Comparisons of Mixed-Phase Icing Cloud Simulations with Experiments Conducted at the NASA Propulsion Systems Laboratory

Comparisons of Mixed-Phase Icing Cloud Simulations with

Experiments Conducted at the NASA Propulsion Systems

Laboratory

Tadas P. Bartkus Ohio Aerospace Institute, Cleveland, Ohio 44135, USA Peter M. Struk NASA Glenn Research Center, Cleveland, OH 44135, USA and Jen-Ching Tsao Ohio Aerospace Institute, Cleveland, Ohio 44135, USA This paper presents the evaluation of a numerical model for simulation of the icing cloud development at N ASA Glenn Research Center’s Propulsion Systems Laboratory (PSL). The model is helping the icing facility and the fundamental ice-crystal icing physics research team to better understand the complex interactions between the test parameters and have greater confidence in the conditions at the test section of the PSL tunnel. The model attempts to explain the observed changes in test conditions by coupling the conservation of mass and energy equations for both the cloud particles and flowing air, while accounting for compressibility and the variable PSL geometry. A subroutine has been added to more accurately simulate the tunnel when water vapor conditions potentially exceed saturation. The model simulation results are compared to experimentally measured values that were taken during the first fundamentals of ice-crystal icing physics tests conducted at PSL in March 2016. The tests simulated ice-crystal and mixed-phase icing that relate to ice accretions within turbofan engines. Experimentally measured air temperature, humidity, total water content, liquid and ice water content, as well as cloud particle size, are compared with model predictions. The model showed good trend agreement with experimentally measured values, but often over-predicted aero-thermodynamic changes. This discrepancy is likely attributed to radial variations that this one-dimensional model does not address. One of the key findings of this work is that greater aero-thermodynamic changes occur when humidity conditions are low. In addition a range of mixed-phase clouds can be achieved by varying only the tunnel humidity conditions, but the range of humidities to generate a mixed-phase cloud becomes smaller when clouds are composed of smaller particles. In general, the model predicted melt fraction well, in particular with clouds composed of larger particle sizes.

Nomenclature 1-D = one-dimensional 2-D = two-dimensional CDP = Cloud Droplet Probe (by Droplet Measurement Technology, Inc.)

CIP = Cloud Imaging Probe (by Droplet Measurement Technology, Inc.)

Cond# = condition number exp = experiment HSI = High Speed Imaging probe (Artium Technologies, Inc.)

Senior Research Associate, Icing Branch, 21000 Brookpark Road, MS 110-3, AIAA Member.

Aerospace Engineer, Icing Branch, 21000 Brookpark Road, MS 11-2, AIAA Senior Member.

Principal Research Scientist, Icing Branch, 21000 Brookpark Road, MS 11-2, AIAA Associate Fellow.

American Institute of Aeronautics and Astronautics IKP2 = Isokinetic Probe, version 2 (Science Engineering Associates, Inc.)

-3 IWC = ice water content, defined as ice water mass per unit of volume (g m ) -3 LWC = liquid water content, defined as liquid water mass per unit of volume (g m ) MVD = median volumetric diameter (μm) P = pressure (kPa or psi) PDI = Phase-Doppler Interferometer probe (Artium Technologies, Inc.)

PSL = NASA Propulsion Systems Laboratory RFTP = rearward-facing temperature probe RH = relative humidity (%) sim = simulation T = temperature or air temperature (K, °C, or °F) Twb = wet-bulb temperature (K, °C, or °F) -3 TWC = total water content, defined as the sum of liquid water mass and ice water mass (g m ) -1 U = air velocity (m s ) Δ = change 𝜂 = melt fraction, defined as liquid water mass divided by the total water content (non-dimensional) -1 𝜔 = mass mixing ratio, defined as water vapor mass divided by the dry air mass (g kg ) Subscripts 0 = total 100%RH = value at 100% relative humidity aero = aerosol e = exit exp = experiment i = inlet (or initial) none = none (does not use aerosol subroutine) off = spray off s = static sim = simulation water = water I. Introduction INCE the 1990’s, several jet engine power-loss events of commercial aircraft have been reported. Mason et al.

discussed how at high altitudes, ice crystals entering the jet engine core can partially melt within the warm engine S and refreeze on internal components, resulting in uncommanded power-loss. To better understand this phenomenon and determine the physical mechanism of ice-crystal icing within jet engines, a series of tests has been collaboratively conducted between NASA Glenn Research Center and the National Research Council of Canada at the Research 2,3 Altitude Test Facility between 2010 and 2012. Most recently, an initial study on the fundamentals of ice-crystal icing physics was conducted at the NASA Propulsion Systems Laboratory (PSL) in March 2016. The long-term goal of these tests is to improve aero-thermodynamic modeling tools that assess icing in current and future N+2/N+3 aircraft.

5-7 During the 2016 icing tests, and on previous occasions, it has been observed that the test conditions, most notably air temperature and humidity, change when the icing cloud is activated. In some cases, measurements of air temperature decreased by as much as several degrees Celsius, accompanied by a measurable increase in water vapor content. It is hypothesized that a thermal exchange occurs between the water particles of the cloud and the air flow, changing the air temperature and vapor content as the cloud and air masses travel down the length of the tunnel to the test section.

8-13 Existing thermal models show little to no effect due to one-way coupling of the ice and water particles with the flowing air, approximating the air mass as an infinite thermal reservoir with unchanging properties. A thermal model, 14,15 previously written by Bartkus et al. attempts to explain the observed changes in test conditions by fully coupling the conservation of mass and energy equations between the icing cloud and the flowing air mass. This paper builds on the previous work by comparing simulations with experiments run during the first fundamentals of ice-crystal icing physics tests. Specifically, experimentally measured air temperature, humidity, total water content, liquid and ice water content, as well as cloud particle size are compared with model predictions.

American Institute of Aeronautics and Astronautics II. Experiment Description The NASA PSL is an altitude jet-engine test facility that generates ice particles using a liquid water spray nozzle 16-18 system whereby the injected water droplets can freeze out as the cloud flows towards the engine/test section. The droplets freeze due to a combination of convective heat transfer and evaporative cooling. Figure 1 depicts the PSL geometry used for the icing tests. The PSL icing tunnel has a 27:1 area contraction ratio from the plenum (Tunnel Inlet) to the test section (Tunnel Exit), with an axial distance of about 8.84 m (29.0 ft). The spray nozzles and spray bar system are located at the tunnel inlet in the plenum, whereas the test section is a 0.91-m (36-in) diameter free-jet exit. PSL has the capability to spray either de-ionized water or city water (i.e. untreated water). City water was used for all icing tests conducted. A spray nozzle configuration was chosen that attempted to maintain the center 0.15 m (6 in) diameter area approximately uniform at the test section, and contain the entire cloud within an approximately 0.61-m (24-in) diameter area. Figure 2 shows a diagram of the spray bar geometry at the PSL tunnel with the two spray nozzle patterns that were primarily used for these tests. The red and blue ellipses represent the approximate initial area coverage for the two different spray patterns.

Of the total eight days of testing, two days were dedicated to airfoil ice accretion testing and six days to characterize the cloud. Figure 3 shows the various test configurations that were used as viewed from the plenum looking aft towards the tunnel exit. The cloud characterization consisted of spray pattern optimization and uniformity characterization (Fig. 3A), mixed-phase generation and measurements (Fig. 3B), particle size measurements (Fig. 3D and 3E), and total water content measurements (Fig. 3F). The two days of ice accretion testing using a NACA 0012 airfoil is shown in Fig. 3C.

Several test conditions were varied including air temperature, pressure, velocity and humidity levels, as well as the cloud total water content, particle inlet temperature, and particle size distribution. The liquid and ice water content 19-21 of the cloud was measured using the Science Engineering Associates Multiwire probe (Fig. 3B), and the total water content was measured using the Science Engineering Associates isokinetic probe, version 2 (IKP2) as seen in Fig. 3F. Particle size measurements were made using the Droplet Measurement Technologies ’ Cloud Droplet Probe (CDP), which can measure droplet sizes ranging from 2 to 50 μm in diameter, and the Cloud Imaging Probe (CIP), for sizes ranging from 15 to 960 μm in diameter, which are seen in Fig. 3D and 3E, respectively. In addition, particle size data were collected using High Speed Imaging (HSI) and Phase-Doppler Interferometer (PDI) probes that are being developed by Artium Technologies, Inc. The HSI and PDI probes measured particle size at a location approximately 0.15 m (6 in) above or below the centerline, the location of which is represented by the crosshair symbols ( ⊕ ) in Fig. 3B-3D. Humidity was measured (both spray-off and spray-on values) by as many as two different rearward-facing humidity probes, depending on the configuration, and ar e labeled as Humidity Probe “A” and Humidity Probe “B” in Fig. 3 . Similarly, air temperature was measured (both spray-off and spray-on values) by a rearward-facing temperature probe (RFTP) and/or a commercial forward-facing total air temperature probe, depending on the configuration. Finally, cloud uniformity was gauged using a tomography system that generates a two-dimensional (2-D), time-averaged, intensity map of the cloud across a portion of the 0.91-m (36-in) duct. The tomography system was used in every test configuration. A combined temperature and humidity traversing probe (Fig. 3A) was also used to measure aero-thermal uniformity at the tunnel exit. Temperature and humidity data presented in this paper were taken from the combined temperature-humidity traversing probe, the RFTP and Humidity Probe “B”. Measurements made with the commercial total air temperature probe and Humidity Probe “A” are believed to be less accurate and are not presented in this paper. A detailed description of the experimental configuration, the data collected, and results of the first fundamental ice-crystal icing physics tests can be found elsewhere.

III. Model Description The one-dimensional (1-D) thermal model simulates an icing wind tunnel by calculating particle and air properties using expressions derived from the conservation of mass, momentum and energy equations. Differential equations for change in air temperature, particle temperature, air mass, particle mass, air velocity, particle velocity, particle melt fraction, air density, and pressure were derived from equations of conservation. Compressible flow expressions were incorporated into the model to relate air temperature, velocity, density and pressure with respect to the changing cross-sectional area of the PSL. The reader is referred to Bartkus et al. for all listed assumptions and equations related to the aero-thermodynamic changes of the cloud and flowing air.

The model simulates both ice particle melting and water droplet freezing, including spontaneous latent heat release during freezing for supercooled droplets. In the model, humidity changes are due to ice particle sublimation or deposition and water droplet evaporation or condensation. The associated latent heating or cooling and mass transfer change the ice water content ( IWC ) and liquid water content ( LWC ) from the moment of injection to the time the cloud reaches the test section. As a result, the model predicts changes in air temperature, humidity, IWC and LWC , as well American Institute of Aeronautics and Astronautics as the particle phase, size, temperature and velocity at the test section. To more accurately simulate spray injection, the model incorporates a distribution of particle sizes (20+ bins) as an initial input condition.

A. Supersaturation Aerosol Condensation Subroutine Experimental test conditions can exist where the water vapor content can reach saturation or even possibly exceed saturation (supersaturation) at the tunnel exit. For example, elevated humidity levels in the tunnel inlet (plenum) can lead to supersaturation at the tunnel exit, due to a rapid static temperature decrease from compressible air expansion in the tunnel contraction. With higher air velocity flow rates, it is possible that the vapor potentially in excess of saturation does not condense (or deposit as ice) quickly enough to not exceed saturation, in particular if the condensation is strictly limited to the rate as physically allowed through diffusion. In the previous version of the model by Bartkus et al., this supersaturation condition was dealt with in one of two ways.

The first method simply condensed (or deposited) on existing particles through diffusion only, when conditions exceeded saturation. Using this method, condensing on existing particles through diffusion was not enough to reduce the vapor content of the supersaturated air down to saturation by the time the humid parcel of air reached the tunnel exit. In general, the total time from the moment a parcel of air became supersaturated part way in the tunnel, to the tunnel exit, was on the order of about 0.1 seconds.

The second method dealt with the potentially supersaturated situation by limiting the humid air to never exceed saturation. In order to not exceed saturation, the algorithm implemented instantaneous condensation. Instantaneous condensation is referred to the amount of vapor mass that must condense in order to not exceed saturation. The resulting condensation rate is greater than what could occur through diffusion alone. In the model, no new particles were created for this approach; all condensation occurred on existing cloud particles only.

A subroutine has been added to the model that attempts to simulate these potentially supersaturated conditions more accurately. In the absence of aerosol particulates, supersaturation can reach several hundred percent in a 25,26 quiescent environment. If the conditions are pristine enough, supersaturation can reach high enough levels that homogeneous condensation can occur. These high supersaturation levels are necessary due to the high activation energy barrier required for this type of phase transition. The PSL icing tunnel, however, used ambient air that was not filtered and likely contained aerosol particulates. Water vapor in excess of saturation readily condenses on aerosol particulates, which is referred to as heterogeneous condensation. In fact, even in nature, the typical limit of supersaturation is no higher than about 1%, as aerosol particles serve as surfaces on which water vapor can condense.

It reasonable that supersaturated air flowing within the icing tunnel condenses not only on existing cloud particles, but also on aerosol particulates that naturally exist in the atmosphere. Thus, this subroutine condenses supersaturated air on existing cloud particles and existing aerosol particulates through diffusion.

A great deal of literature can be found on the chemical content, concentration number and distribution size of airborne aerosol particulates. The particulates are both organic and inorganic in composition. The number and size of particulates is not constant as the values vary with respect to time and location. Research has shown that there can be seasonal variation in these values as there are higher aerosol particulate concentrations during the heating winter months as explained by the increase in fossil fuel combustion. Diurnal variations of particulate number exist as the maximum concentration occurs during peak traffic hours. Measurements of particulate concentrations have been 31 3 measured at times over 100,000 per cubic centimeter of air if near a city background , or as low as 3,100/cm in more rural areas such as the Alps. A study conducted in Pittsburgh, PA, measured a particulate concentration of 3 33 22,000/cm with a mode particulate size of 0.04 μm (40 nm). A distribution of particulate sizes were measured that ranged from 0.003 μm (3 nm) to 2.5 μm (2500 nm) in this Pittsburgh study. For this work, an aerosol particle concentration of 22,000/cm with a uniform aerosol distribution size of 0.04 μm is used. NASA Glenn Research Center is located in Cleveland, OH, and since Pittsburgh is similar in population and nearby in geography, these aerosol particulate values are approximated for our air quality in Cleveland.

B. Aerosol Condensation Algorithm With respect to modeling the icing tunnel and taking aerosol particles into consideration, the aerosol particles are introduced only when conditions exceed saturation. These aerosol particles are handled as a separate bin in the particle distribution of the cloud, and are treated like all the other particles in the cloud at this point. Condensation and deposition on these aerosol particulates occur through diffusion and the rate by which mass condenses on these aerosol particles and other cloud particles follow the same rate as described in Eq. (12) by Bartkus et al. The initial mass of the aerosol particle is the mass of an ice or water particle with a diameter of 0.04 μm (40 nm). The change in velocity of the aerosol particle follows as described in Eq. (13) by Bartkus et al. The initial aerosol particle velocity is approximated to be 99.99% of the air velocity at the moment when the aerosol is introduced. The energy balance equations as expressed in Eq. (14) and Eq. (15) by Bartkus et al., describe the aerosol particulates' change in American Institute of Aeronautics and Astronautics temperature and change in melt fraction, respectively. The initial aerosol temperature is assumed to be the local wet-bulb temperature when the particulates are first introduced. The wet-bulb temperature dictates whether the water O vapor will condense or deposit. If the wet-bulb temperature is above 0 C, then the water vapor will initially condense O as liquid water, and if it is at or below 0 C, then it will initially deposit as ice. It should be understood that the initial aerosol particle merely serves as a seeding agent of condensation (or deposition), and once this phase change occurs on these particulates, the aerosol particulate becomes a water/ice particle for computational purposes. Calling it an aerosol particle merely distinguishes it as an ice/water particle that nucleated due to supersaturation, unlike the other particles which were introduced from the spray nozzles.

IV. Simulation Results and Comparisons with Experimental Data Several parameter subscripts need to be defined, which will apply for the remainder of the document. Subscripts of ‘ 0 ’ and ‘ s ’ refer to a total or static value , respectively. A subscript of ‘ i ’ refers to a parameter at the tunnel inlet, whereas a subscript of ‘ e ’ refers to the tunnel exit. A subscript of ‘ exp ’ refers to experimentally measured data, and ‘ sim ’ refers to model simulation results. Variables with a subscript of ‘ calc’ ’ refer to parameters that are calculated using isentropic relations. A subscript of ‘aero’ refers to a simulation where the aerosol condensation subroutine was used, whereas a subscript of ‘none’ refers to the absence of the aerosol condensation subroutine, and condensation only occurs on existing spray particles.

A. Aerosol Condensation Simulation Sample Results A sample simulation where the aerosol condensation subroutine was implemented due to supersaturation is shown in Fig. 4. For this simulation the total inlet air temperature was set to T = 10.0 C, the total inlet relative humidity 0,i was RH = 76.9%, the total inlet pressure was P = 87.2 kPa, and the exit air velocity was U = 135 m/s. The injected 0,i 0,i e cloud had an inlet liquid water content of LWC = 7.1 g/m , and a particle size distribution where the inlet median i volumetric diameter was MVD = 15 μm. As previously mentioned, an aerosol particle concentration of 22,000/cm i with an initial size of 0.04 μm (40 nm) is used for the aerosol condensation subroutine. Figure 4A shows the change in particle size for two particles along the icing tunnel axis. The solid black l ine represents the cloud’s inlet MVD bin, the dashed red line represents the aerosol particles on which water vapor condensed. An outline of the tunnel radial geometry is shown for reference in Fig. 4A. Figure 4B shows the static relative humidity, RH , as well as the mass s mixing ratio, 𝜔 , (humidity) along the axis of the PSL icing tunnel for the same simulation. The mass mixing ratio initially increased due to cloud evaporation and hit a maximum at 3.0 m, but then decreased in value to the end of the tunnel due to constant condensation during this stretch. The relative humidity increased initially until it plateaued near saturation levels in the plenum. As tunnel geometry decreased, static temperature decreased which caused relative humidity to continue increasing, despite an absolute humidity decrease. RH continued to increase until constant tunnel s geometry was reached, where RH was driven down in value due to condensation and steady static temperature to the s exit of the tunnel. In Fig. 4A, aerosol particles are introduced when the air exceeds saturation at 3.0 m. Both particle size bins increase in size when condensation begins, but the aerosol particle bin increases more due to the high surface area to volume ratio compared to the 15-μm particle bin. The aerosol particle reached 2.8 μm in diameter by the tunnel exit.

B. Aerosol Condensation Simulation Parametric Sweeps Literature shows a wide range in aerosol number density and size distribution that can exist for different times and geographical locations. Therefore, parametric sweeps were run to investigate the sensitivity of these two parameters on condensation when vapor conditions exceed saturation. The conditions run for these parametric sweeps were the same as in Fig. 4, save the parameter that was varied. Figure 5A shows the static relative humidity at the tunnel exit where the aerosol number density was varied. The five number density values were taken from previously mentioned 32 29 30 33 -3 31 references and are 3,100, 6,200, 12,000, 22,000, and 100,000 cm . A fairly sensitive relationship exists between exit relative humidity and the aerosol number density. With a greater number of aerosol particles, a greater total surface area exists on which condensation can occur. Figure 5A suggests that a greater supersaturation value exists at the tunnel exit when the air is cleaner and contains fewer aerosol particulates. Figure 5B shows static relative humidity at the tunnel exit where the aerosol initial diameter was varied. As aerosol size increased, it provided greater surface area on which water vapor can condense, hence the decrease in final RH with increasing aerosol size. This s relationship, however, is not strong for the range of aerosol diameters that were run. The final RH values are provided s next to each point in Fig. 5A and 5B for clarity.

American Institute of Aeronautics and Astronautics C. Simulation Comparison with Supersaturated Experiments Experiments were run to observe behavior when conditions exceed saturation in the PSL tunnel. Four tests were run with elevated RH values in the plenum which calculated to RH values near or above 100% at the tunnel exit. The s target conditions for these four tests were T = 10.0 C, P = 87.2 kPa, and U = 135 m/s. The resulting target static 0,i 0,i e values at the tunnel exit were T = 0.8 C, and P = 78.0 kPa. When the spray nozzles were activated, the inlet cloud s,e 0,i 3 4 values were LWC = 7.1 g/m , and MVD = 15 μm. Struk et al. describes the method used to estimate LWC for the i i effective area in which the cloud exists. Figure 6 shows a view of the PSL tunnel, forward-looking-aft, from the spray bar system for these four test conditions. Figure 6 provides the experimentally measured total relative humidity at the tunnel inlet, RH and the isentropically calculated static relative humidity at the tunnel exit with the spray off, 0,i,exp RH . In Fig. 6A through 6C, the spray nozzles are not activated, but with increasing RH in the plenum a s,e,calc 0,i condensation cloud becomes visible in Fig. 6B and 6C, near the exit of the tunnel. The value of RH , in Fig. 6B s,e,calc and 6C (and 6D), with the spray nozzles turned off are above 100%, and suggest potential supersaturation. The condensation cloud in Fig. 6C is denser and more visible than in Fig. 6B, as is expected since the potential level of supersaturation, and hence condensation, at the tunnel exit is greater in the former. The conditions in Fig. 6D are nearly identical to those in Fig. 6C, except the spray nozzles have been activated. The slight increase in RH , as 0,i,exp noted in the images, is due to experimental tunnel drift and rounding. The condensation cloud has been combined with the cloud ejected from the spray nozzles in Fig. 6D. The presence of a condensation cloud in Fig. 6B and 6C, when the spray nozzles are not activated, suggests that there is heterogeneous condensation occurring on aerosol particles.

Given that homogeneous condensation requires very high levels of supersaturation, it is believed that heterogeneous condensation is occurring in the tunnel. The test conditions labeled in Fig. 6A-6D are condition numbers (Cond#) 101, 102, 103, and 105, respectively.

These RH values at the tunnel exit, as shown in Fig. 6, are calculated using the mass mixing ratio measurement s,e,calc at the inlet in the plenum, 𝜔 , along with the static temperature and static pressure values that result at the tunnel i,exp exit. The measured total temperature and total pressure in the plenum along with the measured velocity at the tunnel exit, are used to calculate static temperature and pressure, using isentropic relations. The value of RH , does not s,e,calc account for any condensation and is simply the static relative humidity value at the tunnel exit given total inlet conditions and an exit velocity, in the absence of a cloud.

Temperature and humidity measurements were made at the tunnel exit at the centerline using a rearward-facing probe (Fig. 3A and Fig. 6). Table 1 shows relevant experimental centerline data regarding Cond# 102, 103 and 105, along with simulation predictions. Column 3 of Table 1 shows the experimentally measured total temperature at the tunnel inlet, T , and Column 4 provides the isentropically calculated static temperature at the tunnel exit with the 0,i,exp spray off, T . Similarly, Columns 5 and 6 show the measured total relative humidity at the tunnel inlet, RH , s,e,calc 0,i,exp and the isentropically calculated static relative humidity at the tunnel exit with the spray off, RH . Column 7 s,e,calc represents a theoretical mass mixing ratio at 100% relative humidity at the tunnel exit using T and P values, s,e,calc s,e,calc and is named 𝜔 100%RH . Column 8 shows the experimentally measured mass mixing ratio in the tunnel plenum, 𝜔 i,exp .

As can be seen, 𝜔 is greater than 𝜔 , and therefore the flowing air has the capacity to become supersaturated.

i,exp 100%RH Column 9 shows the experimentally measured mass mixing ratio at the tunnel exit, 𝜔 , which was lower than 𝜔 e,exp i,exp for all three test conditions. Between Cond# 102 and 103, the amount of decrease in 𝜔 between the exit and inlet was greater for the more highly supersaturated test condition, suggesting more phase change from the gaseous phase to the liquid phase. This was corroborated visually as a stronger condensation cloud was present in Cond# 103. As the spray nozzles were turned on, going from Cond# 103 to Cond# 105, 𝜔 increased slightly. Some of the cloud likely e,exp evaporated in the plenum, pushing the absolute humidity slightly higher, but not by much. Since RH was already 0,i,exp elevated at around 77%, there was not much capacity to evaporate and add to the gas phase before the tunnel began to contract, at which point the supersaturated air began to condense. A complicated interaction exists that includes evaporation, condensation, changing tunnel geometry, aerosol particulates and cloud particles before a final 𝜔 e,exp value is measured for Cond# 105.

American Institute of Aeronautics and Astronautics Table 1. Experimental measurements and simulation predictions for three test conditions where saturation was exceeded.

1 2 3 4 5 6 7 8 9 10 1 1 Cond Spray T 0, i,exp T s,e, calc RH 0, i,exp RH s,e, calc 𝜔 100%RH 𝜔 i ,exp 𝜔 e ,exp 𝜔 e,sim, none 𝜔 e,sim, aero O O # On/Off C C % % g/kg g/kg g/kg g/kg g/kg 102 Off 10.9 1.8 64 107 5.61 6.01 5.99 6.01 6.00 103 Off 10.1 1.1 76 127 5.34 6.87 6.35 6.87 6.79 105 On 10.0 1.0 77 128 5.30 6.81 6.42 7.15 6.94 Columns 10 and 11 show the simulation predictions for the mass mixing ratio at the tunnel exit. Simulation results in Column 10 represent the previous model version where condensation occurred only on existing spray particles and are represented as 𝜔 . For the two test conditions where the spray was not activated (Cond# 102 and 103), 𝜔 sim,none sim,none equaled 𝜔 since no condensation was possible in the absence of spray particles using the previous model algorithm.

i,exp A net gain was predicted in 𝜔 at the exit by 0.34 g/kg when the spray was activated in Cond# 105, whereas the experiment measured a net loss of 𝜔 at the exit by about -0.39 g/kg. Column 11 shows simulation results that includes the aerosol condensation subroutine, where condensation occurred on both aerosol particles and existing spray particles, and is represented as 𝜔 . An aerosol particle concentration of 22,000/cm with a uniform aerosol sim,aero distribution size of 0.04 μm (40 nm) was used for the aerosol condensation subroutine. When simulations were run using the aerosol condensation subroutine, simulation 𝜔 values more closely matched the experimentally measured e values, but there still exists a fair difference. The simulation still predicted a net gain in 𝜔 of 0.13 g/kg at the exit, but by a smaller amount when compared to the simulations that did not use the aerosol condensation subroutine (+0.34 g/kg). It is worth noting that a higher aerosol concentration would have predicted 𝜔 values that more closely e 3 3 matched experimental values. For example, aerosol concentrations of 250,000/cm and 138,000/cm would have produced simulation results that matched 𝜔 for Cond# 103 and 105, respectively. There are other physical e,exp mechanisms that are not considered in this simple aerosol condensation model, such as the effects of charged particles that may enhance diffusion, and hence condensation rates. Nor does the model differentiate in the composition of the aerosol particulates, which may promote the rate of condensation as well. The model simply introduces the aerosol particulates when vapor saturation is exceeded and treats it like a non-charged water droplet or ice particle.

D. Simulation Comparisons to Plenum Relative Humidity Sweeps Several tests were conducted during the first ice-crystal icing physics tests at PSL. This paper will focus on four series of tests that systematically varied plenum relative humidity ( RH ) while maintaining the other aero-thermal 0,i and spray parameters constant. Comparisons between experimentally measured data and simulation results will focus on the aero-thermal changes due to the presence of the cloud and the cloud characteristics at the tunnel exit (experiment test section). The aerosol condensation subroutine was used for simulations where vapor content exceeded saturation.

Specific details of the experimental conditions and measurements can be found in Struk et al. The following paragraphs describe experimental and simulation aspects important to understand for these series of tests and predictions.

These four RH sweeps examined the effects of humidity on the freeze-out of the cloud. By lowering RH , 0,i 0,i evaporation and evaporative cooling were enhanced, which helped promote freeze-out of the cloud. The four humidity sweeps were performed at two different target exit velocities of U = 85 m/s and 135 m/s, and two different inlet e particle size distributions, where MVD = 15 μm and 50 μm. The other facility conditions were kept constant with a i target total air temperature at the inlet of T = 7.2 C, and target total pressure of P = 44.6 kPa. The air velocity in 0,i 0,i the plenum was low (~ 2-3 m/s) and therefore the conditions at the inlet are approximated as total conditions (e.g. T , 0,i P , RH ). The isentropic calculations of static air temperature and static pressure at the tunnel exit, where 0,i 0,i O U e = 85 m/s, were T s,e,off = 3.7 C and P s,e,off = 42.8 kPa, respectively. When U e = 135 m/s, the isentropic values at the O 0 tunnel exit became T = -1.8 C and P = 39.9 kPa. The target inlet water temperature was T = 7.2 C with s,e,off s,e,off water,i a target injected total water content of TWC = 7.0 g/m (or LWC since it is initially a liquid spray). Assuming no i i 3 3 water lost to evaporation, TWC = 6.8 g/m and TWC = 6.5 g/m for U = 85 m/s and U = 135 m/s, respectively. The e e e e volume of air expands as a fixed mass of air travels from the higher pressure inlet to the lower pressure tunnel exit, which is the reason for this decrease in TWC calculation when maintaining this momentary no evaporation assumption.

American Institute of Aeronautics and Astronautics For the tests conducted in these four sweeps, the total wet-bulb temperature, Twb , and static wet-bulb temperature, O Twb s , were within a few degrees of 0 C. Wet-bulb temperature is the temperature of a wet adiabatic surface undergoing evaporation (or sublimation of ice). Twb is the resulting temperature when convective heat transfer is balanced by evaporative cooling of a water particle and is a function of T , RH , and P . When wet-bulb temperatures O are at or below 0 C, a cloud that starts off as liquid water droplets can begin to partially freeze, creating a mixed-phase cloud, and can completely glaciate to an ice-crystal cloud if sufficiently cooled.

Radial variations within the tunnel existed during testing. It can be seen in Fig. 2 that the spray area patterns covered a fraction of the total area, causing radial variations within the tunnel. Figure 7 shows time-averaged tomography data that depict the cloud concentration near the tunnel exit. Figure 7A shows the tomography data for a test where U = 85 m/s and Fig. 7B shows a higher velocity test where U = 135 m/s. The tomography data in Fig. 7 e e suggest that there is some spreading of the cloud as the cloud flows from the tight nozzle pattern shown in Fig. 2 to the tunnel exit. The tomography shows that the radial spread is less for the higher velocity test as the cloud is concentrated more towards the center. In general, the spray pattern configurations created a central core region that contained the cloud, and an annular ring region that was cloud-free, but some mixing between the two regions occurred as the cloud traveled axially down the tunnel.

The five squares in Fig. 7A represent the locations of various instruments. The IKP2, CDP, CIP, and Multiwire instruments were measured at the center labeled as 2 in Fig. 6A. Humidity Probe “A” and the commercial total air temperature probe were located at locations labeled as 1 and 4 respectively, but the data from those instruments were not used in this paper. The rearward-facing temperature probe (RFTP) and Humidity Probe “B” were located at the squares labeled 3 and 5, respectively for the majority of test configurations and these measurements were used for comparison with model predictions. The temperature and humidity probes used for comparison in this paper were offset 0.25 m (9.7 in) from the centerline. Since a radial variation exists, a method was developed to approximate centerline humidity and temperature values from the measurements made by the offset temperature and humidity probes, and is described by Struk et al. The particle size measuring PDI and HSI probes made measurements located approximately 0.15 m (6 in) above or below the centerline at the test section. Centerline values were not extrapolated for the PDI and HSI probes. Aero-thermal and cloud measurements were desired at the center of the tunnel exit because that was the focus area of the ice accretion tests. The 1-D model was used to approximate conditions at the centerline of the tunnel, understanding that these are predictions that did not consider mixing between the central cloud-filled region and the cloud-free annular region of the tunnel.

Several tests were run for any one particular condition with various instruments on different days characterizing different aspects of the cloud. It should be noted that while small variations between tests may have existed using this serial testing technique, all tests for a given test condition were approximated as equal and collectively the measurements from the different instruments characterize the same cloud.

Figures 8 through 11 all have the same structure and compare experimentally collected data with simulation predictions at the tunnel exit for a range of inlet relative humidity values (x-axis). Graph A shows the change in total air temperature at the tunnel exit, ΔT , when the cloud was activated (the reader should note that the left axis of this 0,e graph is negative). The change, Δ , refers to the difference in a condition between the presence of an icing cloud and no cloud (spray-on minus spray-off) at the tunnel exit. Similarly, graphs B and C show the change in mass mixing ratio, Δ 𝜔 , and change in total wet-bulb temperature, ΔT wb , at the tunnel exit as the cloud was activated. Graph D e 0,e of each figure shows the melt fraction of the cloud at the tunnel exit, 𝜂 , where a value of 1 indicates a fully liquid e droplet, and a value of 0 represents a fully glaciated ice particle. Struk et al. describes the method by which the melt fraction was calculated, using experimentally collected IKP2 data and Multiwire data. Graph E in all figures shows the total water content of the cloud at the tunnel exit, TWC , and graph F shows the cloud’s median volumetric diameter e at the tunnel exit, MVD . Simulation predictions are marked as empty, red squares and experimental measurements e are represented by solid, black circles. Graph F displays multiple experimentally measured MVD values. A black, e solid circle refers to the CDP and/or CIP particle size measurements. For clouds composed of smaller particles, the CDP captured the entire spectrum of particle sizes, but clouds composed of larger particles required the combination of both CDP and CIP instrument measurements. A solid, green triangle represents measurements made by the PDI, and the blue crosses represent particle size measurements made by the HSI. A simulation was run for every test condition in the RH sweeps and all simulation results are plotted in every graph, but not every type of experimental measurement was made for every test condition, and therefore will appear as missing points in a graph. For example in Fig. 9E, an experiment total water content measurement with the IKP2 instrument was not made at RH = 10%, 0,i and is obviously not plotted, but the simulation prediction at that test condition is plotted in that graph.

American Institute of Aeronautics and Astronautics 1. RH Sweep #1 Figure 8 compares experimental data with simulation predictions for tests where U e = 85 m/s, and MVD i = 15 μm as inlet relative humidity, RH , varied from 20% to 60%. Experimental data in Fig. 8A show a slightly larger decrease 0,i in ΔT as RH was decreased. More evaporation of the cloud occurred at lower RH levels, increasing the drop in 0,e 0,i 0,i air temperature. This greater amount of evaporation at lower RH levels is reflected in the experimental data in Fig.

0,i 8B, where the change in mass mixing ratio was greater as RH decreased. Simulation predictions shows this trend for 0,i ΔT , and Δ 𝜔 just more pronounced than experimental data. These trends with respect to RH are expected because 0,e e, 0,i there is greater evaporation capacity at low RH conditions. When evaporation occurs, energy is removed from the 0,i air and transferred to the water particles for the phase change transition to vapor. Simulation over-predicted the amount of temperature decrease and humidity increase for this humidity sweep.

Whereas ΔT can be large for some tests, the change in total wet-bulb temperature at the tunnel exit, Δ Twb , 0,e 0,e was not as large. When evaporation occurs, air temperature decreases, but the increase in vapor content acts to balance the Twb value. There is, however, a trend that can be seen in both experimental data and prediction in Fig. 8C such that the change in total wet-bulb temperature at the tunnel exit, Δ Twb increased as RH decreased.

0,e 0,i This trend in Δ Twb with respect to RH is attributed to the energy added from the cloud to the air. For an 0,e 0,i adiabatic system containing an unsaturated parcel of air and water particles (a cloud) that is not in equilibrium, the cloud mass and the air mass will thermally interact until the two masses reach a thermodynamic steady-state. During the interaction, the air will become fully saturated ( RH = 100%) as the water particles evaporate, and the air parcel temperature will decrease to the wet-bulb temperature with the air parcel providing the energy to evaporate the cloud.

Any remaining water particles will reach this steady-state Twb as well. If the initial temperature of the water particles equals the Twb , then Twb will remain constant during the thermal interaction and this temperature is known as the adiabatic saturation temperature. If the initial water particle temperature is greater than the Twb , then the additional water sensible heat will cause a final steady-state Twb that will be greater than the adiabatic saturation temperature.

In addition, if the water droplets undergo a phase change from liquid water to solid ice, then that additional latent heat will increase the final steady-state Twb as well. For the experiments conducted, the inlet water temperature of the O liquid spray was approximately T = 7 C, which was greater than the inlet total wet-bulb temperature for the all water,i O O tests in this sweep which ranged between -4 C < Twb < 3 C, depending on RH . As the water particles cooled to 0,i 0,i the wet-bulb temperature, the sensible energy was transferred to the air. Additional heat was transferred to the air if freezing of the water droplets occurred, with more latent heat being transferred to the air when a greater fraction of the cloud froze. When RH was low, Twb was low, which caused more sensible energy to be transferred to the air, 0,i 0,i and a larger latent energy transfer to the air as a greater fraction of the cloud froze. These larger transfers of sensible and latent heat to the air caused a greater increase in Twb at lower RH test conditions .

0,i For this series of tests, the particles within the cloud are small, with MVD = 15 μm and the largest particles in this i distribution being about 50 μm in diameter. Small particles respond to changes in temperature quickly due to the high surface area to volume ratio. It is for this reason that there is a small window to achieve a mixed-phase cloud when the particles within the cloud are small. For example, a liquid cloud that is exposed to conditions where the wet-bulb O temperature is below 0 C will decrease in temperature until phase change occurs. The smallest particles will reach mixed-phase first, followed by the largest particles and the particles follow this order through complete glaciation.

With clouds composed of small particles this transition can occur quickly, and the greater the difference between Twb and T the more quickly the cloud will go through these stages. The experimental data in Fig. 8D suggest that the water window for a mixed-phase cloud ranges from about RH = 25% (fully glaciated) to about RH = 60% (fully liquid), 0,i 0,i for the given fixed temperature, pressure and air velocity. A series of three consecutive tests showed that the melt fraction can change drastically for a very small RH window, as 𝜂 increased from 0.23 to 0.89 when RH varied 0,i e 0,i from 36% to 40%. The simulation points in Fig. 8D do not match the melt fraction precisely, but the range for mixed-phase matches relatively well.

The simulation predictions of TWC depicted in Fig. 8E are noticeably larger than experimentally measured values.

e With simulation predicting greater amounts of evaporation than what was measured experimentally, the simulation prediction of TWC should be lower than the experimentally measured TWC , however this is not the case. It is likely e e the reason for this discrepancy is that there is mixing between the cloud-filled center region and the cloud-free annular regions of the tunnel. The water content spread out into the cloud-free annular ring of the tunnel as the cloud flowed downstream, which reduced the cloud concentration at the center of the tunnel exit where measurements were made.

The 1-D model does not capture this radial movement and highlights a limitation in the model when radial variations exist. This mixing and spreading of the cloud in the radial direction would also explain the model’s over-predicted changes in temperature and vapor content. For example, a volume of air that has decreased in temperature due to evaporation that partially mixes with a volume of air that has not experienced any evaporative cooling will result in an overall decrease in temperature that is not as great if no mixing was involved. The same concept applies with vapor American Institute of Aeronautics and Astronautics content mixing between a more humid volume of air, which was a result of cloud evaporation, and a less humid volume of air in the annular ring of the tunnel.

Figure 8F shows that experimentally measured MVD of the cloud at the tunnel exit measured to be about MVD = 20 μm, irrespective of RH . This experimentally measured value was consistent between the two particle e 0,i measuring instruments, the CDP and the PDI. With an inlet value of MVD = 15 μm, this measurement would suggest i a slight increase in MVD after portion of the cloud evaporated. Simulation also showed an increase, but not as much as what was measured experimentally with a prediction of MVD = 16 μm. This trend of MVD increasing may appear e to be counterintuitive, but can be explained. During evaporation, the smallest particles completely vaporized, leaving just the largest particles behind. The larger particles also become smaller from evaporative mass losses, but the overall effect mathematically resulted in larger MVD values. In addition, particles that glaciated expanded in volume, which e also played a role in increasing the cloud MVD value at the tunnel exit. This increase in MVD due to evaporation is e e only momentary. If the cloud were allowed to continue evaporation to completion, the MVD will eventually decrease e and reach zero.

2. RH Sweep #2 In RH Sweep #2, the parameters were similar to RH Sweep #1, except the inlet size of the cloud particles were larger. Figure 9 compares experimental data with simulation predictions for tests where U = 85 m/s, and e MVD = 50 μm as RH varied from 0% to 61%. Figures 9A, 9B, and 9C show the same trends and nearly the same i 0,i values in ΔT , Δ 𝜔 , and Δ Twb , as in RH Sweep #1 for both experiment and simulation. Again, simulation over- 0,e e 0,e predicted the amount of temperature decrease, humidity increase, and wet-bulb temperature increase for this humidity sweep. Both experiment and simulation did not show a great deal of sensitivity for ΔT , Δ 𝜔 , and Δ Twb with respect 0,e e 0,e to cloud particle size ( MVD ).

Figure 9D shows that simulation predicted melt fraction very well for the larger particle size cloud. It can be seen that a fully glaciated cloud did not occur even at RH = 0%, which was supported by simulation prediction. Despite 0,i O an inlet total wet-bulb temperature of Twb = -7.5 C, when RH = 0%, the largest particles in the cloud did not have 0,i 0,i enough time to fully glaciate. The largest particles of the cloud at the inlet for this series of tests are estimated to be 330 μm in diameter.

Simulations predicted TWC values that matched experiment more closely for this larger particle size cloud than e for the smaller particle size cloud in RH Sweep #1, and can be seen in Fig. 9E. As some mixing between the cloud region and the cloud-free annular region still occurred, the larger particles may have acted in a more ballistic fashion and may have remained near the center of the cloud or even concentrated towards the centerline of the tunnel. Again, this 1-D model does not address radial particle movement, only axial movement.

Figure 9F shows the MVD at the tunnel exit for various particle size measuring instruments, along with the simulation prediction. Both the PDI and HSI measured values that remained relatively constant for the range of RH 0,i tested. The two instruments did not measure much change from the inlet MVD = 50 μm value, as MVD ranged from i e 47 μm to 57 μm. The combined CDP and CIP measurements showed a significant increase in MVD e from the inlet size, as MVD nearly doubled. The CDP and CIP measurements showed a slightly greater increase in MVD as RH e e 0,i decreased. It is important to note that the CDP and CIP instruments were located at the tunnel centerline, whereas the PDI and HSI were located 0.15 m (6 in) off-center. A concentration of larger particles towards the center may explain this measurement discrepancy. Simulation predictions are bounded by the experimental measurements. Simulation also showed a negative correlation between MVD and RH , much like the combined CDP and CIP measurements.

e 0,i As mentioned earlier, when the conditions are more conducive for greater evaporation, MVD e will increase.

3. RH Sweep #3 For RH Sweep #3, the conditions were similar to RH Sweep #1, except the air velocity at the tunnel exit was greater. Figure 10 shows the comparison between experimental data with simulation predictions for tests where U = e 135 m/s, and MVD = 15 μm as RH varied from 25% to 51%. Experimental measurements of ΔT , as shown in i 0,i 0,e Fig. 10A, for this series of tests was greater than the slower air velocity tests in RH Sweep #1. In this series the model under-predicted the measured ΔT 0,e value. In fact the simulations predicted a smaller drop in total air temperature at the tunnel exit as compared to the slower air velocity test predictions of RH Sweep #1. The reason for this smaller ΔT prediction is because the residence time of the cloud in the tunnel was shorter at this higher velocity, and therefore 0,e allowed for less evaporation time. It is not clear why the experimental ΔT value was greater for these tests. One 0,e possible explanation is that there was less mixing between the core cloud and non-cloud annular regions, maintaining more of the air temperature change in the core where the measurements were taken. Another possibility is that the method by which temperature change was extrapolated to the centerline overestimated the temperature change. It can be seen that squares numbered 3 and 5 in Fig. 7B, the locations of the temperature and humidity probes respectively, American Institute of Aeronautics and Astronautics are in the periphery of the cloud. Extrapolation of temperature change with small measurements may have greater error associated with the approximation.

Comparison of Δ 𝜔 between simulation and experiment can be seen in Fig. 10B. Predictions of Δ 𝜔 , while still e e over-predicting the amount of vapor change compared to experimental data, matched more closely for this series of tests than the simulation-experiment comparison with RH Sweep #1. Compared to RH Sweep #1, both the predicted and experimental Δ 𝜔 value were lower. This trend of smaller Δ 𝜔 values, when compared to the lower air velocity e e series of tests, is expected as residence time of the cloud in the tunnel was shorter at this higher velocity, which allowed for less evaporation time.

Experimental Δ Twb , as shown in Fig. 10C, was slightly negative and relatively insensitive to RH . Predictions 0,e 0,i of Δ Twb did not match experimental measurements well for this humidity sweep. Again, the prediction trend was 0,e that Twb increased with cloud activation, and that there was greater increase in Twb for lower RH , as was the 0,e 0,e 0,i trend in the two previous series of tests.

The comparison between simulation and experiment with respect to melt fraction at the tunnel exit, as seen in Fig. 10D, matched reasonably well at lower RH values, however some discrepancy exists at elevated RH values.

0,i 0,i Experimentally, the range of RH to generate a mixed-phase cloud at the tunnel exit at this higher air speed and small 0,i particle size was fairly small.

The comparisons between simulation and experiment regarding TWC , and MVD are very similar to the e e comparisons for RH Sweep #1, where MVD = 15 μm as well. These comparisons are shown in Fig. 10D and 10F, i respectively. The trends and discrepancies that were discussed in RH Sweep #1 for TWC e , and MVD e apply for this series of tests as well. Slightly less evaporated due to the cloud’s shorter residence time, thus slightly higher TWC e values were measured and predicted as compared to the slower air velocity series of tests. Cloud particle size was only measured for one test, and both particle size measuring instruments measured roughly MVD = 20 μm, whereas e simulation predicted about MVD = 15 μm for the entire RH range.

e 0,i 4. RH Sweep #4 Figure 11 compares experimental data with simulation predictions for RH Sweep #4 where U e = 135 m/s, and MVD = 50 μm as RH varied from 10% to 58%. Figure 11A shows prediction of ΔT matched experiment fairly i 0,i 0,e well, except at the lowest RH setting. The trend of larger decreases of T as RH decreases exists again, due to the 0,i 0,e 0,i greater amount of evaporation that occured, and hence greater amounts of heat transferred from the air.

Simulation for this series of tests over-predicted the amount of evaporation, however the trend was correct, and can be seen in Fig. 11B. Compared to the slower velocity tests with same cloud MVD in RH Sweep #2, the amount i of evaporation was lower for this series of tests, which is expected for the shorter cloud residence time. This reduced amount of evaporation compared to RH Sweep #2 was experimentally measured and predicted.

Experiment and simulation of Δ Twb , is similar to the other high velocity series of tests of RH Sweep #3.

0,e Experiment showed a slight decrease in Twb as the spray cloud was activated, whereas simulation predicted an 0,e increase in Twb 0,e with larger increases at lower RH 0,i conditions.

Figure 11D shows how the model consistently predicted slightly lower melt fraction values than what was measured experimentally. Complete glaciation of the cloud did not occur even at the lowest RH setting. A 0,i combination of large particles within the cloud and short residence time were responsible for the cloud not fully glaciating.

Interesting experimental TWC measurements were made during this series of tests. As mentioned earlier, the inlet e total water content at the tunnel inlet was estimated to be TWC i = 7.0 g/m . Assuming no water lost to evaporation, and accounting for air expansion, TWC = 6.5 g/m was calculated at the tunnel exit. This value represents the e maximum value possible, but Fig. 11E shows TWC measurements greater than 6.5 g/m . One possible explanation e that has been mentioned earlier is that the larger particles may have acted in a more ballistic fashion and may have concentrated towards the centerline of the tunnel. Again, this 1-D model assumes uniform radial distribution and therefore does not model any radial variations.

The prediction and simulation values of MVD e for this series of tests are similar to those of the slower air velocity test where MVD = 50 μm of RH Sweep #2. Figure 11F shows that the MVD measurement by the HSI i e instrument did not change much from the inlet MVD value for the range of RH tested. Similar to RH Sweep #2, the i 0,i combined CDP and CIP measurements showed a sizeable increase in MVD with a slightly larger increase in MVD e e for lower RH test conditions. Simulation similarly predicted a negative correlation between MVD and RH . Again, 0,i e 0,i all predicted MVD values were bounded by the experimental measurements. Since less evaporation occurred at higher e velocity, the predicted MVD e is slightly smaller than the slower air velocity MVD e predictions of RH Sweep #2.

American Institute of Aeronautics and Astronautics V. Discussion A thermal model was written to better understand the observed changes in test conditions at the PSL icing test facility. Modifications to the model were made to address conditions that lead to vapor content exceeding saturation.

Whereas many trends in data were matched with what was measured experimentally, the simulations often over- predicted changes in air temperature and humidity content. These discrepancies are likely attributed to the multiple approximations that were required in order to make comparisons between experimentally collected data with simulation prediction.

Most notable of the approximations made, the 1-D model does not account for the 2-D aspects that were inherent during these experiments. The nozzle configuration that was used during experiments created a non-uniform cross section in the radial direction. In general, the center core of the cross section was cloud-filled, but the annular region was cloud-free. As simulations were run to approximate conditions at the center of the tunnel, where a 1-D model would be most applicable, there was however, evidence of mixing between the two regions. This mass and thermal mixing likely diffused differences in vapor content, cloud concentration and temperature between the cylindrical core and annulus. The 1-D simulations do not capture these 2-D interactions. In addition, three traversing tests showed humidity and temperature measurements varied radially, so a method was required to extrapolate temperature and humidity to the centerline from measurements made offset from the center. This method, as described by Struk et al.

can be a source of experimental temperature and humidity uncertainty, contributing to the comparison discrepancy.

Trajectories of particles may have a 2-D component that is not captured with the 1-D model, and may account for some of the TWC discrepancies between experiment and prediction. As mentioned, cloud particles mix into the annular cloud-free region of the cloud, dispersing the TWC concentration. The geometric contraction in the tunnel, however, may cause the larger particles to concentrate towards the center of the tunnel, as larger particles are more ballistic due to greater inertial forces. The combined effects of larger particles concentrating towards the center of the tunnel, and greater diffusion of smaller particles into the annular region via mixing, may explain some of the experimentally measured data. Tomography showed the cloud spreading further into the annular region, but the IKP2 located at tunnel center measured TWC values at the tunnel exit that were greater than inlet TWC values. This increased TWC measurement occurred for some experiments that were conducted with larger particles, suggesting that some of the cloud converged towards the center of the tunnel. This concentration of larger particles at the tunnel center is supported by particle size measuring instruments as well. For clouds composed of larger particles, instruments located at the tunnel centerline measured larger MVD values, compared to instruments located offset from the centerline. The e method used to estimate TWC , as described by Struk et al., can be another source of experimental TWC uncertainty as well.

VI. Conclusion The purpose of the model developed is to better understand the complex interactions between the test parameters and have greater confidence in the conditions at the test section of the PSL tunnel. As such, the model can be used to help explain the observed changes in test conditions for PSL icing tests by fully coupling the conservation of mass and energy equations between the cloud particles and flowing air mass. The model aided the first fundamental physics of ice-crystal ice accretion tests in 2016 by guiding the development of the test matrix and helping to identify the mixed-phase parameter space. Experimental aero-thermodynamic measurements from this testing were used to compare with model simulations.

A subroutine was added to the model that allowed for the option to condense on existing aerosol particulates when water vapor content in the tunnel exceeds saturation. When compared to experiments where vapor saturation was exceeded, predictions with the aerosol condensation subroutine fared better than without it. The model, however, still over-predicted the amount of vapor at the tunnel exit. The concentration of aerosol particulates used in the model was chosen from literature that best represented the atmospheric conditions at the PSL icing tunnel, since direct measurements of aerosol concentration was not made. When a greater aerosol concentration was used in the model, more vapor condensed out of the air, matching experimental measurements at the tunnel exit more closely.

Comparisons were made between prediction and experimental data for four series of tests where relative humidity was systematically varied such that the cloud at the tunnel exit ranged from liquid to fully glaciated. Two air velocities and two cloud size conditions created the matrix of four series of tests. When compared to experimental data, the model matched the trend of greater temperature and humidity changes at low humidity conditions, and matched the cloud melt fraction fairly well. The model, however, often over-predicted the aero-thermodynamic changes. The two dimensional nature of the injected cloud within the tunnel may be a cause of the comparison discrepancy. In addition, the requirement to extrapolate temperature and humidity measurements to the centerline in order to make a comparison between simulation and prediction may have contributed to the discrepancy as well. The Discussion section, American Institute of Aeronautics and Astronautics highlighted the limitations of using a 1-D model to simulate experiments with 2-D aspects. Despite these limitations, the model provided useful comparisons.

Acknowledgments The authors wish to acknowledge the financial support for this work by the Advanced Aircraft Icing (AAI) Subproject of the NASA Advanced Air Transport Technology Project (AATT) under NASA's, Advanced Air Vehicles Program (AAVP).

References Mason, J. G., Strapp, J. W. , and Chow, P., “The Ice particle Threat to Engines in Flight,” 44th AIAA Aerospace Sciences Meeting and Exhibit , AIAA, Reno, NV, 2006, AIAA-2006-206.

Struk, P., Currie, T., Wright, W. B., Knezevici, D., Fuleki, D., Broeren, A. P. , Vargas, M., and Tsao, J., “Fundamental Ice Crystal Accretion Physics Studies”, SAE International Conference on Aircraft and Engine Icing and Ground Deicing , SAE, Chicago, IL, 2011, SAE-2011-38-0018.

Currie, T. C., Struk, P. M., Tsao, J., Fuleki, D., and Knezevici, D. C ., “Fundamental Study of Mixed -Phase Icing with Application to Ice Crystal Accretion in Aircraft Jet Engines,” 4th AIAA Atmospheric and Space Environments Conference , AIAA, New Orleans, LA, 2012, AIAA 2012-3035.

Struk, P. M., Ratvasky, T. P., Bencic, T. J., Van Zante, J. F., King, M. C., Tsao, J., and Bartkus, T. P., “An Initial Study of the Fundamental of Ice Crystal Icing Physics in the NASA Propulsion Systems Laboratory ,” 9th AIAA Atmospheric and Space Environments Conference , AIAA, Denver, CO, 2017 (submitted for publication).

Struk, P. M., Tsao, J., and Bartkus, T., “Plans and Preliminary Results of Fundamental Studies of Ice Crystal Icing Physics in the NASA Propulsi on Systems Laboratory,” 8th AIAA Atmospheric and Space Environments Conference , AIAA, Washington D.C., 2016, AIAA 2016-3738.

Struk, P. M., Bencic, T., Tsao, J., Fuleki, D, and Knezevici, D. C. “Preparation for Scaling Studies of Ice -Crystal Icing at the N RC Research Altitude Test Facility,” 5th AIAA Atmospheric and Space Environments Conference , AIAA, San Diego, 2013, AIAA 2013-2675.

Fuleki, D., Mahallati, A., Currie, T. C., MaCleod, J. D., and Knezevici, D. C., “Development of a Sensor for Total Temperature and Humidity Measurements under Mixed-Phase and Glaciate d Icing Conditions,” 6th AIAA Atmospheric and Space Environments Conference , AIAA, Atlanta, GA, 2014, AIAA 2014-2751.

Willbanks, C. E., and Schulz, R. J. , “Analytical Study of Icing Simulation for Turbine Engines in Altitude Test Cells,” AEDC - TR-73-144, 1973.

Hindmarsh, J. P., Russell, A. B., and Chen, X. D. , “Experimental and numerical analysis of the temperature transition of a suspended freezing wat er droplet,” International Journal of Heat and Mass Transfer, Vol. 46, No. 7, 2003, pg. 1199-1213.

Strub, M., Jabbour, O., Strub, F., and Bédéarrats, J. P. , “Experimental study and modelling of the crystallization of a water droplet,” International Journal of Refrigeration , Vol. 26, No. 1, 2003, pg. 59-68.

Wright, W. B., Jorgenson, P. C. E., and Veres, J. P., “ Mixed Phase Modeling in GlennICE with Application to Engine Icing ” , AIAA Atmospheric and Space Environments Conference , AIAA, Toronto, ON, 2010, AIAA 2010-7674.

Hauk, T., Roisman, I., and Tropea, C., “Investigation of the Melting Behavior of Ice Particles in an Acoustic Levitator,” th 11 AIAA/ASME Joint Thermophysics and Heat Transfer Conference , AIAA, Atlanta, GA, 2014, AIAA 2014-2261.

Villedieu, P., Trontin, P., and Rémi, C., “Glaciated and mixed - phase ice accretion modeling using ONERA 2D icing suite,” th 6 AIAA Atmospheric and Space Environments Conference , AIAA, Atlanta, GA, 2014, AIAA 2014-2199.

Bartkus, T. P., Struk, P. M., Tsao, J. and Van Zante, J. F. , “Numerical Analysis of Mixed -Phase Icing Cloud Simulations in the NASA Propulsion Systems Laboratory,” 8th AIAA Atmospheric and Space Environments Conference , AIAA, Washington D.C., 2016, AIAA 2016-3739.

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. Vol. 8, No. 1, 2015, pg. 15-32, SAE-2015-01-2155.

Oliver, M. J., "Validation Ice Crystal Icing Engine Test in the Propulsion Systems Laboratory at NASA Glenn Research Center," 6th AIAA Atmospheric and Space Environments Conference , AIAA, Atlanta, GA, 2014, AIAA 2014-2898.

Van Zante, J. F., and Rosine, B. M., "NASA Glenn Propulsion Systems Lab: 2012 Inaugural Ice Crystal Cloud Calibration," 6th AIAA Atmospheric and Space Environments Conference , AIAA, Atlanta, GA, 2014, AIAA 2014-2897.

Van Zante, J. F., Bencic, T. J., and Ratvasky, T. P. , “Update on the NASA Glenn Propulsion Systems Lab Ice Cloud Characterization,” 8th AIAA Atmospheric and Space Environments Conference , AIAA, Washington D.C., 2016, AIAA 2016-3897.

Lilie, L., Emery, E., Strapp, J. W., and Emery, J. "A Multiwire Hot-Wire Device for Measurement of Icing Severity, Total Water Content, Liquid Water Content, and Droplet Diameter,” 43rd AIAA Aerospace Sciences Meeting and Exhibit , AIAA, Reno, NV, 2005, AIAA 2005-859.

Struk, P. M., Rigby, D. L., and Venkataraman, K. "A Thermal Analysis of a Hot-Wire Probe for Icing Applications," 6th AIAA Atmospheric and Space Environments Conference , AIAA, Atlanta, GA, 2014, AIAA 2014-2331.

Rigby, D. L., Struk, P. M., and Bidwell, C. S., "Simulation of fluid flow and collection efficiency for an SEA multi-element probe," 6th AIAA Atmospheric and Space Environments Conference , AIAA, Atlanta, GA, 2014, AIAA 2014-2752.

American Institute of Aeronautics and Astronautics Strapp, J. W., Lilie, L. E., Ratvasky, T. P., Davison, C. R. , and Dumont, C., “Isokinetic TWC Evaporator Probe: Development of the IKP2 and Performance Testing for the HAIC-HIWC Darwin 2014 and Cayenne Field Campaigns," 8th AIAA Atmospheric and Space Environments Conference , AIAA, Washington D.C., 2016, AIAA 2016-4059.

King, M. C., Bachalo, W., and Kurek, A. "Particle Size Measurements from the First Fundamentals of Ice Crystal Icing Physics Test at the NASA Propulsion Systems Laboratory," 9th AIAA Atmospheric and Space Environments Conference , AIAA, Denver, CO, 2017 (submitted for publication).

Bencic, T., Fagan, A., Van Zante, J. F., Kirkegaard, J. P., Rohler, D. P., Maniyedath, A., and Izen, S. H., "Advanced Optical Diagnostics for Ice Crystal Cloud Measurements in the NASA Glenn Propulsion Systems Laboratory," 5th AIAA Atmospheric and Space Environments Conference , AIAA, San Diego, CA, 2013, AIAA 2013-2678.

Hobbs, P. V., Ice Physics , Oxford University Press, Oxford, 1973, Chap. 8.

Schaefer, V., “Continuous cloud chamber for studying small particles in the atmosphere,” Industrial and Engineering Chemistry , Vol. 44, No. 6, pp. 1381-3.

McDonald, J. E. , “Homogeneous Nucleation of Vapor Condensation. I. Thermodynamic Aspects,” American Journal of Physics , Vol. 30, 1962, pp. 870-7.

28 nd Wallace, J. M., and Hobbs, P. V., Atmospheric Science - An Introductory Survey , 2 ed., Academic Press – Elsevier, Burlington, 2006, Chap. 6.

Klejnowski, K., Krasa, A., Rogula-Koz ł owska, W., and Bł asz czak, B., “Number Size Distribution of Ambient particles in a Typical Urban Site: The First Polish Assessment Based on Long-Term (9 Mont hs) Measurements,” The Scientific World Journal , Vol. 2013, Article ID 539568.

Birmili, W., Heinke, K., Pitz, M., Matschullat, J., Widensohler, A., Cyrys, J., Wichmann, H., - E., and Peters, A., “Particle number size distributions in Urban air before and after volatilization,” Atmospheric Chemistry and Physics , Vol. 10, No. 10, 2010 pp. 4643-60.

Boogaard, H., Montagne, D. R., Brandenburg A. P., Meliefste K., and Hoek, G., “Comparison of short -term exposure to particle number, PM10 and soot concentrations on three (sub) urban locations ,” Science of the Total Environment , Vol. 408, No.

20, 2010, pp. 4403 – 11.

Held, A., Zerrath, A., McKeon, U., Fehrenbach, T., Niessner, R., Plass-Dülmer, C., Kaminski, U., Berresheim, H., and Pöschl, U ., “Aerosol size distributions measured in urban, rural and high- alpine air with an electrical low pressure impactor (ELPI),” Atmospheric Environment , Vol. 42, 2008, pp. 8501-12.

Stanier, C. O., Khlystov, A. Y., and Pandis, S. N., “Ambient aerosol size distributions and number concentrations measured during the Pittsburgh Air Quality Study (PAQS),” Atmospheric Environment , Vol. 38, 2004, pp. 3275-84.

Veres, J. P., Jorgenson, P. C. E, Wright, W. B, and Struk, P., “A Model to Assess the Risk of Ice Accretion due to Ice Cr ystal th Ingestion in a Turbofan Engine and its Effects on Performance ,” 4 AIAA Atmospheric and Space Environments Conference , AIAA, New Orleans, LA, 2012, AIAA 2012-3038.

American Institute of Aeronautics and Astronautics Figures Figure 1 . PSL tunnel geometry used for the fundamental ice - crystal icing physics tests with spray nozzles located at the tunnel inlet while the temperature, humidity and water content measuring instruments were located near the tunnel exit .

Figure 2 . The PSL spray bar geometry and the two nozzle patterns used during testing. The red and blue ellipses represent the approximate initial area coverage for the two different spray patterns.

American Institute of Aeronautics and Astronautics Figure 3 . The various test configurations with all probes located near the test section (tunnel exit). Images are all forward - looking - aft towards the test section . The crosshair symbols ( ⊕ ) in images B through D represent the approximate measurement locations of the HSI and PDI probes.

American Institute of Aeronautics and Astronautics Figure 4 . Sample simulation results for a test condition that exceeds saturation part way into the tunnel.

Graph A shows the 15 - μm diameter particle bin for a cloud with an inlet MVD of 15 μm, along with the particle bin that condensed on aerosol particulates. Graph B shows the static relative humidity as well as the mass mixing ratio along the axis of the PSL icing tunnel.

150 140 132.4 132.2 140 131.5 132.3 Final Relative Humidity, % Final Relative Humidity, % 1.0E+03 1.0E+04 1.0E+05 1.0E+06 0 500 1000 1500 B 3 A Initial Aerosol Diameter, nm Aerosol Number Density, #/cm Figure 5 . Parameter sweep where A) shows the static relative humidity at the tunnel exit when the aerosol number density is varied and B) shows how insensitive final static relative humidity is to the initial aerosol diameter.

American Institute of Aeronautics and Astronautics Figure 6 . View of the PSL tunnel, aft looking forward, from the spray bars for four test condition s. In images A through C, the spray nozzles are not activated, but with increasing total RH in the plenum a condensation cloud becomes visible in images B and C. The calculated static RH at the exit in images B and C are above 100%. The conditions in image D are nearly equal to those in image C, however the spray nozzles have been activated in i mage D.

Figure 7 . PSL tomography data taken near the tunnel exit for tests with exit velocities of A) U e = 85 m/s and B) U = 135 m/s.

e American Institute of Aeronautics and Astronautics Figure 8 . Comparison of experimental data with simulation results at the tunnel exit for RH Sweep #1, where U e = 85 m/s and MVD i = 15 μm. Note that the left axis in graph A is negative.

American Institute of Aeronautics and Astronautics Figure 9 . Comparison of experimental data w ith simulation results at the tunnel exit for RH Sweep #2, where U = 85 m/s and MVD = 50 μm. Note that the left axis in graph A is negative.

e i American Institute of Aeronautics and Astronautics Figure 10 . Comparison of experimental data with simulation results at the tunnel exit for RH Sweep #3, where U = 135 m/s and MVD = 15 μm. Note that the left axis in graph A is negative.

e i American Institute of Aeronautics and Astronautics Figure 11 . Comparison of experimental data with simulation results at the tunnel exit for RH Sweep #4, where U e = 135 m/s and MVD i = 50 μm. Note that the left axis in graph A is negative.

American Institute of Aeronautics and Astronautics