NASA Contractor Report 198408

ICOMP-95-17; CMOTI'-95-3

//v-_3_/

J

_ _L: .5_

Y

Center for Modeling of Turbulence and

Transition: Research Briefs_1995

(NASA-CR-198408) CENTER FOR

MODELING CF TURBULENCE AND

TRANSITION: RESEARCH BRIEFS, 1995

(NASA. Lewis Research Center}

108 p

N96-12237

Unclas

G3134 0071422

October 1995

Prepared for

Lewis Research Center

Under Cooperative Agreement NCC3-370

National Aeronautics and

Space Administration
Center for Modeling of Turbulence and Transition i

Research Briefs -1995

Table of Contents

Turbulence Program for Propulsion Systems

Tsan-Hsing Shih .............................. 1

k - e Eddy Viscosity Turbulence Models: Development and Applications

Zhigang "fang ............................... 7

Turbulence Models for CFD Codes

Jiang Zhu .............................. 13

Eddy Diffusivity Models for Heat Transfer

Aamir Shabbir ............................. 33

Computation of Compressible Turbulent Benchmark Flows

William W. Liou ............................ 45

Modeling of Turbulent_ Reacting Flows by PDF Methods

Andrew T. Norris ............................ 57

Application of PDF to Compressible Turbulent Reactive Flows

Andrew T. Hsu ............................. 63

Turbulence Model Development for NPARC

J. Zhu, Z. Yang, and T.-H. Shih ...................... 73

Calculations of Turbulent Reacting Flows in Can Combustors

A.T. Norris and W.W. Liou ........................ 89

Assessment of Turbulence Models in Turbomachinery Flows ...... 95
Center for Modeling of Turbulence and Transition ii

Research Briefs -1995

Appendix A : Industry-Wide Workshop ................ 97

Appendix B : Organization - 1995 ................... 101

Appendix C : List of Member's Publications .............. 103
Center for Modeling of Turbulence and Transition

Research Briefs - 1995

Turbulence Program for Propulsion Systems

Tsan-Hsing Shih

1. Introduction

1.1 Background

• The center for modeling of turbulence and transition (CMOTT) at NASA Lewis

Research Center (LeRC) has been in existence for about four years. In the first

three years, its main activities were developing and validating turbulence and

combustion models for propulsion systems, in an effort to remove the deficiencies

of existing models. Three workshops on computational turbulence modeling were

held at LeRC (1991, 1993, 1994).

• A peer review of the turbulence modeling activities at LeRC was held in September, 1993. Seven peers from GE, P&W, RocketDyne, Cornell University, University of California at Berkeley and NASA Ames conducted the review. The

objective of the review was to assess the turbulence program at LeRC/CMOTT

and to suggest the future direction of turbulence modeling activities for propulsion systems. In September of 1994, a NASA-wide turbulence program peer

review was held at NASA Ames Research Center (ARC) by sixteen peers from

the aerospace industry, universities and other agencies.

1.2 Recommendations from the 1993 peer review

• "LeRC should spend substantial effort being responsive to industry's current

pressing perceived needs; this involves extensive discussion with industry during

every phase of model development, analysis of industry's problems, goal oriented

model development, evaluation of models relative to industry's intended application ... "

• "LeRC has an obligation not only to respond to industry's requests for help, but

to play an autonomous, independent leadership role in providing models of the

highest quality, ..., which can be employed not only by the aircraft gas turbine

and rocket industries but also by other industries ..."

• "In the present financial climate, industry does not have the resources to undertake model development and evaluation. LeRC's help in this regard via the

creation of its turbulence modeling effort, is, therefore, welcome from the industry's standpoint."

• "It is important to work with the industry to evaluate the models and rank-order

them by performance and cost in order to identify the most appropriate models

for particular situations."

• There are many other useful suggestions and comments including collaboration

with industry, joint programs, industry-wide workshops, etc.
2 Tsan-Hsing Shih

1.3 Recommendations from the 1994 NASA-wide turbulence peer review

• NASA needs a right well defined turbulence program in order to help industry

meet its needs.

• The turbulence program needs a certain level of support, coordination, linkage

to industry needs, technology transfer, management and leadership.

• The current progress is not reflective of real problems.

Based on the peers' comments and the current environment of NASA LeRC, we

have planned and carried out a turbulence program for propulsion systems. The

program is briefly explained in the following sections.

2. Program goals at CMOTT

• Develop reliable turbulence and bypass transition models and combustion models

for complex flows in propulsion systems.

• Integrate developed models into deliverable CFD tools for propulsion systems in

collaboration with researchers at NASA LeRC and industry partners.

3. Program approaches

• Develop turbulence and combustion subprograms or modules for CFD users.

• Develop collaboration program with NASA LeRC and industry partners to facilitate technology transfer.

• Conduct turbulence modeling research for propulsion systems:

One-point moment closures for non-reacting flows,

Scalar PDF method for turbulent reacting flows,

Validation of existing and newly developed models.

4. Development of turbulence and combustion subprograms (modules)

4.1 Objective

Build quick and efficient vehicles for turbulence and combustion technology transfer to CFD users in propulsion systems, for example, inlet/nozzle, turbomachinery

and combustion, etc.

4.2 The features of the turbulence module

• It contains various turbulence models from which users can choose the model

appropriate for the flows of interest.

• It is self-contained, i.e., it contains its own solver for the turbulence model equations.

• It can be easily linked to industry CFD codes.

4.3 The features of the PDF combustion module

• It can be easily coupled with any existing industry flow codes.

• It has a novel averaging scheme to reduce memory requirement.

• It contains a general chemistry package.

• There is a parallelized workstation version.
Turbulence Program .for Propulsion Systems 3

4.4 NPARC Turbulence module

NPARC is a CFD code widely used by NASA and industry. A turbulence subprogram has been developed for this code.

• The models already built-in at the present time are:

two mixing length models; the Chien k - e model; and the CMOTT k - 6 model.

• The models to be built-in are:

the CMOTT Reynolds stress algebraic equation model; a Reynolds stress transport equation model; and other models requested by CFD users.

• A robust numerical solver is built-in for the model equations.

4.5 V-STAGE turbomachinery turbulence module

V-STAGE is an advanced turbomachinery code developed at NASA LeRC.

• The built-in turbulence model is the CMOTT k - e model.

5. Collaboration program and technology transfer

5.1 Joint programs with NASA LeRC and industry

• Preliminary programs with engine companies and others have been initiated (GE,

P&W, RocketDyne, Naval Research Laboratories).

• Joint programs with the offices of Advanced Subsonic Technology (AST), High

Speed Research (HSR) and Numerical Propulsion System Simulation (NPSS)

have been developed as ongoing programs with which industry partners will also

be closely involved.

5.2 Industry-wide workshop

The industry-wide workshop will be a regular program held once every other year

to:

• release Lewis turbulence and combustion modules to CFD users; and

• discuss the needs of industry and the state of the art in engineering turbulence

modeling.

6. Models developed at CMOTT

Turbulence and combustion model development and validation are ongoing research activities. The following are a part of work done at CMOTT.

• Isotropic eddy viscosity models:

NASA TM 1056331, 1057682, 1062633, 1067214.

• Reynolds stress and scalar flux algebraic equation models:

NASA TM 1061165, 1065138, 1066447, IJNMF 8.

• Second moment transport equation models:

NASA TM 1053519, 1064691°, 10668111, 105954 z2.

• Multiple-scale models for compressible turbulent flows:

NASA TM 10607213.

• Bypass transition models:
4 Tsan-ttsing Shih

NWTF 14 .

• PDF models for turbulent reacting flows:

NASA TM 10661415, AIAA 16, AIAA 17.

PROGRAM SUMMARY

r TURBULENCE

MODEL. |

_, DEVELOPMENT)

C

M FTURSULE"CE

O "_| COMBUSTION

T \ MODULES

_(DIRECT APPUCATIONS_ ;.orvcotaUng parts

-_-[ IN PROPULSION l_1 Combustor problems

Ls_--s J LTu_mach"_

Reynolds stress models_--_

-'L.°, o.- Jr_, ]

• ! PROGRAMS i

,_NPARC code

Main industw's codes J

IN

ID

T

WORKSHOP] R
Turbulence Program for Propulsion Systems 5

References

1 Shih, T.-H. and Lumley, J. L., "Kolmogorov behavior of near-wall turbulence and

its application in turbulence modeling," NASA TM 105663 (1992), also, Comp.

Fluid Dyn. 1993, Vol. 1, p.43-56.

2 Yang, Z. and Shih, T. -H., "A new time scale based k - e Model for Near Wall

Turbulence," NASA TM 105768, also, AIAA J. 31, 1191-1198 (1993).

3 Yang, Z. and Shih, T. -H., "A Galliean and tensorial invariant k - e model for

near wall turbulence," NASA TM 106263, also, AIAA 93-3105.

4 Shih, T.-H., W.W. Liou, A. Shabbir, Z. Yang and J. Zhu, "A new k- e eddy

viscosity model for high Reynolds number turbulent flows -- model development

and validation,", NASA TM 106721, August, 1994, also, Computers Fluids Vol.

24, No.3 pp.227-238, 1995.

Shih, T.-H. and Lumley, J., "Remarks on turbulent constitutive relations," NASA

TM 106116 (1993), also, Mathl. Comput. Modeling Vol. 18, No. 2, pp. 9-16,

1993.

6 Shih, T.-H., J. Zhu and J. Lumley, "Modeling of wall-bounded complex flows and

free shear flows," NASA TM 106513 (1994).

7 Shih, T.-H., J. Zhu and J. Lumley, "A new Reynolds stress algebraic equation

model," NASA TM 106644 (1994), also, to appear in Comput. Methods Appl.

Mech. Eng.

s Zhu, J. and T.-H. Shih, "Computation of confined coflow jets with three turbulence models," Int. Journal for Numerical Methods in Fluids, Vol. 19, pp.

939-956 (1994).

9 Shih, T.-H. and Lumley, J.L., "Critical comparison of second order closures with

direct numerical simulation of homogeneous turbulence," NASA TM 105351, also,

AIAA J., 31,663-671 (1993)

10 Shih, T.-H., Shabbir, A. and Lumley, J.L., "Realizability in Second Moment

Turbulence Closures Revisited," NASA TM 106496 (1994).

11 Shih, T.-H. and Shabbir, A., "Methods of Ensuring Realizability for NonRealizable Second Order Closures," NASA TM 106681 (1994).

12 Shabbir, A. and Shih, T.-H., "Critical comparison of second order turbulence

models in homogeneous flows," NASA TM 105954 (1994).

13 Liou, W. W., Shih, T. H. and Duncan, B. S., "A multiple-scale model for compressible turbulent flows," NASA TM 106072 (1993), also, Phys. Fluids A, Vol.7

(3), 1995.

14 Yang, Z. and Shih, T.-H., "A k - e model for turbulent and transition boundary

layers," Near-Wall Turbulent Flows, R.M.C. So, C.G. Speziale and B.E. Launder

(Editors), 1993 Elsevier Science Publishors B.V.
6 Tsan-Hsing Shih

15 Noris A.T. and Hsu, A., "Comparison of PDF and Moment Closure Methods in

the Modeling of Turbulent Reacting Flows," NASA TM 106614 (1994)

16 Hsu, A. T., "Study of hydrogen diffusion flames using PDF turbulence model,"

AIAA paper 91-1780 (1991).

17 Hsu, A., Raju, M.S. and Norris, A.T., "Application of a PDF method to compressible turbulent reacting flows," AIAA-94-0781, 1994.
Center for Modeling of Turbulence and Transition

Research Briefs - 1995

7

k-, Eddy Viscosity Turbulence Models:

Development and Applications

Z. Yang

1. Motivation and Objective

The purpose of this research is to develop k - e eddy viscosity models for the

complex flows of engineering interest and to incorporate the resulting models into

general purpose computational fluid dynamics (CFD) codes used to calculate practical flows relevant to aerospace and aero-propulsion systems.

2. Work Accomplished

In the following, work on model development and model applications will be

reported briefly.

2.1 Model development

Because of the limitations of computer speed and storage, DNS (direct numerical

simulation) is limited to flows of moderate Reynolds numbers and simple geometries.

For the present and the foreseeable future, turbulence modeling is the only viable

approach for the calculation of complex turbulent flows of engineering interest. In

turbulence modeling, the k-e eddy viscosity model is the most widely used model

in engineering calculations. Many flows have been calculated by k - e models. A

set of model constants have emerged from these computational experiences, giving

rise to what is commonly referred to as the standard k - e model 1'2. The standard

k - e model can be directly used for free shear flows. For wall bounded flows, it

can be used in conjunction with wall functions which prescribe the flow field at the

equilibrium log layer instead of at the wall. However universal wall functions do

not exist in complex flows and it is thus necessary to develop a form of k - e model

equations which can be integrated down to the wall. Jones and Launder 3 were the

first to propose such a low Reynolds number k - e model for near wall turbulence,

which was then followed by a number of similar k - e models. A critical evaluation

of the pre-1985 models was made by Patel et al.4. More recently proposed models

are found in Lang and Shih 5.

The existing k - e eddy viscosity models for near wall turbulence suffer from some

well known deficiencies. (Some of the existing models may be free from one or two

of these deficiencies.) First, a near wall pseudo-dissipation rate was introduced to

remove the singularity in the dissipation rate equation at the wall. However, the

definition of the near wall pseudo-dissipation rate was quite arbitrary. Second, the

model constants were different from those of the standard k - e model, making

the near wall models less capable of handling flows containing both high Reynolds

number turbulence and near wall turbulence, which is often the case for a real flow

situation. Patel et al.4 listed as the first criterion the ability of the near wall models
8 Z. ]Tang

to be able to predict turbulent free shear flows. Third, the variable y+ was used

in the damping function ], of the eddy viscosity formula. Since the definition of

y+ involves u_, the friction velocity, any model containing y+ can not be used in

flows with separation and reattachment. In addition, y+ may not be well defined

for flows with complex geometries.

All the k - e models for near turbulence reviewed in references 4 and 5 were

constructed based on the standard k - e model, i.e., they took the standard k - e

model form as the base form and added on it modifications due to the near wall

effect. The standard k - e model form itself has the following deficiencies: 1) the

model is not realizable in that physically unrealistic Reynolds stresses could be

obtained from the model in certain flow situations; 2) contrary to the experimental

findings, the predicted spreading rate for a round-jet is larger than that for a planarjet, a phenomenon commonly known as the planar-jet/round-jet anomaly; 3) the

inadequate predicted response to pressure gradient of turbulent boundary layer

flows with adverse pressure gradients is also attributed to the model form of the

standard k - e model, as pointed out by Wilcox 6.

The purpose of our research in the area of eddy viscosity modeling is to propose

a k - e eddy viscosity model which is free from the deficiencies mentioned above. In

addition, we require the proposed model to be free from coordinate parameters in the

model formulation. While coordinate information, distance to the wall for example,

provides a convenient parameter to calibrate the model's damping functions, its

use makes the model not tensorially invariant and thus unable to be incorporated

into computational fluid dynamics (CFD) codes with unstructured grids. The final

product of the present research is aimed at a realizable, tensorially invariant, and

Galilean invariant k - e eddy viscosity model which has the capability for both free

shear flows and wall bounded flows.

The components for such a model have been developed in the past few years. A

variable C, formulation has been proposed, in which C, depends on the solutions of

the mean flow field and the turbulence field. The C_, formulation was derived based

on the realizability constrains, rapid distortion theory (RDT), and the invariants

theory. The Reynolds stresses from the k-e eddy viscosity model using the proposed

Cu formulation are always realizable in that no physically un-admissible Reynolds

stresses will ever be produced. The variations of C_, with flow situations are in the

right direction as well, i.e., the value of Cu is smaller than 0.09 in regions where it

should be. Such a variation of C_, with flow situations make the model calculation

numerically more robust, particularly when the model is used to calculate flows with

large strain rate variations, for example, the flow of shock-wave/turbulent boundary

layer interactions. The detail of the proposed C_, formulation and the corresponding

model performance for benchmark free shear flows and wall bounded flows can be

found in Shih et al.7

The dissipation rate equation in the standard k - e model has a singularity at the

wall, where the turbulent kinetic energy is equal to zero. To remove this singularity,

and yet without introducing an arbitrary pseudo-disspation rate, the k - e model

was reformulated using velocity-scale and time-scale instead of the traditionally
k - e Eddy Viscosity Turbulence Models: Development and Applications 9

used velocity-scale and length-scale. It was further argued that the turbulent time

scale should be bounded from below by the Kolmogorov time scale which is always

non-zero. With such a turbulence time scale, the resulting dissipation rate equation

remains singularity-free everywhere in the flow field. The resulting dissipation rate

equation is computationally very robust. The details of using Kolmogorov behaviors

to aid the k - e eddy viscosity models for near wall turbulence can be found in Yang

and Shih s and Shih and Lumley 9.

It is a challenging task to propose a k - e eddy viscosity model for near wall

turbulence which does not contain coordinate information. The difficulty lies in the

lack of proper parameter(s) to calibrate the near wall effect. In their paper on the

first k - e eddy viscosity model for near wall turbulence, Jones and Launder 3 used

Rt, the turbulence Reynolds number, to calibrate the near wall effect. Rt is a field

quantity, and does not depend on the coordinate system. Thus, a k - e model with

Rt is tensorially invariant. However, the performance of the resulting model is not

very satisfactory, partly due to the fact that the near wall region for R_ is confined

to y+ about 20 while the near region for the damping function f_ extends to the

lower end of the equilibrium log layer, which normally has a y+ value of about 60.

Other researchers have later used y, the distance to the wall, as a parameter to

calibrate the near wall effect. Models using y in the damping functions are found to

give a better performance due to the fact that y grows more gradually in the near

wall region compared with R_. However, any model with coordinate information

will not be tensorially invariant. In addition, the definition of wall distance will be

ambiguous in complex flows with complex geometries.

Recently, Yang and Shih 1° have introduced a new parameter R to calibrate the

k _2/sk and is related to the following near wall effect. R is defined as R = _-_ -- _--/-,

two important physical parameters: the turbulent Reynolds number and the time

scale ratio of the turbulence to the mean flow. R defined above is expressed in

terms of the local field variables and is thus coordinate independent. In the near

wall region, R is found to increase with y+ in a gradual and monotonic manner.

This property makes R an ideal candidate for constructing the damping function.

In Yang and Shih 1°, such a damping function f_ was constructed to model the near

wall effect. The proposed k - e model for near wall turbulence uses the standard

k - e model as the corresponding high Reynolds number model and the resulting

model was validated for turbulent channel flows and turbulent boundary layer flows

with different pressure gradients.

To adequately capture the effect of pressure gradient on turbulent boundary layers, inhomogeneity effect was introduced in the dissipation rate equation. Since the

dissipation rate represents the energy transfer from the large eddies to the small

eddies in a turbulent flow, it is expected that the inhomogeneity in the large eddies

would have an effect on the energy transfer rate. The previously used dissipation

rate model is a homogeneous model in the sense that other than the diffusion term,

the source terms remain the same for both homogeneous flows and inhomogeneous

flows. In the present study, the effect of flow inhomogeneity is modeled directly

as an extra production/destruction source term in the dissipation rate equation.
10 Z. Yah 9

This extra source term is required to vanish when the flow is homogeneous. In the

framework of k - e eddy viscosity model, and for the case of weak inhomogeneity, we

characterize the inhomogeneity by Vk and VS, which represent the inhomogeneity

of the turbulent field and the inhomogeneity of the mean field, respectively. Invariants theory was used to deduce the possible forms of the source term. The resulting

model constants in the source term were further calibrated against a turbulent

boundary layer flow with favorable pressure gradients (the Herring and Norbury

flow it) and a turbulent boundary layer flow with adverse pressure gradients (the

Bradshaw flow12.) The resulting model was found to work well for other boundary layer flows with different pressure gradients. The model description and flow

calculations can be found in Yang and Shih 13.

The round-jet/planar-jet spreading rate anomaly is a long standing issue for turbulence modeling community. This anomaly is not only observed for two equation

eddy viscosity models, it is also observed for second order closure models. This

suggests that the dissipation rate equation needs to be modified in order to resolve

this anomaly, as pointed out by other researchers. The present research confirms

this conjecture and finds that either the extra source term in the dissipation rate

equation due to the inhomogeneity effect reported in reference 13 or the new dissipation rate equation based on the dynamic equation for the fluctuating vorticity

reported in reference 7 can resolve this round-jet/planar-jet spreading rate anomaly.

The computational results are found in the corresponding references.

Now, we have in place all the essential components to construct a realizable, tensorially invariant, and Galilean invariant k- e eddy viscosity model which works well

for both free shear flows and wall bounded flows. The variable C u formulation will

ensure the realizability of the proposed model; the analysis of the Kolmogorov behaviors of turbulence will make the model singularity-free without introducing any

arbitrary pseudo-dissipation rate; a new parameter R has been found which possess

the needed properties to calibrate the damping functions for near wall flows; the

effects of pressure gradients on turbulent boundary layers are captured via an extra

source term in the dissipation rate equation representing the contribution from flow

inhomogeneity; the round-jet/planar-jet spreading rate anomaly can be resolved via

the changes in the model dissipation rate equation, either by the contribution of the

inhomogeneity effect or by the new dissipation rate equation based on the dynamic

equation for the fluctuating vorticity. Now, what needs to be done is to put all the

components together and test the resulting model against a large variety of flows

of practical interest. This is a topic of our current research.

2.2 Model applications

The purpose of engineering turbulence modeling research is to use the proposed

turbulence models to calculate flows of engineering interest. Thus, the ultimate test

for our research on turbulence models lies in the ability of the proposed models to

calculate flows in aerospace and aero-propulsion systems. To this end, it is necessary

to install the proposed models into a CFD code which serves as a numerical test-bed.

In order to isolate the effect of turbulence model, the numerical test-bed should be

a well established code with a well developed numerical scheme. The code should
k - e Eddy Viscosity Turbulence Models: Development and Applications 11

be able to handle a wide range of flow situations and should be widely used in the

aerospace and aero-propulsion community.

In the present study, NPARC is chosen as the numerical test-bed for model applications. NPARC is a general purpose computational fluid dynamics (CFD) code.

It has a large user group in aerospace and aero-propulsion community. The code

solves the Navier-Stokes equations in a conservation form in a general coordinate

system. It can handle complex geometries and different types of boundary conditions. A description of the code and its numerical scheme is given by Cooper and

Sirbaugh 14. Updates of code capabilities and a listing of bibliographies of using

the NPARC code to calculate flows of practical interests can be found in the current version of the NPARC User Guide 15. The NPARC code is actively supported

by NASA Lewis Research Center via the NPARC Alliance, a partnership formed

between NASA Lewis Research Center and Air Force Arnold Engineering Development Center. The NPARC users are represented by the NPARC Association

currently headed by Boeing Company and McDonnell Douglas Company.

Two-equation k - e eddy viscosity turbulence models are incorporated into the

NPARC code via a separate turbulence subprogram written by Dr. Zhu of CMOTT.

(For a description of the philosophy and the structure of the subprogram, please see

Zhu and Shih16.) Different turbulence models have been tested against typical flows

in propulsion systems, including, for example, the ejector nozzle flows, transonic

diffuser flows, and the boat-tail nozzle flows. Details of flow calculations using the

NPARC code and the turbulence subprogram axe going to be presented by Yang et

al. 17 in the upcoming 1995 AIAA Joint Propulsion Conference.

By installing the advanced turbulence models into the NPARC code, NPARC's

turbulence modeling capability is enhanced. Turbulence model development for

NPARC is CMOTT's long term activity in supporting NASA's High Speed Research

(HSR) Program. Research in this area, together with model applications using the

NPARC code, is described in more detail in this research brief in the article entitled

"Turbulence Model Development for NPARC."

3. Future Plans

1) Turbulence modeling: development and applications

We will first finish the research on a realizable, tensorially invaxiant, and Galilean

invaxiant k - e eddy viscosity model. The model is going to be validated for both

simple shear flows, such as free shear flows and boundary layer flows, and the complex flows in aero-propulsion systems. Then, we will proceed with the development

of Reynolds stress models for wall bounded turbulent flows. The final product in the

latter area of research is a realizable, tensorially invariant, and Galilean invariant

second order closure model.

2) Transition modeling

We will be working on the modeling and calculation of transitional flows over a

low-pressure (LP) turbine blade. This research will be part of NASA's Low Pressure

Turbine Flow Physics Program. For such flows, the transition is induced by the high

turbulence level associated with the free-stream and the incoming wake, giving rise
12 Z. Yang

to what is referred to as the bypass transition. The final product of this research

will be a two equation model for transitional flow of bypass transition type.

4. References

1 Launder, B.E. and Spalding, D.B., "The numerical computation of turbulent

flow," Computer Methods in Applied Mechanics and Engineering 3, 269 (1974).

Rodi, W., Turbulence Models and Their Application in Hydraulics, Book Pub. of

International Association for Hydraulic Research, Delft, the Netherlands, 1980.

a Jones, W.P., and Launder, B.E., "The calculation of low-Reynolds number phenomena with a two-equation model of turbulence," Intl. J. of Heat and Mass

Transfer 16, 1119 (1973).

4 Patel, V.C., Rodi, W., and Scheuerer, G., "Turbulence models for near-wall and

low-Reynolds number flows: a review," AIAA J. 23, 1308 (1985).

5 Lang, N.J. and Shih, T.H., "A critical comparison of two equation turbulence

models," NASA TM 105237 (1991).

6 Wilcox, D.C., "Application of low Reynolds number two-equation turbulence

models to high Reynolds number flows," in Near Wall Turb. Flows, Elsevier, 1993.

7 Shih, T.-H., Liou, W.W., Shabbir, A., Yang, Z. and Zhu, J., "A new k-e eddy

viscosity model for high Reynolds number turbulent flows - Model development

and validation," J. of Computers and Fluids, 24, no. 3, 1995.

s Yang, Z. and Shih, T.H., "A new time scale based k - e model for near wall

turbulence," AIAA J. 31, 1191 (1993).

9 Shih, T.-H., and Lumley, J.L., "Kolmogorov behavior of near-wall turbulence and

its application in turbulence modeling," Comp. Fluid Dyn. 1, pp.43-56 (1993).

lo yang, Z. and Shih, T.H., "A Galilean and tensorial invariant k - e model for near

wall turbulence," AIAA Paper 93-3105 (1993).

11 Herring, H.J., and Norbury, J.F., "Some experiments on equilibrium turbulent

boundary layers in favorable pressure gradients," J. Fluid Mech. 27, 541 (1967).

12 Bradshaw, P., "Equilibrium boundary layer in moderate positive pressure gradient," Computation of Turbulent Boundary Layers-1968 AFSOR-IFP Stanford

University, Coles, D. E. and Hirst, E. A. edt., Vol.2, pp.241-248.

la Yang, Z. and Shih, T.H., "Flow inhomogeneity and the dissipation rate equation,"

Accepted for presentation at the 5th Euro. Turb. Conf., July 1994, Siena, Italy.

14 Cooper, G.K. and Sirbaugh, J.R., PARC code: theory and usage, AEDC-TR-89-

15, Arnold Engineering Development Center, 1989.

15 NPARC Alliance Technical Team, A User's Guide to NPARC, Version 2.0, 1994.

16 Zhu, J. and Shih, T.H., "A turbulence module for the NPARC code," to be

presented at the AIAA Joint Propulsion Conference, San Diego, July (1995).

17 Yang, Z., Georgiadis, N., Zhu, J. and Shih, T.H., "Calculations of inlet/nozzle

flows using a new k-e model," to be presented at the AIAA Joint Propulsion

Conference, San Diego, July (1995).
Center ]or Modeling of Turbulence and Transition

Research Brie]s - 1995

13

Turbulence Modules for CFD Codes

Jiang Zhu

1. Motivation and Objectives

It has been long recognized that there is a gap between turbulence model developers and CFD users. The former mainly use simple flows to verify new modeling

concepts and evaluate the resulting models, while the latter are usually reluctant to

implement and test new, more advanced turbulence models unless and until they see

such models showing good performance for a wide range of complex flow situations.

Turbulence model equations also require special treatment to ensure numerical realizability such as the positiveness of the turbulent kinetic energy and its dissipation

rate.

In order to bridge this gap, we develop turbulence modules for industry's CFD

codes. The modules are written in a self-contained manner so that the user can use

any turbulence model in the modules without worrying about how it is implemented

and solved. The input to a module is mean flow variables, boundary and geometric

information which are to be provided by a mean flow solver, and the output of

the module is the turbulent diffusivity and relevant turbulent source terms which

are needed for the mean flow calculation. The interaction between the mean flow

solver and the turbulence module will give the final turbulent flow solution. With

the aid of the modules, we can also take the advantage of the well-established and

sophisticated CFD codes to test turbulence models for a variety of complex flows

which are intractable with the simple research codes.

2. Work Accomplished

2.1 The NPARC Code

The NPARC code has been used extensively by the U.S. aerospace community

to analyze propulsion flows. However until recently, only algebraic models such

as the Thomas and Baldwin-Lomax have been available in NPARC for turbulent

flow simulations. The Chien low Reynolds number k-e model with modifications

for compressibility has been available in the 2D version but was not successfully

installed in the 3D version. Another k-e model (the NPARC 1.0 k-e model) was

installed in both the 2D and 3D versons but has not provided desirable accuracy

and stability.

In order to improve the NPARC's capability to calculate turbulent flows, the

two-equation turbulence model in the NAPRC code (both the 2D and 3D versions)

was modified 1 so that the model is based on the low Reynolds number k-e model of

Chien and no longer on the NPARC 1.0 k-e model. Stability enhancements and a

new inflow boundary condition for the turbulent quantities were also added to the

k-e model. Comparisons of the NPARC solutions obtained using the previous and

new models with experimental data indicated that the Chien k-e model installed in

this work improves the capability of the NPARC to calculate propulsion flows.
14 J. Zhu

2.2 A Turbulence Module for the FAST2D Code

FAST2D is a computer program developed at the CFD group of Prof. Rodi 2

for calculating two-dimensional, incompressible, elliptic flows with complex boundaries. Here "two-dimensional" stands for plane or axisymmetric flows with or without swirling. FAST2D is a finite-volume procedure and uses non-staggered grids

and the cell-center arrangement, i.e., the flow variables are stored at the geometric

centers of control volumes, with grid lines forming the faces of the control volumes.

FAST2D is being used at CMOTT as a research code. The turbulence module,

termed as FAST2DTS-1, is similar in structure to FAST2D. To reduce numerical diffusion while maintaining necessary stability, the second-orde accurate HLPA

scheme 3 is used for the convective terms of the turbulent transport equations. This

scheme is implemented such as to blend two schemes by a parameter _, with limiting values £---0 for the (first-order) upwind and _-1 for the second-order and

bounded solution. Seven turbulence models have been built into the module: the

standard k-c, CMOTT k-c, k-w, CMOTT Reynolds stress algebraic equation model,

Launder-Sharma, Chien and Shih-Lumley, with the first four being high- and the

last three being low-Reynolds number models. The module includes the four most

commonly used boundary conditions - inflow plane, outflow plane, symmetry plane

and solid wall. The solid wall boundary condition is formulated either by the standard wall-function approach or by a low-Reynolds-number procedure, depending

on the turbulence models used. The program also allows the user to introduce

boundary conditions other than these four types. Care has been taken to make the

module user-friendly. This is achieved by providing a few parameters defining the

types and locations of boundaries, which allows the user to choose the boundary

conditions given in the module by simply specifying these parameters. The system

of algebraic difference equations is solved using the alternating direction TDMA of

Thomas.

2.3 A Turbulence Module for the NPARC 2D Code

NPARC is a compressible flow code. Its 2D version is currently restricted to

plane or axisymmetric flows without swirling. The turbulence module (NP2DTS1) 4 written for NPARC 2D has much in common with the above FAST2DTS-1.

Only the following points are worthy of mention: (a) the cell-corner arrangement is

used, i.e., the flow variables are stored at grid nodes rather than at the centers of

control volumes; (b) contrary to NPARC and most other compressible codes, the

non-delta form of equations is used which leads to simple linearization and is more

effective to ensure the positiveness of the turbulent kinetic energy and its dissipation

rate than the delta form; and (c) four turbulence models have been built into the

module: Baldwin-Lomax, Chien, Shih-Lumley and CMOTT realizable models, the

first being the algebraic eddy-viscosity model and the last three the low Reynolds

number k-c models. The NPARC 2D code coupled with the module NP2DTS-1 have

been applied 5 to a series of flows including the flow over a fiat plate, in an ejector

nozzle, in a transonic diffuser, and a boat-tail nozzle flow. Both NP2DTS-1 and

FAST2DTS-1 were released at the 1994 Industry-Wide Workshop on Computational

Turbulence Modeling, held at the Ohio Aerospace Institute.
Turbulence Modules for CFD Codes 15

2.4 A Turbulence Module for the VSTAGE Code

NASA LeRC currently has a very extensive research program on turbomachinery

flow physics that is both experimental and computational. World class experimental facilities exist to provide a multitude of experimental data for both component

design needs and CFD comparisons. CFD turbomachinery tools are also well developed after many years of testing and validation. The VSTAGE code, developed

at Dr. Adamczyk's group of the Lewis Research Academy, is one such tool. At

present it is a proven CFD code which is being used by major U.S. engine companies in their design analysis. In a recent blind testing for a rotor compressor

(designated Rotor 37) organized by the ASME/International Gas Turbine Conference held at the Hague, 1994, the VSTAGE prediction turned out to be among the

best of predictions submitted by a total of 12 groups. However, the best predictions

for Rotor 37 were still less than satisfactory for the purpose of design, and this was

mainly attributed to turbulence modeling. The VSTAGE is equipped only with

the commonly used Baldwin-Lomax mixing length model. The limitations of such

model are quite apparent from the predictions for Rotor 37 case. It is questionable

whether improvements in predicting transonic rotor flows can be made using this

model. Thus a turbulence module is developed for VSTAGE. The reasons of developing this module are twofold: 1) to provide VSTAGE with more general and

advanced turbulence models beyond the algebraic level; and 2) using VSTAGE as

a common numerical platform to assess the performance of different existing and

newly developed turbulence models for turbomachinery flows. The module has the

same numerical framework as that of the NPARC module bat is written exclusively

for VSTAGE. As a result, it has the following unique features: (a) the cylindrical

coordinate system (x, r, 0) is used as the absolute (fixed) reference frame, and correspondingly the cylindrical coordinate components of the flow velocity are used in

the curvilinear transformation of equations; (b) the rotational effect is taken into account by transforming the absolute reference frame to the relative (rotating) frame;

(c) instead of the commonly used metrics such as J, _x, _o, ..., the module uses volumes and areas of control volumes as the major geometric quantities, which is in

line with the definition of variables in VSTAGE; and (d) the data transfer between

VSATGE and the module is via the Fortran common blocks. Three high Reynolds

number k-e models have been implemented, one standard and two CMOTT realizable models. The CMOTT models differ from each other in the dissipation

rate equation, one using the standard equation and the other the newly developed

equation 6. The wall function approaches axe used to formulate the solid wall boundaxy conditions. The code validation is being under way, and the preliminary results

for Rotor 37 test case look promising.

2.5 Calculation of Wall-Bounded Complex Flows

This study 7 concentrates on complex turbulent shear flows which are of great

interest in propulsion systems. These flows are backward-facing step flows, confined coflowing jets, confined swirling coaxial jets, U-duct flows and diffuser flows.

Most of these flows have complex structures. For example, the confined coflowing
16 J. Zhu

•jet combines several types of flow structures, such as the shear layer, jet, recirculation, separation and reattachment. Accurate prediction of these flows is of great

importance for engine design in all its key elements.

The turbulence model used in this study is a recently developed realizable Reynolds

stress algebraic equation model s,9 which is fundamentally different from the traditional algebraic Reynolds stress models. The present model is developed using the

invariance theory in continuum mechanics. This theory leads to a general constitutive relation for the Reynolds stress tensor in terms of the mean deformation rate

tensor and the turbulent velocity and length scales characterized by the turbulent

kinetic energy and its dissipation rate. Realizability is imposed directly on the

constitutive relation for the Reynolds stresses to determine the coefficients in the

relation. As a result, a realizable explicit expression for the Reynolds stresses is

obtained for general three-dimensional turbulent flows. Some model constants were

fine-tuned against a backward-facing step flow and then tested in other flows.

The calculations were performed using the FAST2D code. Grid independent

and low numerical diffusion solutions were obtained by using differencing schemes

of second-order accuracy on sufficiently fine grids. The standard wall function

approach was used for wall boundary conditions. Calculations using the standard

k-e (SKE) model were also carried out for the purpose of comparison.

Diffuser Flows. Two conical diffuser flows were calculated, one with a 8 ° total

angle (Trupp et al.'s case) and the other 10 ° (Fraser'case). In both cases, the flows

undergo strong adverse pressure gradients but remain attached. Although the flow

configuration looks simple, it is not easy to calculate this type of flow accurately,

especially for the boundary layer quantities. Fig.1 shows the variation of calculated

and measured wall friction coefficient Cf with the axial distance x/Ro (Ro is the

inlet duct radius). It is seen that the result of the present model is in good agreement

with the experimental data, while the SKE model overpredicts Cf along almost the

entire length of the diffuser. The calculated and measured displacement thickness

5* are compared in Fig.2. The comparison shows that the SKE model gives a good

prediction in the upstream region, but deviates significantly from the experiment

downstream; the present model prediction is good in the whole region. Fig.3 shows

the comparison of calculated and measured shape factor H. This is the case in

which the worst agreement with the measurement has been found for both models.

Nevertheless, the present model still performs considerably better than does the

SKE model.

U-Duct Flow. This case is the experiment of Monson et al. (1990) conducted

in a 180 ° planar turnaround duct. It features flow with large streamline curvature.

Calculations were compared to the experiment taken at a flow Reynolds number

of 10 6. Fig.4 shows the streamlines computed with the present model. A small

separation region is found at the bend exit. However, the SKE model did not

predict the flow separation. Fig.5 shows the comparison of calculated and measured

Cf along the inner wall. The bend is located between 21.7<s/H_<24.8. Both models

are seen to behave in the same manner and produce large discrepancies in the bend

region. The reason for this may partially due to the use of the wall function which
Turbulence Modules for CFD Codes 17

does not respond to the severe pressure gradient.

Backward-Facing Step Flows. Two backward facing step flows, measured by

Driver and SeegmiUer (1985) and Kim et al. (1978), were calculated. The former

(DS case) has a smaller and the latter (KKJ case) a larger step expansion. The computed and measured reattachment points are compared in Table 1. The calculated

reattachment point from the present model agrees well with the experiments. Fig.6

shows the comparison of the computed and the measured static pressure coefficient

Cp along the bottom wall. The SKE model is seen to predict a premature pressure

rise, which is consistent with its underprediction of the reattachment length, while

the present model captures the pressure rise quite well. Fig.7 shows the comparisons of predicted and measured turbulent stresses uu, vv and _ at the location

x=2 which is in the recirculation region. In the KKJ-case, no reliable experimental

data exist for the turbulent stresses due to the unsteadiness of the flow. However,

the experimental data of the DS-case is considered more reliable because of the

smaller unsteadiness of the flow. As compared with the results of the SKE model

in Fig.7, it is seen that the anisotropic terms in the present model increase u-'-_and

decrease V_, leading to significant improvements in both _-_ and _-_ except in the

near-wall region. On the other hand, the anisotropic terms have little impact on

u--_. The improvement obtained by the present model for _-_ is mainly due to the

reduction in C_, by strain rate.

Table 1. Comparison of the reattachment points

Case measurement SKE PRESENT

DS 6.26 4.99 5.82

KKJ 7+ 0.5 6.35 7.35

Confined Jets. The general features of confined jets are sketched in Fig.8. At

the entrance, two uniform flows, a jet of larger velocity and an ambient stream of

smaller velocity, are discharged into a cylindrical duct of diameter Do. The inlet flow

conditions can be characterized by the Craya-Curtet number Ct. The experiment

shows that recirculation occurs when Ct <0.96. For a given geometry, recirculation as well as adverse pressure gradients can be intensified by reducing the value

of Ct at the entrance. The separation and reattachment points of the predicted

recirculation bubbles are compared with the experimental data in Fig.9. The experiment indicated that as Ct decreases, the separation point moves upstream while

the reattachment point remains practically unchanged. The present model captures

this feature well and predicts both the separation and reattachment points much

better than does the SKE model. The variation of the pressure coefficient Cp along

the duct wall is shown in Fig.10. The pressure distribution is governed by the jet

entrainment as well as the contraction and expansion of the flow caused by the

recirculation bubble. The decrease in the ambient velocity induced by the entrainment gives rise to an adverse pressure gradient, while the contraction of streamlines

produces the opposite effect. These two mechanisms interact more intensely with
18 d. Zhu

each other as C, decreases and cause the pressure to vary little in the region upstream of the center of the recirculation bubble. However, in the downstream part

of the recirculation bubble, the deceleration of the flow sets up an adverse pressure

gradient, the slope of which becomes steeper as (7, decreases. Therefore, the ability

to capture the location of the recirculation center will have a direct impact on the

prediction of the pressure. Regarding the comparison between predictions and experiments, it is seen that although both models predict practically the same total

pressure rises which are in excellent agreement with the measurements, the present

model captures the steep pressure gradients better than does the SKE model for all

of the C, values.

Confined Swirling Coaxial Jets. This is the case experimentally studied

by Roback and Johnson (1983). Fig.ll shows the general features of the flow.

At the inlet, an inner jet and an annular jet are ejected into an enlarged duct.

Besides an annular recirculation bubble due to sudden expansion of the duct, a

centerline recirculation bubble is created by flow swirling. Fig.12 compares the

calculation of the centerline velocity with the experiment. The negative velocity

indicates the central recirculation. It is seen that both models predict the strength

of central recirculation and the front stagnation point quite well, but the present

model predicts the rear stagnation point much better than does the SKE model.

Fig.13 shows the comparison of calculated and measured mean velocity profiles at

x=5.1cm. Both models give reasonably good profiles which are within experimental

scatter, except for the peak values of the axial and radial velocities. Both models

have been found to give nearly the same results in the downstream region, which

can also be seen from Fig.12.

These comparisons show that the present realizable Reynolds stress algebraic

equation model significantly improves the predictive capability of k-e equation based

models, expecially for flows involving massive separations or strong shear layers. In

these situations, the standard eddy viscosity model overpredicts the eddy viscosity

and, hence, fails to accurately predict wall shear stress, separation, recirculation,

etc. We find that the success of the present model in modeling the above mentioned complex flows is largely due to its effective eddy viscosity formulation which

accounts for the effect of mean shear rates. According to the present model, the

effective eddy viscosity will be significantly reduced by the mean strain rate and

maintained at a correct level to mimic the complex flow structures.

2.6 A New k-e Eddy Viscosity Model

The model s has a new dissipation rate equation and a new realizable eddy viscosity formulation. The former is based on the dynamic equation for fluctuating

vorticity and the latter is derived from the realizability analysis. The model contains

the effect of mean rotation on turbulent stresses, and its dissipation rate equation

has an always positive production term, which is expected to enhance the numerical

stability in turbulent flow calculations, especially when using more advanced closure

schemes such as second-order closures.

Comparisons with experimental data show that the model performs better than

the standard k-e model for almost all the flows tested s. The performance of the
Turbulence Modules/or CFD Codes 19

model for complicated recirculating flows is demonstrated through calculations for

two backward-facing step flows, one (DS-case) with smaller and the other (KKJcase) with larger step height, both of which have been extensively used in the

literature to benchmark calculations of separated flows. Sufficiently fine grids, with

201x 109 points in the DS-case and 199x91 points in the KKJ-case, were used to

establish numerical credibility of the solutions. The computational dom_n had a

length of 50 step heights, one fifth of which was placed upstream of the step. The experimental data were used to specify the inflow conditions, the fully-developed flow

conditions were imposed at the outflow boundary, and the standard wall function

approach was used to bridge the viscous sublayer near the wall. The comparison with the standard k-e model and the experiments is shown in Table 2 for the

reattachment length.

Table 2. Comparison of the reattachment point locations

Case measurement standard model present model

DS 6.26 4.99 6.02

KKJ 74- 0.5 6.35 7.5

The comparison of the size of the separation bubble, the skin friction pressure

coefficients shows that the overall performance of the present model is better than

that of the standard model.

2.7 A New Reynolds Stress Algebraic Equation Model

A new Reynolds stress algebraic equation model has been developed using a

truncated constitutive relation 1°. This model differs from the previous one s mainly

in: (a) it has a simpler quadratic form; (b) it shows the proper lack of a rotation

effect on the isotropic turbulence, satisfying the rapid distortion theory; and (c) it

is fully realizable, i.e., it ensures both the positivity of the normal Reynolds stresses

and the Schwarz' inequality between turbulent velocity correlations. The model

were first calibrated using well-studied basic flows such as homogenous shear flow

and the surface flow in the inertial sublayer and then applied to complex flows

including the separated flow over a backward-facing step and the flow in a confined

jet.

Basic flows. Two basic cases were considered; they are Tavoularis and Corrsin's

(1981) homogeneous shear flow and the direct numerical simulation of channel flow

(Kim, 1990). The model gives b12 = -0.156, bll = -b22 = 0.123 for the homogeneous shear flow at U1,2k/e = 6.08 and gives b12 = -0.122, bll = -b22 = 0.14

for the channel flow in the inertial sublayer at U1,2k/e = 3.3. These results show

that the present model gives reasonable anisotropy of Reynolds stresses for both

the homogeneous shear flow and the boundary layer flow compared to the standard

k-e eddy viscosity model which gives bll = b22 = 0 for both the flows and gives

b12 - -0.273 for the homogeneous shear flow and b12 = -0.149 for the boundary

layer flow. Detailed comparisons with the experimental and DNS data are shown in

Table 3 for the homogeneous shear flow and in Table 4 and Fig.14 for the channel
20 J. Zhu

flOW.

Table 3. Anisotropy in the homogeneous shear flow

experiment

b12 -0.142

bll 0.202

b22 -0.145

standard

-0.273

0.

0.

present

-0.156

0.123

-0.123

Table 4. Anisotropy in the channel flow

512

bll

b22

DNS data standard present

-0.145 -0.149 -0.122

0.175 0. 0.14

-0.145 0. -0.14

Backward-facing step flows. Figs.15(a) and (b) show the friction coefficient

Cf at the bottom wall calculated with the SKE model and the present model; also

included in fig.15(a) are the experimental data for the DS-case, but no such data

are available for the KKJ-case. It can be seen that the grid refinement does produce

some differences for the results of the present model, more noticeable in the KKJcase, and this is also the case for the SKE results. This indicates that the solutions

obtained on the coarse grids are not sufficiently close to the grid-independent stage.

Recently, Thangam and Hur (1991) have conducted a highly-resolved calculation

for the KKJ-case. They have found that quadrupling a 166×73 grid leads to only

a minimal improvement. Therefore, the present results on the fine grids can be

considered as grid-independent. For the DS-case, the fine grid computations with

the SKE model and present model required 703 and 691 iterations, and took approximately 7.1 and 9 minutes of CPU time on the Cray YMP computer. In the

following, only the fine grid results are presented.

The wall friction coefficient Cf is a parameter that is very sensitive to the nearwall turbulence modeling. It is Cf that the various low Reynolds number k-e models

tested predict much worse than those using wall functions. However, the influence

of the near-wall turbulence modeling is mainly restricted to the near-wall regions.

It is seen from fig.15(a) that both the SKE model and the present model largely

underpredict the negative peak of Cf, pointing to limited accuracy of the wall

function approach in the recirculation region.

The computed and measured reattachment points are compared in Table 5. They

are determined in the calculation from the point where Cf goes to zero. The

reattachment point is a critical parameter which has often been used to assess the

overall performance of turbulence models as well as numerical procedures. Table 5

clearly demonstrates the significant improvement obtained with the present model.

It is important to mention that this improvement is mainly due to the behavior of

C_, in the present model, and that the anisotropic behavior of the turbulent stresses

only makes a marginal contribution to it.
Turbulence Modules .for CFD Codes 21

Table 5. Comparison of the reattachment points

Case

DS

KKJ

measurement SKE

6.26 4.99 5.80

7::t= 0.5 6.35 7.27

PRESENT

Figs.16(a) and (b) show the comparison of computed and measured static pressure

coefficients Cp along the bottom wall. In both cases, the SKE model is seen to

predict premature pressure rises which is consistent with its underprediction of the

reattachment lengths.

Finally, the comparisons of predicted and measured turbulent stresses itu, vv

and _-_ are shown in figs.17 and 18 at various x-locations. In the KKJ-case, no

experimental data for the turbulent stresses are available in the recirculation region,

and the reattachment point was found in the experiment to move forward and

backward continuously around seven step heights downstream of the step, leaving

an uncertainty of ±0.5 step height for the reattachment length. This also points to

some uncertainty in the measured turbulent quantities in the recovery region. On

the other hand, the experimental data in the DS-case should be considered more

reliable because of the smaller uncertainty of the reattachment location, indicating

a smaller unsteadiness of the flow. The SKE model gives unrealistic results about

normal Reynolds stresses: v-_ > _ at all the locations. In contrast, the present

model gives at least qualitatively correct results due to the non-linear terms which

increase u--_ while decreasing v-_, leading to an overall improvement in both _-_ and

v--g results.

Confined Jets. The predicted axial mean velocity profiles at two Ct numbers

are shown and compared with the experimental data in fig.19, where R and Um

are the radius of the cylinder and the sectional mean velocity, respectively. Both

models are seen to predict very well the upstream evolution of the flow. As for the

downstream development, the results of the present model remain in good agreement

with experiments while the SKE model underpredicts the centerline velocity decay

at all Ct numbers. The variation of the pressure coefficient Cp along the duct wall

is shown in fig.20. The pressure distribution is governed by the jet entrainment

as well as the contraction and expansion of the flow caused by the recirculation

bubble. The decrease in the ambient velocity induced by the entrainment gives rise

to an adverse pressure gradient, while the contraction of streamlines produces the

opposite effect. These two mechanisms interact more intensely with each other as Ct

decreases and cause the pressure to vary little in the region upstream of the center

of the recirculation bubble. However, in the downstream part of the recirculation

bubble, the deceleration of the flow sets up an adverse pressure gradient, the slope

of which becomes steeper as C_ decreases. Therefore, the ability to capture the

location of the recirculation center will have a direct impact on the prediction of

the pressure. Regarding the comparison between predictions and experiments, it is

seen that although both models predict the same total pressure rises which are in

excellent agreement with the measurements, the present model captures the pressure

distribution much better than does the SKE model for all the Ct values.
22 J. Zhu

3. Future Plans

1) To develop a 3D version of the turbulence module for the NPARC code;

2) To develop and validate the turbulence module for the VSTAGE code;

3) To test turbulence models developed at the CMOTT and elsewhere.

4. References

1 Georgiadis, N.J., Chitsomboon, T. and Zhu, J., "Modification of the two-equation

turbulence model in NPARC to a Chien low Reynolds number k-e formulation,"

NASA TM 106710 (1994).

2 Zhu, J., "FAST-2D: A computer program for numerical simulation of two-dimensional

incompressible flows with complex boundaries," Report 690, Institute for Hydromechanics, University of Karlsruhe (1991).

3 Zhu, J., "On the higher-order bounded discretization schemes for finite-volume

computations of incompressible flows," Coraput. Meths. Appl. Mech. Eng. 98,

pp.345-360 (1992).

4 Zhu, J. and Shih, T.H., "A turbulence module for the NPARC code," to be

presented at the AIAA Joint Propulsion Conference, San Diego, July (1995).

Yang, Z., Zhu, J. and Shih, T.H., "Calculations of inlet/nozzle flows using a new

k-e model," to be presented at the AIAA Joint Propulsion Conference, San Diego,

July (1995).

6 Shih, T.-H., Liou, W.W., Shabbir, A., Yang, Z. and Zhu, J., "A new k-e eddy

viscosity model for high Reynolds number turbulent flows - Model development

and validation," NASA TM 106721 (1994), also to appear in Computers and

Fluids.

7 Shih, T.H., Zhu, J. and Lumley, J.L., "Modeling of wall-bounded complex flows

and free shear flows," NASA TM 106513, also presented at the ASME Fluids

Engineering Conference, Symposium on Boundary Layer and Free Shear Flows,

Nevada, June (1994).

s Shih, T.H., Zhu, J. and Lumley, J.L.,"A realizable Reynolds stress algebraic

equation model," NASA TM 105993 (1993).

9 Zhu, J. and Shih, T.H., "A numerical study of confined turbulent jets," J. Fluids

Eng. 116, pp.702-706 (1994).

10 Shih, T.-H., Zhu, J. and Lumley, J.L., "A new Reynolds stress algebraic equation

model," NASA TM 106644 (1994), also to appear in Computer Methods in Applied

Mechanics and Engineering.
Turbulence Modules for CFD Codes 23

5 | I

4- • EXP

. -- PRESENT

0 , , , , I , , , , T T ,

0 5 10 15

X/Ro

Fig.1 Wall friction coefficient (Case of Trupp et al.)

_3

1

1.0

0.8

0.6

0.4

0.2

0.0

0

• EXP

- SKE ,,_

- PRESENT _... -_ -

l I I I | I I i l | I

5 10 15

X/Ro

Fig.2 Displ_cmnen_ thickness (Ca_ of Trupp e_ al.)

"r'2

e' I _

• I"

•__" ., EXP

..... SKE

PRESENT

I • , | t i , 1 ] i a _ f I I _ i I I

2 4. 6 8

X/Ro

Fig.3 Shape factor (Case of Fraser)

10
24 J. Zhu

Fig.4 Streamlines

"g ......... i ......... . .......

0 L,g'r...'_. P.I ........ .I ._e-.'... I ....

0 1 2 0 I 2 -'0.5 0.0 0.5

tOOuu lOOw --lOOuY

Fig.7 Tuzbulent stress profiles

(legend as in Fig.5)

.o

o

o

o

o

(eL

o

5 .... I .... ! .... .! .... I ....

• EXP

/_ ..... SKE

I0 [e_ .... PRESENT

I , , , , l = i | • I , , • •

_5 , , , , | s a l l

15 20 25 30 35 40

S/H

Fi_.5 Friction coefficient along the inner wail

0.2

0.1

0.0

-0.1

......... i ......... I ......... t ........ '_

/'/f "./ " . . .

j DS-case

......... ! ......... I ......... I ........ ,

I0 20 30 40

Fig.6 Pzesstue coefficient along the bottom wall

(legend as in Fig.S)

r,. _otential core ambient reqion

d_0:'_.c_ mixing region Dc

Ul U1

Fig.8 Flow configuzstion amd notstions

0.8

0.6

0.4

0.2

0.0

.... I ' ' | " ' I ....

//°

//

J/

l I t , | , , , , | f t 1 ; J I t I I

0 1 2 2;

x/Do

Fig.9 Sepazation and re.attachment points

(legend as in Fig.5)

4
Turbulence Modules .for CFD Codes 25

.... i .... 'l -'" "'" i .... * .... D " " "

I0

Ct=0.976

5

.... i .... i .... i .... i .... i • • 1=-

Ct=0.152 ,//re "

/

5 i '/Z

-5 .... I .... i .... i .... i .... i ....

0 1 2 3 4 5 6

x/Do

Fig.lO Press_xe coemcienl along the duct wail

Oegend as ia Fig.5)

c-_

¢.) --5

0

0

c 10

1.0 ..:.,_ .... ,...

°i__ • I'tl°

e I

0 e°

_: 0.5 •

+If

-1

'" f "}"

0 1 -0.5 0.5 -0.5 0.5

U (m/s) V (m/s) W (m/s)

Fig.13 Mean velocity pzofiles at T=5.1cm

(legend as ia Fig.I2)

EOo'r _rr.)m

RC-GX_ _'TW_Ew JETS

m

zt-----

Fig.ll Flow co_figazafion

1.5

"E" 0.5

v

:S o.o

-0.5

• ' " ' I .... I " " " ' J .... I ....

• EXP

..... SKE

PRESENT

L_ II •

l , , , I i i i i { ' + i , | , , l i I • • • •

0 10 20 30 40 50

x (cm)

Fig.12 Centerline velocity
26 J. Zhu

3

.

+

>

0

i i i i |

• DNS

Standard model

Present model

I I I I I I ! f I _ I | ! I | I I I I

2

+

I

I I I I ! I I I I _ I I ! I I ! I I I

100 200 300

y+

Fig.14 Direct comparison with DNS data of channel flow
Turbulence Modules for CFD Codes 27

2

o

i,,iD++,,|l+lll,,+l_l,l,,,,l'|_'''lll''

__--4_ _ ---o

(a) DS-case

lfflli,,,li,,,il,|Itill.t,,,,l. Stllt,t.

0 10 20 30 40

1

!o

2L3 ......... | .... , .... | ......... I ' ' 1 ' ' _ ' '__._

' ' (b) KKJ-case

llllIfIflItlfllllll_tllPlllll1111tl+l|!

0 10 20 30 40

X

Fig.15 Friction coefficient along the bottom wall.

present model, fine grid; -- --: present model, coarse grid;

.... SKE, fine grid; ." experiment
28 J. Zhu

o.

(D

0.2

0.1

0.0

-0.1

0

, , , , , , , , , | , , , , l , , o • | , , , , , , , , , | , o , , , • , , ,

/./ . . .

j

/

(a) DS-case

. . . I i i i 2 i I i , i , i I ! I I I , ! , i , , , I ! | t | t ! i i | i .

I0 20 30 40

X

0.4

0.3

o.2

0.1

0.0

-0.1

0

,/_'-"--_ _ ----.

!

(b) KKJ-case

llllll,ll|l,+lllllllltlllllllllllll]lll

10 20 30 40

X

Fig.16 Pressure coefficient along the bottom wall.

--: present model;----: SKi_; e- experiment
Turbulence Modules for CFD Codes 29

3

2

>,

1

2

1

0

3

2

>,

1

0

.... i .... I'

_1,,,,11

.... I .... I '

lOOuu

' '' ' i ' ' ' ' I '

.... I .... l '

)

• _

i i i I | i

I

_/c_,.... ,

0 1 2

100w

........

lll. IJ,,ll,_

,,o ,,,,111,,

f ,_io,,

i_l| llllilll,

X=2

-0.5 0.0 0.5 1.0

- 100uv

Fig.17 Turbulent stress profiles in DS-case.

-: present model;----: SKE; ,: experiment
30 J. Zhu

3

2

1

0

3

2

1

'l''''l''''l''''l

|"

_,'I''''I''''I''''I'

,I .... I,

0 I 2 3 4

''1 .... I''''I''''I

_''1''"1''''1''''1' 1

. ?

,,,I,,,,I,,,,I

0 1 2 3 4-

'll .... i .... i"' 1

X=7.660

b

-&..

X=2.667

-0.5 0.0 0.5 1.0

100uu I OOw -I OOuv

1.5

Fig.18 Turbulent stress profiles in KKJ-case.

----: present model;----: SKE; o: experiment
Turbulence Modules for CFD Codes 31

x/Do=O.375 1.125 1.875 2.625

1.0 _,,,i .... J .... i .........

0.5

0 5 10 15 0 3 6 9 0 3 6 0 3 6

1.0

0.0

' ' ' I ' ' i I '

Ct=O. 152

4O 8O -5 2O 45 -5

_' ' I ' ' I

U/Um

tt t _ I

w

25 -5 10 25

Fig.19 Axial mean velocity profiles.

: present model;----: SKE; e. experiment
32 J. Zhu

I0

5

0

-5

10

-5

10

_ 5

g

o o

-5

10

5

0

-5

10

5

0

-5

, ,,,i,,,, i, ,,, i,,,,i ,,, ,i ,_,,i,,,,i , , , ,

t=0.976

.... I,,,, ,,,, ,,,,I,,, , I,, .... I, ,, ,

r_- ' _" I I ' ' I

= .... I .... I .... I .... I ' ' ' ' I ' ' ' ' I .... I .... .

e Ct=0.714

',,,,I,,, , I,,,,I,, ,, I,,, ,I ,o,,I ,,,,I, i, ,

.... I .... I .... I .... I .... I .... I .... I ....

//__ .......

/ •

, _.____//j • Ct=0.506

,,,,I,,,, l*,,,l,, ,,I,,, ,I,, ,,l,,,ll , i, ,

"' ' ' ' I ' ' ' ' I/_" '......l__ ' ' ' I ' 'J• ' ' I ' '_ ' ' I ' ' i ' ' I ' ' ' ' I _ ' ' '

__.._._/* Ct=0.505

.... I .... I .... | .... I .... i .... I .... I''''.

/

/

,//"

Ct=0.152

,,,,I .... I .... I .... I .... I,_,,I, ,,,I, ,, ,

1 2 3 4 5 6 7 8

x/Do

Fig.20 Pressure coefficient along the duct wall.

_: present model;----: SKE; ." experiment
Center for Modeling of Turbulence and Transition

Research Briefs - 1995

33

Eddy Diffusivity Model for

Turbulent Heat Transfer

Aamir Shabbir

1. Motivation and Objectives

Commonly, the turbulent heat transfer is calculated by assuming a constant value

of the turbulent Prandtl number. A mixing length turbulence model or a two

equation k - e model is usually used in such calculations to calculate the turbulent

eddy viscosity. Obviously this approach has limitations in situations where the

turbulent Prandtl number is not constant. Realizing this, some recent studies have

proposed new eddy diffusivity models in which the transport equations for the scalar

variance and its dissipation rate are solved to calculate the thermal diffusivity.

These include, for instance the work of Nagano and Kim 1 and Youssef, Nagano

and Tagawa 2, Hattouri, Nagano, and Tagawa _. The last study presents an updated

and refined version of Nagano and Kim model. There are several other studies

reported in the literature where the near wall modeling issues are discussed. In

this report only the high Reynolds number version of the models is discussed, and

therefore, reader should consult other papers which deal with this issue (for instance

the review paper by So4). Recently Schwab and Lakshminarayana 5 also propsed

a two equation model for scalar field in which a scalar time scale equation is used

instead of the scalar dissipation rate. They sucessfuUy applied their model to a host

of homogeneous flows.

In last year's research briefs I descirbed my project aimed at developing a new

two equation eddy diffusivity model for calculating turbulent heat transfer. In this

model transport equations for the scalar variance and its dissipation rate are used to

calculate the turbulent eddy diffusivity. As was pointed out then, the present effort

differs from the other work in the following two respects. (1) In the above cited

works, the extension of the scalar dissipation rate equation is based upon the work

of Newman et al. 6 who developed its production/destruction mechanisms analogous

to those of the mechanical dissipation rate equation. The model equation proposed

in the present study is based on the exact transport equation for scalar dissipation

and, its production/destruction mechanisms differ from those proposed in the other

studies. (2) The model coefficient in the the scalar flux constitutive relation used

in the present study is not a constant but is a function of the local invariants.

2. Work Accomplished

2.1 A New eddy viscostiy k- e model for High Reynolds Number Flows

Some of the effort was directed in developing and assessing the performance of

a new two equation k-epsilon model in a joint effort with the other CMOTT researchers. The details of the model development and application can be found in
34 A. Shabbir

the paper by Shih et al.7 The model performs reasonably for a host of benchmark

flows which includes: rotating homogeneous shear flow; various fiat plate boundary

layers; mixing layers; and jets. We, in this section, only show the application of the

model to the rotating homogeneous shear flow, for which a large eddy simulation

was carried out by Bardina et al. s. We will compare the results from the present k-e

model as well as the standard k - e model with the LES data for four different cases

of _/S (which are _/S-O.O, _/S--0.50, 12/S=0.25, and _/8-0.50). The initial

conditions in all these cases correspond to isotropic turbulence and eo/Sko = 0.296.

Figure 1 (a) compares the evolution of turbulence kinetic energy, normalized by its

initial value k0, with the non-dimensional time St for the case of 12/S = 0.0. For

this case both the present and the standard k - e (denoted by ske hereafter) models

show the trends exhibited by LES, with the present model closer to the LES data.

Figure 1 (b) shows the comparisons for the case _/S = 0.25. The LES shows that

the growth rate of the turbulence kinetic energy is increased over the no rotation

rate case. The present model is able to pick up this trend while the ske model does

not. Figures 1 (c) and 1 (d) compare the evolution of turbulence kinetic energy for

two more cases of _2/S - 0.5 and l_/S = -0.5. For the first of these cases the LES

shows that the growth rate of the turbulence kinetic energy decreases over the no

rotation rate case. The present model is able to pick up this trend and although

the agreement between the present model and the LES is not as good as it is for the

other cases, it still is a lot better than the ske model. For the case of _/S = -0.5

the ske model does not show the effect of rotation on turbulence as it gives the same

growth rate of turbulence kinetic energy as it did for the no rotation case, a result

which is already known. On the other hand the present model is in reasonable

agreement with the LES data as it shows the decay of the turbulence kinetic energy

with time.

2.2 A New eddy diffusivity 02 -eo model for High Reynolds Number Flows

The details of the development of the scalar (or thermal) dissipation rate equation

and the scalar flux can be found in the last year's research briefs. Since then the

model has been extended to inhomogeneous flows by modeling the diffusion term

in the disspiation rate equation. A simple gradient diffusion type model is used for

the diffusion term (u--_-_),j. With this the scalar dissipation rate equation is given

by

+ v o,j =( -_ o,j +co oS + Co - r vqT_ ¢ - Ce3(k +eoe (i)

where the model constant a s is to be determined. To do so we consider the log-law

region of a wail, where the turbulent scalar (or heat) flux is constant. In this region

the advection term can be ignored and equation (1) reduces to

0 C_T0E0 E%/f_ E0E

Oy ( aS Oy ) + Col eoS+ Cos _ ¢ - Co3 (k + v_) =0 (2)
Eddy Diffusivity Model for Turbulent Heat Transfer 35

Using the log-law relations, and the fact that in the log-layer production equals

the dissipation rate, the above equation gives the following relation for the model

constat a¢

(Cx/C')k2 (3)

(7¢ =/_1/2/_

where the Von Karman constant K is 0.41 and Ko is 0.4829. We know that in

the log-layer C, = 0.09 and CA = 0.1. With these values we obtain a¢ = 1.8

from (3). The expression for the model constant CA used previoulsy was somewhat

cumbersome. We have now considerably simplified it and is given by equation (7).

The model for the scalar (or thermal) field can now be summarized by the following

equations.

o__ o,TN- _oe

u5 -( o, ,j),j- 2,,,o -

+ co, os+ v - T

ker2_l 2_. k a 2 1/2

ui-'---0= - C), --etr; - / "-'" +-_(r) (a_V_,_ + a3V_,i)e,j

(4)

(s)

(6)

CA -- (2 + 2r + 0.Sr 2)

26 + 3.2zj 2 + 2( 2

(7)

Co2 =0.63, Con =C2-1+r, at=l.O, a 4,=1.8

2.3 Application of Model

Since the expression for CA has been simplified than what was reported previously,

it is instructive to re-calculate the flows which were calculated previously. Figures

1-7 show the evolution of the temperature variance and its dissipation rate for the

seven different cases of Sirivat and Warhaft 1° experiment. We note that the present

model is in reasonable agreement with the experiment, and performs better than

the Hattouri, Nagano and Tagawa model (HNT). The results on the overall are

same as what was reported previously. Same is true for the homogeneous shear flow

with a constant temperature gradient, for which the results are shown in figure 8.

The model was also applied to the flat plate boundary layer whose surface is

heated to a constant temperature. The experiment of Gibson et al. 9 provides one

test case for such a flow. This fl0w constitutes an initial value problem in the streamwise direction and a boundary value in the cross-stream direction. It was calcualted

using a quasi-implicit finite difference scheme. Since only the high Reynolds number

form of the present model is developed so far, it is used here in conjunction with the
36 A. Shabbir

wall functions. For the k and e quantities the standard wall functions were used.

The wall functions used for the thermal field are summarized in Appendix B. (It

is found that with these wall functions the log-layer temperature profile agrees well

with the data but the Stanton number is underestimated by about 7%. Work is

underway to address this aspect of the wall functions.) The equations set consisted

of: mean momentum and mean energy equations; turbulence kinetic energy and its

dissipation equation; temperature variance and its dissipation equation.

The results for the mean temperature are shown in figure 9. We note that the

agreement with the experiment is reasonable. For comparisons purposes results

are also shown for the HNT model which agrees with the experiment very well.

Note that the HNT model was not used with the wall functions but was instead

integrated to the wall.

3. Future Plans

The present model will be extended for integration to the wall. This will necessitate the formulation of damping functions for the model. However, due to the

changing research environment at NASA, most of the next year will be spent on

application of two-equation models to the turbomachinery flows.

Acknowledgements

The help of Dr. Z. Yang in the numerical aspects of the boundary layer calculations is gratefully acknowledged. Helpful discussions with J. Schwab of NASA

Lewis are also acknowledged.

4. References

1 Nagano,Y. and C. Kim. "A two equation model for heat transport in wall turbulent shear flows", J. Heat Transfer, 110, 583-589, (1988).

2 Youssef, M.S., Y. Nagano and M. Tagawa. "A two-equation heat transfer model

for predicting turbulent thermal fields under arbitrary wall thermal conditions",

Int. J. Heat Mass Transfer, 35, 11, 3095-3104, (1992).

a Hattori, H., Ngano, Y., and Tagawa, M. "Analysis of Turbulent Heat Transfer

under Various Thermal Conditions with Two-Equation Models", Department of

Mech. Engg., Nagoya Institute of Technology, Gokiso-cho, Showa-ku, Nagoya

466, Japan.

4 So, R.M.C., S.P.Yuan, and T.P. Sommer. "A hierarchy of near wall closures for

turbulent heat transfer", Trends in Heat, Mass, and Momentum Transfer, Vol. 2,

pp. 203-221.

5 Schwab, J.R. and Lakshminarayana. "Dynamic and Thermal turbulent time scale

modeling for homogeneous shear flows", NASA TM 106635, 1994.

8 Newman, G.l_, B.E. Launder and J.L. Lumley. "Modeling the behaviour of

homogeneous scalar turbulence", J. Fluid Mech., 111, 217-232, (1981).

Shih, T.-H., W.W. Liou, A. Shabbir, Z. Yang and J. Zhu. "A new high Reynolds

number k - e model", Computers Fluids, 24, 3,227-238, 1995.
Eddy Diffusivity Model for Turbulen_ Hea_ Transfer 37

s Bardina, J., J.H. Ferziger and W.C. Reynolds, "Improved Turbulence Models

Based on Laxge-Eddy Simulation of Homogeneous, Incompressible Turbulent

Flows". Stanford Univ. Report TF-19, 1983.

9 Gibson, M.M., Verriopoulos, C.A. and Nagano, Y., "Measurements in the Heated

Turbulent Boundary layer on Mildly Curved Convex Surface", Turbulent Shear

Flows 3 (ed. Bradbury, L.J.S., et al.), Springer-Verlag, 80-89, 1982.

10 Sirivat, A., and Z. Warhaft. "The effect of a passive cross-stream temperature

gradient on the evolution of temperature variance and heat flux in grid turbulence", J. Fluid Mech., 128, 326-346, (1983).

11 Tavoularis, S., and S. Corrsin. "Experiments in nearly homogeneous turbulent

shear flow with a uniform mean temperature gradient", Part 1, J. Fluid Mech.,

104, 311-347, (1981).

12 Launder, B.E. "On the computation of convective heat transfer in complex turbulent flows", J. Heat Transfer, 110,1112-1128, (1988).

Appendix A

The k- e model used in the present study is that of Shih et al. 7 and is summarized

below.

U _ zjOk - e

OE 1]T e 2

+ cl - z

(A1)

(A2)

(A3)

1

vT =C_, k--2 e , C, - 4.0 + AsU*(k/e) (A4)

1 0U_ OUj. 1 (OU_ OU t

(A6)

W =SijSjkSk_/S _,

1

¢ = -_arc cos(v/-6W), As = x/-6Cos¢ (A7)

(AS)

Note that the 12ij is the mean rotation rate as viewed in the rotating reference

frame and wk is the angular velocity at which the reference frame is rotating. The

constants in the e equation are: C_ = 0.42 and C2 = 1.9.
38 A. Shabbir

Appendix B

This appendix gives the wall functions which were used to obtain the boundary

conditions for the thermal field (i.e the values of 0, _, and e0 at the first grid point

off the wall). Note that the standard wall functions were used for the velocity field

and, therefore, are not listed here.

For the mean temperature the standard log-law for a flat plate with constant

surface temperature was used as a boundary condition (see e.g. Launder 12) and is

given as

O+_ AO I

Or -Ko Iny+ ÷ Co (B1)

where AO is the difference between the wall temperature and the local temperature.

From the experiment of Gibson et al. 9 Ko = 0.482 and Co = 3.8. For the turbulent

part of the temperature field the local equilibrium assumption was invoked (i.e.

Po = co) to obtain the following relation for the boundary conditions on _.

m

02

=0/2/0 (82)

D

Once 02 is known the boundary condition on eo can be readily obtained by using

the definition of the time scale ratio.

£0 --r( £ -- _-_)0 2 (83)

Note that for the log-layer the time scale ratio r is 2.0; C t, = 0.09; and Cx = 0.1.
3.0

2.0

1.0

Eddy Diffusivity Model for Turbulent Heat Transfer

' ' ' I ' ' ' I ' ' ' I

dT/dy=10.3 deg C/m, U=6.3 rn/s ,

/

• Expt. Data /

/

Present Model ,'

/

HNT Model /"//_"

,g

_,-

39

00. ' ' ' ' ' ' ' J ' , , '

20.0 60.0 100.0 140.0

x/M

Figure 1. Evolution of the normalized temperature variance and thermal dissipation

rate in the experiment of Sirivat and Warhaft (1983). The temperature variance

is normalized as _ = 02/020 and the thermal dissipation is normalized as e_ =

"-___._L__. _ For this case 0"_o = 0.0128 °C2, Uo = O.145515m/s, and lo = 0.011937m.

a_ouo/to

3.0

2.0

1.0

' ' ' I ' ' ' I _" ' '

dT/dy=4.48 deg C]m, U=6.3 a_a/s

/

• Expt. Data /

- Present Model // /

HNT Model///

/// • •

00_-_ 0 • @-@-0---@-

O0. , , , I , , , I , i , I ..

20.0 60.0 100.0 140.0

xlM

Figure 2. Evolution of the normalized temperature variance and thermal dissipation

rate in the experiment of Sirivat and Warhaft (1983). The temperature variance

is normalized as _' = 02/02o and the thermal dissipation is normalized as e_ =

-- '" . For this case 82o = 0.002287 °C2, Uo = O.145515m/s, and lo = 0.011937m.

O%uo/to
4O

3.0

2.0

1.0

A. Shabbir

' ' _ I ' ' _ I ' f

dT/dy=3.68 degC/m, U=6.3m/s//

/

• Expt. Data /

/

Present Model ,,

HNT Model /// .

' I

_0 • • • ID_-.-IL

00. , t I I I , I I l i , I

20.0 60.0 100.0 140.0

x/M

Figure 3. Evolution of the normalized temperature variance and thermal dissipation

rate in the experiment of Sirivat and Warhaft (1983). The temperature variance

is normalized as 0T_ = 02/020 and the thermal dissipation is normalized as e_ =

_'f For this case 02--'o= 0.001705 °C2, uo = O.145515m/s, and lo = 0.011937m.

e%uo/to

3.0

2.0

1.0

' ' ' I ' ' ' I I ' ' '

dT/dy=8.1 deg C/m, UT--t3.4 m/s

• Expt. Data // •

- Present Model

' HNT Mod_l • •

I i • O/

j/ • O_

i 10

• • -0 -O. - • -O- - - lID00. I i i I i I l I , I , l

20.0 60.0 100.0 140.0

x/M

Figure 4. Evolution of the normalized temperature variance and thermal dissipation

rate in the experiment of Sirivat and Warhaft (1983). The temperature variance

is normalized as _ = 02/020 and the thermal dissipation is normalized as e_ =

"" . For this case 0--_o= 0.009059 °C2, Uo = 0.074701m/s, and lo = 0.011071m.

O%uo/to
2.0

1.0

Eddy Diffusivity Model for Turbulent Heat Transfer

' ' ' I ' ' ' I ' ' ' I

dT/dy=2.24 deg C/m, U=3.4 m/s

• Expt. Data

Present Model

HNT Model •

41

00. , , , I , , I I , , , I

20.0 60.0 100.0 140.0

x/M

Figure 5. Evolution of the normalized temperature variance and thermal dissipation

rate in the experiment of Sirivat and Warhaft (1983). The temperature variance

is normalized as _-_ = 02/020 and the thermal dissipation is normalized as e_ =

"_ . For this case 0-_o = 0.000924 °C'2, Uo = O.074701m/s, and Io = 0.01107Ira.

02ouo/lo

3.0

2.0

1.0

' ' ' t ' ' ' I ' '

dT/dy= 1.81 deg C/m, U-_.4 m/s

/

• Expt. Data /

/

Present Model

/

HNT Mod_l • •

/// •

.

* 1

0.0 , , , I , , , I , o , I

20.0 60.0 100.0 140.0

x/M

Figure 6. Evolution of the normalized temperature variance and thermal dissipation rate in the experiment of Sirivat and Warhaft (1983). The temperature

variance is normalized as _' = 02/02o and the thermal dissipation is normalized

as e_ - -- "_ . For this case /72o = 0.0004471 °C2 uo = 0.07470Ira/s, and

0201LO//O

lo = 0.011071m.
42

3.0

2.0

1.0

A. Shabbir

' ' ' I ' ' ' I ' " ' '

dT/dy=l.78 deg C/m, U=3.4 m[s

/

• Expt. Data //

Present Model /"

HNT Model / /

/// _ -

°_"_-0---- $ - • - at - i- - -a0.0 , , , t , , , I , , , I

20.0 60.0 100.0 140.0

x/M

Figure 7. Evolution of the normalized temperature variance and thermal dissipation rate in the experiment of Sirivat and Warhaft (1983). The temperature

variance is normalized as if2' = 0_/02o and the thermal dissipation is normalized

For this case 0Zo = 0.0004955 °C_, uo = O.074701m/s, and

O%uollo "

l0 = 0.011071m.

3.0

2.5

2.0.

¢q

1.5

1.0

' I ' I ' I ' I ' I '

• Expt.

Present Model

.... HNT Model

t

/

/

/

/

/

/t'/_

05. , I , I , I , t , I ,

6 7 8 9 10 11 12

x/H

Figure 8a. Evolution of the normalized temperature variance in the experiment

of Tavoularis and Corrsin (1981). The temperature variance is normalized with

O-_,._.f= 0.01 °C2.
¢'4

_9

O

G)

tO

0.25

0.15

Eddy Diffusivity Model for Turbulent Heat Transfer

' I ' I ' I ' I '/ I '

• Expt.

Present Model

HNT Model

/

/

/

/

/

/

/

/

/

43

0.05 _ I , I , I , I , I ,

6 7 8 9 10 11 12

x/H

Figure 8b. Evolution of the thermal dissipation rate in the experiment of Tavoularis

and Corrsin (1981).

T+

3O

2O

10

0

f ' I ' I ' I ' I '

1: 0 Experiment

HNT model (no wall func)

--Present Model with wall fu one

I I v I , I a I n

0 200 400 600 800 1000

+

Y

Figure 9. Mean temperature profile in wall units for the constant temperature flat

plate boundary layer for the experiment of Gibson et al. 9.
Page intentionally left blank 
Center for Modeling of Turbulence and Transition

Research Brie]s - 1995

45

Computations of Compressible

Turbulent Benchmark Flows

William W. Liou

1. Motivation and Objective

The development of advanced propulsion systems for high speed aerospace vehicles will require accurate computational models of turbulence that can be used in

the CFD calculation of individual flow components. Before a new model is used in

the calculation of complex flows in propulsion systems, the model has to be assessed

in simple benchmark flows that not only contain the essential physical mechanisms

at work in engine component flows, but are also well-documented. An immediate

objective of the research activity described here is to identify these compressible

benchmark flows and performed a preliminary assessment of models developed locally at CMOTT. The long term objective is to develop second-order closure models

for compressible flows.

2. Work Accomplished

2.1 Model Validation in Compressible Flows

Six different types of two-dimensional flow fields typically encountered in neropropulsion systems were selected. They are: (1) compressible free shear layers; (2)

compressible boundary layers; (3) transonic bump flows; (4) supersonic ramp flows;

(5) oblique shock-wave/turbulent boundary-layer interactions; (6) can combustor.

Sketches of the flows are shown in Fig. 1.

(1) The compressible free shear layer is an integral element in engine component

flows, such as combustor and nozzle flows. In many instances, the turbulence mixing

in the shear layers is a major factor in the operating efficiency of an engine component. Experiments seem to support a state of self-preserving in a fully developed

shear flow, which is describable by local characteristic scales. It is also observed

in experiments that the spreading rate of a compressible turbulent shear layer is

far smaller than an incompressible shear layer with the same velocity and density

ratios. The reduced turbulent mixing influences significantly the combustion efficiency and the noise emitted from high speed jet flows. Therefore, it is important

that turbulence models can predict correctly the compressible free shear flow; in

particular, its spreading rate.

(2) To evaluate model performance in wall bounded compressible flows, models

will be assessed in fully developed compressible turbulent boundary layers. It is generally accepted that for free stream Mach numbers less than about 5, the turbulent

structure in compressible boundary layers is nearly the same as in the corresponding

incompressible (constant density) flows. This is Morkovin's hypothesis. Therefore,

with properly scaling, models that gives good results in incompressible flow may

PRECEDING PP, GE BLA_'tK NOT Fit._,_E[_

tkA_L _ .... z::'L"J:'-'"'_'--Z-_ aL,_'"4'
46 W. W. Liou

also perform well in compressible boundary layer with free steam Mach numbers

less than 5.

The above two cases provide tests for models in a "clean" environment. The

following four cases involve flows with either shock-waves or chemical reactions.

(3) Case 3 resembles the transonic flow over an airfoil, which affects the performance of engine components such as compressors and turbines. The strong

viscous/inviscid interactions in the flow, through a change in the effective bodyshape, pose a critical test for turbulence closure models. Two flows were chosen.

The first one corresponds to the experimental set-up of a flow studied by Delery

et. al 1. The boundary layer developed on the bump mounted on a wind-tunnel

wall remains attached throughout the interaction. The flow speed becomes supersonic atop the bump. The wind-tunnel is choked. This flow was also chosen in a

recent EUROVAL _ effort to evaluate the response of a turbulent boundary layer

developing on an airfoil to shock-waves. There, the case was designated as ONERA BUMP A. The second configuration corresponds to an experiment by Bachalo

and Johnson 3, where a toroidal bump is mounted on a cylinder. The axisymmetric

arrangement eliminates the three-dimensional effects due to the boundary layer on

the side walls of the wind tunnel. The boundary layer separates downstream of the

shock wave. The axisymmetric configuration has made this flow particularly suitable for the study of a flow with separation resulting from strong shock/turbulent

boundary-layer interactions.

(4) Cases 4 and 5 involve interactions between shock-wave and boundary layer

in supersonic flows. The compression ramp flow is a classical problem that encompasses many complex phenomena and is of great practical importance in the

design of engine components, such as inlets and turbomachinery. The four ramp

flows selected were investigated in a series of work by Settles and co-workers 4. The

free stream Mach numbers are less than 3 and, depending upon the ramp angle,

flow separation may occurs at the corner of the ramp. The flows were proposed for

model testing in shock wave/turbulent boundary-layer interactions in the 1980-1981

AFSOR-Stanford Conference on complex turbulent flows.

(5) Case 5 includes another fundamentally important flow that is critical to the

performance of inlets and turbomachinery: the impingement of an oblique incident

shock-wave on a wall 5. A reflected shock is formed to turn the downstream flow

parallel to the wall. Again, if the shocks were sufficiently strong, the turbulent

boundary layer may separate. The challenge for a turbulence model is to predict

the incipient separation of the boundary layer and the subsequent separated flow.

(6) The flow in a engine combustor typically involves highly turbulent reacting

regions. Therefore, a robust turbulence model is just as important as a quality

chemistry model in the design of a high efficiency and low pollutant emission combustor. An axisymmetric can combustor configuration 6 tested in UC Irvine was

examined. The experimental set-up includes an air-swirler and gaseous fuel injection of propane. The fuel and air are not premixed before entering the combustor.

This configuration is representative of a gas turbine engine combustor.

These test cases are well-documented and should provide good data base for
Compressible Turbulent Flow Computations 47

objective evaluations of turbulence models. For instance, some experiments have

reported results obtained using different measurement techniques, i.e. hot wire vs.

LDV. In fact, some flows were also selected by previous large-scale model evaluation

efforts.

It is more efficiently to calculate the free shear layer and the flat plate boundary

layer with a boundary-layer equation solver. When it is necessary to use a NavierStokes solver, the model solution of a flat plate boundary layer is compared with the

solution obtained from the boundary-layer equation solver to validate the model implementation in the Navier-Stokes solver. The COMTUR code 7 is currently used in

the solution of compressible Favre--average Navier-Stokes equations. The FAST2D

code s is used for solving the variable density form of the Navier-Stokes equations.

The chemical reactions were modeled by using a CMOTT scalar probability density

function(PDF) code, LPDF2D.

In the following, a few sample results of this model assessment effort are shown.

More details can be found in NASA TM reports that are in preparation.

Fig. 2 shows the skin friction and the wail pressure distributions for a supersonic

flow over a 24 ramp 4. The incoming free stream Mach number is 2.8. and the unit

Reynolds number is 6.3×107/m. No compressibility corrections were used. The

results show that the standard k - e model and the multiple-scale model 9 give

reasonable predictions for the flow recovery and the flow separation. RNG model l°

and the S&Z 11 model predict a larger separation region than the measurement.

Note that only the linear part of the constitutive relation of the S&Z model was

used in this calculation.

Fig. 3 shows the variations of the skin friction coefficient and the wall pressure for

an oblique shock/turbulent boundary-layer interaction. The incoming free stream

Mach number is 2.89 and the unit Reynolds number is 5.73 × 107/m. The impinging

shock is strong enough to cause the boundary layer to separate. KE1 model 12"13,

which is a modified low Reynolds number model of Shih and Lumley 14, performs

better than the Chien's model _5. The high Reynolds number KE2 _s model and the

standard k - e model gives reasonable predictions for both quantities.

The predicted and measured wail pressure distributions for the ONERA BUMP

A are shown in Fig. 4. KE216 model give a slightly better predication downstream

of the shock wave than the other models. All the model predictions asymptote to a

value higher than the measurement in the fully developed region. This may be due

to three-dimensional effects, such as the boundary layer on the wind-tunnel walls.

The computed and measured wall pressure distribution for the Bachalo and Johnson transonic bump is given in Fig. 5. The incoming free stream is subsonic. A

supersonic pocket is formed atop the bump as the ilow accelerates through the

bump. The combined effect of the shock-wave and the bump geometry results in

flow separation near x/c = 0.7. KE2 model 16 shows the best overall agreement

with the measurement than the three other models tested. KE1 model 12'13, a low

Reynolds number model, also gives better results than the Chien's model.

In Fig. 6, the temperature profiles at four measurement stations in the axial

direction for the UC Irvine axisymmetric can combustor are shown. Compared with
48 W. W. Liou

previous published results 17 obtained by using the KIVA code, the pdf chemistry

model give a better agreement with measurement. The standard k - e model was

used in both calculations.

2.2 Low-Reynolds Number Multiple-Scale Model

Damping functions were introduced into the high-Reynolds number multiplescale model developed earlier 9 so that the model transport equations can be integrated down to the wall. The Kolmogorov behavior of near-wall turbulence proposed by Shih and Lumley 14 were applied to determine the boundary conditions

for the turbulent quantities at the wall. The details of the model can be found in

a NASA TM in preparation. Fig. 7 shows the skin friction distribution for a flat

plate turbulent boundary layer. The free stream Mach number is 2.87. The current

model agrees better with the Van Driest II formula than the Launder and Sharma

model.

3. Future Plans

(1) Examine the effects of turbulence models on chemically reacting flow calculations. Select a combustor of a different set-up, if necessary.

(2) Broaden the validation cases to include realistic geometries found in gas turbine engine. This will provide means of realistically projecting model performance

in engine component flows.

(3) Continue the development of compressible turbulence models.

4. References

1 j. Delery and J. Reisz, "Analyse experimentale d'une interaction choc-couche

limite turbulente a Mach 1.3 (decollment naissant)" ONERA Rapport Technique

No. 42/7078 AY 014, Chatillon, December 1980.

2 EUROVAL-A European Initiative on Validation o/CFD Codes. Edt. W. Haase,

F. Bradsma, E. Elsholz, M. Leschziner and D. Schwanborn, 1992.

3 W. D. Bachalo and D. A. Johnson, "Transonic, turbulent boundary-layer separation generated on an axisymmetric flow model," AIAA J. 24, 437-443 (1986).

4 G. S. Settles, I. E. Vas, and S. M. Bogdonoff, "Details of a shock-separated turbulent boundary layer at a compression corner," AIAA Journal, Vol.12, pp.1709-

1715, 1976.

5 D. C. Reda and J. D. Murphy, "Sidewall boundary-layer influence on shock

wave/turbulent boundary-layer interactions," AIAA Journal, Vol.ll, pp.1367-

1368, 1973.

6 R. Charles, Detailed Data Set: Velocity (U, W, U_,_, W_,_) and Temperature

(T) Measurement in the Axisymmetric Can Combustor (ASCC) for a Parametric Variation in Inlet Conditions, Technical Report UCI-ARTR-87-6, Univ. of

California, Irvine, 1988.
Compressible Turbulent Flow Computations 49

7 p. G. Huang and T. J. Coakley, "An implicit Navier-Stokes code for turbulent

flow modeling," AIAA paper 92-0547.

s j. Zhu, FAST-2D: A Computer Program for Numerical Simulation of Two-Dimensional

Incompressible Flows with Complex Boundaries, Institut Fur Hydromechanik,

Report No. 690, Universitaat Karlsruhe.

9 W. W. Liou, T.-H. Shih, and B.S. Duncan, "A Multiple-Scale Model for Compressible Turbulent Flows," Phys. Fluids, Vol. 7, pp.658-666 (1995).

10 V. Yakhot, S. A. Orszag, S. Thangam, T. B. Gatski, and C. G. Speziale, " "

Phys. Fluids A, Vol 4, 1992.

11 T. H. Shih, J. Zhu, and J. L. Lumely, "A realizable Reynolds stress algebraic

equation model," NASA TM 105993, 1993.

12 Z. Yang, N. Georgiadis, J. Zhu, and T.-H. Shih, "Calculations of Inlet/Nozzle

Flows Using a New k - e Model," AIAA paper 95-2761.

la W. W. Liou and P. G. Huang, "Calculations of Oblique Shock Wave/Turbulent

Boundary-Layer Interactions with New Two-Equation Turbulence Models," Proceedings of the Second Symposium on Transitional and Turbulent Compressible

Flows, Hilton Head, South Carolina, 1995.

14 T.-H. Shih and J. L. Lumley, "Kolmogorov behavior of near-wall turbulence and

its application in turbulence modeling," NASA TM 105663, 1992.

is K.-Y. Chien, "Predictions of channel and boundary-layer flows with a low Reynolds

number turbulence model," AIAA Journal, Vol.20, pp.33-38, 1982.

18 T. H. Shih, W. W. Liou, A. Shabbir, Z. Yang, and J. Zhu, "A New k - e eddy viscosity model for high Reynolds number turbulent flows," Computers and Fluids,

Vol. 24, pp.227-238, 1995.

17 G. J. Micklow, M. R. Harper, and J.M. Deur, "The Effects of Turbulence Modelling on the Numerical Simulation of Confined Swirling Flows," AIAA-93-1976

(1993).
50 W. W. Liou

Model Evaluation

Turbulent Shear Flow

TURBULENT FREE SHEAR LAYER TURBULENT BOUNDARY LAYER

Shock/Turbulent-Boundary-Layer Interactions

(> transonic flow

M<I //_ 1 M<I

BUMP FLOW

supersonic flow

IlVllllllll|ql_pt;I,. ,I.L/IIIIIIIIIIIIIIII

M<I _'- J M<I

I"lls'utlla'*lqnn''z|ilil||lll(L|l|| 1

BUMP-IN-CHANNEL

Iflflll Illlllllllllfllllltfll IIIIIIIllllllllltll _

M<I IIM>I] M<I

I I III I l I 1 4 1 1 I III T llllllllll I

BUMP-IN-CHANNEL

SHOCK REFLECTION RAMP FLOW

Figure 1. Sketches of validation cases
Compressible Turbulent Flow Computations

(a)

51

cl

0.003 I l i

0.002

0.001

0.000

-- k.._model ..--_-----"

,-- 2-scale model /" _"

I---- RNGk-tmodel :" /I-I /"

---- SaZk_mod_ /" / '' /

[] measureme*_t :: [] . / / /

- //

-0.001 , I , I , I

-o_o -o.lo o.oo o.1o

s (m)

0.20

(b)

5.0

4.5

4.0

3.5

3.0

2.5

2.O

I |

1.5

1.0

0.5

-O.2O

-- k._modd

--- 2-scaJe modet .._

.17-

I i I , I

-0.10 0.00 0.10

s (m)

0.20

Figure 2. Comparison of surface measurement for supersonic ramp flow. (a) Skin friction

coefficient; (b) Wall pressure.
52 W. W. Liou

(a)

o_ r _ _, _ I

I i , !s I

i _ is il I

' -'---F, i _',, i i

i ,to i %:._

i

J

i

! i : : :

o° i

t i .......................

-o.oo_

-lO.O -,5.0 0.o 5.0 lO.0

_o

(b)

o

_0

1.0

0,0

-10.0

L__........ i................ i.....t

___i__L__i_/l_._--_;"..--_":.-_-

..... _........._..... i_.__ .....

-_i ', ---- I ........ _......... _................. _ ....

I t

t !............. _...........T........... .,"..........

i " i i '

-5.0 0.0 5.0 10.0

_o

Figure 3. Comparison of surface measurement for oblique shock/turbulent boundary-layer

interaction. (a) Skin friction coefficient; (b) Wall pressure.
Compressible Turbulent Flow Computations 53

O.7O

P]PT 0.50

0._

_10 0.15 0.2O

x(m)

Figure 4. Comparison of wall pressure for ONERA BUMP A.

P/PT

O2

0.5

O.5

0.4

--.-- ¢. ---4. i ]_

• i i

-- _;- .------_: .... _."..... .:...... -__--_..----_.----

i_ _.v_.._., : _ _ : i _ :

--_. -__-_--

_i --..'_St_ _ " i i...... i-- _ i i

0.6 0.7 0.8 o.g 1.0 1.1 1.2 13 1.4 1.5

x/c

Figure 5. Comparison of wall pressure for the transonic shock/turbulent boundary-layer of

Bachalo and Johnson.
54 W. W. Liou

X=I cm

X=4 cm

X=IO cm

X=24 cm

1500.0

1000.0 q

500.0

0.0

3OO0.O

25OO.0

2000.0

1500,0

1000.0

5OO.0

0.0

I t I

, I r I , I t

- %

\_ O'X_,,_

500.0

0.0

3000.0

2.5OO.O

2000.0

1500.0(

1000.0

500.0

0.0

0.0

I ' I I

I , I , I i

Ra d,iI Location, ..o ¢1cm) ,.o

. -

• Measurement

- - - Micklow, Harper, & Deur (1993)

LPDF2D

Figure 6. Comparison of temperature (°K) profiles for the UC Irvine axisymmetric can

combustor.
Compressible Turbulent Flow Computations 55

cl

0.0020

0.0018

0.0016

0.0014

0.0012

0.0010

0.0008

0.0006

0.0

1 I ( I 1

[] ---

l-IVan

LS

Dreist II 1

\ []

20000.0 40000.0 60000.0 80000.0 100000.0

Ree

Figure 7. Comparison of skin friction coefficients for a supersonic flat plate turbulent boundarylayer of M_=2.87.
Page intentionally left blank 
Center ]or Modeling o] Turbulence and Transition

Research Brie]s - 1995

57

Modeling of Turbulent, Reacting

Flows by PDF Methods

A. T. Norris

1. Motivation and Objective

The objective of this work is the development, implementation and validation

of Probability Density Function (PDF) models for turbulent, reacting flows. In

addition, version 1.1 of the LPDF2D code was made available for release, containing

many new features and able to be run in a UNIX environment.

2. Work Performed

In this section, the different areas of PDF methods that were worked on this year

are described. Some validation tests axe also described.

2.1 Introduction.

The basis of Probability Density Function (PDF) methods lies in solving the

transport equation for the joint PDF of quantities of interest, such as velocity,

dissipation, enthalpy and/or composition. The reason why this approach to modeling turbulent reacting flows is attractive is that the chemical reaction is treated

exactly 1 .

Consider the transport equation for species ¢_:

0---f-+ + (1)

where U_ is the velocity, J_ is the scalar flux and S_ is the chemical source term.

Now the chemical source term is a unique function of the set of a species and two

state variables;

S_ = S_(¢1,..., ¢_,p,T), (2)

where p and T are the state variables pressure and temperature respectively. For

all but the simplest of turbulent reactive flows, a solution of Eq.(1) is impossible.

However the problem becomes tractable if the solution of the mean scalar field is

attempted:

Ozi Ozi Oxi + (S_); (3)

where (Q) denotes the mean of Q, and (ui¢_} is the velocity-scalar covariance, an

unknown quantity that needs to be modeled. The other term that needs to be

modeled is the mean reaction source term, (S_). The reason that this needs to be

modeled is because:

S_(¢l, ..., ¢_,p, T), _ S_((¢1), ..., (¢_), (p), (T))). (4)
58 A. T. Norris

It has been shown 1,2 that the two sides in this inequality can differ in both sign

and by orders of magnitude. In addition, the modeling of the reaction source term

would require different formulations for different fuels and mixtures. This would be

clearly impractical.

It is because of the difficulties in modeling the reaction source term that PDF

methods hold such attraction for the solving of turbulent reactive flows.

2.2 Hybrid PDF Method.

The PDF code under development is the LPDF2D solver developed by Hsu et

al3. (1993). This code consists of a particle-based Monte Carlo solver for the

solution of the transport equations for the PDF of the species and enthalpy, coupled

with a finite-volume flow solver for the velocity, turbulence and pressure fields.

Because this method is a combination of PDF and moment closure methods, it

is referred to as a Hybrid PDF scheme. The solver can be used to solve steady,

2D turbulent compressible reacting flows and has been extended to 3D. Work this

year has focused on the refinement of the numerical algorithms, implementation of

several different schemes to calculate reaction rates, porting the code to a UNIX

workstation environment, parallel implementation of the code and some validation

work.

2.3 Numerical Refinement

One of the drawbacks in implementing a PDF method is the extra computer

memory required by the Monte Carlo solver. To help reduce the amount of memory

required, important numerical developments were performed on the PDF solver: a

new time-averaging was developed and implemented, and a convection algorithm

that is accurate even with only a few particles was incorporated in the code.

2.3.1 Time Averaging Scheme.

One of the requirements of the Monte Carlo PDF solver is that it return smooth,

or slowly varying values of the mean scalar to the finite-volume part of the code,

otherwise the flow solver may not converge or could blow-up. The simple way to

achieve this is to increase the number of samples being averaged over. However this

will also increase the computer memory requirements.

To achieve a mean scalar value with a small fluctuation, but with only a few

particles, the technique of time averaging is used. The requirements of a good

time averaging system are that only one set of data need to be stored, that the

influence of samples from an earlier time can be minimized and that computation

of the average is quick and simple. To fulfill these requirements, the time weighted

average scheme was developed. In this method, the weighted time-average of some

mean quantity <¢) at the nth time step is given by:

<8>.- i

w,_ + 1 ((¢)'_ + w"<$)n-1)' (8)

where the over-bar indicates a time averaged quantity and wn is a weighting function,

(6)
Modeling of Turbulent, Reacting Flows by PDF Methods 59

where c,_ is a function of n with a value between 0 and 1. For the case of c,_ - 0,

this scheme reduces to the case of no time averaging. With cn = 1, the average is

formed over all time steps, each equally weighted.

By trial and error it was found that a value of cn = 0.9 provides a good compromise between quick response time and smoothing for the initial stages of a calculation. When the solution of both the scalar and flow fields has stopped exhibiting

transient effects, an increase to cn = 0.999 reduces the scalar fluctuation further,

allowing the solution to converge.

It should be noted here, that the concept of convergence in PDF methods is not

the same as that in pure moment closure schemes. Due to the statistical error

always present, solutions can only converge down to the level fluctuation present in

the mean estimates, rather than to machine accuracy.

2.3.2 Convection Scheme.

Convection in the Monte Carlo PDF code consists of the movement of particles

from one cell node to another due to the action of the mean velocity and turbulent

diffusion. Due to the finite number of particles at each node and the need to move

whole numbers of particles, errors occur in this process. For example, if 5.3 particles

are to be convected, the existing code would convect 5 particles. By increasing the

number of particles at a node, more will be convected and the errors would be

decreased, however more particles means more memory required.

To overcome this problem, we make use of the time-averaging used in the code.

At each step, the number of particles to be convected is calculated. This number

is then rounded to an integer value by a random process, with the probability of

rounding up or down dependent on the value of the fractional remainder. /,From the

previous example, 5.3 has a fractional remainder of 0.3. Thus 5.3 would be rounded

up to 6 with a 30% probability, and rounded down to 5 with a 70% probability.

Because of the time averaging used in the scheme, the mean convection will be 5.3

particles.

Another method would be to save the fractional particles and carry them over to

the next time step. This however has two disadvantages. First the fractional value

needs to be saved, thus using memory. Secondly, this introduces an occilation into

the solution. For example, if 5.1 particles were to be convected, then 5 particles

would be convected except for every 10th step when 6 particles would move. This

period would be damped out by the time averaging to some extent, but not totally.

This small periodic fluctuation could possibly excite some non-physical fluctuations

in the solution, such as vortex shedding. For these reasons, this method was not

implemented.

2.4 Operating Environments.

In order to make the LPDF2D code accessible to a wider range of users, several

changes have been made to the code. The most common computing environment

available these days is the UNIX workstation. Because of this, the LPDF2D code

has been modified to run on workstations. This has involved the removal of machine
60 A. T. Norris

dependent commands, the introduction of a portable random number generator and

the ability to perform calculations in a parallel cluster environment.

2.4.1 Random Number Generator.

To aid the portability of the code, a machine independent random number generator was incorporated into the code. This routine provides a "more random" set of

numbers than most built-in generators, while making use of the 32 bit wrap-around

feature of UNIX machines to speed up the program.

2.5 Chemical Kinetics.

Despite the ability of the PDF method to treat chemical reaction exactly, the

implementation of the numerical chemical kinetics in a Monte Carlo scheme is a

non-trivial matter. First, the solution of the full system of rate equations for the

thermochemistry in turbulent reacting flows requires obtaining solutions for of order 50 chemical species, governed by of order 200 stiff, non-linear rate equations.

At present this task is computationally infeasible, and so reduced mechanisms are

employed. These reduced mechanisms represent the full composition by a few representative species, typically two to five, and the reaction rates of these species are

quasi-global equations derived from the full mechanism. However even the solution

of a reduced mechanism is a computationally expensive task. In a typical PDF

calculation 5, there are on order 100,000 particles and 1,000 time steps, resulting in

the code performing 109 solutions of the reduced mechanism.

A computationally efficient way of reducing this task to manageable proportions

is to use a look-up table. In this table, the reduced mechanism is integrated for

discrete time and composition increments, and stored in a table. Thus integration

is replaced by interpolation, at a considerable saving in computer time. By the use

of adaptive tabulation techniques 4 the size of the tables can also be minimized.

However, because of the many different combinations of fuel and oxidizer, as well

as the differing degrees of complexity of mechanisms, several different options for

chemical reaction have been provided apart from look-up tables. These include

equilibrium chemistry (where the reaction is assumed to proceed infinitely fast), 1

step global reactions and the option to integrate the full mechanism.

In these options, the CHEMKIN 6 package has been used to provide the species

data for the reactions, reducing the amount of input required by the user.

2.6 Parallel Computations.

If the option to use the full mechanism is used, the biggest use of time in the

PDF code is the calculation of reaction rates. Because of this, the LPDF2D code

has been altered to run on a work-station cluster, with PVM 7 message passing. At

present this involves only the calculation of reaction rates in parallel, the rest of the

process being performed in serial. Because of the low volume of message passing,

and the statistical nature of the PDF code, this has proved to be a very effective

method of achieving substantial speedup of the code.

4.0 Validation.

In order to validate the PDF code, and test its performance against that of a
Modeling of Turbulent, Reacting Flows by PDF Methods 61

conventional moment closure code with laminar chemistry, a test case was chosen

and then modeled by both methods. The flow chosen was the CO/H2/N2-air flame

of Ma_ri et al s. This flow was chosen as it has high levels of turbulence in the

reactive region of the flame, thus the turbulent/chemical interactions would have a

significant effect on the flow.

The conclusion of the study was that the PDF method did provide a superior

performance to the moment closure method in predicting turbulent reacting flows.

In Figs. 1 and 2, the temperature contour plots for the two numerical methods are

shown. It can be seen that the PDF method predicts peak temperatures of about

1400K in a narrow band, in good agreement with the experimental data. However

the moment closure results show peak temperatures of 2000K lying in a broad band,

contrary to the experimental results.

Details of this comparison are described by Norris and Hsu 9, and were presented

at the 30th Joint Propulsion Conference in Indianapolis.

5.0 References

1 Pope, S.B., "PDF Methods for Turbulent Reactive Flows," Prog. Energy Cornbust. Sci., 11, 119-192 (1985).

2 Pope, S.B., "A Monte Carlo Method for the PDF Equations of Turbulent Reacting

Flows.," Combust. Sci. Tech., 25, 159-174 (1981).

3 Hsu, A.T., Tsal, Y.L.P. & Raju, M.S., "A PDF Approach for Compressible Turbulent Reacting Flows," Proc. AIAA 31st Aerospace Sciences Meeting and Exhibit, Reno, Nevada, AIAA-93-0087 (1993).

4 Norris, A.T., The Application of PDF Methods to Piloted Diffusion Flames, PhD

Thesis, Cornell University, Ithaca, 1993.

5 Pope, S.B., "Computations of Turbulent Combustion: Progress and Challenges,"

In 23rd Symposium (International) on Combustion, 591-612 (1990).

6 Kee, R.J., Rupley, F.M. & Miller, J.A., "CHEMKIN-II: A Fortran Chemical Kinetics Package for the Analysis of Gas-Phase Chemical Kinetics," Report SAND89-

8009, Sandia National Laboratories, Livermore, CA 94551

7 Geist, A., Beguelin, A., Dongarra, J., Jiang, W., Manchek, R. & Sunderam, V.,

"PVM 3 User's Guide and Reference Manual," Report ORNL TM-12187, Oak

Ridge National Laboratory, Oak Ridge TN 37831

s Masri, A.R., Dibble, R.W. & Barlow, R.S., "The Structure of Turbulent, PilotStabilised Flames of CH4/CO/H2/N2 Fuel Mixtures," Combustion and Flame,

in press (1994).

9 Norris, A.T. & Hsu, A.T., "Comparison of PDF and Moment Closure Methods in

The Modeling of Turbulent Reacting Flows" Proc. AIAA 30th Joint Propulsion

Conference, Indianapolis, Indiana, AIAA-94-3356 (1994).
62 A. T. Norris

.12 I

.11

.lg

.09

.98

X .87

1-4

.#6

, l I i l i i i I l I i I I I "

300.0

--300.0

f

f

Figure I:

coordinate is stretched by factor of five. Contour levels are in degrees Kelvin.

.112 .114 .g6 .118 .ill

X

Contour plots of mean flame temperature produced by the PDF method.

o12 .14

Note that radial

X

-60g. 0 -

• 02 .i4 .06 .g$ .19 .12

X

Figure 2: Contour plots of mean flame temperature produced by the moment closure method. Note that

radial coordinate is stretched by factor of five. Contour levels are in degrees Kelvin.
Center ]or Modeling of Turbulence and Transition

Research Brie]s - 1995

63

Application of PDF to Compressible

Turbulent Reactive Flows

Andrew T. Hsu

1. Motivations and Objectives

In an earlier paper[l], we have introduced a probability density function (pdf)

approach for high speed reactive flows. The extension of this pdf approach to threedimensional applications is reported here. The major obstacle to the application of

a Monte Carlo pdf solver to complex three-dimensional reactive flows is the limitation imposed by the need of a large number of notational particles in a Monte Carlo

simulation to provide a converged, relatively smooth solution free of large statistical

errors, which in general requires enormous amount of both computer memory and

cpu time. To overcome this obstacle, a novel averaging scheme is introduced. This

new averaging scheme allows one to use as few as 50 sample particles per computational cell and yet still provides Monte Carlo solutions that are smooth enough

to be applicable in the finite-volume/Monte-Carlo scheme. Further, the averaging

scheme automatically eliminates the effect of initial conditions_ making it ideal for

large-scale applications.

2. Work Accomplished

In turbulent reactive flow computations, the conventional moment closure models

have difficulty in treating the nonlinear chemical reaction source terms; this problem

is known as the chemical closure problem. One of the methods that can overcome the

chemical closure problem is the probability density function (pdf) method. Several

key advances were made in the development of pdf methods during the past decade.

Most of the earlier advances are primarily in the low-speed combustion area.

In an earlier paper a probability density function turbulence model for compressible reacting flows has been proposed by the present authors [1]. The probability

density function of the species mass fraction and enthalpy is obtained by solving a

pdf evolution equation using a Monte-Carlo scheme. The pdf solution procedure is

coupled w_th a compressible finite volume flow solver which provides the velocity

and pressure fields. A modeled pdf equation for compressible flows, capab!e of capturing shock waves and suitable to the present coupling scheme, has been proposed

and tested. Several 2D supersonic diffusion flames were studied and the results

compared favorably with the available experiment_-data.

In the present paper, the extension of the pdf method to three-dimensional supersonic combustion is reported. The main obstacle to the application of a Monte

Carlo pdf solver to complex three-dimensional reactive flows is the limitation of

computer memory and cpu time. To overcome this obstacle, a novel averaging
64 A. T. Hsu

scheme is introduced. This new averaging scheme allows one with the use as few as

50 sample particles per computational cell to obtain Monte Carlo solutions that are

smooth enough to be supphed to the finite-volume solver as required during each

of the iteration steps of the finite-volume/Monte-Carlo scheme. Furthermore, the

averaging scheme automatically eliminates the effect of initial conditions, making

it ideal for large-scale steady-state applications.

A three-dimensional jet-in-crossflow hydrogen combustion case has been successfully computed using the present scheme. The results are compared with those from

a finite volume reactive flow solver.

2.1 The PDF Method

The modeled pdf equation for the species mass fraction and specific enthalpy used

in the present work can he written as follows:

(pP),,+ (p< > P),j+ = (D,Pj), + M(P)-

where the first two terms on the right hand side of the equation are the modeled terms for turbulent diffusion and molecular mixing; the last term is the term

representing the compressibility effect.

The turbulent diffusion term is simulated using a simple gradient diffusion model

by O'Brien[2] and Pope [3]:

!

- < u i [Yi, h > P = D,P,i.

The use of this gradient diffusion model is consistent with the use of a k - e model

in the flow solver.

The molecular diffusion term is modeled using the continuous mixing model by

Hsu and Chen [4]. This model is an extension of Curl's model. In order to achieve

continuous mixing, we assume that all the particles within a cell participate in

mixing. The extent of the mixing is controlled at the individual particle level. That

is to say, the N particles within a given cell are randomly grouped into N/2 pairs;

the properties of all the particles change according to

Y_(t + St) = AY,_(t) + (1 - A)Y_(t)

Ym(t + 6t) = AY,_(t) + (1 -- A)Y,_(t)

The extent of mixing is controlled at the individual particle level through the parameter A, which is defined as

T

where C' = 2.0. At the hmit 5t -4 O, the above equations becomes

dY_ F

C'_(Y._(t)- Y,_(t)).

T
Application of PDF to Compressible Turbulent Reactive Flows 65

The above equation states that the change of Y,_ due to mixing is proportional to

the difference between Y,_ and Yn, and inversely proportional to the turbulence time

scale, T.

The compressibility effect is modeled by the following expression:

Sp =< p >,_ + < ui >< p >,_ +0.8p < k >< ui >,i -a2pPrM_ + a_peMZt,

which cart be regarded as the convection velocity of a sample particle in the hdirection in the space spanned by h and Yi's. Details of this model are given in Kef.

[1]

2.2 Solution Procedure

A fractional-step Monte-Carlo method developed by Pope [3] is used in solving

the pdf evolution equation. When coupling a pdf solver with RPLUS, a finite

volume flow solver, the species transport equations in the RPLUS code axe no

longer needed. The information we need, at each marching time step, from the flow

solver (RPLUS) includes the mean velocity, pressure, density, and a turbulence time

scale or quantities such as the turbulent kinetic energy and dissipation rate. The

species transport and chemical reactions are simulated by solving the pdf evolution

equation. At every time step, the temperature calculated from the pdf solution is

fed back to the mean flow solver for the computation of heat transfer and pressure.

The Monte Carlo and finite-volume solvers are run in parallel, and information

exchange occurs at every time step until a converged solution is obtained.

2.3 Averaging Scheme

In a Monte Carlo computation, the statistical error is estimated to be of the

order of 1/V_, where N is the number of sample particles per computational cell.

If this error is large, it could lead to an erroneous finite-volume solution when

information from the pdf solver is transferred to the flow solver. Because of this

slow convergence, to obtain a smooth solution from the Monte Carlo solver, one

needs thousands or even tens of thousands of sample particles per cell. In our

previous paper [1] we proposed a combined ensemble-time averaging scheme. That

scheme reduced the required particle per cell, but it requires some more memory

for the time averaging process. For any realistic 3D combustion problems, the

memory requirement of that scheme exceeds the CRAY-YMP capacity. In what

follows we introduce a novel averaging scheme that eliminates the effects of the

initial conditions automatically while requiring no additional memory.

Let an be the combined time-ensemble average at the n th time step, and let b,_

be the ensemble average at the n th time step. A general weighted time average can

be written as

1

- - + l(b + fna -l)

In terms of ensemble average, the above equation can be written as

I fn :

f,_+ 1 (fn + 1)(f,_-1 + 1) bn-l-_ (f,_+ 1)(f,_-i+ 1)(fn-2 + 1)

bn-2+"" 6b_
66 A. T. Hsu

Where fn is a discrete function that can be chosen appropriately to give the desired

weighting for solutions from various time levels. For instance, f,_ - n would give a

equally weighted running average, and the above equation becomes

1

n + 1 (b,, + b,___ +... + b_ + bo)

Fig. 1 shows the weighting distribution as a function of the time step when using

the following function

fn "-- nC

with c = 0.6, 0.8, and 1.0. This weighting procedure gives a fairly smooth Monte

Carlo solution, but it eliminates the effect of only the first few time steps. By

choosing

f. = c(/._1 + 1)

In which case Eq. (9) is equivalent to

1 (b,_ + cb,_-i + c:b,_-2 +'" + c'_-Ibl + c'_bo)

The weighting distribution, as shown in Fig. 2, is shifted drastically towards the

latest solutions. The weakness of this particular formulation is unevenness of the

weight distribution.

Ideally, one would like to eliminate the effect of the initial conditions_ yet have

a fairly even weight distribution for the later solutions. Such a weighting function

can be achieved by replacing the constant c in the equation for f,_ by a variable c,_,

i.e.,

In = + 1)

There axe any number of possibilities in choosing c_. For example, c,_ = 1 - 1/n,

or, = 1 - 1/n::, cn = 1 - 1/e '_, c,_ = a:1/'_, etc. One formulation we have found to be

fairly successful is

1

c_-- I -m

X n

with z > 1. Fig. 3 shows the weight distribution for z -- 1.014, 1.0088, and 1.007.

The advantages of this weighted averaging procedure are that the effect of the initial

condition can be effectively eliminated, and that after a certain number of iterations,

the solutions are equally weighted. For example, in the case of a: = 1.014, the effect

of the solutions from the first 200 time steps are completely eliminated_ and after

about 500 time steps, the solutions are equally weighted in the averaging procedures.

The pdf results reported in what follows axe obtained using this weighting scheme.

2.4 3D Application

In the present work, one of our major objectives is to demonstrate the feasibility

of applying the pdf method to realistic 3D supersonic combustion problems. Unfortunately, there is a lack of experimental data in this area. The case we have chosen
Application of PDF to Compressible Turbulent Reactive Flows 67

to demonstrate the application of the pdf method is a supersonic jet in crossflow.

A 300°K, Mach 1.5 hydrogen jet is vertically injected into a Mach 2 crossflo'_ with

1000°K vitiated air. The jet diameter is 6 mm. Figs. 4-8 are the temperature

and species mass fraction contours from the pdf solution as compared to the solution from a finite volume code (lZPLUS) using laminar chemistry. The grids used

in both calculations are 40 x 15 x 25. The computational domain reaches about

30 jet diameters downstream of the jet. The Monte Carlo simulation of the pdf

equation used 50 sample particles per computational ceU. The coupling between

the pdf solver and the f-mite volume flow solver is done through feeding density,

velocity, and turbulent time scale from the finite volume solution to the pdf solver,

and feeding temperature distribution from the pdf solver to the finite-volume flow

solver.

Fig. 4 shows the central-plane temperature distribution. The results show that

the pdf solution gives a much shorter cold jet core than that from the finite volume

scheme, which is consistent with our previous 2D results. The maximum temperature from both solutions are the same. Fig. 5 shows the temperature contours at

a cross section at about 12 jet diameters downstream of the injection port. The

basic structure of the two solutions are similar, with some differences in details:

The finite volume solution gives a very high temperature region at the bottom wM1

away from the centerline while the pdf solution does not show this. This difference

could be caused by the different boundary conditions used at the wall.

Figs. 6, 7, and 8 show the mass-fraction contours of hydrogen, oxygen, and water

vapor, respectively. The results for the mass fractions are consistent with what

we have observed from the temperature distribution. For instance, the hydrogen

distribution from the finite volume solution shows a maximum of about 0.8 near the

center while the pdf solution there has a lower value of only 0.35 indicating some of

the differences observed between the solutions. In Fig. 8, where it corresponds with

the high temperature at the wall, there is a higher concentration of water vapor

(about 0.35) in the finite-volume solution while the pdf solution shows a maximum

of only O. 12.

Because of the lack of experimental data, it is difficult to judge the validity of

either solution. However, the contours from the Monte Carlo solution is as smooth

as that from the finite volume solution, showing that it is feasible to make large

scale computations on practical problems using the Monte Carlo pdf solver.

2.5 Concluding Remarks

A pdf method for high speed turbulent combustion has been successfully extended

to 3D applications. A novel averaging scheme has been introduced which made the

applications of the present pdf method to reahstic large scale computations feasible.

3. Future Work

3.1 PDF

Parallel computing and NOx prediction will be the two major areas for future

PDF applications. The goal is to make the PDF method an industrial design tool.
68 A.T. Hsu

4. References

1 Hsu, A.T., Tsai, P.Y. and Raju, M.S., " A PDF Approach for Compressible

Turbulent Reacting Flows," AIAA 93-0087, January, 1993, Reno, NV.

2 O'Brien, E.E., "The Probability Density Function Approach to Reacting Turbulent Flows." Turbulent Reacting Flows, Springer-Verlag, Berlin, P. 185, 1980.

3 Pope, S.B., " A Monte Carlo Method for the PDF Equations of Turbulent Reactive Flow." Combustion Science and Technology, vol 25, p.159, 1981.

4 Hsu, A.T. and Chen, J.Y., "A Continuous Mixing Model for PDF Simulations and

its Applications to Combusting Shear Flows." Proceedings of the 8th Symposium

on Turbulent Shear Flows, Munich, September, 1991.
Application of PDF to Compressible Turbulent Reactive Flows 69

0.012

0-010 t

0.008

r-

•-_ 0.006

0.OO4

0.002

0.000

0

C

:.-c2(fn_ , 4- 1 )

C2=0"99 -_

C2= 1.0C2=0"_._

_00 800

time step

Fig.2. Weightdistribution forthe firstI000time

stepsusingEq. (12).

0.0020

0.0016

"_ 0.0012

._

0.0008

0.0O04

C1=1.0

0.0000 ' r

0 _00 8OO

time step

Fig. I. Weightdistribution forthefirstI000time

stepsusingEq. (11).

0.005

0.004

- 0.003

._

0.002

0.001

0.000

0

fn=cn(f._l + 1)

cn=l-1/x n x-'1.007-_

x= 1.00_/" x=l .014 Y//

_00 800

time step

Fig.3. Weightdistribution forthe first1000time

stepsusingEqs.(14)and (15).
70 A. Hsu

TEMPERATURE CONTOURS

PDF SOLUTION

RPLUS SOLUTION

Fig. 4. Center plane temperature distribution for the

3D supersonic jet combustion problem.

ZOO

TEMPERATURE CONTOURS

RPLUS SOLUTION__ PDF SOLUTION

1 _0"" 800

x-12d

Fig..5. Temperature distribution at x=12d for the

3D supersonic jet combustion problem.
Application of PDF to Compressible Turbulent Reactive Flows 71

1-12Mass Frac_on

RPLUS SOLUllON __ PDF SOLUTION

.-....\

Fig. 6. Hydrogen distribution at x=12d for the 3D

supersonic jet combustion problem.

a_ MI_ Fraclkm

RPLtI_SOLUllON PDF $Ott/IION

Fig. 7. Oxygen distribution at x=12d for the 3D

supersonic jet combustion problem.

H_OMu= _m

RPLU_:SO(.U'IIOfl --'_ POFSOLUTION

at

Fig. 8. Water vapor distribution at x=12d for the

3D supersonic jet combustion problem.
Page intentionally left blank 
Center for Modeling of Turbulence and Trar_sition

Research Briefs - 1995

73

Turbulence Model Development for NPARC

J. Zhu, Z. Yang, and T. H. Shih

1. Motivation and Objective

The purpose of this researchis to enhance the turbulence modeling capability of

NPARC. The outcome of this research is a turbulence subprogram which is independent of NPARC but linked to NPARC in code execution. In the solution process,

the turbulence subprogram interacts with the NPARC code, the mean field solver,

to produce the solution. This turbulence subprogram contains a number of turbulence models suitable for aerospace and aero-propulsion system applications. These

models are chosen from the state of the art of turbulence models developed by both

CMOTT researchers and researchers from other groups of the turbulence modeling

community. The turbulence subprogram with advanced turbulence models serves

the following purposes: first, it will enhance the turbulence modeling capability of

NPARC; second it will provide a vehicle for incorporating future turbulence models, and third it can be used to validate turbulence models against complex flows of

industry interest.

NPARC is a general purpose computational fluid dynamics (CFD) code. It solves

Navier-Stokes equations in a conservation form in a general coordinate system. It

can handle complex geometries and different types of boundary conditions. The

NPARC code is widely used in the aerospace and aero-propulsion community. Up

to now, many flows of practical interest have been calculated using the NPARC

code, ranging from an inlet flow in an air-breathing engine to the flows around a

maneuvering airplane. A description of the code and its numerical scheme is given

by Cooper and Sirbaugh 1. Updates of code capabilities and a listing of bibliographies of using the NPARC code to calculate flows of practical interests can be found

in the current version of the NPARC User Guide 2. The NPARC code is actively

supported by NASA Lewis Research Center via the NPARC Alliance, a partnership

formed between NASA Lewis Research Center and Air Force Arnold Engineering

Development Center. The NPARC users are represented by the NPARC Association

currently headed by Boeing Company and McDonnell Douglas Company.

2. Work Accomplished

To enhance the turbulence modeling capabilities of NPARC, a turbulence module

approach is adopted in the present project. In the following, the turbulence module

for NPARC is briefly described. We will then list the turbulence models that are

available in the current version of the turbulence module. NPARC calculations

using the turbulence module for some propulsion flows will also be presented.

The work accomplished is done by J. Zhu, Z. Yang, T. H. Shih at Center for

Modeling of Turbulence and Transition (CMOTT), and by N. Georgiadis at NASA

Lewis Research Center. This project is coordinated by D. R. Reddy of NASA Lewis

Research Center.
74 J. Zhu, Z. Yang, and T. H. Shih

2.1 Introduction to Turbulence Models in NPARC

Most of the flows in aerospace and aero-propulsion systems are turbulent because

of the high Reynolds numbers of the flows. Thus, for a code with a well developed

numerical scheme such as the NPARC code, its turbulence modeling capability is

the pacing item for computational accuracy. Currently, NPARC contains a number of mixing length eddy viscosity models 3,4 and the Chien k - E eddy viscosity

model 5. A mixing length eddy viscosity model has the advantage of being simple in

implementation and fast in execution; however it is a 'one physics/one flow' model,

i.e., the mode is tailered to a particular physics/flow, and its performance deteriorates drastically when used to calculate other flows. In comparison, a k - e eddy

viscosity model is much more general than mixing length eddy viscosity models

and could be used for a much wider range of flows, for example, both free shear

flows and wall-bounded flows. Among the existing low-Reynolds number k - _ eddy

viscosity models, the Chien k - e eddy viscosity model is perhaps the simplest and

gives a reasonable performance for attached boundary layer flows. However, it has

some well known deficiencies. For example, its performance is rather poor for flows

with transition, separation and reattachment, because y+ is used in modeling the

near-wall effect. These deficiencies can be removed and the applicability of k - c

model can be substantially expanded. Recently, a CMOTT k - E eddy viscosity

model has been implemented into the NPARC code. This model does not contain

the friction velocity, unlike the Chien k - _ model, and can thus be used for flows

with separation and reattachment. The detail of the model form and the results of

model performance will be presented in later subsections.

2.2 CMOTT Approach - Turbulence Subprogram for NPARC

In the effort of enhancing the turbulence modeling capabilities of NPARC, we

implement different turbulence models in a stand-alone turbulence subprogram (also

referred to as turbulence module) instead of implementing them in the NPARC code

directly. We feel that this is a preferred approach in a situation that a gap exists

between turbulence model developers and CFD users. The former mainly use simple

flows to verify new modeling concepts and evaluate the resulting models, while the

latter are usually reluctant to implement new turbulence models unless they see

that such models have shown a good performance for a wide range of complex flows

relevant to industry situations. In order to bridge this gap, we develop turbulence

modules for industry's CFD codes. Since the NPARC code is a major industry code

with a well established numerical scheme and a large user group in the aerospace

and aero-propulsion community, we use the NPARC code as our reference CFD

code to write the turbulence modules. With minor modifications, the turbulence

module written for the NPARC code can be adopted for other industry CFD codes.

The turbulence module is written in a self-contained manner so that the user can

use any turbulence model in the module without worrying about how it is implemented and solved. The inputs to the turbulence module are mean flow variables,

boundary and geometric information which are provided by a mean flow solver. The

outputs of the turbulence module are the Reynolds stresses and heat fluxes which

are needed for the mean flow calculation. In the current formulation, the Reynolds
Turbulence Model Development for NPARC 75

stresses and heat fluxes are expressed in terms of turbulent diffusivities and relevant turbulent source terms. The interaction between the mean flow solver and

the turbulence module gives the final turbulent flow solution. With the aid of the

turbulence module, we can take the advantage of the well-established and sophisticated CFD codes to test turbulence models for a variety of complex flows which are

intractable with the simple research codes. The turbulence module, on the other

hand, enhances the turbulence modeling capabilities of existing CFD codes and

provides a vehicle to facilitate the technology transfer from the model development

to model applications.

In the turbulence subprogram for NPARC, the governing transport equations

for turbulence quantities (k and e) are written in conservation form and a finitevolume procedure is used to solve these equations. Comparing with the NPARC

code, the turbulence subprogram has the following features: (a) the cell-corner

arrangement is used, i.e., the flow variables are stored at grid nodes rather than at

the centers of control volumes; (b) contrary to the NPARC code and most other

codes for compressible flows, the non-delta form of equations is used which leads

to a simpler linearization and is more effective to ensure the positiveness of the

turbulent kinetic energy and its dissipation rate (k and e) than the delta form.

To reduce numerical diffusion while maintaining necessary stability, the secondorder accurate HLPA scheme 6 is used for the convective terms of the turbulent

transport equations. This scheme is implemented such as to blend two schemes by

a parameter )_, with limiting values ,k=0 for the (first-order) upwind and _=1 for the

second-order and bounded solution scheme. The turbulence subprogram includes

four most commonly used boundary conditions - the inflow condition, the outflow

condition, the symmetry plane condition, and the condition on the solid wall. The

turbulence subprogram also allows the user to introduce boundary conditions other

than these four types. Care has been taken to make the turbulence module userfriendly. This is achieved by providing a few parameters that the user can define in

model selection and problem specification. The resulting algebraic set of equations

are solved iteratively via a under-relaxation procedure. The user can also control

the relaxation parameters in the problem specific part of the subprogram. The

turbulence module is written in a similar manner as the FAST2D code written by

Zhu 7. A detailed description of the turbulence module for the NPARC code is

going to be presented in a paper by Zhu and Shih s in the upcoming 31st AIAA

Joint Propulsion Conference.

2.3 Turbulence Models in the Subprogram

Currently, four turbulence models are available in the turbulence subprogram.

They are the mixing length eddy viscosity model of Baldwin-Lomax 3, the low

Reynolds number k - e eddy viscosity model of Chien 5, the low Reynolds number k - e eddy viscosity model of Shih and Lumley 9, and a new low Reynolds

number realizable k - e eddy viscosity model developed here at CMOTT.

Because of the increasing computer power available and the increasing demand

for more accurate flow prediction, it is recommended that a k - e eddy viscosity

model is used for calculating flows of engineering interests. The Baldwin-Lomax
76 J. Zhu, Z. Yang, and T. H. Shih

mixing length eddy viscosity model is included in the turbulence subprogram for

two reasons. First, it is chosen as a representative of the mixing length models.

Second, it can be used to generate the initial field for k - e model calculations.

The Baldwin-Lomax mixing-length model, the Chien k - e model, and the ShihLumley k - e model have been published in archival journals. Here, only the newly

developed low Reynolds number realizable k - e eddy viscosity model is described.

The model has the same model form as the Shih-Lumley low Reynolds number

k - e eddy viscosity model, i.e, it has the same form of transport equations for

k and e, it uses the same near-wall boundary conditions for k and e as derived

from the Kolmogorov behavior of the turbulence field near the wall, and it uses

the same parameter Rv - ykl/2/v to construct the damping function f_,. The

boundary conditions based on Kolmogorov behaviors of near wall turbulence makes

the computing robust. The parameter Rv is a better choice for modeling of near

wall turbulence than y+ since it does not involve the skin friction u_. The major

improvement in the present realizable k - e eddy viscosity model compared with the

Shih-Lumley model is that a realizable formulation for C_, is used. In an isotropic

eddy viscosity model, the Reynolds stresses are related to the mean velocity field

by

2 2

-u uj = vr(v ,j + - - (1)

In the framework of k - e two equation eddy viscosity model, the eddy viscosity VT

is further expressed as

=c./. (2)

where the damping function f_ is introduced to account for the effects of wall.

Traditionally, a constant value of C_, = 0.09 is used in equation (2). However,

any two-equation eddy viscosity model with a constant C_ is unrealizable in the

sense that the turbulent component energy given by the model can become negative in certain flow situations. In addition, data from direct numerical simulations

and from experiments suggest that C_, vary with flow situations. For example, to

match the data, different C_ values should be used for homogeneous shear flows

and for turbulent flows in a two dimensional channel. In this study, a variable C_

formulation is used, in which C_ depends on the solutions of the mean flow field

and the turbulence field. The C_ formulation was derived based on the realizability

constrains, rapid distortion theory (RDT), and the invariants theory and takes the

following form

1

C_, = Ao + A,U* k_" (3)

E

In equation (3), k/e represents the characteristic time scale associated with the

turbulent motion. U* is related to the characteristic time scale of the mean field
Turbulence Model Development .for NPARC 77

and it is given by

= &j - -Uk,k

(4)

where S_j is the mean strain tensor, fli*j is the mean rotation tensor in the inertia

frame, and wk is the angular velocity of the rotating reference frame. The mean

rotation tensors in the inertia frame and in the rotating frame are related by

9ti_ = Qij - eijkwk. (5)

In a rotating frame, the mean strain tensor Sij and the mean rotation tensor Q/j

are related to the mean velocity field by

1 1 U (6)

It is worth noting that the effect of frame rotation is included in the present model

directly via its effect on U*, and consequently on C u and VT.

The parameter As in equation (3) is a scalar function formed from the invariants

of the mean velocity field. It has the following expression

As = v_cos¢,

W - SijSjkSm,

S _

1

¢= -garccos(v w)

S = v/_jS_.

(7)

As shown in Reynolds t°, the choice of As guarantees the correct limiting behavior

of the model in the rapid distortion limit.

The parameter A0 in equation (3) is a model constant. In the present study,

Ao = 4.0 is chosen so that the present C u formulation gives a Cu value of 0.09

for the equilibrium log layer regions of two-dimensional turbulent boundary layer

flows. It was based on the experimental data for these flows that a constant value

of C_, = 0.09 was first suggested.

Equations (3) -- (7) give the Cu formulation that we use in the present study.

For such a C u formulation, the resulting Reynolds stresses are always realizable in

that no physically un-admissible Reynolds stresses will be produced. In addition,

the variations of Cu with flow situations are in the right direction, i.e., the value of

Cu is smaller than 0.09 in regions where it should be. Thus, it can be expected that

the current formulation of C_, will improve the performance of k - e eddy viscosity

model considerably, particularly in regions where strain rate is large.

In a paper by Shih et al.11, the above mentioned Cu formulation was used in

conjunction with a new (high-Reynolds number) k - e eddy viscosity model, in

which the dissipation rate equation was derived based on the asymptotic behavior of the fluctuating vorticity equation in the limit of large Reynolds numbers.
78 J. Zhu, Z. ]rang, and T. H. Shih

The paper shows that the C_ formulation helps the model performances for rotating homogeneous shear flows, free shear flows, and wall boundary flows. For wall

bounded flows, the near-wall boundary conditions were provided by the standard

wall function formulation presented in Launder and Spalding 1_.

2.4 Flows Calculations using NPARC and Turbulence Module

A series of flow calculations have been made using the NPARC code to solve the

mean flow field and the turbulence module to solve the turbulence model equations.

The turbulence module was first validated against the turbulent boundary layer

flow over a flat plate with zero pressure gradient. Complex flows of nero-propulsion

interest were then calculated. These include an ejector nozzle flow which provides a

good test for flow mixing; a transonic diffuser flow which tests shock wave/turbulent

boundary layer interactions; and a boat-tail nozzle flow which tests flow separation

and jet-boundary layer interactions. In each case, the numerical convergence and

the grid sensitivity are checked to make sure that the numerical solutions obtained

are creditable. In the following, each flow case is briefly reported. These calculations

are going to be presented in a paper at the upcoming 31st AIAA Joint Propulsion

Conference by Yang et al.la.

2.4.1 Flow over a fiat plate

The purpose of calculating this flow is for code validation. The computational

domain is shown in figure 1. The total grid points are 111 by 81 with 111 grid

points in the x (streamwise) direction and 81 grid points in the y direction. In the

x direction, 14 grid points are located before the leading edge of the flat plate. In

the y direction, grid points are stretched so that the first off wall points have an

average y+ value of about 1. Since the NPARC code is for compressible flows, and

computational results are going to be compared with an incompressible turbulent

boundary layer, a freestream Mach number of 0.2 is used in the present calculation.

Although straightforward on the surface, the calculation of this flow is quite

demanding in that the flow starts from a laminar-like flow, then undergoes a transition process, and then finally settles down to flow of the turbulent boundary layer

type. It is very difficult to capture all three flow regions accurately. In the present

computation, we are mainly interested in the solutions in the turbulent boundary

layer region. Figure 2 shows the skin friction coefficient distribution as a function

of Reo, the Reynolds number based on the momentum thickness of the boundary

layer. Both the model predictions and the experimental results of Wieghardt and

Willmann 14 are shown for comparison. It is seen that all models give similar results

and all model predictions are close to the experimental data. Comparisons were

also made between solutions of the NPARC code with built-in turbulence models

and solutions of the NPARC code with the turbulence module. These solutions are

found to be the same, thus the code with the turbulence module is validated.

2.4.2 Flow in an ejector nozzle

Ejector nozzle flows have important applications in propulsion systems because

ejector nozzle is a candidate device to reduce the jet noise from the engine. The
Turbulence Model Development for NPARC 79

noise reduction is particularly relevant to the High Speed Research (HSR) Program

and the High Speed Civil Transport (HSCT) Program in NASA. A decade or two

ago, ejector nozzle was also studied for its applications in the vertical land/taking

off aircraft.

The geometry for an ejector nozzle is shown in figure 3. The performance of

the ejector nozzle is measured by its mixing properties. To see the mixing effect,

computed velocity and temperature profiles at a few stations along the streamwise

direction are shown in figure 4 and figure 5, respectively, along with the experimental

data of Gilbert and Hill 15 . It is seen that different k - e eddy viscosity models give

similar results and that results from k - e eddy viscosity models are much better

than those from the mixing length eddy viscosity models. The same conclusion was

also reached by Georgiadis and Yoder 16 who calculated this flow using the NPARC

code with its built-in turbulence models.

2.4.3 Flow in transonic diffuser

Transonic diffuser flow is another flow of interest in propulsion systems. In the

present calculation, the experiment by Salmon, Bogar, and Sajben t_ is chosen as

the test case. Depending on the pressure at the nozzle exit, there can be three

different kinds of flow patterns. They are the no-shock case, the weak-shock case,

and the strong-shock case. All these three cases were calculated in the present

study, however, only the results for the strong-shock case, which is also the most

difficult case, are shown here.

Figure 6 shows the diffuser geometry, together with grid information. The pressure distribution along the upper wall of the diffuser is shown in figure 7. Both the

experimental results and the results from different model predictions are shown.

For this flow, the critical parameters that are of interest to the design engineers

are the location of the shock wave and the pressure recovery after the shock wave.

Among the model predictions, the new realizable k - e eddy viscosity model gives

the best performance in these critical parameters. Figure 8 shows the velocity profile at X = 2.88, a station located after the shock wave where the experimental

data are available. The realizable k - e eddy viscosity model gives a much better

prediction for both the overall feature of the field and the size of the separation.

Since the controlling mechanism for this flow is the shock wave/turbulent boundary layer interaction, we can conclude that the new realizable k - e eddy viscosity

model captures this mechanism much better than the k - e eddy viscosity model

with constant Ct, , such as the Chien k - e eddy viscosity model. The Chien k - e

model is itself much better than the Baldwin-Lomax mixing length eddy viscosity

model for this flow.

2.4.4 Flow around a boat-tail nozzle

Boat-tail nozzle is a device commonly found in propulsion systems. A particular

interest of the present study is on the effect of the jet plume on the turbulent

boundary layer formed on the nozzle afterbody. Depending on the flow situations,

the boundary layer can be either separated or attached. Boat-tail nozzle flow also

provides a good test case for the two dimensional version of the NPARC code in
80 J. Zhu, Z. Yang, and T. H. Shih

that it tests the code's capability for axisymmetric flows.

The experimental studies of Mason and Putnam is, Abeyounis and Putnam 19 are

chosen to validate the model calculation. In these studies, the freestream Mach

number is 0.8. The stagnation pressure in the nozzle is such that the jet exits the

nozzle at the local sonic speed, i.e. at Mach number of one. The nozzle geometry

and the computational domain is shown in figure 9. Figure 10 shows the pressure

distribution on the nozzle afterbody from both the computation and the experiment.

It is seen that both the Chien k - e model and the new realizable eddy viscosity

model give similar results, and they are all fairly close to the experimental data.

The k-e eddy viscosity models give improved predictions compared with the mixing

length eddy viscosity model.

3. Future Plans

Turbulence model development for NPARC is a long term activity. What was

reported above represents the first part of this project. Many other activities need

to be carried out in this project. For example, the numerical capability of the

turbulence module needs to be enhanced; more advanced turbulence models need

to be incorporated into the turbulence subprogram; and NPARC calculations with

the turbulence modules need to be carried out for a large range of flows relevant

to propulsion systems. To this end, a research team has been formed. In addition

to the CMOTT members for this project (J. Zhu, Z. Yang, and T.H. Shih), this

research team also includes N. Georgiadis and D.R. Reddy of NASA Lewis Research

Center, and J.R. Sirbaugh at NYMA Inc. The activities of this research team has

been incorporated into the NPARC Development Plan of the NPARC Alliance. The

following are some of the planned highlights of these activities.

3.1 Turbulence Module for 3D Flows

Currently, the turbulence module is written for 2D flows (both planar flows and

axisymmetric flows.) Two dimensional flows provide a basis for code and model

validation. However, a code for 2D flows is only the first step in the coding process

since almost all practical flows in aerospace and aero-propulsion systems are three

dimensional. Thus, the next step is to extend the turbulence module to three

dimensional flows and to test the turbulence module and models accordingly. It is

expected that a 3D version of the turbulence module will be available by Summer

1995.

3.2 Tensorial and Galilean Invariant Eddy Viscosity Model

So far, all the k - e eddy viscosity models in the turbulence module depend

on geometry information in their near wall treatment, e.g. the distance from wall

appears in the damping function. Models with geometry information are undesirable

for flows with very complex geometries because in those situations, the definition of

wall distance can be ambiguous. In addition, models without geometry information

can be easily incorporated into CFD codes using unstructured grids while models

with geometry information can not. To propose such a tensorial invariant model, it

is necessary to use parameters other than Rv or y+ to calibrate the near wall effect.
Turbulence Model Development for NPARC 81

In the paper by Jones and Launder e°, in which the first low-Reynolds number k - e

eddy viscosity model was proposed, Rt, the turbulence Reynolds number, was used.

In a more recent paper by Yang and Shih 21, a new parameter, R = k/(Su), was

introduced to model the near wall effects. Both Rt and R are field variables and are

candidates to construct tensorial invariant models for near wall turbulence. These

two parameters will be tested in the turbulence module for the NPARC code for

complex flow calculations.

3.3 Wall Function Features for 1NPARC

An alternative to the low-Reynolds number turbulence model approach, in which

the integration is carried down to the wall, is to use the wall functions which provide

boundary conditions for the solutions at the equilibrium log layer of a turbulent

boundary layer instead of the wall. The wall function approach has the appeal that

it reduces the near-wall grid resolution requirements considerably, and consequently

reduces the numerical stiffness as well. However, the traditional wall functions were

developed based on the turbulent boundary layer flow over a flat plate at zero

pressure gradient. It does not contain the pressure gradient effect and it can not be

used for flows with separation, a common feature for flows in propulsion systems.

Better wall functions have been developed by Yang and Shih 22 and Shih et al.23

These new wall functions take into account the effect of the pressure gradient. In

addition, the paper by Shih et al.23 also addresses the issue of how to implement

wall functions to calculate flows with separation and reattachment. These new wall

functions are going to be incorporated into the turbulence module for the NPARC

code. We expect that the wall function feature would be available by the end of

1995.

3.3 Anisotropic Algebraic Reynolds Stress Model for NPARC

All the turbulence models mentioned above are isotropic eddy viscosity models

which assume a linear relationship between the Reynolds stresses and the mean

strain rate. This assumption is not valid for complex flows. Anisotropic Reynolds

stress models overcome this deficiency and will give a better performance than the

current isotropic eddy viscosity models for flows where anisotropy is important, for

example, the stagnation point flow and flows with massive separation. Recently, one

such anisotropic algebraic Reynolds stresses model was developed by Shih et al, 24

and has shown improved predictions for a large range of complex flows including

the confined jet flow and the swirling jet flow. This anisotropic algebraic Reynolds

stress model is going to be implemented into the turbulence module for the NPARC

code, and the resulting turbulence module will be used together with the NPARC

code to calculated a large range of practical flows. To implement this model, the

NPARC code for the mean flow solver has to be modified because in addition to

the eddy viscosity term, the anisotropic Reynolds stresses model also brings other

source-like terms into the mean momentum equations. The planned milestone for

this work is the end of FY96.
82 J. Zhu, Z. Yang, and T. H. Shih

4. References

1 Cooper, G.K. and Sirbaugh, J.R., PARC code: theory and usage, AEDC-TR-89-

15, Arnold Engineering Development Center, 1989.

2 NPARC Alliance Technical Team, A User's Guide to NPARC, Version 2.0, 1994.

a Baldwin, B. S. and and Lomax, H., "Thin layer approximation and algebraic

model for separated turbulent flows," AIAA Paper 78-256 (1978).

4 Thomas, P. D., Numerical method for predicting flow characteristics and performance of nonaxisymmetric nozzles, Theory, NASA CR-3147, 1979.

5 Chien, K.Y., "Predictions of channel and boundary-layer flows with a low Reynolds

number turbulence model," AIAA J. 20 pp.33-38, (1982).

6 Zhu, J., "On the higher-order bounded discretization schemes for finite-volume

computations of incompressible flows," Comp. Meths. Appl. Mech. Eng. 98,

pp.345-360 (1992).

7 Zhu, J., FAST-2D: A computer program,for numerical simulation of two-dimensional

incompressible flows with complex boundaries, Rep. No. 690, Institute for Hydromechanics, University of Karlsruhe, 1991.

s Zhu, J. and Shih, T.H., "A turbulence module for the NPARC code," to be

presented at the AIAA Joint Propulsion Conference, San Diego, July (1995).

9 Shih, T.-H., and Lumley, J.L., "Kolmogorov behavior of near-wall turbulence and

its application in turbulence modeling," Comp. Fluid Dyn. 1, pp.43-56 (1993).

10 Reynolds, W.C., Fundamentals of turbulence for turbulence modeling and simulation, Lecture Notes for Von Karman Institute, AGARD Report No. 755, 1987.

11 Shih, T.-H., Liou, W.W., Shabbir, A., Yang, Z. and Zhu, J., "A new k-e eddy

viscosity model for high Reynolds number turbulent flows - Model development

and validation," J. of Computers and Fluids, 24, no. 3, 1995.

12 Launder, B.E. and Spalding, D.B., "The numerical computation of turbulent

flow," Computer Methods in Applied Mechanics and Engineering 3, 269 (1974).

13 Yang, Z., Georgiadis, N., Zhu, J. and Shih, T.H., "Calculations of inlet/nozzle

flows using a new k-e model," to be presented at the AIAA Joint Propulsion

Conference, San Diego, July (1995).

14 Wieghardt, K. and Willmann, W., "On the turbulent friction layer for rising

pressure," NACA TM-1314 (1951).

15 Gilbert, G. B. and Hill, P. G., "Analysis and testing of two-dimensional slot

nozzle ejectors with variable area mixing sections," NASA CR-2251, 1973.

16 Georgiadis, N. J. and Yod.er, D. A., "Use of Navier-Stokes methods for calculation

of high speed nozzle flow fields," AIAA Paper 94-3212, June 1994.

17 Salmon, J. T., Bogar, T. J., and Sajben, M., "Laser Doppler velocimeter measurements in unsteady, separated transonic diffuser flows," AIAA J. 21, 1690

(1983).
Turbulence Model Development for NPARC 83

is Mason, M.L., and Putnam, L.E., "Pitot pressure measurements in flow fields

behind circular-arc nozzles with exhaust jets at subsonic free-stream Mach numbers," NASA TM-80169, December 1979.

19 Abeyounis, W.K., and Putnam, L.E., "Investigation of the flow field surrounding

circular-arc boattail nozzles at subsonic speeds," NASA TP-1633, May 1980.

20 Jones, W.P., and Launder, B.E., "The calculation of low-Reynolds number phenomena with a two-equation model of turbulence," Intl. J. of Heat and Mass

Transfer 16, 1119 (1973).

21 Yang, Z. and Shih, T.H., "A Galilean and tensorial invariant k- e model for near

wall turbulence," AIAA Paper 93-3105 (1993).

22 Yang, Z. and Shih, T.H., "On the wall functions for two equation turbulence

models," presented in 42nd APS/DFD meeting, Albuquerque, New Mexico, 1993.

23 Shih, T.-H., Yang, Z., Zhu, J., and Shabbir, A., "Surface flows in a turbulent

boundary layer", In preparation.

24 Shih, T.-H., Zhu, J. and Lumley, J.L., "A new Reynolds stress algebraic equation

model," NASA TM 106644 (1994), also to appear in Computer Methods in Applied

Mechanics and Engineering.
84 J. Zhu, Z. ]Pang, and T. H. Shih

SUP WALL NO WALL I

Figure 1: Configuration for turbulent boundary layer flow over a fiat plate

0

0

o 4.0

2.0

-- Boldwin_Lomox

.... Chien

.... Shih

-- CMOTT-Cmu

0.0 , I I I ,

0.0 5.0 10.0 15.0

Re_theto/1000

20.0

Figure 2: Skin friction distribution for the turbulent boundary layer flow
Turbulence Model Development .for NPARC 85

Nozzle exit height= 0.122.In.

SECONDARY

__. -,____.... ____-4 1_87in.

Figure 3: Configuration for an ejector nozzle flow

1.0 ' ' • i t • • " t " " ' t ' " '

1

0.5. x=7.0in

-0.5o.... '1.

-1.0 ''

0.0 0.2 0.4 0.6 0.8 1.0

U

• - ' t "U t • • ' I • ' " I • • •

x-!O'_

, o Exp.

_o. _7-_,c. Chien

, , , , ,CMOT'I

0.0 0.2 0.4 0.6 "0.8 1.0

U

Figure 4: Velocity distributions in the ejector nozzle flow

ORIGINAL PAGE IS

OF POOR QUALITY
86 J. Zhu, Z. }rang, and T. H. Shih

.0 ' " ' * ' 0 _ q ! , * , • I _ 0 0 i

0.5

0.0 , .,,,,,i, i,[! ;,, '

-0.5

-1.0

0.95 1.00 1.05 1.10 0.95

To

=10.5in

.,:_ o Exp.

_ Chien

.... [o.,, :::, c_omT',,

1.00 1.05 1.10 1.15

To

Figure 5: Temperature distributions in the ejector nozzle flow

/THROAT HEIGHT (HT) = 44 mm

/

•i/Z/Z/-., /....//i/i//iiliililiiii/ "

,.,.,o. +.

Figure 6: Configuration for an transonic diffuser flow
Turbulence Model Development for NPARC 87

10

0.8

o

._0.6

Q.

0.4-

Pressure on the upper woll

0 EXP

.... Boldwin-Lomox

0 _ .... Chien

CMOTT

0/0' , 0

I ,

• s

0.2

-5.0 0.0 5.0 iO.O

Figure 7: Pressure distribution on the upper wall of the diffuser

1.5

10

0.5

0.0

-0.2

Velocity ot X = 4.6

i i i i i

°°o0"0 - . . .

0 EXP ,' ib

.... Boldwin-Lomax ' /

.... Chlen

-- CMOTT , 0

, °

' 0

o

,\ o

.'¢o

0.0 0.2 0.4 0.6 0.8

Figure 8: Velocity prof-_les after the shock-wave
88 J. Zhu, Z. }rang, and T. H. Shih

Figure 9: Configuration for a boat-tail nozzle flow

0.4

0.2

0.0

-0.2

(_EXP

.... Chien

CMOTT

i

_f

o

o/

/"

-0.4 r i T L

0.0 0.2 0.4 0.6 0.8

Figure 10: Pressure distribution on the after-body of the boat-tail nozzle
Center for Modeling of Turbulence and Transition

Research Briefs - 1995

89

Calculations of Turbulent Reacting

Flows in Can Combustors

A. T. Norris and W. W. Liou

1. Motivation and Objective

A newly developed probability density function (PDF) code is applied in the

numerical simulation of swirling can combustor flows.

Due to enviromental concern over NOx emissions, the design of future gas turbine combustors, operating at high pressures and temperatures, requires a detailed

knowledge of the combustion processes in the different operating conditions of an

engine. Therefore, robust and accurate computational modeling of the comubstor

flow fields are necessary in the design of efficient, low emission combustors.

In modeling turbulent reactive flows, PDF methods have an advantage over the

more traditional moment closure schemes. PDF methods treat the chemical reaction

source term in an exact manner, while moment closure schemes model the mean

reaction rate. For example, the commonly used laminar chemistry approximation

is based on the assumption that the effects of turbulence on the chemistry are

negligible. This, or similar assumptions fail for flows with high turbulence levels

and finite rate chemistry.

2. Hybrid PDF Method

The hybrid PDF scheme used is a combination of moment closure and PDF

methods. The continuity and momentum equations are solved by a moment closure

method while the transport equation for the joint PDF of composition is used to

solve the scalar fields. This combination has the advantage in that the PDF part

of the code can be added to an existing moment closure method with a minimal

amount of coding. Thus existing combustion codes can be extended to include the

PDF model of turbulent chemical reactions without the need for extensive rewriting.

3. Experiment

The flow chosen to test the new code, is the axisymmetric can combustor (ASCC)

tested by UC Irvine.[1] This combustor is cylindrical, with an 80 mm ID and is 320

mm long. A sketch of the combustor is shown in Fig. 1. The ASCC is operated

at 1 atmosphere pressure. The fuel is gaseous propane, which is injected through

an annular cone nozzle with an inner radius of 5.4 mm and an average injection

velocity of 12.4 m/s. Air swirling is achieved by 12 curved blades. The flow rate of

the center swriling air was kept the same as the outer dilution air.

4. Numerical Solutions

The mean velocity field is obtained from the solutions of Navier-Stokes equations

using the FAST-2D code. The code uses a finite volumn scheme with cell-centered
90 A. T. Norris and W.W.Liou

nodes. The pressure -velocity coupling is handled with the SIMPLEC algorithm.

The standard high-Reynolds number k - _ models were used to provide the turbulence mixing characteristics.

The scalar PDF code is similar to that developed by Hsu et al [2] but with a several

important numerical changes. These changes primarily involve the treatment of the

convection and molecular mixing. Details of these changes will be discussed in a

future paper.

5. Initial Conditions

The first flow measurement was reported at 0.5 cm downstream of the nozzle exit

plane. This was found unsatisfactory for use as initial conditions for a numerical

solution of the flow field. Therefore, the numerical initial velocity profiles for the

streamwise, radial, and azimuthal directions were determined by the air and fuel

mass flow rates, the angle of the swirl vanes, and the gaseus propane injection angle

at the nozzle exit. The initial turbulent kinetic energy was set at 1% of the mean

flow energy. The turbulent dissipation rate was then determined by assuming an

equavilent laminar viscosity.

Species initial conditions are simply pure fuel and air. To ignite the flame, we set

the initial species and temperature in the combustor to correspond to a fully burnt

fuel-air mixture at the experimental equivalence ratio.

6. Thermochemistry

Calculations were performed using the equilibrium chemistry assumption. In this

it is assumed that the chemical time scales are much smaller than the turbulent time

scales. In the future, a full mechanism, simplified by the Intrinsic Low-Dimensional

Manifold (ILDM) method of Pope and Maas [3], will be used.

7. Results

The calculated temperature and axial velocity profiles at four different streamwise locations in the combustor are compared with the measurement in Figs. 2 and

3, respectively. For comparison, the results obtained by a moment closure model

for the same combustor are also included. Fig. 2 shows the comparison of the axial

mean velocity profiles. At x--1 cm, the present calculation correctly predicts the

magnitude and the radial location of the peak velocity, which roughly corresponds

the edge of the air swirler. The reciculation zone indicated by the measurement

near the inner edge of the swirler is not predicted well by either the present calculation based on PDF method or the moment closure chemistry model. Therefore,

it is possible that this may be due to the k - e turbulence model, which is known

to underpredict flow recirculation in highly swirling flows. At the next three axial

stations, the PDF calculation does a better job capturing the correct local characteristics of the flow than the momentum closure procedure.

Fig. 3 show the radial temperature distributions at the same four axial locations as in Fig. 2. The PDF method has returned better prediction for the peak

termperature and profiles in the region near the nozzle exit, x=l and 4 cm than
Modeling of Turbulent, Reacting Flows by PDF Methods 91

the momentum closure chemistry model. The slightly higher predicted temperature may be due to either the insufficient turbulent mixing in the recirculation zone

returned by the k - E model or the equilibrium chemistry. The PDF method also

correctly predicted the centerline temperature. At farther downstream the PDF

model also show better temperature profile distribution. This can also be observed

from the relatively uniform temperature coutours in this region, Fig. 4. Note that

the experimenal flow shows higher turbulent mixing, resulting in more uniform temperature near the center region of the combustor. Again, this broad temperature

plateau region may be due to the deficiency of the turbulence models or the equilibrium chemistry. The preliminary results of finite rate chemistry calculations shows

a lower peak temperature at x---1 cm.

8. References

1 R. Charles, "Detailed Data Set: Velocity and temperature measurements in the

axisymmetric can combustor (ASCC) for a parametric variation in inlet conditions" UCI-ARTR-87-6, University of California, Irvine 1988.

2 A. T. Hsu, M. S. Raju and A. T. Norris, "Application of a PDF method to

compressible turbulent reacting flows" AIAA 94-0781, Reno NV 1994.

3 S.B.Pope and U.Maas, "Simplifying chemical Kinetics: Trajectory-generated lowdimensional manifolds" FDA-93-11, Cornell University, Ithaca NY 1993.
92 A. T. Norris and W.W.Liou

DILUTION

AIR /_ "/ J"

SWIRL AIR _

_'_UEL

_w,.LA,._ _ _

DILUTION _ _

AIR _ "_

_ CONE/ANNULAR NOZZLE

Figure 1. A sketch of the UC Irvine axisymmetric can combustor

3000.0

2000.0

1000.0

0.O

-1_.O

I j

0O

5000.0

4000.0

3000.0

2OOO.O

1000.0

0.0

(c)

U (cm/s)

4000.0

3000.0

2000.0

1000.O

0.O

-1000.0

0.O

i t • i

Co)

45OO.O

35OO.O

2500.0

1500.0

- • i ' i , i ,

_--- (d)

0 e • 0 $ Q • 0 I

1 , _ , I , 5_I.O

"° Fi_c)_°m"_° ".o °°

, I , I , I ,

1.O 2.0 3.O 4.0

Figure 2. Radial distribution of the mean streamwise velocity. •

Micklow et. a/[4]; , present.

, measurement[l]; - - -
Modeling of Turbulent, Reacting Flows by PDF Methods 93

2ooo.oI (a)

°°X \t

o.0, , . .... _ ,

2500,0

20(X)..O

1500.0

5OO.O

0.0

i i i

(c)

, , I L, 1 , I ,

3000.0

2500.0

2_00.0

T(°k)

1111111.0

O_

i i I

(b)

, I , I - I ,

o.oRadi'_l Locat,on, 2.0. _(cm)"°

:3OOO.O

2500.0

2000.0

,5_0.0

IOO0.O

500.0

0.0

0.0

i i ' i

(d)

-_ e e

i _ i 1 , i ,

1.o 2.o 3.o

Figure 3. Comparison of Radial distribution of temperature.

Micklow et. a/[4];--, present.

4.0

• , measurement[I]; - - -

x=l cm 4 cm

I I

10 cm

_-----3

24 cm

Level T

8 1700

__ PDF

6

._ .7 ,.

7 1500

6 1300

5 1100

4 9O0

3 700

2 5O0

1 300

Figure 4. Temperature contours
Page intentionally left blank 
Center for Modeling of Turbulence and Transition

Research Briefs - 1995

95

Assessment of Turbulence

Models in Turbomachinery Flows

1. Motivation and Objectives

Turbomachinery flows pose a major challenge for CFD since they are threedimensional and time dependent. NASA Lewis Research Center is playing a lead

role in developing CFD tools which can be used by industry in the analysis of their

designs.

Recently, a blind test case for a rotor compressor (ROTOR 37) was organized by

the ASME/International Gas Turbine Conference held at the Hague (A. Strazisar,

J. Wood, and K. Suderl). The objective was to judge the predictive capabilities

of turbomachinery CFD tools. A total of 12 groups submitted their predictions

for this case, which were then compared with the detailed experimental data. The

results, which are given in the above reference, showed a wide scatter relative to

the experimental data. There were several reasons for this, one of them being the

inadequacies of the turbulence models.

Since different codes using the same turbulence models produced different results

in the ROTOR 37 blind test case, it is obvious that there is a need to assess

the performance of different turbulence models from the same numerical platform.

If different turbulence models are implemented in a single CFD code, then the

differences in the results can be attributed only to the differences in the models.

In this way the impact of different models on the prediction of turbomachinery

flows can be assessed in a systematic manner. In order to answer this question,

CMOTT has started a joint project with the turbomachinery group at the Lewis

Research Academy. Since the Average Passage Equation code, VSTAGE, predicted

the ROTOR 37 blind test case reasonably weU 2, it was selected as the numerical

platform to carry out this task.

The approach taken is to implement several turbulence models, of varying sophistication, in the VSTAGE code in the form of a subroutine. To begin with the

standard k - e model of Launder and Spaulding 3 will be implemented. This will be

followed by implementing an improvement to this model which has been developed

by CMOTT. This improvement has already shown to give better results over the

standard k - e model for a some benchmark flows. Eventually a complete CMOTT

two equation model (Shih, et. al.4) will be implemented.

All of the above models will be used with the wall functions. The standard wall

functions approach in CFD calculations assumes the existence of the law of the

wall, even though this is true for a flat plate boundary layer only. This approach

is popular because it requires a mesh size that is computationatly affordable and

because it, apparently, provides good engineering accuracy. At a later time, and

based on the results obtained using the above models, there may be a need to assess

PRr--GLDir,_G PAGt,_ _JLAr!, [,;OT FiLMEr)
96

the performance of the different wall functions on the prediction of turbomachinery

flows. It should be pointed out that all of the turbulence model assessment will be

done against the ROTOR 37 data.

In the course of this project we may find some deficiencies in the turbulence

models. It is also the objective of this effort to identify and remedy these.

This project will make the following two important contributions. First an assessment of the hierarchy of turbulence models (including those that CMOTT has

developed) in the context of turbomachinery flows. This information will be extremely useful to the turbomachinery community. Second the models which perform the best against the test cases can be incorporated in the Average Passage

Equation codes for both the single stage and the multistage turbomachinery. This

will lead to an improved turbomachinery CFD tool.

2. Work in Progress

We have finished implementing the standard two equation k - E model and its

CMOTT improvement in the VSTAGE code. The preliminary results for the ROTOR 37 test case look promising. We are currently conducting a grid sensitivity

study. The results of this work will be presented in an upcoming report. This work

is being carried out by A. Shabbir, J. Zhu, M. Celestina, J.J. Adamczyk, and T.-H.

Shih.

3. References

I Strazisar,A., Wood, J.,and Suder, K. ASME/International Gas Turbine Conference, The Hague, Netherlands, 1994.

2 Celestina, M. ASME/International Gas Turbine Conference, The Hague, Netherlands, 1994.

3 Launder, B.E., and Spaulding, D.B., "The numerical computation of turbulent

flows", Comp. Meth. Appl. Mech. Engg., 3, 269-289, 1974.

4 Shih, T.-H., W.W. Liou, A. Shabbir, Z. Yang and J. Zhu. "A new high Reynolds

number k - e model", Computers Fluids, 24, 3,227-238, 1995.
Center for Modeling of Turbulence and Transition

Research Briefs - 1995

97

Appendix A

Industry-Wide Workshop on

Computational Turbulence Modeling

CMOTT organized an Industry-wide Workshop on Computational Turbulence Modeling on October 6-7, 1994. The purpose of the workshop was to initiate the transfer of

technology developed at Lewis Research Center (LeRC) to industry and to discuss the

current status and the future needs of turbulence models in industrial CFD. The workshop organizing committee and the titles of presentations are listed below. Details of the

presentations can be found in the NASA CP 10165.

Workshop Organizing Committee

Industries

C. Prakash, General Electric

M. Sindir, Rocketdyne

S. Syed, Pratt & Whitney

Universities

LY. Chen, University of California, Berkeley

LL. Lumley, Comell University

National Aeronautics and Space Administration

L.A. PovineUi, Chairman

IL Mankbadi, Lewis Research Center

D.1L Reddy, Lewis Research Center

P. Richardson, Headquarters

RJ. Shaw, Lewis Research Center

Institute for Computational Mechanics in Propulsion

T. Keith, Ohio Aerospace Institute

A. Shabbir, Center for Modeling of Turbulence and Transition

T.-H. Shih, Center for Modeling of Turbulence and Transition
98 Center ]or Modeling of Turbulence and Transition

Research Briefs - 1995

Workshop Presentations

TURBULENCE PROGRAM FOR PROPULSION SYSTEMS

Tsan-Hsing Shih, Institute for Computational Mechanics in Propulsion and Center for +

Modeling of Turbulence and Transition, NASA Lewis Research Center

TURBULENCE MODEL DEVELOPMENT AND APPLICATION AT LOCKHEED FORT

WORTH COMPANY

Brian R. Smith, CTD Group, Lockhev_i Fort Worth Company

A SUMMARY OF COMPUTATIONAL EXPERIENCE AT GE AIRCRAFt ENGINES

FOR COMPLEX TURBULENT FLOWS IN GAS TURBINES

R. Zerklc and C Prakash, GE Aircraft Engines

THE APPLICABILITY OF TURBULENCE MODELS TO AERODYNAMIC AND PROPULSION

FLOWFIELDS AT McDONNELL DOUGLAS AEROSPACE

Linda D. Kral, John A. Ladd, and Mori Mani, McDonnell Douglas Aerospace

EXPERIENCE WITH k-e TURBULENCE MODELS FOR HEAT TRANSFER

COMPUTATIONS IN ROTATING

Prabhat Tckriwal, GE Corporate Research and Development

TURBULENCE MODELS FOR GAS TURBINE COMBUSTORS

AndrejaBrankovic,CFD Group,Pratt& Whitney

COMBUSTION SYSTEM CFD MODELING AT GE AIRCRAFr ENGINES

D. Burrus and H. Mongia, GE A/rcraft Engines, and A. Tolpadi, S. Corma, and M. Braaten,

GE Corporate Research and Development

RECENT PROGRESS IN THE JOINT VELOCITY-SCALAR PDF METHOD

M.S. Anand, Allison Engine Company

OVERVIEW OF TURBULENCE MODEL DEVELOPMENT AND APPLICATIONS AT

ROCKETDYNE

A.H. Hadid, E.D. Lynch, and/VIM. Sindir, Roekctdyne Division, Roek'well

tntamational

RECENT ADVANCES IN PDF MODELING OF TURBULENT REACTING FLOWS

A.D. Leonard and F. Dai, CFD Research Corporation

EXPERIENCE WITH TURBLrLENCE INTERACTION AND TURBULENCE-CHEMISTRY

MODELS AT FLUENT INC.

D. Choudhury, S.E. Kirn, D.P. Tselepidakis, and M. Missaghi, Fluent Inc.

EXPERIENCES WITH TWO-EQUATION TURBULENCE MODELS

Ashok IC Singhal, Yong G. Lai, and Ram K. Awa, CFD Research Corporation

PROGRESS IN SIMULATING INDUSTRIAL FLOWS USING TWO-EQUATION MODELS:

CAN MORE BE ACHIEVED WITH FURTHER RESEARCH?

Vahd Haroutunian, Fluid Dynamics International, I.nc
Center for Modeling of Turbulence and Transition

Research Briefs - 1995

99

TURBULENCE MODELING NEEDS OF COMMERCIAL CP-'D CODES: COMPLEX FLOWS

IN THE AEROSPACE AND AUTOMOTIVE INDUSTRIES

Bizhan A. Befrui, adapco

TURBULENCE REQI.II:REMENTS OF A COMMERCIAL CFD CODE

J.P. Van Doormaal, C.M. Mueller, and M.J. Raw, Advanced Scientific

Computing Ltd.

SECOND-ORDER CLOSURES FOR COIVIPRF_SSIBLE TURBULENCE

J.L. Immley, S. Savarese, and C.C. Volte, Mechanical and Aerospace Engineering

Department, ComeU University

MODELING OF TURBULF2qT CHEMICAL REACTION

L-Y. Chen, Department of Mechanical Engineering, University of California,

Berkeley

INTRODUCTION TO TURBULENCE SUBPROGRAM

T.-H. Shih and J. Zhu, Institute for Computational Mechanics in Propulsion and Center for

Modeling of Turbulence and Transition, Lewis Research Center

DESCRIFHON OF TURBULENCE SUB-PROGRAM

J. Zhu, Institute for Computational Mechanics in Propulsion, NASA Lewis Rese.arch Center .

OVERVIEW OF PROBABILITY. DENSITY FUNCTION (PDF) MODELING AT LeRC

D_ Reddy, Internal Fluid Mechanics Division, NASA Lewis Research Center

PDF METHODS FOR TURBULENT REACTIVE FLOWS

Andrew T. Hsu, NYMA, Inc., NASA Lewis Research Center

A COMPOSITION JOINT PDF METHOD FOR THE MODELING OF SPRAY FLAMES

M.S. Raju, Nyma, Inc., NASA Lewis Research Center

IMPROVEMENTS AND NEW FEATURES IN THE PDF MODULE

AT. Norris, Institute for Computational Mechanics in Propulsion
Page intentionally left blank 
101

Appendix B

Organization- 1995

Position Chart
102

Center for Modeling of Turbulence and Transition

Research Briefs - 1995

ICOMP Director Dr. Louis A. PovineUi

Deputy Chief

Internal Fluid Mechanics Division

NASA Lewis Research Center

CMOTT Technical Leader Dr. Tsau-Hsing Shih

Senior Research Associate

ICOMP, NASA Lewis Research Center

Advisors Dr. Marvin E. Goldstein

Chief Scientist

NASA Lewis Research Center

Professor John L. Lumley

Sibley School of Mechanical and

Aerospace Engineering

CorneU University

Professor Eli Reshotko

Department of Mechanical and

Aerospace Engineering

Case Western Reserve University
Center for Modeling of Turbulence and Transition

Research Briefs - 1995 103

Listing of current

Names/Term

Liou, William W.

11/1990 - present

Norris, Andrew

4/1993 - present

Shabbir, Aamir

5/1990 - present

Shih, Tsan-Hsing

5/1990 - present

Yang, Zhigang

7/1990- present

Zhu, Jiang

4/1992 - present

Members

Affiliation

ICOMP

ICOMP

ICOMP

ICOMP

ICOMP

ICOMP

Research Areas

Compressible Flow

Modeling, Weakly

Nonlinear Wave Models

PDF Turbulence Modeling,

DNS

Buoyancy Effects on

Turbulence, Turbulence

Modeling

Turbulence Modeling

Modeling of Bypass

Transition

Application of Turbulence

Models in Complex Flows
Page intentionally left blank 
Center for Modeling of Turbulence and Transition

Research Briefs - 1995

105

Appendix C

List of Member's Publications

.Barton, J. M., Rubinstein, R. and Kirtley, K. R. "Nonlinear Reynolds stress model for

turbulent shear flows, " AIAA Paper No. 91-0609, (Jan. 1991).

Barton, J. M. and Rubinstein, R. _, "Nonlinerar algebraic Reynolds stress model for

anisotropic turbulent flows, " llth U. S. National Congress of Applied Mechanics, Tucson,

AZ (May 1990).

Barton, J. M. and Rubinstein, R., "Renormalization group theory and turbulence modeling, " 2nd International Workshop on Chaos and Turbulence, Tsukuba, Japan (Jan.

1992).

Barton, J. M. and Rubinstein, R., "Renormalization group analysis of turbulence- driven

secondary flows, " llth Australian Fluid Mechanics Conference, Hobart, Australia (Dec.

1992).

Brown, S., Leibovich, S., and Yang, Z., "On the linear instability of the Hall-Stewartson

vortex," Theo. and Comp. Fluid Dynamics, Vol 2, no. 1, (1990).

Duncan, B. S., Liou, W. W. and Shih, T.-H., "A multiple-scale turbulence model for

incompressible flow," AIAA 93-0086 (1993).

Duncan, B.S., Lumley, J.L., Shih, T.H. and To, W.M., "A new model for the turbulent

dissipation" International Conference of Fluid Mechanics and Theoretical Physics, Bejing,

China (1992).

Durbin, P.A. and Ya_g, Z., "A transport equation for eddy viscosity," Proc. of the 1992

Summer Program, CTR, Stanford, 1992.

Georgiadis, N.J., Chitsomboon, T. and Zhu, J., "Modification of the two-equation turbulence model in NPARC to a Chien low Reynolds number k-e formulation," NASA TM

106710 (1994).

Hsu, A.T., Raju, M.S., and Norris, A.T., Application of a PDF method to compressible

turbulent reacting flows. In AIAA 32nd Aerospace Sciences Meeting and Exhibit, Reno,

Nevada, pages AIAA-94-0781, 1994.

Hsu, A.T. and Liou, M., "Computational analysis of underexpanded jets in the hypersonic

regime," Journal of Propulsion and Power, 7, No. 2, (1991).

Hsu, A.T., "Progress in the development of PDF turbulence models for combustion," 10th

NASP Symposium, April 23-16, 1991, Monterey, California.

Hsu, A.T., "A study of hydrogen diffusion flames using PDF turbulence model," AIAA

91-1780, AIAA 22nd Fluid Dynamics Conference, June 24-26, 1991, Honolulu, Hawaii.

,/
106 Publications List

Hsu, A.T. and Chen, J.Y., "A continuous mixing Model for PDF simulations and its

applications to combusting shear flows," 8th International Symposium on Turbulence Shear

Flows, Sept. 9-11, 1991, Munich, Germany.

Hsu, A.T., Raju, M.S., and Norris, A.T., Application of a PDF method to compressible

turbulent reacting flows. In AIAA 32nd Aerospace Sciences Meeting and Exhibit, Reno,

Nevada, pages AIAA-94-0781, 1994. Huang, P.G. and Liou, W.W., "Numerical calculations of shock-wave/boundary-layer flow interactions," NASA TM 106694 (1994).

Lang, N.J. and Shih, T.-H., "A critical comparison of two-equation turbulence models,"

NASA TM 105237, (1991).

Liou, M. S. and Steffen, Jr., C. J., "A new flux splitting scheme," NASA TM 104404

(1991).

Liou, W. W. and Shih, T.-H., "Transonic flow predictions with new two--equation turbulence models," AIAA 95-1805.

Liou, W. W. and Shih, T.-H., "Compressible turbulent boundary layer predictions using

a multiple-scale model," NASA TM (in preparation) (1995).

Liou, W. W., Shih, T.-H., and Duncan, B.D., "A multiple-scale model for compressible

turbulent flows," Phys. Fluids, 7 (3), 658-666 (1995).

Liou, W. W., "Linear instability of curved free shear layer," Phys. F!uids, 6 (2), 541-549

(1994).

Liou, W. W. and Shih, T.-H., "On the basic equations for sencond-order modeling of

compressible turbulence," NASA TM 105277 (1991).

Liou, W. W. and Morris, P. J., "Weakly nonlinear models for turbulentmixing in a plane

mixing layer," Phy. Fluids, 4, 2788-2808 (1992).

Liou, W. W. and Morris, P. J., "The eigenvalue specturm of the P_yleigh equation for a

plane shear layer," International Journal of Numerical Methods in Fluids, 15, 1407-1415

(1992).

Liou, W. W., "A new energy transfer model for turbulent free shear flow," NASA TM

105854 (1992).

Mansour, N N. and Shill, T.-H. and Reynolds, W. C., "The effects of rotation on initially

anisotropic homogeneous flows," Physics of fluid A, 3, 10, 2421-2425 (1991).

Michelassi, V. and Shih, T.-H., "LowReynolds number two-equation modeling of turbulent

flows," NASA TM 104368, (1991).

Michelassi, V. and Shih, T.-H., "Elliptic flow computation by low Reynolds number twoequation turbulence models ," NASA TM 105376, (1991).

Moin, P., Shill, T.-H., Driver, D. and Mansour, N., "Direct numerical simulation of a

three-dimensional turbulent boundary layer," Physics of Fluids A 2, 10, 1846-1853 (1990).
Publications List 107

Norris, A.T. and Pope, S.B., Modehng of extinction in turbulent diffusion flames by the

velocity-dissipation-composition PDF method. Combustion and Flame., 100(1/2):211-

220, 1994.

Norris, A.T. and Hsu, A.,T., Comparison of PDF and moment closure methods in the

modeling of turbulent reacting flows. In 30th AIAA/ASME/SAE/ASEE Joint Propulsion

Conference, Indianapolis, Indiana, pages AIAA-94-3356, 1994.

Norris, A.T., "The application of PDF methods to piloted diffusion flames," Ph.D. Thesis,

Cornell University, 1993.

Rubinstein, R. and Barton, J. M., "Renormalization group analysis of the Reynolds stress

transport equations, nonlinear Reynolds stress models and the renormalization group,"

Phys. Fluids A 4, 1759 1766 (1992).

Rubinstein, R. and Barton, J. M., "Nonlinear Reynolds stress models and the renormalization group," Phys. Fluids A 2, t472 (1990).

Rubinstein, R. and Barton, J. M., "Renormalization group analysis of anisotropic diffusion

in turbulent shear flow," Phys. Fluids A 3, 415 (1991).

Rubinstein, R. and Barton, J. M., "Renormalization group analysis of the Reynolds stress

transport equation," Phys. of Fluids A 4, 1579 (1992).

Shabbir, A., and George, W.K. "Experiments on round turbulent buoyant plume", J. Fluid

Mech, 1994, 275, 1-32.

Shabbir, A., Beuther, P.D. and George, W.K. "X-Wire response in turbulent flows of high

intensity turbulence and low mean velocities", 1995, Accepted to Experimental Fluid &

Thermal Science.

Shabbir, A., "Experimental balances for the second moments for a buoyant plume and

their implication on turbulence modeling," Eighth Symposium on Turbulent Shear Flows,

27-1-1 to 27-1-6 (1991).

Shabbir, A. and Shih, T.-H., "Critical comparison of second order turbulence models in

homogeneous flows," Presented at the AIAA Annual meeting in Reno, Jan. 1992.

Shabbir, A. and Taulbee, D.B., Evaluation of Turbulence Models for Predicting Buoyant

Flows. J. Heat Transfer, 1990, Vol 112, No 4, pp 945-951.

Shabbir, A. and Taullbee, D.B., Evaluation of the Coefficients of the Reynolds Stress

Models Using Experiments on Homogeneous Turbulence. A. Shabbir and D.B. Taulbee.

Forum on Turbulent Flows. ASME (1987), FED, 51, pp 67-72.

Shih, T.-H., Zhu, :I. and Lumley, J.L., "A new Reynolds stress algebraic equation

model,'NASA TM 106644 (1994).

Shih, T.-H., Zhu, J. and Lumley, J.L., "Modeling of wall-bounded complex flows and free

shear flows," NASA TM 106513 (1994).
108 Publications List

Shih, T.H., Liou, W.W., Shabbir, A., Yang, Z., and Zhu, J., "A new k - e eddy viscosity

model for high Reynolds number turbulent flows," J. Computers and Fluids, 24, (1995).

Shih, T.-H. and Shabbir, A., "Methods of ensuring realizability for non-realizable second

order closures," NASA TM 106681 (1994).

Shih, T.-H., Shabbir, A. and Lumley, J.L., "Realizability in second moment turbulence

closures revisted,", NASA TM 106469 (1994).

Shih, T.-H., Liou, W.W., Shabbir, A., Yang, Z., and Zhu, J., "A vorticity dynamics based

model for the turbulent dissipation," NASA TM 106177 (1993).

Shih, T.-H. and Lumley, J.L., "Remarks on turbulent constitutive relations," NASA TM

106116 (1993).

Shih, T.-H. and Lumley, J.L., "Kolmogorov behavior of near-wall turbulence and its application in turbulence modeling," NASA TM 105663 (1992)

Shih, T.-H., Chen, J.-Y. and Lumley, J. L., "Second order modeling of free turbulent shear

flows," AIAA Journal, 30, 6, 1553-1560 (1992).

Shih, T.-H. and Lumley, J.L., "A critical comparison of second order closures with direct

numerical simulation of homogeneous turbulence," NASA TM 105351 (1991).

Shih, T.-H. and Hsu, A.-T., "An improved k - e model for near-wall turbulence," AIAA

paper, 91-0611 (1991).

Shih, T.-H. Shili, Shabbir, A. and Lumley, J.L., "Advances in modeling the pressure correlation terms in the second moment equations," The Lumley Symposium: Resent DeveIopments in Turbulence, November, 12-13, 1990, ICASE, NASA Langley Research Center,

Edited by Gatski, T.B., Sarkar, S., Speziale, C.G.

Shih, T.-H., "An improved k-e model for near-wall turbulence and comparison with direct

numerical simulation," NASA TM 103221 (1990).

Shih, T.-H. and Mansour, N. N., "Modeling of near-wall turbulence," Engineering Turbulence Modeling and Experiments, September 24-28, 1990, Dubrovnik, Yugoslavia, Editors:

W.Rodi, E.N. Ganic.

Shih, T.-H. and Reynolds, W. C., "A spectrum model for weakly anisotropic turbulence,"

Physics of Fluids A, 2, 8, (1990).

Steffen, Jr., C. J., "A critical comparison of several low Reynolds number k-epsilon turbulence models for flow over a backward-facing step," NASA TM 106173 (1993).

Steffen, Jr., C. J., "An investigation of DTNS2D for use as an incompressible turbulence

modelling test-bed," NASA TM 105593 (1992).

Steffen, Jr., C. J. and Beard, L. M., "Incompressible Navier Stokes solutions using the

pseudo compressibility technique," Proceeding of the 23rd Annual Pittsburgh Conference

on Modeling and Simulation, in print (1992).
Publications List 109

Van der Vegt, Jaap, "ENO-Osher schems for Euler equations," AIAA Paper 93-0335 (1993)

Van der Vegt, J., "Assessment of flux splitting for viscous compressible flows," AIAA Paper

91-0243 (1991).

Van der Vegt, Jaap, "Overview of the Osher approximate Riemann solver for threedimensional flows," NASA TM, in print (1992).

Yang, Z. and Leibovich, S., "On the lindar marginal stability of the trailing line vortex,"

Submitted for publication.

Yang, Z. and Leibovich, S, "Nonlinear dynamics near the stability margin in rotating pipe

flow," J. Fluid Mech. 233, 329-347, (1991).

Yang, Z. and Leibovich, S., "Unstable viscous wail modes in rotating pipe flow," AIAA

Paper 91-1801, (1991). Also submitted to Physics Fluids A for publication.

Yang, Z. and Shih, T.H., "A k - e modeling of near wall turbulence," Proc. of 4th Intl.

Sym. on CFD, Davis, CA, 1305-1310, (1991). Also available as NASA TM 105238.

Yang, Z. and Shih, T.H., "A k - e calculation of transitional boundary layers," Instability,

IYansition, and Turbulence, M.Y. Hussalni, A. Kumar and C.L. Steett, ed. SpringerVerlag, 1992.

Yang, Z. and Leibovich, S., "On the lindar marginal stability of the trailing line vortex," Presented at the 18th International Congress of Theoretical and Applied Mechanics,

September 1992, Haifa, Israel. Submitted for publication.

Yang, Z. and Shih, T.H. "A new time scale based k - e model for near wall turbulence,"

hIhi J. 31, 1191 (1993).

Yang, Z. and Shih, T.H., "A k - e model for turbulent and transitional boundary layers,"

Proc. Near Wall Turbulent Flows, B.E. Launder, C. Speziale and R.M.C. So, ed. Elsevier,

1993.

Yang, Z. and Shih, T.H., "A Galilean and tnesorial invariant k - e model for near wall

turbulence," AIAA Paper 93-3105 (1993).

Yang, Z. and Shih, T.H., "Flow inhomogeneity and the dissipation rate equation." Accepted for presentation at the 5th European Turbulence Conference, July 1994, Siena,

Italy. In preparation for publication.

Yang, Z., Georgiadis, N., Zhu, J., and Shih, T. H., "Calculations of inlet/nozzle flows using

a new k-e model," to be presented at the AIAA Joint Propulsion Conference, San Diego,

July (1995).

Zhu, J. and Shih, T.H., "Calculations of turbulent separated flows," NASA TM 106154

(1993). also in CFD Journal 3, pp.343-354 (1994).

Zhu, J. and Shih, T.H., "A numerical study of confined turbulent jets," NASA TM 106197

(1993), also in J. Fluids Eng. 116, pp.702-706 (1994).
110 Publications List

Zhu, J. and Shih, T.H., "Computation of confined coflow jets with three turbulence models," AIAA Paper 93-3120, also in Int. J. for Num. Methods in Fluids 19, pp.939-956

(1994).

Zhu, G., Lai, M.-C and Shill, T.-H., "Second-order closure modeling of turbulent buoyant

wall plumes," NASA TM (in print) (1992).

Zhu, J., "FAST-2D: A computer program for numerical simulation of two dimensional incompressible flows with complex boundaries," Rept. No.690, Institute for Hydromechanics,

University of Karlsruhe (1991).
Form Approved

REPORT DOCUMENTATION PAGE OMBNo. 0704-0188

Publicreportingburd_lforthis ¢o41ecZion of informationisestimatedtoaverage1 hourperresponse,includingthe timeforreviewinginstructions,seamhingexi=mgdatasourcce_..,

gatheringandmaJntaittingthedaXaneeded,andco_ing andreviewingthecollectionof _rdotma=ion..S_ _¢_nts regardingthisburdenl_t.imaleer =my=her alkoectofthes

collectionof informality,includingsuggestionsforreducingthisburden,1oWashingtonHeadquartersSerwce$,Dfectorateler InformationOpera_,onsand Reports,1215Jefferson

DavisH_bway,Suite1204,/_li_gton,VA 2220_-4302,andtothe Officeof ManagementandBudget,PaperworkReductionProjec=(0704-0188),Weshk_gton. DE;20503.

1. AGENCY USE ONLY (Leave blank)

4. TITLE AND SUBTITLE

2. REPORT DATE 3. REPORT TYPE AND DATES COVERED

October 1995 ConlzacmrReport

5. FUNDING NUMBERS

Center for Modeling of Turbulence and Transition: Research Briefs---1995

s. Atn'HOR(S)

Center for Modeling of Turbulence and Transition

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

Institute for Computational Mechanics in Propulsion

22800 Cedar Point Road

Cleveland, Ohio 44142

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

National Aeronautics and Space Administration

Lewis Research Center

Cleveland, Ohio 44135-3191

WU-505-90-5K

NCC3-370

8. PERFORMING ORGANIZATION

REPORT NUMBER

E-9932

10. SPONSORING/MONITORING

AGENCY REPORT NUMBER

NASA CR-198408

ICOMP-95-17

CMOTT-95-3

11. SUPPLEMENTARY NOTES

ICOMP Program Director, Louis A. Povinelli, organization code 2600, (216)433-5818.

12a. DISTRIBUTION/AVAILABILITY STATEMENT

Unclassified -Unlimited

Subject Category 34

This publication is available from the NASA Center for Aerospace Information, (301) 621-0390. !

12b. DISTRIBUTION CODE

13. ABSTRACT (Maximum 200 words)

This research brief contains the progress reports of the research staff of the Center for Modeling of Turbulence and

Transition (CMOTT) from July 1993 to July 1995. It also constitutes a progress report to the Institute of Computational

Mechanics in Propulsion located at the Ohio Aerospace Institute and the Lewis Research Center. CMOTT has been in

existence for about four years. In the fast three years, its main activities were to develop and validate turbulence and

combustion models for propulsion systems, in an effort to remove the deficiencies of existing models. Three workshops

on computational turbulence modeling were held at LeRC (1991, 1993, 1994). At present, CMOTT is integrating the

CMOTT developed/improved models into CFD tools which can be used by the propulsion systems community. This

activity has resulted in an increased collaboration with the Lewis CFD researchers.

14. SUBJECT TERMS

Turbulence modeling

17. SECURITY CLASSIFICATION

OF REPORT

Unclassified

18. SECURITY CLASSIFICATION

OF THIS PAGE

Unclassified

NSN 7540-01-280-5500

19. SECURITY CLASSIFICATION

OF ABSTRACT

Unclassified

15. NUMBER OF PAGES

114

16. PRICE CODE

A06

20. LIMITATION OF ABSTRACT

Standard Form 298 (Rev. 2-89)

Pre_d by ANSI Sld. Z39-18

298--102