# Development of a coupled air and particle thermal model for engine icing test facilities.

ABSTRACTThis paper describes a numerical model that simulates the thermal interaction between ice particles, water droplets, and the flowing air applicable during icing wind tunnel tests where there is significant phase-change of the cloud. It has been previously observed that test conditions, most notably temperature and humidity, change when the icing cloud is activated. It is hypothesized that the ice particles and water droplets thermally interact with the flowing air causing the air temperature and humidity to change by the time it reaches the test section. Unlike previous models where the air and particles are uncoupled, this model attempts to explain the observed changes in test conditions by coupling the conservation of mass and energy equations. The model is compared to measurements taken during wind tunnel tests simulating ice-crystal and mixed-phase icing that relate to ice accretions within turbofan engines. The model simulates trends that were experimentally observed, but does not fully explain the measured values. Some possible explanations for this discrepancy are offered. This model, written in MATLAB, is based on fundamental conservation laws and empirical correlations.

CITATION: Bartkus, T, Struk, P., and Tsao, J., "Development of a Coupled Air and Particle Thermal Model for Engine Icing Test Facilities," SAE Int. J. Aerosp. 8(1):2015, doi: 10.4271/2015-01-2155.

INTRODUCTION

Several jet engine power-loss events of commercial aircraft at high altitude have been reported since the 1990's. Mason et al. [1] discussed how power-loss can result from ice crystals entering the engine core, partially melting from the warm engine airflow and refreezing on internal components. According to Mason et al. [2], the accumulation of ice in the engine, from these ice-crystals, can cause engine surge, stall, flameout, and even damage to the blades and vanes from ice shedding.

Due to these numerous power-loss events, potential ice accretion within the engine due to the ingestion of ice crystals is being investigated. To better understand this phenomenon and determine the physical mechanism of engine ice accretion, fundamental tests have been collaboratively conducted by NASA Glenn Research Center and the National Research Council of Canada (NRC) ([3, 4]). NRC's Research Altitude Test Facility (RATFac) is capable of generating ice-crystal and mixed-phase cloud conditions. While conducting icing tests, the total air temperature measured at the test section was observed to decrease when the icing cloud was activated. In some cases, the temperature decrease was considerable, by as much as several degrees Celsius, and was accompanied by a measurable increase in water vapor ([5, 6]). It is hypothesized that the ice particles and water droplets thermally interacted with the flowing air causing the air temperature and humidity to change by the time it reached the test section. Since a difference of a few degrees can mark the difference between an icing event and a non-icing event, it is important to understand the cause of this temperature change so that accurate test conditions can be specified. To that end, the focus of this paper is to model the observed changes in air temperature and water vapor content due to the air-particle enthalpy and mass transfer exchange. In this paper, the term 'particle' is generically used to refer to liquid water, ice or mixed phase water.

Existing thermal models ([7, 8, 9, 10, 11, 12]) show little effect due to coupling of the ice/water particles with the flowing air, or do not couple, approximating the air mass as an infinite thermal reservoir with unchanging properties. The model presented in this paper attempts to explain the observed changes in test conditions by coupling the conservation of mass and energy equations for both the ice/water particle and flowing air mass. The model simulates both ice particle melting and droplet freezing, including spontaneous latent heat release during freezing for supercooled droplets. In the simulation, humidity changes are due to ice particle sublimation/deposition and water droplet evaporation/condensation. The associated latent heating/cooling and mass transfer change the ice water content (IWC) and liquid water content (LWC) from the moment of injection to the time the flowing mass reaches the test section. The conservation of momentum is employed to account for particle velocity. In addition, to more accurately simulate spray injection, the model incorporates a distribution of particle sizes as an initial input condition. As a result, the model predicts changes in air temperature, humidity, IWC and LWC, as well as the particle phase, size, temperature and velocity at the test section.

THERMAL MODEL DESCRIPTION

The thermal model simulates an icing wind tunnel by calculating particle and air properties using the conservation of mass, momentum and energy equations. The ice, water, and humid air mass in an icing tunnel is broken into fundamental control volume units. For a cloud containing uniformly sized particles, the fundamental unit is composed of a single ice or water particle and a parcel of humid air. For an icing cloud containing a spectrum of particle sizes, the fundamental unit contains the particle distribution and corresponding air parcel. The particle(s) and the parcel of air are each treated as separate control volumes where energy and mass transfer between each other. The energy and mass in each fundamental control volume are maintained as constant. Details of the control volumes are described later. The model simulates particles on the icing tunnel centerline, providing thermodynamic conditions at the test section.

Assumptions

The thermal model operates under the following assumptions:

1. Air and particle flow are one dimensional where the geometry of the icing tunnel is constant in the direction of flow.

2. Air and water vapor are treated as ideal gases.

3. Air is well mixed, meaning that temperature, water vapor content and all thermodynamic properties are homogeneous within the air control volume at any time (or axial location).

4. The fundamental control volume is adiabatic and mass is maintained as constant.

5. Supersaturated water vapor content in air does not occur.

6. Particles are evenly spaced within the air from the injection point to the test section.

7. All particles are perfectly spherical.

8. Particle aggregation / coalescences and breakup through collision are negligible.

9. Particle size distribution is characterized as discrete diameters.

10. Particles are injected in the direction of the flow and remain entrained with negligible gravity affects.

11. Temperature is uniform within the particle.

12. Water droplets can be supercooled to a homogenous crystallization temperature that is a function of particle diameter.

13. The latent heat release stage during supercooled freezing occurs instantaneously.

14. Mixed phase particles are spatially homogeneous in water/ice content. This means that evaporation and sublimation occurs from a mixed phase particle surface and the rate is determined by the water/ice content ratio. This applies to the reverse process of condensation and deposition as well.

15. Evaporation, sublimation, condensation and deposition occur at the particle surface at particle temperature.

16. Air velocity and pressure do not change due to heat or mass transfer.

17. Heat transfer at the tunnel walls is negligible.

18. The flow of particles and air is a continuous stream. This means that while following a particular set of particles within a reference control volume, faster moving particles and air that exit the reference control volume are replaced by thermodynamically identical particles and air from a neighboring control volume that is upstream of the reference control volume. The reference control volume is centered on the slowest moving particle in distribution of particles.

To validate Assumption 11, that temperature is uniform within the particle, the Biot number, Bi, is calculated. The Biot number is a ratio of convective heat transfer rate across the particle boundary to the rate of conductive heat transfer within the particle. A Bi value that is below a critical value of 0.1 verifies the use of the lumped capacitance method where particle temperature is uniform within the particle.

Bi = [[[hL.sub.c]]/[[k.sub.p]]] = [hd]/[[6k.sub.p]] < 0.1 for spheres (1)

In the numerator, h is the convective heat transfer coefficient, L is the critical length which equals 1/6th of the particle diameter, d, for a sphere. The particle's thermal conductivity is represented by [k.sub.p]. The largest Bi values will occur with large particles (d = ~200 microns) and high slip velocities (U = ~30 m/s). Using these extreme values along with other appropriate properties, Bi = ~0.015, much less than the critical value needed to validate use of lumped capacitance.

Thermal Model Algorithm

The thermal model was written in MATLAB version R2014a. The model solves the differential equations of mass, momentum and energy transfer. All conservation equations presented in the following model formulation section are solved using MATLAB's built-in ODE45 solver. The system of differential equations is solved by first inputting the initial conditions. Time marching methods are employed by the MATLAB solver such that each iterative solution is solved with numerical relative and absolute tolerances of [10.sup.-12]. Physical accuracy of the solution is dependent on the accuracy of all the property values input into the model. For all simulations presented in this paper, the energy transferred between the air and particle(s) is balanced to the order of [10.sup.-4]. The mass transfer between the air and particle(s) initial and final states is balanced to the order of [10.sup.-15].

Thermal Model Equations: Single Particle Formulation

For an icing cloud of uniform particle size, the fundamental control volume consists of just a single particle control volume and air control volume. The mass of the air control volume, [m.sub.air], is determined by the icing conditions and is the following equation.

[m.sub.air] = [[[[pi]d.sup.3]]/[6]] [[[[rho].sub.p][[rho].sub.air]]/[LWC]] (2)

In Eq. (2), d is the particle diameter, [[rho].sub.p] is the particle density, and [[rho].sub.air] is the air density. The denominator, LWC, can be interchanged with LWC for an ice particle cloud. LWC and LWC represent water content density with units of mass of water per volume of air. The conservation of mass, momentum and energy equations are applied to the two interacting control volumes.

Conservation of Energy

The energy balance for a particle (ice particle, water droplet or mixed phase) thermally interacting with a parcel of air is the balance between the internal energy of the particle and the heat flux across the particle boundary. More specifically, the change in the particle's enthalpy, [H.sub.p], with respect to time, t, is due to the convective heat transfer, [q.sub.conv], and latent energy due to mass transfer, [q.sub.latent], and is given by the following general expression.

[[[partial derivative][H.sub.p]]/[partial derivative]t] = [q.sub.conv] + [q.sub.latent] (3)

Dalton [13] is credited for starting the empirical hydrodynamic approach to the evaporation problem. He stated that the rate of mass loss from a water surface is proportional to the difference in water vapor pressure between the particle surface and the ambient air. He also stated that this rate is proportional to the wind speed. Applying Dalton's formulation to the latent energy term, as well as expanding on the internal energy and convective energy terms of Eq. (3). gives the following expression for an ice particle undergoing sublimation.

[[partial derivative].sub.p][C.sub.p] [[[[pi]d.sup.3]]/[6]] [[partial derivative][T.sub.p]]/[[[partial derivative]t]] = [[pi]d.sup.2]h([T.sub.air] - [T.sub.p]) + [[pi]h.sup.2][h.sub.m][[rho].sub.air][L.ub.subl]([[omega].sub.air] - [[omega].sub.p]) (4)

In Eq. (4). C is the specific heat capacity of the particle, [T.sub.p] is the particle temperature, [T.sub.air] is the ambient air temperature, h is the heat transfer coefficient, [h.sub.m ] is the mass transfer coefficient, [L.sub.subl], is the latent heat of sublimation, [[omega].sub.air] is the water vapor mass fraction in the ambient air, and [[omega].sub.p] is the water vapor mass fraction at the particle surface. The water vapor mass fractions are defined as the following.

[[omega].sub.air] = [[[MW.sub.water]]/[[MW.sub.air]]] [[[P.sub.wv,air]]/[[P.sub.air]]] (5)

[[omega].sub.surf] = [[[MW.sub.water]]/[[MW.sub.air]]] [[[P.sub.wv,surf]]/[[P.sub.air]]] (6)

In Eqns. (5) and (6), MW is molecular weight, [P.sub.air] is air pressure, [P.sub.wv,air] is the ambient vapor pressure, and [P.sub.wv,surf] is vapor pressure at the particle surface (saturated). Eq. (4) applies for a fully frozen ice particle prior to any melting. For an ice particle undergoing melting, the ice/water temperature remains constant and the rate of the phase change is given by the following expression.

[[rho].sub.p][L.sub.melt] [[[[pi]d.sup.3]]/[6]] [[[[partial derivative]][eta].sub.p]/[partial derivative]t] [[pi]d.sup.2]h([T.sub.air] - [T.sub.p]) +

[[pi]d.sup.2][h.sub.m][[rho].sub.air][L.sub.subl/evap]([[omega].sub.air] - [[omega].sub.p]) (7)

In Eq. (7), [L.sub.melt] is the latent heat of melting, [L.sub.subl/evap] is the proportional mix of the latent heat of sublimation and evaporation, and [[eta].sub.p] is the particle melt fraction (or melt ratio). Melt fraction is the ratio of liquid water mass to the total water/ice content (0 is solid ice and 1 is fully liquid). The melt ratio term is used for both melting and freezing cases. The value of [L.sub.subl/evap] is associated with the melt fraction and is given by the following expression.

[L.sub.subl/evap] = (1 - [[eta].sub.p])[L.sub.subl] + ([[eta].sub.p])[L.sub.evap] (8)

For a water droplet undergoing evaporation, the formulation is similar to Eq. (4) except the latent heat of evaporation variable, [L.sub.evap], replaces [L.sub.subl].

[[partial derivative].sub.p][C.sub.p] [[[[pi]d.sup.3]]/[6]] [[partial derivative][T.sub.p]]/[[[partial derivative]t]] = [[pi]d.sup.2]h([T.sub.air] - [T.sub.p]) + [[pi]d.sup.2][h.sub.m][[rho].sub.pair][L.ub.evap]([[omega].sub.air] - [[omega].sub.p]) (9)

Equations (4), (7) and (9) are explicitly defined for ice sublimation, melting and droplet evaporation. These particle energy equations can be used for the reverse thermal processes of ice deposition, freezing and water condensation respectively. The absolute values of the corresponding latent heats of phase change (sublimation : deposition, melt : freeze, and evaporation : condensation) are equal and the direction of the latent energy due to mass transfer is determined by the water vapor mass fraction difference.

The non-dimensional Nusselt number, Nu, and Sherwood number, Sh, are used to determine the heat and mass transfer coefficients in the three particle energy equations and are defined as follows.

Nu = [[hd]/[[k.sub.air]]] (10)

Sh = [[[h.sub.m]d]/[[D.sub.ab]]] (11)

In the above equations, [k.sub.air] is the thermal conductivity of air, and [D.sub.ab] is the diffusivity of water vapor in air. The empirical heat transfer expressions for flow over a sphere [14] is as follows:

Nu = 2 + 0.6[Re.sup.1/2][Pr.sup.1/3] (12)

Using the heat and mass transfer analogy, the Sherwood number is the following:

Sh = 2 + 0.6[Re.sup.1/2][Sc.sup.1/3] (13)

In the above equations, the Reynolds number, Re, Prandtl number, Pr, and Schmidt number, Sc, are defined as follows:

Re = [[[[rho].sub.air]\[v.sub.air] - [v.sub.p]\d]/[[mu].sub.air]] (14)

Pr = [[[C.sub.air][mu].sub.air]]/[[k.sub.air]] (15)

Sc = [[[[mu].sub.air]]/[[[rho].sub.air][D.sub.ab]]] (16)

In the above equations, [C.sub.air] is the specific heat capacity of air, [[mu].sub.air] is the dynamic viscosity and the difference between the air velocity, [v.sub.air], and the particle velocity, [v.sub.p], is also known as the slip velocity.

The energy that transfers across the particle control volume is balanced by an equivalent and opposite amount of air enthalpy change. The change in the air parcel's enthalpy, [H.sub.air], with respect to time is due to the convective heat transfer and the sensible energy change of the water vapor mass that has changed phase, [q.sub.wv,sens]. Water vapor mass that has changed phase can be vapor that has been liberated from the particle by evaporation/sublimation, or it can refer to the water vapor in the moist air that has condensed/deposited onto the particle. The conservation of energy for the air mass is given by the following general expression.

[[partial derivative][H.sub.air]/[[partial derivative]t]] = [q.sub.conv] + [q.sub.wv,sens] (17)

The sensible energy term in Eq. (17) is typically small compared to the convective heat transfer term, but is necessary to account for all energy being transferred. For example, a cool water droplet in warm air will evaporate, producing water vapor whose temperature is at the cool water droplet temperature. The vaporized mass will increase in temperature as it mixes with warm air. This thermal mixing situation is accounted for by the sensible energy term in Eq. (17). Expanding the above equation gives the following energy balance for air:

[MATHEMATICAL EXPRESSION NOT REPRODUCIBLE IN ASCII] (18)

Where [m.sub.wv] is the mass of water vapor that has changed phase, [C.sub.wv] is the specific heat capacity of vapor, and [T.sub.wv] is the vapor temperature.

Conservation of Mass

The change in particle mass, [m.sub.p], due to phase change (evaporation, sublimation, condensation or deposition) for a spherical particle is given by the following mass balance equation.

[[[partial derivative]m.sub.p]]/[[partial derivative]t]] = [[rho].sub.p] [[[partial derivative]]/[[partial derivative]t]] ([[[pi]d.sup.3]/[6]]) = [[pi]d.sup.2][h.sub.m][[rho].sub.air] ([[omega].sub.air - [[omega].sub.p]) (19)

To balance the overall system mass, the change in mass for the air mass is equal to but opposite the change in particle mass and is given by the following expression.

[[[partial derivative]m.sub.air]]/[[partial derivative]t]] = - [[[[partial derivative]m.sub.p]]/[[partial derivative]t]] = [[pi]d.sup.2][h.sub.m][[rho].sub.air] ([[omega].sub.p - [[omega].sub.air]) (20)

Situations can occur where the air is at vapor saturation, yet the driving force for evaporation (or sublimation) still exists such that co > [[omega].sub.air]. For example, a warm water droplet in saturated cold air is a situation where the vapor gradient would favor evaporation. This model assumes that supersaturation of the air does not exist (Assumption 5). A limit function, [phi], is appended to the mass transfer term such that when the air is near saturation, it will regulate the mass transfer rate asymptotically at RH = 100%. The limit function reduces mass transfer rate at a level near saturation defined by the variable [RH.sub.limit] For the thermal model, mass transfer regulation was arbitrarily chosen to begin at [RH.sub.limit] = 99%.

[phi] = [[100-RH]/[100-[RH.sub.limit]]] (21)

The right hand term of the conservation of mass equations (Eq. 091 and Eq. (20)) is multiplied by this limit function when relative humidity is at or above [RH.sub.limit] and [[omega].sub.p] > [[omega].sub.air]. Similarly, the latent energy due to mass transfer term in the particle conservation of energy equations (Eq. (4), Eq. (7), and Eq. (9)) is multiplied by this limit function for the same conditions.

Conservation of Momentum

Unlike the conservation of mass equations where the air and particle together constitute a closed system, the conservation of momentum equation is solved only in reference to the particle. The forces exerted onto the particle are the forces of drag, [F.sub.drag], and gravity, [F.sub.g].

F = [F.sub.drag] + [F.sub.g] [right arrow] [m.sub.p]a = [[rho] sub.p] [[[pi]d.sup.3]/[6]] [[partial derivative]v.sub.p]/[[partial derivative]t] (22)

In the above equation, a is the particle's acceleration due to the sum of all forces acting on the particle, F. The force of gravity is negligible and can be assumed to be zero.

[F.sub.g] = 0 (23)

The drag force can be described by the following equation.

[F.sub.drag] = 1/2 [[rho].sub.air][U.sub.2][AC.sub.D] = 1/2 [[rho].sub.air] [([v.sub.air] - [v.sub.p]).sup.2] [[pi]d.sup.2][C.sub.D] (24)

In Eg. (24). U is the slip velocity and A is the projected area of the spherical particle. Many correlations exist for the coefficient of drag for flow past a sphere. Morrison [15] provides a correlation that spans a wide range of Reynolds numbers and is given below.

[MATHEMATICAL EXPRESSION NOT REPRODUCIBLE IN ASCII] (25)

Substituting and rearranging terms in Eq. (22) provides a simple differential equation for the conservation of momentum.

[[[[partial derivative]v.sub.p]]/[[partial derivative]t]] = 3/4 [[[[rho].sub.air][C.sub.D][([v.sub.air] - [v.sub.p]).sup.2]]/[[[rho].sub.p]d]] (26)

Energy Balance Equations with a Heat Source

The axial air temperature in the RATFac icing tunnel facility at NRC is not constant even in the absence of an icing cloud. A detailed description of the icing facility is given by Knezevici et al. [16]. The air temperature increases from the point of injection to test section. Without the presence of a cloud, cold air is injected into the flow stream and is mixed with warm (and moist) ambient air. The cold air jet and entrained warm air are drawn by a vacuum downstream of the test section. A computational fluid model (CFD) was previously written that approximated the air temperature and flow field profiles along the icing tunnel (these internal calculations are currently unpublished). These profiles were implemented into this paper's model. Figure 1 shows approximate values for the normalized air temperature difference between the test section air temperature, [T.sub.test], and the injection air temperature, [T.sub.inj] at the tunnel centerline, (T -[T.sub.inj] )/([T.sub.test] - [T.sub.inj]).

Figure 1 shows the computational results where the test section air speed is Mach = 0.2 with a total pressure of 43644 Pa. The distance between cold air injection and test section is approximately 2 m. While the values obtained from the fluid model are specific to the conditions run, for the purposes of this paper, the general normalized temperature profile was used as an approximation for different injection and test section temperatures.

A parcel of air - without any particles - heated by a source, [q.sub.source], is given by the following expression.

[[[partial derivative][H.sub.air]]/[[partial derivative]t]] = [m.sub.air][C.sub.air] [[[partial derivative][T.sub.air]]/[[partial derivative]t]] = [q.sub.source] (27)

To convert the above equation from a stationary air mass to one that simulates an air mass moving in the x-direction, the following manipulation can be done.

[MATHEMATICAL EXPRESSION NOT REPRODUCIBLE IN ASCII] (28)

In Eq. (28), x is distance and [partial derivative][T.sub.air] /[partial derivative]x can be derived from the computational fluid model that approximated the temperature profile. [partial derivative][T.sub.air]/[partial derivative]x is the derivative of the polynomial fit in Figure 1. The term [partial derivative]x/[partial derivative]t is the velocity of the reference frame, which is the particle velocity. A constant value of [partial derivative][T.sub.air]/[partial derivative]x = 0 was prescribed after the distance where the normalized temperature surpassed a value of 1.

To approximate this changing temperature air profile in the icing tunnel, the source term is added to the air energy equation.

[[[partial derivative][H.sub.air]]/[[partial derivative]t]] = [q.sub.conv] + [q.sub.wv,sens] + [q.sub.source] (29)

Mass Balance Equations with a Mass Source

As was the case with the icing tunnel temperature profile, the humidity profile changes. At the icing facility, warm humid air is entrained by a cold dry air stream and mixes. A parcel of air - without any particles - gaining water vapor mass from a source is given by the following simple expression.

[[[[partial derivative]m.sub.air]]/[[partial derivative]t]] = [[[[partial derivative]m.sub.wv]]/[[partial derivative]t]] = [m.sub.source] (30)

Again, to convert the above equation from a stationary air mass to one that simulates an air mass moving in the x-direction, the following manipulation can be done.

[MATHEMATICAL EXPRESSION NOT REPRODUCIBLE IN ASCII] (31)

A water vapor profile for the icing facility, in the absence of a cloud, was not previously modelled. As an approximation, the normalized profile determined for the air temperature is used analogously for the water vapor mass profile. Therefore, in the above equation, [[partial derivative]m.sub.wv]/[partial derivative]x can be derived from the water vapor mass profile.

A water vapor mass source term, [m.sub.source], is added to the air mass balance equation to approximate the changing humidity along the length of the icing tunnel.

[[[partial derivative]m.sub.air]]/[[partial derivative]t]] = - [[[[partial derivative]m.sub.p]/[[partial derivative]t]] + [m.sub.source] (32)

Flow Field Profile

The air velocity is not constant from the point of air injection to the test section. A flow field was calculated for the icing tunnel in the previously written computational fluid model. A normalized velocity profile at the icing tunnel centerline is given in Figure 2. The expression in Figure 2 is valid only when the injection velocity, [v.sub.inj], equals the test section velocity, [v.sub.test], with values near Mach =0.2.

A simple differential equation can be written to approximate the changing flow field air velocity.

[MATHEMATICAL EXPRESSION NOT REPRODUCIBLE IN ASCII] (33)

In the above equation, [[partial derivative]v.sub.air]/[partial derivative]x is the derivative of the polynomial fit in Figure 2. The term [partial derivative]x/[partial derivative]t is the velocity of the reference frame, which is the particle velocity. A constant value of [[partial derivative]v.sub.air]/[partial derivative]x = 0 was prescribed after the distance where the normalized velocity surpassed a value of 1.

Supercooling Formulation

Freezing of liquid water can occur in two different processes: (1) normal freezing and (2) supercooled freezing. Normal freezing occurs with tap water, for example, which contains mineral impurities. These act as nucleation sites for freezing to occur at the equilibrium freezing temperature. At atmospheric pressures, city water freezes at 273.15 K. Supercooling is the process of cooling a liquid below its equilibrium freezing point without it becoming a solid. De-ionized water or cloud droplets can be supercooled to their homogeneous crystallization temperature (~ -40[degrees]C) before crystallization occurs.

The process of normal freezing in air can be described in 3 stages:

1. Sensible cooling of liquid water until water temperature reaches 273.15 K.

2. Latent freezing where heat is transferred from the particle at a constant temperature of 273.15 K until it has completely frozen.

3. Sensible cooling of the solid ice particle until the temperature is reduced to the wet-bulb temperature.

The process of supercooling, also outlined by Hindmarsh et al. [8], can be described in 4 stages for a particle and follows Figure 3.

1. Sensible cooling of liquid water to below the equilibrium freezing temperature until homogeneous crystallization temperature is reached.

2. Latent heat release of supercooled water where rapid crystallization of a particle occurs, along with an increase in particle temperature. The particle will partially freeze and will have reached its equilibrium freezing temperature of 273.15 K.

3. Latent freezing where heat is transferred from the partially frozen particle at a constant temperature of 273.15 K until it has completely frozen.

4. Sensible cooling of the solid ice particle until the temperature is reduced to the wet-bulb temperature.

An expression for the homogenous crystallization temperature, [T.sub.hc], was derived from data given by Heverly [17] and is given below.

[T.sub.hc] = 7.2015 ln (d) + 214.64 (34)

During the latent heat release stage, the fraction of the particle that crystallizes will depend on the amount of sensible heat lost during supercooling, [H.sub.sens,super]. The melt ratio of the particle resulting after the latent heat release is found using the following formulas.

[MATHEMATICAL EXPRESSION NOT REPRODUCIBLE IN ASCII] (35)

[[eta].sub.p] = 1 - [[[H.sub.sens,super]]/[[m.sub.p][L.sub.melt]]] (36)

Thermal Model Equations: Particle Distribution Formulation

The single particle formulation (uniform particle size cloud) in the previous sections can be used to approximate the thermal interaction between the air and icing cloud using the jet spray's median volume diameter, MVD. An icing cloud, however, contains a distribution of particle sizes. A more accurate icing tunnel simulation accounts for the thermal interactions between all the differently sized particles in an icing cloud. Smaller particles have a larger surface area to volume ratio compared to larger particles and therefore will have faster thermal and mass transfer responses. The cumulative differences will result in a different solution compared to a uniform particle size cloud simulation.

Ice particle [18] and water droplet [19] distributions were used for simulations in this paper. An ice particle distribution with MVD = 45.5 [micro]m (Figure 12, Grinder Config. A of Reference 18) and a water droplet distribution with MVD = 40 [micro]m (Figure 8 of Reference 19) were used. With respect to Reference 19, the water spray jets used at NASA Glenn's Icing Research Tunnel are the same model as the spray jets used at NRC's RATFac. The fraction of the total water mass (or ice mass) was calculated for every particle size bin.

For the full particle distribution formulation, a control volume is prescribed around every particle of every size. The number of particles of a particular diameter is determined by its respective mass fraction distribution and the water content (or ice content). Since IWC and LWC are generally given as grams per cubic meter, the air control volume is the cubic meter of air containing the mass of water/ice given by IWC and LWC. Figure 4 shows an example of the control volumes of differently sized particles and the air mass control volume that contains all of the particles.

Conservation of Energy

An energy balance is done for every particle size, therefore the conservation of energy equation for particle size i is the following.

[MATHEMATICAL EXPRESSION NOT REPRODUCIBLE IN ASCII] (37)

The total number of particle energy equations is equal to the total number of particle size bins in the distribution. Only one air energy equation is needed, but it contains the sum of all the heat transfers and water vapor (that has changed phase from each particle) sensible energy transfers.

[MATHEMATICAL EXPRESSION NOT REPRODUCIBLE IN ASCII] (38)

In the above expression, n is the number of particle size bins and [#.sub.i] is the number of particles in the [i.sup.th] bin.

Conservation of Mass

A mass balance is done for every particle size, therefore the conservation of mass equation for particle size i is the following.

[MATHEMATICAL EXPRESSION NOT REPRODUCIBLE IN ASCII] (39)

The total number of particle mass equations is equal to the total number of particle size bins in the distribution. Only one conservation of air mass equation is needed, but it contains the sum of all the mass transfers across each particle boundary.

[MATHEMATICAL EXPRESSION NOT REPRODUCIBLE IN ASCII] (40)

These last four energy and mass equations can be applied to an ice distribution, a water distribution, or combined ice and water situation.

Thermodynamic Properties

Several air, water, ice and vapor thermo-physical properties are dependent on temperature, pressure and relative humidity. Equations were fit to published thermo-physical data ([14, 21, 22, 23, 24, 25, 26, 27]) and are listed in the Appendix.

PARAMETRIC ANALYSIS

A parametric study was performed to show the effect of several parameters on the particle-air interaction. The conditions of two baseline simulations are depicted in Table 1. The conditions for Baseline 1 represent water droplets interacting with above freezing air, while Baseline 2 represents ice particles interacting with an above freezing air mass. In both cases, the water is assumed to undergo normal freezing and does not supercool. Three parametric simulations - Sim for short - were run where slip velocity, particle temperature and pressure were individually varied with respect to the initial conditions of Baseline 1. Sim 1 [right arrow] 3 of Table 2 correspond to these individually varied conditions with respect to Baseline 1. The value of the changed parameter is depicted in the final column of Table 2. Three more simulations were run where IWC, relative humidity and air temperature were individually varied with respect to Baseline 2 conditions. Sim 4 [right arrow] 6 of Table 2 correspond to the individually varied parameter with respect to Baseline 2. Unless noted, all simulations do not vary in slip velocity, air temperature, and gas water content in the absence of cloud.

Figures 5, 6, 7, 8 show the parametric analysis results where the air and particle interacted for 0.0235 seconds. This represents the elapsed time for a particle to travel 2.0 m at a fixed particle velocity of 85 m/s. Figure 5 shows the transient particle conditions for Baseline 1 ([B.sub.1]) and its corresponding 3 parametric simulations. Figure 6 shows the transient air conditions for Sim [B.sub.1] and its corresponding parametric simulations. Figure 7 and Figure 8 are similar to Figure 5 and Figure 6 except they represent the particle and air transient data with respect to Baseline 2 ([B.sub.2]) and its 3 corresponding parametric simulations.

In the first baseline simulation, a water droplet underwent evaporation in warm air. According to Figure 5a, the droplet temperature decreased rapidly until it reached wet-bulb temperature, [T.sub.wb], of about 273.8 K. Wet bulb temperature is the temperature of a wet adiabatic surface undergoing evaporation (or sublimation if ice). [T.sub.wb] is the resulting temperature when convective heat transfer is balanced by evaporative cooling and is a function of [T.sub.air], RH, and P. Wet bulb temperature can vary slightly as the three aforementioned parameters change with time (i.e. axially). Since the particle temperature remained above freezing throughout the simulation, the melt ratio remained at 1 (fully liquid droplet). Figure 6 shows that by the end of the simulation, the air temperature decreased by 0.34 K and had an increase in relative humidity by 3.5%.

The slip velocity was increased from 5 m/s to 25 m/s in Sim 1. This parametric change increased heat and mass transfer rates. The droplet temperature reached [T.sub.wb] slightly faster than Sim [B.sub.1] due to this increased heat transfer rate, however [T.sub.wb] was very nearly the same as Sim B,. This is because the initial [T.sub.air], RH, and P were identical in both simulations. Compared to Sim [B.sub.1] a greater fraction of the droplet had vaporized due to the increased mass transfer rate, resulting in a greater change in air temperature, [DELTA][T.sub.air] = -0.45 K, and relative humidity [DELTA]RH = 4.7%.

The initial water droplet temperature was decreased from 278.15 K to 273.15 K in Sim 2. Despite the change, the particle and air profiles were very nearly identical to Sim [B.sub.1] The only difference was that the particle temperature increased until it reached wet-bulb temperature. The nearly identical profiles between Sim 2 and Sim [B.sub.1] can be explained by the fact that the air thermal mass dominated the interaction. The air thermal mass was several magnitudes greater than the water droplet thermal mass. Since the initial air conditions were identical between Sim 2 and Sim [B.sub.1], the resulting profiles were consequently nearly identical.

The pressure was reduced from 88,000 Pa to 28,000 Pa in Sim 3. The reduced pressure significantly increased the particle's evaporative cooling. At the lower pressure, the wet-bulb temperature decreased to below freezing temperatures at about 271.5 K. The particle temperature responded by decreasing rapidly to the freezing temperature of 273.15 K, at which point the particle began to crystalize. During the latent freezing stage, the particle remained at constant temperature. The particle became fully crystalized in 0.0128 s and decreased in temperature until it reached the wet-bulb temperature near 271.5 K. The rate of vaporization decreased after complete crystallization (less steep slope in Figure 5b) because when the temperature decreased to around 271.5 K, the vapor mass fraction difference between the surface and ambient decreased, which reduced the mass transport rate. However, a greater amount of the particle vaporized compared to Sim [B.sub.1] due to the lower pressure, which was also accommodated by a greater air temperature decrease, [DELTA][T.sub.air] = -1.2 K. As a result of the increased vapor in the air coupled by the colder air temperature, the relative humidity increased by ARH = 8.6%.

The initial conditions for the second baseline simulations were similar to the first baseline simulation, except in Sim [B.sub.2] the air temperature was slightly warmer and the particle temperature was below freezing. In Sim [B.sub.2], a fully frozen ice particle initially vaporized (sublimated) in warm air. For the given air conditions, the wet-bulb temperature was greater than freezing, as a result the particle temperature increased to melting temperature, fully melted and continued to increase until it reached a wet-bulb temperature of 275.3 K. According to Figure 7b, the rate of vaporization increased after the increase in water droplet temperature due to the corresponding increase in vapor pressure at the particle surface. As can be seen in Figure 8, the air temperature decreased and relative humidity increased, [DELTA][T.sub.air] = -0.6 K, and [DELTA]RH = 3.4%, respectively.

The ice water content was increased from 1 g/[m.sup.3] to 5 g/[m.sup.3] in Sim 4. Similar to Sim [B.sup.2], the ice particle increased to melt temperature, melted and continued to increase in temperature until it reached wet-bulb temperature of 274.8 K. Due to the larger water content, more mass evaporated, which resulted in a larger decrease in air temperature, [DELTA][T.sub.air] = -2.5 K. The lower wet-bulb temperature compared to Sim [B.sub.2] (274.8 K vs 275.3 K) is explained by the lower air temperature. The decrease in air temperature combined with the increase in water vapor explains the high relative humidity increase, [DELTA]RH = 15.4%.

The relative humidity was increased from 50% to 80% in Sim 5. Again, the ice particle melted and increased in temperature to the wet-bulb temperature of 278.2 K. The wet-bulb temperature was higher than the previous two simulations (Sim [B.sub.2] and Sim 5) due to the higher relative humidity. During the sensible ice heating and melting stages, the vapor content was greater in the ambient air than at the cold particle surface, driving deposition and condensation, respectively. This growth in particle mass is realized in Figure 7b. As the liquid droplet increased to a higher temperature, the vapor concentration at the particle surface became greater than in the ambient air, at which point the droplet began to evaporate. The high relative humidity limited the amount of evaporation, which resulted in a small air temperature drop, [DELTA][T.sub.air] = -0.3 K, and small increase in relative humidity, [DELTA]RH = 1.6%.

The air temperature was reduced to below freezing from 280.15 K to 271.15 K in Sim 6. The ice particle dropped in temperature reaching the wet-bulb temperature of 267.7 K, remaining as a solid ice particle throughout the duration of the simulation. The evaporation rate was greater than Sim [B.sub.2] because the overall vapor difference between the ambient and particle surface was greater in Sim 6 for the duration of the simulation. The melting phase of Sim [B.sub.2] limited the rate of vaporization because the vapor concentration difference was small. Sim 6 did not have any particle phase change limiting the vapor concentration difference, leading to a greater evaporation rate for the majority of the simulation time. The decrease in air temperature and increase in relative humidity were [DELTA][T.sub.air] = -0.3 K, and [DELTA]RH = 3.9% respectively.

Two simulations were run to demonstrate the differences between supercooled freezing (Sim SC) and normal freezing (Sim Norm). The conditions for the two simulations were identical except for the style of freezing and are listed in Table 3. The conditions were such that a water droplet vaporized in cold air with low relative humidity. Results are shown in Figures 9 and 10.

The particle temperature in Sim Norm decreased to freezing temperature and after complete crystallization continued to decrease to wet-bulb temperature of around 235 K. The particle temperature in Sim SC followed a different profile. Initially, sensible cooling reduced the temperature to supercooled temperatures. The particle remained fully liquid until the homogeneous crystallization temperature was reached, at which point latent heat release of the particle occurred instantaneously. The particle temperature instantly increased to 273.15 K along with a corresponding change in phase to a melt ratio of 0.51. The mixed phase particle continued to freeze at constant temperature. Once the fully crystalized state was reached, the ice particle temperature decreased to the wet-bulb temperature of about 235 K, similar in final temperature as Sim Norm.

Aside from the brief sensible cooling from 278.15K to 273.15K, the rate of vaporization was greatest during the freezing stages for both simulations. This is due to the large difference in vapor mass fraction between the surface and ambient air when the particle temperature was 273.15 K. The vaporization rates differed between the two simulations after t = 0.001 s when Sim Norm began to freeze. While the particle temperature decreased as a supercooled liquid in Sim SC, the vapor concentration difference decreased, reducing the mass transport rate from the droplet surface. This can be seen in Figure 9b. Meanwhile, for Sim Norm, the particle maintained a high vapor concentration difference during freezing, favoring a higher mass transport rate. The overall mass that vaporized was greater for the normal freezing simulation for this reason (4.1% vs 3.4%).

The air profiles are depicted in Figure 10. The air temperature increased for both simulations because of the thermal interaction with the significantly warmer particle. [T.sub.air] was greater for Sim Norm during the early part of the simulation because during freezing, the particle maintained a temperature of 273.15K and had greater convective heat transfer to the air than Sim SC. However, the final air temperature was marginally greater for Sim SC because a smaller amount was vaporized which pulled slightly less energy from the air. With a lower final air temperature and greater amount of vaporized water in the air, the final relative humidity was greater for the normal freezing simulation.

Two computational tests were run to identify the differences when particle size varied. Simulations were run for a 10 urn and a 20 urn water droplet. With equal liquid water content, the number of particles was greater for the 10 urn simulation, 1.91 x [10.sup.9] compared to 2.39 x [10.sup.8]. The conditions are listed in Table 4. Results are shown in Figures 11 and 12.

The particles were completely vaporized in a low relative humidity and low pressure environment. The 10 [micro]m particle took nearly a quarter of the time to vaporize compared to the 20 [micro]m particle. The shorter time was due to the greater surface area to volume ratio which enhanced the heat and mass transport rates.

The graphs in Figure 11 and Figure 12 are plotted on a logarithmic scale to highlight the differences that occurred on a small time scale. The final air temperature and relative humidity were identical between the two simulations. This final result was expected because the air conditions and the amount of water content were the same for both simulations. The general trends in all of the graphs are similar except that the 20 urn particle took longer to reach the same final state.

In all of the previous simulations, air temperature, vapor content and air velocity were constant along the length of the icing tunnel, in the absence of an icing cloud. Two simulations were run to demonstrate the differences when varied conditions are introduced. The first simulation, Constant, simulated an icing tunnel where air velocity, air temperature and relative humidity were constant along the length of the tunnel when an ice cloud was absent, while the second simulation, Source, varied these conditions from the injection point to the test section. The variation of these parameters along the tunnel length, while the cloud was absent, followed the CFD calculations that were described in the formulation section. The conditions for both simulations are tabulated in Table 5. The final air temperature, air velocity and specific humidity (also relative humidity) for Sim Source matched Sim Constant in the absence of a cloud. Specific humidity, SH, is the ratio of vapor mass to dry air mass and is provided to show an absolute value of moisture content in the air. The particle velocity changed during Sim Source and was slower overall than Sim Constant, therefore the particle took a longer time to travel the same 2 meters.

The graphs in Figure 13 and Figure 14 are plotted against distance because the heat source, mass source and air velocity centerline profiles were extracted from the previously mentioned CFD calculations. The water droplet in Sim Constant remained a liquid throughout the duration of the simulation, and the air temperature decreased, [DELTA][T.sub.air] = -0.08 K. The particle in Sim Source partially froze due the initially low relative humidity and air temperature. The air temperature decreased, [DELTA][T.sub.air] = -0.06 K, with respect to the target values. In addition, a greater amount of mass vaporized in Sim Source (4.9% vs 3.7%) due to the initially large vapor concentration difference, driving mass transport from the particle.

It is often convenient to characterize a jet spray of water or ice by the median volume diameter, MVD. However, a more accurate representation of a jet spray incorporates the full distribution of particle sizes. A final parametric study was run to illustrate the differences in results between a single-sized particle and a distribution of particle sizes. Conditions are listed in Table 6.

The graphs in Figure 15 and Figure 16 show the results of a uniform particle size cloud simulation, Single, and a full particle size distribution simulation, Multiple. The results of three particle sizes are shown for Sim Multiple. They represent the particle diameter where the cumulative volume of the distribution totaled 10%, 50% and 90%, labeled as DV10, D V50, and DV90 in Figure 15 respectively. The diameters were DV10= 11 [micro]m, DV50 = 40 [micro]m, and DV90 =150 [micro]m where DV50 is by definition the median volume diameter. For direct comparison, the particle diameter of Sim Single matched the particle MVD of Sim Multiple, with a value of 40 [micro]m. All other conditions were identical for the two simulations and are tabulated in Table 6.

The conditions were such that water particles vaporized in a cold and low relative humidity environment. In Sim Single, the particle temperature decreased to freezing and maintained 273.15 K as 32% of the particle mass froze. The convective heat transferred from the warm particle to the cold air was greater than the latent energy due to mass transfer for the duration of Sim Single which increased the air temperature slightly, [DELTA][T.sub.air] = 0.03 K. There was a small increase in relative humidity, [DELTA]RH = 0.9%.

The profiles for the 3 particles in Sim Multiple are different from each other. The DV50 particle profiles in Figure 15 matched Single particle profiles very closely. Since the particles are of equal size and the environmental conditions were similar throughout the simulation, the matching results were expected. The smaller DV10 particle experienced complete freezing and vaporized a greater fraction of its mass, 19% compared to 4% for the DV50 particle. The heat and mass transfer rates were greater for the small particle due to the higher surface area to volume ratio. The larger DV90 particle had a slower response to the environment. Due to the overall small surface area to volume ratio, the particle remained liquid throughout the duration of the simulation and less than 1% of the particle mass vaporized. The cumulative effect of all the differently-sized particles resulted in an overall temperature increase of [DELTA][T.sub.air] = 0.01 K. The initial rise in air temperature occurred because there was convective heat transfer from the warm particle to the cold air. As smaller particles reached the wet-bulb temperature - cooler than the air temperature - overall convective heat was transferred away from the air. The relative humidity gain was greater for Sim Multiple where [DELTA]RH = 2.1%. A greater amount froze in Sim Multiple (Sim Multiple [n.sub.p] = 0.56 vs Sim Single [n.sub.p] = 0.68) where the small particles contributed to this additional freezing. Despite all particles of every size losing mass, the small particles lost more mass compared to the large particles in the distribution and shifted the volume distribution such that MVD increased a small amount from 40 urn to 43.4 urn. Water particles freezing and expanding partially contributed to this increase as well. The MVD in Sim Single decreased marginally from 40 urn to 39.9 [micro]m.

Small particles had a strong influence on the air properties and cloud composition. Despite small differences in T and RH, a more noticeable difference was realized in MVD and melt fraction in the above example. A full particle distribution altered results and more accurately simulated an icing tunnel than a single-sized particle simulation.

MODEL - EXPERIMENT COMPARISON

To test the validity of the thermal model, computational tests were run simulating experiments conducted at the RATFac at NRC during 2012 ([5, 20]). Three simulations were chosen for comparison, (1) a liquid water spray, (2) an ice particle spray (3), and a combined water and ice spray test. Air temperature, humidity, IWC, and LWC measurements at the test section are compared to numerical results. An SAE multiwire was used to measure water content, where the 83mm hotwire measurement was used for LWC (as-measured) and a corrected value of the half-pipe measurement was used for total water content (TWC). IWC was calculated to be the difference between TWC and LWC.

When conducting experiments, it is often easier to measure total conditions (total air temperature, total pressure, etc.). The thermal model in this paper uses static conditions as input. Incompressible flow was assumed when calculating the static values from the experimentally measured total values. Table 7 lists the approximate static conditions that were injected (inj subscript) into the icing tunnel along with the test section target (target subscript) conditions. The experimental values were measured at the centerline of the test section.

For all three simulations, a full distribution of particle sizes was used. The water spray contained de-ionized water, therefore supercooling was simulated. To more closely approximate the flow field, the heat source, vapor mass source and air velocity field profiles were used in the three simulations as well.

Table 8 lists the main changes in air conditions and particle water content (ice and water) that were observed experimentally (exp subscript) and from running the simulation (sim subscript). The change in air temperature, [DELTA][T.sub.air], and change in relative humidity, [DELTA]RH, represent the changes that occurred at the test section when an icing cloud was introduced. The change in gas water content, [DELTA]GWC, is closely related to relative humidity and is listed to show an absolute value for moisture content in the air. The change in liquid water content, [DELTA]LWC, and ice water content, [DELTA]IWC, represent the change in values between the test section and the injection point. The change in overall melt fraction is denoted by [DELTA][eta].

In Scan 877, liquid water was sprayed into a warm environment. Due to the warm air temperature, icing did not occur. The simulation is in agreement with the experiment with respect to the final phase of the cloud, however, the simulation does not fully match the amount of water that evaporated or the amount of air temperature decrease.

In Scan 982, a high concentration of ice was injected into the icing tunnel. Experimentally, 10% ([DELTA][[eta].sub.exp] = 0.10) of the ice cloud melted. In the simulation, the ice reached steady state wet-bulb temperature of approximately 270 K, so no melting occurred. The model simulated only a fraction of mass that vaporized and a portion of the air temperature decrease that had been observed experimentally.

In Scan 1003, a high concentration of ice was injected into the icing tunnel, along with a supplemental amount of liquid water. Melting of the ice particles occurred experimentally. In the model, melting did not occur. The ice (and water) particles reached wet-bulb temperature at approximately 272 K, which is just below melting temperature of the ice. Again, the model underestimated the amount of mass that vaporized and the air temperature decrease.

The model generally predicted the correct trend, but did not fully explain the changes that were observed experimentally at the RATFac icing wind tunnel. In general, the model was able to account for approximately 20% of the air temperature change that was observed experimentally. In addition, the model simulated approximately 20% of the water that vaporized from the icing cloud. It may not be a coincidence that the ~20% values match. When the three simulations were run for a long enough period so that the amount of water vaporized matched the experimental values, the resulting air temperature change was comparable to the experimentally measured values. This suggests that the model may be underestimating the mass transfer rate. A few reasons for the model-experiment discrepancy are offered.

One explanation may be that other sources of water are not accounted for in this model. For example, water/ice film on the tunnel walls have been observed during testing. The water coating was greater especially at lower air velocity tests (such as the tests that were simulated) where the icing cloud has had the chance to disperse more towards the tunnel walls. The long residence time of the water mass on the tunnel walls may have contributed to the observed changes.

The vapor mass source term may be an area of inaccuracy as well. As mentioned in the formulation section, internal calculations were used to approximate the temperature profile along the icing tunnel centerline in the absence of an icing cloud. Similarity was used to approximate the vapor concentration profile. If the profile overestimated the vapor concentration, then the difference between the air and surface vapor concentration was reduced and the overall mass that vaporized was underestimated.

Also, there may have been increased energy and mass transfer that occurred at the injection nozzle that was not accounted for in the model. The expansion of air at the injection nozzle was likely turbulent, which would have enhance heat and mass transfer between the flowing air and the ejected particles.

Finally, there are uncertainties in experimental measurements. However, the values for LWC/IWC, air temperature, and humidity were measured by independent instruments. Since the air temperature change thermally corresponds to the amount of mass that vaporized, it is not likely that experimental error accounts for much of the discrepancy.

CONCLUSION

Changes in air temperature and air moisture content were experimentally measured when an icing cloud was activated at the RATFac icing tunnel at NRC. A thermal model was written to simulate the icing cloud and air interaction in an attempt to explain these observed changes. The model couples the conservation of energy and mass equations between the air and icing cloud particles. A parametric study was conducted to identify parameter impact on results. Due to the thermal mass dominance of the air mass, air temperature, pressure (mass), and vapor content were the largest factors in the air-particle interaction. These three parameters also determine the wet-bulb temperature, the quasi-steady-state temperature vaporizing particles reached.

Three computational tests were run to simulate experiments conducted at NRC. The model was able to simulate the general trends, but it did not fully explain the experimentally measured changes. Several explanations are offered to explain the discrepancy between experimental values and the model simulation. Further work is necessary to ascertain the incongruity. It is likely that all evaporation--such as evaporation from water film on the walls--is not accounted for properly in the model. Future work will investigate experimental uncertainties and their potential impact to the results. It is desired to run simple bench tests, such as a water droplet evaporating in still air, to verify individual components of the model.

The simulation of Scan 1003 is a persuasive example of how a small 1 K difference can mark the difference between an event and a non-event. Ice melting occurred experimentally but the model did not simulate melting because the ice particle temperature remained 1 K below melting temperature. It is widely accepted that partial melting of ice-crystals is necessary for engine icing to occur. An additional 1 K increase in the simulation would have led to melting in the simulation, which, as mentioned, has great implications for engine icing. It is for this reason that a reasonable level of fidelity is needed not only experimentally, but also numerically to understand engine icing and to avoid the problems that can ensue. While the model presented in this paper was specifically written to simulate specific experiments at the RATFac at NRC, it can be modified to simulate the air-particle interaction of other engine icing test facilities.

REFERENCES

[1.] Mason, J., Strap, J., and Chow, P., "The Ice particle Threat to Engines in Flight," presented at 44th AIAA Aerospace Sciences Meeting and Exhibit 2006, USA, January 9-12, doi:10.2514/6.2006-206.

[2.] Mason, I., Chow, P., and Fuleki, D., "Understanding Ice Crystal Accretion and Shedding Phenomenon in Jet Engines Using a Rig Test," presented as ASME Turbo Expo 2010, UK, June 14-18, ASME GT2010-22550, doi:10.1115/1.4002020.

[3.] Struk, P., Currie, T., Wright, W., Knezevici, D. et al., "Fundamental Ice Crystal Accretion Physics Studies", SAE International Technical Paper 2011-38-0018. 2011, doi:10.4271/2011-38-0018.

[4.] Currie, T., Struk, P., Tsao, J., Fuleki, D. et al., "Fundamental Study of Mixed-Phase Icing with Application to Ice Crystal Accretion in Aircraft Jet Engines," presented at 4th AIAA Atmospheric and Space Environments Conference 2012, USA June 25-28, doi:10.2514/6.2012-3035.

[5.] Struk, P., Bencic, T., Tsao, J., Fuleki, D. et al., "Preparation for Scaling Studies of Ice-Crystal Icing at the NRC Research Altitude Test Facility," presented at 5th AIAA Atmospheric and Space Environments Conference 2013, USA, June 24-27, 2013, doi:10.2514/6.2013-2675.

[6.] Fuleki, D., Mahallati, A., Knezevici, D., Currie, T., et al., "Development of a Sensor for Total Temperature and Humidity Measurements under Mixed-Phase and Glaciated Icing Conditions," presented at 6th AIAA Atmospheric and Space Environments Conference 2014, USA, June 16-20, 2014, doi:10.2514/6.2014-2751.

[7.] Willbanks, C., and Schulz, R., "Analytical Study of Icing Simulation for Turbine Engines in Altitude Test Cells," AEDC-TR-73-144, 1973.

[8.] Hindmarsh, J., Russell, A., and Chen, X., "Experimental and numerical analysis of the temperature transition of a suspended freezing water droplet," Internationaljournal of Heat and Mass Transfer 46: 1199-1213, 2003.

[9.] Strub, M., Jabbour, O., Strub, F., and Bedearrats, J., "Experimental study and modelling of the crystallization of a water droplet," International Journal of Refrigeration 26: 59-68, 2003.

[10.] Wright, W., Jorgenson, P., and Veres, J., Mixed Phase Modeling in GlennICE with Application to Engine Icing, presented at AIAA Atmospheric and Space Environments Conference 2010, Canada. August 2-5, 2010, doi:l0.2514/6.2010-7674.

[11.] Hauk, T., Roisman, I., and Tropea, C., "Investigation of the Melting Behavior of Ice Particles in an Acoustic Levitator," presented at AIAA Aviation 2014, USA, June 16-20, 2014, doi:10.2514/6.2014-2261.

[12.] Villedieu, P., Trontin, P., and Remi, C., "Glaciated and mixed-phase ice accretion modeling using ONERA 2D icing suite," presented at AIAA Aviation 2014, USA, June 16-20, 2014, doi:10.2514/6.2014-2199.

[13.] Dalton, J., "Experimental essays on the constitution of mixed gases; on the force of steam or vapor from water and other liquids in different temperatures, both in a Torricellian vacuum and in air; on evaporation and on the expansion of gases by heat," Memoirs and Proceedings of the Manchester Literary and Philosophical Society, 1802, 5, 535-602.

[14.] Incropera, F., and DeWitt, D., "Fundamentals of Heat and Mass Transfer, Fourth Edition," (New York, John Wiley & Sons, 1996) ISBN: 0471304603.

[15.] Morrison, F., "An Introduction to Fluid Mechanics," (New York, Cambridge University Press, 2013) pg. 625, ISBN: 1107003539.

[16.] Knezevici, D. C., Fuleki, D., and MacLeod, J., "Development and Commissioning of the Linear Compressor Cascade Rig for Ice Crystal Research," SAE Technical Paper 2011-38-0079. doi:10.4271/2011-38-0079.

[17.] Heverly, J., "Supercooling and Crystallization," Transactions of American Geophysical Union 30(2): 205-10, 1949.

[18.] Knezevici, D., Fuleki, D., Currie, T., Galeote, B., et al., "Particle Size Effects on Ice Crystal Accretion - Part II," presented at 5th AIAA Atmospheric and Space Environments Conference 2013, USA, June 24-27, 2013, doi:10.2514/6.2013-2676.

[19.] Van Zante, J., Ide, R., Steen, L., and Acosta, W., "NASA Glenn Icing Research Tunnel: 2014 Cloud Calibration Procedure and Results," NASA/TM-2014-218392.

[20.] Struk, P., Bartkus, T., Tsao, J., Currie, T., et al., "Ice Accretion Measurements on an Airfoil and Wedge in Mixed-Phase Conditions," Presentation at SAE 2015 International Conference on Icing of Aircraft, Engines, and Structures, June 2015.

[21.] The Engineering Toolbox, "Water Vapor - Specific Heat," http://www.engineeringtoolbox.com/water-vapor-d_979.html. accessed Jan, 2015.

[22.] C: Pruppacher, H., and Klett, J., "Microphysics of Clouds and Precipitation, Second Edition," (Dordrecht, Springer, 2010), doi :10.1007/978-0-306-48100-0.

[23.] Dorsey, N., "Properties of Ordinary Water-Substance in all its Phases: Water-vapor, Water, and all the Ices," (New York, Reinhold Publishing Corporation, 1940) ISBN: 1114273740.

[24.] Goff, J., and Gratch, S., "Low-pressure properties of water from -160 to 212[degrees]F.," Transactions of the American Society of Heating and Ventilating Engineers 52, 95-121, 1946.

[25.] The Engineering Toolbox, "Ice - Thermal Properties," http://www.engineeringtoolbox.com/ice-thermal-properties-d_576.html. accessed Jan, 2015.

[26.] West, R., and Lide, D., "CRC Handbook of Chemistry and Physics, 70th Edition," (Boca Raton, CRC Press, 1989), ISBN: 0849304709.

[27.] The Engineering Toolbox, "Density of Fluids - Changing Pressure and Temperature," http://www.engineeringtoolbox.com/fluid-density-temperature-pressure-d_309.html. accessed Jan, 2015.

CONTACT INFORMATION

Tadas Bartkus, Ph.D.

Work phone: (216)433-6915

tadas.p.bartkus@nasa.gov

Affiliation: Ohio Aerospace Institute

ACKNOWLEDGEMENTS

The authors wish to acknowledge the financial support from the Atmospheric Environment Safety Technology Project (AEST) of the NASA Aviation Safety Program. The first author is supported currently under a NASA Glenn ARTS contract.

NOMENCLATURE

a - acceleration (m [s.sup-2])

A - area ([m.sup.2])

Bi - Biot number

C - specific heat capacity (J [kg.sup.-1] [K.sup.-1])

[C.sub.d] - coefficient of drag

d - particle diameter (m)

[D.sub.ab] - diffusivity of water vapor in air ([m.sup.2] [s.sup.-1])

F - force (kg m [s.sup.-2])

GWC - gas water content or vapor (g [m.sup.3])

h - heat transfer coefficient (W [m.sup.-2] [K.sup.-1])

H-enthalpy (J)

[h.sub.m] - mass transfer coefficient (m [s.sup.-1])

IWC - ice water content (g [m.sup.3])

k - thermal conductivity (W [m.sup.-1] [K.sup.-1])

L - latent heat of phase change (J [kg.sup.-1])

[L.sub.c] - characteristic length (m)

LWC - liquid water content (g [m.sup.3])

m - mass (kg)

MVD - median volume diameter (m or [micro]m)

MW - molecular weight (g [mol.sup.-1])

Nu - Nusselt number

P - pressure (Pa)

Pr - Prandtl number

q - heat flux (W)

Re - Reynolds number

RH - relative humidity (%)

Sc - Schmidt number

Sh - Sherwood number

SH - Specific humidity (vapor mass/dry air mass)

t - time (s)

T - temperature (K)

U - slip velocity (m [s.sup.-1])

v- velocity (m [s.sup.-1])

x - distance (m)

[partial derivative] - differential

[DELTA] - change

[eta] - melt fraction (liquid water mass/total mass)

[mu] - dynamic viscosity (Pa s)

[rho] - density (kg nr3)

[phi] - evaporation/sublimation limit function

[omega] - water vapor mass fraction

# - total number of particles in a bin

SUBSCRIPTS

air - ambient air

cond - condensation

dep - deposition

drag - drag

dry - dry (absence of water content)

evap - evaporation

exp - experiment

freeze - freezing

g - gravity

hc - homogeneous crystallization

i - [i.sup.th] number

ice - ice

inj - injection point

latent - latent mass transfer energy

limit - RH at which mass transfer is limited

melt - melting

n - number of bins

p - particle (ice, water or mixed phase)

sat - saturation

sens - sensible energy

sim - simulation

source - source (energy or mass)

subl - sublimation

super - supercooling

target - target value desired at test section

test - test section value

water - water

wb - wet-bulb

wv - water vapor

APPENDIX

Several air, water, ice and vapor thermo-physical properties are dependent on temperature, pressure and relative humidity. Functions were generated to approximate the changing properties. Polynomial equations were fit to published thermo-physical data. The equations used for these changing properties are detailed in Table 9.

Tadas P. Bartkus Ohio Aerospace Institute

Peter Struk NASA John Glenn Research Center

Jen-Ching Tsao Ohio Aerospace Institute

Table 9. Thermo-physical properties for ice, water, water vapor and air used in the thermal model. Parameter (Units) Input (Units) [C.sub.air(dry)] - (J [kg.sup.-1] [k.sup.-1]) [14] [C.sub.wv](J T(K) [kg.sup.-1] [k.sup.-1]) ([21]) [C.sub.ice](J T([degrees]c) [kg.sup.-1] [k.sup.-1]) ([22]) [C.sub.water](J T([degrees]C) [kg.sup.-1] [k.sup.-1]) ([22]) (T > 0[degrees]C) [C.sub.water] T([degrees]c) (J [kg.sup.-1] [k.sup.-1])([22]) (T<0[degrees]C) [D.sub.air,wv] T(K), P(Pa) ([m.sup.2] [s.sup.-1]) [23] [K.sub.air] T(K) (W [m.sup.-1] [K.sup.-1]) [14] [L.sub.evap] or T([degrees]c) [L.sub.cond] (J [Kg.sup.-1]) ([22]) [L.sub.melt] - or [L.sub.freeze] (J [Kg.sup.-1]) ([22]) [L.sub.subl] T([degrees]C) or [L.sub.dep] (J [Kg.sup.-1]) ([22]) [P.sub.sat] (Pa) ([24]) T(K) [[mu].sub.air] T(K) (Pa s) [14] [[rho].sub.air T(K), P(Pa), RH(%) (humid)] (Kg [m.sup.-3]) [[rho].sub.ice] T([degrees]C) (Kg [m.sup.-3]) ([25]) [[rho].sub.water] T([degrees]C), P(Pa) (Kg [m.sup.-3]) ([26,27]) Parameter (Units) Value or Expression [C.sub.air(dry)] 1005 (J [kg.sup.-1] [k.sup.-1]) [14] [C.sub.wv](J 0.7459[T.sup.3] x [10.sup.-7] + 0.4839 [kg.sup.-1] [T.sup.2] x [10.sup.4] - 0.2458T+ 1874.4317 [k.sup.-1]) ([21]) [C.sub.ice](J (0.503 +0.00175T) x 4186.8 [kg.sup.-1] [k.sup.-1]) ([22]) [C.sub.water](J [0.9979 + 0.0000031[(T - 35).sup.2] + [kg.sup.-1] 0.0000000038*[(T - 35).sup.4]] x 4186.8 [k.sup.-1]) ([22]) (T > 0[degrees]C) [C.sub.water] [1.0074 + 0.0000829[T.sup.2]] x 4186.8 (J [kg.sup.-1] [k.sup.-1])([22]) (T<0[degrees]C) [D.sub.air,wv] 0.0000219 x [(77273).sup.1.853] x (101325/P) ([m.sup.2] [s.sup.-1]) [23] [K.sub.air] -3.2432[T.sup.2] x [10.sup.-8]+ 9.7809T x (W [m.sup.-1] [10.sup.-5]- 1.4252 x [10.sup.-4] [K.sup.-1]) [14] [L.sub.evap] or (2500.8 - 2.367T+ 0.0016[T.sup.2] - 0.00006 [L.sub.cond] [T.sup.3]) x 1000 (J [Kg.sup.-1]) ([22]) [L.sub.melt] 333,700 or [L.sub.freeze] (J [Kg.sup.-1]) ([22]) [L.sub.subl] (2834.1 -0.29T - 0.004[T.sup.2]) x 1000 or [L.sub.dep] 10^[-7.90298 x (373.16/T - 1) + 5.02808 x (J [Kg.sup.-1]) ([22]) log10(373.16/T) [P.sub.sat] (Pa) ([24]) - 0.00000013816 x(10^(11.344 x (1 -T/373.16))- 1) + 0.0081328 x (10^(-3.49149 x (373.16/T - 1))- 1) + [log.sub.10](1013.246)] [[mu].sub.air] -2.786813[T.sup.2] x [10.sup.-11] + 6.619341T x [10.sup.-8]+ 1.001538 x [10.sup.-6] (Pa s) [14] [[rho].sub.air [(P- Psat(T)*RH/100) x 0.028964 + (humid)] (Psat(T)*RH/100) x 0.018016] / ( 8.3144621T) (Kg [m.sup.-3]) -2.09927"6x 10"'[degrees]-7.4U2r5x 10"8- 1.025397"4x 10s - 6.97301 7"3 x 10"" [[rho].sub.ice] - 0.02371[T.sup.2] - 0.4360T + 916.1152 (Kg [m.sup.-3]) ([25]) (-5.5493[T.sup.6] x [10.sup.11] + 1.5788 [T.sup.5] x [10.sup.-8] - 1.7941[T.sup.4] x [10.sup.-6]+ 1.2051[T.sup.3] x [10.sup.-4] [[rho].sub.water] -9.1638[T.sup.2] x [10.sup.-3] + 0.06367T + (Kg [m.sup.-3]) ([26,27]) 999.8869)/(1 -(P- 101325)/2.15 x [10.sup.-9]) Table 1. Initial conditions for 2 baseline simulations for parametric analysis. Units Baseline 1 Baseline 2 Particle Diameter [micro]m 10 10 Water Content g/[m.sup.3] LWC=1 IWC=1 Pressure Pa 88,000 88,000 Relative Humidity % 50 50 Air Temperature K 278.15 280.15 Particle Temperature K 278.15 271.15 Slip Velocity m/s 5 5 Table 2. Table depicting the parameter changed for Sim 1 [right arrow] 6 with respect to the corresponding baseline simulation. Sim# Parameter Changed Units Value Baseline 1 1 Slip Velocity m/s 25 2 Particle Temperature K 273.15 3 Pressure Pa 28,000 Baseline 2 4 Water Content (IWC) g/[m.sup.3] 5 5 Relative Humidity % 80 6 Air Temperature K 271.15 Table 3. Initial conditions for supercooled and normal freezing simulations. Units Conditions Particle Diameter [micro]m 20 Water Content (LWC) g/[m.sup.3] 1 Pressure Pa 88,000 Relative Humidity % 10 Air Temperature K 235.15 Particle Temperature K 278.15 Slip Velocity m/s 5 Table 4. Initial conditions for two differently-sized particles composing a uniform particle size icing cloud. Units 10 [micro]m Particle Diameter [micro]m 10 Water Content (LWC) g/[m.sup.3] 1 Pressure Pa 28,000 Relative Humidity % 20 Air Temperature K 278.15 Particle Temperature K 278.15 Slip Velocity m/s 5 Number of Particles # 1.91 x [10.sup.9] Evaporation Time s 0.1634 20 [micro]m Particle Diameter 20 Water Content (LWC) 1 Pressure 28,000 Relative Humidity 20 Air Temperature 278.15 Particle Temperature 278.15 Slip Velocity 5 Number of Particles 2.39 x [10.sup.8] Evaporation Time 0.619 Table 5. Conditions for Sim Constant and Sim Source. Units Constant Particle Diameter [micro]m 20 Water Content (LWC) g/[m.sup.3] 1 Pressure Pa 88,000 Initial SH [g.sub.vapor]/[Kg.sub.dry air] 3.58 (Initial RH) % (50) Target SH [g.sub.vapor]/[Kg.sub.dry air] 3.58 (Target RH) % (50) Initial Air Temperature K 280.15 Target Air Temperature K 280.15 Particle Temperature K 280.15 Air Velocity m/s 85 Slip Velocity m/s 5 Time s 0.0235 Distance m 2.0 Source Particle Diameter 20 Water Content (LWC) 1 Pressure 88,000 Initial SH 0.80 (Initial RH) (50) Target SH 3.58 (Target RH) (50) Initial Air Temperature 260.15 Target Air Temperature 280.15 Particle Temperature 280.15 Air Velocity Varied Slip Velocity Varied Time 0.0293 Distance 2.0 Table 6. Initial conditions for a single-sized particle cloud (Single) and a cloud composed of a full distribution of particle sizes (Multiple). Units Value Water Content (LWC) g/[m.sup.3] 1 Pressure Pa 88,000 Relative Humidity % 20 Air Temperature K 268.15 Particle Temperature K 278.15 Slip Velocity m/s 5 MVD (Single) [micro]m 40 DV10 (Multiple) [micro]m 11 DV50 (Multiple) [micro]m 40 DV90 (Multiple) [micro]m 150 Table 7. Conditions for three experimental tests that were simulated by the thermal model. Injected and target values at the test section are given. Listed air conditions are for when the icing cloud is absent. Units Scan 877 [T.sub.air,inj] K 262.0 [T.sub.air,target] K 288.4 [SH.sub.inj] [g.sub.vapor]/[Kg.sub.dry air] 0.07 ([RH.sub.inj]) % (2) [SH.sub.target] [g.sub.vapor]/[Kg.sub.dry air] 4.07 ([RH.sub.target]) % (16.1) [GWC.sub.inj/target] g/[m.sup.3] 2.1 [LWK.sub.inj/target] g/[m.sup.3] 1.0 [IWC.sub.inj/target] g/[m.sup.3] 0 [[eta].sub.inj/target] - 1.0 P Pa 42806 [v.sub.air,inj/target] m/s 86.8 [T.sub.water,inj] K 278.15 [T.sub.ice,inj] K - [MVD.sub.water,inj] [micro]m 40.0 [MVD.sub.ice,inj] [micro]m - Scan 982 Scan 1003 [T.sub.air,inj] 256.1 256.9 [T.sub.air,target] 278.0 277.9 [SH.sub.inj] 0.07 0.07 ([RH.sub.inj]) (5) (4) [SH.sub.target] 2.88 3.81 ([RH.sub.target]) (35.4) (46.5) [GWC.sub.inj/target] 2.4 3.0 [LWK.sub.inj/target] 0 1.9 [IWC.sub.inj/target] 8.4 8.6 [[eta].sub.inj/target] 0.0 0.18 P 66478 65934 [v.sub.air,inj/target] 85.7 84.1 [T.sub.water,inj] - 278.15 [T.sub.ice,inj] 256.15 256.15 [MVD.sub.water,inj] - 40.0 [MVD.sub.ice,inj] 45.5 45.5 Table 8. Comparison of the simulation results with experimental results. Values are given in terms of change at the test section. Units Scan 877 [DELTA][T.sub.air, exp] K -2.6 [DELTA][T.sub.air,sim] K -0.54 [DELTA][SH.sub.exp] [g.sub.vapor]/[Kg.sub.dry air] 1.19 [DELTA][SH.sub.sim] [g.sub.vapor]/[Kg.sub.dry air] 0.23 [DELTA][GWC.sub.exp] g/[m.sup.3] 0.5 [DELTA][GWC.sub.sim] g/[m.sup.3] 0.13 [DELTA][LWC.sub.exp] g/[m.sup.3] -0.39 [DELTA][LWC.sub.sim] g/[m.sup.3] -0.13 [DELTA][IWC.sub.exp] g/[m.sup.3] 0.0 [DELTA][IWC.sub.sim] g/[m.sup.3] 0.0 [DELTA][[eta].sub.exp] - 0.00 [DELTA][[eta].sub.sim] - 0.00 Scan 982 Scan 1003 [DELTA][T.sub.air, exp] -2.9 -4.4 [DELTA][T.sub.air,sim] -0.75 -0.88 [DELTA][SH.sub.exp] 1.31 1.37 [DELTA][SH.sub.sim] 0.15 0.19 [DELTA][GWC.sub.exp] 1.03 1.23 [DELTA][GWC.sub.sim] 0.16 0.21 [DELTA][LWC.sub.exp] 0.76 0.11 [DELTA][LWC.sub.sim] 0.0 -0.09 [DELTA][IWC.sub.exp] -1.79 -1.34 [DELTA][IWC.sub.sim] -0.16 -.12 [DELTA][[eta].sub.exp] 0.10 0.035 [DELTA][[eta].sub.sim] 0.00 -0.004