a
plinth;
. ..

14.. I
. ‘ .:
unanswuu
a...»

 

 

 

m“ a.

.t
¢ 7.». v.
:1... 5311

19.1}. ‘3

x, .
I.‘

”fig-nuvauflm-

an“ :11

'9‘

.
v
I

5].! W.
. t?- I
x 11-...“

s...
are: ‘ ,

#fiigg 3.1.; .33 ..

 

 

 

If
I rrrt>|l

 

to a

006

Michig'a‘nns'tate
University

This is to certify that the
dissertation entitled

NUMERICAL SIMULATION OF
ROUGH-SURFACE AERODYNAMICS

presented by
Xingkai Chi

has been accepted towards fulfillment
of the requirements for the

PhD. degree in Mechanical lingineerim

 

 

\

 

Major Pbbfessor’s Signature

lo: zoo:

 

Date

MSU is an Affinnative Action/Equal Opportunity Institution

PLACE IN RETURN BOX to remove this checkout from your record.
TO AVOID FINES return on or before date due.
MAY BE RECALLED with earlier due date if requested.

 

DATE DUE

DATE DUE

DATE DUE

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

2/05 clam.mws

 

NUMERICAL SIMULATION OF
ROUGH-SURFACE AERODYNAMICS

By

Xingkai Chi

A DISSERTATION

Submitted to
Michigan State University
in partial fulfillment of the requirements
for the degree of

DOCTOR OF PHILOSOPHY

Department of Mechanical Engineering

2005

ABSTRACT

NUMERICAL SIMULATION OF ROUGH-SURFACE AERODYNAMICS
By
Xingkai Chi

Computational fluid dynamics (CFD) simulations of flow over surfaces with
roughness in which the details of the surface geometry must be resolved pose major
challenges. The objective of this study is to address these challenges through two
important engineering problems, where roughness play a critical role - flow over airfoils
with accrued ice and flow and heat transfer over turbine blade surfaces roughened by
erosion and/or deposition.

CFD simulations of iced airfoils face two major challenges. The first is how to
generate high-quality single- and multi-block structured grids for highly convoluted
convex and concave surface geometries with multiple scales. In this study, two methods
were developed for the generation of high-quality grids for such geometries. The method
developed for single-block gn'ds involves generating a grid about the clean airfoil,
carving out a portion of that grid about the airfoil, replacing that portion with a grid that
accounts for the accrued ice geometry, and performing elliptic smoothing. The method
developed for multi-block grids involves a transition-layer grid to ensure jaggedness in
the ice geometry does not propagate into the domain. It also involves a “thick” wrap-
around grid about the ice to ensure grid lines clustered next to solid surfaces do not
propagate as streaks of tightly packed grid lines into the domain along block boundaries.
For multi-block grids, this study also developed blocking topologies that ensure solutions

to multi-block grids converge to steady state as quickly as single-block grids.

The second major challenge in CFD simulations of iced airfoils is not knowing
when it will predict reliably because of uncertainties in the turbulence modeling. In this
study, the effects of turbulence models in predicting lift, drag, and moment coefficients
were examined for airfoils with rime ice (i.e., ice with jaggedness only) and with glaze
ice (i.e., ice with multiple protruding horns and surface jaggedness) as a function of angle
of attack. In this examination, three different CFD codes — WIND, FLUENT, and
PowerFLOW were used to examine a variety of turbulence models, including Spalart-
Allmaras, RNG k-s, shear-stress transport, v2-f, and differential Reynolds stress with and
without non—equilibrium wall functions. The accuracy of the CFD predictions was
evaluated by comparing grid-independent solutions with measured experimental data.
Results obtained show CFD with WIND and FLUENT to predict the aerodynamics of
airfoils with rime ice reliably up to near stall for all turbulence models investigated.
PowerFLOW was able to predict lift well even up to stall because unsteadiness, which
occurs near stall can be accurately simulated. For airfoil with glaze ice, CF D predictions
by all codes and turbulence models were much less satisfactory because the horns induce
considerable separated regions about the airfoil’s leading edge.

The iced airfoils simulations were 2-D. To examine modeling challenges posed
by 3-D roughness, CFD simulations were performed of the flow and heat transfer over a
3-D roughness panel, extracted from a “roughened” turbine blade surface by using a
Cartesian mesh with prismatic grids and h-refinement next to the rough surface and a
second-order accurate in space and in time CFD flow solver that can handle arbitrary
elements. Results obtained show that the friction coefficient can be predicted within

3.5%, but the Stanton number is under predicted from 8% to 15%.

ACKNOWLEDGMENTS

I would like to thank my advisor, Professor Tom Shih for his advice and
guidance, continuous encouragement, and patient support throughout this work. It has
been a very enlightening and fruitful experience to work with Professor Shih. My sincere
thanks are also extended to my co-advisor Professor Z.J. Wang and the other members of
my Ph.D. committee, Professor Harold Schock, Professor Mei Zhuang and Professor
Guowei Wei. I am also grateful for all the help granted to me throughout the course of
this work by the faculty and staff members, colleagues and many friends.

I would like to acknowledge the financial support by Michigan State University
through the Department of Mechanical Engineering. The icing project was supported by
NASA grant NAG 3-2576 and NNCO4GBZSG from NASA — Glenn Research Center.
The 3-D roughness study was supported from the DOE University Turbine Systems
Research (UTSR) Program. Dr. Wang provided the MUSIC solver and VisCART code;
Dr. Shock provided Gridgen license; Hudong Chen from Exa company provided
Powerflow license; Fluent company provided the FIuent-UNS code; Dr. Bons from BYU
provided real turbine roughness experiment data. The author is grateful for these
supports.

Grateful thanks should also go to my parents. Their love and encouragement are
the backbone of my persistence in the endeavor. Most of all I would like to thank my
wife, Xiaoyan. Her endless love is the source of my dedication through the completion

of this research.

IV

Table of Content

 

 

 

 

 

LIST OF FIGURES - - - - -- - - viii

LIST OF TABLES . .......... - xiv

NUMERICAL SIMULATON OF ROUGH-SURF ACE AERODYNAMICS ............ I

Chapter 1. Introduction . - -- - - - ..... - 1

1.1 Aircraft Icing ..................................................................................................... 1

1.2 Turbine Blade with Surface Roughness ............................................................ 4

1.3 Objective ........................................................................................................... 6

Chapter 2. Grid generation for 2-D Iced Airfoils - - - -- - ....... 8

2.1 Issues and Challenges ....................................................................................... 8

2.2 A Method for Generating High-Quality Single-Block Grids ......................... 11

2.3 A Method for Generating High-Quality Multi-Block Grids ........................... 15

2.4 CFD Simulation of Flow over Clean and Iced Airfoils .................................. 21

2.4.1 Problem description, formulation, and CF D flow solver ....................... 21

2.4.2 Results .................................................................................................... 22

2.5 SUMMARY .................................................................................................... 28

Chapter 3. Effects of Blocking Topology on Convergence Rate ................................. 32

3. 1 Introduction ..................................................................................................... 32

3.2 Problem Description, Formulation and CFD Flow Solver ............................. 33

3.3 Multi-Block Grids Evaluated .......................................................................... 34

3.3.1 Multi-Block Grids from a Single-Block Grid ........................................ 35

3.3.2 Different Partitions of a Multi-Block Grid ............................................ 38

3.3.3 Overlapping Chimera Grids ................................................................... 39

3.4 Results on Convergence Rate ......................................................................... 40

3.4.1 Multi-Block Grids from a Single-Block Grid ........................................ 41

3.4.2 Different Partitions of a Multi-Block Grid ............................................ 45

3.4.3 Overlapping Chimera Grids ................................................................... 47

3.5 Summary ......................................................................................................... 49
Chapter 4. Computing Aerodynamic Performance of 2-D Iced Airfoils with

Structured Grids---- - -- - - - - - 50

4. 1 Introduction ..................................................................................................... 50

4.2 Problem Description ....................................................................................... 52

4.3 Formulation and Numerical Method of Solution ............................................ 54

4.4 Results on Grids Generated ............................................................................ 55

4.5 Results on Solution Generated ........................................................................ 62

4.6 Summary ......................................................................................................... 67

Chapter 5. CFD Analysis of the Aerodynamics of a Business-Jet Airfoil with

 

 

 

Leading-Edge Ice Accretion 68
5.1 Introduction ..................................................................................................... 68
5.2 Objective and Approach ................................................................................. 69
5.3 Problem description ........................................................................................ 70
5.4 Formulation and Numerical Method of Solution ............................................ 72
5.5 Results ............................................................................................................. 76

5.5.1 Clean Airfoil .......................................................................................... 76

5.5.2 Airfoil with Glaze Ice ............................................................................ 79

5.5.3 Airfoil with Rime Ice ............................................................................. 83

5.5.4 Prediction of the Velocity Field ............................................................. 88

5.6 Summary ......................................................................................................... 94
Chapter 6. A Comparative Study by Using CFD to Predict Iced Airfoil

Aerodynamics ------------ - - 95
6.1 Introduction ..................................................................................................... 95
6.2 Objective ......................................................................................................... 97
6.3 Problem Description ....................................................................................... 98
6.4 Formulation and Numerical Method of Solution ............................................ 99

6.4.1 WIND and Fluent ................................................................................... 99
6.4.2 PowerFLOW ........................................................................................ 101
6.4.3 Grid Systems for WIND and Fluent .................................................... 102
6.4.4 Grid System for PowerFlow ................................................................ 105

6.5 Results ........................................................................................................... 106
6.5.1 CFD Predictions of 212 Rime Ice Shape ............................................. 107
6.5.2 CFD predictions of 944 glaze ice shape .............................................. 115

6.6 Summary ....................................................................................................... 123
Chapter 7. CFD Study of Turbine Blades with Real Rough Surface.-- - - 125
7.1 Experimental Setup ....................................................................................... 126
7.1.1 Roughness Surface Measurement and Fabrication .............................. 126
7.1.2 Wind Tunnel Facility ........................................................................... 128
7.1.3 Cf Measurement ................................................................................... 129
7.1.4 St Measurement ................................................................................... 130

7.2 Viscous Cartesian Grid Generation Approach .............................................. 131
7.2.1 Adaptive Cartesian Grid Generation .................................................... 131
7.2.2 Cartesian Grid Front Generation and Smoothing ................................ 132
7.2.3 Projection of the Cartesian Front to the Body Surface ........................ 132

7.3 Numerical Method ........................................................................................ 135
7.4 Results and Discussions ................................................................................ 138
7.4.1 Turbulent Flow over a Flat Plate in a Wind Tunnel ............................ 139
7.4.2 Turbulent Flow over Rough Surfaces in a Wind Tunnel ..................... 142

7.5 Conclusions ................................................................................................... 148
7.6 Future Works ................................................................................................ 149
Appendix Summary of the Sparlart-Allmaras Model ------------- 150

 

VI

Reference - - 153

 

VII

LIST OF FIGURES

Figure 2-1. Single-block gridm Streaks of clustered grids extending far into the flow
domain. (a) Close-up view showing source of streaks. (b) Overall view showing
amount of extension .................................................................................................... 9

Figure 2-2. Multi-block grid: grid lines clustered next to walls extend into the flow

domain along block boundaries. ............................................................................... 10
Figure 2-3. Geometry of NLFO414 airfoil with 623 ice shape. ........................................ 12
Figure 2-4. Grid for NLFO414 clean airfoil. (3) Overall grid; (b) Inner grid; (c) Inner grid

near the airfoil’s leading edge; ((1) Inner grid near the airfoil’s trailing edge. .......... 13
Figure 2-5. Equal arc-length concept used to generate smooth grid along the airfoil. ..... 13
Figure 2-6. Single-block grid around NLFO414 with 623 ice. ......................................... 14
Figure 2-7. Comparison of grid before and after smoothing ............................................ 14
Figure 2-8. A multi-block grid around NLFO414 with 623 ice. ....................................... 16
Figure 2-9. “Thick wrap-around” grid. ............................................................................. 17
Figure 2-10. “Transition-layer” grid. ................................................................................ 17
Figure 2-11. A-multi—block grid (MB2) with a “thick” wrap-around grid and transition-

layer grids for NLF 0414 airfoil with 623 ice shape. ................................................. 18
Figure 2-12. A-multi-block grid (MB3) with a “thicker” wrap-around grid and transition-

layer grids for NLFO414 airfoil with 623 ice shape. ................................................. 19
Figure 2-13. Geometry of Commercial airfoil with 145m ice shape. ............................... 20
Figure 2-14. A multi-block grid with "thick wrap-around" and "transition-layer" grids for

B757/767 airfoil with 145m ice. ............................................................................... 20
Figure 2-15. Lift coefficient as a function of attack angle for NLF 0414 clean airfoil. 23

Figure 2-16. Lift coefficient as a function of attack angle for NLFO414 with 623 ice shape.
................................................................................................................................... 24

Figure 2-17. Convergence history for NLFO414 with 623 ice with Single-block grid (SB—
2). .............................................................................................................................. 27

Figure 2-18.Convergence history for NLFO414 with 623 ice with Multi-block gn'd (MB-
3). .............................................................................................................................. 27

VIII

Figure 2-19. Predicted Mach number and velocity vectors near the leading edge of

NLFO414 with 623 ice. ............................................................................................. 29
Figure 2-20. Predicted Mach number and velocity vectors near the leading edge of

B757/767 with 145m ice. .......................................................................................... 30
Figure 2-21 Predicted Mach/Velocity near ice using different grids. ............................... 31
Figure 3-1. (a) NLFO414 airfoil. (b) 623 ice shape. (c) 623 iced airfoil ......................... 33
Figure 3-2. A single-block grid for NLFO414 with 623 ice ............................................. 37

Figure 3-3. MB grids from the SB grid shown in Figure 3-2. Solid lines indicate block
boundaries. ................................................................................................................ 37

Figure 3-4. MB grids with wrap-around grids from the SB grid shown in Fig. 2. Solid

lines indicate block boundaries. ................................................................................ 38
Figure 3-5. A multi-block grid for NLFO414 with 623 ice. .............................................. 39
Figure 3-6. Different partitions of the MB shown in Figure 3-5. ..................................... 39
Figure 3-7. Chimera grids for the NLFO414 airfoil with 623 ice. ................................... 40

Figure 3-8. Convergence history for SB, S-MB-2 and S-MB-3. (Top) or = 22°. (Bottom)
(1 = 52° ..................................................................................................................... 42

Figure 3-9. Convergence history for SB, S-MB-S, and S-MB-6. (Top) or = 22°.
(Buttom) a = 52°. ..................................................................................................... 43

Figure 3-10. Lift and drag coefficients for the SB, S-MB-Z, S-MB-3, S-MB-5, and
experimentally measured ones. ................................................................................. 44

Figure 3-11. Convergence history for Multi-Block grid. (Top) (1 = 22°. (Buttom) ct =

52°. ........................................................................................................................... 46
Figure 3-12. Lift and drag coefficients for the Multi-Block grids and experimentally
measured ones. .......................................................................................................... 47
Figure 3-13. Convergence history of SB, M-MB-l, and Chimera grids for or = 22° ....... 48
Figure 3-14. Lift and drag coefficients for the SB, M-MB-l, M-MB-Z , Chimeral, and
ChimeraZ. .................................................................................................................. 48
Figure 4-1. (Top) Commercial transport airfoil; (Buttom) 145m ice shape near airfoil’s
leading edge. ............................................................................................................. 53
Figure 4-2. A single-block grid (several different views) ................................................. 58

Figure 4-3. Multi-block 12MB1-SO ................................................................................... 58

Figure 4-4. Multi-block 2 (smooth grid): MBl-Sl. ........................................................ 59
Figure 4-5. Multi-block 3: MBl-SZ. ............................................................................... 59
Figure 4-6. Multi—block 4: MB2-82. ............................................................................... 60
Figure 4-7. Convergence history for the residual of CL and CD. See Table 1 for definition
of case numbers ......................................................................................................... 65
Figure 4-8. Mach number with streamlines for or = 6 degrees. ....................................... 66
Figure 4-9. Mach number with streamlines for (1 = 10 degrees. ..................................... 66
Figure 5-1.GLC 305 airfoil with 944 ice shape ................................................................ 71
Figure 5-2. GLC 305 airfoil with 212 ice shape. .............................................................. 71
Figure 5-3.Grid about GLC 305 clean airfoil. .................................................................. 74
Figure 5-4. Grid about GLC 305 airfoil with 212 ice. ...................................................... 75
Figure 5-5. Grid about GLC305 airfoil with 944 ice. ....................................................... 76

Figure 5-6. Computed and measured lift coefficient for clean GLC305 airfoil (M = 0.12,
P = 20.5 psi, Re = 3.5 x 106). CFD: WIND, S-A model. ......................................... 77

Figure 5-7. Computed and measured pressure coefficient for clean GLC305 airfoil (M =
0.12, P = 20.5 psi, Re = 3.5 x 106). AOA = 4. CFD: WIND, S-A model. ............ 78

Figure 5-8 Computed and measured pressure coefficient for clean GLC305 airfoil (M =
0.12, P = 20.5 psi, Re = 3.5 x 106). AOA = 6. CFD: WIND, S-A model. ............. 78

Figure 5-9. Computed and measured lift coefficient for GLC305 airfoil with 944 ice (M =
0.12, P = 20.5 psi, Re = 3.5 x 106). CFD: WIND, S-A model. ................................ 80

Figure 5-10. Computed and measured drag coefficient for GLC airfoil with 944 ice (M =
0.12, P = 20.5 psi, Re = 3.5 x 106). CFD: WIND, S-A model. ............................... 80

Figure 5-11. Computed and measured pressure coefficient for GLC305 airfoil with 944
ice (M = 0.12, P = 20.5 psi, Re = 3.5 x 106). AOA = 4°. CFD: WIND, S-A model
................................................................................................................................... 81

Figure 5-12. Computed and measured pressure coefficient for GLC305 airfoil with 944
ice (M = 0.12, P = 20.5 psi, Re = 3.5 x 106). AOA = 6°. CFD: WIND, S-A model
................................................................................................................................... 81

Figure 5-13. Effects of grid resolution and smoothing on CL and CD. 944 ice shape and
grid (M = 0.12, P = 20.5 psi, Re = 3.5 x 10°). CFD: WIND, S-A model. .............. 82

Figure 5-14. Computed and measured lift coefficient for GLC305 airfoil with 944 ice (M
= 0.12, P = 37.0 psi, Re = 6.0 x 10°). CFD: WIND (S-A) and Fluent. ................... 84

Figure 5-15. Computed and measured drag coefficient for GLC305 airfoil with 944 ice
(M = 0.12, P = 37.0 psi, Re = 6.0 x 10°). CFD: WIND (S-A) and Fluent. ............. 84

Figure 5-16. Computed and measured lift coefficient for GLC305 airfoil with 944 ice at
AOA = 4° and 6° (M = 0.12, P = 37.0 psi, Re = 6.0 x 10°). .................................... 85

Figure 5-17. Computed and measured drag coefficient for GLC305 airfoil with 944 ice
at AOA = 4° and 6° (M = 0.12, P = 37.0 psi, Re = 6.0 x 10°). ................................. 85

Figure 5-18. Computed and measured lift and drag coefficients for GLC305 airfoil with
212 ice (M = 0.12, P = 20.5 psi, Re = 3.5 x 106). CFD: WIND, S-A model. ........ 86

Figure 5-19. Computed and measured lift and drag coefficients for GLC305 airfoil with
212 ice (M = 0.12, P = 20.5 psi, Re = 3.5 x 10°). CFD: WIND, S-A model. ......... 86

Figure 5-20. Computed and measured pressure coefficient for GLC305 airfoil with 212
ice as a function normalized distance along the chord (M = 0.12, P = 20.5 psi, Re =
3.5 x 10°). AOA = 4°. CFD: Wind, S-A model. ...................................................... 87

Figure 5-21. Computed and measured pressure coefficient for GLC305 airfoil with 212
ice as a function normalized distance along the chord (M = 0.12, P = 20.5 psi, Re =

3.5 x 10°). AOA = 6°. CFD: Wind, S-A model. .................................................... 87
Figure 6-1. GLC305 airfoil with 212 rime icem .............................................................. 98
Figure 6-2. GLC305 airfoil with 944 glaze ice. [I] ............................................................ 99

Figure 6-3. Grid used by WIND for GLC305 airfoil. (a) Irmer grid. (b) Outer gn'd. 103

Figure 6-4. Grid used by Fluent for GLC305 airfoil with 944 glaze ice. ....................... 104
Figure 6-5. Grid for GLC 305 airfoil: (8) Inner block grid for 212 rime ice; (b) Close-up
of (a); (c) Inner block grid for 944 glaze ice; (d) Close-up of (c). .......................... 105
Figure 6-6. Cartesian grid method (left) and Variable Resolution regions (right) used in
PowerFLOW. .......................................................................................................... 106
Figure 6-7. 212 rime ice: lift coefficient = f (AOA). ...................................................... 108
Figure 6-8. 212 time ice: drag coefficient = f (AOA). ................................................... 108
Figure 6-9. 212 rime ice: lift coefficient = f (AOA). ...................................................... 110

XI

Figure 6-10. 212 rime ice: drag coefficient = f (AOA) ................................................... 110

Figure 6-11. 212 rime ice: moment coefficient = f (AOA) ............................................. 111
Figure 6-12. 212 rime ice: lift coefficient = f (AOA). .................................................... 113
Figure 6-13. 212 rime ice: drag coefficient = f (AOA) ................................................... 113
Figure 6-14. 212 rime ice: lift coefficient = f (AOA). .................................................... 114
Figure 6-15. 212 rime ice: drag coefficient = f (AOA) ................................................... 115
Figure 6-16. 944 glaze ice: lift coefficient = f (AOA). ................................................... 116
Figure 6-17. 944 glaze ice: drag coefficient = f (AOA). ................................................ 117
Figure 6-18. 944 glaze ice: lift coefficient = f (AOA). ................................................... 118
Figure 6-19. 944 glaze ice: drag coefficient = f (AOA). ................................................ 119
Figure 6-20. 944 glaze ice: moment coefficient = f (AOA). .......................................... 119
Figure 6-21. 944 glaze ice: lift coefficient = f (AOA). ................................................... 120
Figure 6-22. 944 glaze ice: drag coefficient = f (AOA). ................................................ 121
Figure 6-23. 944 glaze ice: lift coefficient = f (AOA). ................................................... 122
Figure 6-24. 944 glaze ice: drag coefficient = f (AOA). ................................................ 122
Figure 7-1. Geometry characteristic of surface #6 .......................................................... 127
Figure 7-2. Geometry characteristic of surface #4 .......................................................... 127
Figure 7-3. Schematic of the Flat Plate Wind Tunnel in Heat Transfer Measurement
Configuration. (Bons GT-2002-30198) .................................................................. 128
Figure 7-4. Dimensions (mm) of the test section with the location of six roughness paneélgs.
Figure 7-5. Geometry of the roughness surface section. ................................................ 129
Figure 7-6. Schematic of floating panel Cf measurement apparatus on AFUL wind tunnel.
................................................................................................................................. 130
Figure 7-7. Geometry of Roughness Surface #4 used in the CFD study. Roughness panel
mirrored in the stream-wise direction. .................................................................... 133
Figure 7-8. Cutting planes showing the viscous adaptive Cartesian grids. .................... 134

XII

Figure 7-9. The surface grid on the lower channel wall showing the refinement near the
rough panels. ........................................................................................................... 134

Figure 7-10. CFD simulation results compared to experimental data and standard
roughness correlations. ........................................................................................... 140

Figure 7-11. Convergence histories for flat panels with S-A and DES approaches ....... 1140

Figure 7-12. Comparison of CFD results for flat panels using S-A and DES approaches
with experimental and correlation data. .................................................................. 142

Figure 7-13. Convergence histories using the S-A model with the fine mesh for Surface

#4 ............................................................................................................................. 143
Figure 7-14. Comparison of Skin fiictions for roughness panels. .................................. 145
Figure 7-15. Comparison of Stanton numbers for roughness panels. ............................. 146
Figure 7-16. Velocity vector plot on a cutting plane near the roughness panel .............. 147
Figure 7-17. Pressure distribution on the rough surface (Surface #6). ........................... 148

X111

LIST OF TABLES

Table 2-1. Summary of Cases Simulated .......................................................................... 25
Table 2-2. Summary of Grids ........................................................................................... 25
Table 3-1. Summary of Multi-Block Grids Studied ......................................................... 35
Table 4-1. Summary of Grids Generated .......................................................................... 55
Table 4-2. Summary of Cases ........................................................................................... 61
Table 4-3. Predicted Drag and Lift Coefficients ............................................................... 61

Table 7-1. Comparison of experimental and computational results for rough Surfaces # 4
and #6. ..................................................................................................................... 144

XIV

Chapter 1. Introduction

Computational fluid dynamics (CFD) has advanced tremendously over the past
three decades. Today, CFD is used to understand the fundamentals of fluid flow
phenomena such as turbulence and two-phase flows through first-principle simulations.
It is also being used to analyze engineering designs in a wide range of applications that
involve fluid flow. Examples of these applications include automotive, aerospace, power
generation, propulsion, materials processing, bioengineering, and meteorology.

Despite the tremendous progress made in CFD and its wide use, CFD is still a
developing science. There are still many problems for which CFD cannot be applied or is
still unreliable. One area, where CFD still needs to be advanced, is flow and heat transfer
over rough surfaces. Two important applications involving rough surfaces are airfoils
roughened by ice accretion and turbine blades roughened by erosion and deposition. This
dissertation addresses advances needed for CFD to attack these two problems in a reliable
fashion.

The remainder of this chapter is organized as follows. In Sections 1.1 and 1.2,
roughness issues associated with the aircraft icing and turbine cooling are reviewed. In

Section 1.3, the objectives of this study are described.

1.1 Aircraft Icing

When an aircraft flies into environments where meteorological conditions can
cause ice to form on its wings, the aircraft’s ability to maintain flight will diminish
quickly with time unless there are ways to eliminate the ice formed. Compared to a

“clean” wing (i.e., a wing without accrued ice), an “iced” wing can have not only reduced

lift, but also stall occurring at markedly lower angles of attackm. This is because ice
buildup occurs mostly on the Wing’s leading edge, whose geometry is critical in
determining lift and stall. If the maximum lift force that can be generated by the aircraft
is less than its weight, then the aircraft will crash. Even if the lift is sufficient to sustain
flight, uneven ice buildup on the wings can produce flight control problems. Thus, it is
critically important to understand the different ice shapes that can form on the wings and
how they affect airfoil’s aerodynamic performance.

The icing problem has two components. One is ice accretion, and the other is
icing effect. Ice accretion is concerned with how ice forms on the wing surface and its
evolving shape as a function of Mach number, angle of attack, airfoil shape, airfoil size,
environmental parameters, and duration under favorable conditions for ice accretion.
Icing effects are concerned with the aerodynamic performance of an wing that has ice of
different sizes and shapes formed on it.

This study is concerned with computational fluid dynamics (CFD) simulations of
icing effects. Though icing effects can also be studied by flight and wind tunnel tests.
CF D is the most cost effective method and has the ability to simulate actual geometry and
flight conditions. However, its accuracy depends on the quality of the grid system in
resolving the geometry of accrued ice and the flow field and the “appropriateness” of the
turbulence model in capturing the relevant flow physics.

Because of the complexity in the geometry of wings roughened by ice accretion
and the flow induced by that roughness, most CFD studies on iced wings have been two-

dimensional (2-D), focusing on grid generation methods that can provide high quality

grids for rough surfaces and turbulence models that can predict lift and drag coefficients
as a function of the angles of attack [31'l7]’“°]'“°].

In order to generate a grid for an iced wing, one must first obtain the ice-surface
geometry. Currently, methods exist that can be used to measure ice wings in 2-D cross
sections, which involve meticulous tracing with subsequent digitization. Methods for
measuring three-dimensional (3-D) ice wings, however, are still being developed.
Traditional optical scanning methods for 3-D surfaces are not usefiil because ice is
transparent. Though a transparent surface can be made opaque by covering it with a paint
or powder, it is extremely difficult to cover a material that can melt with a thin layer that
is of the same thickness everywhere. As a result, reliable data of ice shapes on airfoils
are mostly 2-D. Also, even if the geometry of 3-D ice shapes can be measured precisely,
the computational requirement for computing 3-D iced wings as a function of angles of
attack is prohibitive expensive because the number of grid points needed to resolve the
ice geometry is extremely high.

With limited information on 3-D iced wings and the high computational cost in
simulating 3-D iced wings, there are very few 3-D CFD simulations of iced wings. When
3-D simulations are performed, the interest is on the qualitative features of the flow with
focus on 3-D effects [2]. So far, simulations that focus on iced Wing’s aerodynamic
performance have assumed the flow to be 2-D and steady.

2-D simulations that study how ice shapes affect airfoil performance have been
reported by a number of investigators [21.[3]. Several different types of ice shapes have
been studied. Validity of these studies was assessed by comparing computed results with

measurements from wind-tunnel tests. Results obtained show the 2-D CFD simulations to

provide reasonable predictions for some cases but not for other cases. There are several
reasons for the inconsistent prediction. One reason is that high-quality grid systems
become much harder to generate. The other reason is uncertainty in the turbulence model
used.

In summary, for iced wings, the focus has been on 2-D iced airfoils, which
remains an unsolved problem. There are two major research needs. The first is the need
to develop grid generation methods that can enable the generation of high-quality grids
for highly complicated ice shapes with multiple protruding horns. The second is to assess
the usefulness of existing turbulence models in predicting aerodynamic performance of

the iced airfoils as a function of angles of attack.

1.2 Turbine Blade with Surface Roughness

Land-based electric-power-generation gas turbines operate in harsh environments.
All surfaces such as blades, vanes, end-walls, and hubs that come in contact with the
combustor’s hot gases invariably become rough with service. The degree and the nature
of the roughness — due to mechanisms such as erosion, fuel deposition, pitting, and
spallation of thermal-barrier coatings — depend on the environment fi'om which the air is
ingested, the operating conditions, the effectiveness of cooling management such as film
cooling in maintaining material temperatures within acceptable limits, and the service
time. This material degradation in the form of surface roughness is known to increase
surface skin friction and heat transfer in a significant way. For a given cooling
management, increase in surface heat transfer increases material temperature, which

hastens further material degradation. In order to estimate the service life of turbine

blades, it is necessary to estimate the skin friction and heat transfer augmentation
generated by various types of surface roughness.

The importance of this problem has led a number of investigators to study it by
using both experimental and computational methods. Since the geometry of surface
roughness can vary widely, most investigators study rough surfaces by modeling it with

regularly or irregularly simple shapes. Simple shapes studied include distributed

[37],[38] [39],[40] [41],[42]

cylinders , spherical segments , cones , and pedestals [43]. The
usefulness of these studies is unclear because the connection between artificial and real
roughness has not been made clear. Bons and co-workers [441145] did study real rough
surfaces. Their study showed that tradition sandgrain models do not predict correctly
skin friction or surface heat transfer. This is because some parts of the roughness extend
significantly into the flow when compared to the boundary-layer thickness.

Thus, there is a need to understand in detail how realistic roughness affects skin
friction and surface heat transfer. This could be accomplished by performing simulations
based on the Navier—Stokes equations that resolve the detailed geometry of the roughness
and the flow induced by the roughness. Such direct numerical simulations (DNS) of
surface roughness can be validated by a recent experimental study performed by Bons

and co-workers [45]

, where detailed measurements of Cf and St were made. The
understanding gained from validated DNS type simulations can be used to guide the
development of mathematical models to predict skin friction and surface heat transfer in
CFD simulations that do not resolve the geometric details of roughness.

The main challenge of doing 3-D DNS simulations of surface roughness is the

need to resolve the highly disparate length scales in the physical problem. First, there is

the model length scale, L, in the present case the length of the wind tunnel. Second, there
is a feature length associated with a typical roughness element, 1, which is 3 or 4 orders

less than L. Lastly, there is the thickness of the viscous sub-layer, 6, which is 5-6 orders

less than L. It is critical that all three length scales are resolved in the CFD simulation.

The disparity between L and 5,, can be efficiently captured by one-dimensional grid

clustering near solid walls in the wall normal direction. For the computational grid to
capture rough elements, the grid size along the tangential direction of the rough wall must
be comparable to 1. If one uses a non-adaptive structured grid, and employs 50 layers in
the wall normal direction, the grid size is roughly 1,000 x 800 x 50 (4 x 107 points),
which will overwhelm most computer clusters. A far more efficient computational grid
for this type of geometry is an unstructured adaptive grid. It appears a recent developed
adaptive grid approach named viscous adaptive Cartesian approach [4°] is the most
suitable to tackling the challenge. A finite volume Navier-Stokes solver [47] capable of
handling viscous adaptive Cartesian grids was employed to carry out the computational
simulations with both RAN S [25] (Reynolds-averaged Navier-Stokes) and DES (Detached
Eddy Simulation) approaches to model flow turbulence. The DES [481149145 °] approach is
a hybrid RANS/LES (Large Eddy Simulation) approach designed to capture both the
boundary layer and large separation regions. The reason for using both RANS and DES
approaches is to compare the computational results of both with experimental data to see

which achieves a better agreement.

1.3 Objective

The objective of this study is to address computational issues associated with

CFD simulations of flow over surfaces with roughness in which the detailed geometry of

the roughness is resolved. This will be accomplished by addressing challenges in

performing CFD simulations of two important engineering problems, where roughness

play a critical role — the aircraft icing problem and turbine erosion/deposition problem.

The detailed objectives are as follows:

1.

Develop grid generation techniques and approaches to enable the generation
of high-quality single-block and multi-block structured grids for 2-D rime and
glaze iced airfoils. Rime ice shapes are characterized by surfaces that are
rough and jagged but have no protruding horns. Glaze ice shapes are
characterized roughness and jaggedness as well as having two or more
protruding horns near the airfoil’s leading edge.

Study and propose how best to partition blocks for fastest convergence rate to
steady state for multi-block grid systems.

Study the effects of turbulence models in predicting lift and drag of 2-D lime
and glaze iced airfoils as a function of angle of attack and examine when CFD
can predict iced airfoil aerodynamics reliably.

Study flow over 3-D roughness that represent a “roughened” turbine surface

and explore usefulness of CPD in computing such flows.

Chapter 2. Grid generation for 2-D Iced Airfoils

In this chapter, methods are presented to generating high-quality single and multi-
block structured grids for iced airfoil that can be applied to both rime and glaze ice

shapes.

2.1 Issues and Challenges

The grid system used in a CFD simulation must resolve the geometry and enable
the governing equations with the turbulence model to resolve all of the relevant physics
with minimum grid-induced errors. For airfoils with ice formed on them, this is a
challenging task. From what NASA has learned from previous studies {21131, structured-
grid approach should serve as a good reference or baseline, against which other grid
approaches can be compared and assessed. If structured grids are to be used, then there
are two main issues.

The first main issue is surface preparation. This issue arises because the geometry
of iced airfoils is complicated not just with protruded horns and feathers, but also with
small-scale surface roughness. One system that has been used to prepare ice-surface
geometry with various levels of smoothing is NASA’s Turbo-GRD code [511°]. A
comprehensive system currently being developed and maintained at NASA Glenn
research center is SmaggIce [7]_ SmaggIce 1.0 is a software tool kit that provides
interactive ice surface preparation for grid generation and ice shape characterization. In
this study, SmaggIce 1.0 was used to prepare the ice shapes through its smoothing

routines.

The second main issue is blocking topology. When using structured grid systems,
two choices are possible: single-block grids and multi-block grids. Each approach has
its advantages and disadvantages. For single-block grids, precise controls are needed to
negotiate the series of convex and concave surfaces, while maintaining proper grid aspect
ratio, orthogonality, smoothness, and grid alignment with the flow. This may require
extensive internal blockings that can later be combined. One problem for single-block

grids of iced airfoils that has not been resolved is illustrated in Figure 2-1 [8].

 

Figure 2-1. Single-block gridial Streaks of clustered grids extending far into the
flow domain. (a) Close-up view showing source of streaks. (b) Overall view
showing amount of extension

In this figure, concave regions can be seen to cause grid lines to cluster, forming
streaks that extend far into the flow field. These streaks can degrade the accuracy of
solutions. For multi-block grids, the major problem is illustrated in Figure 2-2. Grid
lines that are highly clustered next to solid walls propagate into the interior of the domain
along block boundaries. Since these clustered grids next to walls typically have very

high aspect ratios (from hundreds to hundreds of thousands), they can induce

considerable errors in the solution. If this problem is not resolved, then it is impossible to

generate a hi gh-quality multi-block structured grid for iced airfoils.

-- -_.—n _— —..
- --'_--'_._.—'-__F:v==._—_r.‘-_l

.' .v.‘ —’=*—— ”Iv-”r: ,_ 4—I—'—’-_-: :_~_-_—_’ :::::

9000' -.--._._.—-—-‘._.—-r,_.

55=§§ s
-. —_-;> ‘

     
 

block
boundary

ll _.~-
‘.l:-"”-o’__,_wii.j‘ ‘
tau-s (’- ' A.
-‘.°,‘_',_'._",- " ‘
*o‘. ‘ r a
.- ‘ \

Figure 2-2. Multi-block grid: grid lines clustered next to walls extend into the flow
domain along block boundaries.

Another problem that plagues both single- and multi— block structured grids is
small-scale surface jaggedness or roughness. If hyperbolic or algebraic methods are used
to generate grids, then these surface discontinuities propagate into the domain interior,
creating a poor quality grid.

The focus of this study is to develop methods to overcome the aforementioned
problems associated with single- and multi- block grids. All methods developed are
rooted in blocking strategy, and they are presented in the next two sections.

Before presenting the methods developed in this study, it is noted that once the
blocking topology has been established, many methods are available to generate high

quality grid systems. In this study, all grids were generated by using a four-boundary

10

technique based on Hermite interpolation. This transfinite-interpolation method has
controls that enable multi-dimensional clustering, enforce orthogonality or any other
specified intersection angle at block boundaries, and ensure C1 continuity across block

boundaries [91410].

2.2 A Method for Generating High-Quality Single-Block Grids

To illustrate the method developed in this study to overcome the problem shown
in Figure 2-1, consider the NLFO414 airfoil and the airfoil with 623 ice shape as shown in
Figure 2-3. The method developed involves the following three steps: 1. Generate a grid
for the clean airfoil irrespective of the ice shape. Since the clean airfoil is smooth, the
resulting grid is smooth throughout the flow domain. A grid thus generated for the
NFL0414 airfoil is shown in Figure 2-4. Note that this grid is made up of two parts. One
is a fine grid next to the airfoil (701 x 91), extending at least 0.6 chord length from the
airfoil in all directions (referred to as the inner grid). The other is a coarser grid (125 x
21) that overlaps the fine grid by 0.1 chord length and extends at least 15 chord lengths
away from the airfoil in all directions. The inner grid is the single-block grid of interest.
While generating this grid, grid lines were clustered next to the airfoil surface so that the
first grid point is within a y+ of unity. Along the airfoil surface, equal arc-length concept
was employed to create as smooth 3 grid as possible (see Figure 2-5). 2. Choose one set
of grid lines near the airfoil surface, but a small distance away from the iced airfoil. All
grid lines between the selected set of grid lines and the iced airfoil are regenerated to
account for the accrued ice shape. Thus, the selected grid lines form two blocks. For the
block that contains the iced airfoil, many sub-blocks are created in order to generate a

high-quality grid. Figure 2-6 shows a grid generated by using this approach with the

11

blocking topologies shown in different colors. With this blocking topology, all grids
generated in all blocks can be combined into a single-block grid. 3. Repeat Step 2 until a
high-quality grid is generated. This may be done in a solution adaptive manner. 4.
Merge the the two blocks separated by the selected grid line to form a new single block.
5. Apply elliptical smoothing to the grid generated above to improve the smoothness of
the final grid.

As shown in Figure 2-6, the method just described can generate good quality
single-block grids without streaks of highly clustered grid lines propagating into the flow
domain. With this approach, only grid lines very close to the iced airfoil are affected by
the ice’s geometric complexities. The method just presented for iced airfoils is based on
an idea proposed by Tai [I ”-

Figure 2-7 shows the graph of the grid before and after elliptical smoothing, from

this figure it can be seen that elliptical smoothing can effectively improve the smoothness

>
fix.

of the grid.

 

NLFO414 airfoil.

623 ice shape neaL/P

airfoil’s leading edge.

.\_-

Figure 2-3. Geometry of NLFO414 airfoil with 623 ice shape.

12

 

   
    
     

   

W%w%\§w“\wwww
$3 °~R®§W\\.\\‘ L :
E.-\R\\$5§m§~\\~x\\\\\\\\~§\

\xmmwmvww
~ n
\"R‘ w$\m ‘ \\\\\‘I\\
(e *RX‘V“ sevwx‘ww mum-xx
‘ ~‘ ~ ‘ ~ \‘ \ ~ *~ ea use» . \\\\\\\\\\“‘\V\\\
w *'§ §\ \®‘\.3~~ “ \avcnmr \“ "\\\\\\\\\N\‘Q VA
&~.~ s‘ “\ W¢$¢m¥kfit~mStu-mitt???“ no“...
- s \e‘s . asses»“‘«w‘rm-swn-w ~~
.~ a .o n
.s

        
  
     

    
  

     

  

         
  

           
      
  
  

 

    
   
  
  
 

  

   

            
    
  
                 

  

   
   
   
  
   
   

 

 

 

 

        

,. 9—535-v
’ ’ o o o‘o ,
.4???»
”’%%%$$'o'o o
.«59»;' '0: ’v '4
’ ’ ’ 9%6 , x», a" u o.
’rw’g’vfiv vow ’lye n
, , , ’ ,9 0, ., , ,...,. mm...
’9’ v“ ’ $4943 , av - 4117/}; momma
’6. ’ , d”; 4 ” , ,, ll.[dwwma’y/"N'flMl/tluuflhrzmz'"
’ ,’ ’ ” I I I, ’1; Iv” ,W/I”, Ifl‘llllllllllav «any
/‘ ”$43446 “7.x “” I’mwmmgmzw
. 11 , , / 1,:
4,4M, ,- az~2zmmflflfigmmm
”. ////4,4 5% «Wm! / xix/Mr.
1 I ' 11' . ’.’//,./. ; AI? /V/l/Ivr’

 

 

 

(C) (d)

Figure 2-4. Grid for NLFO414 clean airfoil. (a) Overall grid; (b) Inner grid; (0) Inner
grid near the airfoil’s leading edge; (d) Inner grid near the airfoil’s trailing edge.

, . . . . . . ilits."‘..‘.,ri‘lu.s‘ltw'lllll.‘.."lllIllllllllllllllllllllll
ngsfifisq“ . . - ' n ' I “‘Illilllllll"“lll"‘ll'l|‘llllll
3N“.- _ . , .. ‘ 1, “‘“CW‘ it"l ‘ 'lilll

.l‘

‘u l l !’
l‘t'dl‘r'r Wit 1 l I ll 1 ‘.

. I ‘
l‘l'l‘l'u‘ ' IiI'IIIIPI'I'I'IIlllll’lllmllli»! l I

 

Figure 2-5. Equal arc-length concept used to generate smooth grid along the
airfoil.

\\\\\‘:\“\“\Z‘~‘\\\\‘““““““““““““ \‘\““‘
\x\ \ \\h\\
\\\\\~‘ Q--.

\\\\‘ \\
“\‘.““ \“u‘. u \
§\§\§\‘ ‘“ “‘ l:‘\““‘ \‘n‘
\x ‘\ ‘
- .-. :_~ \ . \“5‘33
‘1‘“ \

   

4' wzzfl: :21: .. .-.
/ 4 MM" éfrwn”afl u on. . an ”L": an
z/i zap?”uw"zoplun agar” m???" fill": .mfilggwu "In";
" III
4” 4’ millilmuwut‘ m,‘ Zuuumll" ’" u

w, [I /'

 
  
  
 
 
    

   

   

 

II

II I I
ill/ll Il’fl’
”Ii/”7M"
ng/mliiiflrfg'fér’.”

I ,1, r
$5342.:

       
   

lut'llnli" .

.ul”‘.’.will
a

r .1- w
.a‘ “a“ ‘

mm .a
’55:“ \“ \‘fiu

a, m.\‘\\'.‘
\.\
W I.

\

.‘

‘r‘ 'l
“lit .

Figure 2-7. Comparison of grid before and after smoothing

 

2.3 A Method for Generating High-Quality Multi-Block Grids

Even though multi-block grids are generally considered to be superior to single-
block grids, there are problems such as the one illustrated in Figure 2-2. To further
examine problems that may arise from the use of multi-block grids for iced airfoils, a
multi-block grid was generated for the NLFO414 airfoil with the 623 ice. This grid is
shown in Figure 2-8, and two problems can be discerned. The first problem is the
prOpagation of highly clustered grid lines into the domain interior along block
boundaries, which was mentioned in relation to Figure 2-2. This problem will occur
whenever a block boundary next to a wall extends into the interior of the domain. The
second problem shown in Figure 2-8 is that roughness of ice surface propagates into the
interior of the domain. This problem affects both single- and multi-block grids.

To overcome the first problem, a “thick wrap-around” grid is proposed as
illustrated in Figure 2-9. Basically, one layer of grid with clustering next to the wall is
wrapped around the iced airfoil. The thickness of this wrap-around is such that grid
spacing from contiguous blocks will be comparable in size at all block boundaries.

To overcome the second problem mentioned, a “transition-layer” grid is proposed
as illustrated in Figure 2-10. With a transition layer, surface discontinuities are confined.
To demonstrate the usefulness of the “thick wrap-around” and the “transition-layer”
grids, multi-block grids were generated for the NLFO414 airfoil with the 623 ice (Figure
2-3) and the B757/767 airfoil with the 145m ice (Figure 2-13). Figure 2-11 and Figure
2-12 show two different multi-block grids for the NLF 0414 airfoil with the 623 ice. They
differ in the thickness of the wrap-around grid and in the block boundary near the

airfoil’s trailing edge. From these two figures, it can be seen that if the wrap-around grid

15

is of the right thickness, the problem shown in Figure 2-2 can be eliminated. From

Figure 2-12, the usefulness of the transition layer grid is evident.

 

i
III"
.I

h

    
   
  
   
  

.
$118“

that
Jo“

n“
“:5“:
9‘“
I"

‘3
er
.‘

fi/II/ 4" ‘
Until/ea, "1 "
flirt“ 'l‘lill1lllllllliiiIt!"iti’iiiiiiiiii’tifithgz
.1ij:‘ttts'ttitiriil'iimm"fem"
tuptgttilimjm

3:
‘ It" ”I!"
Iiillilll“‘““lilriii""'"tn”

JlIIlI
lrllll
:II I

     
  
 
  

   
 
 
 
  
  

 

ref/”mt

ill/Iii
lllll

  
   
   
 
   
    
  
      

 

%
“iii Mill f

I ", 55:51:“,14'F
’ r], ' t It I, #2 .
//’, 'fit/Iyé’éf ”
,ilm |l||||l|||| /

   
  

  
 

Amplified view 1

  
 

      

4%. /
WW%% 4 v»
‘7, |||ll||ll||l|ll/////{|/W ”If” I
{it //I I I II in III:
7 flit”! if Ill-Ii .IJIm
II II

it .. ////l / ‘9
”WI/ill lilti’l’i’i’lllf'liii
itWI'M/"”I”’I’i"‘iiiii’iit"ii ““ll
[gludllllllll :i””’”gl' 5” I!”
”W55!” I "I" I

  

I ' I
'i’” I‘ll “III
5, I

 

I n 'l
’ Ilium"!!!
’1 I

Figure 2-8. A multi-block grid around NLFO414 with 623 ice.

           

l6

  
  

    

     

  
       

 
  
   
   
  
  

 

      
    

I | IIIIII fill
:Efililiiiiirnlmm} I. I n
-- .“ .l'll “I .
..::.i:IuIII:"Hflm III:
==.~‘=‘III“I |||||
:-=-!~:‘ :llll .
=_-=§\§ .

E §§§§'

= x =

.5- \\

    
     
 

    

‘5) int
\\\\\\\\R§§\‘\\ I

“x \\\\‘ \\ ‘
§§k\\\\\II|\\‘

 
         

   

Wrap-around grid

   

\\\c 7“. fl
‘

.“3
..\\

    

\ “um-

  

    

\\

ll

1

/
,///
y/II

Figure 2-9. “Thick wrap-around” grid.

 
  
   

    
    

 

      
     

          

    
 

 
          
 
 

      

        

        
 
 

 
   
 
     
  

 
   

  
  

  
  
     

  

, W. . §§§§g§sc
lav ,, , ; ~5§§§5§§a=
in Witt? we” ,. 5% --E§“

, MI / //,i/,/_z//,,/anzf,7777,,/.o,my; 1%, ~ ~ .~..‘
l W mm! u/V/I/ flu/[MrW/‘W’ (‘1 ~~g~a--~i
" ”FA/[fl]: /i ”gig; 7“ WIN” 'II’I/y/LQZI/[éf/gi; _ .s#5:5:555:
H H Ml. ' ‘. “7,,”1‘ ". U , I
,i , ii , l W dill efiafl~5a3==
, .//, / fl////it,.,../,l,l, . ~- ~ ~

/ fl ////, [/14], I [HI/”I’M”, '/f l' ' [ill/'y/r [I ---~

.7617 M, WWI/l, W ', Fl ”(if in n

Figure 2-10. “Transition-layer” grid.

 

a. g ,
. 5w?
a a a; 5
aa’aannw’mg—uv a

a 5 up. 1%?

0‘
a; a

 

        

     

 

   

     

  

              

‘ ‘0
‘ ‘ o. .
«v.5...
‘Qfi‘eno conulawéé,
\h‘ . c...n..u.m¥.,.w
. 2.
$fi....va.w7v.7

77.
. 7.x.» X77777”)
: 77w {xx/m/VXXM

: 7X7 /\.,\, X /.vm/..VV,

\

 

th a “thick” wrap-around grid and

2-11. A-multi-biock grid (MB2) wi
transi ion-layer grids for NLFO414 airfoil with 623 ice shape.

Figure

18

                     

    
              

   

:2
“77..., 7. .. 77.... u
47/ IIIooIOOONoNoNnI
, , ””9"";auvonnnunuuuunu. a
277%.»..nmnuuonounuouuuno... :
r v 0 O. O 0 o O~
5 ’0 O O 0.000 QO.OQOOOOOOO A
"6"HON“NNONO”NMNNNMNH”MHMWMOHH§J¢J¢I
I 0000!... 90.6.
zuzaoorooooooouoiiunouuunouuuuououcnuiuoo
I I O‘vIOIDOIOOOO OOOOOOOO 9 O
o 71.. .flio: «01.0100 0100001.! d
.7 wflllflfiaflfi :o/oooooo
11H. il'lllitllllo IOOIOO .
7

77/1

         

         

  

            
    

         

2V

,,/,
’a
I

   

   

II
III
/

   

.74 17/
$1 «.7 z i
2.5;? 7; ’7‘
47 V9
2.53647, 6%” ,

r 6
9:2???”

    

      

    

7. 4 /. / \ / > . 7
777% ,. . , .
a»: 7 , Vow , . a... «a.»
«7, 7 7.. ... 55.53
7 7 7.... 5...;
_ - - ail/fl?
5 22

                 
            

  

/

/. ..
, 7 wkayfl/

        

Figure 2-12. A-multi-block grid (MB3) with a “thicker" wrap-around grid and

transi ion-layer grids for NLFO414 airfoil with 623 ice shape.

19

Figure 2-14 shows a multi-block grid for the B757/767 airfoil with the 145m ice.
The high-quality grid generated for this complicated ice shape with multiple significantly
protruded horns is another testament of the usefulness of the “thick wrap-aroun ” and

“transition-layer” grids.

 

Commercial Transport
B757I767 airfoil

145 ice shape near /

airfoil’s leading edge

Figure 2-13. Geometry of Commercial airfoil with 145m ice shape.

 
   

  

“Hill iii” “110‘s“

 
   

  
   
    
    
 
 

“1.7%“. ”V-

“Hailing”: , 7 ~
ll.“ .7.
will!” . I [/1’ z. ‘0’;
0757577 $3.}:

 

I
'1
’1

 
      
   
       
    

I
Ill".

, 7
”"1. {1' ‘0
I:I'l$'0a','6'
. la}: 0,,

   

Figure 2-14. A multi-block grid with "thick wrap-around" and "transition-layer" grids
for 8757/76? airfoil with 145m ice.

20

2.4 CFD Simulation of Flow over Clean and feed Airfoils

In the previous two sections, we presented two methods for generating high-
quality structured grids over iced airfoils. One is a single-block grid approach for
moderately complicated ice shapes. Another is a multi-block grid for highly complicated
ice shapes. In this section, we show the CFD solutions generated on grid systems
generated by these two methods. In the following, we first describe the flow problem, the

formulation, and the flow solver. Finally, we present the results.

2.4.1 Problem description, formulation, and CFD flow solver

Two clean and two iced airfoils were simulated in order to contrast the effects of
ice accretion. The two clean airfoils investigated are NFLO414 (Figure 2-3) and
B757/767 (Figure 2-13). The two iced airfoils investigated are NLFO414 airfoil with 623
ice (Figure 2—3), and B757/767 airfoil with 145m ice (Figure 2-13). The detailed
geometries of these airfoils and ice shapes can be found in Ref. [1]. For each clean and
iced airfoil, a series of angles of attack were simulated. When the angles of attack were
changed, the grids were not changed. The flow conditions employed for all simulations
were as follows: The approaching Mach number was 0.29. The freestream temperature
was —5° C, and the Reynolds number based on the freestream conditions and the chord
length was 6.4 x 106.

The flow past the clean and iced airfoils was modeled by the ensemble-averaged
conservation equations of mass (continuity), momentum (full compressible Navier-
Stokes), and energy for a thermally and calorically perfect gas. Turbulence was modeled

by the two-equation shear-stress transport (SST) model. Wall functions were not used.

21

Integration of the conservation equations and the SST model were to the wall, where no-
slip and adiabatic wall conditions were imposed.

Though a number of open-source codes are available, one of the most versatile
and well supported code is the WIND code “21’“3]. In this code, the convective terms
were approximated by second—order upwind differencing formulas. Time marching to
steady state was accomplished by an implicit method based on ADI-type approximate

factorization.

2.4.2 Results

A summary of all cases simulated is given in Table 2-1. Summarized in Table 2-2
is the number of grid points used for each simulation. Below, we first describe the
validation of the CFD simulations. Afterwards, we examine the effects of the grids on

the solutions and the convergence rate to steady state

2.4.2.1 Validation

Of the simulations performed, only the NLFO414 airfoil with and without the 623
ice have experimental data, which can be used to validate this CFD study. The
experimental data is on the lift coefficient as a function of the angle of attack, and was
obtained at NASA Glenn’s Icing Research Tunnel and at NASA Langley’s LTPT facility
[15].

A comparison between the CFD predicted and experimentally measured lift
coefficient for the clean NLFO414 is shown in Figure 2-15. This figure shows that CFD
is able to predict the lift coefficient quite well up to stall. After stall, however, the CF D

predictions do not match the experiments. The reason for the mismatch after stall is

22

unclear. But, the excellent agreement before stall, at stall, and slightly after stall is

encouraging.

 

1.8
1.6 -
1.4 r

1.2 - /

1.0 r /x

CL
\2.

0.8 r x/
0.6 ‘ /
0.4 ‘ X

0.2 i x

 

 

0.0 T f T l l T l l T l l I 1

-6-4-2 0 2 4 6 810121416182022
AOA(deg)

 

Figure 2-15. Lift coefficient as a function of attack angle for NLFO414 clean airfoil.

Figure 2-16 shows a comparison between the CFD predicted and experimentally
measured lift coefficient for the NLFO414 airfoil with the 623 ice. In this comparison,
there are two types of ice shapes in the experiments: 3-D (actual 3-D 623 ice shape) and
2-D (i.e., same 2-D ice shape for the entire airfoil span; same as CFD). The 2-D ice
shape selected is the one at the mid-span. From this figure, it can be seen that the
agreement is good at low angles of attack. Though agreement is not as good at higher

angles of attack when compared to the LTPT data, the stall angle was well predicted

23

1.0

 

.

l 7
.
I
.

,4 ‘_ .
g 7' I‘» \
’I!

     
 

 

 

 

l \0.
0.6 “ l I ““X” ”53—“
l i l
d 04 ’ ~X—CFDSB-1—-~
~A~CFD 83-2
0.2 : -<>~ IRT ,
l -B- LTPT (2-D)
0.0 i . “°"L:EI-_(§:|Q)
-o.2 I 1 7 7 1 1 .
-4 -2 O 2 4 6 8 10 12
AOA (deg)

Figure 2-16. Lift coefficient as a function of attack angle for NLFO414 with 623 ice
shape.

The number of grid points in the outer block is 125 x 21. The comparison given
above was based on CFD simulations that used in single-block grid, SB-l, and SB-2 (see

Tables 1 and 2). The generally good agreement obtained gives some confidence to the

CFD simulations performed.

24

 

 

Table 2-1. Summary of Cases Simulated

 

 

 

 

 

 

 

 

 

Airfoil/Ice Shape Grid Type" Angle of Attack
NLFO414/ Clean SB 0'43653371063261533123132561 5'
NLFO414/ 623 SB-1 1.1, 6.2, 10.2
NLFO414/ 623 88-2 0.0, 8.0
NLFO414/ 623 MB-1 0.0
NLFO414/ 623 MB-2 0.0, 6.0
NLFO414/ 623 MB-3 0.0
8757/767/ Clean SB 2.0, 7.0, 14.0
B757/767/ 145m MB 0.0

 

 

 

*Refers to inner block (Figure 24); $8 = single block, MB = multi-block
(MB-1: Figure 2-8, MB-2: Figure 2-11, MB-3: Figure 2-12)

 

Table 2-2. Summary of Grids

 

 

 

 

 

 

 

 

 

 

 

Airfoil/Ice Shape Grid Type Grid Number"
NLFO414/ Clean SB 701 x 91
NLFO414/ 623 88-1 781 x 702
NLFO414/ 623 88-2 136 x 101
NLFO414/ 623 MB-1 131,197
NLFO414/ 623 MB-2 98,898
NLFO414/ 623 MB-3 1 1 1,594
8757/767/ Clean SB 977 x 101
B757/767/ 145m MB 97,070

 

*lnformation is for inner block. The number of grid points in the outer

block is 125 x 21.

 

25

 

 

2.4.2.2 Effects of multi-block grids on convergence rate

Though “thick wrap-around” and “transition-layer” grids introduced in this study
enabled the generation of high-quality multi-block grids, there is one more issue
regarding the use of multi-block grids. This issue has to do with the CFD code. There
are two types of CFD codes. One type codes are written around the multi-block grid
structure. With these codes, computations are performed in one block at a time until all
blocks are computed. This process of computing in one block at a time is repeated until a
converged solution is obtained on all blocks. The other type of CFD codes do not see the
block boundaries when generating solutions. In those codes, each grid point or cell
knows who its neighbors are, and does not use any information on the multi-block grid
structure. These codes effectively see a single block grid, though they may perform
domain decomposition for parallel computing.

In this study, the WH\ID code was used to obtain solutions. The current version of
the WIND code belongs to the first type of CFD codes mentioned. Basically, it does
computations in one block at a time until all blocks are computed, and then repeats until
convergence. For all of the simulations summarized in Table 1, the WIND code
converged very quickly if the inner grid was a single-block grid. Figure 2-17 shows the
convergence history for the case with NLFO414 airfoil, 623 ice shape, and SB-2 grid.
From this figure, it can be seen that the second norm of the residual (L2) drops more than
four orders of magnitude in about 10000 iterations. Each iteration means all blocks have
been computed once. The convergence history shown for this case is fairly typical of all

single-block grids.

26

 

 

 

 

 

 

 

Zone2
a?
:L
(D
O
_J
_8 l l l 4
0 2 4 6 8 10

Iteration (1 03)

Figure 2-17. Convergence history for NLFO414 with 623 ice with Single-block grid
(SB-2).

 

 

 

 

 

—Zone1 ----Zone2 --Zone3
'7 L ---Zone4 ---~Zone5
‘8 r r r 1
0 5 10 15 20 25
Iteration (103)
Figure 2-18.Convergence history for NLFO414 with 623 ice with Multi-block grid
(MB-3).

27

When the inner grid is a multi-block grid, then the convergence has been either
extremely slow or does not converge as shown in Figure 2-18. The culprit is the data
transfer across block boundaries in critical regions such as flow separation. Though the
solutions did not converge for cases with multi-block grids, reasonable flow fields were

obtained as shown below.

2.4.2.3 Flow induced by ice shapes

Figure 2-19 shows the Mach number contours with velocity vectors in the region
about the leading edge for the NLFO414 airfoil with 623 ice predicted by using SB-2 (a
.= 0), MB-2 ((1,: 0), and MB-3 ((1,: O). This comparison shows that predicted flow
features are similar even though the multi-block-grid solutions did not seem to converge.
From this figure, it is clear that the two horns of the ice cause the formation of two large
separated regions. Figure 2-20 shows the very complicated flow induced by three highly
protruded horns (B757/767 with 145m ice at q = 0). From this figure, it can be seen that
horns, especially highly protruded ones, can severely disrupt the flow over the airfoil.
Thus, it is no wonder that lift can drop so significantly and stall can occur at much lower
angles of attack. Figure 2-21 shows the predicted Mach/Velocity near ice by using

different grids, it can be seen these grids will give reasonable flow field prediction.

2.5 SUMMARY

A method was developed to enable the generation of high-quality single-block
grids for moderately complicated ice shapes. "Thick wrap-around” and “transition layer”
grids were introduced to enable the generation of high quality multi-block grids for iced
airfoils with highly complex ice shapes. For multi-block grids, the convergence issue

needs to be addressed for some codes.

28

 

Figure 2-19. Predicted Mach number and velocity vectors near the leading edge
of NLFO414 with 623 ice.

29

 

Figure 2-20. Predicted Mach number and velocity vectors near the leading edge
of B757/767 with 145m ice.

30

 

Figure 2-21 Predicted MachNelocity near ice using different grids.

31

Chapter 3. Effects of Blocking Topology on

Convergence Rate

In this chapter, we examine the effects of blocking topology on convergence rate

to steady state.

3.1 Introduction

On grid generation, in Chapter 2, we presented a number of methods to generate
high-quality single- and multi-block structured grids for complicated ice shapes. For
moderately complicated ice shapes, single-block grids are possible. But, for really
complicated ice shapes, multi-block grids are needed. In Chapter 2, we demonstrated that
a thick wrap-around grid is needed to ensure that grid points clustered next to solid
surfaces do not propagate into the interior of the flow domain. We also suggested using a
transition layer next to solid surfaces that which rough and jagged to confine and smooth
the effects of the irregular geometry on the grid.

Though we were able to generate high quality multi-block grids for highly
complicated ice shapes, convergence rate to steady-state solutions for some multi-block
grids were found to be slow, slower than problems for which single-block grids could be
generated. The objective of this study is twofold. First, compare the convergence rate of
single- and multi-block grids on a consistent basis. Second, examine the effects of
blocking strategy on the convergence rate of multi-block grids.

Blocking strategies evaluated include location of block boundaries in relation to
flow features and thickness of the wrap-around block to resolve the viscous region next to

solid surfaces. Three ways were used to generate different multi-block grids: (1)

32

partitioning an initially single-block grid, (2) re-partitioning an initially multi-block grid,
and (3) overlapping different grids. Method 1 enables a direct comparison between the
convergence rates of single- and multi-block grids since the overall grid is identical.

Methods 2 and 3 enable a comparison of different blocking strategies.

<3 7..

(b)

C} (c)

Figure 3-1. (a) NLFO414 airfoil. (b) 623 ice shape. (0) 623 iced airfoil.

 

3.2 Problem Description, Formulation and CFD Flow Solver

The convergence rate of single- and multi-block grids is evaluated by obtaining
steady-state solutions to the natural laminar flow airfoil (NLFO414) with the 623-ice
shape (see Figure 3-1) at two angles of attack, one low (22°) and one near stall (52°).
The detailed geometry of this airfoil and ice shape can be found in Ref. [1].

This iced airfoil was selected because it is sufficiently complicated geometrically,
but a high-quality single-block grid can still be generated by using a technique presented
in Chapter 2. Having a geometry that has both single- and multi-block grids enables a

direct comparison between the two types of grids.

33

The flow conditions employed for all simulations are as follows. The
approaching Mach number is 0.29. The freestream temperature is —5° C. The Reynolds
number based on the freestream conditions and the chord length is 6.4 x 106.

The flow around the iced airfoil is model by the ensemble-averaged conservation
equations of mass (continuity), momentum (full compressible Navier-Stokes), and energy
for a thermally and calorically perfect gas. Turbulence is modeled by the two-equation
shear-stress transport (SST) model. Wall functions were not used. Integration of the
conservation equations and the SST model were to the wall, where no-slip and adiabatic
wall conditions were imposed.

Though a number of open-source codes are available, one of the most versatile
and well supported codes is the WIND codeml’m] In this code, the convective terms
were approximated by second—order upwind differencing formulas (physical Roe). Time
marching to steady state was accomplished by the Euler implicit method and ADI
approximate factorization. The default setting in WIND was used for all run parameters

in all simulations. The Courant number was set at 1.3.

3.3 Multi-Block Grids Evaluated

Three different types of multi-block grids were evaluated, generated by the
following three methods: (1) partitioning an initially single-block grid, (2) re-partitioning
an initially multi-block grid, and (3) overlapping different grids. All grids investigated

are summarized in Table 3-1.

34

Table 3-1. Summary of Multi-Block Grids Studied

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

Grid Type Description
Mufti-Block Grids from a
Single-Block Grid
SB single-block
S-MB-2 partitioned into multi-block
S-MB-3 partitioned in complex flow regions
S-MB-4 with “thin” wrap-around block
S-MB-S with “thick” wrap-around block
S-MB-6 with “medium” wrap-around block
Different Blocking of a
Multi-Block Grid
M-MB-l block boundaries parallel to flow
M-MB-2 block boundaries perpendicular to flow
M-MB-3 larger block size
M-MB-4 2-cell overlap at block boundary
M-MB-5 3-cell overlap at block boundary
Chimera Grids”
Chimera] big overlap
Chimera2 minimal overlap
'Number of grid points is 85,041 (816 x 101, 125 x 21); l-cell overlap at block
boundary

"Number of grid points: 114,469; l—cell overlap at block boundary unless
otherwise noted
#Number of grid points:100,851

3.3.1 Multi-Block Grids from a Single-Block Grid

Figure 3-2 shows a high-quality single-block grid for the NLF 0414 airfoil with
the 623 ice shape. By single-block, it meant that the grid around the iced airfoil is a
single-block grid. The grid shown in Figure 3-2 has 85,041 grid points with five grid

points within a y+ of 5 about most of the iced airfoil. This single-block grid was

35

generated in two steps. First, a grid was generated for the clean NFL0414 airfoil (i.e., the
airfoil without ice). Next, choose one set of grid lines near the airfoil surface, but a small
distance away from the iced airfoil, and regenerate the grid between the selected set of
grid lines and the iced airfoil surface to account for the ice shape. For further details on
how to generate high-quality single-block grids for complicated geometries, see Refs.
[1 l] and [16].

To study the effects of single- versus multi-block grids, the single block (SB) grid
shown in Figure 3-2 was partitioned into a number of different multi-block (MB) grids as
shown in Figure 3-3 and Figure 3-4 by cutting along different grid lines. In Figure 3-3,
the SB grid in Figure 3-2 was partitioned to form two MB grids, one with seven blocks
(S-MB-2) and another with nine blocks (S-MB-3). S-MB-2 has block boundaries in
regions where the flow patterns are not too complicated. S-MB-3 has block boundaries
in regions with more complicated flow features (e. g., boundaries cut through the
recirculation regions downstream of the two ice horns).

Since wrap-around grids are critical to the generation of high-quality multi-block
grids, Figure 3-4 shows three MB grids generated from the SB grid in Figure 3-2 with
different thicknesses of the wrap-around grid: thin (S-MB-4), thick (S-MB-S), and
medium (S-MB-6).

With all of the MB grids shown in Figure 3-3 and Figure 3-4 being identical to the
SB grid shown in Figure 3-2, the convergence rates of the MB grids can be compared

with the convergence rate of the SB grid on a consistent basis.

36

7
,7 1',
V! l I
.7 ’J-Ilh 1,57}!
. “Him/uh:
. ,7I'I,7Il,l',:',l'
. «' II' N 7' I
""r"7l’l"l
.u/,.I,.I 7'

‘ A'.’

.
~ \
. \‘ \\ \
5.. ,. .~ ~.~2~=
'.~:.~:.~ «:9 ‘
0.;.;:.~;_:g.§‘
-... .c,.~ r: a

I

[I
III

m” I"
”In"

I"

H.

[U I
"lull/"H'H/Ifl/
H",,,mll’
‘1‘ ”lull
I”! 7

mm

’1

I! ’/I

III Ill-
1, ,"r,Zl/‘ ‘

I! / I

1,704,”),

I ’1’, I,” I

"II“ /I I
I

t
I"!

I], I
l
1,0,.

MIMI. r
l I I.
..':""'7'"'I‘:""i
'r
’ l

 

 

 

 

 

 

 

 

 

 

 

S-MB-2 S-MB-3

Figure 3-3. MB grids from the SB grid shown in Figure 3-2. Solid lines indicate
block boundaries.

37

I
J/ //-/?//x / ,/’/l/- l// / / "L”:
\\ / // \‘\ P, J/ \\ / / / fix 7" :. /
‘7 ~=~ ' / WI’ \ .1 \r'\ / /
\- l / [’4 A\\ x/ T]! T l/ '/- /
( ~Gay/’1 E\ K |’\‘~~/"l (0 .kiv

i l ‘\ Va

3“ \,v \) ‘\,‘\

76/ ./ J'- l{ l/

I ll '1 (/
'7' 5 l

/ // Al 1 ///S\\
\ l l

’V‘m.‘q_. 5‘. i T‘ E i f' T“
\ -=\::u‘\:2:;==‘§ \" j}: a fl\ ~~‘~- i. J 4 ix“? “ —. r-\
k “\ \ «\wr‘fi’ \j “~ {A
\\ \

Figure 3-4. MB grids with wrap-around grids from the SB grid shown in Fig. 2.
Solid lines indicate block boundaries.

3.3.2 Different Partitions of a Multi-Block Grid

Figure 3-5 shows a high-quality MB grid generated for the NLFO414 airfoil with
623 ice. This grid has a “thick” wrap-around grid and transition layers. The total number
of grid points is 114,469. There are five grid points within a y+ of 5 over most of the
iced airfoil. This MB grid enables higher resolution of the horn regions of the ice shape
than does the SB grid shown in Figure 3-2.

To study the effects of blocking strategies, the MB grid shown in Figure 3-5,
referred to as M-MB-l, was repartitioned into two other MB grids, M-MB-2 and M-MB-
3, as shown in Figure 3-6. Note that all three MB grids shown in Figure 3-6 have a
“thick” wrap-around grid, which is necessary to form high-quality multi-block grids. In
M-MB-l, block boundaries start from the ice horns and follow the contour of the airfoil.
In M-MB-2, the block boundaries are perpendicular to the flow direction. M-MB-3 differ

from M-MB-2 by having blocks that are larger in direction normal to the airfoil.

38

 

’ / ’/_ ii:// / \\. / :1 A: ”" g ‘ \ ///:' ’ :1////
7:: :2 77; , 7 / «Lily/:5
(‘7‘ LT“ (.‘ ~,\ \ "~.“\7 ‘ "/1
\ l l i! ‘l k
‘7‘ .7 \ ‘7‘ 7 \, , N
1‘. 7 -- T _ 1 S r: K \_‘- 71;: V
‘5‘}- “~‘\\'\:\ ‘\ T ‘5: :g, ‘ / ‘* x- -1:
\ ‘\
\ \\_._,7/ f!
M-MB-1 M-MB-2 M-MB-3

Figure 3-6. Different partitions of the MB shown in Figure 3-5.
In all of the MB grids mentioned so far (Figure 3-3, Figure 3-4, and Figure 3-6),
there is a perfect overlap of one cell between all adjacent blocks to transfer information
between blocks. Since more overlap may improve convergence rate, the M-MB-l grid

was studied with one cell overlap (M-MB-l), two-cell overlap (M-MB-4), and three-cell
overlap (M-MB-S).
3.3.3 Overlapping Chimera Grids

Of the structured grids, the overlapping grid, also known as the Chimera grid, is
the easiest to generate. When applying Chimera grids to iced airfoils, there are four

steps. First, generates a grid for the clean airfoil without ice; this grid is referred to as the

39

main grid. Second, generate one or more grids for the ice shape without regard for the
airfoil. Third, combine the different grids generated at the correct physical locations.
Fourth, cut a hole in the main grid to account for the geometry of the ice shape.

Figure 3-7 shows two Chimera grids employed for the NLFO414 airfoil with the
623 ice. Chimera1 has a large overlapped region between the main grid and the grid

about the ice. Chimera2 has minimal overlap between the two grids.

" Mal/74721571241.
MWflWM’ ‘
Mavmw .

”7317771 .. ’
Chimera1 Chimera2

 

Figure 3-7. Chimera grids for the NLFO414 airfoil with 623 ice.

3.4 Results on Convergence Rate

As noted in the introduction part of this chapter, the objective of this study is to
compare the convergence rate of single- and multi-block structure grids and to examine
the effects of blocking strategy on the convergence rate of multi-block grids. The test
problem is the NLFO414 airfoil with the 623 ice at two angles of attack, or = 22° and a =
52°.

The convergence rate of the different single- and multi-block grids will be
represented by three parameters: (1) the number of iterations needed for the residual to

plateau (denote by N), (2) the value of the residual after it plateaus (denote by Rmin), and

40

(3) oscillatory behavior of the residual after it plateaus (denote by A for amplitude and P

for period of oscillation). The residual R is defined as follows:

R = _11_J|:2(Un+l _Un)2:|

131'
where I is the total number of grid points; J = 4 is the number of flow variables; i
is the index for the grid points; j is the index for the flow variables; U is the solution

n+1

vector; and U — Un is the change in the solution after one iteration.

3.4.1 Multi-Block Grids from a Single-Block Grid

Figure 3-8 shows the convergence history of the SB grid (Figure 3-2) and its MB
versions, S-MB-2 and S-MB-3 (Figure 3-3), for two angles of attack (or = 22° and 52").
From this figure, it can be seen that SB, S-MB-2, and S-MB-3 all require about the same
number of iterations to converge, N z 30,000 when or = 22° and z 50,000 when or =
5.20. When or = 52°, even Rm,n and the behavior of the convergence curve are identical
for SB, S-MB-2, and S-MB-3. When or = 22", there are a few differences. First, Rm,n is
slightly lower for the SB grid at 8 x 10'7 versus 2 x 10'7 for S-MB-2 and S-MB-3.
Second, SB grid has low amplitude oscillations, high-frequency, after R plateaus. But, S-
MB-2 and S-MB-3 have high amplitude, low frequency oscillations (A z 2 x 107; P =
2000 to 3000 iterations). The low-frequency, high amplitude oscillations for S-MB-2 and

S-MB-3 are due to wave reflections across some common block boundaries.

41

 

 

 

 

 

 

 

 

 

_3 ‘ —SB ----- S-MB2 - -*S-MB3
a 4 -
:L " r-
0 '5 L “'0‘ a A
9 V “9“
—6 8 '3‘
-7 L
'8 r r r r
0 1O 20 30 40 50
Iteration (103)
-2
_3 1 —SB S----MB2 -- S-MB3
(<74 7
cl
0 -5
Q d...
_r -6 7
-7 L
l
l
_8 L r r r r
o 10 20 30 40 50
Iteration (103)

Figure 3-8. Convergence history for SB, S-MB-2 and S-MB-3. (Top) or = 22°.
(Bottom) or = 52°

At this point, it is interesting to note that S-MB-2 and S-MB-3 have essentially
the same convergence history even though S-MB-3 has block boundaries that cut
perpendicularly through some complicated flow regions.

Figure 3-9 shows the convergence history of the SB grid (Figure 3-2) and MB
versions with thick and medium wrap-around grids, S-MB-S and S-MB-6 (Figure 3-4).

From this figure, it can be seen that if the wrap-around grid is sufficiently thick (S-MB-

42

5), then the number of iterations needed to converge and the minimum residual achieved
are about the same as those of the SB grid. For the S-MB-6 grid, which has a “medium”
wrap-around grid, Rm is two orders of magnitude higher than that for the SB grid, and
did not converge. For the S-MB-4 grid, which has a “thin” wrap-around grid, the

solution diverged. It is for this reason that results for S—MB-4 are not given.

 

— SB ----- S-MB5 -+- S-MB6

{*‘WENWMEMF‘flK

        
 
   

'

 

 

 

0 10 20 30 40 50
Iteration (103)

 

-— SB ----- S-MB5 -+- S-MBG

E‘Rfiqumwwwm

 

 

 

 

0 10 20 30 40 50
Iteration (1 03)

Figure 3-9. Convergence history for SB, S-MB-5, and S-MB-6. (Top) or = 22°.
(Buttom) a = 52°.

The results presented in this section show that MB grids can converge as quickly

as SB grids with nearly the same minimal residual. It appears that block boundaries

43

perpendicular to the flow direction are preferred over ones that are parallel to the flow
direction, especially in regions where the flow may be complicated. On wrap-around
grid, they must be sufficiently thick to converge quickly. It turns out that high-quality
multi-block only use “thick” wrap-around grids.

Figure 3-10 shows the predicted lift (CL) and drag (CD) coefficients for the SB
and the MB grids studied that generated converged results. From this figure, it can be
seen that SB and MB grids give similar results. Not shown is that the Mach number and

density contours are also similar for all of the SB and MB grids studied if they converge.

 

0.9

O 8 _ @ AOA=2.2 AOA=5.2

0.7 r

CL

0.6 r

0.5 -

 

 

     

 

 

III
1 II.

0.4 ~ ~ -
SB S-MB-2 S-MB-3 S-MB-4 S-MB—5 S-MB-6 Exp(30)

 

0.12

 

0-10 * EAOA=2.2 §AOA=52
0.08 —

0.06

CD

0.04
0.02

  

 

 

 

 

0.00

 

SB S-MB-2 S-MB-3 S-MB—4 S-MB-5 S-MB-6 Exp(3D)

Figure 3-10. Lift and drag coefficients for the SB, S-MB-2, S-MB-3, S-MB-5, and
experimentally measured ones.

44

The predicted lift and drag coefficients, however, do differ somewhat from the
experimentally measured ones reported in Ref. [1]. This discrepancy is due to the
turbulence model used, and not the grid. Chung & Addy [4] showed that the shear-stress-
transport (SST) turbulence model used in this study provided the best prediction on lift

and drag for clean airfoils, but the worst results for iced airfoils.

3.4.2 Different Partitions of a Multi-Block Grid

Figure 3-11 shows the convergence history of the Multi-Block grids (Figure 3-6).
From this figure, it can be seen that M-MB-2 has the best convergence rate of the three
MB grids. When or = 22°, the number of iterations needed to converge, N, is more than
80,000 for M-MB—l and a little over 70,000 for M-MB-Z. When or = 52°, N for M-MB-l
and M-MB-Z are about the same, around 50,000. But, Rmin for M-MB-2 is two orders of
magnitude smaller than that for M-MB-l. M-MB-3 also has a Rm,“ that is as low as that
for M-MB-2, but N is appreciably larger than 50,000. M-MB-2 converges faster than M-
MB-3 because for M-MB-2, the blocks around the ice extended only a short distance
away from the iced airfoil and are surrounded by a SB grid. That SB increased
communication among the blocks.

From the convergence history results for M-MB-4 and M-MB-S, it can be seen
that When or = 22", increased number of cells that are overlapped in the adjacent blocks
— 2, or 3 — did not appear to improve convergence rate; when or = 52°, it did improved
the convergence rate, but its Rm... is two orders of magnitude higher than that of M-MB-2.

Figure 3-12 shows the predicted lift (CL) and drag (CD) coefficients for M-MB-l
to M-MB-S along with the experimentally measured ones. Unlike the S-MB-# grid

results shown in Figure 3-10, there is greater variation in the predicted lift coefficient

45

from the different M-MB-# grids.

Not shown, however, is that the predicted Mach

number and density contours are similar for all of the M-MB-# grids studied. There are

only minor differences in the separated region induced by the lower ice horn.

Similar to SB cases, the predicted lift and drag coefficients differ from the

experimentally measured ones reported in Ref. 1. As noted, the discrepancy is due to the

turbulence model used, and not the grid.

 

     

 

 

 

 

 

 

 

-2
——- M-MB-1 -+- M-MB-2 - - - M-MB-3
-3 j . . ----- M—MB-4 -+- M-MB-5
g '4 " 9 g . ’ YWQW ‘ . '''''''''
o -5 ' «."W. """"""" . .
O ’ a". (- 5W. a ' ~v7vv
.—l '6 _ . MEN," Wm;
-7 e
'8 r r r J r r r
0 10 20 30 40 50 60 70 80
Iteration (103)
-2
-— M-MB-1 --°-- M-MB-2 - - - M-MB-3
-3 . ----- M-MB-4 -+- M-MB-5
A4 _ l
g 5 ....”9..":q MIN ’-
- a....::. 0,7.
8 .....:..:Q;:;\\V\flflfipfi
-6 _ "0'03“... My“, .7
“flow". w fotjl ‘l 7 7‘
-7 7
_8 r r r u r r r
0 10 20 30 40 50 60 70 80
Iteration (103)

Figure 3—1 1. Convergence history for Multi-Block grid. (Top) or = 22°. (Buttom) or

= 52°.

46

0.9

 

E AOA=2.2 Q AOA=5.2

CL

 

 

 

0.12

 

0.10 -
0.08
0.06 ~

Co

0.04
0.02

 

 

 

0.00

M-MB-1 M-MB-2 M-MB—3 M-MB-4 M-MB—5 Exp(3D)

Figure 3-12. Lift and drag coefficients for the Multi-Block grids and experimentally
measured ones.

3.4.3 Overlapping Chimera Grids

Results obtained for the Chimera grids show that Chimera] and Chimera2 grids
give similar results and convergence history. Figure 3-13 compares the convergence
history of the Chimera grids with the SB grid and M-MB-l grids. This figure shows the
Chimera grid to converge almost as quickly as single-block grids with even smaller
minimum residual. Figure 3-14 compares the predicted lift (CL) and drag (CD)
coefficients for SB, M-MB-l, M-MB-Z, Chimera1, and Chimera2. This figure shows
that the predicted lifi and drag coefficients are similar to those predicted by using the

single- and multi-blocks.

47

 

— SB —~— M-MB—1 ----- Chimera1 -— Chimera2

 

 

 

Iteration (103)

Figure 3-13. Convergence history of SB, M-MB-1, and Chimera grids for a = 22°.

 

0.9
B AOA=2.2 AOA=5.2

 

 

 

 

 

 

 

 

 

Chimera1

Figure 3-14. Lift and drag coefficients for the SB, M-MB-1, M-MB-2 , Chimera1,
and Chimera2.

48

3.5 Summary

This study showed that when a single-block grid is partitioned into a multi-block
grid with block boundaries perpendicular to the flow direction, the number of iterations
needed to converge is about the same as that for the single block. The minimum residual,
however, is slightly higher and can oscillate with higher amplitude. When a single-block
grid is partitioned into a multi-block grid with block boundaries parallel to the flow
direction (e. g., wrap-around grids), convergence rate can deteriorate if the wrap-around
grid is too thin. If the wrap-around grid is sufficiently thick, then the convergence rate
can be nearly as good as a single-block grid.

The adverse effect of having block boundaries parallel to the flow direction in
region with flow complexities was confirmed with different partitions of a multi-block
grid. 2 or 3 overlap cells in adjacent blocks may increase convergence rate than using
just 1 overlap cells. The Chimera grids were found to converge as fast as single block

grids with about the same minimum residual.

49

Chapter 4. Computing Aerodynamic Performance of 2-D
Iced Airfoils with Structured Grids

The ice accrued on airfoils can be very complicated geometrically with horns,
feathers, and rough surfaces. Even in 2-D, it is difficult to generate high-quality
structured grids. This chapter examines promising methods and guidelines developed for
the generation of single- and multi-block grids about 2-D iced airfoils.

This examination is based on the commercial transport airfoil with the 145m ice —
one of the most complicated ice shape that can form -— at two high angles of attack, 6° and
10°, with 10° being very close to stall. The flow was modeled by the ensemble-averaged
compressible Navier-Stokes equations, closed by two different turbulence models: the
shear-stress-transport turbulence model and the Spalart-Allmaras model with integration
to the wall. All solutions were generated by using the NPARC WIND code.

Results obtained show grid quality to have significant effects on convergence rate.
With a poor quality grid, the solution could “blow up” or not converge with large
oscillations in the residual. Results obtained also show that solutions obtained on multi-
block grids may differ from those obtained on single blocks under certain flow
conditions, indicating the importance of correct data transfer across block boundary

during the solution procedure.

4.1 Introduction

For iced airfoils — even two-dimensional (2-D) ones — the generation of high-
quality structured grids is a major challenge. In Chapter 2, we presented a number of

methods to generate high-quality single- and multi-block structured grids for complicated

50

2-D ice shapes. We applied this method on grid generation for the NLF 0414 airfoil with
the 623 ice shapem To generate high-quality multi-block grids, we showed that a “thick
wrap-around” grid is needed to ensure that grid points clustered next to solid surfaces do
not propagate into the interior of the flow domain. To minimize and confine the adverse
effects of rough and jagged surfaces on the quality of the grid, we suggested using a
transition layer next to solid surfaces. We also developed a method for generating single-
block grids for iced airfoil based on an idea presented by Tai [1”.

Since multi-block grids can converge slower, We studied the effects of blocking
strategy on the convergence rate to steady-state. We showed that when a single-block
grid is partitioned into a multi-block grid with block boundaries perpendicular to the flow
direction, the number of iterations needed to converge is about the same as that for the
single block. When that same single—block grid is partitioned into a multi-block grid with
block boundaries parallel to the flow direction (e.g., wrap-around grids), convergence
rate can deteriorate if the wrap-around grid is too thin. If the wrap-around grid is
sufficiently thick, then the convergence rate can be nearly as good as single-block grids.

So far, we have only demonstrated their grid generation techniques and blocking
strategies to a moderately complicated ice shape (the 623 ice), one with two large but
relatively widely separated horns. This iced airfoil has also been studied by others using
simpler methods (e.g., Chung, et aim). Thus, the applicability of this method and
conclusions for more complicated ice shapes has not been demonstrated. The objective
of this study is to assess the methods and conclusions of Chi, et al. “6] and Zhu, et al. [17]
by applying them to a much more complicated ice shape that has not been simulated

before because of its complexity (e.g., one with multiple and highly extended horns in

51

close proximity to each other) and to a more complicated flow condition (e. g., one with
large separated regions).

The organization of the remainder of this chapter is as follows. Section 4.2
describes the iced-airfoil problem selected for investigation, and the challenges that it
presents. Section 4.3 summarizes the formulation and the numerical method of solution.
Section 4.4 describes the grids generated and the methods employed. Section 4.5

presents the results obtained by using the grids generated.

4.2 Problem Description

The iced airfoil selected to carry out the objectives of this study is the
commercial—transport airfoil with the 145m ice shape (see Ref. 1 for details of the
geometry). The airfoil and the ice shape about the airfoil’s leading edge are shown in
Figure 4-1. The 145m ice shape is one of the most complicated ice shapes that can form
about an airfoil. This ice shape is considered complicated because it has multiple horns
that are relatively thin, highly extended, and very close to each other. Also, it has the
usual complexities such as feathers and surface roughness and jaggedness.

In this study, the flow conditions employed for all simulations are as follows. The
approaching Mach number is 0.29. The free stream temperature is -11.6° C. The
Reynolds number based on the free stream conditions and the chord length is 6.4 x 10°.

Two angles of attack (a) were simulated, 6 and 10 degrees.

52

 

Figure 4-1. (Top) Commercial transport airfoil; (Buttom) 145m ice shape near
airfoil’s leading edge.

The selected iced airfoil, owing to its geometric complexity, poses a number of
challenges in the generation of structured grids. The first is that it may not be possible to
generate a high-quality single-block grid because of the deep concave regions between
horns. The second is that with highly extended horns being so close to each other, severe
constraints are imposed on the thickness of the wrap-around grid. This constraint is a
concern because we has shown in Chapter 3 that if the wrap-around grid is too thin, then
the solution will either not converge or converge extremely slowly with large oscillations
in the residual as if the solution is unsteady.

The ultimate measure of a high-quality grid is the ability to generate accurate
solutions. For solutions with steady states, another measure is the ability of the grid to
enable the solution algorithm to converge to the steady state solutions in as few iterations
as possible with residuals that plateau near the round off errors. Clearly, solutions that

describe more complicated flow features will be a more stringent test on convergence

53

rate. It is for this reason that the angles of attack selected are quite high. With high
angles of attack, the separated regions will be quite large, extending across several
different blocks of a multi-block grid. Also, unsteadiness effects may occur. The authors

know of no one who has reported a solution for the 145m ice shape yet.

4.3 Formulation and Numerical Method of Solution

The flow past the iced airfoil described in the previous section is modeled by the
ensemble-averaged conservation equations of mass (continuity), momentum (full
compressible Navier-Stokes), and energy for a thermally and calorically perfect gas.
Two turbulence models were used. One is the two-equation shear-stress transport (SST)
model.[27]’[28] The SST model uses the k-w model in the region next to the wall and the k-
2 model in the region away from walls. The other turbulence model used is the one-
equation Spalart-Allmarass (S-A) model. Chung, et all“ have shown the SST model to
perform better for clean airfoils (i.e., airfoils without ice buildup) and the S-A model to
perform better for iced airfoils. For both models, wall functions were not used.
Integration of the conservation equations and the turbulence models were to the wall,
where no-slip and adiabatic wall conditions were imposed.

Though a number of open-source codes are available, one of the most versatile
and well supported code is the WIND code.“2]’“3] Several options are available on
solution algorithm. In this study, the convective terms were approximated by second-
order Roe upwind differencing formulas. When only the steady-state solution is of
interest, time marching to steady state was accomplished by an implicit method based on
ADI-type approximate factorization with local time stepping. When unsteady solutions

are of interest, local time stepping is not used, and sub-iterations in pseudo-time are used

54

to ensure that time-accurate solutions are obtained. Relatively few investigators have

reported unsteady solutions from WIND. Hamed, et a1.“41

presented a study of an
unsteady problem. Their approach was adopted here, which ensures second-order

accuracy in time.

Table 4-1. Summary of Grids Generated

 

Grid Name Description

 

SB single-block grid around ice and airfoil
(Figure 4-2) (995 x 97 + 125 x 23 = 99,390)

 

multi-block grid around ice and airfoil

“(3?1'23) (995 x 31+ 385 x 35 +12 x 33 +14 x 42+ 341x 35 + 781x 51 +125
gm x23 =99,945)

 

 

 

 

MBl-Sl same as MBl-SO except smoothed grid but kept the wrap-around grid
(Figure 4-4) unchanged
MBl-SZ same as MBl-SO except smoothed all grid and allowed a variable
(Figure 4-5) thickness wrap-around grid (102,195 grid points)
MBZ'SZ same as MBl-SZ except partitioned a few more blocks to test ideas
(F 1gure 4-6)

same as MBl-SO except use the blocking strategy shown in Figure
MBZ-SO 4-6

 

 

same as MBl-Sl except use the blocking strategy shown in Figure
MB2-Sl 4-6

 

 

 

4.4 Results on Grids Generated

Seven grids were generated for the iced airfoil shown in Figure 4-1. The grids
generated are summarized in Table 4-1. The reasons for generating these seven grids will
be made clear in the results section on solutions generated. All seven grids were

[91,001

generated by using transfinite interpolation with varying degrees of local elliptic

smoothing.

55

The first of the seven grids, referred to as SB in Table 4-1, is a single-block grid
(Figure 4-2). This grid was generated by using the approach described by Chi, et 3106]
However, considerable amounts of local elliptic smoothing were needed. The generation
of this grid demonstrates that the single-block grid-generation method presented by Chi,

et a1. “6]

can be applied to highly complicated ice shapes.

The second grid, referred to as MBl-SO in Table 4-1, is a multi-block grid (Figure
4-3). It is generated by using the approach described by Chi, et a]. [16] in which no elliptic
smoothing was used. This grid has a transition-layer grid that is merged with a wrap-
around grid. The thickness of the wrap-around grid was nearly constant to ensure that the
first few grid spacings away from solid surfaces are the same around the entire airfoil.
Since the wrap-around grid has essentially the same thickness through out, the maximum
thickness permitted is dictated by the minimum distance between adjacent horns. It is for
this reason that the wrap—around grid for MBl-SO is relatively thin.

The third grid, referred to as MBl-Sl in Table 4-1, is also a multi-block grid
(Figure 4-4). It is identical to MBl-SO in that it also has a wrap-around grid of nearly
constant thickness all around the airfoil. It is different from MBl-SO in being smoother.
Basically, elliptic smoothing was used locally at several locations to generate the
smoothest grid possible for a multi-block grid with a wrap-around grid of constant
thickness.

The fourth grid, referred to as MBl-SZ in Table 4-1, is another multi-block grid
(Figure 4-5). It is identical to MBl-Sl in that it has considerable elliptic smoothing. It is

different from MBl-Sl in that the thickness of the wrap-around grid is no longer kept

constant around the airfoil, removing an unnecessary constraint. By allowing for variable

56

thickness, a much more flexible high-quality multi—block grid can be generated. In
particular, the thickness of the wrap-around grid can now be thick enough to ensure a
converged steady-state solution. Also, grid points can be used more efficiently. For
example, the grid spacing next to solid surfaces needed to resolve recirculating flow
between horns can be much larger than those needed to resolve high-speed boundary-
layer flows on the airfoil surface away from separated regions.

The fifth grid, referred to as MB2-$2 in Table 4-1, is still another multi-block grid
(Figure 4—6). MB2-82 is identical to MBl-SZ in the location of every grid point. The
only difference is that blocking strategy is different. With MBI, the block boundaries
next to the horns are nearly parallel to the main flow direction about the airfoil. With
MB2, the block boundaries are nearly perpendicular to the main flow direction. In
Chapter 3, we showed MB2 type blocking to converge faster than the M31 type
blocking.

The sixth and seventh grids, referred to as MB2-SO and MB2-SI, are identical to
MBl-SO and MBl-Sl, respectively in the location of the grid points. They differ only in
the blocking.

For all grids generated, the first grid point away from all solid surfaces have y+
values that are less than unity. When the SST model is used, there are at least five grids

points within a y+ of five.

57

\
\s“\\‘“‘\\“\t\\m
W“\\\“‘{\\“‘§3\i
Mss‘
fefige‘fififi
e ~

\
cs‘s‘
§

\\
~‘

I . 1| It
“iii" 111311,}:
IIIIIIIIiIgliilIIIIII.I.
Illlllll I

 

  

II‘II‘II‘fi‘I ‘I
ll
lllllllt'tlill‘l‘thh
It“ (“Irma
‘ I

“ll ‘
I'l'lumullllwt‘l‘l‘lé‘w

lllllllll “W“
ll

I
“m“ H II
ll“ |II
““Illjilllj'l‘lll‘l“ lIlIII

    
    

  

 

IIIIIIIIl‘I I
I1 i
IIIIIIIII‘I-..I‘Iut.u

I“
“Minn.“

       
 
 

     

 

      
      
   
   

I

F'
Igure 4-3. Multi-block 12MB1-SO

58

\ .
Make“

‘9‘:

\ ‘\ ‘\\

‘\“\ \
s‘.

\ \\ \\\\ \\\“
WWII:

/,,/
I
I
II

’I
’ :IlgI/l
I
520/4 I

III]

I I

Ila, I
I I

 

   
    
  
  
     
   
  
   

 
 
   
   

[Ill/77107”
,I' I We;
IIIIIIIIII’llllI'I'tIq

 
   

I .C Q
‘9‘
,w
3"

llllllllII "tint...
IIIIIIIIIII‘I‘illl‘IIIIIItt
lll\\\\\\\\
\\\\\\\\\ tit:

llll
\\ t
\\ \\ \\\\°‘
Willi)“

      
 

 

Figure 4-5. Multi-block 3: MB1-SZ.

59

 
    
 

    

 
 

     
   
 
 

      
  
   
  

   
   

  

  

   

§\§\ § \\ \\ .
§§\ §§§§W‘ \\\\\
\ '\ \ \ \ h \‘i \~\\\\\\\\ 111111
$§§§§$§§iiiWiim‘lliiiiiliiiiii‘llllliiit‘illit‘t‘iI‘t‘tt I
\ \V i \\\ 1“ m
§§§§§§§wt
\ \ § § § \\ ‘ \ ‘\ \\ \\\~ “\N \W‘“ \uI\|“““.mm I
\ § § \ § § \\ S \\ ‘2‘ “ ‘0‘“ \‘t‘ “on“ . .
§§$§s NW \ -\“‘w\“\ 3:“ was ‘
\\ \ §%se t «. «m,
N \ \ \ \ x x \
\ N \ N § % §$$
§ § § esss ss
‘- : § §§ e e s
:'¢$§ :s * *
h§§ §S¢$s
- his tzg§§§
; :

\\\\\\\\\\‘
\\\

\\\\\\\\\\\\\\\\\\\\\
\\\\\\\iI \

\‘~\\“‘\\\\\\\‘\ “mm \\
\\\\‘\“s&\\“}:\\\u
. \

\\\\‘

\‘N\\\\\\\ ‘
\\\\\\\\\\\\\\\\

\\\\\\\\\
\\\\\
“‘3: \\\\““

\\\\\\\\\\\\\\\\\\\\\
\\\\
at

w

W
\\\

\
\\
\
\\\\\\\\\\\\\\\\
®\\\\\\\\\\\\\\\

\\\\\\\\\\\

\\\\\‘

M
\\\\\
«s

W
)\\
\\\

\\\\\\\\

s

I '1];
/ 1
I7 2“},
I I 715'13'0?

W

\\\

s“
s

\\\\\\‘\\‘\\
ssss

\\\\\

\\\\\\\\

// I I, III"! I 0,;
/ ///"4/ " Ill/I’ll’ll’llj'uzr.
/ , ’ ”W ’

s

W
¢

/¢///

\\

// //
Wr’gmm
”ll/7”” [’IIIII’M'
I ”I, It”, “In '1 It '
”II "'Ittl’"tl’tstijéjltr.
m

Figure
4-
6. Multi—bIOCk 4. M

60

Table 4-2. Summary of Cases

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

(121a: (1 Grid T323306 A3312” Co/At Converged?
1 6 SB S-A 2nd 5 $38
2 6 MBl-SO S-A 2nd 1.3 no
3 6 MBl-SO ssr 2'"r 1.3 no
4 6 MBl-SO S-A 1st 1.3 yes
5 6 MBl-SO S-A 1st then 2nd 1.3 no
6 6 MBl-SO ssr 1st 1.3 yes
7 6 MBl-SO ssr 1st then 2“d 1.3 no
8 6 MB2-SO S-A 2nd 1.3 no
9 6 MBl-Sl S-A 2nd 1.3 no
10 6 MBZ-Sl S-A 2nd 1.3 no
1 1 6 MBl-S2 S-A 2"(1 1.3 yes
12 6 MB2-32 S-A 2"r 1.3 yes
13 6 MB2-82 ssr 2"d 1.3 no
14 6 MB2-82 S-A 2"d 5 yes
15 6 MBZ-SZ S-A 2nd 5 x 10'“ 8 yes
16 10 MB2-82 S-A 2'“F 1.3 yes
17 10 MBZ-SZ S-A 2nd 1 x 10'6 5 yes
18 10 SB S-A 2nd 5 yes

 

 

 

 

 

 

 

 

*

a = angle of attack in degrees. For steady-state computations, Co = CF L number
used in local time stepping. For unsteady computations, At is the time-step size.

Table 4-3. Predicted Drag and Lift Coefficients

 

 

 

 

 

 

 

 

 

 

 

Case No. Drag Coefficient Lift Coefficient
1 0.0484 0.540
4 0.0827 0.625
6 0.0816 0.620
11 0.0492 0.535
12 0.0460 0.581
14 0.0459 0.581
15 0.0461 0.583
16 0.1172 0.830
17 0.1172 0.830
18 0.1431 0.583

 

 

 

 

 

61

4.5 Results on Solution Generated

By using the seven grids described in the previous section, 18 simulations were
performed by using the WIND code to examine (1) how grid quality affects convergence
rate to steady state (assuming a steady—state solution exists), (2) the difference that can
result from using single and multi-block grids in the WIND code, (3) the effects of
turbulence model on convergence rate, and (4) unsteadiness effects that may be present
for the two high angles of attack examined. Table 4-2 surmnarizes all simulations made,
and Table 4-3 summarizes the predicted lift and drag coefficients. Figure 9 shows the
convergence history on the second norm of the residual, the lift coefficient, and the drag
coefficient for all cases that converged. Figures 10 and 11 show the predicted Mach
number contours and streamlines for selected cases.

From Table 4-2, it can be seen that with the single-block grid, SB, a converged
solution could be generated for both angles of attack (Cases 1 and 18). As shown in
Figure 4-7, the convergence rate was fast and with essentially no oscillations in the
residual as it plateaus. The results fi'om the single-block grid are used to assess the
results obtained by the multi-block grids. This is because the single-block grid is of high
quality and the solution from it is grid independent.

When the multi-block grids, MBl-SO and MB2-SO, a converged solution cannot
be obtained unless the order of Roe differencing is dropped from second order accuracy
to first order accuracy. By dropping to first-order accuracy, which has higher numerical
diffusion, a convergent solution was obtained for both the S-A and the SST turbulence
models (Cases 4 and 6). Though a converged solution was obtained, the lift and drag

coefficients predicted are much higher than those predicted by the single-block (Table

62

4-3). Also, the separated region predicted is much less than that predicted by the single-
block grid (Figure 4-8). Since Cases 4 and 6 gave almost identical results, only Case 4
was plotted in Figure 4-8.

To try to get a converged second-order accurate solution with a multi-block grid
that retains a nearly constant thickness wrap-around grid, MBl-Sl and MB2-SI were
generated. MBl-Sl and MB2-$1 are higher quality grids than MBl-SO and MB2-SO by
being smoother and less skewed. Despite the improved grid quality, a converged solution
cannot be obtained. One reason for this is that the wrap-around grid is too thin. As we
noted in Chapter 3, the thickness of the wrap-around grid strongly affect the convergence
rate.

For the 145 m ice shape, the only way to increase the thickness of the wrap-
around grid is to allow for variable thickness. Multi-block grids MBl-SZ and MB2-82
were generated with variable thickness wrap-around grid. With these two grids,
converged solutions were obtained with S-A model (Cases 11 and 12), but not the SST
model (Case 13). This indicates that SST is more sensitive to the quality of the grid. Of
the two blockings used, MBl-SZ converged faster than MB2-82 (Figure 4-7). Though
MBl—SZ and MBZ-SZ are identical in grid location and only differ in the blocking, there
were an 8% differences in the predicted lift coefficients and 7% difference in the
predicted drag coefficients (see Table 3). Since MBl-SZ grid gave lift and drag
coefficients that were closer to those obtained by SB than did MB2-82, MBl-SZ has a
better blocking than MB2-82. Though the lift and drag coefficients predicted by MBl-SZ
were fairly close to those predicted by SB, there are still some differences in the size of

the separated region just downstream of the ice (Figure 4-8).

63

Cases were also run to compare the convergence rate between the SB and MB2-
S2 grids with a CFL number of 5 (Cases 1 and 14). Figure 4-7 shows that the
convergence rates between the single- and multi-block grids are comparable. By using a
CFL condition of 5 instead of 1.3, the convergence rate can be accelerated by a factor of
about 10 and still obtain the same converged solutions.

Though MBl-SZ and MB2-82 converged to a steady state solution, a simulation
was performed in time-accurate mode (Case 15). The results obtained in this fashion
were essentially the same as the one obtained by using local time-stepping (Cases 12 and
14).

The above discussion was for the 6° angle of attack. Computations were also
performed for the 10° angle of attack (Cases 16, 17, and 18). Results obtained showed
even at this high angle of attack, the flow field is steady since both time-accurate and
non-time-accurate solutions yielded the same results. One disturbing observation is that
the size of the separated region predicted by the multi-block grid differed considerably
from the one predicted by the single-block grid (Figure 4-9). This may be due to how
boundary conditions at the block boundaries are implemented. This anomaly is still

being investigated.

64

— Casel — Case4 Case6 -— Casell — CaselZ -- Casel4 — CaselS

 

(L2)

 

 

 

 

 

 

 

 

 

 

0.20

0.15 I .
0.10 I. .....................................................
0.05 ‘1”, . . ____________
0.00 Al ' .......................................................
-0.05 I _______________________________________________________ l
-o,1o '— ............................................................

.0015 i l T l I
0 10 20 30 40

Iterations (10$)

 

CD

 

 

Figure 4-7. Convergence history for the residual of CL and CD. See Table 4-1 for
definition of case numbers.

65

 

Case 11: MBl-SZ, 2nd-order Case 12: MB2-$2, 2nd—order

Figure 4-8. Mach number with streamlines for a = 6 degrees.

   

Case 16: MB2-$2, 2nd-order Case 18: single-block grid, 2nd-order

Figure 4-9. Mach number with streamlines for a = 10 degrees.

 

66

4.6 Summary

The key findings of this work are as follows: (1) If grid quality is poor (e.g.,
MBl-SO, MBl-Sl), then the solution either blows up or appears as if it is unsteady. If
the grid quality is sufficiently good, then steady-state solution were able to be generated
with low residuals (second norm less than 105). (2) For both angles of attack
investigated (6° and 10°), there were no unsteadiness effects since time-accurate and non-
time-accurate runs yielded identical results. (3) Were able to generate a high quality
single-block grid (SB) by using the method of as shown in Chapter 2 But, selective
elliptic smoothing is needed. (4) Were able to generate high-quality multi-block grids
(MBl-SZ, MB2-82). But, variable thickness wrap-around grids must be used. (5) SB,
MBl-SZ, and MB2-82 grids gave similar results on convergence rate, lift and drag
coefficients, and flow features if the angle of attack is 6°, where separation bubble is not
so big. (6) SB and MB2-82 gave very different results if the angle of attack is 10°. The

reason for this is still being investigated.

67

Chapter 5. CFD Analysis of the Aerodynamics of a
Business-Jet Airfoil with Leading-Edge lce Accretion

For rime ice — where the ice buildup has rough and jagged surfaces but no
protruding horns - this study shows two-dimensional CFD analysis based on the one-
equation Spalart-Almaras (S-A) turbulence model to predict accurately the lift, drag, and
pressure coefficients up to near the stall angle. For glaze ice — where the ice buildup has
two or more protruding horns near the airfoil’s leading edge — CFD predictions were
much less satisfactory because of the large separated region produced by the horns even
at zero angle of attack. This CFD study, based on the WIND and the Fluent codes,
assesses the following turbulence models by comparing predictions with available
experimental data: S-A, standard k-s, shear-stress transport, vz-f, and differential

Reynolds stress.

5.1 Introduction

To date, most CFD studies on 2-D iced airfoils have focused on the following:
grid generation methods, grid resolution, turbulence models, and lift and drag coefficients
as a function of the angles of attack [3]-[7],[l6]-[l8]. Though much has been learned
from these studies, accuracy of CFD predictions is still unclear. For example, what are
the error bounds and confidence level in the computed lift and drag coefficients as a
function of angle of attack? In addition, relatively little emphasis has been made on
understanding the flow field induced by the ice accrued on the airfoil and how well they
are predicted. One reason is that very little experimental data are available that can be

used to assess the accuracy of such CFD predictions. Recently, experimental data for the

68

x-component velocity have been made available for a business-jet airfoil (GLC305) with
a 944 glaze ice shape and is reported in Ref. [20]. The availability of these experimental
data enables a more thorough interrogation of the CFD results generated. That is, in
addition to lift, drag, and pressure coefficients, the detailed flow field can also be
examined with confidence to provide a better understanding of how ice shapes affect

aerodynamics.

5.2 Objective and Approach

With the above backdrop, the main objective of this study was to understand how
well CFD can predict lift, drag, surface pressure, and the velocity field as a function of
the angle of attack for 2-D iced airfoils. The accuracy of CFD predictions was assessed
by comparing computed results with experimentally measured lift, drag, surface pressure,
and velocity field.

Since it is now possible to generate hi gh-quality grids for geometrically
complicated 2-D iced airfoils through the work described in Refs. [l6],[17] and [18], this
study on accurate CFD predictions focuses on the effects of turbulence models, the other
major source of error. Of particular interest is how well state-of-the-art turbulence models
can predict the aerodynamics induced by glaze and rime ice shapes. Glaze and rime ice
shapes, formed under different icing conditions, constitute the two fundamental ice
shapes. Though both ice surfaces are rough and jagged, glaze ice also has horns but rime
ice does not. Thus, these two ice shapes produce very different flow fields, which may
have different requirements on turbulence modeling. For rime ice, which produces only
very small separated regions at all angles of attack except near stall, a simpler turbulence

model may be adequate. For glaze ice, which can produce large separated regions

69

downstream of the horns even at zero angle of attack, one— or two-equations models may
be inadequate. Differential Reynolds stress models that can account for each of the
Reynolds stresses individually may be needed to predict the effects of streamline
curvature and the time-lagged response of turbulence to changes in the mean flow.

In order to examine as many turbulence models as possible, two codes were used
in this study. One is the widely used, open-source code, referred to as WIND, which was
used in the work reported in Refs. [2],[3],[4],[6][l6],[l7], and [18]. The other is Fluent, a
popular commercial CFD code with many turbulence models. WIND and Fluent codes

are described in more detail later.

5.3 Problem description

The airfoil, the ice shapes, and the flow conditions selected for study are those for
which the flow field is sufficiently complicated and for which there are experimental data
that can be used to assess the accuracy of the CFD predictions. The airfoil selected is the
business—jet airfoil (GLC305). The glaze ice selected is the 944 ice shape with two large
protruding horns. The rime ice selected is the 212 ice shape, which has considerable
surface jaggedness but no protruding horns. The airfoil and the ice shapes about the
airfoil’s leading edge are shown in Figure 5-1 and Figure 5-2. Details of the geometry are

given in Ref. [1].

7O

_

 

Figure 5-1.GLC 305 airfoil with 944 ice shape

I I I - ' ‘-

a

Figure 5-2. GLC 305 airfoil with 212 ice shape.

In this study, the freestream Mach number (Ma) was 0.12. Two freestream static
pressures (P = 20.5 psi and 37.0 psi) and two Reynolds numbers based on the freestream
conditions and the chord length (Re = 3.5 x 106 and 6.0 x 106) were investigated. The
angles of attack (AOA) simulated were 0, 4, 6, 7, 8, 9 for the 944 ice and 0, 4, 6, 7, 8, 9,
10, ll, 12 for the 212 ice. The reason for simulating so many angles of attack was to

obtain the details on lift and drag coefficients about the stall angle.

71

5.4 Formulation and Numerical Method of Solution

Two different CFD codes were used to generate solutions for the iced-airfoil
problems described in the previous section. One is a widely used, open-source code,
known as WIND [121.121]. The other is a popular commercial code, Fluent [22]. These codes
were selected because they are highly versatile and contain a wide range of turbulence
models.

For both codes, the flow past the GLC305 airfoil with the 944 and 212 ice is
modeled by the ensemble-averaged conservation equations of mass (continuity),
momentum (full compressible Navier-Stokes), and energy for a thermally and calorically
perfect gas. For most simulations, the one-equation Spalart-Allmaras (S-A) model [24] is
used to mimic the effects of turbulence. For iced airfoils, Chuang & Addy [4] showed the
S-A model to out-perform two-equation turbulence models including the highly regarded
shear-stress transport (SST) model [271428]. For Fluent, the following turbulence models
were investigated: S-A [24], SST, standard k-e [29], Durbin’s v2-f model [32], and a
differential Reynolds stress model [33134]. With Fluent, the near-wall region is always
modeled by the two-layer model of Chen and Patel [3 5 1. In all simulations with both
WIND and Fluent, the conservation equations and the turbulence models were integrated
to the wall, where the no-slip and adiabatic wall conditions were imposed (i.e., wall
functions were not used).

The numerical methods of solution used are as follows. For WIND, the
convective terms were approximated by second-order Roe upwind differencing scheme.
Since only steady-state solutions are of interest, time marching to steady state was

accomplished by an implicit method based on ADI-type approximate factorization with

72

local time stepping. For Fluent, which uses a finite-volume method, fluxes at the cell
faces are interpolated by using second-order upwind differencing. The SIMPLE
algorithm was used to generate steady-state solutions. For this segregated solver, the
convergence criteria used are to ensure that normalized residual is less than 10'6 for the
energy equation and less than 10'3 for all other equations.

To ensure proper comparison between the codes and among the turbulence
models, both WIND and Fluent used essentially the same grid systems as explained
below. For WIND, all grid systems generated consisted of two overlapping single-block
grids. One was a fine grid next to the airfoil, extending 0.6 chord length from the airfoil
in all directions (referred to as the inner grid). The other was a coarser grid (125 x 21)
that overlapped the fine grid by 0.1 chord length and extended 15 chord lengths away
from the airfoil in all directions (referred to as outer grid). The inner grid was the single-
block grid of interest. While generating this grid, grid lines were clustered next to the
airfoil surface so that the first grid point was within a y+ of unity. Along the airfoil
surface, equal arc-length was employed to create a grid as smooth as possible.

Since Fluent does not accept overlapping grids, the grids used by WIND and
Fluent are not exactly the same. The inner grid used by WIND is also used in Fluent.
The outer grid used, however, does differ from the one used by WIND. In Fluent,
quadrilateral cells with size 940*40 were generated in the outer region so the final grid
number has increased compared with the grid number used for WIND code. Finally, all
grids are merged together with total 132,500 nodes and 131,600 quadrilateral cells. Since

the grids closest to the iced airfoils are identical up to 0.6 chord length for both codes and

73

that is where all of the important flow features take place, there is a basis to compare the

two codes.

( l
“II-Illlllll-
III-III...-
_____

‘ .
- _.-———- ----—_ —

 

(C)

Figure 5-3.Grid about GLC 305 clean airfoil.

In this study, five different inner grids were used. The details of the inner grids
are as follows: For the business-jet airfoil without ice (referred to as clean), the inner
grid has 913 x 101 grid points (Figure 5-3). The second grid (referred to as 825) has the
944 ice shape represented by only 25% of the control points (not shown but similar to the
one shown in Figure 5-4). This means that there is smoothing of the jagged ice geometry,
which makes grid generation easier. For this grid system, the inner grid has 941 x 101

grid points. The third grid (referred to as SVl) uses 100% of the control points (Figure

74

5-4). For this grid, the 944 ice shape was smoothed less so that grid generation is more
tedious in having to capture more details of the jagged geometry. For this grid system,
the inner grid also has 941 x 101 grid points. The fourth grid (referred to as SV2) also
uses 100% of the control points. SV2 differs from SVl in having more grid points in the
region next to solid surfaces (989 x 129 grid points). The fifth grid generated is for the
212 ice, and it has 987 x 131 grid points (Figure 5-5).

All grid systems described above were generated by using transfinite
interpolation mum in the manner described in Refs. [l6],[17] and [18] with varying

degrees of local elliptic smoothing.

u‘.‘"4¢
a, fun,

a

o, w. a I out hon, " I
3*2”4”'r7'tz”‘~07‘}7'35“’01 "”345. ~73?
/

 

I ,
an” ‘ u nu
’fr ’0 ”fr,”lnto; '9 an» mm nu” In 1
«9:» Zita/”f” 4:91am, "'7' u ”7 ”3:"- 075' 4517.5 "mun

(c) (d)

Figure 5-4. Grid about GLC 305 airfoil with 212 ice.

75

mum

fl sh,“ ~ . .. .

“\o\ Mun-“at" mu

“'3 \\“ wrm.mm Immu‘fxii‘uw
\‘ "I

\ w
o“
. ““s§3’3“”$“‘“‘ m—zmmu We?”
*txsfil‘xnsa‘: .«sgi- ““

.0, ,,, 0,? 4, ~44, u- M.
r 22’ mm?
'% I ’49 "I "Wpooom
I I” l at”,
" %m£’~"m~uu

 

Figure 5-5. Grid about GLC305 airfoil with 944 ice.
5.5 Results

The main objective of this study is to assess how well CFD can predict lift, drag,
surface pressure, and the velocity field as a function of the angle of attack for a 2-D
airfoil with a glaze or a rime ice shape built upon its leading edge. The focus is on the

effects of turbulence modeling on the predictions.

5.5.1 Clean Airfoil

Figure 5-6, Figure 5-7 and Figure 5-8 show the results obtained by using the
WIND code with the S-A turbulence model for the GLC305 clean airfoil. These figures
show that CFD can predict the lift and surface pressure coefficients quite well for angles
of attack up to near stall. But the stall angle is slightly under-predicted. One reason for
not predicting the lift coefficient correctly near stall is that the flow under those

conditions becomes unsteady. In this study, only steady-state solutions were sought.

76

 

1.4

12 3 —t— Exp-Clean

 
 
  

1.0 ~
---0-~ CFD-Clean
0.8 i
0.6 r
0.4 J

0.2 —

CL

0.0 -
-O.2 1
-0.4 ..

 

 

 

-06 firrIIrIrIr
-6-4-20 246 810121416

AOA (Deg)

Figure 5-6. Computed and measured lift coefficient for clean GLC305 airfoil (M =
0.12, P = 20.5 psi, Re = 3.5 x 10°). CFD: WIND, S-A model.

77

 

 

 

 

 

 

 

 

 

-2.0,
i -------- Exp-Clean
'1'63 —CFD-Clean
-1.2§ ; I .
-0.8§
n- : ............................
0 -0.4: .............
0.0 f - ......... ..............
0.43
0.8—j
-0.1 0.0 0.1 0.2 0.3 0.4 0.5 0.6 0.7 0.8 0.9 1.0

x/c

Figure 5-7. Computed and measured pressure coefficient for clean GLC305 airfoil
(M = 0.12, P = 20.5 psi, Re = 3.5 x 106). AOA = 4. CFD: WIND, S-A model.

 

-2.o _

 

j .. -------- Exp-Clean
'15 : 7 ——CFD-Clean
-1.2 f " ‘

 

 

 

-0.8 f

-0.4 5

CP

 

 

0.0 {

 

0.4 f

0.8 {

 

 

 

 

II]IIIIIIITTTTTFTIITTTTIT—ITIIIIIIIIITITIiTIIIIII[IIII

-0.1 0.0 0.1 0.2 0.3 0.4 0.5 0.6 0.7 0.8 0.9 1.0
x/c

Figure 5-8 Computed and measured pressure coefficient for clean GLC305 airfoil
(M = 0.12, P = 20.5 psi, Re = 3.5 x 10°). AOA = 6. CFD: WIND, S-A model.

78

5.5.2 Airfoil with Glaze Ice

Figure 5-9 to Figure 5-12 show the results obtained by using the WIND code with
the S-A turbulence model for the GLC305 airfoil with the 944 glaze ice. These figures
show CFD to under-predict the lift and drag coefficients at all angles of attack. The
predicted pressure coefficient is shifted towards the leading edge, indicating incorrect
predictions of the separated region induced by the horns.

To ensure that grid resolution or ice-shape smoothing were not the cause of the
inaccuracies, simulations were performed with grids based on smoothed (825) and
unsmoothed (SVl) ice shapes as well as increased grid resolution (SV2). Figure 5-13
shows these effects to be insignificant. Note that SVl is already a very fine mesh, and is
the one used in all remaining simulations of the 944 ice with the WIND and Fluent codes.
Thus, using a smoothed ice shape (e.g., one represented by only 25% of the control
points) may be adequate when predicting lift and drag. This can be significant, since grid

generation is easier and less time consuming with a smoother ice surface.

79

 

 
 
   

 

 

 

0.8
—t—-Exp-lced
0.6 ~
~--0--- CFD-Iced
0.4 ~
0.2 -
.1
0
0.0 -
-0.2 -
-0.4 -
-O,6 r T I r I I I
-6 -4 -2 0 2 4 6 8 1O

AOA (Deg)

Figure 5-9. Computed and measured lift coefficient for GLC305 airfoil with 944 ice
(M = 0.12, P = 20.5 psi, Re = 3.5 x 105). CFD: WIND, S-A model.

0.16

 

Co

0.14 -

0.12 r

0.10 _

0.08 -

0.06 r

0.04 -

0.02 -

 

—I— Exp—Iced

~-0--- CFD-Iced

 

 

 

0.00

-2 0 2 4
AOA(Deg)

Figure 5-10. Computed and measured drag coefficient for GLC airfoil with 944 ice
(M = 0.12, P = 20.5 psi, Re = 3.5 x106). CFD: WIND, S-A model.

80

 

 

 

 

 

 

 

 

111IrrrfirIIIIIII.II|IrIr[rrfljr l

-0.1 0.0 0.1 0.2 0.3 0.4 .5 .6 0.7 0. 0. 1.0

T

 

Figure 5—1 1. Computed and measured pressure coefficient for GLC305 airfoil with
944 ice (M = 0.12, P = 20.5 psi, Re = 3.5 x106). AOA = 4°. CFD: WIND, S-A
model

 

 

 

 

rr—I—rl

l l I

-0.1 0.0 0.1 0. 0. 0.4 0.5 O. 0.7 0.8 0.9 1.0
x/c
Figure 5-12. Computed and measured pressure coefficient for GLC305 airfoil with
944 ice (M = 0.12, P = 20.5 psi, Re = 3.5 x106). AOA = 6°. CFD: WIND, S-A
model

rift illlllllf

81

325 fi SV1 lllllll SV2

 

 

 

AOA (Deg)

825 E SV1 M SV2

 

0.09 ,2, '* *fl'v
0.08
0.07
0.06
0.05
0.04 “5
0.03

l

 

CD

0.02
0.01 .
0.00 -"

 

 

AOA (Deg)

Figure 5-13. Effects of grid resolution and smoothing on CL and CD. 944 ice shape
and grid (M = 0.12, P = 20.5 psi, Re = 3.5 x 10°). CFD: WIND, S-A model.

82

In an attempt to improve predictions, more sophisticated turbulence models were
evaluated by using the Fluent code, which has more turbulence models encoded. Figure
5-14 to Figure 5-17 show that WIND and Fluent gave nearly the same results when the S-
A model was used on the same inner grid. This gives some confidence towards using two
different codes to evaluate a variety of turbulence models. Figure 5-14 to Figure 5-17
also show that S-A gives the best results for the lift and drag coefficients, confirming the
findings of Chung & Addy [4]. SST and the standard k-e models were found to under-
predict lift and to over-predict drag when compared to the other models. Differential
Reynolds stress and vz-f models gave the best results for the drag coefficient at AOA = 4°
but not at AOA = 6°. Also, the lift was severely under-predicted at AOA = 6°. The less
than satisfactory results for the vz-f and the differential Reynolds stress models is that

Fluent uses the one-equation Chen and Patel two-layer model in the near wall region.

5.5.3 Airfoil with Rime Ice

Figure 5-18 to Figure 5-21 show the results obtained by using the WIND code
with the S-A turbulence model for the GLC305 airfoil with the 212 rime ice. These
figures show CFD to predict accurately the lift, drag, and pressure coefficients until near
stall. This shows CFD predictions of airfoils with rime ice, where separated regions are
small, to be quite reliable. Similar to the situation with the clean airfoil, the stall angle is
slightly under-predicted. Again, the reason is that the flow becomes unsteady near stall,

and only steady-state solutions were generated in this study.

83

 

 

 

 

 

-6 -4 -2 0 2 4 6 810
AOA(Deg)

Figure 5-14. Computed and measured lift coefficient for GLC305 airfoil with 944
ice (M = 0.12, P = 37.0 psi, Re = 6.0 x 10°). CFD: WIND (S-A) and Fluent.

 

0.16
—e—Exp
0.14 i ---0---Wind-SA ,0
012 A —+—F-SA of,
-- - F-SST 9'

 
 
 
  

0.10 ‘ ma... F-KE
0.08 _ -->K-- F-RSM
—-- F-V2F

Co

0.06 r
0.04 r
0.02 r

 

 

0,00 l I I I I I I
-6 -4 -2 0 2 4 6 8 10

AOA (Deg)

 

Figure 5-15. Computed and measured drag coefficient for GLC305 airfoil with 944
ice (M = 0.12, P = 37.0 psi, Re = 6.0 x 10°). CFD: WIND (S-A) and Fluent.

84

fiAOA=4E§lAOA=6

 

  
 
  
  

 

 

ts

  

we

 
 

‘ W

 

Figure 5-16. Computed and measured lift coefficient for GLC305 airfoil with 944
ice at AOA = 4° and 6° (M = 0.12, P = 37.0 psi, Re = 6.0 x 10°).

EAOA=4 EAOA=6

 

 

   

w?

stN

F-SA F-SST

  

 

F-KE F-RSM V2F

Figure 5—17. Computed and measured drag coefficient for GLC305 airfoil with
944 ice at AOA = 4° and 6° (M = 0.12, P = 37.0 psi, Re = 6.0 x106).

85

 

1.0

0.8 e

0.6 r

0.4 4

CL

0.2 —

0.0 i

 

 

 

-O.4 TIIIIIWITI
~6-4-202 46 810121416

AOA (Deg)

 

Figure 548. Computed and measured lift and drag coefficients for GLCBOS airfoil
with 212 ice (M = 0.12, P = 20.5 psi, Re = 3.5 x 106). CFD: WIND, S-A model.

 

0.18

0.16 — ,9
---0---CFD a:
0.14 — .f ,
—-t—Exp
0.12 ~

0.10 5

Co

0.08 —
0.06 .
0.04 I
0.02 —

<9 ....................
000 T I l r I I

0 2 4 6 8 10 12 14
AOA (Deg)

Figure 5-19. Computed and measured lift and drag coefficients for GLC305 airfoil
with 212 ice (M = 0.12, P = 20.5 psi, Re = 3.5 x 10°). CFD: WIND, S-A model.

 

 

 

 

86

-2.8

 

-2.4
-2.0 3
-1.6 E
-1.2
-0.8 g
-o.4 g
0.0 g
0.4 i

0.8 j

CP

 

 

1.2

 

 

-O.1

0. 0.3 0.4 0.5 0.6 0.7 0. 0.9 1.0

I

 

ltll’!r)"l’¥_T||l 1 [fl

x/c

Figure 5-20. Computed and measured pressure coefficient for GLC305 airfoil with
212 ice as a function normalized distance along the chord (M = 0.12, P = 20.5 psi,
Re = 3.5 x 106). AOA = 4°. CFD: Wind, S-A model.

-2.8 7

 

-2.4 -,
-2.0 5
-1.6
-1.2 g

0 -0.8-;
-o.4 ;
0.0 E
0.4 g
0.8 g
1.2

 

 

—CFD 1

-~-><-~Exp

 

[\IVV

 

 

-0.1 0.0 0.1 0

l l T I l

. 0.3 o. 0.5 o. 0.7 o. o. 1.0
x/c

Figure 5-21. Computed and measured pressure coefficient for GLC305 airfoil with
212 ice as a function normalized distance along the chord (M = 0.12, P = 20.5 psi,
Re = 3.5 x 106). AOA = 6°. CFD: Wind, S-A model.

87

5.5.4 Prediction of the Velocity Field

Figure 5-22, Figure 5-23, and Figure 5-24 show the predicted contours of the x-
component velocity magnitude for the GLC305 airfoil with 944 glaze ice at AOA equal
to 0°, 4°, and 6°. Also shown are the experimentally measured values reported in Ref.
[20]. Comparing the CFD results with the measured ones shows that CFD incorrectly
predicts the size of the separated region downstream of the horn on the airfoil’s suction
side. This caused the surface pressure to be shifted, which in turn caused the lift and drag
to be under-predicted.

Figure 5-25, Figure 5-26 show the CFD and measured contours of the x-
component velocity magnitude for the GLC305 airfoil with 212 rime ice at AOA equal to
6° and 8°. These figures show the separated region to be predicted correctly though there
are still discrepancies between the CFD and the measured results. Since the lift, drag,
and pressure coefficients were predicted well by CFD, this indicates that predicting the

separated region correctly is paramount for CFD to get good results.

88

c
a
I

bbbooooa-a
hhbhhbbbhb

 

 

E
.5.

ooooaaa
hhbhbhh

' on

 

   

 

 

 

 

UN”
1A
1%
0.1 0',
.2 . as
>\ 04
02
09
02
l I I I I I I T I T I I I T I I I I I I I I I I I r I fiT f1 0'4
0.2 0.3 0.4 0.5 0.6
xlc
(C)
0.1
s3.
0
I I I I I I I I I I I I T I I I I I r I I I I I I I I T I rfi—‘
b o 0.2 0.3 0.4 0.5 o.
xlc

(d)

Figure 5-22. Normalized x-component velocity magnitude for GLC305 airfoil
with 944 ice at AOA = 0° (M = 0.12, P = 20.5 psi, Re = 3.5 x 106). (a) CFD with
WIND and S-A model; (b) Measurement; (c) CFD with WIND and S-A model

(Amplified); (d) Measurement (Amplified).

89

 

 

 

UNhf

14

12

1D

03

03

04

02

on

02

fihfllllI IIII'IIIIIIIIIIIIIIII 01‘

0.1 0.2 0.3 0.4 0.5 0.6 0.7 0.8 0.9
xlc
(a)

UNhf

Ea 13

£03

a 0.8

a 01

1:32..

5 0'2

04

 

  

 

I l l I I I I T I | I I I I T I I T I l I I I I I I I I I I

0.2

 

xlc
(C)

‘é
:

bbbooooaaa
hbbhbbbbhh

 

 

 

| I b I l I Iol.1l l I I012! I I lol-al I I IOI.4I I I Iol.5I I I lol-
xlc

(d)

Figure 5-23. Normalized x-component velocity magnitude for GLC305 airfoil
with 944 ice at AOA = 4° (M = 0.12, P = 20.5 psi, Re = 3.5 x 106). (a) CFD with
WIND and S-A model; (b) Measurement; (0) CFD with WIND and S-A model
(Amplified); (d) Measurement (Amplified).

90

 

 

 

 

 

 

 

 

U/Ulnf

1.4

1.2

1.0

0.8

0.6

0.4

0.2

0.0

-0.2

0.4

UNI"

1.4

1.2

1.0

0.8

0.8

0.4

0.2

-0.0

0.2

0.4
U/Ul'lf

1.4

1.2

1.0

, 0.8

0.8

0.4

0.2

0.0

02

I I I I I I l l I I I I I l I I I I I I I I l I I I I I I I 0'4

. 0.1 0.2 0.3 0.4 0.5 0.6
xlc
(C)

U/Uhf

1.4

1.2

1.0

0.8

0.8

0.4

0.2

0.0

0.2

l—T I I I I I I I I l I I I l I I I I I l I I I I I I I I I l 1 0"

0.1 0.2 0.3 0.4 0.5 o.
xlc

((0
Figure 5-24. Normalized x-component velocity magnitude for GLC305 airfoil
with 944 ice at AOA = 6° (M = 0.12, P = 20.5 psi, Re = 3.5 x 106). (a) CFD with
WIND and S-A model; (b) Measurement; (c) CFD with WIND and S-A model
(Amplified); (d) Measurement (Amplified).

91

C

bbboooo-aaa
hhbhhhhbhb:

 

 

S
:

bbbooooaaa
hhbhbbhbhb

 

E
:

bbbooooaae
bhbhhbbbbb

 

 

0.2 0.3 0.4 0.5 0.6
xlc

YBIIYIIIIIIIIIIIIIIltlrr1tllttr]

(C)

c
E

hhbbhbbbhhg
q

   

bbbooooaa-

 

III!lllllllllllllllllllllllllll

0.2 0.3 o .4
xlc

(d)

   

Figure 5-25. Normalized x-component velocity magnitude for GLC305 airfoil
with 212 ice at AOA = 6° (M = 0.12, P = 20.5 psi, Re = 3.5 x 106). (a) CFD with
WIND and S-A model; (b) Measurement; (0) CFD with WIND and S-A model

(Amplified); (d) Measurement (Amplified).

92

S
:

 
 

rxc'v! ”I? 5.3?!" I‘ ,‘J 13;..{.7J.,2I;.; 214’:.e§_.»;n’..

 
    

bbbooooaaa
bbbhhbhbhb

         

 

 

IIIIIIIIIII I‘I—IIIIWVYIIIlIlllIfir‘II

0213.3 0.4 0.5 0.6 0.7
xlc

 

(@

 

 

IIbIIIIIIIIIIIIIIIIIIIIIIIIIIIIITTIIIIIIII'IIIIIIIII

03 04 05 as OJ on as
xlc

 

 

 

 

o)
UN”
IA
12
on . 33
8 on
>. 0.4
02
03
02
I I I I I I I r I I I I I I I I I I I I I I 1 I77 I I I I I I I 0‘

0.2 0.3 0.4 0.5 0.6
xlc
(0
UN”
11
'12
1s
03
03
04
02
22
o I I I I I I I I I I I I I I I I I I I I I I I I I I I I I I I I I 0"
. 0.4 0.5
xlc
@)

Figure 5-26. Normalized x-component velocity magnitude for GL0305 airfoil
with 212 ice at AOA = 8° (M = 0.12, P = 20.5 psi, Re = 3.5 x106). (a) CFD with
WIND and S-A model; (b) Measurement; (c) CFD with WIND and S-A model

(Amplified); (d) Measurement (Amplified).

93

5.6 Summary

This study showed that if there are no large separated regions (e.g., for the 212
ice), CFD simulation based on the one-equation Spalart-Allmaras turbulence model can
accurately predict the lift, drag, and pressure coefficients. Thus, CFD can predict
aerodynamics of airfoils with rime ice quite adequately for angles of attack up to near
stall. For airfoils with glaze ice (e.g., 944 ice), where the horns produce large separated
regions about the airfoil’s leading edge, CFD predictions are much less satisfactory even
at low angles of attack. For airfoils with glaze ice, this study showed that even the v2-f
and the differential Reynolds stress models do not provide better results than the simple
S-A model, which was found to provide the best results. However, more study is needed
for the v2-f and the differential Reynolds stress models in which the near-wall treatment
is not the two-layer model of Chen and Patel. Comparing the predicted x-component
velocity magnitude with the measured ones show CFD to over-predict the size of the
separated region induced by the horns, and hence incorrectly predicting lift, drag, and
pressure coefficients.

This study also showed that for glaze ice, some smoothing of the ice shape is
acceptable, which makes grid generation easier.

Finally, it is noted that WIND and Fluent provided nearly identical results for lift,
drag, and pressure coefficients when the following were the same: grid, turbulence model

(S-A), and similar order of accuracy for the solution algorithms.

94

Chapter 6. A Comparative Study by Using CFD to Predict
lced Airfoil Aerodynamics

WIND, Fluent, and PowerF LOW were used to predict the lift, drag, and moment
coefficients of a business-j et airfoil with a rime ice (rough and jagged, but no protruding
horns) and with a glaze ice (rough and jagged and has two or more protruding horns) for
angles of attack from zero to and after stall. The performance of the following turbulence
models were examined by comparing predictions with available experimental data:
Spalart-Allmaras (S-A), RNG k-s, shear-stress transport, v2-f, and a differential Reynolds
stress model with and without non-equilibrium wall functions. For steady RANS
simulations, WIND and FLUENT were found to give nearly identical results if the grid
about the iced airfoil, the turbulence model, and the order of accuracy of the numerical
schemes used are the same. The use of wall fimctions was found to be acceptable for the
rime ice configuration and the flow conditions examined. For rime ice, the S-A model
was found to predict accurately until near the stall angle. For glaze ice, the CFD
predictions were much less satisfactory for all turbulence models and codes investigated
because of the large separated region produced by the horns. For unsteady RAN S,
WIND and Fluent did not provide better results. PowerFLOW, based on the Lattice
Boltzmarm method, gave excellent results for the lift coefficient at and near stall for the

rime ice, where the flow is inherently unsteady.

6.1 Introduction

For iced airfoils — even two-dimensional (2-D) ones — the generation of high-

quality structured gn'ds is a major challenge. Chi, et all”) presented a number of

95

methods to generate high-quality single- and multi-block structured grids for complicated
2-D ice shapes. Since multi-block grids can converge slower, Zhu, et allm studied the
effects of blocking strategy on the convergence rate to steady-state. Zhu, et a1. “81
examined the grid-generation and blocking techniques of Chi, et a1. “6] and Zhu, et a1. “8]
by applying them to a much more complicated ice shape, one with multiple, highly
extended, and closely packed horns.

Relatively few studies have focused on the effects of turbulence models. Chung,
et al.[2],[4] used WIND,I(’7]‘[68] an open source CFD code, to compare the performance of
several turbulence models, including the one-equation Spalait-Allmaras (S-A) modelml
and the two-equation shear-stress transport (SST) model.[27]’[28] They found the SST to
perform best with “clean” airfoil (i.e., airfoil without ice) and the S-A model to predict
better for an iced airfoil with horns. Chi, et a1.[69] studied the effects of turbulence models
by using WIND and the popular commercial code, Fluentml With WIND, they
examined the S—A modelm] and the SST model.[27]‘[28] With Fluent, they examined S-A,
SST, standard k-e,[29] VZ-f,[70]'[32] and a differential Reynolds stress model (RSM).[29H34]
They found WIND and Fluent to give essentially the same results if the grid about the
airfoil and the turbulence model used were the same. They also found the S-A model to
give better results than the more complicated SST model for iced airfoils. More
importantly, they found WHNID and Fluent with the S-A model to predict lift and drag
with good accuracy until near stall for airfoils with rime ice (i.e., ice shapes that have
only roughness and jaggedness but no protruding horns). But WIND with S-A and SST
and Fluent with S-A, SST, k-e, vz-f, and RSM predicted lift and drag much less

satisfactorily for airfoils with glaze ice (i.e., ice shapes with two or more protruding horns

96

near the airfoil’s leading edge). However, the conclusions described above were obtained

by evaluating only two angles of attack. Also, Chi, et 31(69]

generated only steady-state
solutions so that unsteadiness in the mean flow that may occur near stall or about the

horns of glaze ice were not considered.

6.2 Objective

The objective of this study is threefold. The first is to confirm that WIND and
Fluent do indeed generate nearly identical results if the grid, the turbulence model used,
and the order of accuracy of the numerical schemes used are the same for a range of
angle of attacks from zero to and after stall, where the flow can change substantially.
This confirmation will allow studies based on WIND to be compared with those based on
Fluent on a sound basis. The second is to examine S-A, SST, v2-f, RSM, and a
renormalization group (RNG) k-s modelnn’m] for a range of angles of attack for an
airfoil with rime ice and with glaze ice, again with and without the use of wall functions.
The third objective is to examine the unsteadiness that may exist in the mean flow by
performing time-accurate simulations (e.g., unsteady RANS or very large-eddy
simulation (VLES)).

Three codes will be used to meet these objectives: WIND, Fluent, and
PowerFLOW. WIND and Fluent are based on finite-volume methods that solve the
ensemble-averaged Navier-Stokes equations, closed by a turbulence model. Unlike
WIND and Fluent, PowerFLOW is based on a method, referred to as the Lattice
Boltzmann method WNW] to solve the Boltzmann equation. The accuracy of CFD
predictions will be assessed by comparing computed results with experimentally

measured lift, drag and moment coefficients.

97

The organization of the remainder of this chapter is as follows. Section 3
summarizes the glaze and rime, iced-airfoil problems studied. Section 4 outlines the
formulation and the numerical method of solution used in the codes and the turbulence
models examined. Section 5 presents the results generated. The key results of this study

are summarized in Section 6.

6.3 Problem Description

The airfoil, the ice shapes, and the flow conditions selected for study are those for
which the flow field is sufficiently complicated and for which there are experimental data
that can be used to assess the accuracy of the CFD predictions. The airfoil selected is the
business-jet airfoil (GLC305 [66]). The rime ice selected is the 212 ice shape, which has
considerable surface jaggedness but no protruding horns. The glaze ice selected is the
944 ice shape with two large protruding horns. The airfoil and the ice shapes about the
airfoil’s leading edge are shown in Figure 6-1 and Figure 6-2. See Ref. [66] for details of

the geometry.

Figure 6-1. GLCBOS airfoil with 212 rime icem

98

9;]

Figure 6-2. GLC305 airfoil with 944 glaze ice. I”

In this study, the freestream Mach number (M) is 0.12. Two freestream static
pressures (P = 20.5 psi and 37 .0 psi) and two Reynolds numbers based on the freestream

conditions and the chord length (Re = 3.5 x 106 and 6.0 x 106) are investigated.

6.4 Formulation and Numerical Method of Solution

Three different CFD codes were used to generate solutions for the iced-airfoil
problems described in the previous section. One is a widely used, open-source code,
known as WIND [671‘I68]. The second code is a popular commercial code, Fluentml The

[73],[74].[75]

third code is another popular commercial code, PowerFlow , which is based on a

fundamentally different method. These codes were selected because they are highly

versatile and contain a wide range of turbulence models.

6.4.1 WIND and Fluent
For the first two codes, the flow past the GLC305 airfoil with the 944 and 212 ice

shape is modeled by the ensemble-averaged conservation equations of mass (continuity),

99

momentum (full compressible Navier-Stokes), and energy for a thermally and calorically
perfect gas.

The turbulence models used for simulations of flow past the GLC305 airfoil with
the 212 ice shape are as follows. For WIND, the one-equation Spalart-Allmaras (S-A)
model was used. For Fluent, S-A turbulence model with and without non-equilibrium
wall functions were used [79]'[8”. More turbulence models were not used for two reasons.
First, Chi, et al. [69] has shown that the S-A model gives excellent results for rime ice until
near stall. Second, Chung, et a1. [2],”) has done an extensive comparative study on
turbulence models with WIND. Here, we want to compare predictions from WIND and
Fluent in which the grid about the airfoil is the same, the turbulence model is the same,
and order of accuracy is the same. If both codes yield identical results, then we can
examine results from WIND and Fluent as if they are from the same code.

The turbulence models used for the simulations of flow past the GLC305 airfoil
with the 944 ice shape are as follows. For WIND, again, the S-A was used. For Fluent,
the following turbulence models were investigated: S-A[24], SSTWHZSMZS], RNG k-
8 {mm}, Durbin’s vz-f modelml’m], and a differential Reynolds stress modell29]’[33]’[34].
With Fluent, the near-wall region is always modeled by enhanced wall treatment, which
uses low-Reynolds number models if the grid is sufficiently fine and wall functions if the
grid next to the wall is coarse.

The numerical methods of solution used are as follows. For WIND, the
convective terms were approximated by second-order Roe upwind differencing. Since
only steady-state solutions are of interest, time marching to steady state was

accomplished by an implicit method based on ADI-type approximate factorization with

100

local time stepping. For Fluent, which uses a finite-volume method, fluxes at the cell
faces are interpolated by using second-order upwind differencing scheme since we found
the accuracy of first-order discretization to be very poor although it generally yields
better convergence than the second-order scheme. The SIMPLE pressure-velocity
coupling algorithm was used to generate steady-state solutions. The convergence criteria
used was to require the scaled residuals to be less than 106 for the energy equation and
less than 10'3 for all other equations.

A limited number of time-accurate solutions were also generated by using WIND
and Fluent. For these unsteady RANS simulations, local-time-stepping were not invoked,

and a converged solution was obtained at each time step.

6.4.2 PowerF LOW

For PowerFLOW, the flow field is simulated by using a different set of equations.
PowerFLOW does not solve the “macroscopic” Navier-Stokes equations, which are the
ones solved in WIND and Fluent. Instead, it uses an extended Boltzmann kinetic
approach [76], and solves the “mesoscopic” equations, known as the Lattice Boltzmann
equation (LBE), that describes the kinetics of flow particles. The basic hydrodynamic
quantities (density, velocity, ...) are obtained through the moment summation of particle

density distribution functions [75]

. In this approach, sub-grid scale contributions to
turbulence are realized through an effective particle relaxation time scale, which can be
determined from the renormalization-group based transport equations (revised RNG k-
8 sub-grid model)[76]. The LBE based description of turbulent fluctuation carries flow

history and upstream information, and contains high order terms to account for the

nonlinearity of the Reynolds stress [741.[77]. A wall-shear stress model is used to reduce the

101

resolution requirement in the near wall region [74].

The unsteady flow solutions are
averaged over a representative time scale to generate mean flow characteristics.
In the numerical implementation, a BGK collision operator with a single

relaxation time approximation is used for the Boltzmann equations [76].

The particle
density distribution functions are cell centered, and the particle advection is solved by
explicit time marching, resulting in an upwinding scheme (due to the linear advection)

[7S],[78]

with second-order accuracy in space and time The RNG k-e equations are also

solved on the same lattice via a modified Lax-Wendroff-like explicit time marching finite

difference scheme [74].

6.4.3 Grid Systems for WIND and Fluent

In order to ensure proper comparison between the codes and among the
turbulence models, two criteria must be satisfied. The first is that the grid used for each
ice shape must be of high quality and provide enough resolution for grid-independent
solution. The second is that all codes and turbulence models must use the same grid for
each ice shape.

WIND and Fluent used essentially the same grid systems as explained below. For
WIND, all grid systems generated consist of two overlapping single-block grids. One is a
fine grid next to the airfoil, extending 0.6 chord length from the airfoil in all directions
(referred to as the inner grid, as shown in Figure 6-3(a)). The other is a coarser grid that
overlaps the fine grid by 0.1 chord length and extends 15 chord lengths away from the
airfoil in all directions (referred to as outer grid, as shown in Figure 6-3(b)). The inner
grid is the most important because that is where the flow physics is the most complicated.

While generating this grid, grid lines were clustered next to the airfoil surface to resolve

102

the turbulent boundary layer flow. Along the airfoil surface, equal arc-length was
employed to create a grid as smooth as possible. The details of the grids used are
described later in this section.

Since Fluent cannot handle overlapped grids, the grids used by WIND and by
Fluent will not be exactly the same. In this study, Fluent used the WIND’s inner grid so
that within 0.6 chord length next to the iced airfoil surface, the grids used by WIND and
by Fluent are identical. Away from 0.6 chord length, the grids used by WIND and Fluent
do differ (contrast Figure 6-3 and Figure 6-4). Since the most important flow physics
about the iced airfoil occur within 0.6 chord length, we do not expect this difference in
the grids beyond 0.6 chord length to be important. As results will show, this is indeed the

case.

    

,.» r~a~1-r4~+~+= , 71>

 

(a) ‘
Figure 6-3. Grid used by WIND for GLC305 airfoil. (a) Inner grid. (b) Outer grid.

With the above backdrop, the actual grids used for the airfoil with the 212 ice are
as follows. For WIND, the outer grid has 125 x 21 grid points, and the inner grid has 987
x 131 grid points (Figure 6-5 (a), (b)). For this high quality grid (smooth and nearly

orthogonal), the y+ of the first grid point away from the iced airfoil surface is within

103

unity for all angles of attack simulated. Also, the first 5 grid points have y+ values within
five. Thus, this is an extremely fine grid. When wall functions are used, the forty grid
lines next to the iced airfoil surface were removed so that the first grid point away has a
y+ between 30 and 50. For Fluent, the single-block grid used has 987 x 171 grid points,
when wall functions are not used. When wall functions are used, the grid has 987 x 140

grid points.

\‘ ‘\l

~ “' YMiir m umrliinllmullmlm lillll

. l‘ lii‘“ autumnal“. “um“mmiiiiidl‘li?
‘“ {mm-mt; .

 

”Mill
' 'mfli'lillnl‘l' “‘

Figure 6-4. Grid used by Fluent for GLC305 airfoil with 944 glaze ice.

The actual grids used for the airfoil with the 944 ice are as follows. For WND,
the outer grid has 125 x 21 grid points, and the inner grid has 941 x 101 grid points (Fig.
5 (c), (d)). Just like the grid for the 212 ice, the y+ of the first grid point away from the
iced airfoil surface is within unity for all angles of attack simulated. Also, the first 5 grid
points have y+ values within five. For Fluent, the single-block grid used has 941 x 141

grid points. For the 944 ice, WIND and Fluent did not use wall functions.

104

‘ l s \ ‘ ‘9
fiesta-etc:
f¢~§¢fio
‘\‘ w‘ “

ll,
4
: c

I a um
”'04 l

1
”"un ” "In " In "‘ I

0. Mn», u , 0114!».
0’55” fiHflJJ' ll ‘ “u'

“fiat“‘mnnmun '
“‘88“‘“ Wilton; “M:\‘|‘I‘.I‘l‘\‘|“w

M“

\ ‘6 n

W‘sggfiogww- --n-m-u--
\ \\ \fl‘oi‘ wu-

~
#-
‘2’.

,7, 49,,qu at", 'Qégnln-
We “w“"W m,
I Wit "It ”0:111:05...- ,.',
45%! M'm”'"'”' mum
I , . "Illlltulpullmu

(C)
Figure 6-5. Grid for GLC 305 airfoil: (a) Inner block grid for 212 rime ice; (b)
Close-up of (a); (c) Inner block grid for 944 glaze ice; (d) Close-up of (c).

 

The grid systems described above for WIND and Fluent were arrived at afier a
grid sensitivity study for the RSM model, which is the most stringent on grid
requirement. All grids were generated by using transfinite interpolation [91,[101 in the
manner described in Refs.[16],[l7], and [18] with varying degrees of local elliptic

smoothing. Gambit® was used to generate the outer grid for Fluent.

6.4.4 Grid System for PowerFlow
For PowerFLOW, a Cartesian grid structure is used for the flow simulation. The
fluid field is discretized into a set of regular cubic lattice cells and the original CAD

geometry (an STL format file for the iced airfoil) is overlaid on the Cartesian mesh to

105

represent the exact fluid/solid interface as shown in Figure 6-6 (left). In order to achieve
computation efficiency, variable resolution (VR) regions are applied as shown in Figure
6-6 (right). Here, each bounding box represents one grid resolution level and VRs
cascade outwards from the fine resolution region toward the coarse resolution region.
Resolutions differ by a factor of two between two adjacent VR regionsm]. The total

number of lattice cells used is 132,042 for the 212 ice and 104,002 for the 944 ice.

 

 

 

 

 

 

 

 

 

             

l .
crt‘ Tl 'r

Figure 6-6. Cartesian grid method (left) and Variable Resolution regions (right)
used in PowerFLOW.

6.5 Results

The main objective of this study is to assess how well CFD can predict lift, drag,
and moment as a function of the angle of attack for a 2-D airfoil with a rime ice shape
and a glaze ice shape. The focus is on the effects of turbulence modeling and code on the
predictions.

Before presenting the results, the abbreviations used in the plots are defined as
follows: Exp denotes measured experimental data from Ref. [1]. WIND, F, and LBM
denote results obtained by WIND, Fluent, and PowerFLOW code, respectively. SA,
RNG, SST, V2F, and RSM denote Spalart-Allmaras, , RNG k-s, shear-stress transport,

vz-f, and differential Reynolds stress model, respectively. EW and NW denote enhanced

106

wall treatment and non—equilibrium wall function options in Fluent. If a result is denoted
by EW or nothing is said about near-wall treatment, then the finest grids were used to
generate that solution and wall functions were not used.

The x-axis of all plots is the angle of attack (AOA). Note that every symbol on
the plot (e.g., +, A, C], ...) indicates a simulation has been done at that angle of attack.
The y-axis labels — CL, CD, and CM — denote the lifi, drag, and moment coefficients,
respectively. The moment coefficient was taken with respect to V4 chord distance from

the leading edge of the clean airfoil.

6.5.1 CFD Predictions of 212 Rime Ice Shape
6.5.1.1 WIND versus Fluent

Figure 6-7, Figure 6-8 show results obtained by using WIND and Fluent for the
business jet airfoil with the 212 ice for angles of attack from zero to past stall. These
figures show that if the grid about the iced airfoil, the turbulence model, and the order of
accuracy of the schemes are the same, then WIND and Fluent give nearly identical results
for the lift and drag coefficients are nearly at all angles of attack. This means that it is
possible to use comparative studies based on WIND and on Fluent to assess turbulence
models (i.e., findings based on one CFD code should also apply to another CFD code).
This is very comforting for all who work in CFD because the results should be the same
if the formulation and the grid are the same. These figures also show that the S-A model
provides excellent results until near stall. Both WIND and Fluent under-predict the stall

angle by about 2 degrees.

107

 

—e— Exp --+-- WIND-SA -~---~>< F-SA-EW

 

 

 

 

 

0.0 I l I l I f l
0 2 4 6 8 10 12 14 16

AOA (Deg)

Figure 6-7. 212 rime ice: lift coefficient = f (AOA).

—e— Exp --+- WIND-SA ---><---- F-SA-EW
0.18

0.16 — 2‘“
0.14 — ,vi'

0.12 e I
0.10 4
0.08 -
0.06 i
0.04 e
0.02 e

0.00 T l 1 l 1
2 4 6 8 1O 12 14

AOA (Deg)

 

Co

 

 

 

 

Figure 6-8. 212 rime ice: drag coefficient = f (AOA).

108

6.5.1.2 S-A versus RNG k-s with and without wall functions

Figure 6-9, Figure 6-10, and Figure 6-11 show the results obtained by using the
Fluent code with the S-A and the RNG k-s models and with different wall treatments.
From these figures, it can be seen that for the business-jet airfoil with the 212 ice, S-A
and RNG k-e models give very similar results for the lift, drag, and moment coefficients.
Also, it can be seen that the non-equilibrium wall functions can be used with reasonable
accuracy. By allowing wall functions to be used, the number of grid points can be

markedly reduced and convergence to steady state can be achieved with few iterations.
Thus, the S-A model with wall functions and the RNG k—e model with wall
functions perform as well as the S-A model that do not use wall functions. The reason
for comparing with the RNG k-e with wall functions is because PowerFLOW uses wall
functions, and only has the RNG k-e model. This study shows that this model is as good

as the S-A model.

109

—e—— Exp -------x F-SA-EW ~23“ F-SA-NW ~~>K~- F-RNG-NW
1.0

0.9 ~
0.8 -
0.7 -
0.6 .
0.5 ~
0.4 4
0.3 ~

 

CL

 

 

 

0.0 l I l 7 T I l
0 2 4 6 8 10 12 14 16

AOA (Deg)

 

Figure 6-9. 212 rime ice: lift coefficient = f (AOA).

—e— Exp ~--x F-SA-EW --A—~ F-SA-NW --->K--- F-RNG-NW
0.14

 

0.12 -
0.101

0.08 a

CD

0.06 —

0.04 .

 

0.02 —

 

 

0.00 T l I I
2 4 6 8 1O 12

AOA (Deg)

 

Figure 6-10. 2.12 rime ice: drag coefficient = f (AOA).

110

Q

—e— Exp m--)<-- F-SA-EW --A--- F-SA-NW ---)K--- F-RNG-NW

 

 

 

 

 

2 4 6 8 10 12 14
AOA (Deg)

Figure 6-11. 212 rime ice: moment coefficient = f (AOA).

6.5.1.3 Steady RAN S versus unsteady RANS and VLES

For the business-jet airfoil with the 212 ice shape, WIND and Fluent predict lift,
drag, and moment coefficient quite well at low angles of attack. It is only when the stall
angle is approached that lift is under predicted. Since the flow is known to become
unsteady as the stall angle is approached because of vortex shedding and unsteady
separation, time-accurate solutions may need to be performed; i.e., instead of steady
RANS simulations, we need to do unsteady RANS simulations. Efforts were made to
perform time-accurate computations with WIND and Fluent. However, for both of these
codes, unsteadiness did not appear to be significant. As a result, we wanted to explore

other codes that can resolve the unsteadiness with greater fidelity. The code that we

111

found success was the PowerFLOW code, which uses the Lattice Boltzmann method
(LBM) to do VLES.

Figure 6-12, Figure 6-13 show how PowerF LOW with RNG k-a sub-grid model
compare with the Fluent RNG k-a model. Figure 6-12 shows PowerFLOW to predict the
lift coefficient significantly better near stall. It also predicts the stall angle better, within
about 1 degree of the experimental data. However, Figure 6-13 shows the drag to be
predicted with less accuracy when compared to Fluent.

The improvement in the lift prediction might be that PowerFLOW is an unsteady
flow solver and is able to capture more flow physics at high AOA via VLES, where flow
field is quite unsteady. The difference in the drag prediction maybe due to the fact that
the two codes use different models to describe the near wall flow region, which may lead
to different surface pressure distributions near the leading edge flow separation region,
and the airfoil drag coefficient, is quite sensitive to such a difference in Cp distributions.

Further examining of the wall models used in these two approaches are needed.

112

—6— Exp ---)K~- F-RNG-NW —-—— LBM

 

 

 

 

0.0 I I I I I I I
0 2 4 6 8 10 12 14 16

AOA (Deg)
Figure 6-12. 212 rime ice: lift coefficient = f (AOA).

 

—6— Exp --—>K~— F-RNG-NW —-— LBM
0.14

 

0.12 ~

0.10 -

0.08 e

Co

0.06 ~

0.04 —

0.02 ~

 

 

 

0.00 r l l r
2 4 6 8 10 12

AOA (Deg)

 

Figure 6-13. 212 rime ice: drag coefficient = f (AOA).

113

 

6.5.1.4 WIND vs. Fluent vs. PowerFLOW

Figure 6-14, Figure 6-15 compares all results generated for the business-j et airfoil
with the 212 ice by using WIND, Fluent, and PowerFLOW. From Figure 6-14, it can be
seen that all codes and models examined predict well at low angles of attack. At angles
of attack near or alter stall, the lift is always under predicted. Only PowerF LOW with
LBM predicts lift well near stall.

Figure 6-15 shows WIND and Fluent to predict drag with reasonable accuracy.

PowerFLOW, however, does worse here. Thus, more researches are needed here.

—e—Exp --+—-WlND-SA wa-SA-EW
--a—--F-SA-NW ------x F-RNG-NW ——— LBM

1.0
0.9 -
0.8 1
0.7 —
0.6 a
0.5 -
0.4 4
0.3 J

 

CL

 

 

 

 

0.0 I I I I I I I
0 2 4 6 8 1O 12 14 16

AOA (Deg)

Figure 6-14. 212 rime ice: lift coefficient = f (AOA).

114

It”.

+Exp --+- WIND-SA x--F-SA-EW

MA F-SA-NW x F-RNG-NW —— LBM
0.20
0.18 —
0.16 —
0.14 ~
0.12 —
0.10 ~
0.08 —
0.06 «
0.04 —
0.02 7:

0.00 l r a l l
2 4 6 8 10 12 14

AOA (Deg)

 

Co

   

 

 

 

Figure 6-15. 212 rime ice: drag coefficient = f (AOA).

6.5.2 CFD predictions of 944 glaze ice shape

So far, we have only examined the capability of CFD in predicting lift, drag, and
moment of an airfoil with rime ice, which are ice shapes that have roughness and
jaggedness but no protruding horns. In this section, we examine how well CFD can
predict the lift, drag, and moment of an airfoil with glaze ice, which has roughness,
jaggedness, and large protruding horns. Once there are two or more horns, the flow
becomes considerably more complicated. For glaze ice, all simulations performed here
by WIND and Fluent did not use wall functions because we wanted to resolve the near-
wall region flow features as accurately as possible. With PowerFLOW, however, wall

fiinctions are still used.

115

6.5.2.1 WIND versus Fluent

Similar to our study on rime ice, we want to make sure that for glaze ice, WIND
and Fluent will provide the same results if the grid about the iced airfoil, the turbulence
model, and the order of accuracy of the schemes used are the same. Figures 16 and 17
show that this is indeed the case.

Figure 6-16 shows the S-A model to under-predict the lift coefficient even at
fairly low angles of attack. This is because when there are horns, large separated regions
form even at zero angle of attack. Figure 6-17 shows the S-A model to under predict

drag at all angles of attack. Only the trend is predicted correctly.

—e— Exp ---o----WIND-SA --+—- F-SA
0.7

 

0.6 ~

0.5 a

0.4 J

CL

0.3 ~

0.2 -

 

0.1 —
E

 

 

0.0 l l 1 l
0 2 4 6 8 1O

AOA (Deg)

 

Figure 6-16. 944 glaze ice: lift coefficient = f (AOA).

116

+ Exp ---0- WIND-SA --+—- F-SA

 

0.16
0.14 . ..o
0.12 i .o

0.10 ~
0.08 ~

CD

0.06 4
0.04 -

 

0.02 ~
0.00 -

 

 

-0.02 r w l
-2 0 2 4 6 8 10

AOA (Deg)

 

—I
d

Figure 6-17. 944 glaze ice: drag coefficient = f (AOA).

6.5.2.2 Comparison among Fluent’s different turbulent models

Since the S-A model is only a one-equation model, perhaps more advanced
models such as RNG k-c, SST, vz-f, and RSM can yield better results. In Chapter 0, we
did not find this to be the case. But, we only evaluated these models at two angles of
attack. In this chapter, these models are evaluated at angles of attack from zero to after
stall. Enhanced wall treatment option is selected for all Fluent simulations.

The results generated are shown in Figure 6-18, Figure 6-19, and Figure 6-20.
These figures show that different turbulence models give quite different prediction, which

implies turbulence modeling is the key part for further improvement.

117

 

Figure 6-18 shows the S-A and the RNG k-e turbulence models to give the best
CL predictions and SST turbulence model to be second. Surprisingly, vz-f and RSM
models did not provide satisfactory results on lift.

Figure 6-19 shows the v2-f and the RSM models to give the best results on the
drag coefficient. While v2-f model is a little better than the RSM model. S-A, SST, and
RNG k-s models performed poorly.

Figure 6-20 shows predictions of the moment coefficient. Predictions by vz-f and
RSM turbulence models match experimental data much better than the results by S-A,
SST, RNG k-s turbulence models. Here, it is noted that for these steady RANS
simulations, vz-f and RSM are very hard to converge. S-A, SST, and RNG k-e converge

relatively easily.

—e— Exp --+—- F-SA ----><--- F-SST
—-A—~ F-RNG —--->l<- F-RSM —-— F-V2F

0.7

 

0.6 -

0.5 -

 

 

 

 

AOA (Deg)

Figure 6-18. 944 glaze ice: lift coefficient = f (AOA).

118

—e— Exp ---l-- F-SA ~--~x- F-SST
--A--- F-RNG --)K-- F-RSM —— F-V2F

 

0.12

0.10 —

0.08 -

Co

0.06 A

0.04 —

 

 

 

 

AOA (Deg)

Figure 6-19. 944 glaze ice: drag coefficient = f (AOA).

+Exp -—+—- F-SA -><---- F-SST
--A--- F-RNG ---x-- F-RSM —— F-V2F
0.08
0.07 -
0.06 —
0.05 ~
0.04 -
0.03 -
0.02
0.01 —
0.00 -
.001 _ if:
002 l l .

 

CM

 

 

 

 

AOA (Deg)

Figure 6-20. 944 glaze ice: moment coefficient = f (AOA).

119

6.5.2.3 Comparison F—RNG with LBM

Figure 6-21 shows time-averaged results from time-accurate simulations
performed by using PowerF LOW. From this figure, it can be seen that lift is predicted a
little better. Figure 6-22 shows PowerFLOW to predict drag reasonably well at low
angles of attack, but slightly worse than RNG at high angles of attack. The inability of
PowerFLOW to perform better for glaze ice may be due to its wall function treatment of

the near-wall region.

—e— Exp ---A--- F-RNG —El— LBM
0.7

 

0.6 -

0.5 ~

0.4 ~

CL

0.3 ~

0.2 ~

 

 

 

 

AOA (Deg)

Figure 6-21. 944 glaze ice: lift coefficient = f (AOA).

120

—e— Exp “A“ F-RNG —B—- LBM

 

0.12
0.10 i

0.08 ~

CD

0.06 *

 

 

 

 

-2 0 2 4 6 8
AOA (Deg)

Figure 6-22. 944 glaze ice: drag coefficient = f (AOA).

6.5.2.4 Comparison of all results of 944 ice shape

Figure 6-23, Figure 6-24 compares all results generated for the business-jet airfoil
with the 944 glaze ice by using Fluent and PowerFLOW. These figures show that

PowerFLOW with the LBM method predicts lifi better. The v2-f and the RSM models

predict drag better.

121

—e—Exp ——+—- F-SA ------x- F-SST —-A—-- F-RNG
--->K- F-RSM —— F-V2F —EI— LBM

0.7

 

0.6 -

0.5 ~

0.4 ~

CL

0.3 ~

0.2 ~

0.1 -

 

 

 

 

AOA (Deg)
Figure 6-23. 944 glaze ice: lift coefficient = f (AOA).

—e——Exp ——+—- F-SA ----x---iF-SST ~A—~-F-RNG
-x- F-RSM —-—— F-V2F —B— LBM

0.12

 

0.10 -

0.08 a

CD

0.06 -

0.04 i

 

 

 

 

-2 0 2 4 6 8
AOA (Deg)

Figure 6—24. 944 glaze ice: drag coefficient = f (AOA).

122

 

 

6.5.2.5 Grid sensitivity study

Grid sensitivity studies has been given for AOA = 4, 6 of both ice shape. Studies
show that CL is nearly not changed after the refinement of the grid, but CD, CM may

change 1%.

6.6 Summary

For the 212 rime ice and the 944 glaze ice, if the grid used about the iced airfoil,
the turbulence model, and the order of accuracy of the numerical schemes used are the
same, then the results obtained are essentially identical whether one uses WIND, Fluent,
or PowerFLOW. Thus, government and commercial CFD codes are now like LDV and
PIV system that you can buy or license, but still must use correctly to get meaningful
results.

For the 212 rime ice, WIND, Fluent, and PowerFLOW all gave excellent results
for lift except at angles of attack near stall. At angles of attack near and after stall,
PowerFLOW gave the best results because the unsteady mean was resolved by VLES.
On drag, WIND and Fluent provided excellent agreement with experimental data. If the
Reynolds number of the flow is high so that grid lines at y+ of 30 to 50 are still very
close to the iced airfoil surface (i.e., the key features of the ice geometry is still resolved
by the coarser mesh), then the use of non-equilibrium wall functions were found to yield
results similar to those from low Reynolds turbulence models.

For the 944 glaze ice with two large protruding horns in addition to surface
roughness and jaggedness, CFD yielded much less satisfactory results. Among the

RANS turbulence models, S-A gives the best lift predictions followed by RNG k-e. vz-f

123

and RSM give much less satisfactory results on lift, but provide the best drag prediction.
PowerFLOW with the LBM gives the best lift prediction by resolving the unsteadiness
that may occur, but prediction on drags could be improved. For PowerFLOW, a low
Reynolds number near-wall model is needed.

Further studies are needed to investigate the unsteady solver influence and wall

model influence on flow prediction of iced airfoils.

124

Chapter 7. CFD Study of Turbine Blades with Real
Rough Surface

The main objective of this research is to directly compute the skin friction (Cf) and
heat transfer (St) coefficients on real rough surfaces using a state-of-the-art unstructured
adaptive grid-based finite volume method. Recent experiments with real roughness

panels by Bons [45]

were computationally simulated in this study. Computational results
were compared with experimental data to assess the simulation accuracy. A RANS
(Reynolds-Averaged Navier-Stokes) approach based on the Spalart-Allmaras turbulence
model and a DES (Detached Eddy Simulation) approach were employed for the
computations, and grid refinement studies were conducted to assess the effects of grid
resolution. In two cases with rough surfaces, the RANS approach was capable of
accurately predicting Cf (within 3.5%) while under-predicting St by 8-15%. The DES
approach was able to predict Cf and St for smooth plate but failed in the cases with real
roughness. The cause will be further investigated.

This chapter is organized as follows: First, the experimental setup is reviewed to
set the stage for computational simulations. Second, the viscous adaptive Cartesian grid
generation approach is briefly described. Third, the basic features of the finite-volume
Navier-Stokes solver with a RANS Spalart-Allmaras (S-A) model and the DES approach
are presented. Forth, several validation cases are presented and direct simulations

including the surface roughness are performed, and the results are compared with the

experimental data. Finally, conclusions from the present study are summarized.

125

 

7.1 Experimental Setup

7.1.1 Roughness Surface Measurement and Fabrication
To prepare for our current study, many land-based turbine components were
assembled fiom four turbine manufacturers: General Electric, Solar Turbines, Siemens-
Westinghouse, and Honeywell Corporation. They are selected by the manufacturers as
representative of general surface conditions of land-based gas turbine. 3-D surface
measurements by contact stylus measurement system were made on these components to
get the “real” roughness surface. Of the many 3-D maps obtained [44], six maps were
selected for our study. They included one pitted surface, two coated/spalled surfaces, one
fuel deposit surface, and two erosion/deposit surfaces. Surface #4 and #6, out of the six
are selected in this study. Surface #4 is an example of surface with fuel deposits that are
elliptical in shape and aligned with the streamwise direction. Surface #6 is representative
of combined erosion and deposits surface with smaller, more jagged roughness elements
than surface #4. Nikuradse [5” classified roughness into three regimes based on k+:
aerodynarnically smooth (k+<5), transitionally rough (5<k+<70), and completely rough
(k+>70). The surface data were properly scaled to make sure the scaled model and the
actual parts have the same roughness regime as defined by k+. Finally, plastic roughness
models were fabricated by a StrataSys Inc. GeniSys Xi 3-D printer. Figure 7-1 shows the
geometry characteristic of surface #6. Figure 7-2 shows the geometry characteristic of

surface #4.

126

 

Figure 7-2. Geometry characteristic of surface #4.

127

7.1.2 Wind Tunnel Facility

For a detailed description of the wind tunnel facility, refer to Ref. [45]. A very
brief introduction is given here. Figure 7-3 is the schematic of the flat plate wind tunnel
used in Bons’ experiment. Figure 7-4 shows the dimensions of the test section with the
location of six roughness panels. The leading edge of the boundary layer starts fiom the
boundary layer suction point. The cross-section area of the roughness panel section is
240 mm by 380 mm. The leading edge of the roughness panel section is located 1040
mm from the boundary layer suction point. Typically, six individual roughness panels
(140 mm length x 120 mm width) are installed in a 280 mm stream-wise gap in the
bottom wall. Figure 7-5 shows the geometry of the roughness surface section after six

roughness panels are arranged together there. The tunnel then continues 620 mm beyond

the trailing edge of the roughness panels.

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

HOt Air from Jet
l Main Blower Turbulence
Grid 7
Conditioning
Plenum Rough Plate
__ . 71 .
Turbulence =\L._
Blower Boundary
Layer
0 . =1: . Bleed
I: i“
77//////// 7///// ////////// //7///////// /7///.

Figure 7-3. Schematic of the Flat Plate Wind Tunnel in Heat Transfer
Measurement Configuration. (Bons GT-2002-30198)

128

  
 
  

Position ofthe
infrared camera

--——_.

   

Six individual
roughness panels

Figure 7-4. Dimensions (mm) of the test section with the location of six roughness
panels.

 

Figure 7-5. Geometry of the roughness surface section.
7.1 .3 Cf Measurement
It is difficult to make precise drag measurements over rough surfaces. Acharya et

al. [52] suggested a force-balance method rather than boundary layer momentum balance

129

and log-region curve fitting methods which are velocity based Cf measurement methods.
In Bons’ experimental study, a hanging element balance (As shown in Figure 7-6) was

used to obtain Cf.

7.1.4 St Measurement

To measure St, a FLIR Therrnacam SC 3000 infrared camera system was
mounted with lens fit into a hole in the ceiling of the trmnel. The camera field of view is
roughly 70x90 mm. The limited field of view was centered at a distance of 1200 mm
from the leading edge of the tunnel floor. The surface temperatures measured from the
camera were area-averaged to obtain the representative surface temperature history to
calculate St. The St was then determined from this surface temperature history using the

method of Schultz and Jones [53] which is based on Duhamel’s superposition method.

I? III

 

 

      

 

 

 

 

Top Wall of Tunnel
Steel
‘/ Suspension
Wires
Uinf
Upstream fl
Overlapping .
Shim Capacrtec
Roughness Panels on Metal Support Capacitance
/ Meter
. I
Plexrglass Plexiglass
Upstream Downstream
Pressure Pressure

Figure 7-6. Schematic of floating panel Cf measurement apparatus on AF UL wind
tunneL

130

7.2 Viscous Cartesian Grid Generation Approach

The first step in a CFD simulation is to define or import the geometry, and
generate a computational grid. In a viscous Cartesian grid method, a volume grid is first
generated before a surface grid is produced through projections. A unique advantage of
the method is that “dirty” geometries may be automatically handled without geometry

I47]

repair The generation of a viscous Cartesian grid can be accomplished in the

following steps:

7.2.1 Adaptive Cartesian Grid Generation

Two meshing parameters, dmin and dmax are specified first. They represent the
minimum and maximum sizes of Cartesian grid cells to be generated. One of the popular
data structures for adaptive Cartesian grids is the Octree. A more flexible data structure
is the so called 2N tree, which supports anisotropic subdivisions. The adaptive Cartesian
grid is generated by recursively subdividing a single coarse root Cartesian cell. Since the
root grid cell must cover the entire computational domain, the surface geometry is
contained in the root cell. The size of the Cartesian cells intersecting the geometry is
controlled by two parameters, disT and disN. Parameter disN controls the Cartesian cell
size in the geometry normal direction, whereas disT specifies the Cartesian cell size in
the geometry tangential direction. The ratio disT/disN determines the maximum cell
aspect ratio in the Cartesian grid. The recursive sub-division process stops when all the
Cartesian cells intersecting the geometries satisfy the length scale requirement. For the
sake of solution accuracy, it is very important to ensure that the Cartesian grid is smooth.
In the present study, the sizes of any two neighboring cells in any coordinate direction

cannot differ by a factor exceeding 2.

131

7.2.2 Cartesian Grid Front Generation and Smoothing

In order to “insert” a viscous layer grid between the Cartesian grid and the body
surface, Cartesian cells intersected by the geometry must be removed, leaving an empty
space between the Cartesian grid and the body surface. All the Cartesian cells intersected
by the geometry can be determined efficiently using a tree-based search algorithm. In
addition, the intersected cells also serve to divide cells “outside” the geometry from the
cells “inside” the geometry. Depending on whether the problem is external or internal,
cells “inside” or “outside” the geometry must be removed. The 2N tree is not only used to
record the recursive cell subdivision process, it is also used to perform efficient
intersection operations with the geometry. For example, if a (coarse) Cartesian cell does
not intersect a geometric entity, all of the child cells from the Cartesian cell must not
intersect the geometric entity.

Once the Cartesian cells intersected by the geometry, and cells outside the
computational domain are removed, we are left with a “volume” Cartesian grid. The
boundary faces of this volume Cartesian grid form the so-called Cartesian front. Before
this front is “projected” to the geometry, it is smoothed with a Laplacian smoother to
produce a smoother front. To prevent the smoothed Cartesian front from intersecting the
body geometry, Cartesian cells which are within a certain distance of the body are also

removed.

7.2.3 Projection of the Cartesian Front to the Body Surface
After the smoothed front in the Cartesian grid is obtained, each node in the front
needs to be connected to the body surface to form a single layer of viscous grids. After

the front is projected to the boundary geometric entities, a "water-tight" surface grid is

132

 

 

generated on the boundary. The “foot prints” of the layer grids on the body surface have
the same topology (or connectivity) as the Cartesian fi'ont. With this assumption, the
viscous layer grids are naturally “blended” with the adaptive Cartesian grid, eliminating
the need of cell-cutting currently adopted by many Cartesian grid generators. By
connecting each point on the Cartesian front and the corresponding projected point on the
boundary, we obtain a single layer of prism grids. This single layer can be sub-divided
into multiple layers with proper grid clustering near the geometry to resolve a viscous

boundary layer.

 

Figure 7-7. Geometry of Roughness Surface #4 used in the CFD study.
Roughness panel mirrored in the stream-wise direction.

An example of viscous adaptive Cartesian grid for 2 panels of Surface # 4 (shown

in Figure 7-7) is displayed here. Figure 7-8 shows two cutting planes across the

133

computational grid. Note that the viscous layer grid is used to resolve the turbulent
boundary layer. The surface grid generated from front projection is displayed in Figure
7-9. It is observed that the grid cells near the roughness panels are adaptively refined to

resolve the roughness elements.

 

Figure 7-8. Cutting planes showing the viscous adaptive Cartesian grids.

 

Figure 7-9. The surface grid on the lower channel wall showing the refinement
near the rough panels.

134

7.3 Numerical Method

A flow solver capable of handling arbitrary polyhedrons was developed to
uniformly handle the adaptive Cartesian and the viscous layer grids [47]. The so-called
hanging node problem actually disappears because of the use of a cell-centered finite-
volume method supporting arbitrary grid cells. A Cartesian face with a hanging node is
actually treated as four separate faces. The hanging nodes are, in fact, not visible to the
flow solver. This simple treatment is not only accurate, but fully conservative as well.

The Reynolds-averaged Navier-Stokes equations can be written in the following

integral form:

6Q _ _
lEdVarim mars—0 (1)

where Q is the vector of conserved variables, F and Fv are inviscid and viscous
flux vectors, respectively. The integration of Eq. (1) in an arbitrary control volume, V,
gives:

d .
'(7QtLI/t +2}?fo =ZFv,fo (2)
f f

where Q. is the vector of cell-averaged conserved variables, Ff and Fv; are the

numerical inviscid and viscous flux vectors through face f, and Sf is the face area. The
overbar will be dropped from here on. The key is then how to compute the inviscid and
viscous fluxes through any given face. Here the standard Godunov-type finite volume
approach is employed. Using a linear least-squares reconstruction algorithm, a cell-wise
linear distribution can be built for each solution variable (in the present study the

primitive variables). To compute the inviscid flux, an approximate Riemann solver such

135

 

 

as Roe flux difference splitting [54] is used given the reconstructed solutions at both sides
of a face. To handle steep gradients or discontinuities, a limiter due to Venkatakrishna

[57] is used. The viscous flux is computed using a simple and robust approach presented

in Reference [55]

without a separate viscous reconstruction.

Although explicit schemes are easy to implement, and are often useful for steady-
state, inviscid flow problems, implicit schemes are found to be much more effective for
viscous flow problems with highly clustered computational grids. An efficient block LU-
SGS (Lower-Upper Symmetric Gauss-Seidel) implicit scheme [58] has been developed for
time integration on arbitrary grids. This block LU-SGS (BLU-SGS) scheme takes much
less memory than a fully (linearized) implicit scheme, while having essentially the same
or better convergence rate than a fully implicit scheme. The BLU-SGS scheme can be
used to integrate Eq. (2) with both first or second order accuracy. For steady flow
computations, the backward Euler approach is employed, i.e.,

{1+1 ____Q__i +1 n+IS
At ————V+ZF" Sf: 73va (3)
For time accurate computations, we employ a very robust second-order backward

difference scheme

 

n+1 n n—l
-3 . + .
4,-Q 21ft! Q1 Vi + 2f ijLle = Eijv’ZIISf (4)

To further speed convergence, local time steps are used in Eq. (3). Both Eq. (3)
and Eq. (4) are then solved with the BLU-SGS approach. Multiple sub-iterations are

utilized to solve Eq. (4) to improve time accuracy.

136

 

 

 

To simulate flow turbulence, a RANS Spalart-Allmaras (S-A) model and a DES
approach are employed. They are briefly described next. The S-A one-equation model
[241126] solves a single partial differential equation for the variable V which is related to
the turbulent viscosity. The differential equation is derived by using empiricism and
arguments of dimensional analysis, Galilean invariance and selected dependence on the

molecular viscosity. The model includes a wall destruction term that reduces the

turbulent viscosity in the log layer and laminar sub-layer. The equation can be written in

the following form:

~

~ 2
£=Cb1SV—Cw1fw[i] +l-[V.((V+I7)Vl7)+cb2(VV)2]
0’ d a (5)
The turbulent viscosity is determined via,
13 V
Vt=va1, fv1=—3——3—, 15— (6)
I + cvl V

where v is the molecular viscosity. Using S to denote the magnitude of the

vorticity, the modified vorticity is defined as

 

~ {7 Z
SES+——-—f2, f2=1— (7)
K2612 v v 1+3fvl

where d is the distance to the closest wall. The wall destruction function is

definedas
1 6 1/6 ~
+c V
fw=g[——6 W63] , g=r+cw2(r6—r), re. 2 2 (8)
g +Cw3 SK (1

The closure coefficients are given by:

137

 

 

em = 0.1355, 0' = 2/3, cb2 = 0.622, x = 0.41,

c 1+c
Cw] = -b%'+ 0b.? , CW2 = 0.3, CW3 = 2, CV1 = 7.1

K

 

For the DES approach, a new length scale is defined for DES, i.e.,
c7 = min(d,CDEsA) (9)
where CDES is a constant, and A is the measure of local mesh spacing, taken to be

the maximum distance from the current cell centroid to the centroids of its neighbors.

Then the distance to the closest wall (I in the S-A model is replaced with the new length
scale (7 to obtain DES. The purpose of using this new length is that in boundary layers,
A far exceeds (1 and the standard S-A model rules since 3 = d. The model comes with its
experience base and is fair accuracy. Away from walls, we have (7 = C 053A and the

model turns into a simple one equation Sub—Grid-Scale (SGS) model, close to
Smagorinsky’s in the sense that both make the “mixing length” proportional to A. On the
other hand, the approach retains the full sensitivity to the RANS model’s predictions of
boundary separation. The model constant CnEs was calibrated at 0.65 through the study

of isotropic turbulence [58].

7.4 Results and Discussions

The main objective of this study is to assess how well CFD can predict Cf, St for
real rough surfaces by comparing computational results with experimental results and
correlation formulas. First, turbulent flow over a smooth flat plate in a wind tunnel is
computed to validate the finite volume flow solver. Following the validation,
computations with rough surfaces in the wind tunnel are performed, and comparisons are

made to assess the adequacy of the CFD results.

138

.s-

 

 

7.4.1 Turbulent Flow over a Flat Plate in a Wind Tunnel

This case serves two purposes. One purpose is to validate the S-A model and DES
approach in the flow solver. The other is to see how well CFD results agree with
correlation formulas and experimental data for smooth walls. The computational domain
is the wind tunnel shown in Figure 7-4. To further reduce the computational cost, only
1/6 of the span is included in the computational domain. Symmetry or slip wall boundary
conditions are used on the two end walls in the span-wise direction. Therefore the flow is
essentially two-dimensional. Two viscous adaptive Cartesian meshes were generated for
this case. The coarse mesh has 37,888 cells while the fine mesh has 58,368 cells. The
average y+ value of the first cell from the wall is 0.8 for the coarse mesh and 0.3 for the
fine mesh. Simulations were carried out using the S-A models on both meshes. Figure
7-10 shows the computed average Cf and St over the area where roughness panels are
located and the experimental and correlation values. Standard flat plat correlations for Cf
and St [62] are calculated as follows:

c 0.026
f = 1
Rex”

 

050/:
_ Pr, + Josef [5 Pr+ 51n(5 Pr+ 1) — 14 Pr,]

 

 

5:

Note that the computed cf and St on the coarse grid agree quite well with those on
the fine grid indicating the grid resolution is adequate for the present simulation. The
computed Cf and St are also in good agreement with the experiment and correlation
results. The computational results suggest that the S-A model is capable of predicting

both Cf and St for the flat plate case.

139

.4 u.—

 

0.0040
0.0035 i
0.0030 -
0.0025 -

 

+ + XCf
+31

0.0020 J + +

 

 

 

0.0015 1
0.0010 1
0.0005 ~

 

 

0.0000 1 1 1

Corase Mesh Fine Exp Correlation
(y+=0.8) mesh(y+=0.3)

 

 

Figure 7-10. CFD simulation results compared to experimental data and standard
roughness correlations.

0.10

 

 

— Cf S—A Steady
- - - - St S-A Steady
— Cf DES‘I
—— St DES1
—— Cf DESZ
— — - St DESZ

 

 

 

 

 

 

\-
“q
A
Q‘
-—~-—-———————--.--—P-———-_-—q

 

 

 

0 1000 2000 3000
Iterations

Figure 7-11. Convergence histories for flat panels with S-A and DES approaches

140

Next, both S-A and DES approaches were studied and compared on the fine
mesh. For the S-A model, a local time stepping strategy with a CF L of 50 was employed.
For the DES approach, the simulation must be carried out in a time-accurate manner.
Therefore to study temporal convergence, two cases were performed with different time
steps. In the case of DES1, the time step is 2.457e-4 (which corresponds to a CFL
number of 1000 for the smallest cell), while the time step is doubled for the case of
DES2. The convergence histories of Cf and St from all three simulations are plotted in
Figure 7-11 in terms of the number of iterations or time steps. In terms of physical time,
both DES cases showed identical convergence histories, indicating the simulation is time-
step independent. Note that convergence in Cf and St was achieved in about 2000
iterations for the S-A model, and in about 200 time steps for the DES approach. The
convergence for the DES approach is a lot faster because of the larger time step in the
viscous boundary layer and the use of multiple sub-iterations to achieve time accuracy.
Figure 7-12 shows the computed average Cf and St with comparison to experimental data
and correlation formula. It is observed that the DES approach predicted a slightly lower

value of Cf and St than the S-A model.

141

 

0.0040

 

0.0035 1 x
0.0030 ~

0.0025 -

 

+ x Cf
0.0020 1 + +
T + + St

 

 

 

0.0015 «

0.0010 ~

0.0005 r

 

 

0.0000

 

S-A steady DES1 DESZ Em Correlation

Figure 7-12. Comparison of CFD results for flat panels using S-A and DES
approaches with experimental and correlation data.

7.4.2 Turbulent Flow over Rough Surfaces in a Wind Tunnel

A set of studies using the flat panel indicated that the computational grid should
have a y+ value near or less than unity for the cells near a no-slip wall. Since the smaller
the y+ value, the more CPU time is needed to achieve solution convergence. For all the
simulations with roughness panels, the average y+ value of the first grid layer near the
wall is near unity.

As mentioned earlier, two rough surfaces (Surface #4 and #6 from Reference
[45]) were employed in the present study. In addition, two different grids were generated
for each surface. For example, for Surface #4, the coarse grid has 364,005 total cells with
108,709 hexahedrons, 234,952 prism cells and 20,344 polyhedral cells, while for the fine
grid it has a total of 1,260,051 cells with 400,662 hexahedrons, 782,340 prisms, and
77,049 polyhedral cells. For Surface #6, the coarse and fine grids have 873,221 and

1,601,430 cells respectively. The fine grid essentially doubles the grid resolution near the

142

.—

 

 

roughness panels while maintaining an average y+ of 1 in the wall normal direction. The
coarse mesh has 64 cells in the span-wise direction with a grid resolution of about 1 mm,
while the fine mesh has a grid resolution of about 0.5 mm. The fine grids have about 50-
80 layers in the tunnel height direction. If structured grids were used for this
configuration, the fine grid then would have about 25-40 million cells. Using the viscous
adaptive Cartesian grid approach, the number of grid cells can be reduced by over an

order of magnitude. All the grids look similar to those shown in Figure 7-8 and Figure

 

 

 

 

 

 

 

 

 

7-9.
1.000
—Cf
0.100 ~
0, - ----- St
0
..l
0.010 L
0.001 1 . 1 .
0 2000 4000 6000 8000 10000

Iterations

Figure 7-13. Convergence histories using the S-A model with the fine mesh for Surface
#4.

All the simulations using the S-A model were carried out with a local time
stepping strategy and all the DES runs were in time accurate mode with a 2nd-order
backward difference formula. The flow convergence is monitored by the history of the
average of and St over the roughness panels. For example, Figure 7-13 shows the
convergence history of the average Cf and St using the S-A model on the fine grid for

surface #4. Note that the solution is converged after a few thousands of iterations.

143

 

 

 

Table 7-1 summarizes the computational results with comparison to
experimentally measured Cf and St. The computed average Cf and St on the coarse and
fine meshes still show large discrepancies, indicating the coarse grid is just too coarse. In
order to demonstrate grid independency, an even finer grid is necessary. On the fine grid,
the results are more encouraging. With both rough surfaces, the computed Cf number is
within 3.5% of the experimental data. The difference between the computed and
measured St number is larger at about 15% for Surface #4 and 8% for Surface #6 on the
fine grids. This difference may be due to several factors. One factor is that a constant
wall temperature was used in the computation. In the actual experiment, the wall
temperature is not constant. Another factor is insufficient grid resolution. It appears that
St increases with grid refinement. Further investigation with finer grids and non-constant

wall temperature will be carried out to find the reason. The computed Cf and St using

DES are too small, and the cause will be investigated.

Table 7-1. Comparison of experimental and computational results for rough

 

 

 

 

 

Surfaces # 4 and #6.
Surface Surface Surface Surface Surface Surface Surface
#4 #4 #4 #4 #6 #6 #6
(Coarse) (Fine) (Exp.) (DES ) (Coarse) (Fine) (Exp.)
Cf 0.0128 0.00970 0.00937 0.00329 0.0113 0.0100 0.0103
St 0.00255 0.00260 0.00308 0.00061 0.00259 0.00284 0.00308

 

 

 

 

 

 

 

 

 

the experimental data in Figure 7-14 and Figure 7-15 at Rex

correlation formulas are given below.

144

The computational results are also compared with four roughness correlations and

= 900,000. These

 

 

cf = [1.4 + 3.7 log(x / k, )1"2 from White[60]

cf = [2.87 + l.5810g(x / k, )1‘25 from Schlichting[6l]
cf = [3.476 + 0.707 ln(x / 1191446 from Mills[62]
cf = 0.168 [111(845/11, )1‘2 from Kays and Crawford[63]

The k, value in these correlations were computed based on A, which was

tabulated in Table l of Ref. [45].

The dashed line in Figure 7-14 is the of number of the smooth panel as a

 

reference. As shown in Figure 7-14, the computed Cf matches the experimental data very
well. The Schlichting correlation also gives good prediction. All the correlations appear

to bound the experimental and CFD data. This may suggest that CFD can be used as an

 

effective tool to predict Cf for “real” rough surfaces.

 

 

 

 

 

 

 

0.012 . Er
El 0
O
0.01 1 ¥ . 1
A
0.008 ~ A
Cf 0.006 ~
—+—— CFD Cf
+ Exp Cf
”004 ‘ ,1 .......... + .......... ,. .......... ,. ——-B—White (ks)
x 1 Schlichting(ks)
+Mills (ks)
0.002 ‘ —e—Kays (ks)
-----4— Smooth-Reference
I
0 1 fl l

 

Surface1 Surfacez Surface3 Surface4 Surface5 Surface6
Figure 7-14. Comparison of Skin frictions for roughness panels.

145

The St results are compared in Figure 12. In this graph, three correlations are used

 

 

for comparison
0.50f
S , = 0 2 from Kays and Crawford [63]
Pr,+,/0.5cf (ki' Prm/ C)
0.5cf .
S , = 0 2 from Drpprey and Sabersky [64]
1+ 1/0.5cf (5.19ki‘ Pro'“—8.5)
0.5cf .
S from Wassel and Mrlls [65]

 

' Pr, +,/0.5c, (4811:“2 Pr°-‘“—7.65)

In all the correlation formulas, the Cf value predicted with the Schlichiting
correlation is selected as the reference since it has the best agreement with experimental
and CFD data. As shown in Figure 7-15, CFD gives the lowest St prediction. Comparing

with the experimental data, all correlations predict higher rough surface St numbers.

 

 

 

 

0.005
E]
0.0045 4 CI
0.004 1
0 0035 A A
0.003 4 0 0
+
St 0.0025 « +
+ ---------- + ---------- + ---------- + ---------- 4» --------- -+
0.002 ~ ——+—CFD St
—0— Exp St
0.0015 —
—B— Kays&Crfrd (ks) [C=1]
0.001 1 x Kays&Crfrd (ks) [C=0.35]
—-A— Mills(ks)
0.0005 1 —e— Dipprey 81 Sabersky (ks)
0 - --+—- - Smooth-Reference

 

 

 

 

Surface1 Surface2 Surface3 Surface4 Surfaces Surface6

Figure 7-15. Comparison of Stanton numbers for roughness panels.

146

 

 

Since CFD can provide more detailed flow field data, a few “flow pictures” are
shown here to give the reader some ideas on the flow characteristics. A velocity vector
plot showing the flow near the rough panel is displayed in Figure 7-16. It is observed
that very complex separated flow regions exist near the rough surfaces. The surface
pressure distribution near the rough panel of Surface #6 on the fine grid is shown in

Figure 7-17. Clearly on the rough surfaces, the pressure drag is a dominant force in the

 

overall drag. In fact, over 75% of the total drag is due to pressure in both cases.

 

Figure 7-16. Velocity vector plot on a cutting plane near the roughness panel.

147

 

Figure 7-17. Pressure distribution on the rough surface (Surface #6).

7.5 Conclusions

In the present study, the skin friction (Cf) and heat transfer (St) coefficients on real
rough surfaces were directly computed using a state-of-the-art unstructured adaptive grid-
based finite volume flow solver. Both RANS and DES approaches were employed for
the computations. Computational results were compared to experimental data to assess
the simulation accuracy. Based on the present study, the following conclusions can be
drawn:

The unstructured adaptive grid generation method is very efficient in resolving
disparate length scales. It is estimated that the number of cells generated is more than an
order of magnitude less than that with a structured grid. The rough surfaces can be

handled by the grid generator with minimum user interference.

148

 

In the flat panel case, both the S-A and DES approaches are capable of predicting
Cf and St, thus validating the implementation.

With proper grid resolution, the S-A model is able to accurately predict Cf for both
rough surfaces (within 3.5% of experimental data). The computational predictions for St
showed 8-15% differences from the experimental data. This could be attributed to the
constant wall temperature used in the computational simulation, or insufficient grid
resolutions. Further investigation will be carried out to understand the reason for the
discrepancy.

The DES approach failed to predict either cf or St for both rough surfaces. The

reason will be investigated in the future.

7.6 Future Works

In this study, we performed CFD simulations on a flat plate by both S-A
turbulence model and DES approach. We have also finish simulation for two roughness
panels with surface #4 and surface #6 by S-A model. Surface #4 is characterized in
deposit and surface #6 is characterized in erosion. We plan to continue our study by
finishing the other four roughness panels. They are characterized by pitted surface and
spalled surface. So we can verify if CFD can predict the Cf and St precisely for all typical
turbine roughness. After that, DES turbulence model will be applied to these six
roughness surface to study the effects of turbulence modeling.

In addition, finer computational grids will be used to firrther assess grid

independence for all the cases with rough surfaces.

149

 

Appendix Summary of the Sparlart-Allmaras Model

The Spalart-Allmaras model is a one-equation model that solves a modeled
transport equation for the kinematic eddy (turbulent) viscosity. It was designed
specifically for aerospace applications involving wall-bounded flows and has been shown
to give good results for boundary layers subjected to adverse pressure gradients. It is also
popular for turbomachinery applications.

By Spalart-Allmaras model, the Reynolds-averaged Navier-Stokers equations and

a transport equation for the turbulence model are solved. The Reynolds stresses are given

by — uiuj = 2v,S,-j. The kinematic turbulent viscosity v, is given by

~ 1 _
V1=vara fvl =“3—r3—1 Z=—

where v is the molecular viscosity, fvl is the viscous damping function, i7 obeys

the transport equation as follows:

D17 ~~ 1 8 ~ 3‘7
_=cb1[1—f,2]Sv+—[—{(V+V) }+Cb2(
Dt 0' axj 6x]-

c i7
— [Cwlfw — 'Lzl’fIZ “2142 + ftlAUz
K

a

an

)2]

The term cb1[1— £21397 is the production term, in this term,

~

~ v
S=S+m3fvza fvzzl— Z

1'i’xiffvl

 

where S is the scalar measure of the deformation tensor. According to the
original model proposed by Sparlart and Alhnaras, S is based on the magnitude of the
vorticity, and d is the distance to the closest wall. The magnitude of the vorticity is

defined as follows:

150

While Qij is the mean rate-of—rotation tensor and it is defined by

Q. .1. %_5i

“:2 6x,- 6x]-

~

v . . . . . .
The term Cwlfw (5)2 rs the wall destructron term, 1n thrs term, the functron fw rs

 

defined as:
1+ 6 l "’
fw=8—_6 CW; 6 whereg=r+cw2(r6—r),ra~: 2
g +Cw3 Sk d

The term C_bzl fag—)2 is the wall production term, the ftz function is defined as

follows:
7.2 = 43 “pl—c.4272)

The trip function f,l is defined as follows:

 

2
w
ft] = ctlgt “IX-€12 Alt/2 [‘12 + 8:251:21)

where d t is the distance from the field point to the trip, which is on a wall, w, is

the wall vorticity at the trip, and AU is the difference between the velocity at the field

point and that at the trip. Then, g, E min(0.l,AU / thx,) where A x, is the grid spacing

along the wall at the trip.

151

The constants used above are Cb1=0.1335, 0:2/ 3 , Cbz =0.622 , K=O.41,
cw1=cb1/K2+(l+cb2)/0', sz =0.3, CW3 =2, Cv1=7.1, ct] =1, 6,2 =2, c,3 =1.2,

6’4 = 0.5 .

In our real roughness study, we don’t consider the trip function, so the SA
transport equation is simplified as follows:

D17 ~~ l 6 ~ 617
—— = 0615" +-[—-{(V+ V)--} +Cbz(
Dt 0' axj ax]

9::

6x-

2 _ 22
1)] Cwlfwid]

While Fluent company use a little bit different definition about the S , for details,

please refer to Fluent manual [22].

152

 

 

 

[1]

[2]

[3]

[4]

[5]

[61

[7]

[8]

[9]

[10]

[11]

Reference

Addy, H.E., “Ice Accretions and Icing Effects for Modern Airfoils,”
NASA/TP-2000-210031, April 2000.

Chung, J., Choo, Y.K., Reehorst, A., Potapczuk, and Slater, J., “Navier-
Stokes Analysis of the Flowfield Characteristics of an Ice Contaminated
Aircraft Wing,” AIAA 99-0375.

Shim, J., Chung, J., and Lee, K.D., “A Computational Investigation of Ice
Geometry Effects on Airfoil Performances,” AIAA Paper 2001-0540, Jan.
2001.

Chung, J.J. and Addy, HE, “A Numerical Evaluation of Icing Effects on a
Natural Laminar Flow Airfoil,” AIAA Paper 2000-0096, Jan. 2000.

Choo, Y.K., Slater, J.W., Henderson, T., Bidwell, C., Braun, D., and Chung
J.J., “User Manual for Beta Version of Turbo-GRD: A Software System for
Interactive Two-Dimensional Boundary/Field Grid Generation,
Modification, and Refinement,” NASA TM-1998-206631, 1998.

Chung, J., Reehorst, A.L., Choo, Y.K., and Potapczuk, “Effects of Airfoil
Ice Shape Smoothing on the aerodynamic Performance,” AIAA Paper 98-
3242, July 1998.

Vickerman, M., Choo, Y., Braun, D., Baez, M., and Gnepp, S., “SmaggIce:
Surface Modeling and Grid Generation for Iced Airfoils, Phase 1 Results,”
AIAA Paper 2000-0235, Jan. 2000.

Thompson, 0.8. and Soni, B., “ICEGZD (v2.0) — An Integrated Software
Package for Automated Prediction of Flow Fields for Single-Element
Airfoils with Ice Accretion,”NASA/CR-2001-210965, June 2001.

Shih, T.I-P., Bailey, R.T., Nguyen, H.L., and Roelke, R.J., “Algebraic Grid
Generation for Complex Geometries,” International Journal for Numerical
Methods in Fluids, Vol. 13, 1991, pp. 1-31.

Steinthorsson, E., Shih, T.I-P., and Roelke, R.J., “Enhancing Control of
Grid Distribution in Algebraic Grid Generation,” International Journal for
Numerical Methods in Fluids, Vol. 15, 1992, pp. 297-311.

Tai, T.C., “Single-Block Structured Body-Conforming Grid for Complex
Geometries,” Numerical Grid Generation in Computational Field
Simulations, edited by B.K. Soni, J.F. Thompson, J. Hauser, and P. R.
Eiseman, ISGG, Mississippi State University, 2000, pp. 91-100.

153

 

 

 

[12]

[13]

[14]

[15]

[15]

[17]

[18]

[19]

[20]

[21]

[22]

[2 3]

[24]

[25]

[26]

Power, G.D. and Underwood, M.L., ‘WIND 2.0: Progress on an
Applications-Oriented CFD Code,” AIAA Paper 99-3212, June 1999.

Michal, T. and Oser, M., “Improving Zonal Coupling Accuracy and
Robustness in the WIND Code,” AIAA Paper 2001 -0222, Jan. 2001.

Harned, A., Das, K., and Basu, 0., “Numerical Simulation of Unsteady
Flow in Resonance Tube,” AIAA Paper 2002-1118, January 2002

Addy, H., Chung, J., “A Wind Tunnel Study of Icing on a Naturally Laminar
Flow Airfoil,” AIAA Paper 2000-0095, January 200

Chi, X., Zhu, 8., Shih, T.I-P., Slater, J.W., Addy, HE, and Choo, Y.K.,
“Computing Aerodynamic Performance of 2D Iced Airfoils: Blocking
Topology and Grid Generation,” AIAA Paper 2002-0381, January 2002.

Zhu, B., Chi, X., Shih, T.I-P., and Slater, J.W., “Computing Aerodynamic
Performance of 2D Iced Airfoils: Blocking Strategies and Convergence

Zhu, 8., Chi, X., Shih, T.I-P., Slater, J.W., Addy, HE, and Choo, Y.,
“Computing Aerodynamic Performance of 2-D Iced Airfoils with Structured
Grids,” AIAA Paper 2003-1071, January 2003.

Shim, J., Chung, J., and Lee, K.D., “A Computational Investigation of Ice
Geometry Effects on Airfoil Performances,” AIAA Paper 2001-0540, Jan.
2001.

Broeren, A.P., Addy, HE, and Bragg. M.B., “Flowfield Measurements
about an Airfoil with Leading-Edge Ice Shapes,” AIAA Paper 2004-0559,
Jan.2004.

Michal, T. and Oser, M., “Improving Zonal Coupling Accuracy and
Robustness in the WIND Code,” AIAA Paper 2001-0222, January 2001.

FLUENT, Software Package, Ver. 6.1.18, Lebanon, New Hampshire,
03766, URL:http://www.fluent.com/software/fluentlindex.htm [cited 1
February 2005].

http://www.fluent.com/software/fluent/index.htm.

Spalart, PR. and AIImaras, S.R., “A One-Equation Turbulence Model for
Aerodynamic Flows,” AIAA Paper 1992-0439, January 1992.

Spalart, P.R., AIImaras, S.R., “A one-equation turbulence model for
aerodynamic flows,” AIAA-92-0439

Spalart, P.R., AIImaras, SR, 1994. “A one-equation turbulence model for
aerodynamic flows,” La Recherche A erospatiale 1, 5—21

154

[27]

[28]

[29]

I30]

[31]

[32]

[331

I34]

[35]

[36]

[37]

[33]

I39]

Menter, FR, 1991, “Performance of Popular Turbulence Models for
Attached and Separated Adverse Pressure Gradient Flows,” AIAA J., Vol.
30, No. 8. Pp. 2066-2071.

Menter, FR, 1993, “Zonal Two-Equation k-w Turbulence Models for
Aerodynamic Flows,” AIAA Paper 93-2906.

Launder, BE. and Spalding, 0.8., “The Numerical Computation of
Turbulent Flows,” Computer Methods in Applied Mechanics and
Engineering, Vol. 3, 1974, pp. 269-289.

Launder, B.E., “Second-Moment Closure: Present... and Future?”
International J. Heat Fluid Flow, Vol. 10(4), 1989, pp. 282-300.

Launder, BE. and Spalding, D.B., “The Numerical Computation of
Turbulent Flows,” Computer Methods in Applied Mechanics and
Engineering, Vol. 3, 1974, pp. 269-289.

Durbin, PA, “A Reynolds Stress Model for Near-Wall Turbulence,”
Journal of Fluid Mechanics. Pp. 465-498, Vol. 249, 1993.

Fu, 8., Launder, BE, and Leschziner, M.A., “Modeling Strongly Swirling
Recirculating Jet Flow with Reynolds-Stress Transport Closures,” Sixth
Symposium on Turbulent Shear Flows, Toulouse, France, 1987.

Gibson, M. M. and Launder, B. E., “Ground Effects on Pressure
Fluctuations in the Atmospheric Boundary Layer,” J. of Fluid Mechanics,
Vol. 86, 1978. pp. 491-511.

Chen, H.C. and Patel, V.C., “Near-Wall Turbulence Models for Complex
Flows Including Separation," AIAA J., Vol. 26, No. 6, 1988, pp. 641-648.

Shih, T.I-P., Bailey, R.T., Nguyen, H.L., and Roelke, R.J., “Algebraic Grid
Generation for Complex Geometries,” lntemational J. for Numerical
Methods in Fluids, Vol. 13, 1991, pp. 1-31.

Goldstein, R., Eckert, E., Chiang, H., and Elovic, E., “Effect of Surface
Roughness on Film Cooling Performance,” ASME Journal of Engineering
for Gas Turbines and Power, January 1985, Vol. 107, pp. 111-116.

Pinson, MW. and Wang, T., “Effects of Leading Edge Roughness on Fluid
Flow and Heat Transfer in the Transitional Boundary Layer over a Flat
Plate,” International Journal of Heat and Mass Transfer, 1997, Vol. 40, No.
12, pp. 2813-2823.

Taylor, R.P., Scaggs, W.F., and Coleman, H.W., 1988, “Measurement and
Prediction of the Effects of Non-uniforrn Surface Roughness on Turbulent

155

[40]

I41]

[42]

[43]

[44]

[4 5]

I46]

[47]

[43]

[49]

Flow Friction Coefficients,” ASME Journal of Fluids Engineering, Vol. 110,
Dec 1988. pp. 380-384.

Hosni, M.H., Coleman, H.W., and Taylor, RP, 1991, “Measurements and
Calculations of Rough-Wall Heat Transfer in the Turbulent Boundary
Layer,” lntemational Journal of Heat and Mass Transfer, Vol. 34, No. 4/5,
pp. 1067-1082.

Scaggs, W.F., Taylor, RR, and Coleman, H.W., “Measurement and
Prediction of Rough Wall Effects on Friction Factor - Uniform Roughness
Results,” ASME Journal of Fluids Engineering, Vol. 110, Dec 1988, pp.
385-391.

Bogard, D.G., Schmidt, D.L., and Tabbita, M., “Characterization and
Laboratory Simulation of Turbine Airfoil Surface Roughness and
Associated Heat Transfer,” Journal of Turbomachinery, Vol. 120, No. 2,
April 1998, pp. 337-342.

Barlow, ON. and Kim, Y.W., 1995, “Effect of Surface Roughness on Local
Heat Transfer and Film Cooling Effectiveness,” presented at the June,
1995 ASME lntemational Gas Turbine Exposition in Houston, Texas,
ASME paper #95-GT-14.

Bons, J.P., Taylor, R., McClain, S., Rivir, R.B., “The Many Faces of
Turbine Surface Roughness,” Journal of Turbomachinery, V0. 123, No. 4,
October 2001, pp. 739-748

Jeffrey P. Bons, “St and Cf augmentation for real turbine roughness with
elevated freestream turbulence,” Proceeding of TURBOEXPO 2002:
International gas turbine conference June 3-6, 2002 in Amsterdam, The
Netherlands.

Wang, Z.J. and Srinivasan, K., 2002, ”An Adaptive Cartesian Grid
Generation Method for 'Dirty’ Geometry," lntemational Journal for
Numerical Methods in Fluids, vol. 39, pp. 703-717.

Wang, Z.J. and Chen, R.F., 2002, “Anisotropic Solution-Adaptive Viscous
Cartesian Grid Method for Turbulent Flow Simulation,” AIAA Journal, vol.
40. pp. 1969-1978.

Spalart, P.R., Jou, W.-H., Strelets, M., AIImaras, SR, 1997. “Comments
on the feasibility of LES for wings, and on a hybrid RANS/LES approach,”
In: Liu, C., Liu, Z. (Eds), First AFOSR International Conference on
DNS/LES, 4-8 August, Ruston, LA, Advances in DNS/LES, Greyden
Press, Columbus, OH.

P. R. Sparlart, “Strategies for turbulence modeling and simulations,”
lntemational Journal of Heat and Fluid Flow 21(2000) p252-263

156

[50]

[51]

[52]

[53]

[54]

[55]

I55]

[57]

[581

I591

[60]

[61]

[62]

Spalart, “Trends in turbulence treatments,” AIAA Paper 2000-2306, Fluids
2000, Denver.

Nikuradse, J., 1933, “Laws for Flows in Rough Pipes,” VDI-Forchungsheft
361, series B, Vol. 4. (English Translation NACA TM 1292, 1950)

Acharya, M., Bomsterin, J., Escudier, M., 1986, “Turbulent Boundary
Layers on Rough Surfaces,” Experiments in Fluids, No 4, pp. 33—47

Schultz, BL. and Jones, T.V., 1973, “Heat-transfer Measurements in
Short-duration Hypersonic Facilities,” Advisory Group for Aerospace
Research and Developments, No. 165, NATO

Roe, P. L., “Approximate Reimann Solvers, Parameter Vectors, and
Difference Schemes,” Journal of Computational Physics, Vol. 43, 1983,
p.357.

Wang, Z.J., “A Quadtree-Based Adaptive Cartesian/Quad Grid Flow
Solver for Navier-Stokers Equations,” Computers and Fluids, vol. 27, No.
4, 1998. pp. 529-549.

Wang, Z.J., “Proof of Concept - An Automatic CFD Computing
Environment with a Cartesian/Prism Grid Generator, Grid Adaptor and
Flow Solver,” Proceedings of 4th lntemational Meshing Roundtable,
Sandia National Labs, Albuquerque, NM, 1995, pp.87-99.

Venkatakrishnan, V., “Convergence to Steady State Solutions of the Euler
Equations on Unstructured Grids with Limiters,” Journal of Computational
Physics,” Vol. 118, 1995, pp. 120-130.

Chen, R.F., and Wang, Z.J., “Fast, Block Lower-Upper Symmetric Gauss-
Seidel Scheme for Arbitrary Grids,” AIAA Journal, Vol. 38, No. 12, 2000,
pp. 2238-2245

Shur, M., Spalart, P.R., Strelets, M., Travin, A., 1999. “Detached-eddy
simulation of an airfoil at high angle of attack,” Proceedings of the Fourth
International Symposium on Engineering Turbulence Modeling and
Measurements, 24-26 May, Corsica, Elsevier, Amsterdam.

White, F.M., Viscous Fluid Flow, 2nd Edition, 1991, McGraw-Hill, New
York.

Schilichting, H. Boundary layer Theory, 7th Editon, 1979, McGraw-Hill,
New York

Mills, A.F., Heat transfer, 2st Edition, Prentice-Hall, New Jersey. ISBN 0-
13-947624-5 .

157

[63]

[54]

[6 5]

[56]

[671

[68]

[69]

[70]

[71]

[72]

[73]

[74]

[75]

[76]

Kays, W.M., Crawford, M.E., Convective Heat and Mass Transfer, 3rd
Edition, 1993, McGraw-Hill, New York

Dippery, BF, and Sabersky, R.H, “Heat and Momentum Transfer in
Smooth and Rough Tubes at Various Prandal Numbers,” International
Journal of Heat and Mass Transfer, Vol. 6, pp.329-353

Wassel, A.T. and Mills, A.F., “Calculating of Variable Property Turbulent
friction and Heat Transfer in Rough Pipes,” Journal of Heat Transfer, Vol.
101 August 1979. pp.469-474

Addy, H.E., “Ice Accretions and Icing Effects for Modern Airfoils,”
NASA/TP-2000-210031, April 2000.

Mani, M., “A Structured and Hybrid- unstructured Grid Euler and Navier-
Stokes Solver for General Geometry,” AIAA Paper 2004-0524, Jan. 2004.

Nelson, C., Lankford, D., Nichols, R., “Recent Improvements to the
WIND(-US) Code at AEDC,” AIAA Paper 2004-0527, Jan. 2004.

Chi, X., Zhu, 8., Shih, T.I-P., Addy, HE, and Choo, Y.K., “CFD Analysis
of the Aerodynamics of a Business-Jet Airfoil with Leading-Edge Ice
Accretion,” AIAA Paper 2004-0560, Aerospace Sciences Meeting, Reno,
January 2004

Durbin, P.A., “Separated Flow Computations with the k-s-v2 Model,” AIAA
Journal, pp. 659-664, Vol. 33(4), 1995

Yakhot, V. and Orszag S. A., “Renorrnalization group analysis of
turbulence. I. Basic theory,” J. Sci. Comput., Vol. 1, pp. 3-51, 1986.

Smith, L. M. and Woodruff S. L., “Renorrnalization-group analysis of
turbulence”, Annu. Rev. Fluid. Mech. Vol. 30, pp. 275-310, 1998.

PowerFLOW, Software Package, Ver. 3.4, Burlington, MA, 01803,
URL:http://www.exa.com[cited 1 February 2005]

Li, Y., Shock, R., Zhang, R., Chen, H., and Shih, T.I-P. “Simulation of Flow
over an Iced Airfoil by Using a Lattice-Boltzmann Method,” AIAA Paper
2005-1103, January 2005.

Chen, S. and DooIen, G. D.,“Lattice Boltzmann Methods for Fluid Flows,”
“Lattice Boltzmann Method for Fluid Flows,” Ann. Rev. Fluid Mech., Vol.
30, 1997. PP. 329-364.

Chen, H., Kandasamy, S., Orszag, 8., Shock, R. , Succi, S., and Yakhot,

V., “Extended Boltzmann kinetic Equation for turbulent flows.” Science
Vol. 301, 2003. pp. 633-636.

158

[77]

[78]

[79]

[801

[81]

Chen, H., Orszag S. Staroselsky I., and Succi S., “ Expanded Analogy
between Boltzmann Kinetic Theory of Fluid and Turbulence”, J. Fluid
Mech., Vol. 519, 2004, pp. 307-314.

Chen, H., Teixeira, C., and Molvig, K. “Digital physics approach to
computational fluid dynamics: some basic theoretical features” Int. J. Mod.
Phys. C Vol. 8, 1997, pp. 675-684.

Jayatilleke, C., “The lnfuence of Prandtl Number and Surface Roughness
on the Resistance of the Laminar Sublayer to Momentum and Heat
Transfer,” Prog. Heat Mass Transfer, Vol. 1, 1969, pp. 193-321.

Kim, SE. and Choudhury, D., “A Near-Wall Treatment Using Wall
Functions Sensitized to Pressure Gradient,” ASME FED Vol. 217,
Separated and Complex Flows, ASME, 1995.

Viegas, J.R., Rubesin, M.W., and Horstman, 0.0., “On the Use of Wall
Functions as Boundary Conditions for Two-Dimensional Separated
Compressible Flows,” AIAA-85-O180, AIAA 23rd Aerospace Sciences
Meeting, Reno, Nevada, 1985.

159

llllllllllllllllllllllllIllIIIIIIIIlllllllllllllll