Analysis of Experimental Ice Accretion Data and Assessment of a Thermodynamic Model During Ice Crystal Icing

SAE 2019-01-2016

Analysis of Experimental Ice Accretion Data and Assessment of a Thermodynamic

Model During Ice Crystal Icing

Tadas P. Bartkus and Jen-Ching Tsao

Ohio Aerospace Institute

Peter M. Struk

NASA John Glenn Research Center

Abstract Introduction

This paper analyzes ice crystal icing accretion data and evaluates a Since the 1990s, there have been numerous reports of turbofan engine thermodynamic ice crystal icing model, which has been previously power-loss or damage events that have been attributed to the ingestion presented, to describe the possible mechanisms of icing within the core of ice crystals. These events have been observed at altitudes at or above of a turbofan jet engine. The model functions between two distinct ice the upper limit at which water droplets can naturally exist as liquid.

accretions based on a surface energy balance: freeze-dominated icing These events typically have occurred in deep convective updraft and melt-dominated icing. Freeze-dominated icing occurs when liquid systems and have included engine stall, rollback, flameout, and engine water (from melted ice crystals) freezes and accretes on a surface along component damage. Mason et al. [1] hypothesized that ice crystals with the existing ice of the impinging water and ice mass. This ingested into the engine begin to undergo partial melting within the freeze-dominated icing is characterized as having strong adhesion to compressor system and then, as a mixed-phase water mass, accrete on the surface. The amount of ice accretion is partially dictated by a freeze surfaces within the engine core.

fraction, which is the fraction of impinging liquid water that freezes.

Melt-dominated icing occurs as unmelted ice on a surface accumulates. This threat of engine icing has spawned a substantial research effort in This melt-dominated icing is characterized by weakly bonded surface understanding the fundamental physics of ice crystal icing. To better adhesion. The amount of ice accumulation is partially dictated by a understand the physical mechanisms of ice crystal icing, experiments melt fraction, which is the fraction of impinging ice crystals that melts. have been conducted at the National Aeronautics and Space Experimentally observed ice growth rates suggest that only a small Administration (NASA) Glenn Research Center. Experimental studies fraction of the impinging ice remains on the surface, implying a mass on the fundamentals of ice crystal icing physics were conducted at the loss mechanism such as splash, runback, bounce, or erosion. The NASA Propulsion Systems Laboratory (PSL) in 2016 [2] and 2018 [3].

fraction of mass loss must be determined in conjunction with the Other fundamental physics studies conducted elsewhere [4-12] have fraction of freezing liquid water or fraction of melting ice on an icing investigated altitude scaling [8], mixed-phase sticking efficiency [9], surface for a given ice growth rate. This mass loss parameter, however, particle size effects [10, 11], and accretion angle effects [12] related to along with the freeze fraction and melt fraction, are the only ice crystal icing.

experimental parameters that are currently not measured directly.

Using icing growth rates from ice crystal icing experiments, a In 2013, Currie et al. [8] of the National Research Council of Canada methodology that has been previously proposed is used to determine conducted a series of ice crystal accretion tests in development of these unknown parameters. This work takes ice accretion data from scaling laws. Accretions were conducted on an axis-symmetric tests conducted by the National Aeronautics and Space Administration hemisphere (44.5 mm diameter) attached to a conical afterbody. Tests (NASA) at the Glenn Research Center in 2018 that examined the were conducted at air speeds of Mach = 0.25, at absolute pressures of fundamental physics of ice crystal icing. This paper continues 34.5 kPa and 69 kPa, and warmer than freezing temperatures, where evaluation of the thermodynamic model from a previous effort, with the wet-bulb temperature ranged from approximately 0 ° C to 6 ° C. The additions to the model that account for sub-freezing temperatures that cloud median volumetric diameter, MVD , was measured to be 57 μm have been observed at the leading edge of the airfoil during icing. The [13]. Currie reported ice accretions that reached a steady-state size predicted temperatures were generally in good agreement with during a continuous exposure to a mixed-phase icing cloud. The measured temperatures. Other key findings include the total wet-bulb steady-state ice shape was hypothesized to be a balance between temperature being a good first order indicator of whether icing is accretion and erosion. Currie developed a semi-empirical model of the freeze-dominated (sub-freezing values) or melt-dominated (above ice accretions that introduced a concept called sticking efficiency. This freezing). Maximum sticking efficiency values, the fraction of sticking efficiency, 𝜂 ¸ was defined as the fraction of an impinging ௦௧௜௖௞ impinging mass that adheres to a surface, was calculated to be about mixed-phase water mass flux that was retained on the surface. The 0.2, and retained this maximum value for a range of melt ratios (0.3 to model treated the accretion process as a strictly physical sticking 0.65 and possibly higher), which is defined as the ratio of liquid water phenomenon, ignoring heat transfer, phase change, and runback.

content to total water content. Higher air velocities reduced the Currie began to quantify an icing severity chart that identifies regions maximum sticking efficiency and shifted the icing regime to higher of maximum sticking efficiency with respect to the impinging cloud melt ratio values. Finally, the leading edge ice accretion angle was melt ratio, 𝜂 . The melt ratio is defined as the ratio of liquid water ெோ found to be related to ice growth (lower growth rates for smaller content to the total water content of the impinging cloud. For the angles) and melt ratio (smaller melt ratios resulted in smaller angles, conditions run, Currie reported maximum sticking efficiency in the likely due to erosion effects).

melt ratio range of 0.10 – 0.25 where the sticking efficiency ranged from 0.27 to 0.40 at the stagnation point, depending on the total water Page 1 of 17 content of the impinging cloud. Sticking efficiency dropped off considered in work by Currie et al. [8, 9]. The work by Tsao stated that there are two distinct types of ice accretions based on an overall energy significantly at 𝜂 values smaller and greater than the mentioned ெோ range, creating an icing severity plateau. balance at the accretion site. The first type is freeze-dominated icing, which occurs when liquid in a mixed-phase cloud freezes onto a In 2014, Currie et al. [9] expanded the research efforts of ice crystal surface along with ice present in the cloud. This freeze-dominated icing by running ice accretion test on three different airfoil geometries. icing is characterized as having strong adhesion to the surface. The The airfoils included the aforementioned axis-symmetric hemisphere, amount of ice accretion is partially dictated by a freeze fraction, n , a crowned cylinder, and an axis-symmetric cone, where all three which is the fraction of impinging liquid water that freezes.

forebodies had a 44.5 mm diameter. Tests were conducted at flow Melt-dominated icing occurs as unmelted ice on a surface accumulates.

This melt-dominated icing is characterized by weakly bonded surface speeds of Mach 0.25 and 0.40 and absolute pressures of 34.5 kPa and 69 kPa, with cloud MVD of 57 μm [13]. Variations existed between adhesion. The amount of ice accumulation is partially dictated by a different airfoil geometries, however, sticking efficiency was greatest melt fraction, m , which is the fraction of impinging ice crystals that melts. To determine the mass fraction values of n or m , a mass loss in the melt ratio range of 0.1 – 0.2 and reached stagnation point 𝜂 0 0 ௦௧௜௖௞ values around 0.4 - 0.5. Abrupt decreases in 𝜂 value were again parameter, n loss , must be determined. This mass loss fraction, along ௦௧௜௖௞ with the freeze fraction or melt fraction, are the only parameters that reported outside of the 0.1 - 0.2 melt ratio range, creating an icing are currently not measured directly.

severity plateau. Currie reported that increased Mach number decreased the icing severity map, by creating a narrower icing severity Using reported icing growth rates from published ice crystal icing plateau. Some of Currie’s data suggests that at the higher Mach experiments, a methodology was proposed by Bartkus et al. [17] to number, the height of the plateau is shallower too (possibly, according to Currie). Sticking efficiency decreased at oblique particle determine these unknown parameters. The paper built on the previously proposed model by Tsao et al. [16] by adding a transient impingement angles (away from the stagnation point). Currie also conduction term between the airfoil and ice accretion to explain ice reported that a steady-state ice shape was reached for all tests except for the highest total water content and slowest Mach number cases, growth behavior at the onset of experimental tests that was observed to be different from steady-state ice growth that occurred later in the where ice growth did not reach a limit. This is due to the sticking efficiency remaining non-zero at all oblique impingement angles. No test run. Several key findings were reported in the paper by Bartkus et significant effects due to pressure or total water content on sticking al. The work suggested that mass loss fractions can exceed n loss = 0.90 for steady ice growth periods. The ice accretion rate and mass loss efficiency were found.

value for this steady ice growth period were measured and calculated for times between 120 and 180 s of mixed-phase ice cloud exposure on The French National Aerospace Center, ONERA, has made significant a NACA 0012 airfoil. Lower mass loss fraction values were calculated efforts to better understand ice crystal icing. They have developed a during the initial transient period. Due to the additional melt from numerical icing tool, IGLOO2D, to address concerns raised by mixed conduction when using an initially warmer-than-freezing airfoil, a wet, phase icing. [12, 14]. The accretion model, based on the extension of sticky surface was likely the physical mechanism that allowed for more the Messinger model [15] that accounts for the presence of ice crystals of the incoming cloud to be captured, reducing n loss . All initial transient among the impinging cloud particles, was evaluated against results periods in the experiments that were evaluated were determined to be from comprehensive icing wind tunnel tests conducted by the French melt-dominated icing, where transient measurements and calculations agency. Accretion tests were conducted by Baumert et al. [12] utilizing were taken within the first 20 s of mixed-phase ice cloud exposure. In a NACA 0012 airfoil and a cylindrical airfoil. The cylindrical airfoil addition, the paper noted that the total wet bulb temperature, Twb , had a diameter of 60 mm, which corresponded to the maximum which is the balance between convective and evaporative energy thickness of the NACA 0012 airfoil (500 mm chord length). Tests were fluxes, was a good indicator for determining freeze-dominated or conducted with an air velocity of 40 m/s (Mach = 0.12), and cloud melt-dominated icing in which conduction was negligible. Generally, MVD of 80 μm. Baumert was able to independently control air when Twb is below freezing, freeze-dominated icing exists, and when temperature and melt ratio. The paper reported several key findings.

Twb 0 is above freezing, melt-dominated icing exists. Another finding With respect to sticking efficiency and the icing severity map, Baumert stated that for conditions that were sufficiently cold and freeze fraction reported wide icing severity plateaus for both airfoils, where 𝜂 ெோ values were n 0 = 1, then the icing surface temperature at the stagnation ranged from about 0.2 to 0.6, with 𝜂 values between 0.3 and 0.4.

௦௧௜௖௞ point could continue to decrease below 0 °C until it reached a Accretion tests with the cylindrical airfoil extended the range to thermodynamic equilibrium temperature. This ice temperature for 𝜂 = 1, at wet-bulb temperatures of -5 ° C and -15 ° C, where icing ெோ these sufficiently cold cases was suggested to be the total wet-bulb transitioned to supercooled liquid ice accretion. As reported by temperature, but the paper did not provide any derivations or Baumert, at mixed phase conditions below the freezing point, no right calculations to support this statement. A new addition to the model in end of the icing severity map exists. Baumert reported that 𝜂 ௦௧௜௖௞ this work addresses this icing surface temperature. Using ice growth decreases with decreasing wet-bulb temperature, but is independent of rates, airfoil thermocouple measurements and test conditions from a total water content. Baumert observed time dependent accretion rates, companion paper by Struk et al. [3], the thermodynamic model will be which can be attributed to the leading edge ice accretion angle, φ .

further assessed, providing values to the unknown parameters of the Investigating the accretion angle, Baumert reported lower 𝜂 values ௦௧௜௖௞ melt-dominated or freeze-dominated ice accretions. These at lower φ values. This can be attributed to a reduction in collection fundamental icing physics experiments were conducted at the NASA efficiency around the stagnation point with narrower angles, and not PSL icing wind tunnel in June 2018 using a NACA 0012 airfoil. These benefitting from re-impingement after particle bounce/break-up that a icing tests and this work are part of NASA’s Advanced Air Transport wider angle ice accretion would experience. Finally, Baumert reported Technology (AATT) Project roadmap to improve understanding of the reduced φ values at lower 𝜂 values, likely due to erosion effects of ெோ ice growth physics and expand engine aero-thermodynamic modeling the more glaciated cloud.

capability to predictively assess the onset and growth of ice in current and future engines during flight.

In 2014, Tsao et al. [16] proposed a thermodynamic model to describe the possible mechanisms for ice crystal icing on surfaces within the core of a jet engine. This thermodynamic model included factors not Page 2 of 17 Surface Energy Balance Equations

Objectives

The main objective of this work is to evaluate the ice crystal icing Tsao et al. [16] identified two distinct mechanisms for ice crystal icing thermodynamic model proposed by Tsao et al. [16] by solving for the growth: freeze-dominated icing and melt-dominated icing. The conservation of energy expressions for the icing surface for both unknown parameters of n loss and n 0 or m 0 . These parameters will be determined by utilizing experimental ice accretion data from tests mechanisms are taken from Tsao [16] and Bartkus [17] and are conducted by NASA in 2018 [3]. The model will be evaluated utilizing reproduced in the following sections.

reported ice growth rates at the midspan of the NACA 0012 airfoil along with the measured conditions at the airfoil leading edge.

Freeze-Dominated Regime Additions were made to the model to determine the thermodynamic The governing energy conservation law at the icing surface for freeze- equilibrium icing surface temperature for tests where conditions were dominated icing is shown in Eq. (1). The sign convention assumes that sufficiently cold to produce freeze fraction values of n = 1. This new the net energy transfer from the right-hand side of the equation is addition will be evaluated utilizing temperature data measured by a positive. Therefore, the left-hand side of the expression is also positive thermocouple located at the midspan and leading edge of the airfoil.

and represents the energy that is available for freezing, and is the latent " " heat of fusion surface energy for freezing, 𝑞 . In Eq. (1), 𝑞 is ௙௥௘௘௭௘ ௘௩௔௣ A final objective is to investigate sticking efficiency for the NASA the evaporative heat transfer flux (heat transferred away from icing 2018 experimental ice accretion data, and how 𝜂 varies with the " ௦௧௜௖௞ surface = positive), 𝑞 is the convective heat transfer flux (heat to ௖௢௡௩ controlled parameters in the experiments, in particular, 𝜂 . In ெோ " icing surface = positive), 𝑞 is the kinetic energy transfer flux ௞௜௡௘௧௜௖ addition, this paper will examine the roll that the ice accretion angle, " (energy into icing surface = positive), and 𝑞 is the conductive heat ௖௢௡ௗ φ , played in the accretion process. These findings will be compared transfer flux (heat to icing surface = positive). These five surface with previous work by Currie [8, 9] and Baumert [12].

balance energy terms are described in greater detail in the next paragraphs.

Ice Crystal Icing Thermodynamic Model

Description

" " " " " 𝑞 = 𝑞 − 𝑞 − 𝑞 − 𝑞 ( 1 ) ௙௥௘௘௭௘ ௘௩௔௣ ௖௢௡௩ ௞௜௡௘௧௜௖ ௖௢௡ௗ A thermodynamic model for ice crystal icing within the core of jet engines was proposed by Tsao et al. [16]. Bartkus et al. [17] built on The evaporative heat transfer flux at the icing surface is shown in Eq.

" the previously proposed model by Tsao by adding a transient (2). In the expression, 𝑚 ̇ is the evaporative mass flux and 𝐿 is the ௘ ௩ conduction term to explain initial ice accretion growth rates that latent heat of vaporization. It should be noted that ice on the surface differed from ice accretion growth rates that occurred later in the ice will sublimate, and water will evaporate, so the value of 𝐿 will be ௩ cloud spray. Bartkus et al. [17] described the assumptions and surface dependent on the mixture (quality) of ice and liquid water. The energy balances that compose the model. They are repeated in the evaporative mass flux term can be expressed in terms of an evaporative following sections for thoroughness. In this current study, a surface mass transfer coefficient, ℎ , total temperature, 𝑇 , static temperature, ௠ ଴ temperature energy balance is added to the thermodynamic model to 𝑇 , total pressure, 𝑝 , static pressure, 𝑝 , saturation vapor pressure of ௦ ଴ ௦ address conditions where the icing surface temperature can become water in air, 𝑝 , and the saturation vapor pressure of water at the icing ௩ , ௦ sub-freezing in temperature. In addition to the surface temperature surface, 𝑝 .

௩ , ௦௨௥௙ energy balance, expressions for sticking efficiency, along with other mass loss factors are provided in the following sections.

೛ ೛ ೡ , ೞೠೝ೑ ೛ ೡ , ೞ బ ି ∙ ∙ ோு ೞ ೅ ೅ ೛ " " ೞ బ ೞ 𝑞 = 𝑚 ̇ 𝐿 = ℎ ቆ ቇ 𝐿 ( 2 ) ೛ ௘௩௔௣ ௘ ௩ ௠ ௩ భ ೛ ೡ , ೞೠೝ೑ బ Model Assumptions ∙ ି బ . లమమ ೅ ೅ బ ೞ The following assumptions are listed for assessing this thermodynamic The convective heat transfer flux at the icing surface is shown in Eq.

model: (3). In the expression, ℎ is the convective heat transfer coefficient, 𝑈 ௖ is the air velocity, 𝐶𝑝 is the specific heat capacity of air, and 𝑇 ௔௜௥ ௦௨௥௙  Steady icing cloud flow and air flow relative to the is the icing surface temperature. The temperature recovery factor, the accretion process.

recovery of energy from static temperature as flow decelerates in the  Accretion growth rates are taken at the stagnation point of boundary layer, can be approximated to be unity (1.0) at the leading an airfoil.

edge. For this reason, the sum of the first two terms within the  All water mass comes from the impinging water or ice parentheses of Eq. (3) is equivalent to the total air temperature.

cloud, (i.e. no water flow from neighboring surface control volumes).

మ ௎  All impinging water mass is at the freezing temperature of " 𝑞 = ℎ ቀ 𝑇 + − 𝑇 ቁ (3) ௖௢௡௩ ௖ ௦ ௦௨௥௙ ଶ஼௣ ೌ ೔ೝ water (0 °C).

 Coefficients of heat and mass transfer, as initially measured on a non-iced airfoil surface, remain constant despite The kinetic energy transfer flux at the icing surface is shown in Eq. (4).

changing geometries as ice accretes on the airfoil. " In the expression, 𝑚 ̇ is the impinging mass flux that sticks ௦௧௜௖௞ (accretes), and 𝑈 is the velocity of the impinging mass (which is approximated to be equal to the air velocity). The composition of " 𝑚 ̇ will be more fully described in the mass balance equation ௦௧௜௖௞ Page 3 of 17 " " " " section. It should be noted that in this form, the equation neglects any ௤ ି ௤ ି ௤ ି ௤ ೐ೡೌ೛ ೎೚೙ೡ ೖ೔೙೐೟೔೎ ೎೚೙೏ 𝑛 = (8) ଴ " kinetic energy transfer from the impinging mass flux that did not stick. ௠ ̇ ∙ 𝐿 ೔೘೛ , ೗೔೜ 𝑓 ∙ ൫ 1− ೙ ൯ ೗೚ೞೞ మ ௎ " " The individual energy fluxes in the numerator of Eq. (8) can be 𝑞 = 𝑚 ̇ ∙ (4) ௞௜௡௘௧௜௖ ௦௧௜௖௞ ଶ substituted with the respective expression from Eq. (2) through Eq. (5).

It should be noted that liquid mass must initially be present for freeze-dominated icing to occur. In addition, it should be noted that the The conduction heat transfer flux to the icing surface is shown in Eq.

kinetic energy transfer flux is a function of 𝑛 (it will be shown in (5). In the expression, 𝑘 is the thermal conductivity of the airfoil ௪௔௟௟ ௟௢௦௦ " surface, and the 𝑑𝑇 ⁄ 𝑑𝑛 ො term refers to the change in temperature in the mass balance section that 𝑚 ̇ is a function of 𝑛 ) and ௟௢௦௦ ௦௧௜௖௞ the normal direction within the airfoil at the ice and airfoil interface therefore Eq. (8) contains two unknown parameters.

( 𝑛 ො = 0 ). Heat transfer from the wall to the water/ice mix is positive in value, according the sign notation in Eq. (5).

Melt-Dominated Regime " ௗ் According to the sign convention in Eq. (1), the value of 𝑞 will " ௙௥௘௘௭௘ 𝑞 = 𝑘 ቃ (5) ௖௢௡ௗ ௪௔௟௟ ௗ௡ ො ௪௔௟௟ be positive in the freeze-dominated regime. If the sum of the terms on the right side of Eq. (1) is negative, then melt-dominated freezing " occurs. This is explicitly expressed in Eq. (9), where 𝑞 is the latent It should be noted that the conduction term as presented in Eq. (5) is in ௠௘௟௧ heat of fusion surface energy for melting (the energy available for steady-state form. Bartkus et al. [17] provided a derivation for a one melting). According to the sign convention in Eq. (9), the value of dimensional (1D) transient conduction term in which an airfoil is " initially at a different temperature than the surrounding icing cloud 𝑞 is positive for melt-dominated icing.

௠௘௟௧ environment. Approximating the stagnation point of the airfoil as a 1D plane of thickness 𝐿 , where the stagnation point temperature at 𝑛 ො = 0 " " " " " 𝑞 = 𝑞 + 𝑞 + 𝑞 − 𝑞 (9) ௠௘௟௧ ௖௢௡௩ ௞௜௡௘௧௜௖ ௖௢௡ௗ ௘௩௔௣ is fixed at the freezing temperature of 0 ºC, the internal wall at 𝑛 ො = 𝐿 as adiabatic, and the initial airfoil temperature as uniform, 𝑇 , the ௪௔௟௟ airfoil temperature at any time, 𝑡 , and any location within the airfoil, The individual energy fluxes on the right-hand side of Eq. (9) are 𝑛 ො , can be determined by the expression in Eq. (6). The infinite series identical as expressed in the previous Freeze-Dominated Regime " solution in Eq. (6) is solved by using Separation of Variables with section. The value of 𝑞 can be expressed in terms of the impinging ௠௘௟௧ Fourier Series methods.

" ice mass flux rate, 𝑚 ̇ , the latent heat of fusion, the ice mass melt ௜௠௣ , ௜௖௘ fraction, 𝑚 , and 𝑛 , and is shown in Eq. (10). The impinging ice ଴ ௟௢௦௦ మ మ ( మೕషభ ) ഏ ೖ mass flux can be expressed in terms of total water content, melt ratio, ೢೌ ೗೗ ష ቈ ቉ ೟ ഏ 𝑛 ො మ రಽ ഐ ಴೛ ೢೌ ೗೗ ೢೌ ೗೗ ௦௜௡ ቂ ( ଶ௝ିଵ ) ቃ ∙ ௘ collection efficiency at the stagnation line, and particle mass velocity.

ସ் ೢೌ ೗೗ , ೔ మಽ ஶ ∑ 𝑇 ( 𝑛 ො , 𝑡 ) = (6) ௝ୀଵ గ ଶ௝ିଵ " " ( ) 𝑞 = 𝑚 ̇ ∙ 1 − 𝑛 ∙ 𝑚 ∙ 𝐿 ௠௘௟௧ ௜௠௣ , ௜௖௘ ௟௢௦௦ ଴ ௙ In Eq. (6), 𝑗 represents the positive integers in the summation of the ( ) ( ) = 𝑇𝑊𝐶 ∙ 1 − 𝜂 ∙ 𝛽 ∙ 𝑈 ∙ 1 − 𝑛 ∙ 𝑚 ∙ 𝐿 (10) ெோ ଴ ௟௢௦௦ ଴ ௙ infinite series solution, 𝜌 represents the airfoil density, and 𝐶𝑝 ௪௔௟௟ ௪௔௟௟ represent the airfoil specific heat. The thermal conductivity heat flux into the water/ice mass for any time, 𝑡 , can be determined by solving Melt-dominated icing focuses on the fraction of ice that can melt, 𝑚 .

଴ for 𝑇 and calculating the spatial change in temperature at the ice and Equations (9) and (10) can be arranged for 𝑚 , and is shown in Eq.

଴ wall interface ( 𝑑𝑇 ⁄ 𝑑𝑛 ො at 𝑛 ො = 0). Since the value of 𝑑𝑇 ⁄ 𝑑𝑛 ො changes (11).

with time, the thermal conductivity heat flux term is transient for the " " " " situation described in Eq. (6). The average conductive heat flux, ௤ ା௤ ା௤ ି ௤ ೎೚೙ೡ ೐ೡೌ೛ ೖ೔೙೐೟೔೎ ೎೚೙೏ 𝑚 = (11) " ଴ " 𝑚 ̇ ∙ ௅ 𝑞 ത , can be taken to approximate the conductive heat transferred 𝑖𝑚𝑝 , 𝑖𝑐𝑒 ೑ ∙ ( భష 𝑛 ) ௖௢௡ௗ 𝑙𝑜𝑠𝑠 between the airfoil and icing surface over a period of time It should be noted that ice mass must initially be present for The latent heat of fusion surface energy for freezing, is shown in Eq.

melt-dominated icing to occur. Similar to the freeze mass fraction " (7) in terms of freeze fraction, 𝑛 . In the main expression, 𝑚 ̇ is ଴ ௜௠௣ , ௟௜௤ expression in Eq. (8), the kinetic energy transfer flux is a function of the mass flux of liquid water, 𝐿 is the latent heat of fusion, and 𝑛 ௙ ௟௢௦௦ 𝑛 , and therefore Eq. (11) contains two unknown parameters.

௟௢௦௦ is the fraction of ice mass lost due to bounce and erosion. The value of 𝑚 ̇ can be expressed in terms of total water content ( 𝑇𝑊𝐶 ) , melt ௜௠௣ , ௟௜௤ Mass Balance Equations Using Stagnation Icing ratio, collection efficiency at the stagnation line ( 𝛽 ), and particle mass ଴ Growth Rates velocity.

" " All parameters in the expressions for the thermodynamic model can be 𝑞 = 𝑚 ̇ ∙ ( 1 − 𝑛 ) ∙ 𝑛 ∙ 𝐿 ௟௢௦௦ ଴ ௙ ௙௥௘௘௭௘ ௜௠௣ , ௟௜௤ measured experimentally except for  , n , and n in Eq. (8) and m 0 loss 0 0 = ( 𝑇𝑊𝐶 ∙ 𝜂 ) ∙ 𝛽 ∙ 𝑈 ∙ ( 1 − 𝑛 ) ∙ 𝑛 ∙ 𝐿 ெோ ଴ ௟௢௦௦ ଴ ௙ in Eq. (11). Values of  are dependent on particle size, velocity, and (7) airfoil geometry and can be approximated using simulation [18]. For this analysis, the collection efficiency value will be approximated to Freeze-dominated icing focuses on the fraction of liquid water that can be unity,  = 1 at the stagnation point. This leaves n as the only 0 loss freeze, 𝑛 . Equations (1) and (7) can be combined and arranged to ଴ unknown for determining either n or m . The values of n , and n or 0 0 loss 0 solve for 𝑛 , and is shown in Eq. (8).

଴ m 0 , can be determined utilizing experimentally measured ice growth rates within mass balance equations. Equation (12) expresses the mass Page 4 of 17 balance for freeze-dominated icing, utilizing the mass flux that ൫ ଵିி ൯ ̇ accreted (i.e. stick). The ice growth rate, 𝑡 , is explicitly written ೗೚ೞೞ , ೟೚೟ೌ೗ ௙௥௘௘௭௘ 𝑛 = 1 − (17) ௟௢௦௦ ఎ ∙ ௡ ା ( ଵିఎ ) ಾೃ బ ಾೃ to distinguish this as the freeze-dominated icing mass balance. Both ” ” 𝑚 ̇ and 𝑚 ̇ are broken down into more measurable ௜௠௣ , ௟௜௤ ௜௠௣ , ௜௖௘ In a similar fashion, combining Eqs. (13), (14) and (15), then re- components. Equation (13) expresses the mass balance for arranging for 𝑛 provides the correlation for melt-dominated icing, ௟௢௦௦ ̇ melt-dominated icing, where the ice growth rate, 𝑡 , is explicitly ௠௘௟௧ and is shown in Eq. (18).

written to distinguish this as the melt-dominated icing mass balance.

The units of ice growth rates are in terms of distance over time.

” ௠ ̇ ∙ ൫ ଵିி ൯ ೔೘೛ , ೟೚೟ೌ೗ ೗೚ೞೞ , ೟೚೟ೌ೗ 𝑛 = 1 − (18) ௟௢௦௦ ” ௠ ̇ ∙ ( ଵି௠ ) బ ೔೘೛ , ೔೎೐ ” ” ” ̇ 𝑚 ̇ = 𝑡 ∙ 𝜌 = ൫ 𝑚 ̇ ∙ 𝑛 + 𝑚 ̇ ൯ ∙ ( 1 − 𝑛 ) ௦௧௜௖௞ ௙௥௘௘௭௘ ௜௖௘ ௜௠௣ , ௟௜௤ ଴ ௜௠௣ , ௜௖௘ ௟௢௦௦ Again, the numerator of Eq. (18) is the impinging mass flux that sticks, = 𝑇𝑊𝐶 ∙ 𝛽 ∙ 𝑈 ∙ ( 𝜂 ∙ 𝑛 + ( 1 − 𝜂 ) ) ∙ ( 1 − 𝑛 ) (12) ଴ ெோ ଴ ெோ ௟௢௦௦ while the denominator is the unmelted ice mass flux rate. Yet again, subtracting the quotient from “1” gives the fraction of ice mass flux that did not stick, or the fraction of ice mass that is lost due to bounce combined with any erosion. Eq. (18) can be also be formulated in terms ” ” ̇ ( ) ( ) 𝑚 ̇ = 𝑡 ∙ 𝜌 = ቀ 𝑚 ̇ ∙ 1 − 𝑚 ቁ ∙ 1 − 𝑛 ௦௧௜௖௞ ௠௘௟௧ ௜௖௘ ௜௠௣ , ௜௖௘ ଴ ௟௢௦௦ of 𝜂 , and is shown in Eq. (19).

ெோ = 𝑇𝑊𝐶 ∙ 𝛽 ∙ 𝑈 ∙ ൫ ( 1 − 𝜂 ) ( 1 − 𝑚 ) ൯ ∙ ( 1 − 𝑛 ) (13) ଴ ெோ ଴ ௟௢௦௦ ൫ ଵିி ൯ ೗೚ೞೞ , ೟೚೟ೌ೗ 𝑛 = 1 − (19) ௟௢௦௦ ( ଵିఎ ) ∙ ( ଵି௠ ) In Eq. (12) and (13), 𝜌 is the density of accreted ice (approximated ಾೃ బ ௜௖௘ to be 916 kg/m for this work). In Eq. (12), all of the ice mass along Solving for the unknown values of n , and n or m , provides with water mass that freezes, after bounce and erosion losses, loss 0 0 information on the amount of mass flux that is lost due to splash and contributes to the ice growth. In Eq. (13), ice mass that has not melted, runback. For freeze-dominated icing, the fraction of impinging mass after losses, contributes to ice growth. For freeze-dominated icing, Eq.

flux that does not stick due to splash and runback, 𝐹 , can be (8) and Eq. (12) are solved together, providing values for 𝑛 and 𝑛 .

௟௢௦௦ , ௦௣௟ / ௥௕ ௟௢௦௦ ଴ For melt-dominated icing, Eq. (11) and Eq. (13) are solved together, determined from the amount of impinging liquid mass flux that does not freeze, and is shown in Eq. (20).

providing values for 𝑛 and 𝑚 . The unknown values cannot be ௟௢௦௦ ଴ solved for directly and must be determined by using an iterative ” ̇ solving method. ௠ ∙ ( ଵି௡ ) ೔೘೛ , ೗೔೜ బ 𝐹 = 1 − (20) ௟௢௦௦ , ௦௣௟ / ௥௕ ” ௠ ̇ ೔೘೛ , ೟೚೟ೌ೗ Currie [8, 9] proposed a straightforward mass balance that defined a sticking efficiency, 𝑛 , which was simply the fraction of the For melt-dominated icing, the fraction of impinging mass flux that ௦௧௜௖௞ impinging mixed-phase mass flux that managed to stick to the surface. does not stick due to splash and runback, is determined using the It is defined in Eq. (14) amount of impinging liquid mass flux along with the ice mass flux that melted, and is shown in Eq. (21).

̇ ௧ ∙ ఘ ೔೎೐ 𝜂 = (14) ௦௧௜௖௞ ” ” ” ௠ ̇ ೔೘೛ , ೟೚೟ೌ೗ ௠ ̇ ା ௠ ̇ ∙ ௠ బ ೔೘೛ , ೗೔೜ ೔೘೛ , ೔೎೐ 𝐹 = 1 − (21) ௟௢௦௦ , ௦௣௟ / ௥௕ ” ௠ ̇ ೔೘೛ , ೟೚೟ೌ೗ ” Again, in Eq. (14), 𝑚 ̇ is the total impinging mass flux, which ௜௠௣ , ௧௢௧௔௟ ” ” is also the sum of 𝑚 ̇ and 𝑚 ̇ . Currie did not differentiate ௜௠௣ , ௟௜௤ ௜௠௣ , ௜௖௘ Icing Surface Equilibrium Temperature Derivation ̇ between freeze-dominated or melt-dominated icing, and therefore 𝑡 is simply the general ice growth rate. Currie defined the complementary For Eqs. (1) through Eq. (13), the thermodynamic model centered on total mass loss fraction as well, 𝐹 , which is simply the fraction ௟௢௦௦ , ௧௢௧௔௟ the energy available for freezing or melting (latent energy), where the of impinging mass flux that did not stick, and is defined in Eq. (15).

icing surface temperature was 𝑇 = 0 °C. Bartkus et al. [17] stated ௦௨௥௙ that for conditions that were sufficiently cold and freeze fraction values 𝐹 = 1 − 𝜂 (15) ௟௢௦௦ , ௧௢௧௔௟ ௦௧௜௖௞ were calculated to be n 0 = 1, then the icing surface temperature at the stagnation point could continue to decrease below 0 °C until it reached It should be noted that 𝑛 and 𝐹 represent different losses.

௟௢௦௦ ௟௢௦௦ , ௧௢௧௔௟ a thermodynamic equilibrium temperature. The derivation of that While 𝐹 is the fraction of impinging total mass that is lost, ௟௢௦௦ , ௧௢௧௔௟ equilibrium temperature is provided in Eq. (22).

𝑛 is the fraction of ice mass that is lost due to bounce and erosion.

௟௢௦௦ These two loss terms are related. Combining Eqs. (12), (14) and (15), " " " " " " 𝑞 + 𝑞 = 𝑞 − 𝑞 − 𝑞 − 𝑞 (22) then re-arranging for 𝑛 provides the correlation for ௦௘௡௦ , ௜௖௘ ௙௥௘௘௭௘ ௘௩௔௣ ௖௢௡௩ ௞௜௡௘௧௜௖ ௖௢௡ௗ ௟௢௦௦ freeze-dominated icing, and is shown in Eq. (16).

Equation (22) expands the expression in Eq. (1) by including the ” ௠ ̇ ∙ ൫ ଵିி ൯ ೗೚ೞೞ , ೟೚೟ೌ೗ ೔೘೛ , ೟೚೟ೌ೗ " 𝑛 = 1 − (16) ௟௢௦௦ ” ” sensible heat surface energy for frozen ice, 𝑞 , and allowing ௦௘௡௦ , ௜௖௘ ௠ ̇ ∙ ௡ ା௠ ̇ బ ೔೘೛ , ೗೔೜ ೔೘೛ , ೔೎೐ cooling of the iced surface below 0 °C. The other surface energy terms in Eq. (22) are the same as previously defined in Eqs. (1) through The numerator of Eq. (16) is simply the impinging mass flux that " Eq. (7). The 𝑞 term can be expanded further and is shown in sticks, while the denominator is the ice mass flux rate, which is the ௦௘௡௦ , ௜௖௘ Eq. (23).

sum of the impinging liquid water mass flux that froze and impinging ice mass flux rate. Subtracting the quotient from “1” gives the fraction of ice mass flux that did not stick. Or put another way, it is the fraction " ” 𝑞 = 𝑚 ̇ ∙ 𝐶𝑝 ∙ ൫ 𝑇 − 𝑇 ൯ (23) ௦௘௡௦ , ௜௖௘ ௦௧௜௖௞ ௜௖௘ ௦௨௥௙ ௙௥௘௘௭௘ of ice mass that is lost due to bounce and erosion. Eq. (16) can be also be formulated in terms of 𝜂 , and is shown in Eq. (17).

ெோ Page 5 of 17 In Eq. (23), 𝐶𝑝 is the specific heat of the ice which is approximated particular case of a series. A total of 16 icing tests were conducted for ௜௖௘ these four series. The tunnel exit target velocity, U , plenum to be constant with temperature, and 𝑇 is the temperature at targ ௙௥௘௘௭௘ ” pressure, p PL,targ , and plenum temperature, T PL,targ , represent the target which water freezes (i.e. 0 °C). In Eq. (24), 𝑚 ̇ is expanded further ௦௧௜௖௞ values for each Test Condition Series. The plenum relative humidity as defined by freeze-dominated icing mass balance equation. As a values, RH PL , are as measured in the tunnel plenum. As can be seen reminder, sub-freezing surface temperatures can only occur in from Table 1, tests were conducted at 3 different target velocities freeze-dominated icing. In addition, 𝑛 = 1 by definition for this ଴ (85 m/s, 135 m/s, and 185 m/s), and at two target total pressures sensible heat cooling, and therefore portions of the expression can be (44.8 kPa and 87.5 kPa). The target air flow speeds equate to Mach ” reduced. In the final formulation in Eq. (24), 𝑚 ̇ is expressed 𝑖𝑚𝑝 , 𝑡𝑜𝑡𝑎𝑙 values of 0.25, 0.40, and 0.56 respectively. Not listed in the table, but in terms of total water content, collection efficiency at the stagnation important to mention, the injected particle size distribution was held line, particle mass velocity, and 𝑛 .

௟௢௦௦ constant for each test with an approximate initial median volumetric diameter, MVD , of 20 μm. In addition, the injected TWC was held constant for each test, with a target value of 2.0 g/m . The column " 𝑞 = ௦௘௡௦ , ௜௖௘ labeled Escort File # (or Esc# for short) in Table 1 is the test run ” ” number system used during testing. It should be noted that since flow ( ) ൫ 𝑚 ̇ ∙ 𝑛 + 𝑚 ̇ ൯ ∙ 1 − 𝑛 ∙ 𝐶𝑝 ∙ ൫ 𝑇 − 𝑇 ൯ ௜௠௣ , ௟௜௤ ଴ ௜௠௣ , ௜௖௘ ௟௢௦௦ ௜௖௘ ௦௨௥௙ ௙௥௘௘௭௘ velocity in the plenum is slow, plenum values can be approximated to ” ( ) = ൫ 𝑚 ̇ ൯ ∙ 1 − 𝑛 ∙ 𝐶𝑝 ∙ ൫ 𝑇 − 𝑇 ൯ be total conditions.

௟௢௦௦ ௜௖௘ ௦௨௥௙ ௙௥௘௘௭௘ ௜௠௣ , ௧௢௧௔௟ Table 1. Test conditions for the four RH sweeps (Test Condition Series).

( ) = 𝑇𝑊𝐶 ∙ 𝛽 ∙ 𝑈 ∙ 1 − 𝑛 ∙ 𝐶𝑝 ∙ ( 𝑇 − 𝑇 ) (24) ଴ ௟௢௦௦ ௜௖௘ ௦௨௥௙ ௙௥௘௘௭௘ TCS U p T Case RH Escort targ PL,targ PL,targ PL kPa °C File With the solved freeze fraction value of n = 1 for these sufficiently # m/s (psia) (°F) % # cold cases, the unknown values become n and 𝑇 . The surface A 0.4 60 loss ௦௨௥௙ " " " " 44.8 7.2 B 34.0 228 energy terms of 𝑞 , 𝑞 , 𝑞 , and 𝑞 in Eq. (22) are all ௖௢௡௩ ௘௩௔௣ ௦௘௡௦ , ௜௖௘ ௖௢௡ௗ 1 85 (6.5) (45) C 36.9 72 functions of 𝑇 . Therefore, the unknown values of n and 𝑇 loss ௦௨௥௙ ௦௨௥௙ D 40.8 235 must be determined by iteratively solving Eq. (12) and Eq. (24), where A 20.6 170 experimental growth rate data is used with respect to Eq. (12). The 44.8 7.2 B 25.2 593 2 135 resulting solved value of 𝑇 is the thermodynamic equilibrium icing ௦௨௥௙ (6.5) (45) C 30.2 171 surface temperature. D 34.5 169 A 28.1 187 44.8 7.2 B 29.5 611

Model Assessment with Experimental Data

3 185 (6.5) (45) C 31.4 214 D 33.2 205 The following sections present an assessment of the thermodynamic A 5.2 618 model utilizing experimental data from fundamental icing physics B 10.7 249 87.5 7.2 4 135 testing conducted at the NASA PSL in 2018. The thermodynamic (12.7) (45) C 15.0 243 model is evaluated using data measured at various velocities, D 19.5 253 pressures, humidities, temperatures, and melt ratios. Struk et al. [3] provides a description of the test facilities and experiments in an All plenum values in Table 1 refer to conditions upstream of the tunnel accompanying paper. The experiments generated a set of ice shapes spray bars. It should be noted that while the conditions upstream of the spray bar in the tunnel plenum were maintained constant, as it has on a NACA 0012 airfoil under well-characterized conditions. At the center of the icing experiments were four sets of relative humidity ( RH ) been previously reported [2-5], conditions downstream at the test section change with the activation of the spray cloud. Modeling efforts sweeps. The model will be evaluated utilizing accretion data on the [19-23] show that a thermodynamic interaction between the icing 0.267 m (10.5 in) chord airfoil, where tests were conducted with a 0º angle of attack (AOA). The assessment will focus on times after which cloud and flowing air results in changes in total air temperature, 𝑇 , ଴ and humidity at the test section, compared with pre-spray conditions.

the initial temperature difference between the airfoil leading edge temperature and ice temperature became negligible. According to In addition, changes occur in total water content, melt ratio and cloud Bartkus [17] this initial conductive heat transfer between the airfoil and particle size at the test section, compared with the initial spray conditions in the plenum. For the NASA 2018 experiments, whereas accreting ice occurs primarily during the first 20 s of airfoil exposure to the icing cloud. This assessment will look at accretions beyond the only humidity content was varied at the tunnel inlet (plenum), other initial transient time period. The following will present sections on the parameters changed, most notably 𝑇 , 𝑇𝑤𝑏 , 𝜂 , 𝑇𝑊𝐶 , and the cloud ଴ ଴ ெோ experimental conditions, experimental analysis, and finally the particle size distribution (and MVD ). Due to the highly-coupled thermodynamic model evaluation. parameter set, the ideal experimental situation where one variable is isolated and varied independently was not possible.

Experimental Conditions Table 2 shows the measured and calculated conditions at the test section (i.e. just forward of the airfoil leading edge) after the icing Ice shapes were generated across a series of four different flow cloud was activated, for all 16 icing tests. See Struk et al. [3] for a conditions where the relative humidity in the plenum of PSL was the description of how these conditions were measured and calculated. In primary parameter varied. Table 1 provides the initial test conditions Table 2, the humidity parameter, MMR , is the mass mixing ratio, which for the four RH test sweeps. Four icing tests (labeled as Case A through is the ratio of vapor mass to dry air mass. Not listed in Table 2 is that Case D) were conducted for each RH test sweep. For convenience, the cloud MVD just forward of the airfoil leading edge was measured each RH test sweep will be referred to as Test Condition Series (TCS-1 to be approximately 35 μm for every test [3, 24]. The value varied by through TCS-4) for the remainder of the paper, and may be coupled just a few microns from case to case. It can be seen from Table 2 that with an individual case letter (i.e. TCS-1.A) for quick reference to a Page 6 of 17 within each TCS, the values of 𝑇 , 𝑇𝑤𝑏 , 𝜂 , MMR, and 𝑇𝑊𝐶 ଴ ଴ ெோ increased as RH PL was increased. It should be noted that total conditions are provided in Table 1, but several of the energy expressions provided utilize static conditions.

Table 2. Conditions as measured and calculated just forward of the airfoil leading edge as the icing cloud was activated.

Test U T p MMR Twb TWC 𝜂 0 0 0 ெோ Case (m/s) (°C) (kPa) (g/kg) (°C) (g/m³) (-) Test Condition Series 1 A 84 -1.0 44.7 3.9 -5.5 1.2 0.03 B 83 1.8 44.7 7.2 -0.8 2.1 0.37

a

C 83 2.1 44.8 7.5 -0.4 2.1 0.57 D 84 2.2 44.7 7.9 0.0 2.2 0.93 Test Condition Series 2 A 133 0.5 44.8 5.9 -2.6 3.2 0.14 B 133 1.1 44.7 6.3 -1.9 3.4 0.21 C 133 2.0 44.8 6.8 -1.1 3.6 0.28 D 133 2.1 44.8 7.1 -0.8 3.8 0.51 Test Condition Series 3 A 182 3.0 44.8 6.2 -1.4 4.1 0.29 B 182 3.6 44.7 6.5 -1.0 4.2 0.31 C 182 3.6 44.8 6.6 -0.8 4.2 0.41 D 182 3.9 44.8 6.7 -0.6 4.3 0.53 Test Condition Series 4

b

A 133 2.4 87.5 2.6 -1.4 2.9 0.32 Figure 1. Image of ice accretion from the a) span-view and b) profile-view.

B 133 2.5 87.5 2.9 -0.9 3.2 0.40 C 133 2.5 87.5 3.1 -0.7 3.3 0.65 Figure 2 through Fig. 5 show profile-view ice accretions of the four D 133 3.0 87.5 3.3 -0.1 3.5 0.91 Test Condition Series. The ice shapes at 240 s, and at the end of test (generally near 600 s) are shown for each case. Two tests were Experimental Analysis terminated short of 600 s, as no accretion was observed on the mid- section of the airfoil. These two tests are noted in the captions. Three Analysis of the accreted ice was performed for each icing test and is of the images show ice accretion that occurred on the airfoil extension utilized for this paper to better understand the icing process. Videos of (denoted by red ovals). These ice accretions on the extension were the ice accretion were recorded at two perpendicular angles. Figure 1 dismissed from ice accretion analysis since they did not occur on the shows images of the ice accretion from both camera angles. Figure 1a mid-section of the airfoil. Cases A through Case D for each Test shows the accretion along the span (span-view), Fig. 1b shows the Condition Series (TCS-1 through TCS-4) along with the time in two-dimensional accretion from the side of the airfoil (profile-view). seconds from the start of spray activation are shown in the bottom left In Fig. 1a, the light gray mid-section region (~10.5 cm in width) is the corner of each image. The earlier image (left) also lists the melt ratio main icing target area (a hollow titanium alloy shell), while the dark in the bottom corner as well for reader convenience. These images are regions on either side of the mid-section are the airfoil extensions shown to provide the reader a general idea of ice shape and growth. In (solid aluminum). Leading edge ice growth rates were measured from general, the figures show that as melt ratio increased for each Test both perspectives, however, only the span-view ice growth rate at the Condition Series, the ice shape grew in size and also changed from a midspan is utilized for this analysis. The span-view camera was able pointy wedge-like shape to a blunt shape. Precursors of double-horn to capture the entire ice growth, including the initial growth, while the shapes are noticeable and more prevalent in the higher melt ratio profile-view was unable to capture the initial ice growth conditions, and higher velocity conditions. A full double-horn ice (approximately the initial 2.5 mm) due to blind spots that result from shape formed in Case D of TCS-4. Some of the ice shapes became optical perspective effects. Profile-view images were backlit to quite irregular beyond 240 seconds of ice cloud exposure.

provide good contrast between the iced airfoil and background for imaging analysis purposes. Ice shape analysis also included measuring the leading edge ice accretion angle. This was measured utilizing the profile-view images and generating two line segments along the leading edge. Each line segment extended about 4 mm on opposite sides of the leading edge. The two line segments produced the leading edge ice accretion angle, φ , and can be seen in Fig. 1b. Thermocouple (T/C) data that measured airfoil surface temperatures was also analyzed. Several thermocouples were affixed to the surface of the airfoil, however the leading edge thermocouple at the airfoil midspan is the primary thermocouple of interest for this analysis.

Page 7 of 17 3 cm 3 cm A – 145 s – 𝜂 = 0.03 A – 198 s A – 632 s A – 240 s – 𝜂 = 0.29 ெோ ெோ B – 240 s – 𝜂 = 0.37 B – 610 s B – 607 s B – 240 s – 𝜂 = 0.31 ெோ ெ ோ C – 240 s – 𝜂 = 0.57 C – 240 s – 𝜂 = 0.41 C – 623 s C – 709 s ெோ ெோ D – 240 s – 𝜂 = 0.93 D – 240 s – 𝜂 = 0.53 D – 624 s D – 619 s ெோ ெோ Figure 2. Profile ice shapes for TCS-1 at 240 s and at the end of test (generally Figure 4. Profile ice shapes for TCS-3 at 240 s and at the end of test (generally near 600 s). The red circles show accretion on the airfoil extension. The Case near 600 s). The Case (A through D) along with time of ice exposure is shown (A through D) along with time of ice exposure is shown in the bottom left of in the bottom left of every image.

every image.

3 cm 3 cm A – 240 s – 𝜂 = 0.32 A – 613 s A – 501 s A – 240 s – 𝜂 = 0.14 ெோ ெோ B – 615 s B – 240 s – 𝜂 = 0.21 B – 240 s – 𝜂 = 0.40 B – 608 s ெோ ெோ C – 607 s C – 608 s C – 240 s – 𝜂 = 0.28 C – 240 s – 𝜂 = 0.65 ெோ ெோ D – 636 s D – 240 s – 𝜂 = 0.51 D – 240 s – 𝜂 = 0.91 D – 608 s ெோ ெோ Figure 3. Profile ice shapes for TCS-2 at 240 s and at the end of test (generally Figure 5. Profile ice shapes for TCS-4 at 240 s and at the end of test (generally near 600 s. The red circle shows accretion on the airfoil extension. The Case (A near 600 s). The Case (A through D) along with time of ice exposure is shown through D) along with time of ice exposure is shown in the bottom left of every in the bottom left of every image.

image.

Page 8 of 17 Figure 6. Analyzed data for TCS-1 (Case A through Case D). Analysis includes Figure 7. Analyzed data for TCS-2 (Case A through Case D). Analysis includes midspan ice thickness (with linear fits), leading edge ice accretion angle, and midspan ice thickness (with linear fits), leading edge ice accretion angle, and leading edge thermocouple temperature data.

leading edge thermocouple temperature data.

Page 9 of 17 Figure 8. Analyzed data for TCS-3 (Case A through D). Analysis includes Figure 9. Analyzed data for TCS-4 (Case A through D). Analysis includes midspan ice thickness (with linear fits), leading edge ice accretion angle, and midspan ice thickness (with linear fits), leading edge ice accretion angle, and leading edge thermocouple temperature data. leading edge thermocouple temperature data.

Page 10 of 17 Figure 6 through Fig. 9 show analyzed data for Case A through D for and colder cases (lowest T and Twb cases), the T/C temperature 0 0 each Test Condition Series (TCS-1 through TCS-4). Analysis includes decreased below 0 °C. There are two likely reasons for these midspan ice thickness (from span-view profile images), leading edge sub-freezing readings. The first reason for sub-freezing thermocouple ice accretion, angle and leading edge thermocouple data with respect readings, is that for the coldest T and Twb cases, the energy surface 0 0 to ice cloud activation time. Several linear fits are fitted to the ice balance equations will show that the icing surface should decrease thickness for each case. The slope of each line represents the below 0 °C. These surface energy equilibrium temperature equations approximate ice growth rate (mm/s) for the respective time period. were previously discussed, and will be evaluated later in the paper. The second explanation is due to conduction from colder sections of the There are several observations and points to be made for Figs. 6 – 9. airfoil, downstream of the leading edge, and generally relates to higher Firstly, of all the cases shown in Figs. 6 – 9, only case TCS-1.A did velocity cases. The temperature recovery factor aft of the leading edge not accrete any ice, presumably because the melt ratio for the test was in the boundary layer is expected to be less than one. Due to lower too low to accrete ice. A second point, the leading edge ice thickness temperature recovery factors, airfoil metal temperatures will be colder for TCS-4.D (the double-horn test) in Fig. 9 was only measured for the downstream of the leading edge. Eventually conduction of heat from initial period (ice thickness = 0 mm), prior to the onset of the the leading edge towards those lower metal temperature regions drive double-horn ice growth. The view of the leading edge (inside the the leading edge T/C temperature down. This is most prevalent for the double-horn) was obscured by a cloud debris field for a majority of the highest velocity condition (TCS-3), where temperature recovery will test, and therefore was not measured (along with the leading edge be the lowest aft of the leading edge. The hypothesis of metal accretion angle). A single measurement of leading edge ice thickness conduction within the airfoil will be evaluated later in the paper as and angle was made when the cloud was turned off and a clean view well.

of the accretion was imaged at the end of TCS-4.D. Next, aside for two cases (TCS-1.B and TCS-1.C) in Fig. 6, each test experienced a Sticking Efficiency Analysis period of no growth or slow growth at or near the start of ice cloud activation. This period of leading edge no/slow growth lasted It is instructive to look at ice accretions at several times, namely at anywhere from 60 s to 240 s in some cases. It was observed that ice 60 s, 240 s, and 600 s, to see how sticking efficiency changes with accretion further downstream on the airfoil grew forward during these ̇ ̇ ̇ time. The corresponding ice growth rates – 𝑡 , 𝑡 , and 𝑡 ଺଴௦ ଶସ଴௦ ଺଴଴ leading edge no/slow growth periods. As that downstream ice reached respectively – for all 16 cases are shown in Table 3. The same growth near the leading front plane of the airfoil, the leading edge ice growth rate may apply for multiple periods for any individual case.

experienced an increase in rate. A final observation from these figures is that the ice did not reach a maximum ice size for any case where ice Table 3. Ice growth rates for all cases as measured at 60 s, 240 s, and 600 s.

accretion occurred.

̇ ̇ ̇ 𝑡 𝑡 𝑡 ଺଴ ௦ ଶସ଴ ௦ ଺଴଴ ௦ Case Leading edge ice accretion angle was measured for periods where the (mm/s) (mm/s) (mm/s) leading edge was visible for Fig. 6 through Fig. 9. Utilizing the Test Condition Series 1 profile-view, the ice accretion only became visible once it reached A 0.000 0.000 N/A approximately 2.5 mm in thickness, as optical perspective effects (as discussed earlier) created a blind region. Therefore, it was not possible B 0.033 0.033 0.033 to measure this angle for the beginning of each test, sometimes up to C 0.043 0.043 0.043 several hundred seconds for some cases. The ice shape near the leading D 0.000 0.020 0.043 edge became irregular at times, and therefore the accretion angle was Test Condition Series 2 not measured for these periods. In general, the ice shape became more irregular as 𝜂 and U increased. Precursors to double-horn shapes A 0.000 0.004 0.004 ெோ (and fully developed double-horn shapes) had φ values greater than B 0.000 0.013 0.013 180 ° . For TCS-1 and TCS-2, some trends are noticeable regarding the C 0.013 0.013 0.069 leading edge accretion angle, φ . First is that φ reaches a near D 0.015 0.015 0.081 steady-state angle towards the end of each test. Also, φ increased in Test Condition Series 3 value as 𝜂 increased for TCS-1 and TCS-2. The lower φ values are ெோ likely due to erosion effects at lower 𝜂 values. Comparing φ with ெோ A 0.000 0.035 0.013 leading edge growth rate, one can see that larger leading edge ice angle B 0.000 0.056 0.056 values are accompanied by higher ice growth rates. TCS-2.C illustrates C 0.002 0.056 0.056 this coupled trend well. During the slow growth period, φ reaches its lowest value (perhaps the decreasing φ value is indicative of some D 0.028 0.069 0.069 erosion effects during this slow growth period), then as φ increased, Test Condition Series 4 ice growth rate increases. Both the angle and growth rate reach a steady A 0.005 0.041 0.041 value beyond the 360 s mark in the test. Another trend noticeable is B 0.010 0.064 0.064 that larger φ values occurred at higher velocities. Finally, larger φ values occurred for the higher pressure cases (TCS-4) compared to the C 0.000 0.063 0.063 lower pressure cases (TCS-2).

D 0.000 N/A N/A Thermocouple temperature data at the airfoil midspan leading edge was plotted for each case in Fig. 6 through Fig 9. In general, the pre-spray T/C temperature for all 16 cases measured about 6 °C above freezing. This is the approximately the pre-spray air total temperature.

The T/C responded quickly and reached 0 °C within seconds of cloud activation. This reading of 0 °C is generally expected for mixed-phase icing tests. For some tests, generally the highest velocity cases (TCS-3) Page 11 of 17 Sticking efficiency is plotted against melt ratio for the four test cases 0.4 at 60 s, 240 s and 600 s in Fig. 10.There are several key points to take TCS-1 away from Fig. 10. First is that the sticking efficiency was the lowest a - 60 s TCS-2 for the earliest time period (Fig. 10a - 60 s), while the highest 𝜂 ௦௧௜௖௞ TCS-3 values occurred towards the end of the test (Fig. 10c - 600 s). This low TCS-4 sticking efficiency in Fig. 10a corresponds with the no/slow growth 0.2 that occurred during the first few minutes for most tests, which was highlighted in the previous experimental analysis section. The two points of TCS-1 (Case B and Case C) with higher 𝜂 values during ௦௧௜௖௞ the early time period correspond to points that did not experience a no/slow growth. As noted earlier, when ice downstream of the leading Sticking Efficiency (non-dim) edge grew forward reaching the leading edge plane, the leading edge 0.0 ice growth rates increased. Fig. 10b and 10c illustrate this higher 0.0 0.2 0.4 0.6 0.8 1.0 Melt Ratio (non-dim) growth rates with higher sticking efficiency values. A plateau of icing severity (where severity relates to higher 𝜂 values) begins to ௦௧௜௖௞ 0.4 emerge at later times for each Test Condition Series. According to TCS-1 b - 240 s Fig. 10c, the 𝜂 plateau value depends on the TCS. The plateau ௦௧௜௖௞ TCS-2 decreases with increasing velocity where 𝜂 plateau values are ௦௧௜௖௞ TCS-3 approximately 0.21, 0.15, and 0.07 for TCS-1, TCS-2, and TCS-3 TCS-4 respectively. Again, the measured velocities for these three Test 0.2 Condition Series are U = 84 m/s, 133 m/s, and 182 m/s respectively.

Fig. 10c shows that approximately the same plateau 𝜂 values were ௦௧௜௖௞ measured for the two different total pressure cases. Again, the total pressures for TCS-2 and TCS-4 are p = 44.8 kPa and 87.5 kPa, respectively. It is not clear if an upper 𝜂 limit exists with ெோ Sticking Efficiency (non-dim) 0.0 mixed-phase icing for the conditions run in these experiments. The 0.0 0.2 0.4 0.6 0.8 1.0 point plotted with the highest melt ratio (TCS-1.D, 𝜂 = 0.92) may ெோ Melt Ratio (non-dim) have been experiencing supercooled liquid icing. A sufficient number of tests were not conducted to identify if a maximum 𝜂 mixed-phase 0.4 ெோ icing limit exists. Despite that lone point, the plateau, in general ranged TCS-1 c - 600 s from approximately 𝜂 = 0.30 to at least 𝜂 = 0.65. Fig. 10c shows ெோ ெோ TCS-2 that the onset of icing is pushed to the right (towards higher melt ratio TCS-3 values) with increasing velocity, although the uncertainty in melt ratio TCS-4 measurements precludes this from being definitive. Should this trend 0.2 be correct, the shift to more severe icing with respect to melt ratio may be explained by greater amounts of erosion occurring at higher velocities at lower 𝜂 values.

ெோ Sticking Efficiency (non-dim) 0.0 0.0 0.2 0.4 0.6 0.8 1.0 Melt Ratio (non-dim) Figure 10. Sticking efficiency vs melt ratio at for all 16 cases at a) 60 s, b) 240 s, and c) 600 s.

It should be noted that the case where no ice accretion was observed and terminated after 198 s of cloud exposure (TCS-1.A) is plotted in all three plots. It is momentarily assumed that no ice would grow at any time under those conditions, and is plotted in the graphs to illustrate that there exists a lower melt ratio boundary where icing will not occur. The fully developed double-horn test (TCS-4.D) is not plotted in Fig 10b and 10c as ice growth rates at the leading edge of the airfoil could not be measured.

Thermodynamic Model Assessment As mentioned earlier, this assessment focuses on time periods after the initial temperature transient period. The thermodynamic model will be evaluated at 240 s as this was generally a good compromise with ice accretion size. Zero or small ice accretion rates occurred early in most tests, while large irregular shapes that did not conform to the dry airfoil shape often occurred at later times in a test. Not many interesting results can be extracted from zero or very slow growth rates. Also, when ice shapes become excessively large and irregular, the Page 12 of 17 coefficients of heat and mass transfer used in the model (which Mass Loss Fraction Analysis assumes an approximate dry airfoil geometry) become less applicable.

For these reasons, the model will be evaluated at 240 s.

Running the thermodynamic model to determine the values of n loss , and n 0 or m 0 , also allowed for the determination of what fraction of the Freeze Fraction and Melt Fraction Analysis impinging ice is lost due to splash and runback and what fraction is lost by bounce and erosion. Figure 12 shows loss fractions plotted The thermodynamic ice crystal icing model was run for 15 of the 16 against melt ratio. The large, solid symbols represent the total mass cases (the fully developed double-horn case TCS-4.D was not fraction lost as defined by Eq. (15). This again is simply the fraction evaluated) using the test conditions provided in Table 2, along with the of the impinging mass flux that did not stick. The smaller empty ̇ ice growth rate at 240 s, 𝑡 , presented in Table 3. Figure 11 shows symbols represent the fraction of the impinging mass flux that did not ଶସ଴௦ stick due to splash and runback. These values are determined by the calculated freeze fraction or melt fraction as determined from the model, plotted against the measured total wet-bulb temperature just Eq. (20) if the icing event was freeze-dominated, and by Eq. (21) if the upstream of the airfoil. The large, solid symbols in Fig. 11 represent icing event was melt-dominated. The difference between 𝐹 ௟௢௦௦ , ௧௢௧௔௟ cases that experienced freeze-dominated icing, while the smaller and 𝐹 represents the fraction that was lost due to bounce and ௟௢௦௦ , ௦௣௟ / ௥௕ empty symbols represent melt-dominated cases. The color of the erosion. It can be seen in Fig. 12 that 𝐹 increases with ௟௢௦௦ , ௦௣௟ / ௥௕ symbol is grouped with the Test Condition Series. It can be seen that increasing 𝜂 . This is generally intuitive as there is more liquid to be ெோ the there is a transition between melt-dominated and freeze-dominated lost at higher 𝜂 values. With only small variations in the trend, ெோ icing around -0.5 to -1.0 °C (depending on condition). In general, Twb Fig. 12 suggests that 𝜂 is the most domininat factor is determining ெோ is a good first-order indicator for determining the type of icing that will 𝐹 as compared to any effects that U and p 0 contribute. The ௟௢௦௦ , ௦௣௟ / ௥௕ occur. Sub-freezing Twb 0 values generally indicate freeze-dominated 𝐹 values for the highest velocity case (green empty diamond) ௟௢௦௦ , ௦௣௟ / ௥௕ icing, while above freezing values indicate melt-dominated icing.

is consistently slightly greater than the other cases and may play as a Melt-dominated icing shifts slightly into the sub-freezing Twb domain second-order factor in determining the fraction lost due to splash and due to the transfer of kinetic energy from a fast moving particle to a runback. This is intuitive as higher velocities contribute more to the stationary ice accretion surface. It can be seen that higher velocity kinetic energy flux, which promotes more melt.

cases pushed further into the sub-freezing Twb 0 domain, up to -0.8 °C for one of the highest velocity cases (TCS-3.B). A slight trend is noticeable with the freeze-dominated (large, solid) points. As Twb 1.0 decreases, the freeze fraction increases. This is expected for two reasons. The first is that the colder the condition, the more readily the 0.8 liquid water will freeze into ice. Secondly, since it was not possible to isolate conditions, as humidity in the tunnel plenum was decreased, not 0.6 only did Twb decrease, but so did 𝜂 . With small amounts of liquid ெோ TCS-1: total TCS-2: total water available to freeze at lower 𝜂 values, the greater the freeze ெோ TCS-3: total fraction will be as well. It is noted that while there are not a large 0.4 TCS-4: total number of melt-dominated icing test points to show on Fig. 11, an TCS-1: spl/rb inflection of the freeze-dominated points around -0.5 °C is expected.

0.2 TCS-2: spl/rb The same reasoning mentioned for the freeze-dominated trend applies Loss Fractions (non-dim) TCS-3: spl/rb to the melt-dominated regime. As Twb increases, the melt fraction will 0 TCS-4: spl/rb 0.0 increase as ice will more readily melt at warmer conditions. Similarly, 0.0 0.2 0.4 0.6 0.8 1.0 due to the nature of many parameters being coupled in this icing tunnel, Melt Ratio (non-dim) namely 𝜂 and Twb 0 , a warmer environment will be coupled with a ெோ Figure 12. Total mass loss fraction and mass loss fraction due to splash and higher 𝜂 value, and with less ice to melt, the greater the melt fraction ெோ runback, plotted against melt ratio.

value will be.

It should be noted that the overlap between the two loss fractions for the highest 𝜂 case (TCS-1.D) suggests that there is not enough ice ெோ 1.0 in the impinging ice to create the amount of ice growth measured. This suggests either some supercooling liquid water accretion occurred, or 0.8 that the measurement and subsequent calculation of 𝜂 was too high.

ெோ For example the points would be equal in loss fraction if 𝜂 was ெோ calculated to be 0.90 instead of 𝜂 = 0.93.

ெோ 0.6 TCS-1: Freeze Frac TCS-2: Freeze Frac Icing Surface Equilibrium Temperature Analysis TCS-3: Freeze Frac 0.4 TCS-4: Freeze Frac TCS-1: Melt Frac A subroutine was added to the thermodynamic ice crystal icing model TCS-2: Melt Frac to determine the icing surface equilibrium temperature, T . The icing surf 0.2 TCS-3: Melt Frac surface temperature is calculated to be 0 °C when the freeze fraction TCS-4: Melt Frac or melt fraction equals a value less than unity (if m = 1, there is no 0 Freeze or Melt Fraction (non-dim) 0.0 icing as all ice has melted). If n = 1, conditions are sufficiently cold -6 -4 -2 0 2 to affect the sensible heat at the accretion surface. The expressions O Total Wet-Bulb Temperature ( C) derived in the Icing Surface Equilibrium Temperature Derivation Figure 11. Freeze fraction or melt fraction as determined from the model, section are used to find T surf . For this exercise, T surf was calculated at plotted against total wet-bulb temperature.

60 s, 240 s, and 600 s for all possible cases. These calculated values are compared to corresponding airfoil leading edge thermocouple Page 13 of 17 decrease in temperature recovery in the boundary layer aft of the temperature data, T and are shown in Table 4. An important T/C leading edge is even more significant. It is this reasoning that explains distinction to note is that T surf is a calculated value at the leading edge the continuously decreasing leading edge T/C temperature trend.

of the icing surface, while T is a measured value located at the T/C Finally, moving to the TCS-4, sub-freezing temperatures are measured leading edge of the airfoil. These two temperatures are more closely and calculated for the two coldest cases. The temperature trends related when the ice thickness at the leading edge is small. Therefore, between prediction and experiment as seen in TCS-4.A and TCS-4.B the comparisons are less likely to agree later in an icing test when the can be explained similarly as described for the TCS-3 series.

ice accretion is large. The comparisons, however, are still made at these later times to illustrate some points.

Evidence is provided for the hypothesis that heat conducting from the leading edge to colder regions on the metal airfoil contributed to Table 4. Comparison of calculated icing surface temperature (model) with experimentally measured thermocouple temperature data (exp.) at 60 s, 240 s, reduced leading edge thermocouple readings later into the test.

and 600 s, for all applicable test cases.

Thermocouples located downstream of the leading edge at the airfoil midspan corroborate this hypothesis. Figure 13 shows the location of T , T , T , T , T , T , T/C surf T/C surf T/C surf the leading edge thermocouple (T/C_03) and the location of the next (exp.) (model) (exp.) (model) (exp.) (model) closest thermocouple aft of the leading edge (T/C_06). The location of Case 60 s 60 s 240 s 240 s 600 s 600 s T/C_06 is in a region where the temperature recovery factor in the (°C) (°C) (°C) (°C) (°C) (°C) boundary layer is expected to be less than unity. In general, the Test Condition Series 1 1/2 recovery factor is expected to be approximately Pr for laminar flow A -4.0 -5.1 -4.8 -5.1 N/A N/A 1/3 and Pr for turbulent flow [25], where Pr is the Prandtl number ( Pr = B 0.0 0.0 0.0 0.0 0.0 0.0 0.71 for air). Figure 14 illustrates this reduced temperature C 0.0 0.0 0.0 0.0 0.0 0.0 downstream of the airfoil leading edge. Figure 14 recreates the ice thickness, leading edge angle, and thermocouple data of TCS-4.A, but D 0.0 0.0 0.0 0.0 0.0 0.0 also includes the downstream thermocouple reading (T/C_06) along Test Condition Series 2 with the leading edge thermocouple reading (T/C_03). Fig. 14 shows A -0.5 -2.5 -0.6 -2.5 -0.6 -2.5 that T/C_06 reads colder than T/C_03 shortly after the activation of the B 0.0 -1.8 -0.2 -1.7 -1.1 -1.7 cloud and. This lower temperature will tend to eventually conduct heat C 0.0 0.0 0.0 0.0 -0.1 0.0 away from leading edge of the airfoil, and in this example reducing the D 0.0 0.0 0.0 0.0 0.0 0.0 leading edge to sub-freezing temperatures.

Test Condition Series 3 A 0.0 -1.2 -0.1 0.0 -2.3 0.0 B 0.0 -0.8 -0.7 0.0 -2.7 0.0 C 0.0 -0.7 -0.8 0.0 -2.7 0.0 D 0.0 0.0 -0.5 0.0 -1.9 0.0 Test Condition Series 4 A 0.0 -1.8 -0.2 0.0 -1.6 0.0 B 0.0 -1.3 0.0 0.0 -0.4 0.0 C 0.0 0.0 0.0 0.0 -0.1 0.0 D 0.0 0.0 0.0 N/A 0.0 N/A The following will assess the predicted temperature with T/C data by stepping through each Test Condition Series. Beginning with TCS-1, only Case A shows calculated and measured sub-freezing temperatures. Even though no ice accreted on the surface, sub-zero values were measured. The airfoil, even if just barely wet from the mostly ice impinging cloud, acted as a wetted surface, and registered Figure 13. Locations of the leading edge thermocouple and the next closest thermocouple aft of the leading edge of the NACA 0012 airfoil.

values near the total wet-bulb temperature. Indeed, for this case T surf was calculated to be the total wet-bulb temperature. Cases B through D for TCS-1 were measured and calculated to be 0 °C at all times.

There is good agreement for the measured and calculated values for this test series. Continuing to TCS-2, the two coldest cases (TCS-2.A and TCS-2.B) show calculated and measured sub-freezing temperatures. There is fair agreement in temperature value for this test series as well. Next is the highest velocity series, TCS-3. Some discrepancies occur in this test series. With respect to calculated model values, the three coldest cases (TCS-3.A through TCS-3.C) show sub- freezing temperatures at 60 s. However, as ice growth rate increased from no/slow growth rates to higher growth rates, the addition of particle kinetic energy to the icing surface increased the surface temperature at later times (240 s and 600 s). All four T/C temperatures initially start at 0 °C, but continue to decrease in temperature as time Figure 14. Recreation of the ice thickness, leading edge angle, and progresses for each test. It is believed that heat is conducted in the thermocouple data of TCS-4.A, but also shows airfoil thermocouple readings metal from the airfoil leading edge to colder metal regions downstream downstream of the leading edge, T/C_06 along with the leading edge on the airfoil, due to reduced temperature recovery factors aft of the thermocouple reading, T/C_03.

leading edge. At higher velocities, such as is the case with TCS-3, the Page 14 of 17 that occurred with a greater fraction of impinging particles being

Discussion

glaciated at lower 𝜂 values. The accretion angle plays an important ெோ role in ice growth. As reported by Currie [9], most tests reached a There are many important points to discuss regarding the results of this steady-state ice size, except where the 𝜂 value remained non-zero ௦௧௜௖௞ paper. The first is the explanation for the time dependent ice growth for all oblique impingement angles, in which case the ice grew rates that were observed and measured for the majority of the tests.

indefinitely. All ice accretions grew without reaching a size limit Most tests underwent a no/slow ice growth period for the first several during the ~600 s test runs reported in this paper (aided by the minutes of the test. It was observed that ice accretion, downstream of secondary icing front). This indefinite ice growth is in agreement with the leading edge grew forward toward the leading edge plane. This Baumert.

second front created a larger leading edge surface, which led to a greater leading edge ice growth rate, and subsequently a higher

Summary/Conclusions

sticking efficiency. For some tests, the leading edge thermocouple, after reaching a freezing temperature, began to slowly decrease to sub-zero temperatures. It is hypothesized that the conduction of heat There were several key findings in this paper. First, the thermodynamic from the airfoil leading edge to colder regions on the metal airfoil, ice crystal icing model showed that Twb is a good predictor for which is colder due to a reduced temperature recovery factor, determining the type of icing that will accrete. To the first order, contributed to the decreasing thermocouple reading. In addition, it is freeze-dominated icing occurs for sub-freezing Twb values, while likely that this local colder region downstream of the airfoil leading melt-dominated icing occurs for Twb above 0 °C values. The kinetic edge aided in accreting ice in this region, which allowed for the energy surface flux provides additional energy to the icing surface, downstream ice to grow forward towards the leading edge plane. These which slightly shifts the melt-dominated regime to small sub-freezing are phenomena that the thermodynamic model does not capture since Twb values. The greater the kinetic energy flux (i.e. higher velocity), it focuses on the leading edge. These comprehensive explanations, the greater the shift into the sub-freezing Twb 0 domain. Also, the lower however, are necessary to understand the changing ice growth rates the Twb , the greater the fraction of liquid water that will freeze. If that were observed that are ultimately used in the model.

sufficiently cold, the freeze fraction reaches n 0 = 1 and the icing surface temperature can shift from mixed phase temperatures of 0 °C In ice crystal icing, there is a general acceptance that a minimum and to fully glaciated sub-freezing temperatures. Additions to the model maximum melt ratio exists for accretion to occur. The minimum limit predicted this sensible energy change. The predicted icing surface is described as the limit below which too little melt occurs preventing temperatures were in good agreement with experimental temperature the ice from sticking, and the ice crystals simply bounce off the surface values that were measured with a thermocouple located at the airfoil without accreting. The maximum limit is described as the limit above leading edge. Some disagreement in temperatures were observed late which there is too much melt and the impinging ice and water mixture in testing for some test cases. Colder regions downstream of the washes away without accreting. Currie [8, 9] suggests that in between leading edge, due to a reduced temperature recovery factor, likely these limits there exists a plateau region where aggressive growth is conducted heat away from the airfoil leading edge, producing possible. The sticking efficiency data as seen in Fig. 10c supports the sub-freezing thermocouple readings late into some tests. This colder existence of a minimum melt ratio limit and plateau region. However, region downstream of the leading edge aided initial ice growth in the the maximum melt ratio limit was not reached for these tests. The local cold region, which grew forward towards the leading edge. This maximum melt ratio limit requires that the liquid water does not secondary front then supported further ice growth at the leading edge.

supercool. In the present experiments it is possible that a portion of the This dynamic helps explain the observed time-dependent leading edge liquid water was supercooling. The 𝜂 range of the plateau as ெோ ice accretion rates. These two-dimensional dynamics point out the measured in this paper (~0.3 to 0.65) exceeds Currie’s [8, 9] melt ratio shortcomings of the thermodynamic model, as it focuses on the range (~0.10 to 0.25), and is more in line with the 𝜂 range that stagnation point, however the model did aid in understanding complex ெோ Baumert [12] reported (~0.2 to 0.6). It should be noted, however, that icing dynamics.

𝜂 values reported in this paper are corrected for several uncertainties ெோ [2, 3], and are calculated differently than by Currie [8,9] and by This paper calculated sticking efficiency, and supports the existence of Baumert [12]. Also, Baumert likely saw supercooled icing with 𝜂 ெோ a minimum melt ratio limit and plateau region. However, the values at or near unity, similar to what was likely observed in this maximum melt ratio limit was not reached for the conditions run. The work. The height of the icing plateau in this paper was measured to be plateau region range from 𝜂 ~0.3 to 0.65, and perhaps even to higher ெோ lower ( 𝜂 = 0.2) than both Currie ( 𝜂 = 0.3 to 0.5) and Baumert ௦௧௜௖௞ ௦௧௜௖௞ melt ratio values. This 𝜂 range is in agreement with Baumert [12] ெோ ( 𝜂 = 0.3 to 0.4). Different airfoil geometries are likely responsible ௦௧௜௖௞ and extends wider than Currie [8, 9]. The plateau region reached a for part of this difference. In addition, while the 𝜂 values reported ௦௧௜௖௞ 𝜂 value of 0.2, which was lower than other reported 𝜂 values.

௦௧௜௖௞ ௦௧௜௖௞ by Baumert also used a NACA 0012 airfoil (about twice the chord as This difference can in part be explained by the higher velocities that the NACA 0012 used in this paper), Baumert’s air flow velocity were run, which this paper showed that the plateau became shallower (40 m/s) was significantly lower than what was utilized in this testing at higher velocities. The higher velocities also pushed the accretion (84 m/s, 133 m/s and 182 m/s). Part of the discrepancy in 𝜂 value ௦௧௜௖௞ range further towards higher 𝜂 values, likely because higher ெோ can be attributed with the finding in this paper that the value of 𝜂 ௦௧௜௖௞ velocities resulted in greater erosion for more glaciated clouds.

increases with decreasing velocity. This correlation between velocity and sticking efficiency is “possibly” in agreement with Currie [9]. The Finally, the ice accretion leading edge angle was investigated for the finding that total pressure has little effect on 𝜂 is in agreement with ௦௧௜௖௞ icing tests. A key finding is that an increase in φ is accompanied by an the work by both Currie and Baumert.

increase in ice growth rate, and therefore a higher 𝜂 value. This ௦௧௜௖௞ angle played in important role when downstream ice grew forward This work examined leading edge ice accretion, and many of the toward the leading edge, which increased φ , which ultimately findings associated with φ are in agreement with Baumert et al. [12].

increased the leading edge growth rate. In addition, as 𝜂 decreased, ெோ These findings include that a decrease in φ is accompanied by a φ decreased as well, likely due to the higher erosion that occurred at decrease in ice growth rate, and therefore a lower 𝜂 value. In ௦௧௜௖௞ lower 𝜂 values. These ice accretion angle findings are in agreement ெோ addition, as 𝜂 decreased, φ decreased as well, likely due to erosion ெோ with Baumert [12].

Page 15 of 17 Structures , SAE International, Prague, CZ, 2015, doi:

References

10.4271/2015-01-2107.

14. Trontin, P., Knotagiannis, A., Blanchard, G., and Villedieu, P., 1. Mason, J. G., Strapp, J. W., and Chow, P., “The Ice Particle “Description and assessment of the new ONERA 2D icing suite Threat to Engines in Flight,” 44th AIAA Aerospace Sciences IGLOO2D.,” 9th AIAA Atmospheric and Space Environments Meeting and Exhibit , AIAA, Reno, NV, 2006, doi: Conference , AIAA, Denver, CO, 2017, doi: 10.2514/6.2017- 10.2514/6.2006-206.

3417.

2. Struk, P. M., Ratvasky, T. P., Bencic, T., Van Zante, J. F., King, 15. Messinger, B. L., “Equilibrium Temperature of an Unheated Icing M. C., Tsao, J.-C., and Bartkus, T. P., “An Initial Study of the Surface as a Function of Air Speed,” Journal of Aeronautical Fundamentals of Ice Crystal Icing Physics in the NASA Science 20: 29-42, 1953, doi: 10.2514/8.2520.

Propulsion Systems Laboratory,” 9th AIAA Atmospheric and 16. Tsao, J-C., Struk, P. M., and Oliver, M. J., “Possible Mechanisms Space Environments Conference , AIAA, Denver, CO, 2017, doi: for Turbofan Engine Ice Crystal Icing at High Altitude,” 6th AIAA 10.2514/6.2017-4242.

Atmospheric and Space Environments Conference , AIAA, 3. Struk, P., Agui, J., Bartkus, T., and Tsao, J-C., “Ice-crystal icing Atlanta, GA, 2014, doi: 10.2514/6.2014-3044.

accretion studies at the NASA Propulsion Systems Laboratory,” 17. Bartkus, T. P., Struk, P. M., and Tsao, J-C., "Evaluation of a SAE 2019 International Conference on Icing of Aircraft, Engines, Thermodynamic Ice Crystal Icing Model Using Experimental Ice and Structures , SAE International, Minneapolis, MN, 2019 Accretion Data", 2018 AIAA Atmospheric and Space (submitted for publication).

Environments Conference , AIAA, Atlanta, GA, 2018, doi: 4. Struk, P. M., King, M. C., Bartkus, T. P., Tsao, J-C., Fuleki, D., 10.2514/6.2018-4129.

Neuteboom, M., and Chalmers, J. L., “Ice Crystal Icing Physics 18. Rigby, D. L., Struk, P. M., and Bidwell, C. S., “Simulation of Study Using a NACA 0012 Airfoil at the National Research Fluid Flow and Collection Efficiency for an SEA Multi-Element Council of Canada’s Research Altitude Test Facility,” 2018 AIAA Probe,” 6th AIAA Atmospheric and Space Environments Atmospheric and Space Environments Conference , AIAA, Conference , AIAA, Atlanta, GA, 2014, doi: Atlanta, GA, 2018, doi:10.2514/6.2018-4224.

10.2514/6.2014 2752.

5. Struk, P. M., Bartkus, T. P., Tsao, J. C., Currie, T., and Fuleki, 19. Bartkus, T. P., Struk, P. M., and Tsao, J-C., “Comparisons of D., “Ice Accretion Measurements on an Airfoil and Wedge in Mixed-Phase Icing Cloud Simulations with Experiments Mixed-Phase Conditions,” SAE 2015 International Conference Conducted at the NASA Propulsion Systems Laboratory,” 9th on Icing of Aircraft, Engines, and Structures , SAE International, AIAA Atmospheric and Space Environments Conference , AIAA, Prague, CZ, 2015, doi: 10.4271/2015-01-2116.

Denver, CO, 2017, doi: 10.2514/6.2017-4243.

6. Struk, P. M., Broeren, A. P., Tsao, J-C., Vargas, M., Wright, W.

20. Bartkus, T. P., Struk, P. M., Tsao, J. C., and Van Zante, J. F., B., Currie, T., Knezevici, D., and Fuleki, D., “Fundamental Ice “Numerical Analysis of Mixed-Phase Icing Cloud Simulations in Crystal Accretion Physics Studies,” SAE 2011 International the NASA Propulsion Systems Laboratory,” 8th AIAA Conference on Aircraft and Engine Icing and Ground Deicing , Atmospheric and Space Environments Conference , AIAA, NASA/TM-2012-217429, 2011, doi: 10.4271/2011-38-0018.

Washington D.C., 2016, doi: 10.2514/6.2016-3739.

7. Currie, T. C., Struk, P. M., Tsao, J., Fuleki, D., and Knezevici, D.

21. Bartkus, T. P., Struk, P. M., and Tsao, J. C., “Development of a C. “Fundamental Study of Mixed-Phase Icing with Application to Coupled Air and Particle Thermal Model for Engine Icing Test Ice Crystal Accretion in Aircraft Jet Engines,” 4th Atmospheric Facilities,” SAE International Journal of Aerospace, Vol. 8, No.

and Space Environments Conference , AIAA, New Orleans, LA, 1, 2015, pp. 15-32, doi: 10.4271/2015-01-2155.

2012, doi: 10.2514/6.2012-3035.

22. Agui, J. H., Struk, P. M., and Bartkus, T. P., “Total Temperature 8. Currie, T. C., Fuleki, D., Knezevici, D. C., and MacLeod, J. D.

Measurements Using a Rearward Facing Probe in Supercooled “Altitude Scaling of Ice Crystal Accretion,” 5th AIAA Liquid Droplet and Ice Crystal Clouds,” 2018 AIAA Atmospheric Atmospheric and Space Environments Conference , AIAA, San and Space Environments Conference , AIAA, Atlanta, GA, 2018, Diego, CA, 2013, doi: 10.2514/6.2013-2677.

doi:10.2514/6.2018-3970.

9. Currie, T. C., Fuleki, D., and Mahallati, A. “Experimental Studies 23. Agui, J., Struk, P., and Bartkus, T., “Total Temperature of Mixed-Phase Sticking Efficiency for Ice Crystal Accretion in Measurements in Icing Cloud Flows using a Rearward Facing Jet Engines,” 6th AIAA Atmospheric and Space Environments Prove,” SAE 2019 International Conference on Icing of Aircraft, Conference , AIAA, Atlanta, GA, 2014, doi: Engines, and Structures , SAE International, Minneapolis, MN, 10.2514/6.2014-3049.

2019 (submitted for publication).

10. Knezevici, D. C., Fuleki, D., Currie, T. C., Galeote, B., Chalmers, 24. King, M. C., Manin, J., Van Zante, J. F., Timko, E. N., and Struk, J., and MacLeod, J. D. “Particle Size Effects on Ice Crystal P. M., “Particle Size Calibration Testing in the NASA Propulsion Accretion - Part II,” 5th AIAA Atmospheric and Space System Laboratory,” 2018 AIAA Atmospheric and Space Environments Conference , AIAA, San Diego, CA, 2013, doi: Environments Conference , AIAA, Atlanta, GA, 2018, doi: 10.2514/6.2013-2676.

10.2514/6.2018-3971.

11. Knezevici, D. C., Fuleki, D., Currie, T. C., and MacLeod, J. D.

25. Schlichting, H., “Boundary-Layer Theory, Seventh Edition,” “Particle Size Effects on Ice Crystal Accretion,” 4th AIAA (New York, McGraw-Hill, 1979), pg. 335, 714, doi: Atmospheric and Space Environments Conference , AIAA, New 10.1002/zamm.19800600419.

Orleans, LA, 2012, doi: 10.2514/6.2012-3039.

12. Baumert, A., Bansmer, S., Trontin, P., and Villedieu, P.,

“Experimental and numerical investigation on aircraft icing at Contact Information

mixed phase conditions,” International Journal of Heat and Mass Transfer 123: 957-978, 2018, doi: Tadas P. Bartkus, Ph.D.

10.1016/j.ijheatmasstransfer.2018.02.008.

Work phone: (216) 433-6915 13. Currie, T., Fuleki, D., and Davison, C., “Simulation of Ice Particle E-mail: tadas.p.bartkus@nasa.gov Melting in the NRCC RATFac Mixed-Phase Icing Tunnel,” SAE Affiliation: Ohio Aerospace Institute 2015 International Conference on Icing of Aircraft, Engines, and Page 16 of 17 RH relative humidity (%)

Acknowledgments

t time (s) ̇ 𝒕 ice thickness growth rate, general (mm/s) The authors wish to acknowledge the Advanced Aircraft Icing (AAI) sub -project of the NASA Advanced Air Transport Technology Project ̇ 𝒕 ice thickness growth rate, freeze-dominated (mm/s) 𝒇𝒓𝒆𝒆𝒛𝒆 (AATT) , under NASA's Advanced Air Vehicles Program (AAVP), for ̇ 𝒕 ice thickness growth rate, melt-dominated (mm/s) 𝒎𝒆𝒍𝒕 financial support for this work .

T temperature (°C or K) 𝑻 temperature of freezing water (0 °C or 273.15 K) 𝒇𝒓𝒆𝒆𝒛𝒆

Nomenclature

𝑻 temperature as measured by a thermocouple (°C or K) 𝑻 / 𝑪 Twb wet-bulb temperature (°C or K) Cp specific heat capacity (J/kg/K) TWC total water content, sum of liquid and ice water contents F loss,total fraction of impinging mass flux that does not accrete (g/m³) (dimensionless) U air velocity (m/s) F fraction of impinging mass flux that does not accrete due loss,spl/rb to splash and runback (dimensionless)  0 collection efficiency at the stagnation line (dimensionless) h c convective heat transfer coefficient (W/m²/K) h mass transfer coefficient (kg/m²/s) η stick fraction of impinging mass flux retained (sticks) on the m surface (dimensionless) j positive integers in the summation of the heat equation infinite series solution (dimensionless)  density k thermal conductivity (W/m/K) φ leading edge ice accretion angle L thickness of airfoil shell (m) L f latent heat of fusion (freezing or melting) (J/g) Subscripts L latent heat of vaporization (J/g) v 0 total conditions 𝜼 melt ratio, ratio of liquid water content to total of liquid 𝑴𝑹 i initial state + ice water content (dimensionless) ice ice MVD median volumetric diameter, in reference to the particle PL plenum of tunnel size spray distribution (microns) s static conditions m surface melting fraction of ice water at stagnation (dimensionless) surf icing surface " 𝒎 ̇ evaporative mass flux (kg/m²/s) targ target 𝒆 " wall wall, airfoil surface 𝒎 ̇ mass flux of impinging solid ice water (kg/m²/s) 𝒊𝒎𝒑 , 𝒊𝒄𝒆 " 𝒎 ̇ mass flux of impinging liquid water time (kg/m²/s) 𝒊𝒎𝒑 , 𝒍𝒊𝒒 " ̇ 𝒎 total impinging water mass flux (ice and liquid) (kg/m²/s) 𝒊𝒎𝒑 , 𝒕𝒐𝒕𝒂𝒍

Abbreviations

n surface freezing fraction of liquid water at stagnation (dimensionless) Escort number Esc # n loss fractional mass loss of ice due to bounce and erosion (dimensionless) NASA National Aeronautics and 𝒏 ෝ in the normal direction; distance into the airfoil leading Space Administration edge wall, in the tunnel axial direction (m) p pressure (Pa) PSL Propulsion Systems p v , s saturation vapor pressure of water in atmosphere, at static Laboratory temperature (Pa) TCS Test Condition Series p v , surf saturation vapor pressure of water over icing surface (Pa) Pr Prandtl number (dimensionless) " T/C Thermocouple 𝒒 conductive heat transfer surface flux (W/m²) 𝒄𝒐𝒏𝒅 " 𝒒 ഥ conductive heat transfer surface flux, averaged over a 𝒄𝒐𝒏𝒅 T/C_03 Thermocouple at leading period of time (W/m ) edge " 𝒒 convective heat transfer surface flux (W/m²) 𝒄𝒐𝒏𝒗 " 𝒒 evaporative heat transfer surface flux (W/m²) T/C_06 Thermocouple just aft of 𝒆𝒗𝒂𝒑 " leading edge 𝒒 latent heat of fusion surface energy for freeze-dominated 𝒇𝒓𝒆𝒆𝒛𝒆 icing (W/m²) " 𝒒 kinetic energy surface flux (W/m²) 𝒌𝒊𝒏𝒆𝒕𝒊𝒄 " 𝒒 latent heat of fusion surface energy for melt-dominated 𝒎𝒆𝒍𝒕 icing (W/m²) " 𝒒 sensible heat surface flux for frozen ice (W/m²) 𝒔𝒆𝒏𝒔 , 𝒊𝒄𝒆 Page 17 of 17