# Experimental and numerical investigation of accretion inception and heat transfer physics in ice crystal icing

Yasir Malik, Lokman Bennani, Stephan Bansmer, Pierre Trontin, Philippe Villedieu

## ► To cite this version:

Yasir Malik, Lokman Bennani, Stephan Bansmer, Pierre Trontin, Philippe Villedieu. Experimental and numerical investigation of accretion inception and heat transfer physics in ice crystal icing. International Journal of Heat and Mass Transfer, 2023, 214, pp.124364. 10.1016/j.ijheatmasstransfer.2023.124364. hal-04192890

HAL Id: hal-04192890

https://hal.science/hal-04192890

Submitted on 10 Jun 2024

HAL is a multi-disciplinary open access archive for the deposit and dissemination of scientific research documents, whether they are published or not. The documents may come from teaching and research institutions in France or abroad, or from public or private research centers.

L'archive ouverte pluridisciplinaire HAL, est destinée au dépôt et à la diffusion de documents scientifiques de niveau recherche, publiés ou non, émanant des établissements d'enseignement et de recherche français ou étrangers, des laboratoires publics ou privés.

<figure><img src="image_1.png" /></figure>

Distributed under a Creative Commons Attribution 4.0 International License
<figure><img src="image_1.png" /></figure>

# Experimental and numerical investigation of accretion inception and heat transfer physics in ice crystal icing

<figure><img src="image_2.png" /></figure>

Yasir A. Malik<sup>a,1,*</sup>, Lokman Bennani<sup>b,2,**</sup>, Stephan Bansmer<sup>c,3</sup>, Pierre Trontin<sup>d,4</sup>, Philippe Villedieu<sup>b,5</sup>

## ARTICLE INFO

**Article history:**

Received 22 February 2023

Revised 5 May 2023

Accepted 29 May 2023

Available online 11 June 2023

**Keywords:**

Aircraft icing

Engine icing

Ice crystal icing

Ice inception

Multiphase flow

Ice accretion

## ABSTRACT

Earlier research highlighted severe detrimental effects related to aircraft performance due to ice crystal icing. As part of current research impact of glaciated icing cloud on heated substrates is investigated. To that end, a series of experiments were carried out at the icing wind tunnel of Technische Universitt Braunschweig to investigate the physics of ice crystal accretion on heated substrates for various operating conditions. Dedicated transparent and metallic heatable substrates were designed for macro and microscopic investigation of ice accretion from a qualitative visualization and heat transfer physics perspective. In addition to that numerical simulations were performed at ONERA. Qualitative observations from dedicated experiments on ice accretion initiation phase resulted in advancement of numerical model to accurately capture the different stages leading up to ice accretion as well as the necessary conditions for growth of slushy and glaciated ice layers. Furthermore, a dedicated test matrix was defined to explicitly study the influence of dominant parameters such as flow velocity, heat flux, wet bulb temperature and ice water content on ice accretion process. The experiments showcased a strong influence of increasing flow velocity and ice water content yielding shorter duration required to accrete an ice layer on a heatable substrate. Further experimental investigation reflected that upon increasing the heating power of the test article the icing cloud had to overcome a larger temperature gradient resulting in longer duration required for accretion. It was found that on one hand, the test run with negative wet bulb temperature required a heating source from the substrate as the necessary condition for ice accretion resulting in glaciated ice layers. On the other hand, for positive wet bulb temperature cases natural melting of ice layer was sufficient to induce ice accretion generating slushy ice layers and the heating source from the substrate had little to no influence on the overall ice accretion growth. Numerical simulations were also performed for same conditions and were able to correctly capture the trends and orders of magnitude in comparison with experimental results. The experimental findings presented in this paper lead to development of an accretion solver based on enthalpy approach which considers accretion process as a homogeneous mixture of crystals and liquid water and not by a mere superposition of ice layer and liquid water film. The findings helped calibrate, validate and advance the numerical model for ice crystal icing and found to be more representative of accretion process than unsteady triple layer approach. The findings ensure better predictive capability resulting in improved flight safety and performance criterion.

© 2023 The Authors. Published by Elsevier Ltd.

This is an open access article under the CC BY license (http://creativecommons.org/licenses/by/4.0/)

<footnote><sup>a</sup> Institute of Fluid Mechanics, Technische Universität Braunschweig, Hermann-Blenk-Str. 37, Braunschweig, Germany</footnote>

<footnote><sup>b</sup> ONERA / DMPE, Université de Toulouse, Toulouse F-31055, France</footnote>

<footnote><sup>c</sup> Dipartimento di Science e Tecnologie Aerospaziali, Politecnico di Milano, Via La Masa 34, Milano 20156, Italy</footnote>

<footnote><sup>d</sup> Univ Lyon, Univ Claude Bernard Lyon 1, CNRS, Ecole Centrale de Lyon, INSA Lyon, LMFA, UMR5509, Villeurbanne 69622, France</footnote>

<footnote>* Principal Corresponding (Lead) Author.</footnote>

<footnote>** Corresponding author.</footnote>

<footnote>E-mail addresses: y.malik@tu-braunschweig.de (Y.A. Malik), Lokman.Bennani@onera.fr (L. Bennani), pierre.trontin@univ-lyon1.fr (P. Trontin), philippe.villedieu@onera.fr (P. Villedieu).</footnote>

<footnote n="1">1 Research Scientist</footnote>

<footnote n="2">2 Research Engineer</footnote>

<footnote n="3">3 Professor</footnote>

<footnote n="4">4 Professor</footnote>

<footnote n="5">5 Research Director</footnote>
<figure><img src="image_2.png" /></figure>

# Experimental and numerical investigation of accretion inception and heat transfer physics in ice crystal icing

<figure><img src="image_3.png" /></figure>

Yasir A. Malik<sup>a,1,*</sup>, Lokman Bennani<sup>b,2,**</sup>, Stephan Bansmer<sup>c,3</sup>, Pierre Trontin<sup>d,4</sup>,
Philippe Villedieu<sup>b,5</sup>

## ARTICLE INFO

**Article history:**

Received 22 February 2023

Revised 5 May 2023

Accepted 29 May 2023

Available online 11 June 2023

**Keywords:**

Aircraft icing

Engine icing

Ice crystal icing

Ice inception

Multiphase flow

Ice accretion

## ABSTRACT

Earlier research highlighted severe detrimental effects related to aircraft performance due to ice crystal icing. As part of current research impact of glaciated icing cloud on heated substrates is investigated. To that end, a series of experiments were carried out at the icing wind tunnel of Technische Universitt Braunschweig to investigate the physics of ice crystal accretion on heated substrates for various operating conditions. Dedicated transparent and metallic heatable substrates were designed for macro and microscopic investigation of ice accretion from a qualitative visualization and heat transfer physics perspective. In addition to that numerical simulations were performed at ONERA. Qualitative observations from dedicated experiments on ice accretion initiation phase resulted in advancement of numerical model to accurately capture the different stages leading up to ice accretion as well as the necessary conditions for growth of slushy and glaciated ice layers. Furthermore, a dedicated test matrix was defined to explicitly study the influence of dominant parameters such as flow velocity, heat flux, wet bulb temperature and ice water content on ice accretion process. The experiments showcased a strong influence of increasing flow velocity and ice water content yielding shorter duration required to accrete an ice layer on a heatable substrate. Further experimental investigation reflected that upon increasing the heating power of the test article the icing cloud had to overcome a larger temperature gradient resulting in longer duration required for accretion. It was found that on one hand, the test run with negative wet bulb temperature required a heating source from the substrate as the necessary condition for ice accretion resulting in glaciated ice layers. On the other hand, for positive wet bulb temperature cases natural melting of ice layer was sufficient to induce ice accretion generating slushy ice layers and the heating source from the substrate had little to no influence on the overall ice accretion growth. Numerical simulations were also performed for same conditions and were able to correctly capture the trends and orders of magnitude in comparison with experimental results. The experimental findings presented in this paper lead to development of an accretion solver based on enthalpy approach which considers accretion process as a homogeneous mixture of crystals and liquid water and not by a mere superposition of ice layer and liquid water film. The findings helped calibrate, validate and advance the numerical model for ice crystal icing and found to be more representative of accretion process than unsteady triple layer approach. The findings ensure better predictive capability resulting in improved flight safety and performance criterion.

© 2023 The Authors. Published by Elsevier Ltd.

This is an open access article under the CC BY license (http://creativecommons.org/licenses/by/4.0/)

<footnote><sup>a</sup> Institute of Fluid Mechanics, Technische Universität Braunschweig, Hermann-Blenk-Str. 37, Braunschweig, Germany</footnote>

<footnote><sup>b</sup> ONERA / DMPE, Université de Toulouse, Toulouse F-31055, France</footnote>

<footnote><sup>c</sup> Dipartimento di Science e Tecnologie Aerospaziali, Politecnico di Milano, Via La Masa 34, Milano 20156, Italy</footnote>

<footnote><sup>d</sup> Univ Lyon, Univ Claude Bernard Lyon 1, CNRS, Ecole Centrale de Lyon, INSA Lyon, LMFA, UMR5509, Villeurbanne 69622, France</footnote>

<footnote>* Principal Corresponding (Lead) Author.</footnote>

<footnote>** Corresponding author.</footnote>

<footnote>E-mail addresses: y.malik@tu-braunschweig.de (Y.A. Malik), Lokman.Bennani@onera.fr (L. Bennani), pierre.trontin@univ-lyon1.fr (P. Trontin), philippe.villedieu@onera.fr (P. Villedieu).</footnote>

<footnote n="1">1 Research Scientist</footnote>

<footnote n="2">2 Research Engineer</footnote>

<footnote n="3">3 Professor</footnote>

<footnote n="4">4 Professor</footnote>

<footnote n="5">5 Research Director</footnote>
## Nomenclature

<table><tr><td>TUBS</td><td>Technische Universität Braunschweig</td></tr><tr><td>AoA</td><td>Angle of attack of the test article</td></tr><tr><td>ICI</td><td>Ice Crystal Icing</td></tr><tr><td>IGLOO2D</td><td>ONERA's two-dimensional icing simulation tool</td></tr><tr><td>α<sub>l</sub></td><td>Liquid volume fraction</td></tr><tr><td>ρ</td><td>Density [kg.m<sup>-3</sup>]</td></tr><tr><td>e</td><td>Accretion mean energy per unit mass [J.kg<sup>-1</sup>]</td></tr><tr><td>h</td><td>Accretion thickness [m]</td></tr><tr><td>u</td><td>Velocity [m.s<sup>-1</sup>]</td></tr><tr><td>ṁ<sub>drop</sub></td><td>Droplet deposition mass flow rate [kg.m<sup>-2</sup>.s<sup>-1</sup>]</td></tr><tr><td>ṁ<sub>c,l</sub></td><td>Deposition mass flow rate of liquid water brought by crystals [kg.m<sup>-2</sup>.s<sup>-1</sup>]</td></tr><tr><td>ṁ<sub>c,s</sub></td><td>Deposition mass flow rate of ice brought by crystals [kg.m<sup>-2</sup>.s<sup>-1</sup>]</td></tr><tr><td>ṁ<sub>evs</sub></td><td>Evaporation/sublimation mass flow rate [kg.m<sup>-2</sup>.s<sup>-1</sup>]</td></tr><tr><td>ṁ<sub>ero</sub></td><td>Erosion mass flow rate [kg.m<sup>-2</sup>.s<sup>-1</sup>]</td></tr><tr><td>T<sub>drop</sub></td><td>Temperature of the deposited droplets [K]</td></tr><tr><td>T<sub>c</sub></td><td>Temperature of the deposited crystals [K]</td></tr><tr><td>u<sub>drop</sub></td><td>Velocity of the deposited droplets [m.s<sup>-1</sup>]</td></tr><tr><td>u<sub>c</sub></td><td>Velocity of the deposited crystals [m.s<sup>-1</sup>]</td></tr><tr><td>Φ<sub>0</sub><sup>-</sup></td><td>Heat flux transferred by the domain Ω<sub>2</sub> to the domain Ω<sub>1</sub> [W.m<sup>-2</sup>]</td></tr><tr><td>h<sub>tc</sub></td><td>Heat transfer coefficient [W.m<sup>-2</sup>.K<sup>-1</sup>]</td></tr><tr><td>T<sub>f</sub></td><td>Recovery temperature [K]</td></tr><tr><td>η<sub>s</sub></td><td>Enthalpy of the ice (per unit mass) [J.kg<sup>-1</sup>]</td></tr><tr><td>η<sub>l</sub></td><td>Enthalpy of liquid water (per unit mass) [J.kg<sup>-1</sup>]</td></tr><tr><td>η<sub>v</sub></td><td>Enthalpy of steam (per unit mass) [J.kg<sup>-1</sup>]</td></tr><tr><td>c<sub>p,s</sub></td><td>Specific heat of ice [J.kg<sup>-1</sup>.K<sup>-1</sup>]</td></tr><tr><td>c<sub>p,l</sub></td><td>Specific heat liquid water [J.kg<sup>-1</sup>.K<sup>-1</sup>]</td></tr><tr><td>c<sub>p,v</sub></td><td>Specific heat of steam [J.kg<sup>-1</sup>.K<sup>-1</sup>]</td></tr><tr><td>T<sub>m</sub></td><td>Melting temperature [K]</td></tr><tr><td>T<sub>b</sub></td><td>Boiling temperature [K]</td></tr><tr><td>L<sub>f</sub></td><td>Latent heat of fusion [J.kg<sup>-1</sup>]</td></tr><tr><td>L<sub>v</sub></td><td>Latent heat of vaporization [J.kg<sup>-1</sup>]</td></tr><tr><td>T</td><td>Temperature of the ice/liquid water mixture [K]</td></tr><tr><td>λ</td><td>Thermal conductivity [W.K<sup>-1</sup>.m<sup>-1</sup>]</td></tr><tr><td>Θ</td><td>Temperature in the flat plate [K]</td></tr><tr><td>Φ<sub>0</sub><sup>+</sup></td><td>Heat flux transferred by the domain Ω<sub>1</sub> to the domain Ω<sub>2</sub> [W.m<sup>-2</sup>]</td></tr><tr><td>T<sup>+</sup></td><td>Temperature of the ice/liquid water mixture at the interface Γ [K]</td></tr><tr><td>Θ<sup>-</sup></td><td>Temperature of the flat plate at the interface Γ [K]</td></tr><tr><td>Γ</td><td>Interface between the flat plate and the ice/liquid water mixture</td></tr><tr><td>Δt</td><td>Time step [s]</td></tr><tr><td>ṁ<sub>ero</sub></td><td>Erosion mass rate [kg.s<sup>-1</sup>.m<sup>-2</sup>]</td></tr><tr><td>u<sub>c</sub><sup>t</sup></td><td>Tangential component of the ice crystal impact velocity [m.s<sup>-1</sup>]</td></tr><tr><td>ε<sub>s</sub></td><td>Sticking efficiency</td></tr></table>

# 1. Introduction

Icing has been a major hazard in aviation since the beginning. One of the secondary icing hazards that was long overlooked is ice crystal icing. It mostly occurs in strongly convective weather, where high concentrations of moisture are lifted to high altitudes. This leads to formation of very small ice crystals as small as 10µm. In glaciated conditions, where there are only ice crystals and no supercooled droplets, they just bounce off the aircrafts cold surfaces. Due to no ice accretion on the airframe and little turbulence,

the conditions associated with ICI often appear benign to pilots [1]. Unlike supercooled droplet icing, ice crystal icing happens inside the aircrafts engine and on external heated surfaces, such as heated stagnation pressure probes. It is believed that in engines, the ice crystals partially melt due to the increasing air temperature in the engines compressor. This results in mixed phase icing conditions, where both ice crystals and supercooled droplets or partially melted ice crystals are present. In turn this allows ice crystals to stick to the surfaces, that have cooled down to the freezing point upon impact [2]. The phenomenon of jet engine icing is pictorially illustrated in Fig. 1.

Icing in the engine core flow path can lead to un-commanded loss of thrust (rollback), flameout and damage to the engine blades due to ice shedding [1]. Icing on heated stagnation pressure probes can lead to biased pressure measurements and in turn, an erroneous display of airspeed and the loss of the autopilot. In 2009, such an event led to the loss of an Air France flight of an Airbus A330 over the Atlantic Ocean [4]. Engine power loss events were first attributed to ice crystals, when a transport aircraft, equipped with Ice Detectors, encountered power loss in glaciated conditions in 2002 [5]. Since then, there have been multiple projects to investigate the processes associated with icing in glaciated and mixed phase conditions. To counter these challenges, icing wind tunnels around the world have been upgraded with production capabilities that include glaciated and mixed phase icing conditions.

The HAIC project (High Altitude Ice Crystals) worked on numerical and test capabilities as well as ice particle detection technologies in order to lay the groundwork for future aircraft and de-icing technology certification procedures [6]. Its successor, MUSIC-haic, aims to build on the results, fill gaps in knowledge, provide reliable experimental data and develop new and better predictive capabilities. The goal is to further the physical understanding of ice crystal icing and develop numerical capabilities to tackle icing challenges. The goal is to provide missing data for ice crystal accretion on heatable substrates in order to improve, calibrate and validate ice crystal accretion models. This work is a continuation of the work previously performed in the HAIC project as well as by CNRC and NASA independently in recent years which was essentially focused on ice crystal accretion on unheated surfaces.

At TUBS, ice crystal icing experiments were carried out with flow velocities up to 40 m/s for which influence of erosion is limited. Ice crystals were produced using cloud chamber technology. The experimental results obtained from the different setups helped further the understanding of icing physics and ice accretion phenomenon. The analysis of the results reveals intriguing insights not only qualitatively (such as different stages of ice accretion inception and growth) but also quantitatively (temperature evolution). The large database of accretion experiments has been exploited to calibrate, validate and advance numerical models.

At ONERA, modelling and numerical simulations were conducted. A classic approach to simulate ice accretion is given by Messinger's model. This method has been used to predict ice shapes in various droplet icing conditions [7–10] and has been extended to ice crystal icing conditions [11–13]. In addition, it has served as a basis for more complex numerical modelling approaches [14–18]. Nevertheless, Messinger's approach is restricted to steady icing conditions. One possible way to model unsteady cases is to use multi-layer approaches [19–22]. This kind of methodology has been applied to simulate electrothermal ice protection systems in de-icing mode and in supercooled droplet icing conditions [23,24]. A multi-layer approach was also developed by Bucknell et al. to study ice crystal accretion on warm substrates [25]. However, the experimental observations presented in this paper show that a multilayered approach does not always accurately describe the ice growth process, especially during the initiation phase of accretion (which is one the main topics of the present
## Nomenclature

<table><tr><td>TUBS</td><td>Technische Universität Braunschweig</td></tr><tr><td>AoA</td><td>Angle of attack of the test article</td></tr><tr><td>ICI</td><td>Ice Crystal Icing</td></tr><tr><td>IGLOO2D</td><td>ONERA's two-dimensional icing simulation tool</td></tr><tr><td>α<sub>l</sub></td><td>Liquid volume fraction</td></tr><tr><td>ρ</td><td>Density [kg.m<sup>-3</sup>]</td></tr><tr><td>e</td><td>Accretion mean energy per unit mass [J.kg<sup>-1</sup>]</td></tr><tr><td>h</td><td>Accretion thickness [m]</td></tr><tr><td>u</td><td>Velocity [m.s<sup>-1</sup>]</td></tr><tr><td>ṁ<sub>drop</sub></td><td>Droplet deposition mass flow rate [kg.m<sup>-2</sup>.s<sup>-1</sup>]</td></tr><tr><td>ṁ<sub>c,l</sub></td><td>Deposition mass flow rate of liquid water brought by crystals [kg.m<sup>-2</sup>.s<sup>-1</sup>]</td></tr><tr><td>ṁ<sub>c,s</sub></td><td>Deposition mass flow rate of ice brought by crystals [kg.m<sup>-2</sup>.s<sup>-1</sup>]</td></tr><tr><td>ṁ<sub>evs</sub></td><td>Evaporation/sublimation mass flow rate [kg.m<sup>-2</sup>.s<sup>-1</sup>]</td></tr><tr><td>ṁ<sub>ero</sub></td><td>Erosion mass flow rate [kg.m<sup>-2</sup>.s<sup>-1</sup>]</td></tr><tr><td>T<sub>drop</sub></td><td>Temperature of the deposited droplets [K]</td></tr><tr><td>T<sub>c</sub></td><td>Temperature of the deposited crystals [K]</td></tr><tr><td>u<sub>drop</sub></td><td>Velocity of the deposited droplets [m.s<sup>-1</sup>]</td></tr><tr><td>u<sub>c</sub></td><td>Velocity of the deposited crystals [m.s<sup>-1</sup>]</td></tr><tr><td>Φ<sub>0</sub><sup>-</sup></td><td>Heat flux transferred by the domain Ω<sub>2</sub> to the domain Ω<sub>1</sub> [W.m<sup>-2</sup>]</td></tr><tr><td>h<sub>tc</sub></td><td>Heat transfer coefficient [W.m<sup>-2</sup>.K<sup>-1</sup>]</td></tr><tr><td>T<sub>f</sub></td><td>Recovery temperature [K]</td></tr><tr><td>η<sub>s</sub></td><td>Enthalpy of the ice (per unit mass) [J.kg<sup>-1</sup>]</td></tr><tr><td>η<sub>l</sub></td><td>Enthalpy of liquid water (per unit mass) [J.kg<sup>-1</sup>]</td></tr><tr><td>η<sub>v</sub></td><td>Enthalpy of steam (per unit mass) [J.kg<sup>-1</sup>]</td></tr><tr><td>c<sub>p,s</sub></td><td>Specific heat of ice [J.kg<sup>-1</sup>.K<sup>-1</sup>]</td></tr><tr><td>c<sub>p,l</sub></td><td>Specific heat liquid water [J.kg<sup>-1</sup>.K<sup>-1</sup>]</td></tr><tr><td>c<sub>p,v</sub></td><td>Specific heat of steam [J.kg<sup>-1</sup>.K<sup>-1</sup>]</td></tr><tr><td>T<sub>m</sub></td><td>Melting temperature [K]</td></tr><tr><td>T<sub>b</sub></td><td>Boiling temperature [K]</td></tr><tr><td>L<sub>f</sub></td><td>Latent heat of fusion [J.kg<sup>-1</sup>]</td></tr><tr><td>L<sub>v</sub></td><td>Latent heat of vaporization [J.kg<sup>-1</sup>]</td></tr><tr><td>T</td><td>Temperature of the ice/liquid water mixture [K]</td></tr><tr><td>λ</td><td>Thermal conductivity [W.K<sup>-1</sup>.m<sup>-1</sup>]</td></tr><tr><td>Θ</td><td>Temperature in the flat plate [K]</td></tr><tr><td>Φ<sub>0</sub><sup>+</sup></td><td>Heat flux transferred by the domain Ω<sub>1</sub> to the domain Ω<sub>2</sub> [W.m<sup>-2</sup>]</td></tr><tr><td>T<sup>+</sup></td><td>Temperature of the ice/liquid water mixture at the interface Γ [K]</td></tr><tr><td>Θ<sup>-</sup></td><td>Temperature of the flat plate at the interface Γ [K]</td></tr><tr><td>Γ</td><td>Interface between the flat plate and the ice/liquid water mixture</td></tr><tr><td>Δt</td><td>Time step [s]</td></tr><tr><td>ṁ<sub>ero</sub></td><td>Erosion mass rate [kg.s<sup>-1</sup>.m<sup>-2</sup>]</td></tr><tr><td>u<sub>c</sub><sup>t</sup></td><td>Tangential component of the ice crystal impact velocity [m.s<sup>-1</sup>]</td></tr><tr><td>ε<sub>s</sub></td><td>Sticking efficiency</td></tr></table>

## 1. Introduction

Icing has been a major hazard in aviation since the beginning. One of the secondary icing hazards that was long overlooked is ice crystal icing. It mostly occurs in strongly convective weather, where high concentrations of moisture are lifted to high altitudes. This leads to formation of very small ice crystals as small as 10µm. In glaciated conditions, where there are only ice crystals and no supercooled droplets, they just bounce off the aircrafts cold surfaces. Due to no ice accretion on the airframe and little turbulence,

the conditions associated with ICI often appear benign to pilots [1]. Unlike supercooled droplet icing, ice crystal icing happens inside the aircrafts engine and on external heated surfaces, such as heated stagnation pressure probes. It is believed that in engines, the ice crystals partially melt due to the increasing air temperature in the engines compressor. This results in mixed phase icing conditions, where both ice crystals and supercooled droplets or partially melted ice crystals are present. In turn this allows ice crystals to stick to the surfaces, that have cooled down to the freezing point upon impact [2]. The phenomenon of jet engine icing is pictorially illustrated in Fig. 1.

Icing in the engine core flow path can lead to un-commanded loss of thrust (rollback), flameout and damage to the engine blades due to ice shedding [1]. Icing on heated stagnation pressure probes can lead to biased pressure measurements and in turn, an erroneous display of airspeed and the loss of the autopilot. In 2009, such an event led to the loss of an Air France flight of an Airbus A330 over the Atlantic Ocean [4]. Engine power loss events were first attributed to ice crystals, when a transport aircraft, equipped with Ice Detectors, encountered power loss in glaciated conditions in 2002 [5]. Since then, there have been multiple projects to investigate the processes associated with icing in glaciated and mixed phase conditions. To counter these challenges, icing wind tunnels around the world have been upgraded with production capabilities that include glaciated and mixed phase icing conditions.

The HAIC project (High Altitude Ice Crystals) worked on numerical and test capabilities as well as ice particle detection technologies in order to lay the groundwork for future aircraft and de-icing technology certification procedures [6]. Its successor, MUSIC-haic, aims to build on the results, fill gaps in knowledge, provide reliable experimental data and develop new and better predictive capabilities. The goal is to further the physical understanding of ice crystal icing and develop numerical capabilities to tackle icing challenges. The goal is to provide missing data for ice crystal accretion on heatable substrates in order to improve, calibrate and validate ice crystal accretion models. This work is a continuation of the work previously performed in the HAIC project as well as by CNRC and NASA independently in recent years which was essentially focused on ice crystal accretion on unheated surfaces.

At TUBS, ice crystal icing experiments were carried out with flow velocities up to 40 m/s for which influence of erosion is limited. Ice crystals were produced using cloud chamber technology. The experimental results obtained from the different setups helped further the understanding of icing physics and ice accretion phenomenon. The analysis of the results reveals intriguing insights not only qualitatively (such as different stages of ice accretion inception and growth) but also quantitatively (temperature evolution). The large database of accretion experiments has been exploited to calibrate, validate and advance numerical models.

At ONERA, modelling and numerical simulations were conducted. A classic approach to simulate ice accretion is given by Messinger's model. This method has been used to predict ice shapes in various droplet icing conditions [7–10] and has been extended to ice crystal icing conditions [11–13]. In addition, it has served as a basis for more complex numerical modelling approaches [14–18]. Nevertheless, Messinger's approach is restricted to steady icing conditions. One possible way to model unsteady cases is to use multi-layer approaches [19–22]. This kind of methodology has been applied to simulate electrothermal ice protection systems in de-icing mode and in supercooled droplet icing conditions [23,24]. A multi-layer approach was also developed by Bucknell et al. to study ice crystal accretion on warm substrates [25]. However, the experimental observations presented in this paper show that a multilayered approach does not always accurately describe the ice growth process, especially during the initiation phase of accretion (which is one the main topics of the present
<figure><img src="image_3.png" /><figcaption>Fig. 1. Illustration of jet engine accretion in ice crystal icing conditions (reproduced with permission from [3]).</figcaption></figure>

work). Indeed, it does not take into account that accretion is ac-
tually constituted by a homogeneous mixture of crystals and liq-
uid water and not by a superposition of ice layer and liquid water
film. Therefore, in this work an unsteady extension to Messinger's
model, based on the enthalpy method [26], and assuming that ac-
cretion is a mixture of ice and liquid water is used.

1.1. Motivation

Foremost a high speed imaging technique is utilized for mi-
croscopic investigation of ice accretion inception phase. This lays
down the foundation of better understanding of initiation phase
of ice crystal accretion. The experimental findings were meant to
develop an accretion solver based on enthalpy approach which is
expected to be more representative of an accretion process than
unsteady triple layer approach. Next the origins of ice accretion are
briefly presented for different types of ice layers generated in the
TUBS IWT. This should help further the understanding of the nec-
essary conditions required for ice accretion to initiate on heated
and non-heated substrates for positive and negative wet bulb con-
ditions. Furthermore, heat transfer physics is investigated for ice
crystal icing on heated substrate under a wide range of parametric
conditions. The results are to be used for calibration, advancement
and validation of the presented numerical model. Trends and order
of magnitudes are to be compared between experimental runs and
numerical simulations.

2. TU Braunschweig icing wind tunnel

To further the understanding of heat transfer physics of ice
crystal accretion, two different experimental setups were designed
and developed at the Technische Universität Braunschweig in order
to support the improvement, calibration and validation of numeri-
cal models at ONERA.

First some details are shared with regards to Technische Uni-
versität Braunschweig icing wind tunnel with particular focus on
ice crystal generation and conveyance mechanism. In the next sec-

tions special focus is laid on the different experimental setups fol-
lowed up with analysis and key highlights originating from the ex-
perimental investigation towards advancement of numerical tools
for ice crystal icing on heatable substrates.

2.1. *Experimental facility*

Pictorial description of the Technische Universität Braunschweig icing wind tunnel is provided in Fig. 2. Wind tunnel is marked with position (1) and is coloured in blue. Position (8) highlights the ice crystal generation and conveyance system (IGS). The ice crystals are injected in the wind tunnel through a pipe at position (3). At the same position there is a spray bar system with multiple supercooled droplet injector nozzles. Depending on the settings, different sizes of droplets can be produced, allowing for supercooled large droplets (SLD) and small droplets (with LWC ranging from 0.1 g/m³ to 2 g/m³) enabling mixed phase icing conditions. Control room which allows for controlling the tunnel systems is marked with position (2). Position (4) through (7) marks the tunnel refrigeration system allowing for accretion experiments with sustained tunnel cooling.

The Technische Universität Braunschweig Icing Wind Tunnel hosts 0.5 0.5 m² test section cross-sectional area. It can reach flow velocities up to 40 m/s and allows for static air temperature in the range of -25 to +30 [27]. However, it needs to be pointed out here that the upper limit of the achievable static air temperature during ice crystal experiments is limited to +3 as beyond that the ice crystal conveyance system gets clogged up due to warm ambient conditions.

2.2. *Ice crystal generation system*

One of the unique aspects of the Technische Universität Braun-
schweig Icing Wind Tunnel is its ice crystal generation system.
Cloud chamber technology is employed to simulate the natural oc-
curring ice crystal generation process as observed in atmospheric
clouds. This approach involves feeding atomized water droplets
<figure><img src="image_4.png" /><figcaption>Fig. 1. Illustration of jet engine accretion in ice crystal icing conditions (reproduced with permission from [3]).</figcaption></figure>

work). Indeed, it does not take into account that accretion is ac-
tually constituted by a homogeneous mixture of crystals and liq-
uid water and not by a superposition of ice layer and liquid water
film. Therefore, in this work an unsteady extension to Messinger's
model, based on the enthalpy method [26], and assuming that ac-
cretion is a mixture of ice and liquid water is used.

1.1. Motivation

Foremost a high speed imaging technique is utilized for mi-
croscopic investigation of ice accretion inception phase. This lays
down the foundation of better understanding of initiation phase
of ice crystal accretion. The experimental findings were meant to
develop an accretion solver based on enthalpy approach which is
expected to be more representative of an accretion process than
unsteady triple layer approach. Next the origins of ice accretion are
briefly presented for different types of ice layers generated in the
TUBS IWT. This should help further the understanding of the nec-
essary conditions required for ice accretion to initiate on heated
and non-heated substrates for positive and negative wet bulb con-
ditions. Furthermore, heat transfer physics is investigated for ice
crystal icing on heated substrate under a wide range of parametric
conditions. The results are to be used for calibration, advancement
and validation of the presented numerical model. Trends and order
of magnitudes are to be compared between experimental runs and
numerical simulations.

2. TU Braunschweig icing wind tunnel

To further the understanding of heat transfer physics of ice
crystal accretion, two different experimental setups were designed
and developed at the Technische Universität Braunschweig in order
to support the improvement, calibration and validation of numeri-
cal models at ONERA.

First some details are shared with regards to Technische Uni-
versität Braunschweig icing wind tunnel with particular focus on
ice crystal generation and conveyance mechanism. In the next sec-

tions special focus is laid on the different experimental setups fol-
lowed up with analysis and key highlights originating from the ex-
perimental investigation towards advancement of numerical tools
for ice crystal icing on heatable substrates.

2.1. *Experimental facility*

Pictorial description of the Technische Universität Braunschweig icing wind tunnel is provided in Fig. 2. Wind tunnel is marked with position (1) and is coloured in blue. Position (8) highlights the ice crystal generation and conveyance system (IGS). The ice crystals are injected in the wind tunnel through a pipe at position (3). At the same position there is a spray bar system with multiple supercooled droplet injector nozzles. Depending on the settings, different sizes of droplets can be produced, allowing for supercooled large droplets (SLD) and small droplets (with LWC ranging from 0.1 g/m³ to 2 g/m³) enabling mixed phase icing conditions. Control room which allows for controlling the tunnel systems is marked with position (2). Position (4) through (7) marks the tunnel refrigeration system allowing for accretion experiments with sustained tunnel cooling.

The Technische Universität Braunschweig Icing Wind Tunnel hosts 0.5 0.5 m² test section cross-sectional area. It can reach flow velocities up to 40 m/s and allows for static air temperature in the range of -25 to +30 [27]. However, it needs to be pointed out here that the upper limit of the achievable static air temperature during ice crystal experiments is limited to +3 as beyond that the ice crystal conveyance system gets clogged up due to warm ambient conditions.

2.2. *Ice crystal generation system*

One of the unique aspects of the Technische Universität Braun-
schweig Icing Wind Tunnel is its ice crystal generation system.
Cloud chamber technology is employed to simulate the natural oc-
curring ice crystal generation process as observed in atmospheric
clouds. This approach involves feeding atomized water droplets
<figure><img src="image_4.png" /><figcaption>**Fig. 2.** Schematics of Icing Wind Tunnel at Technische Universität Braunschweig [27].</figcaption></figure>

<figure><img src="image_5.png" /><figcaption>**Fig. 3.** Ice Crystal Generation System (Left). Sketch of Technische Universität Braunschweig Ice Crystal Generation Mechanism (Right) [27].</figcaption></figure>

through an air stream into a chamber, in which the conveying air stream is mixed with another cooled airstream. This leads to droplets being frozen and then growing to naturally shaped ice particles, using the ambient vapour and liquid droplets [2,28]. Once the ice crystals reach a critical size they settle down in the chest freezer. The cloud chamber-based employing freeze-out technique for ice crystal production system is depicted in Fig. 3. In terms of production capabilities, the ice crystal generation system produces approximately 25 kg of snow per day [29].

For long experimental runs, ice particles need to be stored inside ultra-high cooling refrigerants operating up to -60 In case of ice crystal storage for long duration, particles undergo sintering, leading to formation of larger clumps of ice particles, that need to be broken down manually before feeding the particle mass to the conveyance system. Therefore, freshly produced snow is mixed with stored snow to ensure good quality of snow being fed to the test section throughout experimental runs.

The conveyance system consists of a volumetric screw conveyor, a sieving machine and an airstream that transports the sieved particles, through an after-cooler into the icing wind tunnel to avoid pre-melting of ice crystals through the pipes as they approach the test section. More details on ice crystal shapes, sizes and cloud uniformity has already been presented by Bansmer et al. [27] and Malik et al. [30].

### 3. Ice accretion inception and origins on a heated wall

To investigate the initiation of ice accretion on a heated wall, a dedicated transparent flat plate setup was designed to study inception of ice accretion from interface and frontal perspective.

Fig. 4(a) shows the pictorial illustration of interface layer perspective to study ice accretion inception. Fig. 4(b) shows the frontal view of the transparent visualization test article, which is composed of an Indium tin oxide coated (600 nm etched to the surface) sapphire glass (3 mm thick) which is electrically heatable. The transparent coating is then further coated only along the edges with silver to enable electrical cabling and ensure the flow of electrical current through the glass. The ITO coated sapphire glass is designed in such a way that it has a resistance of 25 Ohms and maximum operational voltage of 100 Volts, thereby enabling a maximum heating power of 400 Watts and maximum heat flux of approximately 22000 W/m². To avoid any heat transfer loss to the side of the flat plate, a 2 mm thin Teflon body is designed to ensure heat flux dominant direction towards the impinging ice crystals.

Fig. 5 illustrates the setup for visualization of ice accretion inception. The figure displays a mirror that was attached to the aluminum body through a rail setup that enabled the mirror to be moved along the span of the flat plate. The impingement of ice
<figure><img src="image_5.png" /><figcaption>**Fig. 2.** Schematics of Icing Wind Tunnel at Technische Universität Braunschweig [27].</figcaption></figure>

<figure><img src="image_6.png" /><figcaption>**Fig. 3.** Ice Crystal Generation System (Left). Sketch of Technische Universität Braunschweig Ice Crystal Generation Mechanism (Right) [27].</figcaption></figure>

through an air stream into a chamber, in which the conveying air stream is mixed with another cooled airstream. This leads to droplets being frozen and then growing to naturally shaped ice particles, using the ambient vapour and liquid droplets [2,28]. Once the ice crystals reach a critical size they settle down in the chest freezer. The cloud chamber-based employing freeze-out technique for ice crystal production system is depicted in Fig. 3. In terms of production capabilities, the ice crystal generation system produces approximately 25 kg of snow per day [29].

For long experimental runs, ice particles need to be stored inside ultra-high cooling refrigerants operating up to -60 In case of ice crystal storage for long duration, particles undergo sintering, leading to formation of larger clumps of ice particles, that need to be broken down manually before feeding the particle mass to the conveyance system. Therefore, freshly produced snow is mixed with stored snow to ensure good quality of snow being fed to the test section throughout experimental runs.

The conveyance system consists of a volumetric screw conveyor, a sieving machine and an airstream that transports the sieved particles, through an after-cooler into the icing wind tunnel to avoid pre-melting of ice crystals through the pipes as they approach the test section. More details on ice crystal shapes, sizes and cloud uniformity has already been presented by Bansmer et al. [27] and Malik et al. [30].

### 3. Ice accretion inception and origins on a heated wall

To investigate the initiation of ice accretion on a heated wall, a dedicated transparent flat plate setup was designed to study inception of ice accretion from interface and frontal perspective.

Fig. 4(a) shows the pictorial illustration of interface layer perspective to study ice accretion inception. Fig. 4(b) shows the frontal view of the transparent visualization test article, which is composed of an Indium tin oxide coated (600 nm etched to the surface) sapphire glass (3 mm thick) which is electrically heatable. The transparent coating is then further coated only along the edges with silver to enable electrical cabling and ensure the flow of electrical current through the glass. The ITO coated sapphire glass is designed in such a way that it has a resistance of 25 Ohms and maximum operational voltage of 100 Volts, thereby enabling a maximum heating power of 400 Watts and maximum heat flux of approximately 22000 W/m². To avoid any heat transfer loss to the side of the flat plate, a 2 mm thin Teflon body is designed to ensure heat flux dominant direction towards the impinging ice crystals.

Fig. 5 illustrates the setup for visualization of ice accretion inception. The figure displays a mirror that was attached to the aluminum body through a rail setup that enabled the mirror to be moved along the span of the flat plate. The impingement of ice
<figure><img src="image_6.png" /><figcaption>**Fig. 4.** (a) Pictorial illustration of interface layer perspective to study inception phase, (b) Frontal view of transparent ITO coated sapphire model with a mirror rail setup in the back.</figcaption></figure>

<figure><img src="image_7.png" /><figcaption>**Fig. 5.** On the left - Top view of the visualization setup, On the right - Side view of the visualization setup.</figcaption></figure>

crystals and the corresponding melted residue on the heated surface is first projected on to the first surface mirror and then reflected into the Ultra-High-Speed Cam (Photron SAZ Cam 2100K) placed outside right next to the test section.

To ensure the right contrast and uniform background for obtaining sharp boundaries of the impinging ice crystals on the heatable ITO coated glass multiple considerations were accounted for. High powered flicker free light sources were used to illuminate the target i.e., beamers (16000 lm) and focused light source (approximately 2000 lm). For uniform spreading of light, sand blasted glasses were produced that were positioned such that it allowed for optimum uniform brightness. Additionally, a scheimpflug adapter was affixed to the lens and sensor of the high-speed camera in such a way that it enabled capturing the desired contrast, brightness and sharp boundaries of the impinging ice crystals. In addition to the high-speed camera focused on the mirror from the back of the model, an additional High-Speed Cam (Photron SA 1.1) was situated such that it captured the ice crystal inception from the front which provided additional interesting insights explained in Section 3.1. Detailed overview of the experimental setup is presented in Fig. 6.

To capture the qualitative description of the ice accretion inception on the transparent target two additional DSLR cameras (Canon 80D) were placed from the side (at an oblique angle) and from the top view. Special care was taken in terms of selection of the trans-

parent glass and Sapphire was chosen as it not only provides excellent mechanical and thermally conductive properties but also high optical transmission properties making it ideal for IR applications. This allowed to capture the surface temperature of the transparent glass by placing an IR camera focused on the flat plate setup from the bottom window.

### 3.1. Inception of ice layer - Optical observations

To study the ice initiation process a series of experiments were performed to investigate the different stages that led to the buildup of ice layer on a heated substrate. This not only provided detailed qualitative visual description and physical understanding of ice accretion inception phase but also aided in advancement of numerical tools.

Fig. 7 demonstrates how the ice accretion initiates (from interface layer perspective) on a heated substrate as recorded from the ultra-high-speed camera. The impingement of glaciated icing cloud onto a dry and warm substrate lead to heat transfer taking place between the ice crystals and the heated wall. The impact of ice crystals lead to formation of secondary fragments and residue that stick to the surface and start melting to form sessile bubbles or sessile droplets that exist pre-dominantly in liquid phase. The incoming ice crystals impinging on the existing sessile bubbles lead to splashing and intake/absorption of incoming impacting
<figure><img src="image_7.png" /><figcaption>**Fig. 4.** (a) Pictorial illustration of interface layer perspective to study inception phase, (b) Frontal view of transparent ITO coated sapphire model with a mirror rail setup in the back.</figcaption></figure>

<figure><img src="image_8.png" /><figcaption>**Fig. 5.** On the left - Top view of the visualization setup, On the right - Side view of the visualization setup.</figcaption></figure>

crystals and the corresponding melted residue on the heated surface is first projected on to the first surface mirror and then reflected into the Ultra-High-Speed Cam (Photron SAZ Cam 2100K) placed outside right next to the test section.

To ensure the right contrast and uniform background for obtaining sharp boundaries of the impinging ice crystals on the heatable ITO coated glass multiple considerations were accounted for. High powered flicker free light sources were used to illuminate the target i.e., beamers (16000 lm) and focused light source (approximately 2000 lm). For uniform spreading of light, sand blasted glasses were produced that were positioned such that it allowed for optimum uniform brightness. Additionally, a scheimpflug adapter was affixed to the lens and sensor of the high-speed camera in such a way that it enabled capturing the desired contrast, brightness and sharp boundaries of the impinging ice crystals. In addition to the high-speed camera focused on the mirror from the back of the model, an additional High-Speed Cam (Photron SA 1.1) was situated such that it captured the ice crystal inception from the front which provided additional interesting insights explained in Section 3.1. Detailed overview of the experimental setup is presented in Fig. 6.

To capture the qualitative description of the ice accretion inception on the transparent target two additional DSLR cameras (Canon 80D) were placed from the side (at an oblique angle) and from the top view. Special care was taken in terms of selection of the trans-

parent glass and Sapphire was chosen as it not only provides excellent mechanical and thermally conductive properties but also high optical transmission properties making it ideal for IR applications. This allowed to capture the surface temperature of the transparent glass by placing an IR camera focused on the flat plate setup from the bottom window.

### 3.1. Inception of ice layer - Optical observations

To study the ice initiation process a series of experiments were performed to investigate the different stages that led to the buildup of ice layer on a heated substrate. This not only provided detailed qualitative visual description and physical understanding of ice accretion inception phase but also aided in advancement of numerical tools.

Fig. 7 demonstrates how the ice accretion initiates (from interface layer perspective) on a heated substrate as recorded from the ultra-high-speed camera. The impingement of glaciated icing cloud onto a dry and warm substrate lead to heat transfer taking place between the ice crystals and the heated wall. The impact of ice crystals lead to formation of secondary fragments and residue that stick to the surface and start melting to form sessile bubbles or sessile droplets that exist pre-dominantly in liquid phase. The incoming ice crystals impinging on the existing sessile bubbles lead to splashing and intake/absorption of incoming impacting
<figure><img src="image_8.png" /><figcaption>**Fig. 6.** Pictorial representation of the actualized setup for visualization of ice accretion inception.</figcaption></figure>

<figure><img src="image_9.png" /><figcaption>**Fig. 7.** Ice Accretion Inception from the interface layer perspective (Flow Velocity = 30 m/s, Heat Flux = 1700 W/m², IWC = 5 g/m³, T<sub>wb</sub> = -5 °C).</figcaption></figure>

ice crystals which in turn lead to the volumetric enhancement of the sessile drops. Over time with impact of incoming ice crystals lead to formation of semi-frozen state (mixed-phase in composition) which upon percolation with neighbouring semi-frozen particles lead to the formation of ice layer.

The inception of ice accretion with spatial resolution of 512×512 and temporal resolution of 20,000 frames per second revealed that as the number of sessile drops started to increase on the surface it leads to more splashing and erosion of the liquid drops as a result of impact of incoming glaciated ice crystals. It appears that the particle sticking is a dynamic process and that the sticking efficiency varies greatly depending on whether the ice crystals impact on dry heated spot compared to wet spot covered by sessile drops where in the case of latter the sticking efficiency seem to be higher. Therefore, during the inception phase as soon as the heated

surface is covered by sessile drops the sticking efficiency starts to rise leading to rapid formation of an ice layer.

**Fig. 8** shows visual description of ice accretion initiation as perceived from the high-speed camera from the frontal view. It is visually observed that the ice crystals upon impact lead to formation of secondary fragments and these secondary fragments radially flying away from the surface under the influence of air flow and incoming ice crystals impact the heated substrate again and get entrapped in the sessile drops to aid the process of ice accretion. The accretion inception stages as perceived from the interface layer perspective were supported from the frontal view as well, leading all the way up to formation of ice layer. The surface temperature evolution captured through the IR camera also reveals that when the surface is heated i.e., the initial surface temperature is above freezing point, the impinging ice crystals leads to decrease of sub-
<figure><img src="image_9.png" /><figcaption>**Fig. 6.** Pictorial representation of the actualized setup for visualization of ice accretion inception.</figcaption></figure>

<figure><img src="image_10.png" /><figcaption>**Fig. 7.** Ice Accretion Inception from the interface layer perspective (Flow Velocity = 30 m/s, Heat Flux = 1700 W/m², IWC = 5 g/m³, T<sub>wb</sub> = -5 °C).</figcaption></figure>

ice crystals which in turn lead to the volumetric enhancement of the sessile drops. Over time with impact of incoming ice crystals lead to formation of semi-frozen state (mixed-phase in composition) which upon percolation with neighbouring semi-frozen particles lead to the formation of ice layer.

The inception of ice accretion with spatial resolution of 512×512 and temporal resolution of 20,000 frames per second revealed that as the number of sessile drops started to increase on the surface it leads to more splashing and erosion of the liquid drops as a result of impact of incoming glaciated ice crystals. It appears that the particle sticking is a dynamic process and that the sticking efficiency varies greatly depending on whether the ice crystals impact on dry heated spot compared to wet spot covered by sessile drops where in the case of latter the sticking efficiency seem to be higher. Therefore, during the inception phase as soon as the heated

surface is covered by sessile drops the sticking efficiency starts to rise leading to rapid formation of an ice layer.

**Fig. 8** shows visual description of ice accretion initiation as perceived from the high-speed camera from the frontal view. It is visually observed that the ice crystals upon impact lead to formation of secondary fragments and these secondary fragments radially flying away from the surface under the influence of air flow and incoming ice crystals impact the heated substrate again and get entrapped in the sessile drops to aid the process of ice accretion. The accretion inception stages as perceived from the interface layer perspective were supported from the frontal view as well, leading all the way up to formation of ice layer. The surface temperature evolution captured through the IR camera also reveals that when the surface is heated i.e., the initial surface temperature is above freezing point, the impinging ice crystals leads to decrease of sub-
<figure><img src="image_10.png" /><figcaption>Fig. 8. Ice Accretion Inception from the front view perspective (Flow Velocity = 30 m/s, Heat Flux = 1700 W/m², IWC = 5 g/m³, T<sub>wb</sub>= -5 °C).</figcaption></figure>

strate temperature which upon approaching the freezing tempera-
ture results in accretion on the substrate surface.

The key observations retrieved from the visualization setup
serve as a major input for further advancement of numerical and
theoretical icing models. The experimental results show that the
thermal behaviour of the heated flat plate in ice crystal icing con-
ditions is unsteady. Therefore the unsteady triple layer approach
[22,25] is a possible candidate for the description of this config-
uration. In this approach, the ice formation and melting physics
is described by the superposition of at most three single phase
layers of either liquid water or ice. However, the stage preceding
to formation of ice layer is not a complete liquid layer but rather
percolation of semi-frozen particles (mixed-phase in its composi-
tion) bridging together to form an ice layer. Hence a triple deck
approach, such as the one proposed in [22], would not provide
a good description of the experimentally observed ice accretion
physics. In order to make a first step towards the description of
this ice formation phenomena, an extension to Messinger's classic
steady state model is considered and is described in the following
section. It is an enthalpy based approach that introduces the liquid
volume fraction and hence is able to represent locally the presence
of both the liquid and solid phase.

3.2. Origins of ice accretion - Qualitative experimental observations

Origins of ice accretion for positive and negative tunnel wet
d bulb temperatures was investigated experimentally resulting in for-
mation of glaciated and slushy ice layers which is briefly discussed
here. The test article and setup used to study the influence of wet
bulb temperature and substrate heating on the origins and extent
of ice accretion growth is described in Section 4.1. The wet bulb
temperature is described in the literature as the lowest tempera-
ture to which air can be cooled by the evaporation of water into
the air at a constant pressure [31]. For icing conditions this trans-
lates into ice particle melting at T<sub>wb</sub> > 0, whereas T<sub>wb</sub> < 0 affects
film freezing [2]. The wet bulb temperature can be calculated from
the air temperature and the relative humidity using the Stull for-
mula [32].

To study the origins of ice accretion, three different heating
modes were actualized and analyzed wherein substrate heaters
were turned off at different instances of time resulting in an in-
fluence on ice accretion growth. Mode 1 and 2 were necessary to
allow for continuous heating and turning off heaters as soon as
the icing cloud impacts the test article respectively. Mode 3 on
the other hand was to consolidate the findings by allowing the
heaters to remain switched on for 20 seconds after the icing cloud
impacted the test article. The corresponding visual description of
ice accretion growth for the set of experiments carried out in the
TUBS IWT is presented in Fig. 9. The summary of analysis and
key discoveries originating from the set of experiments performed
in three different heating modes are presented in tabular form in
Fig. 10.

It is observed that for accretion to initiate liquid water needed
to be present in the particle impingement area. This source of
water is either produced in the case of negative wet bulb tem-
perature conditions through heat exchange between impinging ice
crystals and heating from the substrate and in case of positive
wet bulb temperature conditions through dominant natural melt-
ing of ice crystals. In case of heated wall, the ice accretion growth
is attributed to melt water produced through melting of ice crys-
tals owing to substrate heating and liquid water capillary rise
through the ice layer. In case of non-heated wall, the accretion
growth is aided by natural melting of the outermost surface of ice
layer as well as partial melting of incoming impinging ice crys-
tals.

It is also very interesting to observe that for positive wet
d bulb temperature cases the maximum ice accretion thickness did
not vary greatly for cases with and without continuous heat-
ing mode as ice accretion growth is primarily driven by natu-
ral melting of ice layers. The nature of ice layers and extent of
maximum accretion thickness was significantly different for the
different type of ice layers, wherein positive wet bulb tempera-
tures lead to significantly larger ice layers with slushy consistence
more vulnerable to shedding unlike strongly adhered glaciated ice
layers generated in test cases with negative wet bulb tempera-
ture.
<figure><img src="image_11.png" /><figcaption>**Fig. 8.** Ice Accretion Inception from the front view perspective (Flow Velocity = 30 m/s, Heat Flux = 1700 W/m², IWC = 5 g/m³, T<sub>wb</sub> = -5 °C).</figcaption></figure>

strate temperature which upon approaching the freezing tempera-
ture results in accretion on the substrate surface.

The key observations retrieved from the visualization setup
serve as a major input for further advancement of numerical and
theoretical icing models. The experimental results show that the
thermal behaviour of the heated flat plate in ice crystal icing con-
ditions is unsteady. Therefore the unsteady triple layer approach
[22,25] is a possible candidate for the description of this config-
uration. In this approach, the ice formation and melting physics
is described by the superposition of at most three single phase
layers of either liquid water or ice. However, the stage preceding
to formation of ice layer is not a complete liquid layer but rather
percolation of semi-frozen particles (mixed-phase in its composi-
tion) bridging together to form an ice layer. Hence a triple deck
approach, such as the one proposed in [22], would not provide
a good description of the experimentally observed ice accretion
physics. In order to make a first step towards the description of
this ice formation phenomena, an extension to Messinger's classic
steady state model is considered and is described in the following
section. It is an enthalpy based approach that introduces the liquid
volume fraction and hence is able to represent locally the presence
of both the liquid and solid phase.

3.2. Origins of ice accretion - Qualitative experimental observations

Origins of ice accretion for positive and negative tunnel wet
d bulb temperatures was investigated experimentally resulting in for-
mation of glaciated and slushy ice layers which is briefly discussed
here. The test article and setup used to study the influence of wet
bulb temperature and substrate heating on the origins and extent
of ice accretion growth is described in Section 4.1. The wet bulb
temperature is described in the literature as the lowest tempera-
ture to which air can be cooled by the evaporation of water into
the air at a constant pressure [31]. For icing conditions this trans-
lates into ice particle melting at T<sub>wb</sub> > 0, whereas T<sub>wb</sub> < 0 affects
film freezing [2]. The wet bulb temperature can be calculated from
the air temperature and the relative humidity using the Stull for-
mula [32].

To study the origins of ice accretion, three different heating
modes were actualized and analyzed wherein substrate heaters
were turned off at different instances of time resulting in an in-
fluence on ice accretion growth. Mode 1 and 2 were necessary to
allow for continuous heating and turning off heaters as soon as
the icing cloud impacts the test article respectively. Mode 3 on
the other hand was to consolidate the findings by allowing the
heaters to remain switched on for 20 seconds after the icing cloud
impacted the test article. The corresponding visual description of
ice accretion growth for the set of experiments carried out in the
TUBS IWT is presented in Fig. 9. The summary of analysis and
key discoveries originating from the set of experiments performed
in three different heating modes are presented in tabular form in
Fig. 10.

It is observed that for accretion to initiate liquid water needed
to be present in the particle impingement area. This source of
water is either produced in the case of negative wet bulb tem-
perature conditions through heat exchange between impinging ice
crystals and heating from the substrate and in case of positive
wet bulb temperature conditions through dominant natural melt-
ing of ice crystals. In case of heated wall, the ice accretion growth
is attributed to melt water produced through melting of ice crys-
tals owing to substrate heating and liquid water capillary rise
through the ice layer. In case of non-heated wall, the accretion
growth is aided by natural melting of the outermost surface of ice
layer as well as partial melting of incoming impinging ice crys-
tals.

It is also very interesting to observe that for positive wet
d bulb temperature cases the maximum ice accretion thickness did
not vary greatly for cases with and without continuous heat-
ing mode as ice accretion growth is primarily driven by natu-
ral melting of ice layers. The nature of ice layers and extent of
maximum accretion thickness was significantly different for the
different type of ice layers, wherein positive wet bulb tempera-
tures lead to significantly larger ice layers with slushy consistence
more vulnerable to shedding unlike strongly adhered glaciated ice
layers generated in test cases with negative wet bulb tempera-
ture.
<figure><img src="image_11.png" /><figcaption>Fig. 9. Influence of Heater turn off on Ice Accretion growth for different tunnel wet bulb temperatures.</figcaption></figure>

<table><caption>Fig. 10. Necessary conditions for ice accretion for negative and positive wet bulb temperature conditions.</caption><thead><tr><th>Mode</th><th>T<sub>wb</sub> > 0</th><th>T<sub>wb</sub> < 0</th></tr></thead><tbody><tr><td>Mode 1: Continuous Heating</td><td>Ice accretion (combined effect of natural melting and heating from the substrate)</td><td>Ice accretion (resulting from the heated substrate)</td></tr><tr><td>Mode 2: Heaters turned off as soon as the ice cloud is activated</td><td><ul><li>Ice accretion initiated by thermal inertia inside substrate</li><li>Ice accretion does not stop<ul><li>Liquid water necessary for the ice particles to stick comes from the melting of the accreted ice layer since T<sub>wb</sub> > 0</li></ul></li></ul></td><td><ul><li>Ice accretion initiated by thermal inertia inside substrate</li><li>After a while, ice accretion stops as soon as the thin layer of ice is formed<ul><li>Non-availability of source of liquid water for ice accretion</li></ul></li></ul></td></tr><tr><td>Mode 3: Heaters turned off when the ice cloud is activated + 20 seconds</td><td><ul><li>More ice accretion than in mode 2<ul><li>Imbibition highlighting (liquid water capillary rise)</li></ul></li></ul></td><td><ul><li>More ice accretion than in mode 2<ul><li>Imbibition highlighting (liquid water capillary rise)</li></ul></li><li>After a while, accretion stops<ul><li>No more liquid water available for accretion</li></ul></li></ul></td></tr></tbody></table>

## 4. Heat transfer physics on a metallic heatable substrate

In this section, first the experimental setup is described which provides a detailed overview of investigation of heat transfer physics for ice crystal icing on a heatable substrate applicable for wide range of parametric conditions. Next the numerical setup is presented which provides a holistic overview of the modeling strategy adopted to capture the desired trends and order of magnitude for a considerably large parametric space.

### 4.1. Experimental setup

To study heat transfer physics and ice accretion phenomenon on a heatable metallic substrate a dedicated test article was designed for testing in the icing wind tunnel of Technische Universität Braunschweig as shown in Fig. 11.

The flat plate model is composed of multiple layers with base support structure (a) responsible for providing structural strength and transportation of cables to outside of the test section. To avoid any heat leakage to the sides Teflon strips are used and for insu-
<figure><img src="image_12.png" /><figcaption>Fig. 9. Influence of Heater turn off on Ice Accretion growth for different tunnel wet bulb temperatures.</figcaption></figure>

<table><caption>Fig. 10. Necessary conditions for ice accretion for negative and positive wet bulb temperature conditions.</caption><thead><tr><th>Mode</th><th>T<sub>wb</sub> > 0</th><th>T<sub>wb</sub> < 0</th></tr></thead><tbody><tr><td>Mode 1: Continuous Heating</td><td>Ice accretion (combined effect of natural melting and heating from the substrate)</td><td>Ice accretion (resulting from the heated substrate)</td></tr><tr><td>Mode 2: Heaters turned off as soon as the ice cloud is activated</td><td><ul><li>Ice accretion initiated by thermal inertia inside substrate</li><li>Ice accretion does not stop<ul><li>Liquid water necessary for the ice particles to stick comes from the melting of the accreted ice layer since T<sub>wb</sub> > 0</li></ul></li></ul></td><td><ul><li>Ice accretion initiated by thermal inertia inside substrate</li><li>After a while, ice accretion stops as soon as the thin layer of ice is formed<ul><li>Non-availability of source of liquid water for ice accretion</li></ul></li></ul></td></tr><tr><td>Mode 3: Heaters turned off when the ice cloud is activated + 20 seconds</td><td><ul><li>More ice accretion than in mode 2<ul><li>Imbibition highlighting (liquid water capillary rise)</li></ul></li></ul></td><td><ul><li>More ice accretion than in mode 2<ul><li>Imbibition highlighting (liquid water capillary rise)</li></ul></li><li>After a while, accretion stops<ul><li>No more liquid water available for accretion</li></ul></li></ul></td></tr></tbody></table>

## 4. Heat transfer physics on a metallic heatable substrate

In this section, first the experimental setup is described which provides a detailed overview of investigation of heat transfer physics for ice crystal icing on a heatable substrate applicable for wide range of parametric conditions. Next the numerical setup is presented which provides a holistic overview of the modeling strategy adopted to capture the desired trends and order of magnitude for a considerably large parametric space.

### 4.1. Experimental setup

To study heat transfer physics and ice accretion phenomenon on a heatable metallic substrate a dedicated test article was designed for testing in the icing wind tunnel of Technische Universität Braunschweig as shown in Fig. 11.

The flat plate model is composed of multiple layers with base support structure (a) responsible for providing structural strength and transportation of cables to outside of the test section. To avoid any heat leakage to the sides Teflon strips are used and for insu-
<figure><img src="image_12.png" /><figcaption>**Fig. 11.** Construction of test article with heaters, temperature and heat flux sensors.</figcaption></figure>

lation in the bottom direction, Multitherm Insulation 550M (b) is placed which has very high insulation properties (thermal conductivity coefficient = 0.023 W/m.K). This ensures heat flux only in one direction i.e., towards the impinging ice crystals. Grooves (c) are installed for the heat flux sensors as well as temperature sensors to have them completely flushed with the surface and ensure a smooth flat finish.

To provide the heating capability, a total of six heat foils are glued to the aluminium plate (d). The heat foils are placed such that they are able to provide uniform heating to the flat plate model, each heat foil has a resistance of 14.4 Ohms, a maximum operational voltage of 28 Volts, able to provide a total heating power of 326 Watts and maximum heat flux of around 15000 W/m². To account for the evolution of temperature and heat flux as a result of impinging ice crystals leading to ice accretion, a total of 34 temperature sensors (PT1000 sensors) and 5 heat flux sensors (Differential-Temperature Thermopile) were installed 0.3 mm to the ice crystal target surface (e). Cameras from oblique angle were installed to obtain qualitative description of ice accretion growth under varying parametric conditions. For generation of ice layers under different regimes (Slushy and Glaciated ice layers) focus was set on influential parameters such as heat flux, tunnel temperature, ice water content and flow velocity. The experiments were performed in such a way that the surface is first allowed to reach thermal equilibrium and employment of a shield that covers the model allowing for stabilization of icing cloud (for 30 seconds) inside the test section.

For experiments 08 to 10 presented in Fig. 12, repeatability of icing experiments in TUBS IWT was analyzed with regards to temperature evolution upon impact of icing cloud on heated flat plate test article. The experimental conditions for these runs are heating power of 50 Watts, tunnel temperature of -10 ice water content of 5 g/m³ and flow velocity of 40 m/s with test article oriented orthogonal to the impinging icing cloud. The excellent repeatability of icing experiments for same experimental conditions ensures high confidence interval and allows for robust investigation of accretion phenomenon resulting from impact of glaciated ice crystals on heated substrates.

As can be seen in Fig. 12, the test article reaches thermal equilibrium and at time $t = 70$ seconds, the icing shield is removed which allows impact of ice crystal on heated flat plate. The ice accretion initiates over the model within the first 20 seconds of impact of ice crystals and this is shown in the form of sharp temperature gradient, as surface temperature reaches close to the freezing

<figure><img src="image_13.png" /><figcaption>**Fig. 12.** Heat Transfer Physics - Temperature evolution pre and post-impact icing cloud.</figcaption></figure>

point. It is important to point out here that as soon as the icing cloud impacts the heated substrate it leads to the transient phase which encompasses the different stages of ice inception that lead to the formation of a coherent ice layer as described in Section 3.1.

To monitor temperature evolution for a wide array of temperature sensors on the flat plate model, we can clearly see in Fig. 12 that the ice accretion initiates as soon as the surface temperature reaches close to the freezing point which is in line with the findings of the previous studies [33].

An overview of the experimental setup is illustrated in Fig. 13, Cameras (1) and (3) are used to obtain qualitative description of ice accretion growth from top and side oblique perspective. IR Camera is illustrated as (2) used to capture the surface temperature evolution. Light sources (4) are placed from the top and bottom window to illuminate the test article and allow for better observation of ice accretion initiation and growth during the runs. The setup allows for accurately measuring the surface temperature evolution as well as qualitative description of ice accretion growth process.

## 4.2. Numerical setup

### 4.2.1. Simulation framework

In this work, modelling and simulation are performed in a two-dimensional framework. Fig. 14 shows the simulation frame-
<figure><img src="image_13.png" /><figcaption>**Fig. 11.** Construction of test article with heaters, temperature and heat flux sensors.</figcaption></figure>

lation in the bottom direction, Multitherm Insulation 550M (b) is placed which has very high insulation properties (thermal conductivity coefficient = 0.023 W/m.K). This ensures heat flux only in one direction i.e., towards the impinging ice crystals. Grooves (c) are installed for the heat flux sensors as well as temperature sensors to have them completely flushed with the surface and ensure a smooth flat finish.

To provide the heating capability, a total of six heat foils are glued to the aluminium plate (d). The heat foils are placed such that they are able to provide uniform heating to the flat plate model, each heat foil has a resistance of 14.4 Ohms, a maximum operational voltage of 28 Volts, able to provide a total heating power of 326 Watts and maximum heat flux of around 15000 W/m². To account for the evolution of temperature and heat flux as a result of impinging ice crystals leading to ice accretion, a total of 34 temperature sensors (PT1000 sensors) and 5 heat flux sensors (Differential-Temperature Thermopile) were installed 0.3 mm to the ice crystal target surface (e). Cameras from oblique angle were installed to obtain qualitative description of ice accretion growth under varying parametric conditions. For generation of ice layers under different regimes (Slushy and Glaciated ice layers) focus was set on influential parameters such as heat flux, tunnel temperature, ice water content and flow velocity. The experiments were performed in such a way that the surface is first allowed to reach thermal equilibrium and employment of a shield that covers the model allowing for stabilization of icing cloud (for 30 seconds) inside the test section.

For experiments 08 to 10 presented in Fig. 12, repeatability of icing experiments in TUBS IWT was analyzed with regards to temperature evolution upon impact of icing cloud on heated flat plate test article. The experimental conditions for these runs are heating power of 50 Watts, tunnel temperature of -10 ice water content of 5 g/m³ and flow velocity of 40 m/s with test article oriented orthogonal to the impinging icing cloud. The excellent repeatability of icing experiments for same experimental conditions ensures high confidence interval and allows for robust investigation of accretion phenomenon resulting from impact of glaciated ice crystals on heated substrates.

As can be seen in Fig. 12, the test article reaches thermal equilibrium and at time $t = 70$ seconds, the icing shield is removed which allows impact of ice crystal on heated flat plate. The ice accretion initiates over the model within the first 20 seconds of impact of ice crystals and this is shown in the form of sharp temperature gradient, as surface temperature reaches close to the freezing

<figure><img src="image_14.png" /><figcaption>**Fig. 12.** Heat Transfer Physics - Temperature evolution pre and post-impact icing cloud.</figcaption></figure>

point. It is important to point out here that as soon as the icing cloud impacts the heated substrate it leads to the transient phase which encompasses the different stages of ice inception that lead to the formation of a coherent ice layer as described in Section 3.1.

To monitor temperature evolution for a wide array of temperature sensors on the flat plate model, we can clearly see in Fig. 12 that the ice accretion initiates as soon as the surface temperature reaches close to the freezing point which is in line with the findings of the previous studies [33].

An overview of the experimental setup is illustrated in Fig. 13, Cameras (1) and (3) are used to obtain qualitative description of ice accretion growth from top and side oblique perspective. IR Camera is illustrated as (2) used to capture the surface temperature evolution. Light sources (4) are placed from the top and bottom window to illuminate the test article and allow for better observation of ice accretion initiation and growth during the runs. The setup allows for accurately measuring the surface temperature evolution as well as qualitative description of ice accretion growth process.

## 4.2. Numerical setup

### 4.2.1. Simulation framework

In this work, modelling and simulation are performed in a two-dimensional framework. Fig. 14 shows the simulation frame-
<figure><img src="image_14.png" /><figcaption>Fig. 13. Illustration of equipment used to observe ice accretion on heated flat plate.</figcaption></figure>

<figure><img src="image_15.png" /><figcaption>Fig. 14. IGLOO2D workflow. $\rho_e$, $u_e$, $p_e$ and $T_e$ are respectively the density, velocity, pressure and temperature outside the boundary layer that develops around the airfoil. $\dot{m}_{imp}$ is the mass flow rate of impacting particles and $T_{imp}^{tot}$ is their total temperature. $\tau$, $h_{tc}$ and $T_{rec}$ are respectively the aerodynamic shear stress at the wall, the heat transfer coefficient and the recovery temperature.</figcaption></figure>

work provided by IGLOO2D, ONERA's two-dimensional icing simulation tool. A detailed description of the ONERA icing tool IGLOO2D is given in [34]. In this work, the Euler2D inviscid flow solver was used in combination with the CLICET Prandlt boundary layer solver to perform the computation of the aerodynamic field. The lagrangian solver TrajL2D is used to compute the ice crystals' trajectories. Data concerning the heat transfer coefficient, the collection efficiency and the heating power calibration can be found in Appendix A.

To perform unsteady computations, as is the case in this study, IGLOO2D typically uses an unsteady triple layer ice accretion solver that is coupled with an unsteady heat conduction solver. The coupling is performed with an optimized Schwarz algorithm [23]. However, as stated in the introduction and the previous section, based on the observations from the experiments (Section 3.1), it seems that the triple deck approach proposed in [22] does not accurately describe what is observed at early times with the formation of a semi frozen layer containing both liquid and solid water on the wall. The superposition of the three layers is replaced here by a single mixed phase layer characterized by its liquid water volume fraction $\alpha_l$. Therefore, a new solver based on an enthalpy method was developed and will be presented in a following section.

### 4.2.2. Geometry and test matrix

Fig. 15 shows the flat plate geometry with its dimensions and the thermal conductivity of each layer. A temperature probe located mid-span and in the middle of the thermally conductive adhesive layer is used to compare the numerical results with the experimental data. The physical properties of the materials composing the flat plate are given in Table 1.

### 4.2.3. Accretion model and heat conduction in the flat plate

As illustrated in Fig. 16, the problem of ice accretion on a heated plate can be represented using two separate domains. The first one, noted $\Omega_1$, represents the ice/liquid water mixture. The second one, noted $\Omega_2$, represents the heated plate. Each domain is modelled with a specific set of equations which are coupled at the interface $\Gamma$.

The outer boundary of $\Omega_1$, noted $\Gamma_{ext}$, is in contact with the aerodynamic flow. Hence, that is where the convective heat transfer and evaporation/sublimation occurs. Moreover, ice crystals and droplets may impinge that boundary. The equations solved for the liquid/solid mixture describe ice accretion, that is to say heat and mass transfer, in domain $\Omega_1$. They are given by:

$$ \frac{\partial}{\partial t}(\rho h) + \nabla \cdot (\rho h \mathbf{u}) = \dot{m}_{drop} + \dot{m}_{c,l} + \dot{m}_{c,s} - \dot{m}_{evs} - \dot{m}_{ero} \quad (1a) $$
<figure><img src="image_15.png" /><figcaption>Fig. 13. Illustration of equipment used to observe ice accretion on heated flat plate.</figcaption></figure>

<figure><img src="image_16.png" /><figcaption>Fig. 14. IGLOO2D workflow. $\rho_e$, $u_e$, $p_e$ and $T_e$ are respectively the density, velocity, pressure and temperature outside the boundary layer that develops around the airfoil. $\dot{m}_{imp}$ is the mass flow rate of impacting particles and $T_{imp}^{tot}$ is their total temperature. $\tau$, $h_{tc}$ and $T_{rec}$ are respectively the aerodynamic shear stress at the wall, the heat transfer coefficient and the recovery temperature.</figcaption></figure>

work provided by IGLOO2D, ONERA's two-dimensional icing simulation tool. A detailed description of the ONERA icing tool IGLOO2D is given in [34]. In this work, the Euler2D inviscid flow solver was used in combination with the CLICET Prandlt boundary layer solver to perform the computation of the aerodynamic field. The lagrangian solver TrajL2D is used to compute the ice crystals' trajectories. Data concerning the heat transfer coefficient, the collection efficiency and the heating power calibration can be found in Appendix A.

To perform unsteady computations, as is the case in this study, IGLOO2D typically uses an unsteady triple layer ice accretion solver that is coupled with an unsteady heat conduction solver. The coupling is performed with an optimized Schwarz algorithm [23]. However, as stated in the introduction and the previous section, based on the observations from the experiments (Section 3.1), it seems that the triple deck approach proposed in [22] does not accurately describe what is observed at early times with the formation of a semi frozen layer containing both liquid and solid water on the wall. The superposition of the three layers is replaced here by a single mixed phase layer characterized by its liquid water volume fraction $\alpha_l$. Therefore, a new solver based on an enthalpy method was developed and will be presented in a following section.

### 4.2.2. Geometry and test matrix

Fig. 15 shows the flat plate geometry with its dimensions and the thermal conductivity of each layer. A temperature probe located mid-span and in the middle of the thermally conductive adhesive layer is used to compare the numerical results with the experimental data. The physical properties of the materials composing the flat plate are given in Table 1.

### 4.2.3. Accretion model and heat conduction in the flat plate

As illustrated in Fig. 16, the problem of ice accretion on a heated plate can be represented using two separate domains. The first one, noted $\Omega_1$, represents the ice/liquid water mixture. The second one, noted $\Omega_2$, represents the heated plate. Each domain is modelled with a specific set of equations which are coupled at the interface $\Gamma$.

The outer boundary of $\Omega_1$, noted $\Gamma_{ext}$, is in contact with the aerodynamic flow. Hence, that is where the convective heat transfer and evaporation/sublimation occurs. Moreover, ice crystals and droplets may impinge that boundary. The equations solved for the liquid/solid mixture describe ice accretion, that is to say heat and mass transfer, in domain $\Omega_1$. They are given by:

$$ \frac{\partial}{\partial t}(\rho h) + \nabla \cdot (\rho h \mathbf{u}) = \dot{m}_{drop} + \dot{m}_{c,l} + \dot{m}_{c,s} - \dot{m}_{evs} - \dot{m}_{ero} \quad (1a) $$
<figure><img src="image_16.png" /><figcaption>Fig. 15. Flat plate setup with dimensions and thermal conductivities.</figcaption></figure>

<table><caption>**Table 1**

Properties of the materials composing the flat plate.</caption><thead><tr><th>Material</th><th>ρ [kg.m<sup>-3</sup>]</th><th>c<sub>p</sub> [J.kg<sup>-1</sup>.K<sup>-1</sup>]</th><th>λ [W.m<sup>-1</sup>.K<sup>-1</sup>]</th></tr></thead><tbody><tr><td>Aluminum</td><td>2700</td><td>900</td><td>237</td></tr><tr><td>Teflon</td><td>2200</td><td>1000</td><td>0.25</td></tr><tr><td>Thermally conductive adhesive</td><td>2500</td><td>900</td><td>5.8</td></tr><tr><td>Kapton heaters</td><td>2700</td><td>900</td><td>0.8</td></tr><tr><td>Multitherm insulation</td><td>2800</td><td>787</td><td>0.023</td></tr></tbody></table>

<figure><img src="image_17.png" /><figcaption>Fig. 16. Illustration of the physical problem and definition of the ice/liquid water and airfoil domains and their coupling interface $\Gamma$.</figcaption></figure>

$$ \frac{\partial}{\partial t}(\rho h e) + \nabla \cdot (\rho h e \mathbf{u}) = \dot{m}_{\text{drop}} \left[ \eta_l(T_{\text{drop}}) + \frac{u_{\text{drop}}^2}{2} \right] + \dot{m}_{c,l} \left[ \eta_l(T_c) + \frac{u_c^2}{2} \right] + \dot{m}_{c,s} \left[ \eta_s(T_c) + \frac{u_c^2}{2} \right] - \dot{m}_{\text{evs}}\eta_v(T) - \dot{m}_{\text{ero}}e + h_{tc}(T_r - T) + \Phi_0^{-} \quad (1b) $$

The system of Eqs. (1) can be derived from the 3D equations modeling the total mass conservation (equation for $\rho h$) and total energy conservation (equation for $\rho e$) and by integrating these equations with respect to the z variable. The divergence operator $\nabla$ is defined with respect to the wall surface coordinate system. $h$ is the mixture thickness, $\mathbf{u}$ its travel velocity and $e$ its energy per unit mass (or enthalpy since they are very close for both water and ice). $\dot{m}_{\text{drop}}$ is the deposition mass flow rate of droplets

whereas $\dot{m}_{c,l}$ and $\dot{m}_{c,s}$ are respectively the deposition mass flow rates of liquid water and ice brought by crystals. The mass flow rate for evaporation and sublimation is given by $\dot{m}_{\text{evs}}$ and the mass flow rate of erosion is given by $\dot{m}_{\text{ero}}$. Temperatures of the deposited droplets and ice crystals are respectively given by $T_{\text{drop}}$ and $T_c$. The same convention is used for the velocity magnitudes $u_{\text{drop}}$ and $u_c$. $\Phi_0$ is the heat flux transferred by the domain $\Omega_2$ to the domain $\Omega_1$. Finally, $T$ is the mean temperature of the accretion layer.

The travel velocity $\mathbf{u}$ is computed using a thin film lubrication model in the case where the liquid volume fraction $\alpha_l$ is equal to 1 (purely liquid state). If the liquid volume fraction is strictly less than 1, this velocity is set to $0 \text{ m.s}^{-1}$. This model assumes that in the presence of ice crystals in the accretion the liquid water is entirely trapped in the pores and cannot be transported by aerodynamic or capillary forces.

Choosing the solid ice at the melting temperature ($T_m = 273 \text{ K}$) as the reference state for the energies, the enthalpies $\eta_s$, $\eta_l$ and $\eta_v$ for the solid ice, liquid water and steam can be written:

$$ \eta_s(T) = c_{p,s}(T - T_m) \quad (2a) $$
<figure><img src="image_17.png" /><figcaption>Fig. 15. Flat plate setup with dimensions and thermal conductivities.</figcaption></figure>

<table><caption>**Table 1**

Properties of the materials composing the flat plate.</caption><thead><tr><th>Material</th><th>ρ [kg.m<sup>-3</sup>]</th><th>c<sub>p</sub> [J.kg<sup>-1</sup>.K<sup>-1</sup>]</th><th>λ [W.m<sup>-1</sup>.K<sup>-1</sup>]</th></tr></thead><tbody><tr><td>Aluminum</td><td>2700</td><td>900</td><td>237</td></tr><tr><td>Teflon</td><td>2200</td><td>1000</td><td>0.25</td></tr><tr><td>Thermally conductive adhesive</td><td>2500</td><td>900</td><td>5.8</td></tr><tr><td>Kapton heaters</td><td>2700</td><td>900</td><td>0.8</td></tr><tr><td>Multitherm insulation</td><td>2800</td><td>787</td><td>0.023</td></tr></tbody></table>

<figure><img src="image_18.png" /><figcaption>Fig. 16. Illustration of the physical problem and definition of the ice/liquid water and airfoil domains and their coupling interface $\Gamma$.</figcaption></figure>

$$ \frac{\partial}{\partial t}(\rho h e) + \nabla \cdot (\rho h e \mathbf{u}) = \dot{m}_{\text{drop}} \left[ \eta_l(T_{\text{drop}}) + \frac{u_{\text{drop}}^2}{2} \right] \\ + \dot{m}_{c,l} \left[ \eta_l(T_c) + \frac{u_c^2}{2} \right] \\ + \dot{m}_{c,s} \left[ \eta_s(T_c) + \frac{u_c^2}{2} \right] \\ - \dot{m}_{\text{evs}}\eta_v(T) - \dot{m}_{\text{ero}}e \\ + h_{tc}(T_r - T) \\ + \Phi_0^{-} \quad (1b) $$

The system of Eqs. (1) can be derived from the 3D equations modeling the total mass conservation (equation for $\rho h$) and total energy conservation (equation for $\rho e$) and by integrating these equations with respect to the z variable. The divergence operator $\nabla$ is defined with respect to the wall surface coordinate system. $h$ is the mixture thickness, $\mathbf{u}$ its travel velocity and $e$ its energy per unit mass (or enthalpy since they are very close for both water and ice). $\dot{m}_{\text{drop}}$ is the deposition mass flow rate of droplets

whereas $\dot{m}_{c,l}$ and $\dot{m}_{c,s}$ are respectively the deposition mass flow rates of liquid water and ice brought by crystals. The mass flow rate for evaporation and sublimation is given by $\dot{m}_{\text{evs}}$ and the mass flow rate of erosion is given by $\dot{m}_{\text{ero}}$. Temperatures of the deposited droplets and ice crystals are respectively given by $T_{\text{drop}}$ and $T_c$. The same convention is used for the velocity magnitudes $u_{\text{drop}}$ and $u_c$. $\Phi_0$ is the heat flux transferred by the domain $\Omega_2$ to the domain $\Omega_1$. Finally, $T$ is the mean temperature of the accretion layer.

The travel velocity $\mathbf{u}$ is computed using a thin film lubrication model in the case where the liquid volume fraction $\alpha_l$ is equal to 1 (purely liquid state). If the liquid volume fraction is strictly less than 1, this velocity is set to $0 \text{ m.s}^{-1}$. This model assumes that in the presence of ice crystals in the accretion the liquid water is entirely trapped in the pores and cannot be transported by aerodynamic or capillary forces.

Choosing the solid ice at the melting temperature ($T_m = 273 \text{ K}$) as the reference state for the energies, the enthalpies $\eta_s$, $\eta_l$ and $\eta_v$ for the solid ice, liquid water and steam can be written:

$$ \eta_s(T) = c_{p,s}(T - T_m) \quad (2a) $$
<table><caption>**Table 2**

Physical properties of water.</caption><thead><tr><th>ρ [kg.m<sup>-3</sup>]</th><th>c<sub>p,s</sub> [J.kg<sup>-1</sup>.K<sup>-1</sup>]</th><th>c<sub>p,l</sub> [J.kg<sup>-1</sup>.K<sup>-1</sup>]</th><th>c<sub>p,v</sub> [J.kg<sup>-1</sup>.K<sup>-1</sup>]</th><th>L<sub>f</sub> (T<sub>m</sub>) [J.kg<sup>-1</sup>]</th><th>L<sub>v</sub> (T<sub>b</sub>) [J.kg<sup>-1</sup>]</th></tr></thead><tbody><tr><td>997</td><td>2060</td><td>4185</td><td>1850</td><td>333000</td><td>2257000</td></tr></tbody></table>

$$ \eta_l(T) = c_{p,l}(T - T_m) + L_f(T_m) \quad (2b) $$

$$ \eta_v(T) = c_{p,v}(T - T_b) + c_{p,l}(T_b - T_m) + L_f(T_m) + L_v(T_b) \quad (2c) $$

where $c_{p,l}$ and $c_{p,s}$ are respectively the heat capacities for liquid water and solid ice. The latent heats for melting and vaporization are given by $L_f$ and $L_v$. Their values are summarized in Table 2. In Eq. (1b) and neglecting the difference between energy and enthalpy for both ice and liquid water, the temperature $T$ and the energy $e$ are related by the following relationship:

$$ \begin{array}{llll} T = T_m + \frac{e}{c_{p,s}} & \text{and} & \alpha_l = 0 & \text{if} & e \le 0 \\ T = T_m & \text{and} & \alpha_l = \frac{e}{L_f(T_m)} & \text{if} & 0 \le \frac{e}{L_f(T_m)} \le 1 \\ T = T_m + \frac{L_f(T_m)}{c_{p,l}} \left( \frac{e}{L_f(T_m)} - 1 \right) & \text{and} & \alpha_l = 1 & \text{if} & 1 \le \frac{e}{L_f(T_m)} \end{array} \quad (3) $$

The thermal conduction inside the flat plate (domain $\Omega_1$ in Fig. 16) is described using the heat equation and is given by:

$$ \frac{\partial}{\partial t} (\rho c_p \Theta) - \nabla \cdot (\lambda \nabla \Theta) = S \quad (4) $$

where $\rho$, $c_p$ and $\lambda$ are respectively the density, the heat capacity and the thermal conductivity of the material composing the wall. The source term $S$ allows taking into account heater mats that can be embedded into the plate. In domain $\Omega_2$, $\Theta$ is used to denote the temperature field.

The boundary conditions are given by:

$$ \lambda \nabla \Theta \cdot n = \Phi_0^+ \quad \text{on } \Gamma \quad (5) $$

Equation (5) stands for the boundary conditions at the bottom of the inner boundary and at the interface with the liquid/solid mixture. $\Phi_0^+$ is the heat flux transferred by the domain $\Omega_1$ to the domain $\Omega_2$.

ONERA's icing simulation tool IGLOO2D is used to compute the aerodynamics, droplet and crystal trajectories (which are used to compute the impacting mass flow rates) as well as the heat transfer coefficient [34]. Moreover, the evaporation mass rate is computed using the classic Chilton-Colburn analogy. To complete the model, one needs to specify an erosion model and a sticking efficiency model.

### 4.2.4. Discretization and coupling procedure

In order to numerically solve the governing equations in the ice/liquid water mixture and in the flat plate, a discretization method as well as a coupling algorithm are required.

Equations (4), (5) and (1) are coupled to ensure the continuity of both the temperature and the heat flux at the interface $\Gamma$ between the liquid/solid mixture and the heated wall:

$$ T^+ = \Theta^- \quad (6a) $$

$$ \Phi_0^+ = -\Phi_0^- \quad (6b) $$

To do this, a coupling procedure derived from [23] is used. The idea is to perform sub-iterations within a time step, here noted with the superscript $k$, until equations (6) are satisfied to a given

tolerance. To this end, the heat fluxes at the interface are expressed in the following way:

$$ \Phi_0^{-(k)} = \omega(\Theta^{-(k-1)} - T^{+(k)}) + \Phi_0^{+(k-1)} \quad (7a) $$

$$ \Phi_0^{+(k)} = \omega(T^{+(k)} - \Theta^{-(k)}) - \Phi_0^{-(k)} \quad (7b) $$

where $\omega$ is a numerical coupling coefficient. In this work its value is fixed to 10000. The sub-stepping iterations are initialised by using the values of the temperatures and fluxes at time step $n$.

At each sub-iteration $k$, a discrete form of Eqs. (1a) and (1b) is solved. They are discretized using a finite volume approach. The considered surface is discretized into control volumes (with only one cell in the direction normal to the surface) and on any given control volume one obtains the following discrete equations:

$$ \rho \frac{h^{(k)} - h^n}{\Delta t} = \dot{m}_{drop} + \dot{m}_{c,l} + \dot{m}_{c,s} - \dot{m}_{evs}^n - \dot{m}_{ero}^n + \frac{F_{m,L}^n - F_{m,R}^n}{\Delta x} \quad (8a) $$

$$ \begin{aligned} \rho \frac{h^{(k)} e^{(k)} - h^n e^n}{\Delta t} = & \dot{m}_{drop} \left[ \eta_l(T_{drop}) + \frac{u_{drop}^2}{2} \right] \\ & + \dot{m}_{c,l} \left[ \eta_l(T_c) + \frac{u_c^2}{2} \right] \\ & + \dot{m}_{c,s} \left[ \eta_s(T_c) + \frac{u_c^2}{2} \right] \\ & - \dot{m}_{evs}^n \eta_v(T^n) - \dot{m}_{ero}^n e^n \\ & + h_{tc} (T_r - T^{(k)}) + \frac{F_{e,L}^n - F_{e,R}^n}{\Delta x} \\ & + \Phi_0^{-(k)} \end{aligned} \quad (8b) $$

Where $h^n$ is the average thickness of the accretion layer in a given cell at time step $n$, $e^n$ the average energy in the cell and $T^n$ the average temperature in the cell. $\Delta x$ is the length of the cell and $\Delta t$ is the time step. $F_{L}^m, F_{R}^m, F_{L}^e, F_{R}^e$ are the numerical mass and energy fluxes at the left and right boundaries of the cell. They are computed using a first order upwind scheme.

On the other hand the heat equation in the plate is discretized using a $\mathbb{P}_1$ finite element method. This yields the following global linear system:

$$ M \frac{\Theta^{(k)} - \Theta^n}{\Delta t} + A \Theta^{(k)} = R \quad (9) $$

where $M$ is the “mass” matrix, $A$ the “stiffness” matrix and $R$ the right hand side which typically includes the heat provided by the heater mats.

To close the coupling methodology, a relationship is needed between $T^+$ and the corresponding average temperature of the ice layer $T$ which is the only available information on the temperature profile in the accretion shape, since Eqs. (1) has been obtained by eliminating the variable $z$. To do so, a uniform temperature (in the normal direction) is assumed in the whole accretion thickness. This leads to $T^{+, (k)} = T^{(k)}$.

On the other hand, domain $\Omega_2$ is meshed with triangles using gmsh [35]. This domain was meshed using 23000 elements. Eq. (4) is then solved numerically using Lagrange $\mathbb{P}_1$ finite elements. The meshes of domains $\Omega_1$ and $\Omega_2$ are constructed so as to be coincident.
<table><caption>**Table 2**

Physical properties of water.</caption><thead><tr><th>ρ [kg.m<sup>-3</sup>]</th><th>c<sub>p,s</sub> [J.kg<sup>-1</sup>.K<sup>-1</sup>]</th><th>c<sub>p,l</sub> [J.kg<sup>-1</sup>.K<sup>-1</sup>]</th><th>c<sub>p,v</sub> [J.kg<sup>-1</sup>.K<sup>-1</sup>]</th><th>L<sub>f</sub> (T<sub>m</sub>) [J.kg<sup>-1</sup>]</th><th>L<sub>v</sub> (T<sub>b</sub>) [J.kg<sup>-1</sup>]</th></tr></thead><tbody><tr><td>997</td><td>2060</td><td>4185</td><td>1850</td><td>333000</td><td>2257000</td></tr></tbody></table>

$$ \eta_l(T) = c_{p,l}(T - T_m) + L_f(T_m) \quad (2b) $$

$$ \eta_v(T) = c_{p,v}(T - T_b) + c_{p,l}(T_b - T_m) + L_f(T_m) + L_v(T_b) \quad (2c) $$

where $c_{p,l}$ and $c_{p,s}$ are respectively the heat capacities for liquid water and solid ice. The latent heats for melting and vaporization are given by $L_f$ and $L_v$. Their values are summarized in Table 2. In Eq. (1b) and neglecting the difference between energy and enthalpy for both ice and liquid water, the temperature $T$ and the energy $e$ are related by the following relationship:

$$ \begin{array}{llll} T = T_m + \frac{e}{c_{p,s}} & \text{and} & \alpha_l = 0 & \text{if} & e \le 0 \\ T = T_m & \text{and} & \alpha_l = \frac{e}{L_f(T_m)} & \text{if} & 0 \le \frac{e}{L_f(T_m)} \le 1 \\ T = T_m + \frac{L_f(T_m)}{c_{p,l}} \left( \frac{e}{L_f(T_m)} - 1 \right) & \text{and} & \alpha_l = 1 & \text{if} & 1 \le \frac{e}{L_f(T_m)} \end{array} \quad (3) $$

The thermal conduction inside the flat plate (domain $\Omega_1$ in Fig. 16) is described using the heat equation and is given by:

$$ \frac{\partial}{\partial t} (\rho c_p \Theta) - \nabla \cdot (\lambda \nabla \Theta) = S \quad (4) $$

where $\rho$, $c_p$ and $\lambda$ are respectively the density, the heat capacity and the thermal conductivity of the material composing the wall. The source term $S$ allows taking into account heater mats that can be embedded into the plate. In domain $\Omega_2$, $\Theta$ is used to denote the temperature field.

The boundary conditions are given by:

$$ \lambda \nabla \Theta \cdot n = \Phi_0^+ \quad \text{on } \Gamma \quad (5) $$

Equation (5) stands for the boundary conditions at the bottom of the inner boundary and at the interface with the liquid/solid mixture. $\Phi_0^+$ is the heat flux transferred by the domain $\Omega_1$ to the domain $\Omega_2$.

ONERA's icing simulation tool IGLOO2D is used to compute the aerodynamics, droplet and crystal trajectories (which are used to compute the impacting mass flow rates) as well as the heat transfer coefficient [34]. Moreover, the evaporation mass rate is computed using the classic Chilton-Colburn analogy. To complete the model, one needs to specify an erosion model and a sticking efficiency model.

### 4.2.4. Discretization and coupling procedure

In order to numerically solve the governing equations in the ice/liquid water mixture and in the flat plate, a discretization method as well as a coupling algorithm are required.

Equations (4), (5) and (1) are coupled to ensure the continuity of both the temperature and the heat flux at the interface $\Gamma$ between the liquid/solid mixture and the heated wall:

$$ T^+ = \Theta^- \quad (6a) $$

$$ \Phi_0^+ = -\Phi_0^- \quad (6b) $$

To do this, a coupling procedure derived from [23] is used. The idea is to perform sub-iterations within a time step, here noted with the superscript $k$, until equations (6) are satisfied to a given

tolerance. To this end, the heat fluxes at the interface are expressed in the following way:

$$ \Phi_0^{-(k)} = \omega(\Theta^{-(k-1)} - T^{+(k)}) + \Phi_0^{+(k-1)} \quad (7a) $$

$$ \Phi_0^{+(k)} = \omega(T^{+(k)} - \Theta^{-(k)}) - \Phi_0^{-(k)} \quad (7b) $$

where $\omega$ is a numerical coupling coefficient. In this work its value is fixed to 10000. The sub-stepping iterations are initialised by using the values of the temperatures and fluxes at time step $n$.

At each sub-iteration $k$, a discrete form of Eqs. (1a) and (1b) is solved. They are discretized using a finite volume approach. The considered surface is discretized into control volumes (with only one cell in the direction normal to the surface) and on any given control volume one obtains the following discrete equations:

$$ \rho \frac{h^{(k)} - h^n}{\Delta t} = \dot{m}_{drop} + \dot{m}_{c,l} + \dot{m}_{c,s} - \dot{m}_{evs}^n - \dot{m}_{ero}^n + \frac{F_{m,L}^n - F_{m,R}^n}{\Delta x} \quad (8a) $$

$$ \begin{aligned} \rho \frac{h^{(k)} e^{(k)} - h^n e^n}{\Delta t} = & \dot{m}_{drop} \left[ \eta_l(T_{drop}) + \frac{u_{drop}^2}{2} \right] \\ & + \dot{m}_{c,l} \left[ \eta_l(T_c) + \frac{u_c^2}{2} \right] \\ & + \dot{m}_{c,s} \left[ \eta_s(T_c) + \frac{u_c^2}{2} \right] \\ & - \dot{m}_{evs}^n \eta_v(T^n) - \dot{m}_{ero}^n e^n \\ & + h_{tc} (T_r - T^{(k)}) + \frac{F_{e,L}^n - F_{e,R}^n}{\Delta x} \\ & + \Phi_0^{-(k)} \end{aligned} \quad (8b) $$

Where $h^n$ is the average thickness of the accretion layer in a given cell at time step $n$, $e^n$ the average energy in the cell and $T^n$ the average temperature in the cell. $\Delta x$ is the length of the cell and $\Delta t$ is the time step. $F_{L}^m, F_{R}^m, F_{L}^e, F_{R}^e$ are the numerical mass and energy fluxes at the left and right boundaries of the cell. They are computed using a first order upwind scheme.

On the other hand the heat equation in the plate is discretized using a $\mathbb{P}_1$ finite element method. This yields the following global linear system:

$$ M \frac{\Theta^{(k)} - \Theta^n}{\Delta t} + A \Theta^{(k)} = R \quad (9) $$

where $M$ is the “mass” matrix, $A$ the “stiffness” matrix and $R$ the right hand side which typically includes the heat provided by the heater mats.

To close the coupling methodology, a relationship is needed between $T^+$ and the corresponding average temperature of the ice layer $T$ which is the only available information on the temperature profile in the accretion shape, since Eqs. (1) has been obtained by eliminating the variable $z$. To do so, a uniform temperature (in the normal direction) is assumed in the whole accretion thickness. This leads to $T^{+(k)} = T^{(k)}$.

On the other hand, domain $\Omega_2$ is meshed with triangles using gmsh [35]. This domain was meshed using 23000 elements. Eq. (4) is then solved numerically using Lagrange $\mathbb{P}_1$ finite elements. The meshes of domains $\Omega_1$ and $\Omega_2$ are constructed so as to be coincident.
<figure><img src="image_18.png" /><figcaption>(a) Test 07</figcaption></figure>

$$Q=100W - T_{\infty}=268.15K$$
$$IWC=5g.m^{-3} - U_{\infty}=20m.s^{-1}$$

<figure><img src="image_19.png" /><figcaption>(b) Test 08</figcaption></figure>

$$Q=100W - T_{\infty}=268.15K$$

$$IWC=5g.m^{-3} - U_{\infty}=30m.s^{-1}$$

<figure><img src="image_20.png" /><figcaption>(c) Test 04</figcaption></figure>

$$Q=100W - T_{\infty}=268.15K$$
$$IWC=5g.m^{-3} - U_{\infty}=40m.s^{-1}$$

<figure><img src="image_21.png" /><figcaption>(d) Plot of sticking efficiency function $F_2$</figcaption></figure>

<caption>Fig. 17. Calibration results for the cases relative to the influence of velocity and sticking efficiency function $F_2$.</caption>

### 4.2.5. Erosion model

The erosion model used in this work is based on previous studies [12,36]. The model links the erosion mass rate to the impacting ice crystal mass rate through a factor $\epsilon_{ero}$. In the present work, the crystals are assumed fully glaciated, yielding:

$$\dot{m}_{ero} = \epsilon_{ero}\dot{m}_{c,s}^{imp} \qquad (10)$$

The factor $\epsilon_{ero}$ is given by:

$$\epsilon_{ero} = E \left( \frac{u_c^t}{V_0} \right)^2 \frac{\alpha_{l0}}{\alpha_{l0} - \min(\alpha_l, \alpha_{l0})} \qquad (11)$$

where $E = 0.314466$, $V_0 = 84.5m.s^{-1}$ and $\alpha_{l0} = 0.6$. $u_c^t$ is the tangential component of the ice crystal impact velocity. It is computed at every time step to take into account the change in geometry of the growing ice shape.

### 4.2.6. Sticking efficiency model

The sticking efficiency model relates the deposition mass rate of ice crystals to the impacting mass rate of ice crystals by using a factor $\epsilon_s$, which yields:

$$\dot{m}_{c,s} = \epsilon_s \dot{m}_{c,s}^{imp} \qquad (12)$$

$$\dot{m}_{c,l} = \epsilon_s \dot{m}_{c,l}^{imp} \qquad (13)$$

In order to compute $\epsilon_s$ several functions are combined. First, the modelling of the sticking of partially melted ice crystals is based on previous work by Trontin et al. [12] and relies on the introduction of a function $F_1$.

$$F_1(\zeta_m) = (K_1 - 2)\zeta_m^3 + (3 - 2K_1)\zeta_m^2 + K_1\zeta_m \qquad (14)$$

where $\zeta_m$ is the melt ratio of the ice crystals and $K_1 = 2.5$.

On the other hand, the sticking of glaciated particles is modelled using a function noted $F_2$ that is given by:

$$F_2(\alpha_l) = 0.02 + 0.75\alpha_l + \frac{1}{2} \left( K_2 \left( \frac{u_0}{u_c} \right)^{0.5} - 3 \right) \alpha_l^2 + \frac{1}{3} \left( 2.25 - K_2 \left( \frac{u_0}{u_c} \right)^{0.5} \right) \alpha_l^3 \qquad (15)$$

where $K_2$ is an empirical parameter equal to 1.89, $u_0 = 40m.s^{-1}$, $\alpha_l$ is the liquid volume fraction, $u_c$ is the magnitude of the crystal impact velocity. This empirical relation was obtained in the following way. First, for liquid volume fractions ranging from 0 to 0.3
<figure><img src="image_19.png" /></figure>

$$IWC = 5g.m^{-3} - U_{\infty} = 20m.s^{-1}$$

<figure><img src="image_20.png" /></figure>

$$IWC = 5g.m^{-3} - U_{\infty} = 30m.s^{-1}$$

<figure><img src="image_21.png" /><figcaption>Fig. 17. Calibration results for the cases relative to the influence of velocity and sticking efficiency function $F_2$.</figcaption></figure>

<figure><img src="image_22.png" /></figure>

### 4.2.5. Erosion model

The erosion model used in this work is based on previous studies [12,36]. The model links the erosion mass rate to the impacting ice crystal mass rate through a factor $\epsilon_{ero}$. In the present work, the crystals are assumed fully glaciated, yielding:

$$\dot{m}_{ero} = \epsilon_{ero}\dot{m}_{c,s}^{imp} \quad (10)$$

The factor $\epsilon_{ero}$ is given by:

$$\epsilon_{ero} = E \left( \frac{u_c^t}{V_0} \right)^2 \frac{\alpha_{l0}}{\alpha_{l0} - \min(\alpha_l, \alpha_{l0})} \quad (11)$$

where $E = 0.314466$, $V_0 = 84.5m.s^{-1}$ and $\alpha_{l0} = 0.6$. $u_c^t$ is the tangential component of the ice crystal impact velocity. It is computed at every time step to take into account the change in geometry of the growing ice shape.

### 4.2.6. Sticking efficiency model

The sticking efficiency model relates the deposition mass rate of ice crystals to the impacting mass rate of ice crystals by using a factor $\epsilon_s$, which yields:

$$\dot{m}_{c,s} = \epsilon_s \dot{m}_{c,s}^{imp} \quad (12)$$

$$\dot{m}_{c,l} = \epsilon_s \dot{m}_{c,l}^{imp} \quad (13)$$

In order to compute $\epsilon_s$ several functions are combined. First, the modelling of the sticking of partially melted ice crystals is based on previous work by Trontin et al. [12] and relies on the introduction of a function $F_1$.

$$F_1(\zeta_m) = (K_1 - 2)\zeta_m^3 + (3 - 2K_1)\zeta_m^2 + K_1\zeta_m \quad (14)$$

where $\zeta_m$ is the melt ratio of the ice crystals and $K_1 = 2.5$.

On the other hand, the sticking of glaciated particles is modelled using a function noted $F_2$ that is given by:

$$F_2(\alpha_l) = 0.02 + 0.75\alpha_l + \frac{1}{2} \left( K_2 \left( \frac{u_0}{u_c} \right)^{0.5} - 3 \right) \alpha_l^2 + \frac{1}{3} \left( 2.25 - K_2 \left( \frac{u_0}{u_c} \right)^{0.5} \right) \alpha_l^3 \quad (15)$$

where $K_2$ is an empirical parameter equal to 1.89, $u_0 = 40m.s^{-1}$, $\alpha_l$ is the liquid volume fraction, $u_c$ is the magnitude of the crystal impact velocity. This empirical relation was obtained in the following way. First, for liquid volume fractions ranging from 0 to 0.3
<figure><img src="image_22.png" /><figcaption>(a) Test 04</figcaption></figure>

<figure><img src="image_23.png" /><figcaption>(c) Test 03</figcaption></figure>

<figure><img src="image_24.png" /><figcaption>(b) Test 05</figcaption></figure>

<figure><img src="image_25.png" /><figcaption>(d) Test 06</figcaption></figure>

<caption>Fig. 18. Comparison between simulation and experimental data for the cases relative to the influence of IWC.</caption>

the data of Baumert et al. [36] provides values for the sticking efficiency. On the other hand, for liquid volume fractions close to 1, the experimental database presented in this paper can be used to estimate the value of the sticking efficiency. Indeed, as shown in Fig., when the ice crystals start impacting the heated plate, the temperature decreases from its thermal equilibrium value in dry conditions to its icing steady state value. During this transient phase, ice accretion has not yet begun and $\alpha_l = 1$. The sticking efficiency will drive the slope of this temperature decrease as it controls the cooling power brought by the ice crystals. Therefore one may estimate the value of the sticking efficiency by calibrating its value so as to obtain a good match for the slope of the temperature curve at the moment when the ice crystals start to impinge the flat plate. In addition, tests 07, 08 and 04 investigate the influence of velocity on the thermal and ice build up behaviour of the flat plate configuration. This parameter has an important influence on the thermal dynamics and hence the sticking efficiency at $\alpha_l = 1$ is assumed to depend on it. To sum up, tests 07, 08 and 04 are used to estimate a value for the sticking efficiency at $\alpha_l = 1$ and for three different impact velocities. Once values have been obtained for the sticking efficiency at $\alpha_l = 1$, a polynomial fit is performed so as to also take into account the data of Baumert et al. [36] for lower values of the liquid volume fraction. Note that the model also includes a residual value of the sticking efficiency equal to 0.02. This value was estimated thanks to experiments dedicated to studying the impact of ice particles on a dry solid substrate [37].

Both functions are combined by taking the maximum value between both:

$$ \epsilon_s = \max(F_1, F_2) \qquad (16) $$

Note that in this work the crystals are assumed to be fully glaciated. Hence, the part of the sticking efficiency model relative to the partial melt ratio of the crystals, $F_1$, will be always equal to 0.

Figs. 17a, 17b and 17c show the results obtained by simulation for tests 07, 08 and 04 where the farfield air velocity was varied from 20m/s to 40m/s. As stated previously, the experimental data from these runs was used to calibrate the sticking efficiency model. This is shown in particular with respect to the slope of the curves at $t = 0$s. The resulting sticking efficiency function is shown in Fig. 17d.

## 5. Results and discussion

A comparison of experimental and numerical results focusing on qualitative observation of ice accretion thickness and temper-
<figure><img src="image_23.png" /></figure>

<figure><img src="image_24.png" /><figcaption>Fig. 18. Comparison between simulation and experimental data for the cases relative to the influence of IWC.</figcaption></figure>

the data of Baumert et al. [36] provides values for the sticking efficiency. On the other hand, for liquid volume fractions close to 1, the experimental database presented in this paper can be used to estimate the value of the sticking efficiency. Indeed, as shown in Fig., when the ice crystals start impacting the heated plate, the temperature decreases from its thermal equilibrium value in dry conditions to its icing steady state value. During this transient phase, ice accretion has not yet begun and $\alpha_l = 1$. The sticking efficiency will drive the slope of this temperature decrease as it controls the cooling power brought by the ice crystals. Therefore one may estimate the value of the sticking efficiency by calibrating its value so as to obtain a good match for the slope of the temperature curve at the moment when the ice crystals start to impinge the flat plate. In addition, tests 07, 08 and 04 investigate the influence of velocity on the thermal and ice build up behaviour of the flat plate configuration. This parameter has an important influence on the thermal dynamics and hence the sticking efficiency at $\alpha_l = 1$ is assumed to depend on it. To sum up, tests 07, 08 and 04 are used to estimate a value for the sticking efficiency at $\alpha_l = 1$ and for three different impact velocities. Once values have been obtained for the sticking efficiency at $\alpha_l = 1$, a polynomial fit is performed so as to also take into account the data of Baumert et al. [36] for lower values of the liquid volume fraction. Note that the model also includes a residual value of the sticking efficiency equal to 0.02. This value was estimated thanks to experiments dedicated to studying the impact of ice particles on a dry solid substrate [37].

<figure><img src="image_25.png" /></figure>

<figure><img src="image_26.png" /></figure>

Both functions are combined by taking the maximum value between both:

$$ \epsilon_s = \max(F_1, F_2) \qquad (16) $$

Note that in this work the crystals are assumed to be fully glaciated. Hence, the part of the sticking efficiency model relative to the partial melt ratio of the crystals, $F_1$, will be always equal to 0.

Figs. 17a, 17b and 17c show the results obtained by simulation for tests 07, 08 and 04 where the farfield air velocity was varied from 20m/s to 40m/s. As stated previously, the experimental data from these runs was used to calibrate the sticking efficiency model. This is shown in particular with respect to the slope of the curves at $t = 0$s. The resulting sticking efficiency function is shown in Fig. 17d.

## 5. Results and discussion

A comparison of experimental and numerical results focusing on qualitative observation of ice accretion thickness and temper-
<figure><img src="image_26.png" /></figure>

<figure><img src="image_27.png" /><figcaption>**Fig. 19.** Comparison between simulation and experimental data for the cases relative to the influence of heating power.</figcaption></figure>

<figure><img src="image_28.png" /></figure>

<figure><img src="image_29.png" /></figure>

ature evolution concerning heat transfer physics is presented and analyzed next.

## 5.1. Influence of ice water content

Fig. 18 shows the results obtained for tests 04, 05, 03 and 06. This group of tests investigates the influence of IWC. The simulations are in quite good agreement with the experimental data. However, the numerical simulations show a much higher effect than what is observed experimentally. Indeed, the numerical simulations predict that the temperature drop occurs faster when the IWC is increased. This can be explained in the following way: as $\alpha_l = 1$ during the temperature drop phase, the sticking efficiency is maximal. Increasing the IWC hence means increasing the cooling power induced by the impacting crystals and therefore higher slopes. The fact that this is not observed experimentally may be due to a shading effect linked to a part of the crystals bouncing back from the flat plate's surface and not taken into account in the model.

It is important to highlight here that in Fig. 18 for the experimental runs the temperature continues to drop before the impact of icing cloud whereas in numerical simulations it has reached thermal equilibrium for the same duration of the run. This is because in the case of experiments for some cases to allow for a stabilized icing cloud in the test section a shield is placed in front of

the test article for the first 40 seconds after the turn on of the icing cloud. This ensures a fully developed and stabilized icing cloud in the test section as previously presented by Bansmer et al. [27]. The consequence of the use of the shield is that it results in formation of flow structures in the wake of shield that results in slight reduction (i.e. cooling effect) of surface temperature upon impact on the test article in experiments unlike numerical simulations which do not employ such a shield since the icing cloud is considered stabilized from the beginning.

## 5.2. Influence of heating power

Fig. 19 shows the results obtained for tests 28, 29, 30 and 31. For this group of tests the tunnel temperature is fixed to 263.15K and the power is varied from 50W to 200W. The numerical simulations predict that when the heating power is increased, the dry air equilibrium temperature is increased as well as the time required to reach the ice crystal icing equilibrium temperature. This trend is consistent with the experimental measurements. On a more quantitative side, for low heating powers such as 50W and 100W, the results are in quite good agreement with the experimental data. The temperature drop rate is overestimated by the experiments but the maximum difference between simulation and experiment is of about 5K. On the other hand for heating powers of 150W and 200W the discrepancy starts to increase.
<figure><img src="image_27.png" /></figure>

<figure><img src="image_28.png" /><figcaption>Fig. 19. Comparison between simulation and experimental data for the cases relative to the influence of heating power.</figcaption></figure>

<figure><img src="image_29.png" /></figure>

<figure><img src="image_30.png" /></figure>

ature evolution concerning heat transfer physics is presented and analyzed next.

## 5.1. Influence of ice water content

Fig. 18 shows the results obtained for tests 04, 05, 03 and 06. This group of tests investigates the influence of IWC. The simulations are in quite good agreement with the experimental data. However, the numerical simulations show a much higher effect than what is observed experimentally. Indeed, the numerical simulations predict that the temperature drop occurs faster when the IWC is increased. This can be explained in the following way: as $\alpha_l = 1$ during the temperature drop phase, the sticking efficiency is maximal. Increasing the IWC hence means increasing the cooling power induced by the impacting crystals and therefore higher slopes. The fact that this is not observed experimentally may be due to a shading effect linked to a part of the crystals bouncing back from the flat plate's surface and not taken into account in the model.

It is important to highlight here that in Fig. 18 for the experimental runs the temperature continues to drop before the impact of icing cloud whereas in numerical simulations it has reached thermal equilibrium for the same duration of the run. This is because in the case of experiments for some cases to allow for a stabilized icing cloud in the test section a shield is placed in front of

the test article for the first 40 seconds after the turn on of the icing cloud. This ensures a fully developed and stabilized icing cloud in the test section as previously presented by Bansmer et al. [27]. The consequence of the use of the shield is that it results in formation of flow structures in the wake of shield that results in slight reduction (i.e. cooling effect) of surface temperature upon impact on the test article in experiments unlike numerical simulations which do not employ such a shield since the icing cloud is considered stabilized from the beginning.

## 5.2. Influence of heating power

Fig. 19 shows the results obtained for tests 28, 29, 30 and 31. For this group of tests the tunnel temperature is fixed to 263.15K and the power is varied from 50W to 200W. The numerical simulations predict that when the heating power is increased, the dry air equilibrium temperature is increased as well as the time required to reach the ice crystal icing equilibrium temperature. This trend is consistent with the experimental measurements. On a more quantitative side, for low heating powers such as 50W and 100W, the results are in quite good agreement with the experimental data. The temperature drop rate is overestimated by the experiments but the maximum difference between simulation and experiment is of about 5K. On the other hand for heating powers of 150W and 200W the discrepancy starts to increase.
<figure><img src="image_30.png" /><figcaption>(a) Test 25</figcaption></figure>

<figure><img src="image_31.png" /><figcaption>(b) Test 26</figcaption></figure>

<figure><img src="image_32.png" /><figcaption>(c) Test 27</figcaption></figure>

<figure><img src="image_33.png" /><figcaption>(d) Test 32</figcaption></figure>

<caption>**Fig. 20.** Comparison between simulation and experimental data for the cases relative to the influence of tunnel temperature.</caption>

## 5.3. Influence of tunnel temperature

Fig. 20 shows the results obtained for tests 25, 26, 27 and 32. For this group of tests the tunnel temperature was varied from 273.15K to 253.15K. In order to always start with the same surface temperature, the heating power was varied from 100W to 158W. A reasonably good agreement is observed for tests 25, 26 and 27. For test 32, the experimental measurement shows that the temperature first decreases to reach approximately 280K at $t = 40s$, after which the temperature then starts to increase.

## 5.4. Qualitative comparison

Figs. 21a and 22b show a side view of the flat plate with ice accreted on top. The direction of airflow is from top to bottom. As stated in Section 4.2, the simulations presented in this work are two-dimensional. Hence, interpreting the numerical results shown in Fig. 21b and 22b as cuts of the corresponding experimental views in the spanwise direction (or, equivalently, in the direction parallel to the airflow) enables a qualitative comparison with experimental observations.

Fig. 21 shows a comparison between experiment and numerical simulation of the ice shape for Test 04. After 100s of ice crystal impact the numerically predicted ice shape is approximately 3mm thick. The temperature in the ice is equal to $T_m = 273.15K$ and the liquid volume fraction is of 0.35. The predicted ice shape covers the whole heated area of the plate and is slightly thicker at the stagnation point which is consistent with the experiment.

Fig. 22 shows a comparison between experiment and numerical simulation of the ice shape for Test 28. After 100s of ice crystal impact the numerically predicted ice shape is approximately 2.5mm thick at the middle of the plate. The thickness quickly reduces to 0mm towards the sides of the plate. At the middle of the plate, the temperature in the ice is equal to $T_m = 273.15K$ and the liquid volume fraction is of 0.15. Towards the side of the plate the temperature goes slightly below 273.15K and at these locations the liquid volume fraction is of 0. In this case, the predicted ice shape only partly covers the heated area of the plate and is thicker at the stagnation point. This is consistent with the experimental visualisation of the ice shape.

For these two cases, the numerical simulations are able to capture the correct trend. That is to say there is less ice for Test 28 than for Test 04. In both cases, the temperature field in the flat
<figure><img src="image_31.png" /><figcaption>(a) Test 25</figcaption></figure>

$$Q=100W - T_{\infty}=273.15K$$
$$IWC=5g.m^{-3} - U_{\infty}=40m.s^{-1}$$

<figure><img src="image_32.png" /><figcaption>(b) Test 26</figcaption></figure>

$$Q=114.5W - T_{\infty}=268.15K$$
$$IWC=5g.m^{-3} - U_{\infty}=40m.s^{-1}$$

<figure><img src="image_33.png" /><figcaption>(c) Test 27</figcaption></figure>

$$Q=130W - T_{\infty}=263.15K$$
$$IWC=5g.m^{-3} - U_{\infty}=40m.s^{-1}$$

<figure><img src="image_34.png" /><figcaption>(d) Test 32</figcaption></figure>

$$Q=145W - T_{\infty}=258.15K$$
$$IWC=5g.m^{-3} - U_{\infty}=40m.s^{-1}$$

<caption>**Fig. 20.** Comparison between simulation and experimental data for the cases relative to the influence of tunnel temperature.</caption>

## 5.3. Influence of tunnel temperature

Fig. 20 shows the results obtained for tests 25, 26, 27 and 32. For this group of tests the tunnel temperature was varied from 273.15K to 253.15K. In order to always start with the same surface temperature, the heating power was varied from 100W to 158W. A reasonably good agreement is observed for tests 25, 26 and 27. For test 32, the experimental measurement shows that the temperature first decreases to reach approximately 280K at $t = 40s$, after which the temperature then starts to increase.

## 5.4. Qualitative comparison

Figs. 21a and 22b show a side view of the flat plate with ice accreted on top. The direction of airflow is from top to bottom. As stated in Section 4.2, the simulations presented in this work are two-dimensional. Hence, interpreting the numerical results shown in Fig. 21b and 22b as cuts of the corresponding experimental views in the spanwise direction (or, equivalently, in the direction parallel to the airflow) enables a qualitative comparison with experimental observations.

Fig. 21 shows a comparison between experiment and numerical simulation of the ice shape for Test 04. After 100s of ice crystal impact the numerically predicted ice shape is approximately 3mm thick. The temperature in the ice is equal to $T_m = 273.15K$ and the liquid volume fraction is of 0.35. The predicted ice shape covers the whole heated area of the plate and is slightly thicker at the stagnation point which is consistent with the experiment.

Fig. 22 shows a comparison between experiment and numerical simulation of the ice shape for Test 28. After 100s of ice crystal impact the numerically predicted ice shape is approximately 2.5mm thick at the middle of the plate. The thickness quickly reduces to 0mm towards the sides of the plate. At the middle of the plate, the temperature in the ice is equal to $T_m = 273.15K$ and the liquid volume fraction is of 0.15. Towards the side of the plate the temperature goes slightly below 273.15K and at these locations the liquid volume fraction is of 0. In this case, the predicted ice shape only partly covers the heated area of the plate and is thicker at the stagnation point. This is consistent with the experimental visualisation of the ice shape.

For these two cases, the numerical simulations are able to capture the correct trend. That is to say there is less ice for Test 28 than for Test 04. In both cases, the temperature field in the flat
<figure><img src="image_34.png" /><figcaption>(a) Experimental visualisation of the ice shape for Test 04 at $t = 100s$</figcaption></figure>

<figure><img src="image_35.png" /><figcaption>(b) Numerical prediction of the ice shape and temperature field for Test 04 at $t = 100s$</figcaption></figure>

<figure><img src="image_36.png" /><figcaption>(a) Experimental visualisation of the ice shape for Test 28 at $t = 100s$</figcaption></figure>

<figure><img src="image_37.png" /><figcaption>(b) Numerical prediction of the ice shape and temperature field for Test 28 at $t = 100s$</figcaption></figure>

<caption>**Fig. 22.** Qualitative comparison between experiment and numerical simulation of the ice shape for Test 28.</caption>

plate almost does not diffuse beyond the rectangular region located just on top of the heaters, reflecting the specific design of the heated flat plate.

## 5.5. Discussion

As shown in the previous sections, the numerical simulations are able to capture the correct trends and orders of magnitude with respect to the experimental results. A common aspect to all cases is that the numerical simulations predict a steady state temperature equal to the melting temperature. This is due to the fact that the accretion layer is assumed to be a uniform mixture. Therefore the temperature field in a given cell is uniform and no a posteriori reconstruction of a temperature profile is performed. In the steady state, the numerical simulation predicts that the accretion layer is in mixed phase ($0 < \alpha_l < 1$) and therefore that the temperature has to be equal to the melting temperature. However, the measured temperature is slightly higher than the melting temperature and in some case starts to increase after the temperature drop phase.

For example, in Fig. 20(d) the results from the experiments show a slight rise in the temperature once the ice layer is formed on the surface and steady state phase is reached. For the same ice water content, the heating power from the substrate is so high in this case that it results in formation of large amount of melt water at the interface (i.e. between the substrate and ice layer). This additional heat generates a temperature gradient in the water-ice layer. When ice starts accreting, the ice layer insulates the surface from incoming ice particles and cold airstream. Due to the constant heat flux from the surface all ice in contact with the surface melts and ultimately the surface temperature rises due to heat up of the liquid water present in the interface since the ice layer is no longer in direct contact with the substrate. Once this interfacial liquid water reaches a critical threshold it can also result in ice shedding as previously modeled by Kintea et al. [38] based on wedge icing tests performed by Currie et al. [39].

Hence, the model could be improved by going beyond the uniform temperature distribution in the simulation control volumes. Indeed, the resolution of the temperature field in the direction normal to the wall could help better predict the experimentally observed thermal dynamics.
<figure><img src="image_35.png" /><figcaption>(a) Experimental visualisation of the ice shape for Test 04 at $t = 100s$</figcaption></figure>

<figure><img src="image_36.png" /><figcaption>(b) Numerical prediction of the ice shape and temperature field for Test 04 at $t = 100s$</figcaption></figure>

<figure><img src="image_37.png" /><figcaption>(a) Experimental visualisation of the ice shape for Test 28 at $t = 100s$</figcaption></figure>

<figure><img src="image_38.png" /><figcaption>(b) Numerical prediction of the ice shape and temperature field for Test 28 at $t = 100s$</figcaption></figure>

<caption>**Fig. 22.** Qualitative comparison between experiment and numerical simulation of the ice shape for Test 28.</caption>

plate almost does not diffuse beyond the rectangular region located just on top of the heaters, reflecting the specific design of the heated flat plate.

## 5.5. Discussion

As shown in the previous sections, the numerical simulations are able to capture the correct trends and orders of magnitude with respect to the experimental results. A common aspect to all cases is that the numerical simulations predict a steady state temperature equal to the melting temperature. This is due to the fact that the accretion layer is assumed to be a uniform mixture. Therefore the temperature field in a given cell is uniform and no a posteriori reconstruction of a temperature profile is performed. In the steady state, the numerical simulation predicts that the accretion layer is in mixed phase ($0 < \alpha_l < 1$) and therefore that the temperature has to be equal to the melting temperature. However, the measured temperature is slightly higher than the melting temperature and in some case starts to increase after the temperature drop phase.

For example, in Fig. 20(d) the results from the experiments show a slight rise in the temperature once the ice layer is formed on the surface and steady state phase is reached. For the same ice water content, the heating power from the substrate is so high in this case that it results in formation of large amount of melt water at the interface (i.e. between the substrate and ice layer). This additional heat generates a temperature gradient in the water-ice layer. When ice starts accreting, the ice layer insulates the surface from incoming ice particles and cold airstream. Due to the constant heat flux from the surface all ice in contact with the surface melts and ultimately the surface temperature rises due to heat up of the liquid water present in the interface since the ice layer is no longer in direct contact with the substrate. Once this interfacial liquid water reaches a critical threshold it can also result in ice shedding as previously modeled by Kintea et al. [38] based on wedge icing tests performed by Currie et al. [39].

Hence, the model could be improved by going beyond the uniform temperature distribution in the simulation control volumes. Indeed, the resolution of the temperature field in the direction normal to the wall could help better predict the experimentally observed thermal dynamics.
In a similar way, a non-uniform liquid volume fraction distribu-
tion could help account for phenomena such as imbibition and also
help improve the sticking efficiency model. The model could also
be improved regarding the treatment of mass transport under the
effect of aerodynamic and capillary forces. Indeed, the model used
in this work assumes that in the presence of ice crystals in the ac-
cretion the liquid water is entirely trapped in the pores and cannot
be transported by aerodynamic or capillary forces. This assumption
is not correct for high liquid fractions and the model may be im-
proved in the future to account for this phenomenon.

Finally, when the liquid fraction is very high, the water/ice mix-
ture behaves as a very viscous non-Newtonian fluid that can be
carried downstream by the airflow. The model can also be adapted
in the future to take into account this physical mechanism.

6. Conclusion

This research work is an effort towards better understanding of
origins, necessary conditions for ice accretion, accretion inception
and icing physics. To that end, dedicated transparent and metallic
test articles were designed and tested in the icing wind tunnel of
TU Braunschweig. Foremost, a dedicated test campaign was con-
ducted utilizing a transparent ITO coated sapphire glass flat plate
model with heating capability for microscopic investigation of ice
accretion inception phase. High speed cameras and high powered
flicker-free light sources were used to capture the different stages
of ice accretion inception from interface and frontal perspective.
It was observed that ice crystals upon impact fragment and un-
dergo sticking and melting to form sessile droplets which with in-
take of incoming crystals from the icing cloud result in formation
of semi-frozen particles. Finally percolation and bridging of semi-
frozen particles leads to formation of ice layer. The semi-frozen
particle in its constitution is comprised of both liquid and solid
phase and served as a key assumption for the development of nu-
merical model presented in detail as part of this research work.

Furthermore, a dedicated test matrix was designed to study
the necessary conditions for origins of ice accretion. To have an
in depth overview of origins of ice accretion, a heatable flat plate
model was utilized whereby the heating power was varied as well
as the instance at which heaters were turned off. It was found that
for T<sub>wb</sub> < 0 with unheated substrate ice accretion did not initiate
because no liquid water was formed on the substrate for sticking
of ice particles on the surface. Substrate heating was found to be
necessary for ice accretion on the substrate leading to formation
of glaciated ice layer. On the other hand, for T<sub>wb</sub> > 0 heated and
unheated substrate both lead to ice accretion due to natural melt-
ing of ice crystals. Substrate heating was not necessary whereas
natural melting of ice crystals was the driving force for ice accre-
tion initiation leading to formation of slushy ice layer. In addition
to that parametric study was performed both experimentally and
numerically to study the influence of dominant parameters on ice
accretion phenomenon. It was interesting to observe that increas-
ing the flow velocity lead to a drop in ice accretion growth due
to influence of erosion. The time it took for the ice to accrete on
the surface drastically reduced upon increasing the flow velocity,
ice water content and wet bulb temperature. Increasing the heat-
ing power lead to longer duration required for the ice to accrete
on the surface as the icing cloud had to overcome a larger tem-
perature gradient. It was observed that the icing severity increased
with warmer wet bulb conditions. Increasing ice water content also
resulted in less time required for the ice to accrete on the substrate
as the larger mass flux of ice was able to overcome the large tem-
perature gradient in shorter amount of time.

In terms of numerical simulation and modelling, an extended
Messinger approach based on the enthalpy method was proposed.
This framework enables to capture the unsteady feature of the ex-

periments. It is also more representative of the observed accretion
process than an unsteady triple-layer approach because it allows
to take naturally into account the fact that the accretion is consti-
tuted by a homogeneous mixture of crystals and liquid water and
not by a superposition of ice layer and liquid water film. Within
this framework, the experimental database enabled to extend the
existing sticking efficiency model to the case of heated substrates.
The model was calibrated using three cases so as to take into ac-
count the influence of ice crystal velocity. The other cases of the
database were used to assess the predictive performance of the
model. The numerical simulations showed that the model was able
to capture the correct trends and orders of magnitude. Therefore
the proposed model provides an interesting way forward for ice
crystal icing simulations. On the other hand, some improvements
to account for additional physical mechanisms were discussed and
are part of ongoing research. To summarize in a nutshell, the pre-
sented work helps further the understanding of ice crystal icing
physics for heated and non-heated scenarios and contributes to
advancement of numerical tools towards improved predictive ca-
pability.

Declaration of Competing Interest

The authors declare the following financial interests/personal
relationships which may be considered as potential competing in-
terests:

Yasir A. Malik, Lokman Bennani, Stephan Bansmer, Pierre Tron-
tin, Philippe Villedieu reports financial support was provided by
European Commission.

CRediT authorship contribution statement

**Yasir A. Malik:** Conceptualization, Methodology, Software, Val-
idation, Formal analysis, Investigation, Resources, Data curation,
Writing – original draft, Writing – review & editing, Visualiza-
tion. **Lokman Bennani:** Methodology, Software, Validation, Formal
analysis, Investigation, Resources, Writing – original draft, Writ-
ing – review & editing. **Stephan Bansmer:** Conceptualization, Su-
pervision, Writing – review & editing. **Pierre Trontin:** Methodol-
ogy, Software, Validation, Formal analysis, Investigation, Resources,
Writing – original draft, Writing – review & editing. **Philippe
Villedieu:** Conceptualization, Methodology, Writing – review &
editing, Project administration, Funding acquisition.

Data availability

Data will be made available on request.

Acknowledgments

This research work is carried out and received funding under
the framework of European Union's Horizon 2020 research and in-
novation program with grant agreement No. 767560. The authors
would like to express gratitude to Stephan Sattler and Reinhard
Kerbstadt for sharing valuable input and expertise on model as-
sembly and electronics. Special thanks to Alex Vorgias who as-
sisted as part of experimental setup and campaign runs. A special
note of gratitude to Elmar Bonaccurso from Airbus CRT for lending
the High speed Camera (Photron SAZ Cam 2100K) utilized as part
of the experimental setup.

Appendix A. Heat transfer coefficient, collection efficiency and heating power calibration

In order to simulate ice accretion on the heated flat plate, one
must first compute the heat transfer coefficient and the ice crys-
tal collection efficiency. In this work the heat transfer coefficient
In a similar way, a non-uniform liquid volume fraction distribu-
tion could help account for phenomena such as imbibition and also
help improve the sticking efficiency model. The model could also
be improved regarding the treatment of mass transport under the
effect of aerodynamic and capillary forces. Indeed, the model used
in this work assumes that in the presence of ice crystals in the ac-
cretion the liquid water is entirely trapped in the pores and cannot
be transported by aerodynamic or capillary forces. This assumption
is not correct for high liquid fractions and the model may be im-
proved in the future to account for this phenomenon.

Finally, when the liquid fraction is very high, the water/ice mix-
ture behaves as a very viscous non-Newtonian fluid that can be
carried downstream by the airflow. The model can also be adapted
in the future to take into account this physical mechanism.

6. Conclusion

This research work is an effort towards better understanding of
origins, necessary conditions for ice accretion, accretion inception
and icing physics. To that end, dedicated transparent and metallic
test articles were designed and tested in the icing wind tunnel of
TU Braunschweig. Foremost, a dedicated test campaign was con-
ducted utilizing a transparent ITO coated sapphire glass flat plate
model with heating capability for microscopic investigation of ice
accretion inception phase. High speed cameras and high powered
flicker-free light sources were used to capture the different stages
of ice accretion inception from interface and frontal perspective.
It was observed that ice crystals upon impact fragment and un-
dergo sticking and melting to form sessile droplets which with in-
take of incoming crystals from the icing cloud result in formation
of semi-frozen particles. Finally percolation and bridging of semi-
frozen particles leads to formation of ice layer. The semi-frozen
particle in its constitution is comprised of both liquid and solid
phase and served as a key assumption for the development of nu-
merical model presented in detail as part of this research work.

Furthermore, a dedicated test matrix was designed to study
the necessary conditions for origins of ice accretion. To have an
in depth overview of origins of ice accretion, a heatable flat plate
model was utilized whereby the heating power was varied as well
as the instance at which heaters were turned off. It was found that
for T<sub>wb</sub> < 0 with unheated substrate ice accretion did not initiate
because no liquid water was formed on the substrate for sticking
of ice particles on the surface. Substrate heating was found to be
necessary for ice accretion on the substrate leading to formation
of glaciated ice layer. On the other hand, for T<sub>wb</sub> > 0 heated and
unheated substrate both lead to ice accretion due to natural melt-
ing of ice crystals. Substrate heating was not necessary whereas
natural melting of ice crystals was the driving force for ice accre-
tion initiation leading to formation of slushy ice layer. In addition
to that parametric study was performed both experimentally and
numerically to study the influence of dominant parameters on ice
accretion phenomenon. It was interesting to observe that increas-
ing the flow velocity lead to a drop in ice accretion growth due
to influence of erosion. The time it took for the ice to accrete on
the surface drastically reduced upon increasing the flow velocity,
ice water content and wet bulb temperature. Increasing the heat-
ing power lead to longer duration required for the ice to accrete
on the surface as the icing cloud had to overcome a larger tem-
perature gradient. It was observed that the icing severity increased
with warmer wet bulb conditions. Increasing ice water content also
resulted in less time required for the ice to accrete on the substrate
as the larger mass flux of ice was able to overcome the large tem-
perature gradient in shorter amount of time.

In terms of numerical simulation and modelling, an extended
Messinger approach based on the enthalpy method was proposed.
This framework enables to capture the unsteady feature of the ex-

periments. It is also more representative of the observed accretion
process than an unsteady triple-layer approach because it allows
to take naturally into account the fact that the accretion is consti-
tuted by a homogeneous mixture of crystals and liquid water and
not by a superposition of ice layer and liquid water film. Within
this framework, the experimental database enabled to extend the
existing sticking efficiency model to the case of heated substrates.
The model was calibrated using three cases so as to take into ac-
count the influence of ice crystal velocity. The other cases of the
database were used to assess the predictive performance of the
model. The numerical simulations showed that the model was able
to capture the correct trends and orders of magnitude. Therefore
the proposed model provides an interesting way forward for ice
crystal icing simulations. On the other hand, some improvements
to account for additional physical mechanisms were discussed and
are part of ongoing research. To summarize in a nutshell, the pre-
sented work helps further the understanding of ice crystal icing
physics for heated and non-heated scenarios and contributes to
advancement of numerical tools towards improved predictive ca-
pability.

Declaration of Competing Interest

The authors declare the following financial interests/personal
relationships which may be considered as potential competing in-
terests:

Yasir A. Malik, Lokman Bennani, Stephan Bansmer, Pierre Tron-
tin, Philippe Villedieu reports financial support was provided by
European Commission.

CRediT authorship contribution statement

**Yasir A. Malik:** Conceptualization, Methodology, Software, Val-
idation, Formal analysis, Investigation, Resources, Data curation,
Writing – original draft, Writing – review & editing, Visualiza-
tion. **Lokman Bennani:** Methodology, Software, Validation, Formal
analysis, Investigation, Resources, Writing – original draft, Writ-
ing – review & editing. **Stephan Bansmer:** Conceptualization, Su-
pervision, Writing – review & editing. **Pierre Trontin:** Methodol-
ogy, Software, Validation, Formal analysis, Investigation, Resources,
Writing – original draft, Writing – review & editing. **Philippe
Villedieu:** Conceptualization, Methodology, Writing – review &
editing, Project administration, Funding acquisition.

Data availability

Data will be made available on request.

Acknowledgments

This research work is carried out and received funding under
the framework of European Union's Horizon 2020 research and in-
novation program with grant agreement No. 767560. The authors
would like to express gratitude to Stephan Sattler and Reinhard
Kerbstadt for sharing valuable input and expertise on model as-
sembly and electronics. Special thanks to Alex Vorgias who as-
sisted as part of experimental setup and campaign runs. A special
note of gratitude to Elmar Bonaccurso from Airbus CRT for lending
the High speed Camera (Photron SAZ Cam 2100K) utilized as part
of the experimental setup.

Appendix A. Heat transfer coefficient, collection efficiency and heating power calibration

In order to simulate ice accretion on the heated flat plate, one
must first compute the heat transfer coefficient and the ice crys-
tal collection efficiency. In this work the heat transfer coefficient
<figure><img src="image_38.png" /></figure>

<figure><img src="image_39.png" /><figcaption>Fig. A.1. Heat transfer coefficients computed with the boundary layer solver CLICET and collection efficiency computed with the lagrangian trajectory solver TrajL2D (for tests 07, 08 and 04). s is the curvilinear abscissa who's reference point is the middle of the interface Γ.</figcaption></figure>

was computed using the Prandtl boundary layer solver CLICET. Fig. A.23a shows the heat transfer coefficient for tests 07, 08 and 04. The ice crystal trajectories were computed using a lagrangian solver (TrajL2D). The computed collection efficiency is shown in Fig. A.23b for tests 07, 08 and 04.

In addition, for every simulation, the heating power was reduced by 10% to account for possible heat losses due to three dimensional effects which are not captured by the present two dimensional framework.

## Appendix B. Mesh convergence analysis

In all simulations presented in Section 5, the mesh used to solve the heat equation inside the flat plate was composed of 25000 triangular elements. This mesh is sufficiently refined to provide converged results. Indeed, Figs. B.24 and B.25 show a comparison between results obtained with the base mesh (25000 elements) and a finer mesh (90000 elements). With the exception of very slight differences at the stagnation point for the ice shapes, the results are identical.

<figure><img src="image_40.png" /><figcaption>Fig. B.1. Comparison of temperature drop obtained with the base mesh and the refined mesh.</figcaption></figure>
<figure><img src="image_39.png" /><figcaption>Fig. A.1. Heat transfer coefficients computed with the boundary layer solver CLICET and collection efficiency computed with the lagrangian trajectory solver TrajL2D (for tests 07, 08 and 04). s is the curvilinear abscissa who's reference point is the middle of the interface Γ.</figcaption></figure>

was computed using the Prandtl boundary layer solver CLICET. Fig. A.23a shows the heat transfer coefficient for tests 07, 08 and 04. The ice crystal trajectories were computed using a lagrangian solver (TrajL2D). The computed collection efficiency is shown in Fig. A.23b for tests 07, 08 and 04.

In addition, for every simulation, the heating power was reduced by 10% to account for possible heat losses due to three dimensional effects which are not captured by the present two dimensional framework.

## Appendix B. Mesh convergence analysis

In all simulations presented in Section 5, the mesh used to solve the heat equation inside the flat plate was composed of 25000 triangular elements. This mesh is sufficiently refined to provide converged results. Indeed, Figs. B.24 and B.25 show a comparison between results obtained with the base mesh (25000 elements) and a finer mesh (90000 elements). With the exception of very slight differences at the stagnation point for the ice shapes, the results are identical.

<figure><img src="image_40.png" /><figcaption>Fig. B.1. Comparison of temperature drop obtained with the base mesh and the refined mesh.</figcaption></figure>
<figure><img src="image_41.png" /><figcaption>Fig. B.2. Comparison of ice shapes predicted at $t = 100$s with the base mesh and the refined mesh.</figcaption></figure>

References

[1] J. Mason, Engine power loss in ice crystal conditions (2007).

[2] A. Baumert, Experimental and Numerical Studies on Ice Crystal Icing of Civil Aircraft, TU Braunschweig-Niedersächsisches Forschungszentrum für Luftfahrt, 2019.

[3] V. Charton, Modélisation de l'accrétion de glace dans les turboréacteurs en conditions cristaux, Université de Toulouse, 2020 Ph.D. thesis.

[4] Bureau d'Enquêtes et d'Analyses pour la sécurité de l'avion civile, Final report on the accident on 1st June 2009 to the airbus a330-203 registered f-gzcp operated by air france flight af 447 rio de janeiro-paris, 2012,

[5] J. Mason, W. Strapp, P. Chow, The ice particle threat to engines in flight, in: 44th AIAA aerospace sciences meeting and exhibit, 2006, p. 206.

[6] Community research and development information service. haic project information., (https://cordis.europa.eu/project/id/314314).

[7] B. Messinger, Equilibrium temperature of an unheated icing surface as a function of air speed, J. Aerospace Sci. 20 (1953) 29.

[8] W.B. Wright, R. Gent, D. Guffond, DRA/NASA/ONERA collaboration on icing research, Technical Report, 1997.

[9] E. Radenac, Validation of a 3D ice accretion tool on swept wings of the SUNSET2 program, in: 8th AIAA Atmospheric and Space Environments Conference, 2016, p. 3735.

[10] W.B. Wright, C.E. Porter, A revised validation process for ice accretion codes, in: 9th AIAA Atmospheric and Space Environments Conference, 2017, p. 3415.

[11] P. Villedieu, P. Trontin, R. Chauvin, Glaciated and mixed phase ice accretion modeling using ONERA 2D icing suite, in: 6th AIAA Atmospheric and Space Environments Conference, AIAA-2014-2199, Atlanta, USA, 2014.

[12] P. Trontin, G. Blanchard, P. Villedieu, A comprehensive numerical model for mixed-phase and glaciated icing conditions, in: 8th AIAA Atmospheric and Space Environments Conference, AIAA-2016-3742, 2016, p. 3742.

[13] P. Trontin, P. Villedieu, A comprehensive accretion model for glaciated icing conditions, Int. J. Multiphase Flow 108 (2018) 105-123.

[14] H. Beaugendre, F. Morency, W.G. Habashi, Fensap-ice's three-dimensional in-flight ice accretion module: ice3d, J. Aircr. 40 (2) (2003) 239-247.

[15] C.N. Aliaga, M.S. Aubé, G.S. Baruzzi, W.G. Habashi, Fensap-ice-unsteady: unified in-flight icing simulation methodology for aircraft, rotorcraft, and jet engines, J. Aircr. 48 (1) (2011) 119-126.

[16] S. Nilamdeen, V.S. Rao, D. Switchenko, J. Selvanayagam, I. Ozcer, G.S. Baruzzi, Numerical simulation of ice crystal accretion inside an engine core stator, Technical Report, SAE Technical Paper, 2019.

[17] Q. Yang, X. Guo, W. Dong, A. Wang, Ice accretion and aerodynamic effects on a turbofan engine nacelle under takeoff conditions, Aerosp. Sci. Technol. 126 (2022) 107571.

[18] C. Laurent, M. Bouyges, V. Charton, L. Bennani, J.-M. Senoner, Ice crystals accretion capabilities of oneras 3d icing suite, in: AIAA AVIATION 2022 Forum, 2022, p. 3537.

[19] T.G. Myers, Extension to the messinger model for aircraft icing, AIAA J. 39 (2) (2001) 211-218.

[20] T. Myers, S. Mitchell, G. Muchatibaya, M. Myers, A cubic heat balance integral method for one-dimensional melting of a finite thickness layer, Int. J. Heat Mass. Transf. 50 (25) (2007) 5305-5317.

[21] G. Gori, M. Zocca, M. Garabelli, A. Guardone, G. Quaranta, Polimice: a simulation framework for three-dimensional ice accretion, Appl. Math. Comput. 267 (2015) 96-107.

[22] R. Chauvin, L. Bennani, P. Trontin, P. Villedieu, An implicit time marching galerkin method for the simulation of icing phenomena with a triple layer model, Finite Elem. Anal. Des. 150 (2018) 20-33.

[23] L. Bennani, P. Trontin, R. Chauvin, P. Villedieu, A non-overlapping optimized schwarz method for the heat equation with non linear boundary conditions and with applications to de-icing, Comput. Math. Appl. 80 (6) (2020) 1500-1522.

[24] L. Bennani, P. Trontin, E. Radenac, Numerical simulation of an electrothermal ice protection system in anti-icing and deicing mode, Aerospace 10 (1) (2023) 75.

[25] A. Bucknell, M. McGilvray, D.R. Gillespie, G. Jones, B. Collier, A thermodynamic model for ice crystal accretion in aircraft engines: emm-c, Int. J. Heat Mass Transf. 174 (2021) 121270.

[26] V. Alexiades, A. Solomon, Mathematical Modeling of Melting and Freezing Processes, Hemisphere, 1993.

[27] S.E. Bansmer, A. Baumert, S. Sattler, I. Knop, D. Leroy, A. Schwarzenboeck, T. Jurkat-Witschas, C. Voigt, H. Pervier, B. Esposito, Design, construction and commissioning of the braunschweig icing wind tunnel, Atmos. Meas. Tech. 11 (6) (2018) 3221-3249.

[28] M. Breiling, M. Bacher, S. Sokratov, F.G. Best, Method and Device for Producing Snow, 2016, US Patent 9,429,348.

[29] S.E. Bansmer, A. Baumert, Ice crystal generation and conveyance system documentation, technical report, 2015.

[30] Y.A. Malik, L. Bennani, A. Vorgias, P. Trontin, P. Villedieu, Experimental and numerical investigation of ice crystal icing on a heatable naca0012 airfoil, in: AIAA AVIATION 2022 Forum, 2022, p. 3534.

[31] A. Razak, Industrial Gas Turbines: Performance and Operability, Elsevier, 2007.

[32] R. Stull, Wet-bulb temperature from relative humidity and air temperature, J. Appl. Meteorol. Climatol. 50 (11) (2011) 2267-2269.

[33] J. Löwe, D. Kintea, A. Baumert, S. Bansmer, I.V. Roisman, C. Tropea, Inception of ice accretion by ice crystal impact, in: Journal of Physics: Conference Series, volume 745, IOP Publishing, 2016, p. 032013.

[34] P. Trontin, A. Kontogiannis, G. Blanchard, P. Villedieu, Description and assessment of the new onera 2d icing suite igloo2d, in: 9th AIAA Atmospheric and Space Environments Conference, AIAA-2017-3417, Denver, USA, 2017.

[35] C. Geuzaine, J.-F. Remacle, Gmsh: a 3-d finite element mesh generator with built-in pre-and post-processing facilities, Int. J. Numer. Methods Eng. 79 (11) (2009) 1309-1331.

[36] A. Baumert, S. Bansmer, P. Trontin, P. Villedieu, Experimental and numerical investigations on aircraft icing at mixed phase conditions, Int. J. Heat Mass. Transf. 123 (2018) 957-978.

[37] L. Reitter, H. Lohmann, M. Schremb, I. Roisman, J. Hussong, C. Tropea, Impact of an ice particle onto a dry rigid substrate: dynamic sintering of a residual ice cone, Cold Reg. Sci. Technol. 194 (2022) 103416.

[38] D. Kintea, I. Roisman, C. Tropea, Transport processes in a wet granular ice layer: model for ice accretion and shedding, Int. J. Heat Mass Transf. 97 (2016) 461-472.

[39] T. Currie, P. Struk, J.-C. Tsao, D. Fuleki, D. Knezevici, Fundamental study of mixed-phase icing with application to ice crystal accretion in aircraft jet engines, in: 4th AIAA atmospheric and space environments conference, 2012, p. 3035.
<figure><img src="image_41.png" /><figcaption>Fig. B.2. Comparison of ice shapes predicted at $t = 100s$ with the base mesh and the refined mesh.</figcaption></figure>

References

[1] J. Mason, Engine power loss in ice crystal conditions (2007).

[2] A. Baumert, Experimental and Numerical Studies on Ice Crystal Icing of Civil Aircraft, TU Braunschweig-Niedersächsisches Forschungszentrum für Luftfahrt, 2019.

[3] V. Charton, Modélisation de l'accrétion de glace dans les turboréacteurs en conditions cristaux, Université de Toulouse, 2020 Ph.D. thesis.

[4] Bureau d'Enquêtes et d'Analyses pour la sécurité de l'avion civile, Final report on the accident on 1st june 2009 to the airbus a330-203 registered f-gzcp operated by air france flight af 447 rio de janeiro-paris, 2012,

[5] J. Mason, W. Strapp, P. Chow, The ice particle threat to engines in flight, in: 44th AIAA aerospace sciences meeting and exhibit, 2006, p. 206.

[6] Community research and development information service. haic project information., (https://cordis.europa.eu/project/id/314314).

[7] B. Messinger, Equilibrium temperature of an unheated icing surface as a function of air speed, J. Aerospace Sci. 20 (1953) 29.

[8] W.B. Wright, R. Gent, D. Guffond, DRA/NASA/ONERA collaboration on icing research, Technical Report, 1997.

[9] E. Radenac, Validation of a 3D ice accretion tool on swept wings of the SUNSET2 program, in: 8th AIAA Atmospheric and Space Environments Conference, 2016, p. 3735.

[10] W.B. Wright, C.E. Porter, A revised validation process for ice accretion codes, in: 9th AIAA Atmospheric and Space Environments Conference, 2017, p. 3415.

[11] P. Villedieu, P. Trontin, R. Chauvin, Glaciated and mixed phase ice accretion modeling using ONERA 2D icing suite, in: 6th AIAA Atmospheric and Space Environments Conference, AIAA-2014-2199, Atlanta, USA, 2014.

[12] P. Trontin, G. Blanchard, P. Villedieu, A comprehensive numerical model for mixed-phase and glaciated icing conditions, in: 8th AIAA Atmospheric and Space Environments Conference, AIAA-2016-3742, 2016, p. 3742.

[13] P. Trontin, P. Villedieu, A comprehensive accretion model for glaciated icing conditions, Int. J. Multiphase Flow 108 (2018) 105-123.

[14] H. Beaugendre, F. Morency, W.G. Habashi, Fensap-ice's three-dimensional in-flight ice accretion module: ice3d, J. Aircr. 40 (2) (2003) 239-247.

[15] C.N. Aliaga, M.S. Aubé, G.S. Baruzzi, W.G. Habashi, Fensap-ice-unsteady: unified in-flight icing simulation methodology for aircraft, rotorcraft, and jet engines, J. Aircr. 48 (1) (2011) 119-126.

[16] S. Nilamdeen, V.S. Rao, D. Switchenko, J. Selvanayagam, I. Ozcer, G.S. Baruzzi, Numerical simulation of ice crystal accretion inside an engine core stator, Technical Report, SAE Technical Paper, 2019.

[17] Q. Yang, X. Guo, W. Dong, A. Wang, Ice accretion and aerodynamic effects on a turbofan engine nacelle under takeoff conditions, Aerosp. Sci. Technol. 126 (2022) 107571.

[18] C. Laurent, M. Bouyges, V. Charton, L. Bennani, J.-M. Senoner, Ice crystals accretion capabilities of oneras 3d icing suite, in: AIAA AVIATION 2022 Forum, 2022, p. 3537.

[19] T.G. Myers, Extension to the messinger model for aircraft icing, AIAA J. 39 (2) (2001) 211-218.

[20] T. Myers, S. Mitchell, G. Muchatibaya, M. Myers, A cubic heat balance integral method for one-dimensional melting of a finite thickness layer, Int. J. Heat Mass. Transf. 50 (25) (2007) 5305-5317.

[21] G. Gori, M. Zocca, M. Garabelli, A. Guardone, G. Quaranta, Polimice: a simulation framework for three-dimensional ice accretion, Appl. Math. Comput. 267 (2015) 96-107.

[22] R. Chauvin, L. Bennani, P. Trontin, P. Villedieu, An implicit time marching galerkin method for the simulation of icing phenomena with a triple layer model, Finite Elem. Anal. Des. 150 (2018) 20-33.

[23] L. Bennani, P. Trontin, R. Chauvin, P. Villedieu, A non-overlapping optimized schwarz method for the heat equation with non linear boundary conditions and with applications to de-icing, Comput. Math. Appl. 80 (6) (2020) 1500-1522.

[24] L. Bennani, P. Trontin, E. Radenac, Numerical simulation of an electrothermal ice protection system in anti-icing and deicing mode, Aerospace 10 (1) (2023) 75.

[25] A. Bucknell, M. McGilvray, D.R. Gillespie, G. Jones, B. Collier, A thermodynamic model for ice crystal accretion in aircraft engines: emm-c, Int. J. Heat Mass Transf. 174 (2021) 121270.

[26] V. Alexiades, A. Solomon, Mathematical Modeling of Melting and Freezing Processes, Hemisphere, 1993.

[27] S.E. Bansmer, A. Baumert, S. Sattler, I. Knop, D. Leroy, A. Schwarzenboeck, T. Jurkat-Witschas, C. Voigt, H. Pervier, B. Esposito, Design, construction and commissioning of the braunschweig icing wind tunnel, Atmos. Meas. Tech. 11 (6) (2018) 3221-3249.

[28] M. Breiling, M. Bacher, S. Sokratov, F.G. Best, Method and Device for Producing Snow, 2016, US Patent 9,429,348.

[29] S.E. Bansmer, A. Baumert, Ice crystal generation and conveyance system documentation, technical report, 2015.

[30] Y.A. Malik, L. Bennani, A. Vorgias, P. Trontin, P. Villedieu, Experimental and numerical investigation of ice crystal icing on a heatable naca0012 airfoil, in: AIAA AVIATION 2022 Forum, 2022, p. 3534.

[31] A. Razak, Industrial Gas Turbines: Performance and Operability, Elsevier, 2007.

[32] R. Stull, Wet-bulb temperature from relative humidity and air temperature, J. Appl. Meteorol. Climatol. 50 (11) (2011) 2267-2269.

[33] J. Löwe, D. Kintea, A. Baumert, S. Bansmer, I.V. Roisman, C. Tropea, Inception of ice accretion by ice crystal impact, in: Journal of Physics: Conference Series, volume 745, IOP Publishing, 2016, p. 032013.

[34] P. Trontin, A. Kontogiannis, G. Blanchard, P. Villedieu, Description and assessment of the new onera 2d icing suite igloo2d, in: 9th AIAA Atmospheric and Space Environments Conference, AIAA-2017-3417, Denver, USA, 2017.

[35] C. Geuzaine, J.-F. Remacle, Gmsh: a 3-d finite element mesh generator with built-in pre-and post-processing facilities, Int. J. Numer. Methods Eng. 79 (11) (2009) 1309-1331.

[36] A. Baumert, S. Bansmer, P. Trontin, P. Villedieu, Experimental and numerical investigations on aircraft icing at mixed phase conditions, Int. J. Heat Mass. Transf. 123 (2018) 957-978.

[37] L. Reitter, H. Lohmann, M. Schremb, I. Roisman, J. Hussong, C. Tropea, Impact of an ice particle onto a dry rigid substrate: dynamic sintering of a residual ice cone, Cold Reg. Sci. Technol. 194 (2022) 103416.

[38] D. Kintea, I. Roisman, C. Tropea, Transport processes in a wet granular ice layer: model for ice accretion and shedding, Int. J. Heat Mass Transf. 97 (2016) 461-472.

[39] T. Currie, P. Struk, J.-C. Tsao, D. Fuleki, D. Knezevici, Fundamental study of mixed-phase icing with application to ice crystal accretion in aircraft jet engines, in: 4th AIAA atmospheric and space environments conference, 2012, p. 3035.