William B. Wright

ASRC Aerospace Corporation, Cleveland, Ohio

Philip C.E. Jorgenson and Joseph P. Veres

Glenn Research Center, Cleveland, Ohio

Mixed Phase Modeling in GlennICE With

Application to Engine Icing

NASA/TM—2011-216978

February 2011

AIAA–2010–7674
NASA STI Program . . . in Profi le

Since its founding, NASA has been dedicated to the 

advancement of aeronautics and space science. The 

NASA Scientifi c and Technical Information (STI) 

program plays a key part in helping NASA maintain 

this important role.

The NASA STI Program operates under the auspices 

of the Agency Chief Information Offi cer. It collects, 

organizes, provides for archiving, and disseminates 

NASA’s STI. The NASA STI program provides access 

to the NASA Aeronautics and Space Database and 

its public interface, the NASA Technical Reports 

Server, thus providing one of the largest collections 

of aeronautical and space science STI in the world. 

Results are published in both non-NASA channels 

and by NASA in the NASA STI Report Series, which 

includes the following report types:

• TECHNICAL PUBLICATION. Reports of 

completed research or a major signifi cant phase 

of research that present the results of NASA 

programs and include extensive data or theoretical 

analysis. Includes compilations of signifi cant 

scientifi c and technical data and information 

deemed to be of continuing reference value. 

NASA counterpart of peer-reviewed formal 

professional papers but has less stringent 

limitations on manuscript length and extent of 

graphic presentations.

• TECHNICAL MEMORANDUM. Scientifi c 

and technical fi ndings that are preliminary or 

of specialized interest, e.g., quick release 

reports, working papers, and bibliographies that 

contain minimal annotation. Does not contain 

extensive analysis.

• CONTRACTOR REPORT. Scientifi c and 

technical fi ndings by NASA-sponsored 

contractors and grantees.

• CONFERENCE PUBLICATION. Collected 

papers from scientifi c and technical 

conferences, symposia, seminars, or other 

meetings sponsored or cosponsored by NASA.

• SPECIAL PUBLICATION. Scientifi c, 

technical, or historical information from 

NASA programs, projects, and missions, often 

concerned with subjects having substantial 

public interest.

• TECHNICAL TRANSLATION. Englishlanguage translations of foreign scientifi c and 

technical material pertinent to NASA’s mission.

Specialized services also include creating custom 

thesauri, building customized databases, organizing 

and publishing research results.

For more information about the NASA STI 

program, see the following:

• Access the NASA STI program home page at 

http://www.sti.nasa.gov

• E-mail your question via the Internet to help@

sti.nasa.gov

• Fax your question to the NASA STI Help Desk 

at 443–757–5803

• Telephone the NASA STI Help Desk at

 443–757–5802

• Write to:

 NASA Center for AeroSpace Information (CASI)

 7115 Standard Drive

 Hanover, MD 21076–1320
William B. Wright

ASRC Aerospace Corporation, Cleveland, Ohio

Philip C.E. Jorgenson and Joseph P. Veres

Glenn Research Center, Cleveland, Ohio

Mixed Phase Modeling in GlennICE With

Application to Engine Icing

NASA/TM—2011-216978

February 2011

AIAA–2010–7674

National Aeronautics and

Space Administration

Glenn Research Center

Cleveland, Ohio 44135

Prepared for the

Atmospheric and Space Environments Conference

sponsored by the American Institute of Aeronautics and Astronautics

Toronto, Ontario, Canada August 2–5, 2010
Acknowledgments

The authors would like to thank the members of the NASA Glenn Icing Branch for their feedback on this report, in particular 

Jen-Ching Tsao, Colin Bidwell, Gene Addy, and Mark Potapczuk. William B. Wright would also like to thank NASA for their 

continued support under contract.

Available from

NASA Center for Aerospace Information

7115 Standard Drive

Hanover, MD 21076–1320

National Technical Information Service

5301 Shawnee Road

Alexandria, VA 22312

Available electronically at http://www.sti.nasa.gov

Trade names and trademarks are used in this report for identifi cation 

only. Their usage does not constitute an offi cial endorsement, 

either expressed or implied, by the National Aeronautics and 

Space Administration.

Level of Review: This material has been technically reviewed by technical management. 

This report is a formal draft or working 

paper, intended to solicit comments and 

ideas from a technical peer group.

This report contains preliminary fi ndings, 

subject to revision as analysis proceeds.
NASA/TM—2011-216978 1 

Mixed Phase Modeling in GlennICE With 

Application to Engine Icing 

William B. Wright 

ASRC Aerospace Corporation 

Cleveland, Ohio 44135 

Philip C.E. Jorgenson and Joseph P. Veres

National Aeronautics and Space Administration 

Glenn Research Center 

Cleveland, Ohio 44135 

Abstract

A capability for modeling ice crystals and mixed phase icing has been added to GlennICE. 

Modifications have been made to the particle trajectory algorithm and energy balance to model this 

behavior. This capability has been added as part of a larger effort to model ice crystal ingestion in aircraft 

engines. Comparisons have been made to four mixed phase ice accretions performed in the Cox icing 

tunnel in order to calibrate an ice erosion model. A sample ice ingestion case was performed using the 

Energy Efficient Engine (E3

) model in order to illustrate current capabilities. Engine performance 

characteristics were supplied using the Numerical Propulsion System Simulation (NPSS) model for this 

test case. 

Introduction 

Mason (Ref. 1) describes a situation where an aircraft engine can encounter rollbacks and flameouts 

at high altitude conditions due to ice crystal ingestion. Numerous in-fight encounters had been observed. 

It was hypothesized that the cause of the incidents was the ingestion of a high volume of ice crystals into 

the engine. More so, this report theorized that the ice crystals could be tolerated if the crystals remained 

small solid particles. However, this paper suggested that the ice crystals first melt upon contact with warm 

engine parts. Subsequent ice particles are then more likely to stick to the wetted surface, allowing ice 

accumulation to occur. The engine abnormalities then occur when this ice buildup sheds. Venkataramani 

(Ref. 2), Veillard (Ref. 3), and Rios (Ref. 4) describe previous attempts to model the engine icing 

problem. This report will describe an overview of the GlennICE and NPSS software. It will then describe 

the modifications made to GlennICE in order to model ice crystal scenarios. Included in these 

modifications is an empirical relationship for ice particle induced mass loss (ice erosion) that used to 

govern mixed-phase conditions. Data for this relationship was taken from a report by Al-Khalil (Ref. 5). 

Finally, a sample case was performed for ice crystal ingestion using the results from the NPSS model to 

provide inputs to GlennICE. 

Nomenclature 

Cd drag coefficient 

Cp heat capacity 

Da,i diffusivity of ice in air 

d particle size 

F force 

H enthalpy 

IWC ice water content 

h heat transfer coefficient 
NASA/TM—2011-216978 2 

k thermal conductivity 

L characteristic length 

Lsub heat of sublimation 

LWC liquid water content 

M Molecular weight 

m mass flux 

Nu Nusselt number 

P pressure 

Pv vapor pressure 

q heat flux 

r coefficient of restitution 

Re Reynolds number 

Sc Schmidt number 

Sh Sherwood number 

Stm mass Stanton number 

Stmom momentum Stanton number 

T temperature 

t time 

V velocity 

Vc critical velocity 

 collection efficiency 

 strain rate 

 mass transfer parameter 

 viscosity 

 freezing fraction

 density 

 shear stress 

m dimensionless mass parameter 

mom dimensionless momentum parameter 

t dimensionless thermal parameter 

 mass fraction 

Subscripts: 

a air 

conv convection 

ev evaporation 

f freezing 

film film 

heat heat added 

i ice 

KE kinetic energy 

m melting 

o original 

rb runback 

s surface 

sens sensible 

sub sublimation 

t threshold 

w water 

∞ ambient 
NASA/TM—2011-216978 3 

GlennICE Overview 

The GlennICE (Ref. 6) project was initiated in 1999 as a means to develop a consistent set of 

standards and practices for creation, use and maintenance of NASA administered icing computational 

tools. One goal of this process is to provide a mechanism for a more team-based approach to software 

development. Current development of icing codes at NASA has been more separated and individualbased. The plan has been to establish a common computational framework for integration of NASA icing 

analytical tools, and to create a set of requirements both for the development of these tools and for their 

integration into the common framework. The features of the resulting system would include providing a 

common user interface for all tools and allowing user selection of analysis modules. A modular design 

would allow more streamlined substitution of new or improved software elements, and the defined 

process would include documented standards for rapid integration of new software. This process also 

requires a version control system designed to allow both control of approved software and flexibility for 

rapid prototyping, as well as documentation that allows rapid identification of software design and use. 

NPSS Overview 

The Numerical Propulsion System Simulation (NPSS) code (Refs. 7 to 9) is used to model the 

thermodynamic cycle of complete aircraft engines and therefore can be effective at predicting the engine 

system operation. The development of this tool is a joint effort between NASA Glenn Research Center, 

the Department of Defense, and Unites States aeropropulsion industry. It is used to model the physical 

interactions between components of a turbine engine reducing the need for full-scale engine tests and 

component experiments ultimately impacting the time it takes to bring an aircraft engine to market. 

At its core the NPSS is an object-oriented code that brings together multiple disciplines such as 

aerodynamics, structures, and heat transfer which supports multi-fidelity (numerical zooming) from zerodimensional to three-dimensional engine component codes. Zooming allows a designer to vary the level 

of detail of analysis throughout the engine based upon the physical processes being studied. It is 

important to consider the aircraft engine as a system of components interacting with each other not simply 

in isolation. In the NPSS these component models, potentially at various levels of fidelity, are allowed to 

interact with each other to see the effects they have on the overall engine performance. A simulation of a 

complete aircraft engine will provide sufficient detail of the individual components and their interactions 

with other engine components over a range of flight operating conditions. This allows for an examination 

of many design options before costly tests are conducted. 

GlennICE Modifications 

Three modifications have been made to GlennICE for handling ice particle impacts. First, there are 

the modifications to the particle trajectory routine that track the energy transfer to and from a particle. 

Second, there is a model for mass loss due to erosion in a mixed phase environment. Finally, there are 

changes to the mass and energy balance whereby the incoming ice particles can add to the ice mass at the 

surface. In addition to these modifications, an analysis was performed to assess if coupling between the 

air and ice phases is necessary. 

There are two separate modifications to the particle trajectory equations in order to handle ice 

crystals. First, a separate drag equation is needed. The drag equation for liquid droplets includes an 

increase in drag due to particle deformation whereas ice crystals remain spherical. It is also assumes that 

ice crystals do not breakup due to shear. In addition, drag correlations were added for ice particles that are 

cylindrical as well as disk-shaped. The second modification was to track the change in particle 

temperature (and phase change if necessary) through the air stream. Previously, GlennICE did not track 

droplet temperature as they were assumed to be in equilibrium with the air. However, the mixed phase 

capabilities have been added to analyze engine ice events. In this scenario, the air, water drops and ice 

crystals will increase in temperature when they enter the engine. Since the air will respond more rapidly 
NASA/TM—2011-216978 4 

to this increase, there will exist a driving force for energy transfer between the particles and the 

surrounding air stream. 

Particle Drag Equations 

The following equations for particle drag were obtained by curve-fitting drag data for spheres, 

cylinders, and disks. The reference used was Perry’s Chemical Engineers’ Handbook (Ref. 10), but 

several references contain the same data. 

For spheres 

for Re <100

1 Re

6 0 4

Re

24



Cd   .  (1)

for Re >100

1 Re

6 0 3

Re

24



Cd   .  (2)

For cylinders 

4 194Re for Re < 0.01 0.931 .  Cd  (3)

13 344Re for 0.01< Re < 20 0.691 .  Cd  (4)

Cd  3.2344  0.334logRe for 20 < Re <1000 (5)

Cd  0.93 for Re 1000 (6)

For disks 

for Re < 40

1 Re

6 0 4

Re

24



Cd   .  (7)

1 5 10 Re 0 001Re 1 8176 for 40 < Re <1000 6 2  . *  .  .  Cd (8)

Cd 1.3 forRe 1000 (9)

Particle Energy Balance 

The following analysis will present the energy equation for a spherical ice particle. The particle 

energy equation is simply the increase or decrease in enthalpy of the particle balanced with the losses due 

to convection and sublimation or 

qconv qsub/evap t

H  



 (10)

Prior to or after any phase change by melting of the ice crystal, this equation can be given as 
NASA/TM—2011-216978 5 

   a s a m a s s i p i h d T T d h L

t

d T C         



 

 sub

2 2

3

6 , (11) 

for spheres 

After a particle has started to change phase, Lsub is replaced by Levap in the equation above. It will be 

more convenient to solve this equation in non-dimensional form, hence 

a k

hd Nu  (12) 

a i

m

D

h d

,

Sh  (13)

a

p a a

k

C   , Pr (14)

a a i

a

 D ,

 Sc  (15)

P

P

M

M v

a

w   (16)

a

p i i

T k

C d

12

2 

  , (17)

The resulting energy equation then becomes 

    s

T p i

a s

T

s

C

L T T

dt

dT

  

    

,

sub

Sc

Pr

2

Sh

2

Nu (18) 

where 

For spherical particles 

1 2 1 3 Nu  2  0.6Re Pr (19)

1 2 1 3 Sh  2  0.6Re Sc (20)

where 

a

a Va Vi d



  Re  (21)

For cylindrical particles 
NASA/TM—2011-216978 6 

1 4 2 3

0 8 0 625

5

1 2 1 3

Pr

0 4 1

2 82 10

Re 1

Nu 0 3 0 62Re Pr





















 



 





















 







 

.

. * . .

. .

(22)

1 4 2 3

0 8 0 625

5

1 2 1 3

Sc

0 4 1

2 82 10

Re 1

Sh 0 3 0 62Re Sc





















 



 





















 



 

 

.

. * . .

. .

(23)

For discs 

Nu 0 664Re Pr (Re < 200,000) 1 2 1 3  . (24) 

  (Re > 200,000)

1 2 443Re Pr 1

0 037Re Pr Nu 0 1 2 3

0 8

    .

.

.

. (25) 

Sh 0 664Re Sc (Re < 200,000) 1 2 1 3  . (26) 

  (Re > 200,000)

1 2 443Re Sc 1

0 037Re Sc Sh 0 1 2 3

0 8

    .

.

.

. (27) 

Reduction in particle size for a sphere is given by the mass balance 

d

d D d

t

a s i a a i

     













 





 , 2

3

Sh

6 (28) 

or after integration 

d d t  o   2 2 (29)

     

   s

i

4Sh aDa,i (30)

During phase change, the ice/water temperature remains fixed and the rate of melting of ice is given 

by 

i f h d    Ta Tm d hmL a m t

d L        



   evap

2 2

3

6 (31) 
NASA/TM—2011-216978 7 

By fixing the crystal temperature at the melting point, the melting rate can be calculated directly from 

this equation. 

Coupling Between Particles and Air 

In previous icing analysis, the airflow and particle dynamics were calculated separately, which 

implies that coupling of the equations are unnecessary. In this regime, it should be clarified whether or 

not coupling is needed. Crowe (Ref. 11) estimated mass coupling to be important if 

1

St



m

Z (32)

where 

a

Z

  IWC (33)

L

mVi m

 St  (34)

    

  

a a i s

i m D

d

Sh , 6

2

(35)

If we try to maximize the coupling, we need the maximum ice water content, the smallest particle 

size, the lowest ice velocity, the longest compressor length and the highest mass transfer rate. Let IWC = 

9 g/m3

, d = 15 μm, Vice = 40 m/s, P = 27,000 Pa, Sh = 5, L = 1.4 m and (s–) = 0.02 then a = 0.4, m = 

0.044 and 0 018

St  .

m

Z , thus mass coupling is unimportant. 

Momentum coupling is important when 

1

1 St



 mom

Z (36) 

where 

L

vVi mom

 St  (37)

a

d v

d



   18

2

(38)

thus using the previous values, 

v  0.000744 (39)

Stmom  0.021 (40)

Z  0.0225 0 022

1 0 021

0 0225 . .

.  

(41)
NASA/TM—2011-216978 8 

so momentum coupling can be ignored. 

Energy coupling is given by 

0 1

St

sub .

,

 C T

Z L

m p i

(42)

using the previous values so energy coupling is unimportant as well. Based on this analysis, the one-way 

coupling performed in icing software such as GlennIce is still justified. 

Coefficient of Restitution 

When a solid particle hits a surface, the resulting collision is considered elastic or inelastic. If the 

collision is elastic, none of the incoming energy is absorbed by the surface. For inelastic collisions, some 

of the energy is absorbed by the surface and some is used to determine the velocity of the particle after the 

collision. This energy ratio is called the coefficient of restitution. If this coefficient is equal to one, the 

collision is considered elastic. If this coefficient is equal to zero, all of the energy is absorbed by the 

surface. In this case, all of the mass must remain on the surface as well. Therefore, in the ice crystal 

environment being modeled, the coefficient of restitution can be used to determine the fraction of 

impinging ice crystals that remain on the surface, similar to the splashing model used for supercooled 

large droplets. Higa et al. (Ref. 12) measured the coefficient of restitution for an ice sphere with a 3 cm 

diameter at velocities up to 0.7 m/s and a wide range of temperatures. They curve-fit the restitution 

coefficient to the following equation: 















 















  Vc

V

Vc

V

r

log

(43)

They reported values of the critical velocity (Vc) after which particle collisions were no longer quasielastic at velocities of 25, 35, and 61 cm/s at 269, 261, and 245 K, respectively. The decrease in 

restitution coefficient was caused by fragmentation of the ice spheres. Based on this equation, most ice 

particle collisions at the surface will have essentially a zero restitution coefficient, meaning that almost all 

of the energy of the collision is transferred to the surface and the ice crystals are expected to fragment. 

Initial build up of an ice layer is then caused in part by the kinetic energy of the collision being used to 

melt the surface ice fragments. GlennICE assumes that the fraction of ice crystals present on the surface 

will be (1-r) as this is the fraction of energy transferred to the surface during the collision with the 

surface. However, this does not necessarily mean that the ice crystals will accumulate. In order for 

accretion to occur, a film thickness must be present. Otherwise, the ice fragments would be swept away 

by the flow field. In the limited mixed phase experiments presented below, the ice crystals do not 

accumulate even in a glaze ice condition where a water film is present. However in this case the water 

film is limited to a very small region near the stagnation zone. In the current model, the energy and mass 

balances are calculated assuming the ice crystals remain on the surface. For control volumes that do not 

have a liquid film, the ice crystal mass is removed prior to the geometry modification routines. As more 

experimental data becomes available, this assumption can be revisited. Erosion effects, as described 

below, are accounted for in all regimes. 

Ice Erosion 

Ice erosion occurs when ice particles hit an ice shape that has already accreted. The force of impact 

can dislodge portions of the ice shape, especially individual ice feathers. Al-Khalil (Ref. 5) performed 

tests in the Cox icing tunnel in order to document the mixed phase capabilities of that facility. The model 
NASA/TM—2011-216978 9 

was a 36 in. chord NACA0012 airfoil. All cases were run for 10 min. at 100 kts, = 0° and a water MVD 

= 20 μm. Results from Al-Khalil’s report show that ice impacts can cause erosion of an ice shape in a 

mixed phase test. Four cases were selected for calibration with the GlennICE erosion model. The cases 

are given in Table 1. 

The cases at the lower temperature (cases 3 and 4) were unaffected by the ice crystals as the ice 

thickness compares well with values from GlennICE without any erosion effects. This suggests that the 

low air speed of these conditions is not high enough to cause erosion. At the warmer condition (cases 14 

and 15), the introduction of an IWC of 0.7 g/m3

 caused a 20 percent reduction in the ice shape as shown 

in Figure 1. The amount of ice loss was calculated using the utility program THICK. This result seems to 

indicate that the lower bonding strength of the ice at the warmer condition is enough to cause the onset of 

ice erosion. Rist (Ref. 13) measured the deformation and fracture of 1.5 cm ice particles over a wide 

range of temperatures. He measured the strain rate at failure to the corresponding stress applied and 

produced the following relationship: 

  



 





      RT

Q exp 14 4 2

1 3 .  (44) 

In this equation, Q is the activation energy that was measured as 69 kJ/mol in the experiment and R is 

the gas constant. For the cases listed, the warmer condition would require only half as much force to 

produce the same strain rate in the ice as the colder temperature cases. The impacting force of each 

particle is proportional to the mass of the particle (and thus proportional to the diameter cubed), the 

velocity and the particle frequency. Without performing a stress analysis of each ice shape there is 

insufficient data to build a proper correlation due to the sparseness of the test matrix. However if it is 

assumed that the cases without erosion effect are at the erosion threshold, the erosion effect can be 

estimated. First, the impact stress was assumed proportional to the force of the ice particle impacts. This 

force was determined by multiplying the particle momentum with the frequency: 





































 



 







  3

1

3 IWC

6 i

i d

V F d V (45) 

The term in the second bracket is the particle frequency. This equation would make the stress 

proportional to d2

, V2, and IWC1/3. This would make the applied stress the same for the cases where IWC 

= 0.7 g/m3

. However, since the temperature is higher in case 15 than case 3, twice as much strain rate is 

produced in the ice, causing a 20 percent reduction in ice size. The model in GlennICE first assumes that 

case 4 from the Al-Khalil data is at the threshold where erosion becomes important. For this condition and 

for conditions with lower shear stresses due to particle impact, no erosion takes place. For conditions with 

higher impact energy, the ratio of strain rate can be estimated by 











































 





 



 































 



 



  



t

t t

t

t

RT

Q

RT

Q

V

V

d

d

exp 14

exp 14

IWC

IWC 2 2 1 3



 (46) 

The erosion threshold values are dt = 150 μm, Vt = 100 kts, IWCt = 0.7 g/m3, and Tt = 267 K which is 

the average ice temperature in the impingement region for this condition. If it is further assumed that the 

erosion rate is proportional to the increase in strain rate, the erosion rate can then be related back to the 

incoming conditions. Comparing case 4 with case 15 from Al-Khalil’s data gives: 
NASA/TM—2011-216978 10 

3 5

8 314 10 267

69

exp 14

8 314 10 273

69

exp 14

0 7

0 7

100

100

150

200

3

3 2 2 1 3

.

. * *

. * *

.

. 

































 



 





 



 





 



 



 



 



 



  







t 

 (47) 

Given that case 4 does not appear to have any erosion effects while case 15 shows a 20 percent loss in 

mass due to erosion leads to the following empirical equation: 















  

 erosion  0 08 1

t 

 . (48) 

where t  is the strain rate at the erosion threshold. Obviously, much more data would be necessary to 

determine a more accurate value for the erosion threshold as well as to verify the dependence of erosion 

on the icing parameters. 

Surface Energy Balance 

The mass and energy balance add terms for the added mass of ice, the kinetic impact of the ice 

particles, and the melting of ice crystals. For a glaze surface, 

mimp,w mimp,i m f mrb mev msub mfilm              (49)

where mimp,i  1 erosion 1 riIWCV (50)

r is the restitution coefficient calculated earlier. 

The energy equation becomes 

0 qf/m  qheat  qKE,a  qKE,w  qKE,i  qconv  qevap  qsub  qsens  (51)

where qf/m is the heat removed by additional freezing or heat added by melting of ice. As described below, 

GlennICE assumes that only one of these processes can take place in a particular control volume. 

The first term that is added due to ice crystal impacts is the kinetic energy transfer of the impacting 

crystals, 

q  m V   1 r

2

1 2

KE,i imp,i imp  (52) 

where r is the restitution coefficient calculated earlier. The second term added is the melting term given 

above. The final term that is changed is the sensible heat transfer from the drop or particle. Previously, the 

drops were assumed to be at the ambient air temperature. Currently, the particle temperature is calculated 

as shown earlier and the resulting temperature at impact is transferred to the icing balance. In addition, 

any phase change or changes in particle and droplet size are accounted for in these equations. 

In the solution procedure, the terms other than melting and freezing are summed. Based on the sign of 

that term, the remaining energy is used to either freeze additional surface water or used to melt incoming 

ice. In a situation where more energy is carried away allowing additional freezing, then the freezing 

fraction of the incoming ice crystals is equal to one and the additional freezing fraction of the incoming 

water is calculated from the energy balance. In a situation where more energy is available to melt the 

incoming ice, then the freezing fraction of the incoming water is zero and the fraction of incoming ice that 
NASA/TM—2011-216978 11 

is melted is again calculated from the energy equation. In any situation where the sum of the other terms 

is exactly zero, then the freeing fraction of the incoming water is zero and the freezing fraction of the 

incoming ice is equal to one. The solution procedure determines the mass of ice or water directly and does 

not rely on calculating a freezing fraction. 

Results 

The energy efficient engine, or E3, was designed at NASA Glenn for the intention to lay an advanced 

technology foundation for a new generation of turbofan engines. It was a cooperative program with 

industry aimed at developing and demonstrating advanced component and systems technologies for 

engines. Since the geometry is public domain, it was selected for initial tests of the computer model. 

Figure 2 shows a top view of this engine while Figure 3 shows the cross-section. 

The NPSS code was used to fly the E3

 engine over a typical flight envelope to determine locations in 

the engine where ice accretion might occur. The flight envelope provided the NPSS code with the altitude 

and throttle settings for the engine performance analysis. The NPSS code was used in a 0-D 

(thermodynamic cycle analysis) mode to determine locations in the compression system where ice 

accretion might be supported. Figure 4 illustrates the total temperature variation over a portion of the 

flight envelope for the E3

 at engine station 2.2. Based on icing event statistics of commuter and large 

transport aircraft reported by Mason (Ref. 1), the E3

 was flown over the flight envelope at various altitude 

and temperatures (at and above standard atmospheric conditions). The representative flight envelope 

shown in Figure 4 includes two cruise altitudes and one descent profile. During the two cruise portions of 

the flight, the engine operates at reduced thrust after attaining the desired altitude. At the end of the 

second cruising altitude, the engine is spooled back to 5 percent maximum thrust in preparation for 

descent. The aircraft then descends from maximum cruising altitude to prepare for landing at a reduced 

thrust. An engine operating point at 22 R above standard atmospheric conditions within the 39000 ft 

cruise portion of the flight envelope was selected for testing GlennICE. Figure 5 shows the fan tip 

performance map of the engine at the elevated temperature. The fan operates well away from the surge 

line over the flight profile. However as indicated in Figure 6, the fan hub and quarter stage compression 

system operates near the surge line and with geometry changes due to ice accretion has the potential to 

push the compression system toward surge thus compromising the engine performance and causing a 

safety concern. 

The initial analysis performed for this report focuses on the region from the inlet through the lowpressure compressor (LPC) to the entrance of the high-pressure compressor (HPC). The system modeled 

in NPSS using a 0-D analysis was used to obtain temperature, pressure, and flow rates at various engine 

stations. A diagram of the region analyzed is shown in Figure 7 and the results of the NPSS analysis at the 

selected location in the flight envelope are shown in Table 2. 

A zero-order analysis only provides values at specific locations such as the conditions before and 

after the fan for example. It does not include the dimensions of the system, hence the use of the term zeroorder. While NPSS is capable of performing a higher fidelity model, the zero-order results were sufficient 

for the purposed of testing the initial model. In order to use this information in GlennICE, a grid was 

created in GRIDGEN and the values from Table 2 were used to supply velocities, pressures, and 

temperatures at the grid points. Since a flow analysis was not performed using this grid, it can be 

somewhat coarser than that used for Navier-Stokes analysis since its only purpose is for tracking particle 

trajectories. There were 10 grid blocks used with a total of 19,000 grid points in order to interpolate the 

NPSS results to a grid that GlennICE could use. Spherical ice crystals of 200 μm MVD were then tracked 

through this artificial flow field to find collection efficiencies on the two splitters. This is shown in 

Figure 8. 

The secondary splitter had minimal impingement and no ice growth, so results for that geometry are 

not shown. An icing time of 5 min and an IWC of 8 g/m3

 were used to generate the results for this case. 

Velocity, temperature and pressure were taken as the conditions at the engine inlet as shown in Table 1. 

Results from the energy balance and mass balance are shown in Figures 9 and 10, respectively. Figure 11 
NASA/TM—2011-216978 12 

shows what the ice shape would look like without erosion effects. Based on the preliminary ice erosion 

model used, this ice shape would be completely eroded by the impact of other ice crystals due to the much 

higher velocities inside the engine. It should be emphasized that this result has not been experimentally 

verified. 

Future Work 

This model must be considered at the early stages of development at this juncture. Much work 

remains in order to produce a model ready for release. In particular, experimental data is needed to verify 

erosion effects and the deposition rate of ice crystals when a water film is present. Just as early research in 

the supercooled droplet regime provided references of droplet splashing at low velocities that did not 

coincide with results obtained in the NASA icing tunnel at higher velocities, the models presented in this 

report for coefficient of restitution, particle breakup and ice erosion at low velocities need to be confirmed 

at the higher velocities inside the engine. In addition, while particles impinge on the lower wall of the 

input geometry (shown in Fig. 7), GlennICE needs to be modified to allow deposition and accretion on 

this surface. Experimental data on actual engine geometries in an ice crystal environment is also 

necessary for validation. It would also be useful to test the model with two-dimensional flow solutions 

from NPSS. Once a viable model is ready, it will also be implemented into LEWICE3D so that the full 

geometry can be analyzed. 

Conclusions 

A preliminary model for multiphase ice accretion was developed for GlennICE. An analysis of typical 

ice crystal inputs showed that one-way coupling could still be used in the icing model. The effects of 

temperature change, phase change and sublimation of the ice crystal were added to the trajectory model. 

In addition, non-spherical ice crystal capability was also added. A factor for the coefficient of restitution 

was also added. The energy equation was revised to include the effects of ice crystals, including a 

preliminary model for ice erosion. Erosion results were calibrated to a mixed phase icing test performed 

by Al-Khalil in the Cox icing facility. Icing results were produced for a test case using the E3

 geometry. 

The results show that while ice could form on the splitter lip based on the energy balance, it would be 

quickly eroded by other ice particle impacts. However, this result has not been experimentally validated. 

References 

1. Mason, J.G., Strapp, J.W., and Chow, P., “The Ice Particle Threat to Engines in Flight,” 

AIAA–2006–206. 

2. Venkataramani, K.S., Plybon, R.C., Holm, R.G., and Krupp, B.R., “Aircraft Engine Icing Model,” 

AIAA–2008–440. 

3. Veillard, X., Habashi, W.G., Aube, M.S., and Baruzzi, G.S., “FENSAP-ICE: Ice Accretion in Multistage Jet Engines,” AIAA–2009–4158. 

4. Rios, M.A., and Cho, Y.I., “Analysis of Ice Crystal Ingestion as a Source of Ice Accretion inside 

Turbofans,” AIAA–2008–4365. 

5. Al-Khalil, K., “Assessment of Effects of Mixed Phase Icing Conditions on Thermal Protection 

Systems,” DOT/FAA/AR-03/48, May, 2003. 

6. Wright, W.B., Potapczuk, M.G., and Levinson, L.H., “Comparison of LEWICE and GlennICE in the 

SLD Regime,” AIAA–2008–0439. 

7. Lytle, J.K., “The Numerical Propulsion System Simulation: An Overview,” NASA/TM—2000-

209915. 

8. Jones, Scott M., “An Introduction to Thermodynamic Performance Analysis of Aircraft Gas Turbine 

Engine Cycles Using the Numerical Propulsion System Simulation Code,” NASA/TM—2007-

214690. 
NASA/TM—2011-216978 13 

9. Turner, Mark G.; Reed, John A.; Ryder, Robert; Veres, Joseph P., “Multi-Fidelity Simulation of a 

Turbofan Engine With Results Zoomed Into Mini-Maps for a Zero-D Cycle Simulation,” 

NASA/TM—2004-213076. 

10. Green, D.W., Perry’s Chemical Engineer’s Handbook, 6th Edition, Sect. 5, pg. 54. 

11. Higa, M., Arakawa, M., Maeno, N., “Measuremens of Restitution Coefficients of Ice at Low 

Temperatures,” Planet. Space Sci, vol. 44, no. 9, pp. 917–925, 1996. 

12. Crowe, C., Sommerfeld, M. and Tsuji, Y., Multiphase Flows with Droplets and Particles, CRC Press, 

1998. 

13. Rist, M.A., and Murrell, S.A.F., “Ice Triaxial Deformation and Fracture, J. of Glaciology,” v. 40, 

no. 135, pp. 305–18, 1994. 

TABLE 1.—MIXED PHASE ICING RUNS 

Case no. To, 

°F 

MVD, ice,

μm 

IWC, 

g/m3

LWC, 

g/m3

3 12 150 0.3 0.7 

4 12 150 0.7 0.3 

14 22 N/A 0 0.7 

15 22 200 0.7 0.7 

TABLE 2.—FLOW CONDITIONS FROM NPSS ZERO-ORDER ANALYSIS 

Station T, 

R 

P, 

psia 

Mach V, 

mph 

W, 

lbm/s 

1 462.2 4.27 0.78 776.4 1313.3 

2 462.2 4.27 0.57 580.2 1313.3 

2T 462.2 4.27 0.57 580.2 1017.5 

2H 462.2 4.27 0.57 580.2 295.8 

21 513.3 5.94 0.48 520.0 224.0 

22 530.7 6.67 0.48 534.6 202.8 

23 530.7 6.67 0.48 534.6 96.6 

24 530.7 6.67 0.21 240.0 106.2 

25 530.7 6.67 0.40 439.8 106.2 

13 527.4 6.47 0.46 502.5 716.3 

14 527.4 6.47 0.48 526.1 716.3 

15 527.4 6.50 0.48 531.4 813.0 
NASA/TM—2011-216978 14 

Figure 1.—Erosion effect with mixed phase inputs. 

Figure 2.—Top view of energy efficient engine. 
NASA/TM—2011-216978 15 

Figure 3.—Cross-section view of energy efficient engine. 

Figure 4.—Total temperature range at engine station 2.2 of the E3 over a prescribed 

flight envelope at 22 R above standard atmospheric conditions. Engine operating 

point in cruise at 39000 ft was selected for testing GlennICE. 

Altitude, ft

Tt, R @ Station 2.2

10000 15000 20000 25000 30000 35000 40000

480

500

520

540

560

Cruise 1

Cruise 2

Descent

dT = +22 R

100% max thrust

5% max thrust

Location on Flight Envelope

of Icing Analysis
NASA/TM—2011-216978 16 

Figure 5.—Fan tip performance map showing corrected flow rate versus total pressure ratio, 

for a range of corrected speed lines. Key engine operating points along the vehicle trajectory 

are superimposed onto the fan map. Corrected speed lines are expressed as a percent of 

the design rotational speed. 

Figure 6.—Fan hub and quarter stage performance map showing corrected flow rate versus 

total pressure ratio, for a range of corrected speed lines. Key engine operating points along 

the vehicle trajectory are superimposed onto the fan map. Corrected speed lines are 

expressed as a percent of the design rotational speed. 

Flow Rate, corrected, lbm/sec

Pressure Ratio

500 1000 1500 1

1.1

1.2

1.3

1.4

1.5

1.6

1.7

1.8

1.9

2 Design 35K .80-M 100%-RPM 100%-Thrust

Design 39K .78-M 100%-RPM 100%-Thrust

Cruise 2 39K .78-M 93%-RPM 80%-Thrust

Cruise 1 35K .78-M 92%-RPM 75%-Thrust

Descent 39K .78-M 65%-RPM 5%-Thrust

Descent 25K .55-M 56%-RPM 5%-Thrust

Descent 20K .49-M 52%-RPM 5%-Thrust

Descent 15K .45-M 50%-RPM 5%-Thrust

Descent 10K .40-M 48%-RPM 5%-Thrust

105%

RPM

100%

RPM

9

90

80

%

7

%

0

60

%

%

50%

Surge Line

85%

95%

RPM

40%

10K

1

25K

20K

5K

39K

30%

Flow Rate, corrected, lbm/sec

Pressure Ratio

0 50 100 150 200 250 300 350 1

1.1

1.2

1.3

1.4

1.5

1.6

1.7

1.8

1.9

Design 35K .80-M 100%-RPM 100%-Thrust

Design 39K .78-M 100%-RPM 100%-Thrust

Cruise 2 39K .78-M 93%-RPM 80%-Thrust

Cruise 1 35K .78-M 92%-RPM 75%-Thrust

Descent 39K .78-M 65%-RPM 5%-Thrust

Descent 25K .55-M 56%-RPM 5%-Thrust

Descent 20K .49-M 52%-RPM 5%-Thrust

Descent 15K .45-M 50%-RPM 5%-Thrust

Descent 10K .40-M 48%-RPM 5%-Thrust

105% RPM

100% RPM

95% RPM

90% RPM

85% RPM

80% RPM

70% RPM

60% RPM

5

40

0%

30%

%

39K

25K

20K

1

10

5K

K

Surge Line
NASA/TM—2011-216978 17 

Figure 7.—Engine region analyzed with NPSS and GlennICE. 

Figure 8.—Collection efficiency for E3 test case. 
NASA/TM—2011-216978 18 

Figure 9.—Energy balance terms for E3 test case. 

Figure 10.—Mass balance terms for E3 test case. 
NASA/TM—2011-216978 19 

Figure 11.—Ice shape for E3 test case. 
REPORT DOCUMENTATION PAGE Form Approved 

OMB No. 0704-0188 

The public reporting burden for this collection of information is estimated to average 1 hour per response, including the time for reviewing instructions, searching existing data sources, gathering and maintaining the 

data needed, and completing and reviewing the collection of information. Send comments regarding this burden estimate or any other aspect of this collection of information, including suggestions for reducing this 

burden, to Department of Defense, Washington Headquarters Services, Directorate for Information Operations and Reports (0704-0188), 1215 Jefferson Davis Highway, Suite 1204, Arlington, VA 22202-4302. 

Respondents should be aware that notwithstanding any other provision of law, no person shall be subject to any penalty for failing to comply with a collection of information if it does not display a currently valid OMB 

control number. 

PLEASE DO NOT RETURN YOUR FORM TO THE ABOVE ADDRESS. 

1. REPORT DATE (DD-MM-YYYY)

01-02-2011 

2. REPORT TYPE

Technical Memorandum 

3. DATES COVERED (From - To)

4. TITLE AND SUBTITLE

Mixed Phase Modeling in GlennICE With Application to Engine Icing 

5a. CONTRACT NUMBER

5b. GRANT NUMBER

5c. PROGRAM ELEMENT NUMBER

6. AUTHOR(S)

Wright, William, B.; Jorgenson, Philip, C.E.; Veres, Joseph, P. 

5d. PROJECT NUMBER

5e. TASK NUMBER

5f. WORK UNIT NUMBER

WBS 609866.02.07.03.04 

7. PERFORMING ORGANIZATION NAME(S) AND ADDRESS(ES)

National Aeronautics and Space Administration 

John H. Glenn Research Center at Lewis Field 

Cleveland, Ohio 44135-3191 

8. PERFORMING ORGANIZATION

 REPORT NUMBER

E-17618 

9. SPONSORING/MONITORING AGENCY NAME(S) AND ADDRESS(ES)

National Aeronautics and Space Administration 

Washington, DC 20546-0001 

10. SPONSORING/MONITOR'S

 ACRONYM(S)

NASA 

11. SPONSORING/MONITORING

 REPORT NUMBER

NASA/TM-2011-216978 

12. DISTRIBUTION/AVAILABILITY STATEMENT

Unclassified-Unlimited 

Subject Category: 03 

Available electronically at http://www.sti.nasa.gov 

This publication is available from the NASA Center for AeroSpace Information, 443-757-5802

13. SUPPLEMENTARY NOTES

14. ABSTRACT

A capability for modeling ice crystals and mixed phase icing has been added to GlennICE. Modifications have been made to the particle 

trajectory algorithm and energy balance to model this behavior. This capability has been added as part of a larger effort to model ice crystal 

ingestion in aircraft engines. Comparisons have been made to four mixed phase ice accretions performed in the Cox icing tunnel in order to 

calibrate an ice erosion model. A sample ice ingestion case was performed using the Energy Efficient Engine (E3) model in order to 

illustrate current capabilities. Engine performance characteristics were supplied using the Numerical Propulsion System Simulation (NPSS) 

model for this test case. 

15. SUBJECT TERMS

Turbojet engines; Aircraft icing 

16. SECURITY CLASSIFICATION OF: 17. LIMITATION OF

 ABSTRACT

UU 

18. NUMBER

 OF 

 PAGES

25 

19a. NAME OF RESPONSIBLE PERSON

STI Help Desk (email:help@sti.nasa.gov) 

a. REPORT

U 

b. ABSTRACT

U 

c. THIS 

PAGE

U 

19b. TELEPHONE NUMBER (include area code)

443-757-5802 

Standard Form 298 (Rev. 8-98)

Prescribed by ANSI Std. Z39-18