The primary goals are directed toward the development of a numerical method for computing flow about ice accretion shapes and determining the influence of these shapes on flow degradation. It is expedient to investigate various aspects of icing independently in order to assess their contribution to the overall icing phenomena. The specific aspects to be examined include the water droplet trajectories with collection efficiencies and phase change on the surface, the flowfield about specified shapes including lift, drag, and heat transfer distribution, and surface roughness effects. The configurations computed were models of ice accretion shapes formed on a circular cylinder in the NASA Lewis Icing Research Tunnel. An existing Navier-Stokes program was modified to compute the flowfield over four shapes (2, 5, and 15 minute models of glaze ice, and a 15 minute accumulation of rime ice)

Scott, J. N.

English

NASA Technical Reports Server

NUMERICAL SIMULATION OF THE FLOWFIELD
OVER ICE ACCRETION SHAPES
SEMIANNUAL PROGRESS REPORT
30 October 1985 - 30 April 1986
Prepared For:
NASA Lewis Research Center
21000 Brookpark Road
Cleveland, Ohio 44135
Project Monitor: Dr. Robert J. Shaw
NASA Grant No. NAG 3-665
Submitted By:
/ James N/ Scott, Ph.D.
Principal Investigator
University of Dayton
Research Institute
300 College Park
Dayton, Ohio 45469
(N&SA-CR-176960V N U M E E I C A L SIMOLATION OF
THE FIOWFIELD O V E R ICE iCCBElJOU SHAPES
Semiannual Progress Report, 20 Apr. .1985 -
30 Apr. 1986 (Daytcn Oniv. , Ohic.). 19- P
CSCL 20D G3/34
N86-30093
Unclas
43244
https://ntrs.nasa.gov/search.jsp?R=19860020621 2020-03-20T14:40:51+00:00Z
NUMERICAL SIMULATION OF THE FLOWFIELD
OVER ICE ACCRETION SHAPES
1. INTRODUCTION
This is the first Semi-Annual Status Report submitted on
Grant NAG 3-665. It includes the progress made during the period
of October 1985 to April 1986. The NASA Technical Officer for this
grant is Dr. Robert J. Shaw, NASA Lewis Research Center, Cleveland,
Ohio.
2. R&D STATUS REPORT
The grant program schedule consists of five components, i.e.,
computer code revision, grid generation, computations, extraction
of data, and report preparation. During this reporting period the
first two components were accomplished and progress made on the
computations. In addition, an abstract was prepared and submitted
to the AIAA for potential presentation of these results at the
Aerospace Sciences Meeting, January 12-15, 1987. Details of the
program status follow in the next sections.
2.1 Technical Effort
In this section the progress of the technical effort
will be discussed. The primary goals of this program are directed
toward the development of a numerical method for computing flow
about ice accretion shapes and determining the influence of these
shapes on flow degradation. In pursuing these goals, it is
expedient to investigate various aspects of icing independently in
order to assess their contribution to the overall icing phenomena.
The specific aspects to be examined include the water
droplet trajectories with collection efficiencies and phase change
on the surface, the flowfield about specified shapes including
lift, drag, and heat transfer distribution, and surface roughness
effects. In treating these issues it must be kept in mind that
they will ultimately be coupled together in a fashion which will
permit accurate prediction of the complete icing process as well as
the influence on the surrounding flowfield. The following
paragraphs describe the progress in examining these issues.
The first issue to address is the unsteady nature of the
phenomenon. Although the ice shape configuration changes with
_L|
time, the rate of growth is slow (10 ft/sec) compared to the
flight speed (200 ft/sec), therefore, a quasi-steady analysis is
appropriate. This approach involves computation of Navier-Stokes
solutions for a sequence of shapes starting from the original
configuration, computing the growth rate, advancing the
configuration over a fixed time interval (e.g., 50 sec),
recomputing a new growth rate, etc., until the final period is
attained.
A primary key to this approach is the accurate computa-
tion of the heat transfer over a roughened surface of unusual
geometry. Although the Navier-Stokes equations conceptually
2
possess the ability to simulate this process, numerical
confirmation has not yet been achieved. Fortunately, a series of
experimental heat transfer tests have been accomplished
which assist in the validation of the computational fluid dynamic
(CFD) methods. The purpose of this phase of the investigation is
the computation of a series of these test cases to validate the
numerical simulation.
The configurations computed were models of ice accretion
shapes formed on a circular cylinder in the NASA Lewis Icing
Research Tunnel (IRT). These shapes were 2-, 5-, and 15-minute
models of glaze ice and a 15-minute accumulation of rime ice. An
existing Navier-Stokes program was modified to compute the
flowfield over these four shapes.
2.1.1 Governing equations
The governing equations are obtained by adding
body forces to the Navier-Stokes equations to account for the drag
of the droplet particles and surface roughness. In addition, a
continuity and momentum equation are required to develop the
trajectory equations for the droplets.
Air Equations
|£ + V.p_V = 0
DV
g)ppf-Y
Water Droplet Equations
3p
Ft K K
!|£* yp.v yp . g + f
where
C pS
£ =
 2m |V - V | (V - V ) = droplet drag/unit mass
The water droplets are assumed to be uniform spheres of diameter,
d. Cloud physics experiments indicate d to vary between 5 and 50
microns depending upon the air temperature and liquid water content
(see Aircraft Icing, NASA CP2086). Often an average value of 20
microns is used in ice accretion studies. Therefore,
and
2m
 2PW
CD = CD(Red) (see Figure 1)
ORIGINAL .
OF POOR QUALITY
,;rm4U-ai
,;;! \ Ml' : ' :
til'
I I
j .
e \
I
M-~H
r >•' - »" -".n.9 -' - 6.1... .' - 6 .•!,, , 2 v tf d .
Figure 1. Drag Coefficient for Spheres as a Function
of the Reynolds Number.
A curve fit of this classic plot produces the
following:
oh
Red < 1 CD = p£- (Stokes flow)
1 < Re, < 400
d -D
 (Ra<i)0.6«
400 < Red < 3 x 105 CD = 0.5
These values will be used in the computation of the water droplet
trajectories.
2.1.2 Boundary conditions
Far field. At the far field boundary, the
following conditions are prescribed:
V , T , P
— CD * 00 'O ' ' 00
p = LWC = liquid water content of cloud
V = V
- -oo
Surface . At the surface, the following
conditions prevail:
V = 0, T speci f ied = 0
~" """ W O II
2.1.3 Grid generation
The grid was developed using the hyperbolic grid
generator of Steger and Barth. This technique appears ideally
suited for computing grids for the irregular shapes of icing
configurations. This method first requires a detailed description
of the surface geometry. The program produces an orthogonal, body
oriented grid (which is preferable in CFD computations) and
clusters the grid points near the surface. The hyperbolic method
has little control over the final exterior grid shape, which is
perfectly acceptable for external flow problems (but unacceptable
for internal flows). An example grid is shown in Figure 2.
2.1.3.1 Solution of the water droplet
equations
An efficient means for solving the
water droplet trajectory equations was developed during this
reporting period. Droplet trajectory equations in the past have
been solved using the Lagrangian method. In this investigation it
was proposed to use the Eulerian method for two reasons. First to
make the droplet system of equations compatible with the airflow
equations, and secondly to include the variation of the liquid
water content (p ) throughout the flowfield. The Lagrangian method
presumes the water content to be constant.
The droplet equations are similar in
form to the airflow equations with the exception that the stress
tensor is zero. This difference, however, changes the mathematical
character of the partial differential equations from elliptic to
hyperbolic as will be shown in the following paragraph. This
observation means that the hyperbolic equations may be marched in
space using a Parabolized Navier-Stokes code (PNS). An advantage
results in that the droplet equations may be solved in only a few
seconds on a VAX computer.
2.1.3.2 Eigenvalue analysis of marching
technique
The purpose of this section is to
demonstrate the applicability of space-marching techniques to solve
the governing fluid dynamic equations for water droplets motion.
To correctly apply a space-marching technique the eigenvalues of
the governing equations must be real (indicating the character of
the equations is hyperbolic).
The governing equations in divergence
form are written as:
+ F = Hy
ORJQiNAl PAGE- SS
OE POOR QUALITY
I C E G R I D
Z/RAD
Figure 2. Grid for 15-Minute Glaze Case,
where E, F, and H are defined by
pu
2
pu
puv
F =
pv
puv
_PV 2_
H =
0~
- fx
_-fy_
In order to perform the eigenvalue
analysis the E and F vector are factored such that
f -
3 F
W
where
U
P
u
v
The A and B matrices are found to be
u
0
0
V
0
0
P
pu
0
0
pv
0
0
0
pu
P
0
pv
B =
Applying the above factorization can be rewritten as
U + A~1 B U = A 1H
-x - - -y
To determine the character of this equation set it is necessary to
find the eigenvalues of the matrix A B. A B was found to be
V
u
0
0
111
2
u
V
u
0
£_
u
0
1
u
-1
The characteristic equation of the
matrix A B is determined from
det[A-1B - XI]
or in simplified form
The eigenvalues of the matrix A B are defined as the roots of the
characteristic equation. Therefore, the eigenvalues of A B are
Since the eigenvalues are real, the character of the governing
fluid dynamic equations is hyperbolic. Therefore, the governing
equations can be solved with a hyperbolic space-marching scheme.
2.1.4 Solution procedure
The airflow was computed using a well documented
1 4Navier-Stokes code adapted from the original program developed by
Shang. The MacCormack algorithm is utilized and the program is
vectorized for the CRAY computer. Also, the hyperbolic grid
generator of Steger and Barth was used as discussed in the previous
section. Since this portion of the project is current state of the
art in CFD, little difficulty was encountered.
The computation of the droplet trajectories is
computed using technology developed for the PNS solving schemes.
The droplet equations are hyperbolic since the stress-tensor
vanishes for this model. (No collisions between droplets are
considered in the formulation and hence no pressure or shear terms
are present.) The governing equations are transformed into a £-n
computational domain and take on the following form:
pp
Pp
0
P
TJ
( V P -
u)
v)
TJ
-
VPU
where the contravariant velocity components are represented as
and
U = y u - x v
n p n P
v
The upstream boundary conditions are
Pp(«)
up(-)
v (.)
LWC
V
= 0
Given the air velocity components (u,v) from the solution of the
airflow equations, the droplet equations are marched in the in-
direction inward towards the body. After attaining the body
surface the values of u ,v , and p are recorded from which the
local collection efficiency (B) may be determined.
U COt 9
LWCW v w
10
where tan e = (HW = slope of the wall.
T W
A dimensional analysis of the water droplet
equations indicates that the collection efficiency is a function of
the geometry and only one flow parameter, T'; where
V T
1
 R
This single parameter is a type of Reynolds number, taking into
account all flow variables of the problem, and is very useful for
evaluating collection efficiencies. Figure 3 shows the sensitivity
of 3 to T' for various flows over circular cylinder.
2.2 Results
To calibrate the numerical procedure the heat transfer
distribution over a smooth circular cylinder under laminar condi-
tions was computed for which there exists an exact solution
(Frossling ). Figure 4 shows a comparison of Nusselt number over
a cylinder for the Navier-Stokes solution with the Frossling
solution. Conditions for this case are Re, = 138,000 and a
diameter of 2 inches. Excellent agreement may be observed, proving
the validity of the numerics.
Figure 5 shows the computed droplet velocity vector field
around a 2-inch diameter cylinder in potential flow with a free
stream velocity of 130 fps.
Preliminary computations of the flowfield about each of
the ice shapes have been performed. Refinement of these computa-
tions are currently being examined in order to achieve the required
accuracy.
11
oII
h
m
•
o
II
.3
Q)
c
•H
>H
O
O
ra
0)
cd
H
W
O
•H
05
to
cu
o
•H
O
•H
C
O
•H
-P
O
O
O
<u
-p
^
a
o
o
oo
bO
•H
AON3IDUJ3 aALLoanoo
12
8500
400
300
200
100
FR6SSLING
0 COMPUTED
R« - 138.666
Uoo - 136 FT/SEC
0 - 2.6 INCHES
©
10 20 30 40 50 60
THETA DEGREES
Figure 4. Comparison of Computed Nusselt Number
Distribution with Frossling Solution.
13
COMPUTED DROPLET VELOCITY VECTORS
ao-
>* 2.0-
1.0-
00
-3J5 - -2JO -IA -IJO -QA 0.0
X
Figure 5. Computed Droplet Velocity Vectors
Over a Cylinder T' = 0.^76.
References
1. Shaw, R. J., "Progress Toward the Development of an Aircraft
Icing Analysis Capability," NASA TM 83562, Jan. 1984.
2. Potapczuk, M. G. and Gerhart, P. M., "Progress in development
of a Navier-Stokes Solver for Evaluation of Iced Airfoil
Performance, AIAA Paper 85-0410, Jan. 1985.
3. Van Fossen, G. J. Simoneau, R. J., Olsen, W. A., and Shaw, R.
J. , "Heat .Transfer Distributions Around Nominal Ice Accretion
Shapes Formed on a Cylinder n the NASA Lewis Icing Research
Tunnel, AIAA paper 84-0017, Reno, Nevada, Jan. 1984.
4. Pais, M. and Singh, S. N., "Determination of the Local Heat
Transfer Coefficients of Three Simulated Smooth Ice Formation
Characteristics on a Cylinder," AIAA Paper, Snowmass,
Colorado, Aug. 1985.
5. Smith, M. E., Arimilli, R. V., and Keshock, E. G.,
"Measurement of Local Convective Heat Transfer Coefficients
of Four Ice Accretion Shapes," NASA Contract Report 174680,
May 1984.
6. Bragg, M. B. and Coirier, W. J., "Detailed Measurements of
the Flow Field in the Vicinity of an Airfoil with Glaze Ice,"
AIAA Paper 85-0409, Jan. 1985.
7. Langmuir, 1. and Blodgett, K. B., "A Mathematical Investiga-
tion of Water Droplet Trajectories," AAF-ATSC Tech. Report
5418, Feb. 1946.
8. Bragg, M. B. , "Rime Ice Accretion and Its Effect on Airfoil
Performance," NASA Contractor Report 165599, March 1982.
9. Bragg, M. B., Zaguli, R. J., and Gregorek, G. M., "Wind
Tunnel Evaluation of Airfoil Performance Using Simulated Ice
Shapes," NASA Contractor Report 167960, Nov. 1982.
10. Zaguli, R. G., Bragg, M. B., and Gregorek, G. M., "Results of
an Experimental Program Investigating the Effects of
Simulated Ice on the Performance of the NACA 63A415 Airfoil
with Flap," NASA Contractor Report 168188, Jan. 1984.
11. Messinger, B. L., "Equilibrium Temperature of an Unheated
Icing Surface as a Function of Airspeed," Journal of
Aeronautical Sciences, Vol. 20, Jan. 1953-
12. Hankey, W. L. and Kirchner, R., "Ice Accretion of Wing
Leading Edges," AFFDL TM-79-85-FXM, June 1979.
13- MacArthur, C. D., "Numerical Simulation of Airfoil Ice
Accretion," AIAA Paper 83-0112, Jan. 1983-
15
14. Scott, J. -N. and Hankey, W. L. , "Numerical Simulation of Cold
Flow in an,Axisymmeric Centerbody Combustor," AIAA Journal,
Vol. 23. No. 5, May 1985.
15. Frossling, N., NACA Tech Memo H32.
16
APPENDIX
Description and Use of HGRID
Hyperbolic Grid Generator
HGRID has a nicely structured modular form. Though it is not
completely clear to the uninformed observer just what each subroutine does,
the important routines are easily identifiable. Much of the code should be
considered a "black box". The routines which you will most likely want to
change are the main program, BODIS, and INITIA. These sections of the
program involve setting up the body, defining important parameters such as
loop sizes and output of the grid. A list of the primary subroutines and
their purpose is shown below in order of appearance in the code.
MAIN - Output of the completed grid.
SARC - Called at each arc by STEP, sets up integration parameters
based on input parameters and scaling factors.
BODIS - Reads in or generates the body coordinates, called by
INITIA.
INITIA - Reads in or assigns the smoothing and integration input
parameters, sets up stretching and scaling factors based
on the input parameters, called by MAIN.
METRIC - Generates the coordinate transformation metrics for the
integration, called at each arc by STEP.
STEP - Marches the grid generation out from the body. Calculates
X, Y coordinates at next arc, called by MAIN.
VOLUME - Calculates the "volume" or area of each cell formed by the
present arc and where the new arc should lie based on
variable scaling factors. Called by STEP.
The basic procedure in the process of grid generation is as follows:
1) Body coordinates and input parameters are read, and the
body coordinates are written to the output file.
2) STEP is called and it calls routines to set up:
integration parameters, cell areas and the right hand side
of the PDE to be solved.
3) The PDE is solved by routines called by STEP, and within
STEP the X, Y values for the next arc are calculated.
17
4) In MAIN, the X, Y values of the new arc are written to
tape, and STEP is called again.
The sequence of calculating the next arc out from the body and then
writing it to tape eliminates the need to store X and Y in two dimensional
arrays. This cuts down on the memory required by HGRID.
Changes to HGRID will most likely be changes in the body coordinates,
loop sizes and final grid form. These changes will all take place in MAIN,
BODIS and INITIA. The other input parameters defined in INITIA should be
considered constants that need not be changed. They have been found to work
well for many grid generation cases. Two factors, however, that you may wish
to control are DSETA and ESCAL. DSETA is the initial step size out from the
body in eta, the transformed coordinate normal to the body. ESCAL controls
the scale factor that is used to calculate the cell areas or volumes.
Hint: In the event that only the upper or lower half of a grid is needed
due to symmetry conditions, you must supply HGRID with the whole set of body
coordinates anyway. For HGRID to function, it needs to be given a complete,
closed body. You can cut the mesh in half after the complete grid has been
generated.
18


Numerical simulation of the flowfield over ice accretion shapes

Abstract

Similar works

Full text

Available Versions

NASA Technical Reports Server