NASA Technical Memorandum 105243

ICOMP-91-17; CMOTT-91-07

Center for Modeling of Turbulence

and Transition (CMOTT)

Research Briefs 1990

October 1991

NASA

LEWIS RESEARCH CENTER

iCOMP

CASE WESTERN

RESERVE UNIVM[TY

NASA
Page intentionally left blank
14
Page intentionally left blank
Center for Modeling of Turbulence and Transition

Research Briefs - 1990

Table of Contents

Page

Preface . . . . . . . . . . . . . . . . . . . . . . . . . . . vii

The Study of PDF Turbulence Models in Combustion

Andrew T. Hsu . . . . . . . . . . . . . . . . . . . . . . . . 1

Turbulence Modeling

Tsan-Hsing Shih . . . . . . . . . . . . . . . . . . . . . . . . 4

Experiments and Modeling

Aamir Shabbir . . . . . . . . . . . . . . . . . . . . . . . . 9

Modeling of Compressible Turbulent Shear Flows

William "I . Liou . . . . . . . . . . . . . . . . . . . . . . 13

RNG in Turbulence and Modeling of Bypass Transition

Zhigang Yang . . . . . . . . . . . . . . . . . . . . . . . . 16

Turbulence Modeling in Supersonic Combusting Flows

Tativit Chitsomboon . . . . . . . . . . . . . . . . . . . . . 19

Development of New Flux Splitting Schemes

Meng-Sing Liou and Christopher J. Steffen, Jr. . . . . . . . . . 22

Progress of Simulations for Reacting Shear Layers

Sheng-Tao Yu . . . . . . . . . . . . . . . . . . . . . . . . 24

Appendix A : Organization - 1990 . . . . . . . . . . . . . 26

Appendix B : Seminars and Meetings . . . . . . . . . . . . 29

B.1 CMOTT Seminar Series . . . . . . . . . . . . . . 29

B.2 CMOTT Technical Meetings . . . . . . . . . . . . 48

Appendix C : Collection of Publications . . . . . . . . . . 52

V
Page intentionally left blank
Center for Mod,!ing of Turbulence and Transition

Research Briefs - 1990

Preface

The Center for Modeling of Turbulence and Transition (CMOTT), a cooperative turbulence research team, was formally established in May of 1990

as a focal group within the Institute of Computational Mechanics for Propulsion (ICOMP). The location of CMOTT is shown in the organization structure chart in Appendix A. The objectives of the CMOTT are to develop, validate and implement the models for turbulence and boundary-layer transition

in practical engineering flows. The flows of interest are three-dimensional,

incompressible and compressible flows with chemistry. The schemes being

studied include the two-equation (e.g. k-E) and algebraic Reynolds-stress

models, the full Reynolds-stress (or second moment closure) models, the

probability density function (pdf) models, the Renormalization Group Theory (RNG) and Direct Interaction Approximation (DIA), the Large Eddy

Simulation (LES) and Direct Numerical Simulation (DNS).

Currently, CMOTT has eight formal members working on various aspects of turbulence and transition modeling in collaboration with NASALewis scientists and Case Western Reserve University (CWRU) faculty members. The CMOTT members have been actively involved in international and

national turbulence research activities through meetings, seminars, workshops and exchange-visitors. Since June of 1990, a CMOTT seminar series

has been conducted with speakers invited from within and outside of the

NASA Lewis Research Center, including foreign speakers. In 1991, a new

series of biweekly CMOTT technical meetings was initiated for informal discussions regarding special issues in turbulence and transition modeling. The

CMOTT research activity is advised by a group consisting of Professor J.L.

Lumley (Cornell University), Dr. M. Goldstein (NASA/LeRC) and Professor

E. Reshotko (CWRU ).

This research brief contains the progress reports of the CMOTT Research staff from May 1990 to May 1991. It is intended to be an informational

report of the CMOTT activities as well as an annual report to ICOMP and

NASA. The current CMOTT roster and its organization are listed in the

Appendix A. Listed in Appendix 13.1 are the visiting members and their

seminar abstracts. Appendix 13.2 gives the scientific and technical issues discussed in biweekly CMOTT meetings. Journal and conference publications

by CMOTT members are grouped in Appendix C.

V11
Center for Modeling of Turbulence and 7'rangitin::

Research. Briefs - 1990

Starting in 1991, NASA Technical Memoranda authored by members of

the CMOTT staff will be given a specific number to identify them as CMOTT

reports. These manuscripts will be made available for early dissemination of

completed research results by the CMOTT staff.

Finally, we express our thanks to one of the CMOTT members, Dr.

William W. Liou, who carefully assembled the material and provided editorial assistance.

Louis Povinelli

Meng-Sing Liou

Tsan-Hsing Shih

viii
Center for Modeling of Turbulence and 7'ransition

Research Briefs - 1990

The Study of PDF Turbulence Models in Combustion

Andrew T. Hsu

1. Motivation and Objectives

1.1 Motivation

The accurate prediction of turbulent combustion is still beyond reach

for today's computation techniques. It is the consensus of the combustion

profession that the predictions of chemically reacting flow were poor if conventional turbulence models were used. The main difficulty lies in the fact

that the reaction rate is highly non-linear, and the use of averaged temperature, pressure and density produces excessively large errors. The probability.density function (pdf) method is the only alternative at present time

that uses local instant values of the temperature, density, etc., in predicting

chemical reaction rate, and thus is the only viable approach for turbulent

combustion calculations.

1.2 Objectives

The present work aims at the development and implementation of the

the pdf turbulence models in solving realistic combustion problems. The fact

that the pdf equation has a very large dimensionality renders finite difference

schemes extremely demanding on computer memories and thus impractical,

if not entirely impossible. A logical alternative is the Monte Carlo scheme,

which has been used extensively in statistical physics. However, the use of

Monte Carlo scheme to solve both the flowfield and the chemical reaction is

very time consuming. Further more, since CFD has reached a certain degree

of maturity as well as popularity, it seems less beneficial to abandon CFD

completely and opt for Monte Carlo schemes. Therefore, we propose the

use of a combined CFD and Monte Carlo scheme in the present study. The

scheme would use the conventional flow solvers when calculating the flowfield

properties such as velocity, pressure, etc., while the chemical reaction part

would be solved using Monte Carlo solvers.
2. Works Accomplished

2.1 Code development.

A parabolic code with k - E turbulence models have been developed

in the past months. Three different k - E models have been tested with

satisfactory numerical resuits.

A grid dependent Monte Carlo scheme is being explored. This scheme

discretize the pdf equation on a given grid and write, for parabolic flows:

P2+dT>j = aiP. j+i + Oj Py ,j + -(jp.,j- 1

(1)

and we require aj+Qj+^6 =1 (2)

Using a very simple test case of a convection / diffusion process with two

scalars, it was found that the previous scheme does not conserve mass fractions due to re-contamination. It is found that in order to conserve the

mass fractions absolutely, one needs to add further restriction to the scheme,

namely

aj + ^^j = aj-1 + 7.7+ 1

(3)

A new algorithm was devised and is currently being tested. Again using

the simple test case of two scalars with assumed constant coefficients in the

pdf equation, the new algorithm is shown to conserve the mass fractions

perfectly in cases of uniform flows or pure diffusion problems. Deficiencies

such as directional bias and re - contamination that were found in the previous

algorithm are completely eliminated.

2.2 Applications.

The code developed has been validated by solving a heated turbulent

jet. The temperature is treated as a conserved passive scalar and solved

using the pdf Monte Carlo simulation while the flow field is obtained using

a conventional CFD solver. The mean temperature profile and RMS of the

temperature fluctuation were compared with experimental data.

As a first application to combustion problems, the non-premixed flame

of hydrogen and fluorine is being studied. A comparison between primary

results from the present study and experimental data show that the present

scheme predicts the mean flame temperature accurately.

2
3. Future Plans

1. Further investigate the case of hydrogen-fluorine reaction.

2. Study finite rate calculation of the same non-premixed flame.

3. Study the interaction between mixing and chemical reaction.

4. Study compressibility effects.

4. Publications

1. Hsu, A.T., "The Study of PDF Turbulence Models in Combustion,"

9th National Aero-Space Plane Technology Symposium, November 1-2,

1990.

2. Hsu, A.T., " On Recontamination and Directional Bias Problem in

Monte Carlo Simulation of PDF Turbulence Models," NASA CFD Conference, April 12-14, 1991, Moffett Field, California.

3. Hsu, A.T., "Progress in the Development of PDF Turbulence Models for

Combustion," 10th National Aero-Space Plane Technology Symposium,

April 23-25, 1991, Monterey, California.

4. Hsu, A.T., "The Study of PDF Turbulence Model in Nonequilibrium

Hydrogen Diffusion Flames" AIAA Paper 91-1780, Honolulu, Hawaii,

June, 1991.

3
Center for Modeling of Turbulence and Transition

Research Briefs - 1990

TI irbulence Modeling

Tsan-Hsing Shili

1. Motivation and Objectives

(1) Examine the performance of existing two-equation eddy viscosity models and develop better models for the near-wall turbulence using direct

numerical simulations of plane channel and boundary layer flours.

(2) Use the asymptotic near-wall behavior of turbulence to examine the

problems of current second-order closure models and develop new models

with the correct near-wall behavior.

(3) Use Rapid Distortion Theory to analytically study the effects of mean

deformation (especially due to pure rotation) on turbulence, obtain

analytical solutions for the spectrum tensor, Reynolds stress tensor,

anisotropy tensor and its invariants, which can be used in the turbulence model development.

(4) Explore the potential of the renormalization group (RNG) theory in

turbulence modeling.

(5) Modeling of compressible turbulent flows.

(6) Modeling of bypass transition.

2. Work Accomplished and Ongoing Work

2.1 k-e model

The k-e model is still the most widely used model for computing engineering flows. We have examined the near-wall behavior of various eddy

viscosity models proposed by different researchers, and have studied the

near-wall behavior of the terms in the k-equation budget. We found that

the modeled eddy viscosity in many existing k-e models does not possess

correct near-wall behavior and the pressure transport term in the k-equation

is not modeled appropriately. Based on the near-wall asymptotic behavior of

the eddy viscosity and the pressure transport term in the k-equation, a new

set of improved closure models has been obtained. In addition, a modeled

equation for the dissipation rate is derived more rationally. This work is

reported in NASA TM 103221 ICOMP-90-16111.

In addition, all the existing two-equation models (except Jones & Laun4
der model, which unfortunately does not work well even for some simple

flows) have an "unacceptable" wall distance parameter (y+ ) in their eddy

viscosity damping function fμ (y + ). This will result in an unphysical zero

eddy viscosity near the separation region. In addition, y + can not be well

defined in many flows with complex geometry. To remove this deficiency,

Dr. V. Michelassi, Dr. A. Hsu and I proposed two new eddy viscosity damping functions, and both of them are independent of the wall distance. The

new models have been satisfactorily tested in channel and boundary layer

flows. This work is reported in two papers: AIAA-91-0611 121 and NASA

TM/ICOMP/CMOTT[sl

2.2 Second order modeling of near-wall turbulence

The main emphasis is on developing a near-wall turbulence model for

the velocity pressure gradient correlation and the dissipation tensor in the

Reynolds-stress equation. A modeled dissipation rate equation is also derived

more rationally. Near a wall, a reduction in velocity fluctuations normal to

the wall becomes significant. Because of this wall effect, the viscous diffusion

term in the Reynolds-stress equations becomes the leading term and it must

be properly balanced by the other terms. We have used this as a model

constraint for developing a model for the pressure and dissipation terms. To

test the models, a fully developed channel flow and boundary layer flows

are chosen as the test flows, for which direct numerical simulations and experimental data are available for comparison. The modeled Reynolds stress

equations for the channel flow are steady one-dimensional, and for boundary layer flows are steady two-dimensional. Therefore model testing will be

very accurate. This part of work ['] is reported in the paper: Proceedings

of the International Symposium on Engineering Turbulence Modeling and

Measurements and NASA TM 103222 ICOMP-90-0017.

2.3 Second order modeling of a three - dimensional boundary layer

A study of three-dimensional effects on turbulent boundary layer was

achieved by the direct numerical simulation of a fully developed turbulent

channel flow subjected to transverse pressure gradient (see Physics of Fluids, Vol.2 N0.10, 1990, pp. 1846-1853). The time evolution of the flow was

studied. The results show that, in agreement with experimental data, the

Reynolds stresses are reduced with increasing three-dimensionality and that,

near the wall, a lag develops between the stress and the strain rate. In addition, we found that the turbulent kinetic energy also decreased. To model

these three-dimensional effects on the turbulence, we have tried different

two-equation models and second order closure models. None of the current

5
closure models can predict the reductions in the shear stress and turbulent

kinetic energy observed in direct numerical simulations. Detailed studies

of the Reynolds-stresses budgets were carried out. One of the preliminary

conclusions from these budget studies is that the velocity pressure-gradient

term in the normal stress equation W) plays a dominant role in the reduction of shear stress and kinetic energy. These budgets have been used

to guide the development of better models for three dimensional turbulent

boundary layer flows. This work N was presented in the American Physical

Society Forty-Third Annual Meeting, November, 1990.

2.4 The effect of rotation on turbulence

In addition to the above studies of second order closure models, we have

carried out some RDT analysis on simple homogeneous turbulent flows. An

order of magnitude analysis shoes that under the condition of S(q2 )1E »

Rt , the equations for turbulent velocity fluctuations can be approximated

by a linear set of equations, and if S(q2 )/E » R3/4 , then the turbulent velocity equations can be further approximated by an inviscid linear equation.

Therefore, RDT can be used to analytically study some very basic turbulent

flows, such as, homogeneous shear flows, irrotational strain flows and pure rotational flows. This work focuses on the effect of rapid rotation on turbulence

using RDT. We have obtained analytical expressions for velocity, the spectrum tensor, Reynolds-stress, the anisotropy tensor and its invariants. The

solutions show that the turbulence is strongly affected by the rapid rotation.

Using RDT, we can calculate the rapid pressure-stain term exactly and we

can obtain very useful information for developing corresponding turbulence

models. See the report [ '] for this work.

2.5 Renormalization Group Theory (RNG) in turbulence modeling

RNG method has been introduced to the turbulence modeling mainly in

the Large Eddy Simulation (LES) of turbulence with a subgrid scale model.

One also attempted to use it to develop Reynolds-averaged turbulence model

equations, for example, k-E model equations. However, we found that there

are a few fundamental concepts and important procedures used in the derivation of those model equations which are not clear and well justified. Dr. Z.

Yang and I are working on this subject and try to explore the potential of

RNG in the turbulence modeling.

2.6 Modeling of compressible turbulent flows

The turbulence models for compressible flows are of great interest in hypersonic flows and turbulent combustion. The modeling scheme greatly de6
pends on the averaging schemes (i.e., conventional average, density weighted

average and mixed average) used in the turbulence equations. We start with

the analysis of the turbulent equations derived from the different averaging

schemes to see what kind of averaging scheme is most convenient for both

turbulent modeling and applications in CFD. We concentrate on the second

order closure model (i.e. Reynolds stress model) and two-equation model.

Dr. W. Liou and I are working on this subject. See Reference N for the first

report on averaging schemes for compressible flows.

2.7 Modeling of bypass transition

Most common transition phenomena occurred in engineering flows are

bypass transition. A few papers on modeling of transition with turbulence

models show that the bypass transition can possibly be modeled with the

modified turbulence models developed solely for turbulent flows. However,

most of the work in this direction was based on the parabolic two-equation

models. We expect that the bypass transition phenomena will be more appropriately described by the elliptical equations. Then, the prediction of

normal stresses becomes important. Because of the inability of modeling

normal stresses with the two-equation models, we are pursuing the elliptical Reynolds stress model equations for the bypass transition studies. Dr.

Z. Yang and I are working on the improvement of our previous near-wall

Reynolds stress model for the purpose of modeling bypass transition.

2.8 Modeling of scalar turbulence:

Modeling of scalar turbulence is of great importance in turbulent heat

transfer. Eddy viscosity models often fail in the prediction of heat transfer in many shear flows. We have developed a set of second order closure

models based on the joint realizability (between velocity and scalar) and the

experiments. Dr. A. Shabbir and I are working on this subject. A paperN

was presented in the Lumley Symposium: Recent developments in turbulence,

November, 1990.

3. Publications:

1. Shill, T.-H., 1990, "An Improved k-e Model for Near-Wall Turbulence

and Comparison with Direct Numerical Simulation," NASA TM 103221

ICOMP-90-16.

2. Shih, T.-H. and Hsu, A.T., 1991, "An Improved k-e Model for Near-Wall

Turbulence," AIAA-91-0611.

7
3. Michelassi, V. and Shih, T.-H., 1991, "Low Reynolds Number TwoEquation Modeling of 'Eirbulent Flows," NASA TM 104368, ICOMP91-06, CMOTT-91-01.

4. Shih, T.-H. and Mansour, N.N., 1990, "Modeling of Near-Wall Turbulence," Proceeding of the International Symposium on Engineering Turbulence Modeling and Measurements, September, 1990, Dubrovnic, Yugoslavia, Editors: W. Rodi, E.N. Gami c. or, NASA TM 103222 ICOMP90-0017.

5. Shih, T.-H., 1990, "Modeling of 3D Turbulent Boundary Layer Flows,"

American Physical Society Forty-Third Annual Meeting, 1990, Cornell,

U. Ithaca, New York.

6. Shih, T.-H., 1991, "Rapid Distortion Theory on Homogenous Turbulence with Rapid Rotation," CMOTT Report.

7. Liou, W.W. and Shih, T.-H., 1991, "On the Basic Equations for the

Second-order Modeling of Compressible Turbulence," CMOTT-91-06.

8. Shih, T.-H. and Shabbir, A., 1990, "Advances in Modeling the Pressure

Correlation Terms in the Second Order Moment Equations," the Lumley Symposium: Recent developments in turbulence, November, 1990,

ICASE, NASA Langley Research Center, Edited by T.B. Gatski, S.Sarkar

and C. G. Speziale

9. Shih, T.-H., 1990, "Advancements in Engineering Turbulence Modeling," 9th NASP Technology Symposium, Paper-105, November, 1990,

Orlando FL.

10. Shih, T.-H., Chen, J.-Y. and Lumley, J.L., 1991, "Second Order Modeling of Boundary Free Turbulent Shear Flows," AIAA 91-1779.

8
Center for Modeling of Turhulenr:c and Transition

Research Briefs - 1990

Experiments and Modeling

Aamir Shabbir

1. Motivation and Objectives

The usual approach in establishing the correctness and accuracy of turbulence models is to numerically solve the modeled differential equations

and then compare the results with the experiment. However, in the case

of a discrepancy, this procedure does not pinpoint where in the model the

drawback lies. It is also possible that the model overcompensates one physical phenomenon and undercompensates the other so that the net result is

a good agreement between the two. Therefore a more desirable approach is

to directly compare the individual terms in the equations with their models. To achieve this objective primarily physical experiments have been used

to carry out the second moment budgets. These can then be used to analyze and assess various models and closure assumptions and seek improvements/modifications where models prove deficient.

2. Work Accomplished

2.1 Evaluation and Development of Turbulence Models for Pressure Correlations.

A direct comparison between the pressure strain and pressure temperaturegradient correlations and their closure models is carried out. The flows used

include both physical and numerical experiments on homogeneous shear flows

and physical experiments on buoyant plumes. Models considered include

both the linear and the more elaborate non-linear ones. It is found that

the non-linear models provide a much better agreement with these experiments than the linear ones. A new model for the slow part of the pressure

temperature-gradient correlation is also derived using joint realizability concept.

2.2 On the Ratio of Mechanical to Thermal Time Scales in Turbulent Flows.

The ratio of these time scales is very often employed in the two equation

turbulence models. The study of Beguier et al (1978) recommended this

9
ratio to be around 2.0. The current analysis using the buoyant plume and

homogeneous shear flow experiments shows that this value is about 3.0. It

is shown that this departure from the commonly used value is a consequence

of the local equilibrium assumption being not satisfied by these experiments.

2.3 Turbulent Buoyant Transport - A Comparative Study between

Models and Experiments.

The more popular gradient diffusion type models for the turbulent transport (third moments) were found to underestimate the experiment by an

order of magnitude. More complex models (Andre et al 1976, Lumley et

al 1978), based on the simplification of exact transport equations of third

moments, do a much better job in reproducing the third moments although

the results are still less than satisfactory.

2.4 Experimental Balances of Second Moment Equations for a

Buoyant Plume.

Despite large volume of work on second moment closure, there is very

little experimental information available about the budgets of the second

moments. Part of this reason stems from our inability, at present, to measure the pressure correlations. Experimental budgets for Reynolds stresses

and heat fluxes have been carried out for a boundary free shear flow (round

plume) and the pressure correlations are obtained as the closing terms in

these budgets. These budgets show how different terms in the equations are

distributed across the flow and can be used to analyze some of the modeling

assumptions. For example they show that the assumption of local equilibrium is not justified for bulk of the flow field - an idea fundamental to the

algebraic stress models.

2.5 X-wire Response in Turbulent Flows of High Intensity Turbulence and Low Mean Velocities.

This work is based on an experimental study, which was carried out at

SUNY/Buffalo, of angular response of an x-wire, at low velocities (0.25m/s

to lm/s). It is found that the k-factor in the modified Cosine Law is strongly

velocity dependent. The implications of this on multi component turbulence

measurements are explored. Expressions are also derived for evaluating when

the cross-flow errors begin to affect x-wire measurements.

10
2.6 Modeling of Turbomachinery Flows using the Average Passage

Approach.

Turbomachineray flows are turbulent and unsteady and numerical calculation of a flow in a multistage machine, at present, is not possible. However,

the effects of periodic unsteadiness can be accounted through the models for

deterministic stresses which arise in the average passage equation set (Adamczyk 1985). Exact equations governing the transport of these stresses have

been derived and a two equation model is being developed and tested at

present. The model uses ideas from turbulence modeling such as the gradient diffusion type hypotheses. This work is being performed in collaboration

with J. Adamczyk of the Lewis Research Academy, at the NASA Lewis Research center.

3. Future Plans

3.1 Study the effect of buoyancy on turbulence by computating flows using

turbulence models. In addition to environmental flows, such a work also

has industrial applications e.g. cooling of nuclear reactors and electronic

components, and "geyser" formation in fuel tanks in microgravity.

3.2 Seek improvements in the models for turbulent transport. In general

the transport is not too important in most of the turbulent flows but in

some applications, e.g. geophysical flows, the modeling of the transport

could be critical in the success of a computation.

3.3 Seek improvements in the existing two equation models by incorporating

newer models for pressure correlations etc.

3.4 Assess the models for deterministic stresses in a multistage turbomachinery environment.

4. Publications

1. Pressure Correlations in the Reynolds Stress and Heat Flux Equations

- A Comparison between Experiment and Models. A. Shabbir. APS

Bulletin, Vol. 35, No. 10, Nov. 90, Abstract KC4.

2. Evaluation of Turbulence Models for Predicting Buoyant Flows. A.

Shabbir and D.B. Taulbee. J. Heat Transfer, 1990, Vol 112, No 4, pp

945-951.

3. Experiments on Round Turbulent Buoyant Plumes. A. Shabbir and
W.K. George. Under review for publication in J. Fluid Mechanics.

4. Advances in Modeling Pressure Correlation Terms in the Second Moment Equations. T.-H. Shih and A. Shabbir. Presented at the Symposium honoring J. Lumley's 60th birthday, November 90, NASA Langley

Research Center.

5. X-wire Response in Turbulent Flows of High Intensity Turbulence and

Low Mean Velocities. A. Shabbir, P.D. Beuther, and W.K. George.

Submitted to Experimental Thermal and Fluid Science.

5. References

Ada-Tnczyk, J. J. "Model Equation for Simulating Flows in Multistage

Turbomachinery", ASME Paper No. 85-GT-226. (1985)

Andre J.C., G. De Moor, P. Lacarrere and R. du Vachat "Turbulence

Approximation for Inhomogeneous Flows: Part I. The Clipping Approximation", J. Atmos. Sci., Vol. 33, pp. 476-481, (1976)

Beguier, C., 1. Dekeyser and B. E. Launder "Ratio of Scalar and Velocity

Dissipation Time Scales in Shear Flow Turbulence" ," Phys. Fluids, Vol.

21, pp. 307-310 (1978).

Lumley, J. L., 0. Zeman and J. Siess "The influence of Buoyancy on

Turbulent Transport", J. Fluid Mech., Vol. 84, pp. 581-597 (1978).

12
Center for Modeling of Turbulence and Transition

Research Briefs - 1990

Modeling of Compressible Turbulent Shear Flows

William W. Lion

1. Motivation and Objectives

Despite all the recent development in computer technologies and numerical algorithms, full numerical simulations of turbulent flows are feasible only

at moderate Reynolds numbers and for flows with relative simple geometries.

Turbulence models provide alternatives in the pressing need for the prediction

of turbulent flows and,in fact, have become an important pacing factor for

the successful development of computational fluid dynamics (CFD). With

the advent of supercomputers, however, it has become more affordable to

apply second order closure models in the prediction of flows with complex

effects; such as strong curvature, three-dimensionality and compressibility.

The main goal of this research is to develop new second order moment closures for compressible turbulence. It has been shown that the models based

on the extension of those developed originally for incompressible flows fail

to predict adequately turbulent flows at high Mach numbers. In this attempt, the compressibility effects will be explicitly considered. A successful

development of these models that take into account directly the compressibility effects may have a range of technological implications in the design of

supersonic and hypersonic vehicles.

2. Work Accomplished

During this early stage of the task, the goal is to obtain an objective yet

comprehensive understanding of the development and the current status of

compressible turbulence modeling. Due to the variable density effects in compressible flows, density correlation terms appear in the governing equations

for the mean flow, if the conventional ensemble averaging technique is applied. These terms do not exist in incompressible flows and need to be modeled. On the other hand, the mass-weighted-averaging or Favre procedure

generates a set of mean equations that have the similar forms as they are for

incompressible flows. One may then incline to use the incompressible analog

in compressible flow calculations. A simple test, however, should show that a

direct application of the exact incompressible models fails. One of the main

effects that is excluded in incompressible models is the finite propagation

13
speed of disturbances. In compressible flows, modulation of flow properties

occurs only within Mach cones of influence, with acoustic time delay. This

introduces additional scales for the transport properties. Caution also needs

to be exercised in comparing Favre-averaged calculations with experiments,

since the differences between Favre-averaged and measured quantities may

not be negligible at high Mach numbers.

Compressible turbulence modeling is still in its infancy. This appears

to be true both theoretically and computationally. Recently, a new concept

called dilatation dissipation was proposed. Dilatation dissipation, as opposed to the solenoidal dissipation in incompressible flows, accounts for the

viscous dissipation of turbulent kinetic energy due to volume fluctuations.

Models accommodating this effect show the importance of this additional

drain of turbulent kinet-ic energy in order to obtain adequate model predictions. Dilatation dissipation appears to be among the direct consequencies

of compressibility effects.

Another school of thought on compressibility effects focuses on the

changes of turbulence structures at high Mach numbers. Note that these

structures are identified by using conditional sampling techniques in experiments. Due to the communicability problem between interacting elements,

structures that are highly efficient in extracting energy from the mean flow at

low Mach numbers no longer prevail as the Mach number increases. They are

replaced by structures that are less sensitive to the Mach number. This selective amplifying behavior is describable by quasi-linear theory, which view

the turbulence energetics as physical manifestations of ongoing nonlinear

instability in turbulent shear flows.

The above mentioned matters are described in detail in an ICOMP/

CMOTT report [1] that is in preparation. Equations for the second order

moments and the mean flotiv as a result of the application of different averages

will also be given. Modeling methodologies used in compressible flow calculations will be reviewed. The evaluation process performed during the present

stage of the research has identified avenues that will be pursued during the

next period of this task.

3. ]Future Plans

(1) Develop second order models for compressible turbulence based on ensemble averages. This may be assisted by first developing k — E types of

models to identify important mechanisms.

(2) Develop unconventional models that incorporate explicitly the characteristics of the structures of compressible turbulence.

14
The developed models will be applied to certain benchmark flows with and/or

without chemical reactions.

4. Publications

1. Liou, W. W. and Shih, T.-H., "On the Basic Equations for the Secondorder Modeling of Compressible Turbulence," CMOTT-91-06, 1991.

2. Liou, W. W. and Morris, P. J., "An Comparison of Numerical Methods

for the Rayleigh Equation in Unbounded Domains," CMOTT-91-05,

1991.

15
Center for Modeling of Turbulence and Transition

Research Briefs - 1990

RNG in Turbulence and Modeling of Bypass 'Transition

Z higang Yang

1. Motivation and Objectives

Since I joined CMOTT on July 1990, I have been working on two

research projects. The first project concerns the Renormalization Group

(RNG) analysis of turbulence and the second project is on the calculation of

bypass transition through turbulence modeling. In addition, the preparation

of two papers on work performed at was completed.

Application of RNG in turbulence was proposed by Yakhot and Orszag

in 1986. RNG is a process which eliminates systematically the small scales,

and represents the effect of those eliminated small scales on the uneliminated

large scales as the changes in the transport properties. It is because of this

property of RNG that Yakhot and Orszag suggested that RNG could be

used as a model builder in turbulence modeling. They also presented a

k — c model in their 1986 paper. However, this paper is lengthy, with many

unstated assumptions. Our aim is first to understand, and to validate the

RNG approach in turbulence through an independent study. We will then

study the possibility of constructing RNG based turbulence models, and try

to proceed to do the turbulence modeling through RNG in parallel with the

classical approach. We will also compare the numerical predictions made by

RNG models and by classical models against data from Direct Numerical

Simulation and against experimental data from different benchmark cases.

In a quiescent environment, the transition is initiated by the instability

of the laminar boundary layer to Tollmien-Schlichting waves. These waves

are amplified with streamwise distance and eventually breakdown into turbulent spots, which are precursors of turbulent boundary layer. While in an

environment with high freestream turbulence, the transition is found to be a

bypass one in which turbulent spots are formed without Tollmien-Schlichting

wave amplification. The formation of turbulent spot is a random process,

and flow within a turbulent spot is almost fully turbulent. This suggests the

possibility of using turbulence modeling to describe and predict the bypass

transition. There have been some works in this direction, primarily using

different versions of two equation models. Bypass transition is predicted, as

the level of the freestream turbulence is increased. However, it is found that

the predicted transition is much sharper than that observed in the experi16
ment. In addition, the predicted transition depends on the description of the

initial profiles to an certain extent. The works we propose to do are twofold.

1) We will be using a low Reynolds number version of Reynolds stress model

rather than the two equation model. This would bring in more physics, and

hopefully would fit better the complicated flows such as bypass transition.

2) We will be using an elliptic solver rather than a parabolic solver for this

boundary layer transition. This way, we will be able to include the effect of

the leading edge. The testing case will be flow passing a flat plate. Both the

zero pressure gradient case and the non-zero pressure gradient case will be

tested.

2. Work Accomplished

1. Nonlinear dynamics near the stability margin in rotating pipe flow.

The nonlinear evolution of inarginally unstable wave packets are studied

in rotating pipe flow. These flows depend on two control parameters, which

may be taken to be the axial Reynolds number Re and the rotation rate q.

Marginal stability is realized on a curve in the (Re, q) plane, and we explore

the entire marginal stability boundary. As the floe- passes through any point

on the marginal stability curve, it undergoes a supercritical Hopf bifurcation

and the steady base flow is replaced by a traveling wave. The envelope of

the wave system is governed by a complex Ginzburg-Landau equation. The

Ginzburg-Landau equation admits Stokes waves, which correspond to standing modulations of the linear traveling wavetrain, as well as traveling wave

modulations of the linear wavetrain. Bands of wavenumbers are identified in

which the nonlinear modulated waves are subject to a side-band instability.

This work[1] was reported in APS/DFD meeting in November 1990. A

paper for this work has been submitted to JFM for consideration of publication. The paper is co-authored with S. Leibovich of Cornell University. This

work was supported by AFOSR-89-0346.

2. Unstable viscous wall modes in rotating pipe flow. (AIAA Paper No.

91-1801)

Linear stability of flow in rotating pipe is studied. These flows depends

on two parameters, which can be taken as the axial Reynolds number Re and

the rotating rate q. In the region of Re » 1 and q = 0(1), the most unstable

modes are wall modes. The wall modes are found to satisfy a simpler set of

equations containing two parameters rather than four parameters as in the

full linear stability problem. The set of equations is solved numerically and

asymptotically over a wide range of the parameters. In the limit of Re —> oc,

17
the eigenvalue reaches the inviscid limit and the eigenfunction shows a two

layer structure. The eigenfunction reaches the inviscid limit over the main

part of the domain, while near the wall of the pipe, the eigenfunction is

represented by a viscous solution of the boundary layer type.

This work[2] is to be presented at the AiAA 22nd Fluid Dynamics,

Plasma Physics and Lasers Conferences, June 24-26, 1991. The paper is coauthored with S. Leibovich of Cornell University. This work was supported

by AFOSR-89-0346.

3. RNG in turbulence modeling.

This work is done in collaboration with Dr. T.H. Shih of CMOTT. In

this work, we carry out an independent study of the work done in the paper

by Yakhot and Orszag up to their derivation of k — E equation. Many of their

results are repeated. However, we also found some discrepances, and some

unclaimed assumptions in their derivations.

3. Future Plans

1. Currently, we are testing and improving the low Reynolds number

version of Reynolds stress model proposed by Shih and Mansour (1990) in

the simple shear flows, such as channel flow and boundary layer flow. This

model is going to be used in the calculation of bypass transition.

2. We will carry an independent derivation of k — E equation using RNG,

and compare it with the one presented by Yakhot and Orszag. We will also

compare the prediction of RNG k — E model with the other k — E models, in

both the high Reynolds number case and the low Reynolds number case.

4. Publications

[1] Yang, Z. and Lelbovich, S. 1990 "Nonlinear dynamics near the stability

margin in rotating pipe flow", Submitted to J. Fluid Meth.

[2] Yang, Z. and Leibovich, S. 1990 "Unstable viscous wall modes in rotating

pipe flow", AIAA Paper 91-1801.

5. References

Yakhot, V. and Orszag, S.A. 1986 Renormalization group analysis of

turbulence. J. Sci. Comput. Vol. 1, No. 1, 3-51.

Shih, T.H. and Mansour, N.N. 1990 Modeling of near wall turbulence.

NASA TM-103222, ICOMP-90-0017.

18
Center for Modeling of Turbulence and Transition

Research Briefs - 1990

Turbulence Modeling in Supersonic Combusting Flows

Tawit Chitsomboon

1. Motivation and Objectives

To support the National Aerospace Plane project, the RPLUS3D CFD

code has been developed at the NASA Lewis Research Center. The code

has the capability to solve three dimensional flowfields with finite rate combustion of hydrogen and air. The combustion processes of the hydrogenair system are simulated by an 18-reaction path, 8-species chemical kinetic

mechanism. The code uses a Lower-Upper (LU) decomposition numerical

algorithm as its basis, making it a very efficient and robust code. Except

for the Jacobia.n matrix for the implicit chemistry source terms, there is no

inversion of a matrix even though it uses a fully implicit numerical algorithm.

The main purpose of this work is to incorporate a k -E (two equation)

turbulence model into the RPLUS31) code.

2. Work Accomplished

Since February 1990, when this work was started, some of the more

important accomplishments are categorized as follows:

1) Add a k-e turbulence model: The model selected is in a high Reynolds

number form. The low Reynolds number form could not be used economically in the case of a three dimensional flow with chemical reaction since it

demands too much in computer resources. The addition was designed to be

as modular as possible but some interactions with the main code are needed

in order to be more efficient. The first test case tried was a Mach 0.5 floeover a flat plate. The velocity profile compare very well with the log-law.

profile. The friction coefficient also compares well with the Van Driest correlation. More validations will be performed for other flows such as free shear

layer and jet flows.

2) Improved accuracy and convergence rate: According to a stability

analysis of a model equation, it is shown that the RPLUS3D code excessively added artificial damping to the right hand side of the algorithm while

at the same time it overestimated the spectral radii on the left hand side.

These excessive additions would not give an optimum convergent rate. Mod19
ifications were made such that true directional spectral radii were added to

the left hand side and the artificial damping terms were reduced to optimum

values in accord with the stability analysis.

A test run was made of a Mach 4 flow of air over a 10 degree compression

ramp. It was found that the modified code converged to machine zero about

five times faster than the original code while at the same time it somewhat

improved the shock resolution.

3) Added consistent damping terms at block interfaces: It had been observed that at the block interface of a multi-block grid, wiggles developed.

This happened because the damping terms at the interface were not consistent with those at the interior points.

4) Validated RPLUS31) with laminar flows over flat plates: The test

cases run were for Mach numbers of 0.1, 0.3 and 0.5. Results of all cases

seemed to be good except for the region of high curvature of the velocity

profile near the edge of the boundary layer.

5) Add two dimensional capability to the code: It turned out that this

is not a trivial task especially for a finite volume code like RPLUS31). It is

nice that now the code can solve a 2D flow without having to carry the 3rd

direction along as a redundancy. There is, of course, no need to maintain a

separate 2D code.

6) Implemented a local time stepping capability: The original code always ran at a fixed time step of 1 second. For most flows this corresponds to

using a very large CFL number which may not be conducive to a fast convergence rate. With the local time stepping, it was found that an optimum

CFL number was, in agreement with other investigators, around 5 to 7.

7) Implemented implicit boundary conditions: This addition enhanced

the convergence rate by about 30 percent at the expense of a more complex

code and an increase of about 20 percent in CPU time per iteration step.

There seems to be, then, no net advantage of the implicit boundary condition

except maybe in the area of robustness.

8) Changed input file: Instead of having to scan a whole subroutine to

set up a problem, a user now can change numbers in a small file of length

about one page. To start up a run from a previous run one now needs only to

change a parameter in this input file without having to recompile the input

subroutine as before.

9) Solved two 3D hypermixing flow fields of W.Hingst and D.Davis: This

work was performed in collaboration with Dr. A.C. Taylor of Old Dominion

University. My task was to set up the program for these particular problems

20
and to generate the grid for one of the problems.

3. Future Plans

The work plan for the year 1991 consists of both basic model implementations and practical applications of the code :

• Continue to validate the baseline k-E model.

• Add compressibility effects to the base line model.

• Apply the base line model to re-solve the 3D flowfields mentioned in the

previous section.

21
Center for Modeling of Turbulence and Transition

Research Briefs - 1990

Development of New Flux Splitting Schemes

Meng-S. Liou and Christopher J. Steffen, Jr.

1. Motivation and Objectives

Maximizing both accuracy and efficiency has been the primary objective

in designing a numerical algorithm for computational fluid dynamics (CFD).

This is especially important for solution of complex 3D systems of NavierStokes equations which often include turbulence modeling and chemistry

effects. Recently, upwind schemes have beer. well received for both their

capability of resolving discontinuities and their sound theoretical basis in

characteristic theory for hyperbolic systems. With this in mind, we present

two new flux splitting techniques for upwind differencing.

2. Work Accomplished

The first method is based upon High-Order Polynomial Expansions

(HOPE) of the mass flux vector [1]. The present splitting results in positive

and negative mass flux components that vanish at M=O. Thus the error in

the Van Leer scheme which results in the diffusion of the boundary layer is

eliminated. We also introduce several choices for splitting the pressure and

examine their effects on the solution.

The second new flux splitting is based on the Advection Upwind Splitting Method (AUSM for short) [2]. In Navier-Stokes calculations, the diffusion error present in Van Leer's flux splitting scheme corrupts the velocity

vector near the wall. In the AUSM, a proper splitting of the advective velocity component leads to an accurate resolution of the interface fluxes. The

interface velocity is defined using the Mach number polynomial expansion

in the mass flux, then the convective fluxes follow directly. Again, several

choices of pressure splitting axe possible among which a simple Mach number splitting according to characteristics appears to be the best in terms of

accuracy. The scheme has yielded results whose accuracy rivals, and in some

cases surpasses that of Roe's method, at reduced complexity and computational effort. The calculation of the hypersonic conical flow demonstrates

the accuracy of the splittings in resolving the flow in the presence of strong

gradients. The second series of tests involving the 2D inviscid flow over a

NACA 0012 airfoil demonstrate the ability of the AUSM to resolve the shock

22
discontinuity at transonic speed and the level of entropy generation at the

stagnation point.

In the third case we calculate a series of supersonic flows over a circular cylinder. The Roe splitting in all conditions and grids tested yielded

anomalous solutions (sometimes referred to as the carbuncle phenomenon),

which could appear as non-symmetric, protuberant, or indented contours.

The AUSM gave expected solutions in all calculations.

The fourth test deals with a 2D shock wave/boundary layer interaction.

This provides an opportunity to accurately resolve a laminar separation region and to compare the ability to resolve a non grid-aligned shock with

other methods.

3. Future Plans

Future plans are primarily concerned with the AUSM. A detailed stability analysis for this new technique will be useful. The idea of splitting the

ddvective velocity opens up a whole family of potential schemes. Therefore, a

comprehensive study of the interaction between various pressure and advection velocity splitting methods is necessary to optimize both accuracy and

efficiency. Additionally, a 2D turbulent calculation would be a good test of

the scheme's ability to solve a coupled system of K — e equations.

4. Publications

1. Liou, M.-S. and Steffen, C.J.Jr., "High-Order Polynomial Expansions

(HOPE) for Flux-Vector Splitting," (to be presented at the ICES191

Conference, August 11-16, 1991)

2. Liou, M.-S. and Steffen, C.J.Jr., "Development of a New Flux Splitting

Scheme," (to be presented at the AIAA Tenth CFD Conference, June

24-26,1991)

23
Centcr for Modeling of Turbulence and Transition

Research Briefs - 1990

Progress of Simulations for Reacting Shear Layers

Sheng-Tao Yu

1. Motivation and Objectives

In the past six months, the effort was devoted to the development of

a high speed, chemically reactive shear layer test rig. The purpose of the

experiment is to study the mixing of oxidizer and fuel streams in reacting

shear layers for various density, velocity, and Mach number. The primary

goal is to understand the effects of the compressibility upon mixing and

combustion in a fundamental way. Therefore, a two-dimensional shear layer

is highly desirable for its simplicity to quantify the compressibility effects.

The facility consists of a two-stream wind tunnel with two independent gas supplies. After passing through flow-management devices located

upstream, each gas stream expands to its predetermined Mach number by

means of a contoured center body and tunnel walls. Various combinations

of flow conditions of high-speed stream and low-speed stream allows for the

systematic study of mixing and reactions of compressible shear layers.

2. Work Accomplished

The RPLUS 2D code is used to calculate the flow fields of different

sections of the test rig. The emphasis was on the supersonic nozzle design,

the vitiation process for the hot air stream and the overall thermodynamic

conditions of the test matrix.

The k — e turbulence model with wall function has been successfully

implemented in the RPLUS code. The k and e equations are solved simultaneously and the the LU scheme is used to make it compatible with the flow

solver. The coupling between the flow solver and the k — E solver depends

on the turbulence viscosity only, and the k — e solver is separated from the

flow solver to reduce the complexity. The newly developed code has been

used for the compressible shear layer calculations. Many cases of the compressible free shear layer with various convective Mach numbers and density

ratios have been simulated using the compressible k — e solver. The results

are summarized in two technical papers.' , ' Currently, the k — E solver is a

standard feature in the RPLUS 2D code and the code has been distributed

to the industry and universities through NASP group. Locally, Duncan and

Tsai are using the k — e solver for their research work.

24
3. Future Plans

Physical phenomena of the reactive free shear layer can not be adequately described by the k — E model coupled with the Reynolds averaged

flow equations. The properties of the vortical flows which dominate the

whole flow field of free shear layers can only be illustrated by time marching

numerical method with accurate spatial resolution. Traditionally, the spectrum methods were used for this kind of applications. However, very limited

success has been reported for compressible flows using spectrum methods.

On the other hand, recent development show promising results using high

order central differencing and Essentially Non-Oscillatory (ENO) schemes.

One developed by Lele of Stanford university using high order compact differencing is especially interesting. Future work includes Direct Numerical

Simulation (DNS) of Navier Stokes equations for the chemically reactive

flow and application to free shear layers.

In additional to the above mentioned work, I will serve as a consultant

for the k — E solver in the RPLUS code.

4. Publications

1. S. T. Yu, J. S. Shuen, and Y-L P. Tsai, "Three Dimensional Calculations

of Supersonic Combustion Using a LU Scheme," to appear in the J. of

Comput. Phys.

2. S.T. Yu, C.L. Chang, and C.L. Merkle, "Solar Rocket Plume/Mirror

Interactions, " submitted to the J. of Spacecraft and Rockets.

3. S.T. Yu, B.J. McBride, K.C. Hsieh, and J.S. Shuen, "Hypersonic Flows

Simulations with Equilibrium or Finite Rate Chemistry," submitted to

Computers & Fluids for publication

4. S.T. Yu, "A Convenient Way to Convert 2D CFD Codes to Axisymmetric Ones," submitted to J. of Propulsion and Power as a technical

note.

5. S.T. Yu, C.D. Chang, and C.J. Marek, "Modern CFD applications for

the Design of a Reacting Shear Layer Facility," presented at the AIAA

Science Meeting, 1991.

6. S.T. Yu, C.D. Chang, and C.J. Marek, "Simulation of Free Shear Layers

Using a Compressible k-E Model," accepted for presentation at the AIAA

Propulsion Conference, 1991.

25
Center for Modeling of Turbulence and Transition

Research Briefs - 1990

Appendix A

Organization - 1990

Position Chart

ICOMP '` Research Acad.

LM Aer OVirectorate

26
Center for Modeling of Turbulence and Transition

Research Briefs - 1990

Manager Dr. Louis A. Povinelli

Deputy Chief

Internal Fluid Mechanics Division

NASA Lewis Research Center

Coordinator Dr. Meng-Sing Liou

Senior Scientist

Internal Fluid Mechanics Division

NASA Lewis Research Center

Technical Leader Dr. Tsan-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

Cornell University

Professor Eli Reshotko

Department of Mechanical and

Aerospace Engineering

Case Western Reserve University

27
Center f.;r Modeling of Turbulence and Transition

Research Briefs - 1990

Shabbir, Aamir

5/90 - present

Shih, Tsan-Hsing

5/90 - present

Steffen, Christopher J. Jr

10/90 - present

Yang, Zhigang

7/90 - present

Yu, Sheng-Tao

3/90 - present

Research Areas

ICOMP Algorithms,

Code Development

Sverdrup Tech., Inc. PDF Turbulence Modeling

NASA Intern Turbulence Modeling,

Univ. of Michigan CFD

NASA LeRc CFD Algorithms,

High-speed Flow

ICOMP Compressible Flow

Modeling, Weakly

Nonlinear Wave Models

NRC Buoyancy Effects on

Turbulence

ICOMP Turbulence Modeling

NASA LeRc Upwind Algorithms

Incompressible Flow

ICOMP RNG in Turbulence

Modeling, Modeling of

Bypass Transition

Sverdrup Tech., Inc. Modeling of Chemical

Reacting Flows,

Direct Numerical Simulations

Member Listing

Names/Term

Chitsomboon, Tawit

5/90 - present

Hsu, Andrew T.

5/90 - present

Lang, Nancy J.

5/90 - 8/90

Liou, Meng-Sing

5/90 - present

Liou, William W.

11/90 - present

28
Center for Modeling of Turbulence and Transition

Research Briefs - 1990

Appendix B

Seminars and Technical Meetings

B.1 CMOTT Seminar Seminars

The purpose of these seminars is to exchange ideas and opinions about

the latest developments and current state of turbulence and transition research. The speakers are invited from within and outside of the NASA

LeRC, including foreign speakers. This seminar series complements the informal CMOTT technical group meetings.

The abstracts of the seminars are given below.

29
I'-,-(DR COh^

INSTITUTE FOR COMPUTATIONAL ,,, ^. Q.oi mn O

z ICOMP

s Obd N 1 z'% MECHANICS IN PROPULSION

(COMP SEMINAR SERIES

Measurements in a Circular Jet of Helium into

Quiescent Air*

by

John L Lumley

Cornell University

Monday, June 4, 1990

9:00-10:30 a.m.

ERB Building 5, Room 119

The ability to model an inert flow with density fluctuations allows one to model a

diffusion flame with fast chemistry, and gives insight into one aspect of the

dynamics of compressible flows. With this in view, these measurements were

made to calibrate a second order model for such flows. Measurements using

shuttle-mounted hotwire and Way-Libby probes are described. For comparison,

a jet of air into air was also measured. Moments up to fourth order were

computed. The data are used to evaluate various modeling assumptions.

Attempts are made to explain the greater spreading rate of the helium jet.

*From the Ph.D. thesis of N. R. Panchapakesan

Contact: Charles Feller, PABX 3-6681

30
ICOMP ^ CENTER FOR MODELING OF

\fs..m =^-^ ,^^ TURBULENCE AND TRANSITION

CMOTT SEMINAR SERIES

ENGINEERING TURBULENCE MODELING

Present and Future?

Tsan-Hsing Shih

Center for Modeling of Turbulence and Transition

Wednesday, July 11, 1990

3.00-4:00 PM

ERB Building 5, Room 119

A summary of the present position of eddy-viscosity models (e.g.

k — E) and second-order closure models (Reynolds stress models

is presented. Typical examples (comparisons between model predictions and experiments) show their abilities as well as their limitations.

Development of more advanced and complex schemes is discussed.

The inclusion of such models into a CFD commercial code is feasible,

but need intensive work - a cooperative effort!

Contact: T.-H. Shih, PABX 3-6680

31
^^Jti6 

` DR CphPG'r

N NASA o

ICOMP r CENTER FOR MODELING OF

\ ^ TURBULENCE AND TRANSITION

CMOTT SEMINAR SERIES

FUNDAMENTAL IN-HOUSE EXPERIMENTS

TO SUPPORT THE LEWIS TURBULENCE

MODELING PROGRAM

Edward J. Rice

Inlet, Duct and Nozzle Flow Physics Branch

Wednesday, July 25, 1990

3.00 PM

Building 5, Room 119

The fundamental experiments conducted by members of the Inlet, Duct and Nozzle

Flow Physics Branch of IFMD have been primarily for support of the shear flow control

effort within the Division. Although this support effort is expected to continue, some work

could be diverted (and in time expanded) to support the Turbulence Modeling Program

currently underway within ICOMP. This new emphasis can be accomplished through an

interactive effort between the Numerical Analysts and the Experimentalists. The current

and planned experiments will be summarized in this talk and include: unsteady flow around

airfoils (stationary and oscillating), 2D rapid diffusers (backward facing ramp), aspirated

backward facing step, circular and rectangular jets (subsonic and supersonic), dual stream

supersonic shear layer (annular geometry), and boundary layer transition. A swirl generator

within the plenum of the CW-17, ERB rig allows the addition of swirling flow to any of

the jet or shear layer experiments. The available and planned experimental instrumentation

include: single, X-wire and three-wire hot wire anemometry, single and two element corona

probe, multiple microphone and pressure transducer channels, Schlieren and laser-sheet flow

visualization, and conventional and fiber-optic two-component LDA systems. The new 16

channel anemometry system allows simultaneous measurement with 8 X-wires. Some sample

data will be presented to illustrate the current capability.

This is intended to be a group discussion type of workshop. The presentation will take

30 minutes and will be followed by 30 minutes of discussion between the numerical and

experimental participants. The presentation will be general in nature intended mainly to

acquaint the numerical analyst with the experimental capability which can support the

numerical program through a cooperative effort.

For additional information contact: T.-H. Shih, PABX 3-6680

32
F OR COtipG'^

NS i N ^^ .Q^ mnu OZ

oICOMP^ CENTER FOR MODELING OF

TURBULENCE AND TRANSITION

J1 )^^^bd N l S^ S

CMOTT SEMINAR SERIES

CFD-RELATED RESEARCH AT THE

UNIVERSITY OF BRUSSEL

AND TURBULENCE RESEARCH ACTIVITIES IN

EUROPE

Prof. Charles Hirsch

Vrije Universiteit Brussel

Wednesday, July 25, 1990

1.00-2.00 PM

Building 5, Room 119

For additional information contact: L. Povinelli, PABX 3-5818

33
1J,^ ^pR CphP\

NSA

y uw Matra mnai

z ICOMP CENTER FOR MODELING OF

TURBULENCE AND TRANSITION

CMOTT SEMINAR SERIES

CHAOTIC WANDERINGS THROUGH A LAND

OF TURBULENCE MODELS

Russ Claus

Aerothermochemistry Branch

Wednesday, August S, 1990

3.00 PM

Building 5, Room 119

A turbulence model user details his pursuit of the mythological "correct

answer". The talk opens with some failed attempts to calculate seperated

flows using a two equation turbulence model. Following this, additional

attempts to achieve the "correct answer" through a series of Direct Numerical Simulations and Large Eddy Simulations will be described and the

limitations of these approaches are highlighted. Finally, a re-examination

of two-equation and a second order closure calculations of a jet in crossflow

will be discussed.

A brief discussion of some turbulence modeling efforts being supported

under the NASA SBIR program will also be described. This includes RNG

modeling with Orszag and Yakhot, and PDF modeling for compressible

flows with Kollman and Farshchi.

For additional information contact: T.-H. Shih, PABX 3-6680

34
F OR CpMAG'j

ICOMPo CENTER FOR MODELING OF

^!^_ . TURBULENCE AND TRANSITION

Nl

CMOTT SEMINAR SERIES

A COMPARATIVE ANALYSIS OF TWO

EQUATION TURBULENCE MODELS

Nancy Lang

CFD Branch

Wednesday, August 22, 1990

3.00 PIN1

Building 5, Room 119

Several two equation models have been proposed and tested against

benchmark flows by various researchers. For each study, different numerical

methods or codes were used to obtain the results which were an improvement or success over some other model. However, these comparisons may

be overshadowed by the different numerical schemes used to obtain the results. With this in mind, several existing two equation turbulence models,

including k-c and k-T models, are implemented into a common flow solver

code for near wall turbulent flows. Calculations are carried out for low

Reynolds number, two dimensional, fully developed channel and boundary

laver flows. The accuracy of the different models is established by comparing the turbulent kinetic energy, mean velocity, and shear stress profiles

with the direct numerical simulations and experimental data.

For additional information contact: T.-H. Sbih, PABX 3-6650

35
( pa CC),AG

z ICOMP CENTER FOR MODELING OF

Oyd 

N1 

57^ TURBULENCE AND TRANSITION

CMOTT SEMINAR SERIES

ON LENGTH SCALE EQUATIONS IN

TWO-EQUATION TURBULENCE MODELS

Micha Wolfshtein

Faculty of Aerospace Engineering

Technion, Israel Institute of Technology

Haifa, Israel.

Wednesday, August 29, 1990

3.00 PM

Building 5, Room 119

Two-equation turbulence models show a great survivability between the very successful

mixing length models and the Reynolds stress models. Yet the models are often critisized,

mainly on the wide scatter of the computed results and on the difficulties encountered in the

derivation of a reliable scale equation. The problem has been difficult to resolve due to:

(i) The large resources required for developing solvers and to run test cases on computers;

(ii) The theoretical difficulties to derive turbulence models from the Navier Stokes equations.

The demands from a "good" turbulence model axe difficult to satisfy, and are often conflicting

with one another. This point will be illustrated in a discussion of some possible approaches

to this problem. In particular we shall refer to well established models like the dissipation

or length scale models, as well as newer models like the volume of turbulence or time scale

models. A generalised two-equation turbulence model will be used to demonstrate a possible

approach for the improvement of two-equation models. Finally, a fourth order boundary layer

solver for the generalized two equations model will be presented. The solver can handle both

compressible and incompressible flows, with any two-equation model, with or without wall

functions. Some results will be presented, for a flat plate boundary layer and for unseparated

diffuser flows.

For additional information contact: T.-H. Shih, PABX 3-6680

36
/` pR COG

NASA i

M Q^O^ mR0 =

COMP? CENTER FOR MODELING OF

TURBULENCE AND TRANSITION

CMOTT SEMINAR SERIES

EFFECT OF ACOUSTIC EXCITATION ON

STALLED FLOWS OVER AN AIRFOIL

Khairul Zaman

Inlet, Duct and Nozzle Flow Physics Branch

Wednesday, September 5, 1990

3.00 PM

Building 5, Room 119

Experimental results on the subject are to be summarized focussing attention on post-stalled flows, i.e. flows that are fully separated from near the

leading edge of the airfoil. The excitation results in a tendency towards

reattachement, which is accompanied by an improved airfoil performance,

although the flow may still remain fully seperated. It is observed that with

increasing excitation amplitude, the effect becomes more pronounced but

shifts to a Strouhal number which is much lower than that expected from

linear, inviscid instability of the seperated shear layer. In addition, some

results from a recent experiment on supersonic jets will be briefly reviewed

emphasizing need for collaboration between experiment and computations.

For additional information contact: T.-H. Shih, PABX 3-6630

37
"-"-(UR CpM^

,, IlVR QNO CaR1 Z0

=ICOMP? CENTER FOR MODELING OF .^..c w^wQn Sn TURBULENCE AND TRANSITION

Ni S

CMOTT SEMINAR SERIES

COMPUTATION OF ELLIPTIC FLOWS USING

LOW-REYNOLDS-NUMBER TWO-EQUATION

MODELS

Vittorio Michelassi

Energy Engineering Department

University of Florence, Italy

Wednesday, September 19, 1990

3.00 PM

Building 5, Room 119

Two new low-Reynolds number forms of the k — E model will be presented.

These exhibit better stability and stiffness characteristics as compared to

the previous formulations. Model generality is improved by formulating

the damping functions so that they do not depend on the wall distance.

The proposed formulations are compared with eight other low Reynolds

number two-equation turbulence models by computing the fully developed

channel flow and the incompressible flow past a hill. Results are compared

with available direct numerical simulation and experimental data. The flow

solver is based on the approximate factorization technique and the artificial

compressibility method requiring no (or very little) numerical damping. A

simple linearization technique for the turbulence model source terms based

on Taylor series expansion ensures implicit algorithm stability for all the

models tested. Both numerical accuracy and computational efficiency are

discussed.

For additional information contact: T.-H. Shih, PABX 3-6680

38
♦J^^ (OR 

CDMoG'1

N11W Q.ol mfl^ 0

z ICOMP r CENTER FOR MODELING OF

°. r. \s, TURBULENCE AND TRANSITION

CMOTT SEMINAR SERIES

Evaluation of Models for the Pressure Terms

in the Second Moment Equations

Aamir Shabbir

Wednesday, October 10, 1990

3.00 PM

Building 5, Room 119

The usual approach in establishing the correctness and accuracy of turbulence models is to numerically solve the modeled differential equations and

then compare the results with the experiment. However, in the case of a

discrepancy, this procedure does not pinpoint where the drawback of the

model lies. It is also possible that the model overcompensates one physical

phenomenon and undercompensates another resulting in good agreement

between the prediction and experiment. A more desirable approach is to

individually compare each term in the equations with its model. This talk

will focus on such a comparison for the pressure correlations appearing in

the second moment equations. At present these correlations can not be measured but can be obtained by balancing the turbulent flux and Reynolds

stress equations. Using this approach, pressure correlations will be directly

compared with their models, which range from simple linear to more elaborate nonlinear ones. Also results will be shown for a newly developed

model for the return to isotropy part of the pressure temperature-gradient

correlation. The data used is from the homogeneous shear flow and the

buoyant plume experiments as well as from the direct numerical simulation

of homogeneous shear flow.

For additional information contact: T.-H. Shih, PABX 3-6680

39
COMP CENTER FOR MODELING OF

TURBULENCE AND TRANSITION

CMOTT SEMINAR SERIES

HEAT TRANSFER WITH INHOMOGENEOUS

FREE STREAM TURBULENCE

S.N.B. Murthy

School of Mechanical Engineering

Purdue University, West Lafayette, IN 47907

Thursday, October 25, 1990

2:00 PM

Building 5, Room 119

The presence of free stream turbulence (FST) in a wall-bounded flow with

heat transfer presents several interaction complexities. These depend upon

(a) the state of the boundary layer, laminar, fully turbulent or transitional,

(b) the nature of FST, especially its "peakiness" or inhomogeneity, and (c)

other complications such as pressure gradient, geometry and cooling. An

investigation is in progress on the possible application of (a) large eddy interaction hypothesis (based on Lumley's rational description of turbulence)

and (b) spectral analogy between heat and turbulence kinetic energy to the

simplest case of a flat wall, fully developed turbulent boundary layer with

zero pressure gradient when there is heat transfer in the presence of homogeneous FST. An extension of the approach to the case of other types of

boundary layers with inhomogeneous FST will be discussed.

For additional information contact: T.-H. Shih, PABX 3-66SO

40
L ICOMP^ CENTER. FOR MODELING OF

\^,nsaK ^n ry^' TURBULENCE AND TRANSITION

CMOTT SEMINAR SERIES

PRESSURE-STRAIN MODELLING AND

THE RETURN TO ISOTROPY

M.M. Gibson

Imperial College of Science Technology,

London, England.

Wednesday, November 14, 1990

1.00 PM

Building 5, Room 119

Calculations of complex shear layers give the best results when the constant in Rotta's return-to-isotropy model is given values greater than those

deduced from the measurements in grid turbulence. New data from grid

turbulence are presented and assessed in the light of previous experiments.

The results show that homogeneous turbulence decay is associated with

increasing Reynolds number. It is argued that this finding has important

implications for Reynolds-stress Modelling. Recent theoratical and experimental studies of the analogous model for pressure scrambling in the scalarflux equations reveal even larger discrepancies in them "Monin constant".

The implications for the future of second-moment modelling are discussed.

For additional information contact: T.-H. Shih, PABX 3-6680

41
INSTITUTE FOR COMPUTATIONAL N uH v.ci mnl Z

zICOMP

MECHANICS IN PROPULSION

NI

1COMP SEMINAR SERIES

CURRENT TRENDS IN TURBULENCE MODELLING

by

W. Rodi

University of Karlsruhe

Institute for Hydromechanics

Thursday, November 29, 1990

1:00-2:00 p.m.

ERB Building 5, Room 119

A brief review is given of recent work in the area of modelling turbulence in near-wall regions and by

Reynolds-stress-equation models. Various low-Reynolds-number versions of the k-s model and their

damping functions are examined with the aid of results from direct numerical simulations and they are

compared with respect to their performance in calculating boundary layers under adverse and favorable

pressure gradients. A two-layer model is presented in which near-wall regions are resolved with a oneequation model and the core region with the standard k-e model. Various applications of this model

are shown. The ability of the various models to simulate laminar-turbulent transition in boundary layers

is discussed. Recent applications of a basic Reynolds-stress-equation model to complex flows of

practical interest are presented. Finally, some recent proposals for improved Reynolds-stress-equation

models are outlined and an outlook on possible future turbulence model developments is given.

Contact: Charles Feller, PABX 3-6681

42
G

J '^ tiNASA o

^, lfV9 op[MUI mrtOi Z

zICOMP^ CENTER FOR MODELING OF

TURBULENCE AND TRANSITION

COMPUTATION OF TRANSITIONAL AND TURBULENT

FLOWS IN COMPLEX GEOMETRIES

George Em Karniadakis

Mechanical and Aerospace Engineering,

Program in Applied and Computational Mathematics,

Princeton university, Princeton, NJ.

Friday, March 22, 1991

2:00 - 3:00 PNI

Building 5, Room 119

General, high-order numerical schemes are formulated appropriate for simulation of transitional and turbulent incompressible flows in complex geometries. In particular, the new

schemes are based on mined explicit/implicit stiffly stable time-stepping methods and explicit treatment of the pressure boundary condition that lead to enhanced stabilit y and

arbitrarily high-order time-accuracy dictated entirel y by the employed integration rule. H y

-brid spectral element methods are then used for the spatial discretization of the variable

properties governing equations in three-space dimensions. Special Neumann/viscous sponge

type boundary conditions are developed for open flows. Large or sub-rid scales in the high

Reynolds number regime are modeled through renormalization group ( R\G) techniques.

Simulations are than performed to stud y the transitional and full y turbulent stages of

spatially developing flows. Here, we consider the flow over a backward- facing step, and the

three-dimensional flo%v past a circular c ylinder. For the first flow, the secondary instability is

first investigated and the three-dimensional transitional states are computed through direct

simulation (DNS). Tra n sport and sub-rid RNG models are then employed to simulate the

high Reynolds number flow and heat transfer. Comparisons with experimental data in both

regimes are presented. For the second, three-dimensional equilibria are accurately resolved

via DNS; the series of bifurcations followed is simulated until the cylinder wake becomes

turbulent. Our results <_ug est a successive period doubling in the temporal response of the

flow, which eventuall y leads to a disordered state and renders the wake turbulent.

For additional information contact: T.-H. Shill, PABN 3-66SO

43
FOR COME+Gr'l

ti (AS O taw Q^ mn^ Y

` ICOMP

CENTER FOR MODELING OF

TURBULENCE AND TRANSITION

A DYNAMIC MODEL FOR LARGE EDDY SIMULATION

OF TURBULENT FLOWS

Parviz Moin

Stanford University

and NASA Ames Research Center

Tuesday, March 26, 1991

10:00 - 11:00 AM

Room 176, Sverdrup Building

One major drawback of the subgrid scale models currently used in large

eddy simulations is their inability to correctly represent with a single universal constant different turbulent flows in rotating or sheared flows, near

solid walls or in transitional regimes. A new subgrid scale model will be

presented which alleviates many of these drawbacks. The model coefficient

is inputted dynamically rather than inputted a priori. The basic idea in

the derivation of the model is the utilization of the spectral information

which is computed directly. This rather rich spectral information is not

available in methods based on Reynolds averaged equations. The subgrid

scale stresses obtained using the proposed model vanish in laminar flow and

at a solid boundary. The results of large eddy simulation of transitional and

turbulent channel flow that use the proposed model are in good agreement

with the direct numerical simulation data. The same model was applied

to the decay of isotropic turbulence with excellent agreement with the experimental data. The model has been extended to compressible flows. A

dynamic subgrid scale turbulent Prandtl number was derived. Its dependence on molecular Prandtl number, direction of scalar gradient, and the

distance from the wall are in accordance with direct simulation.

For additional information contact: T.-H. Shih, PABX 3-6680

44
J1(, 

f DR COhf•L 

N IIVY O[..o

ICOMP^

., CENTER FOR MODELING OF

'°d, d N S, rs TURBULENCE AND TRANSITION

CMOTT SEMINAR SERIES

TURBULENCE AND DETERMINISTIC CHAOS

Robert G. Deissler

Lewis Research Academy

NASA Lewis Research Center

Wednesday, April 10, 1991

3:00 - 4:00 PM

Room 119, Building 5

Several turbulent and nonturbulent solutions of the Navier-Stokes equations are obtained.

The unaveraged equations are used numerically in conjunction with tools and concepts from

nonlinear dynamics, including time series, phase portraits, Poincare sections, largest Liapunov exponents, power spectra, and strange attractors.

Initially neighboring solutions for a low-Reynolds-number fully developed turbulence are

compared. The turbulence, which is fully resolved, is sustained by a nonrandom timeindependent external force. The solutions, on the average, separate exponentially with time,

having a positive Liapunov exponent. Thus, the turbulence is characterized as chaotic.

In a search for solutions which contrast with the turbulent ones, the Reynolds number (or

strength of the forcing) is reduced. Several qualitatively different flows are noted. These are,

respectively, fully chaotic, complex periodic, weakly chaotic, simple periodic, and fixed-point.

Of these we classify only the fully chaotic flows as turbulent_ Those flows have both a positive

Liapunov exponent and Poincare sections without pattern. By contrast, the weakly chaotic

flows, although having positive Liapunov exponents, have some pattern in their Poincare

sections. The fixed-point and periodic flows are nonturbulent, since turbulence, as generally

understood, is both time-dependent and aperiodic.

For additional information contact: TAI. Shih, PABX 3-6680

45
J``'

_ NASA 0

z ICOMP=

CENTER FOR MODELING OF

TURBULENCE AND TRANSITION

CMOTT SEMINAR SERIES

A CRITICAL ASSESSMENT OF TURBULENCE MODELING

Charles G. Speziale

ICASE, NASA Langley Research Center

Monday, April 22, 1991

1:30 - 2:30 PM

Room 175, Sverdrup Building

A variety of turbulence models including zero, one and two-equation models as well as second-order closures will be reviewed. Based on comparisons

with results from physical and numerical experiments on homogeneous turbulence, a strong case will be made for the superior predictive capabilities

of second-order closure models. It will be shown how some significant improvements in second order closure models have been recently achieved

by means of invariance arguments coupled with a dynamical systems approach. Several applications will be discussed including recent extensions

to high speed compressible flows that were developed in connection with

the National Aerospace Plane Project.

For additional information contact: T.-H. Shih, PABX 3-6680

46
J i

^~ N/SA o

IC 0 MP CENTER FOR MODELING OF

TURBULENCE AND TRANSITION

CMOTT SEMINAR SERIES

A CRITICAL-LAYER THEORY

FOR BOUNDARY-LAYER TRANSITION

Reda R. Mankbadi

Lewis Research Academy

NASA Lewis Research Center

Tuesday, April 30, 1991

10:00 - 11:00 AM

Room 119, Building 5

An asymptotic critical-layer theory is developed to study nonlinear interactions of a triad

of instability waves leading to boundary-layer transition. This triad consists of a spatially

growing plane fundamental wave and a pair of symmetrical, subharmonic oblique waves. The

spatial development of the waves is determined by nonlinear viscous flow in the critical layer.

The theory successfully captures not only the linear and parametric resonance stages, but also

the later fully interactive regime of the transition process including the saturation and decay

stages. The analysis is fully nonlinear, in that all the nonlinearly generated waves that were

not originally present, are accounted for in the analysis. The three-dimensional nonlinear

modifications of the mean flow are also considered. The theory applies to both ribboninduced and naturally occurring transition. The analytically obtained amplitude equations

are highly accurate but simple enough to be used for practical transition predictions. Results

presented explain experimental observations and reveal novel features of the phenomena.

For additional information contact: T.-H. Shih, PABX 3-6680

47
Center for Modeling of Turbulence and Transition

Research Briefs - 1990

13.2 CMOTT Technical Meetings

These are CMOTT's biweekly group meetings. In each meeting, one

speaker from either CMOTT or other local modeling efforts at the NASA

LeRc presents new ideas or work under progress. The informal seminars were

intended not noly to keep the members informed of the latest development

of local turbulence and transition modeling research but also to increase

interactions between group members and other researchers at the NASA

LeRc.

The schedule of these biweekly events were widely disseminated before

each series began. The following is the meeting schedule during the reporting

period.

48
(OR CO.".

I- _ w^sn l

.o.,.a m..a o

[COMP CENTER FOR MODELING OF

TURBULENCE AND TRANSITION

'^ebd 

NI S^^?

Date: December 19, 1990

To: CMOTT Members and SVR and IFMD Staff

From: William W. Lion (6682)

Subject: CMOTT Biweekly Meeting

The following is a tentative schedule for CMOTT biweekly get-together from Dec.

19, 1990 to Feb. 27, 1991. In each meeting, one member of our group will give an informal

presentation on the subject she/he is interested in or currently working on. This will not

only give us opportunities to interact with each other, but also keep us informed about

exciting approaches in modeling turbulence other than our own specialities. However, the

format of the meeting after the current period is open to suggestions.

The meeting will start at 4:00 p.m. in Room 228, Sverdrup Building.

Dec. 19, 1990 Zhigang Yang (61-2123)

Understanding RNG in Turbulence: A Preliminary

Report

Jan. 2, 1991 Tsan-Hsing Shih (6680)

Realizability Concept and Its Applications in

Turbulence Modeling

Jan. 16, 1991 Wai-Ming To (5937)

Spectral Methods

Jan. 30, 1991 Le Tran (61-6701)

Stagnation Point Turbulence and Heat Transfer

Feb. 13, 1991 William W. Liou (6682)

Weakly Nonlinear Models for Turbulent Free Shear Flows

Feb. 27, 1991 Chris Steffen (8508)

Optimizing Laminar Accuracy

49
J1^ 

` OR COM,oGl

~, .a..a m a C5

X ICOMP= CENTER FOR MODELING OF

TURBULENCE AND TRANSITION

Date: March 13, 1991

To: CMOTT Members and SVR and IFMD Staff

From: William W. Liou (6682)

Subject: CMOTT Biweekly Meeting

The following is a tentative schedule for the CMOTT biweekly get-together from

March 20, 1991 to May 1, 1991.

This round of biweekly meetings will be more of presentations of progress and hurdles

in pursuing the subjects than presentations of results. Therefore, active participation is

expected from the attendants. Note that these meetings are different from the CMOTT

Seminar Series, whcih axe mainly formal presentations.

We would also appreciate some contributions from you. Subjects related to either the

theoretical, experimental or computational aspects of turbulence and transition modeling

are welcomed. Those who are willing to share their experience in these areas can contact

me or Dr. T.-H. Shih at 6680 for further arrangement.

The meeting will start at 4:00 p.m. in Room 228, Sverdrup Building.

Mar. 20, 1991 K. Zaman (5888)

Low Frequency Fluctuations in Separated Flows

- Experimental Evidence

Apr. 3, 1991 K. Kirtley (61-6659)

Let's Talk RNG

Apr. 17, 1991 T.-H. Shih (6680)

Lumley's New Formulation of Dissipation Equation

May 1, 1991 K. Kao (5965)

Stability of Compressible Couette Flows

50
l

cw o«a mnn O

COMP' CENTER FOR MODELING OF

TURBULENCE AND TRANSITION

Date: April 25, 1991

To: CMOTT Members and SVR and IFMD Staff

Froin: William W. Liou (6682)

Subject: CMOTT Biweekly Meeting

The following is a tentative schedule for the CMOTT biweekly get-together from

May 1, 1991 to June 12, 1991. Please note that Dr. Shih's and Dr. Kao's talks have been

rescheduled.

The presentations will be informal and active participation is expected from the attendants. These meetings complement the CMOTT Seminar Series, which are mainly formal

presentations.

We would also appreciate some contributions from you. Subjects related to either the

theoretical, experimental or computational aspects of turbulence and transition modeling

are welcomed. Those who are willing to share their experience in these areas can contact

me or Dr. T.-H. Shih at 6680 for further arrangement.

The meeting will start at 4:00 p.m. in Room 228, Sverdrup Building.

May 1, 1991 T.-H. Shill (6680)

Lumley's New Formulation of Dissipation Equation

May 15, 1991 A. Shabbir (5927)

On Some Closure Assumptions Regarding Time Scales,

Local Equilibrium and Pressure Diffusion

May 29, 1991 Z. Yang (61-2123)

Near Wall k — e Modeling: Another Crusade

June 12, 1991 K. Kao (5965)

Stability of Compressible Couette Flows

51
Center for Modeling of Turbulence and Transition

Research Briefs - 1990

Appendix C

Collection of Publications

52
THE STUDY OF PDF TURBULENCE MODELS

IN COMBUSTION

Andrew T. Hsu

Sverdrup Technology, Inc.

& Center for Turbulence Modeling

NASA Lewis Research Center

Cleveland, Ohio

9th National Aero-Space Plane

Technology Symposium

November 1-2, 1990

Paper Number 107

53
I. Introduction

(U) In combustion computations, it is well known that the predictions

of chemical reaction rates (the source terms in the species conservation equations) are poor if conventional turbulence models are used. The main difficult y lies in the fact that the reaction rate is highly non-linear, and the

use of averaged temperature, pressure and density produces excessively large

errors. Moment closure models for the source terms have attained only limited success. The probability density function (pdf) method seems to be

the only alternative at the present time that uses local instantaneous val-.

ues of the temperature. density, etc., in predicting chemical reaction rates.

and thus is the only viable approach for more accurate turbulent combustion

calculations.

(U) The fact that the pdf equation has a ver y large dimensionalit y renders finite difference schemes extremel y demanding on computer memories

and thus impractical. if not entirel y impossible. A logical alternative is the

Monte Carlo scheme. which has been used extensively in statistical physics.

The evolution equations for the joint pdf of the velocity and species mass

fraction have been successfully solved using Monte Carlo schemes, see, e.g.,

Pope[1]. However, since CFD has reached a certain degree of maturity as

well as acceptance, it seems, at least from the stand point of practical applications. that the use of a combined CFD and Monte Carlo scheme is more

beneficial. Therefore, in the present study a scheme is chosen that uses a

conventional CFD flow solver in calculating the flowfield properties such as

velocity. pressure. etc.. .chile the chemical reaction part is solved using a

Monte Carlo scheme.

(U) A combined CFD-Monte Carlo computer algorithm has been developed recentl y. As a first calibration of the Monte Carlo solver recently

developed in this Nvork. the discharge of a heated turbulent plane jet into

quiescent air is studied. Experimental data for this problem shows that when

the temperature difference between the jet and the surrounding air is small,

buovancv effects can be neglected and the temperature can be treated as a

passive scalar. The fact that jet flows have a self-sillular solution lends convenience in the modeling studv. Furthermore, the existence of experimental

data for turbulent shear stress and temperature variance (temperature fluctuation) make the case ideal for the testing of pdf models wherein these

54
values can be directl y evaluated.

(U) The following section presents the methodologies used in this study,

and a discussion on the numerical results and their comparison with experimental data follows.

II. Method

2.1 Governing equations

(U) The Favre averaged momentum and species transport equations can

be written as

pdtui + pujdJ u i = — diP + UdJTi) — dju+^uJ

pd1 Y; + pic j dj Y = p

(9,( DOJI — I.„u”) + c);(Dd;pOjY,") + pfvi

where ui is the velocity and 3" is the crass fraction. The major problem encountered in the numerical study of turbulent combustion lies in the modeling

of the chemical source term. pwi , known as the chemistry closure problem.

The Source term wi is an exponential function of temperature, and it is well

know that the use of averaged temperature, T, in evaluating pwi can cause

egregious errors. Therefore the accurate prediction of turbulent combustion

using conventional turbulence models becomes ver y difficult. This motivated

the use of pdf methods. If the pdf, P. is given, then the mean of the source

term can be evaluated exactly:

p wi = J ... I pwi(Yi, .... Y,^, T, p ) P( Yc, - - -, Yn, T, p ) dY1 ...dY'l,dTdp

To solve for the pdf, we need the following:

2.2 Evolution equation for the pdf

(U) The evolution equation for the probability density function of the

mass fractions, temperature, and density, P(Y,, ..., Y,, T, p), can be written

as(lj

;v

p91P + p b'a'p + p E dV, { wiN)l, ---, ON)P}

i-1

N N a

_ 0"(p

 C tea I lr^i P) — p

V,G.+G (< E iJ I 'Nk > P)

55
where the terms represent mean convection, chemical reactions. turbulent

convection. and molecular mixing, respectively; P is the density-weighted

joint pdf:

P=PP,

e is the scalar dissipation:

E i; _ —PDi;c;l'„

Vi 's are variables such as the mass fractions, density, and temperature, and

<, > denotes the mathematical expectation of a function.

(U) The left hand side of the above equation can be evaluated exactly

and requires no modeling; the right hand side terms contain the conditional

average of the Favre velocity fluctuation and the conditional average of the

scalar dissipation and require modeling.

2.3 Finite difference solution for the NS equations

(U) In order to simplify the problem so as to concentrate on the stud y of

the pdf models, we have confined our numerical procedures to parabolic flows.

A parabolic NS solver with a k — e turbulence model has been developed.

The k — f model has been tested by solving standard cases such as flat plate

boundary lavers and free shear lavers with satisfactor y numerical results.

2.4 Monte Carlo scheme for the pdf equation

(U) A grid dependent Monte Carlo scheme has been employed, primarily

following Chen and Kolimann[5[. The task here is to construct an ensemble

of sample points, each sample has its own distinct properties such as temperature. mass fraction. etc.: these properties change with time or location

such that the probabilit y function of the ensemble evolves according to the

evolution equation. Consequently, the pdf of' the ensemble is the desired

approximate solution for the pdf equation.

(U) Suppose N samples are assigned to a grid cell in the flow domain,

the pdf for the ensemble of N samples can then be written as

N

where o, is the scalar function value carried by the j` h sample, e.g., the mass

fractions. etc.

56
(U) In order to have the pdf P' evolve according to the pdf evolution

equation. «v e discretize the equation on a given grid and write, for parabolic

flows:

P'r+dx.i = ai P'x.j+i + 3jP*=,j

 + 'iP•=.i -1

and require

aj +;3; +7; = 1

2.5 Recontamination in a Monte Carlo scheme

(U) Using a very simple test case of a convection/ diffusion process with

two scalars. it was found that the previous scheme(5] does not conserve mass

fractions due to re-contamination. It is found that in order to conserve the

mass fractions absolutel y, one needs to add further restriction to the scheme,

namely

A new computer algorithm was devised and tested. This algorithm uses a

few extra arrays to store informations form the previous time level and thus

eliminates the repetition in the sampling process. With the simple test case

of two scalars with assumed constant coefficients in the pdf equation, the

new algorithm is shown to conserve the mass fractions perfectly. Deficiencies

such as directional bias and re-contamination that were found in the previous

algorithm are completely eliminated.

(U) It is not yet known whether one can indeed devise a scheme that

satisfies relation (*) in general flows, where the coefficients of the pdf equation

are variable and are calculated from the flow velocities, turbulence time scale,

etc.

57
III. Results

(U) The sketch of a heated plane jet is given on this page. The flow

domain is divided into 30 cells in the cross direction. and 100 samples are

assigned to each cell. The numerical results from the present study are

compared with experimental data in the following figures.

P _ , — I

Figure 1. Sketch of a heated turbulent free jet.

58
(U) The numerical solution for the ensemble averaged pdf's of the temperature: y/y0.5 = 0. where yo.5 is the location where T = 0.5To. In order to

interpret the data, the ensemble averaged pdfs are integrated to obtain the

approximate continuous pdf: P = P`(OT)/AT. The result of this integration is shown as vertical bars. The normal distribution is plotted as a dotted

line for comparison.

Fig. 2 probability density distribution

0.8 F

Present results 0.7 — — — Normal distrib.

0.6

a 0.5

0.4

a

0.3

0.2

0.1

0.0 ^I III !"I 1 i I 1 11111111

NIT1111I p

-3 -2 -1 0 1 2 3

(T- < T > )/(T

59
(U) The numerical solution for the ensemble averaged pdf's of the temperature: y/yo.s = 0.44, where yO _ 5 is the location where T = MTo. In order

to interpret the data, the ensemble averaged pdfs are integrated to obtain

the approximate continuous pdf: P = P"(OT)IAT. The result of this integration is shown as vertical bars. The normal distribution is plotted as a

dotted line for comparison.

Fig.3 probability density distribution

0.8

Present results 0. 7 -- — Normal distrib.

0.6

a 0.5

\ 0.4 , -

a

0.3 '

0.2

0.1 ^J

0.0 I!^ ^I I I IIl n I I I I !11111 I III ^^ -

-3 —2 —1 0 1 2 3

(T—<T>)/g

60
(U) The computed velocity profile is compared with experimental data

from Foerthmann(7]. The good agreement between the numerical solution

and experiments indicates that the parabolic flow solver emploved in the

present study is fairl y accurate and can serve its purpose.

Fig.4 Velocity distribution

1.0

0.8

0 0.6

0.4

0.2

0 0

0 Present

o Experimentsl7]

0

0

8 0

0

0

0

0 p

0 00

0

8

0.0 0.5 1.0 1.5 2.0 2.5

y\y. 5 U

61
(U) As a calibration for the k — e turbulence model. the turbulent

shear stress from the present stud y is compared with measurements given

b y Brandbury[S], and by Gutmark and NVygnanski[9]. Except for a slight

over prediction towards the edge of the jet, the numerical solution generally

agrees well with the experimental data. This comparison is very crucial for

pdf calculations because it is an indicator of the accuracy of the numerically

predicted turbulent kinetic energy, k. and dissipation, e, which are used to

determine turbulent time scale, t, from t = k1 e. An accurate prediction for

k and e is essential for the correct modeling of molecular mixing in the pdf

calculation.

0.030

0.025

0.020

0.015

0.010

0.005

Fig.5 Turbulent shear stress

Present

O Exp. Ref. [8]

o O° o

c Exp. Ref. [9]

O O

o O

0

0

0

00

0

O 0

0.000

0.0 0.5 1.0 1.5 2.0 2.5

Y/Y. 5Uo

62
(U) The results for mean temperature from the \Monte Carlo solution of

the pdf equation for the temperature compared with experimental data from

Ref. 10-12.

1.0

0.8

0 0.6

0.4

0.2

00

Fig.6 Mean temperature

r

G

•

C? q Present

Q7 o Ref. [10]

Q^q n o Ref. [11]

aoe o Ref. [12]

q . 0

gob

^o 

I,- O

d°

.^

00

^ 8

.p

° of

0.0 0.5 1.0 1.5 2.0 2.5

y/y.5To

63
• . • . •

•, 0 o (D D e

° O ° ° OOo o •

.e o° p m

O

0 O

O

O

O

O

• Present

O Ref. [121

o Ref. [13]

e

e

°

O

O

O

(U) The standard variance of the temperature, or the temperature fluctuation, compared with experimental data from Ref. 12 and 13. Because of the

limited number of sample points (N=100) used in the computation, the numerical results are scattered: however, these results agree very well with the

experimental observations. Particularly important is the agreement between

the predicted and the measured temperature fluctuations, for in case of finite

rate chemistry calculations, the instantaneous temperature will be used in

determining the reaction rate, and the correct prediction of the temperature

distribution will ensure the accuracy in reaction rate calculations.

Fig.7 -Standard variance of temperature 0.40

0.35

0

0.30

0.25

°

A

0.20

0.15

V

0.10

0.05

0 00

0.0 0.5 1.0 1.5 2.0 2.5

y/y. 5 To

64

L
Acknowledgement

The author wish to thank Prof. S.B. Pope of Cornell University and Dr.

J.-Y. Chen of Sandia National Labs for their help. This work is supported

by the NASA Lewis Research Center under Contract NAS3-25266.

References

1. Pope, S.P., Prog. Energy Coinbst. Sci. 1985. Vol-11. pp.119-192

2. Kollmann, W., Theoret. Comput. Fluid Dynamics, 90 1. pp. 249-285

3. Clien, J.-Y. and Kolhnann, NV. Combust. Sci. and Tech. 1989, Vol.

64, pp. 315-346

4. Shuen, J.-S. and Yoon, S. AIAA J. Vol. 27 No. 12, 1989, pp. 1752-1760

5. Chen, J.-Y. and Kollmann, W. Combust. and Flame. 1990, pp.75-99

6. Brabbs, T.A. and Musiak, J.D., Private Communications

7. Schlichting, H., "Boundary Layer Theory." McGraw-Hill, 1989.

8. Brandbury, L.J.S., JFM, Vol. 23, 1965, pp. 31-64.

9. Gutmark, E. and Wygnanski, I., JFNI, Vol. 73, 1976, pp. 465-495.

10. Browne, L..I.S., Antonia., R..A., and Chambers, , A.J., JFM, Vol. 149,

1984, pp. 355-373.

11. Uberoi, M.S. and Singli, P.I., Physics of Fluids, Vol. 18, No.17, 1.975,

pp 764-769.

12. Jenkins, P.E. ASME J. of Eng. for Power., Vol. 98, 1976, pp.501-505.

13. Antonia, R.A., et al., Int. J. of Heat Mass Transfer, Vol. 26, No. 1,

1983, pp. 41-48.

65
On Recontamination and Directional-Bias Problems

in Monte Carlo Simulation of PDF Turbulence Models

Andrew T. Hsu

NASA Lewis Research Center, Cleveland, Ohio 44135

Phone 216/826-6648

Turbulent combustion can not be simulated adequately by conventional moment closure turbulence models.

The difficulty lies in the fact that the reaction rate is in general an exponential function of the temperature, and

the higher order correlations in the conventional moment closure models of the chemical source term can not be

neglected, making the applications of such models impractical. The probability density function (pdf) method

offers an attractive alternative: in a pdf model, the chemical source terms are closed and do not require additional

models.

The partial differential equation for the probability density function, P, can be written as

N

Pat e +Pb,L%p +P t9v,; {w;(01'....tGN)P)

;=t

N N

_ -aQ P < 0;^;(< (+jI VGk > P)

=t;=t

where the terms represent the time derivative, mean convection, chemical reaction, turbulent convection, and

molecular mixing, respectively. The fact that the pdf equation has a very large dimensionality renders finite

difference schemes extremely demanding on computer memory and CPU time and thus impractical, if not entirely

impossible. A logical alternative is the Monte Carlo scheme, wherein the number of computer operations increases

only linearly with the increase of number of independent variables, as compared to the exponential increase in a

conventional finite difference scheme.

A grid dependent Monte Carlo scheme following that of J.Y. Chen and W. Kollmann has been studied in the

present work. In dealing with the convection and diffusion of the pdf, the pdf equation is discretized on a given

grid, e.g.,

Ps+dti=a^P=i+t 1QiP=i+7iPzi—t

where

aj + Q, + 7j = 1

However, if this is the-only restriction satisfied by the numerical algorithm, the mass fractions may not be conserved due to re-contamination, and directional-bias also appears. These phenomena are illustrated in Figure 1:

Consider a mixing layer; use white balls to represent contaminants in the upper stream and black balls to represent

contaminants in the lower stream. As the two streams move toward right, the location of the white balls and black

balls are interchanged randomly to simulate convection and diffusion. From Figure 1, it is clear that directionalbias caused by recontamination caused the center of the mixing layer to drift downward. ( Directional-bias can

be partially corrected by changing sweeping directions. ) One also notices that after the first marching step, the

conservation law is violated, reflected in the Figure as missing white or black balls.

It is found that in order to conserve the mass fractions absolutely, one needs to add further restriction to the

scheme, namely

aj + 7i = aj— t + 7i+1

A new algorithm was devised that satisfies this restriction in the case of pure diffusion or uniform flow problems.

Using the same example, it is shown that absolute conservation can be achieved. This result is shown in Figure 2.

One can see than the diffusion process is symmetric, and the problem of directional- bias is eliminated.

Although for non-uniform flows absolute conservation seems impossible, the present scheme has reduced the

error considerably.

66
01.1

RECTION

W

ECTION

F;e_ure 1. Illustration of directional-bias and recontariination

problems in Monte Carlo simulation.

Fiqure 2. Solution obtai^ea by using additional co:;ltrai-,t

in the solution procedure.

67
PROGRESS IN THE DEVELOPMENT OF

PDF TURBULENCE MODELS FOR COMBUSTION

Andrew T. Hsu

Sverdrup Technology, Inc.

NASA Lewis Research Center

Cleveland, Ohio

10th National Aero-Space Plane

Technology Symposium

April 23-26, 1991

Paper Number 213

68
I. Introduction

(U) A combined Monte Carlo—CFD algorithm has been developed recently at NASA Lewis Research Center for turbulent reacting flows. In this algorithm, conventional CFD schemes are employed to obtain the velocity field

and other velocity related turbulent quantities, and a Monte Carlo scheme

is used to solve the evolution equation for the probability density function

(pdf) of species mass fraction and temperature.

(U) In combustion computations, the predictions of chemical reaction

rates (the source terms in the species conservation equations) are poor if

conventional turbulence models are used. The main difficulty lies in the fact

that the reaction rate is highly non-linear, and the use of averaged temperature produces excessively large errors. Moment closure models for the source

terms have attained only limited success. The probability density function

(pdf) method seems to be the only alternative at the present time that uses

local instantaneous values of the temperature, density, etc., in predicting

chemical reaction rates, and thus may be the only viable approach for more

accurate turbulent combustion calculations. The closure problem and the

need for pdf has been discussed in detail in a previous paper[1].

(U) Assumed pdf are useful in simple problems; however, for more general combustion problems, the solution of a evolution equation for the pdf is

necessary.

H. Method

(U) Conventional CFD flow solvers are quite common and needed no

explanation here. We will concentrate on the derivation and solution of the

pdf evolution equation in this section.

2.1 PDF evolution equation of mass fraction.

(U) For simplicity, the pdf equation for a single scalar will be derived

here. The extension to multi species is trivial. The species transport equations ( for a single species ) can be written as

patY = —puk akY + pak(DakY) + pw (0.1)

where uk is the velocity, Y is the mass fraction, and w is the chemical source

term; summation for repeated indices is understood.

69
Let P(Y) be the probability density function, and define 0 = exp(%Ac),

then the characteristic function of the pdf can be written as the ensemble

average of 0:

< ¢ >= f P(Y)OdY,

and the pdf can be written as the inverse Fourier transform of the characteristic function:

P(Y) 2r J < 0 > exp(—iaY)dY.

Now differentiate 0 with respect to time,

ato = i.lexp(1*AY)(atY),

and replace atY in the above equation with the right hand side of eq. (1),

we obtain

PatO = iAe_xp (%aY) [— PU kakY + ak(pDakY) + Pw],

which can be rearranged as

at(PO) + ak (Puk^) = iA Oak(PDakY) + taopw

Take ensemble average of the above equation gives

at < PO > +ak < puk O >= t*A < pwo > +ia < Oak(PDakY) >

The inverse Fourier transform of the above equation gives the evolution equation for the probability density function:

at(PP) + ak(P < uk > P) + 0Y(PwP)

_ —ak(P < uklY > P) — ay(< ak (pDakY)lY > P)

where the terms represent mean convection, chemical reactions, turbulent

convection, and molecular mixing, respectively.

For multi species and temperature, the pdf is written as P(0 1 ,..., 0„),

where tki , i = I,—, n represent the scalars including mass fractions and temperature. A similar derivation gives

N

Pate + PvQaaP + P a^; {w;( 1 , ..., 4 N)P}

N N

_ —a

"

(P < vale; > P) — P^1:a^:^;(< C:jIOk > P).

70
(U) The left hand side of the above equation can be evaluated exactly

and requires no modeling; the right hand side terms contain the conditional

average of the Favre velocity fluctuation and the conditional average of the

scalar dissipation and require modeling.

2.2 Modeled pdf equation

(U) The first term that needs modeling is the turbulent diffusion term,

for which the following gradient model is used:

— < uk10i >P"'Dta^P,

where Dt is the turbulent diffusion coefficient, which is set to be equal to the

eddy viscosity, i.e., the turbulent Schmidt number is equal to one.

(U) The next term needs modeling is the molecular mixing term. A

coalescence/ dispersion model is used for this term, which has the following

general form:

N N

—P a2,,t; (< E;, I Ok > P)

CD f P(0')P(;b1')T (iP j Of , 011 ) dO ld^" — P

J = D,n P

where T is the transition probability. For more detail description of the

molecular mixing model, see Ref 2. The modeled the pdf equation is then

N

AP + PijAP + P aOi { w+(01, ..., ON)P}

= a,(Dt aaP) + D,nP

2.3 Solution of the pdf equation

(U) A fractional step method is used to solve the pdf equation in a

Monte Carlo simulation. Three processes, namely, (1) convection/ diffusion,

(2) molecular mixing, and (3) chemical reaction, are simulated consecutively:

N

(at + p a^:{w,(^1, ..., ^N)})( at — D,n)(at + pvaa — aa (Dt aa )P = o.

71
Present

O E.P. [3]

J.0

0.8

0 0.6

0.4

0.2

00

This equation is solved using a Monte Carlo simulation, i.e., the pdf is

represented by an ensemble of samples. These sample points move in the

hyperspace of (X 1, x2, 23, 01, ..., ;b„) according to the above equation.

III. Results

(U) A heated turbulent plane jet had been calculated and results reported at last NASP symposium. The case has since been recalculated using

modified turbulence model, and improved results are shown here.

Velocity distribution

1.0

0.8

0.6

0.4

0.2

00

0.0 0.5 1.0 1.5 2.0 2.5

y/y 5U°

0.030

Present

0.025 < 0 E.P. [a].[s]

o•

0.020 i

e °

0.015 -

o

0.010

0.005 -

0

0.000 y

0.0 0.5 1.0 1.5 2.0 2.5

y/y 5Uo

Mean temperature

Present

°oae E.P. [6,1,3]

^e0

0

13

0

°ae

0.0 0.5 1.0 1.5 2.0 2.5

y/y 5To

Standard tianance of temperature

U. <r v

0.35

0 0.30

0.25

In

Iz 0.20

ca 0.15

v 0.10

0.05

0 00

0.0 0.5 1.0 1.5 2.0 2.5

y/y 5 T°

Present

o A o E.P. [8,9]

0 0 •o

- c

o < o

_

- 0 +

°

0

°

0

72
U 

961 N^+ 4% F2

rature

u 2 y or n2 r 91 n2

Velocity distribution

—0.10 0.00 0.10

y/x

7

Turbulent kinetic energy

6

5

4

y

3

2

1

0

—0.10 0.00 0.10

y/s

1.0

0.8 i

N

0.6

j 0.4

0.2

0.0

(U) A hydrogen—fluorine diffusion flame was calculated. The flowfield

and flow conditions, as well as the calculated velocity and turbulent kinetic

energy distribution, are given on this page.

73
(U) Figures show the calculated heat release of the H2-F2 flame (temperature distribution) compared to experimental data, and the species mass

fraction from pdf calculation.

Fig.1 Temverature distribution.: H2+F2=2HF

1.0

0.8

a

0.6

a 0.4

0.05

0.04

0.03

0

0.02

Z

0.01

0 Exp. Hermanson

Dimotakis

?resent PDF solution

0

O

0

0

0

—0.4 0.0 0.4

yob T

Fig.2 Mass fraction: H2+F2=21-fF

-F2----

UHF.

/ A

0.2

00

—0.8

0 00

—0.4 0.0 0.4

y/6T

74
p -0. 15

-0.42

(U) Computed pdf of temperature distribution in the H2-F2 flame.

Fig-2 probability mass function.

0.3(

---- Y/6t=-0.42

Y/bt=-0.15

0.25

0.20

a

0. 15

a

,

,

0.0 0.2 0.4 0.6 0.8

AT/ATQd f

0.10

0.05

0 00

1.0

75
Acknowledgement

This work is supported by the NASA Lewis Research Center under Contract NAS3-25266.

References

1. Hsu, A.T. "A Numerical Study of Combusting Shear Layers Using PDF

Turbulence Model," to be presented at the 8th Symposium on Turbulence Shear Flows; Munich, Germany, Sept. 9-11, 1991.

2. Hsu, A.T. "The Study of PDF Turbulence Models in Combustion," 9th

NASP Symposium, Orlando, Florida, Nov. 1-2, 1990.

3. Schlichting, H., "Boundary Layer Theory." McGraw-Hill, 1989.

4. Brandbury, L.J.S., JFM, Vol. 23, 1965, pp. 31-64.

5. Gutmark, E. and Wygnanski, I., JFM, Vol. 73, 1976, pp. 465-495.

6. Browne, L.J.S., Antonia, RA., and Chambers, , A.J., JFM, Vol. 149,

1984, pp. 355-373.

7. Uberoi, M.S. and Singh, P.I., Physics of Fluids, Vol. 18, No.17, 1975,

pp 764-769.

8. Jenkins, P.E. ASME J. of Eng. for Power., Vol. 98, 1976, pp.501-505.

9. Antonia, R.A., et al., Int. J. of Heat Mass Transfer, Vol. 26, No. 1,

1983, pp. 41-48.

10. Herma.nson, J.C. and Dimotadis, P.E., JFM. (1989), vol. 199. pp.333-

375.

76
A STUDY OF HYDROGEN DIFFUSION FLAMES USING

PDF TURBULENCE MODEL

Andrew T. 11su'

Sverdrup Technology, Inc.

NASA Lewis Research Center, Cleveland, Ohio 44145

Abstract

The application of probability density function (pdf) turbulence models is addressed in this work. For

the purpose of accurate prediction of turbulent combustion, an algorithm that combines a conventional CFD

flow solver with the Monte Carlo simulation of the pdf evolution equation has been developed. The algorithm

has been validated using experimental data for a heated turbulent plane jet. The study of H2—F2 diffusion

flames has been carried out using this algorithm. Numerical results compared favorably with experimental

data. The computations show that the flame center shifts as the equivalence ratio changes, and that for the

same equivalence ratio, similarity solutions for flames exist.

Introduction

It is the consensus of the combustion profession

that the prediction of chemical reaction rate (the

source term in a species conservation equation) is

poor if a conventional turbulence model is used. The

main difficulty lies in the fact that the reaction rate

is highly non-linear, and the use of averaged temperature, pressure and density produces excessively

large errors. Moment closure models for the source

terms have attained only. limited success because the

assumptions for such models to be valid often can

not be satisfied. The probability density function

(pdf) method seems to be the only alternative at the

present time that uses local instantaneous values of

the temperature, density, etc., in predicting the chemical reaction rate. and thus is the only viable approach

for more accurate turbulent combustion calculations.

Two main lines are being followed in pdf methods:

one uses an assumed shape for the pdf, the second

solves a pdf evolution equation; the present paper addresses only the latter. There has been significant

progress in the stud y of pdf turbulence models in low

speed flows in the past decades. These developments

were summarized in Refs. [1], [2] and [3]. In spite

of this progress. pdf turbulence modeling remains a

nascent discipline with many unresolved issues.

The fact that the pdf equation has a very large

dimensionality renders the use of finite difference

schemes extremely demanding on computer memories

and thus impractical, if not entirely impossible. A logical alternative is the \Monte Carlo scheme, which has

been used extensively in statistical physics. The evolution equations for the joint pdf of the velocity and

species mass fraction have been successfully solved

using Monte Carlo schemes, see, e.g., Pope[l]. Ilowever, since CFD has reached a certain degree of maturity as well as acceptance, it seems, at least from

the standpoint of practical applications, that the use

of a combined CFD and Monte Carlo scheme is more

beneficial. Therefore, in the present study a scheme

is chosen that uses a conventional CFD algorithm to

solve the Navier-Stokes equations and provide flowfield properties such as velocity, pressure, etc_, and the

chemical reactions are calculated by using a Monte

Carlo scheme to solve a pdf evolution equation.

The combined CFD-pdf solver has been developed recently and validated using non-reacting flow

data for a heated turbulent plane jet. The algorithm

has been further tested in the numerical study of an

H2 —F2 diffusion flame. This diffusion flame was studied experimentally by \fungal & Dimotakis[4] and

Hermanson & Dimotakis[5]. The numerical results

from the present study are compared with these experimental data.

Theory

Governing emiations for reacting flows.

Flows with chemical reaction are governed by

the continuity equation, momentum equations, energy equation. and species transport equations:

atP+ (di pui = 0

Pat u; + puj aj ui = -ai p + M ai Tij

Pat h+Puj aj — at p — uj aj p (1)

q, + Q

Patl k + Puj Oj 1"k = pa, (Daj Yk ) + wk

k = 1,2 ......V.

where u i is the velocity, }k is the mass fraction, and

wk is the chemical source term. (In addition to these

equations, one also needs the equation of state.) For

turbulent flows, we substitute

ui = ui + ui+

(2) Yi=Y"+Y,

into the above equation and take an ensemble average. In the process of averaging, new unknown quantities in the form of correlations appear, e.g., uiuj',

u'Yk, etc. These quantities can all be modeled using

conventional turbulence models, such as two equation

models or second order closure models. A problem

unique to flows with chemical reaction is the average

77
of the chemical source term, pwk , The major difficulty lies in the fact that w k is a highly nonlinear

(usually exponential) function of the temperature. It

is well know that the use of averaged temperature, T,

in evaluating pwi can cause egregious errors. One may

consider the effect of the temperature fluctuation, T',

by applying a moment closure model to the term pwi;

however, such a closure model results in an infinite series that converges only when T. — T and T' «T. In

many combustion, problems, these two conditions are

violated: in fact, we often have T. » T and T' ^- T

instead. In view of the above, the prospect of an accurate prediction of turbulent combustion using conventional turbulence models seems dismal. The above

fact motivated the use of pdf methods.

Pdf evolution equation

Given a set of m random variables, 0 1 , 02 , • • • , ,pm , and the joint probability density function, P(01 , 02, • • -, 0",; T,, y, Z, t), the mean of

any random function, P01, 02 , • • • Om-1 x, y, z, t), can

be calculated as

f(x, y ,=,f) = f ... f f Pd0i d02 ... dV)-- (3)

For simple ;"lows, one could assume the shape of

the pdf and compute the source term, 77w-7k , based on

the assumed pdf using the above integral. For more

general problems, one needs to solve a pdf evolution

equation for a more accurate pdf distribution. In the

latter case, the evaluation of the source term, pwk , is

no longer necessary because the mean values of the

temperature, species mass fractions, which we previously would obtain by solving the transport equations

(Eq. 1), can now be evaluated directly from the pdf

using the above integral.

The pdf evolution equation can be derived from

the transport equations (Eq. 1) in many different

ways, using a Dirac delta function, a characteristic

function. or a characteristic functional. [1-3] The evolution equation of a single point probability density

function of scalar random variables w 1 , r[,,, can be

written as

m

PatP+pi5'aaP+P {wi(V i,...,ZGm)P}

i.l

= —aa(P < va1¢k( x ) = Ok > P)

—P av Y; (< Eij 10k ( y) _k > P) (4)

i=1 j=1

where the terms represent the rate of time change,

mean convection, chemical reaction, turbulent convection, and molecular mixing, respectively; P is the

density-weighted joint pdf:

P = PP/P,

ei; is the scalar dissipation:

E ij = Daa0iaa-Oj,

(where D is the diffusion coefficient), and < xly >

denotes the mathematical expectation of a random

function x conditioned upon y.

The left hand side of eq. (4) can be evaluated

exactly and requires no modeling; the right hand side

terms contain the conditional expectation of the velocity fluctuation and the conditional expectation of

the scalar dissipation, which are new unknowns and

require modeling.

Closure models for pdf equation

The first term that needs modeling is the turbulent convection term, for which one can use the

following gradient model:

— < v" alV k > P = D=BaP,

where Di is the turbulent diffusion coefficient, which

is set to be equal to the eddy viscosity, i.e., the turbulent Schmidt number is equal to one.

The next term needs modeling is the molecular

mixing term. A coalescence/dispersion model is used

for this term, which has the following general form[6]:

N N

i=1 j=1

= Co f f 

F,( ,)P( „

)7 (t^^P^,;b")dO'd,p" — P

M (P)

where T is the transition probability. A new mixing model continuous in time is recently developed by

Hsu and Chen[7]. For more detail description, of the

molecular mixing models, see Refs. 6 and 7.

The modeled pdf equation is then

m

pa'P+paaaaP+P a0,1wiN1,... (Wpl

i=1

=aj (Do%P) + Al (P), (5)

which can be solved using a Monte Carlo simulation.

Numerical Methods

In the present study, a combined finite

difference—Monte Carlo solver is developed. For the

velocity field, the Navier-Stokes equations and a k — c

turbulence model are solved using a finite difference

method; for scalar variables such as the temperature,

mass fractions, etc., the pdf evolution equation is

solved using a Monte Carlo scheme. The CFD flow

solver provides the Monte Carlo solver with the mean

velocity and a turbulence time scale r, where r = k/e,

and the Monte Carlo solver provides the mean flow

solver with the density, P.

The solution'of the Navier-Stokes equations

78
Since the primary concern of the present work is

pdf modeling, we choose the simplest possible solver

for the N-S equations. :lie equation is transformed

into a general coordinate system using the following

general relation:

a=f

Y = W, n)

and

ar = ac — Uf all

Yv

1

a — Yn an

The transformed momentum equation for steady

flows is

1 1

Ofu

u

+(v — 6 Yf) 1 (9,

= —Pafp+ 1

YU 

a" ( VT av u ) , Y?

where the viscous term involves gradient in the ( -

direction is neglected since Nve are only interested in

shear flows in the present study.

First order upwind difference is used in the

direction so that a marchine scheme for steady

parabolic flows can be used. To ensure stability, a

flux splitting scheme is used in the r7 -direction.

The solution of the Pdf equation

A fractional step Monte Carlo method, as will

be described in the following, is used to solve the pdf

equation. We first discretize the time derivative in eq.

(5) using finite difference: (9,P = ( P1nt — P")/Of

then the pdf evolution equation (eq. 5) can be written

as

m

Pn +1 _ {1— ,tUQatr—^t ^CjL',w,(t3 1,....t^Im)

;-1

+Ot aQ D, aQ + P"

Using approximate factorization, the above equation

is recast as

m

n+1 •/, ))

;.1

x(1+A- 0f)(1 —Atpva, +AtatrD,,% P"

+O(At)

(1 + AtQ1 + AtM)(1 + AtR)P" +O(Ot),

where C denotes the convection operator, m the

molecular mixing, and R chemical reactions. With

the above expression, we can three processes consecutively:

(1) convection:

(2) molecular mixing:

P" - (1 + At.1f) P' ,

(3) chemical reaction:

Pn+1 = (I+At R)P"

In the present study, only steady flows are considered.

so at is replaced by of and a marching scheme in the

^-direction is employed.

In a Monte Carlo simulation, the continuous pdf

is replaced by N x M delta functions,

P'(w1,02," ,1Jm.S,Y,Z,t)

1N

= N b(^1 — oinl(t))

n=1

X 6 (02 — 9(n)(t)) ... b (O.n — 0(n)(t)),

where each product of the m delta functions represents one event of an ensemble of N sample events.

An event can be thought of as a fluid particle. and

the evolution of P' entails the movement of the particles in the physical space as well as the phase space

(V)—space). The movement of the particles is, of

course. governed by the pdf evolution equation. and

is simulated in the following three steps.

Step 1: convection

Replace P by P' in the pdf evolution equation

and transform it into the l —Y coordinate system, the

convection process for a steady flow can be written as

PUaf P' = — P( v — uYf)yna,P. + I av(PDs(9.,P*),

Discretizing the above equation using a finite difference scheme, we can write

P;.i = aP;_1,] +1 + 3P;_l.i +7P, -l.i-=

(Here again we used a one sided difference so that a

marching scheme could be used.) The above equation

states that if we divide the flow field into cells. then

the pdf at point (i, j) can be written as a linear combination of the pdf at neighboring points, To simulate

this process with a Monte Carlo scheme, we move the

sample particles between cells according to the above

equation. For instance. the particles of cell (i, j) will

be obtained by choosing randomly ON particles form

cell (i — 1, j+ 1), pN particles from cell (i— 1, j), and

7tV particles from cell (i — 1, j — 1). In order for the

total number of particles not to change, we require

that a+i3+-r = 1.

It is worth noting that this method can easily

be extended to elliptic flows and applied to general

curvilinear coordinate systems.

Step 2: molecular mixing;

79
The molecular mixing process is simulated by the

following binary interaction model:

puBE P

=CD f f P( ,P )P(0")T(0JO', ^p")d;b'dV," — P

where T is the transition probability. By assigning

various functions to T, we would have different mixing models. In the present study, the modified Curl

model by Janicka et al. [8] is used. In this model, the

transition probability is given as

for t/i' <;b<t/"

0 otherwise.

with CD = 6.0.

In the Monte Carlo simulation, the above model

is realized in the following manner: Divide the flow

domain into small cells, each containing N sample

particles_ Given a small time interval At and a turbulent time scale r, select randomly N,,,r pairs of particles, where

Nmr = O.J^,T.^'.

and let a pair, say, m and n. mix as follows

6„(t + At) = AOm(t) + (1 — A)0,(t)

0m(t + At) = Ao„(t) + (I — .-1)0m(t)

where A = 0.51;, with l a random variable uniformly

distributed on the interval [0.1]. The remaining N —

2N,,, = particles remain unchanged:

0n(t + At) = an(t)

The turbulent time scale is supplied by the finite

difference solution of a k — c turbulence model: r =

k/e.

This model admits the non-physical jump condition and does not produce the correct long time

behavior for decay problems in homogeneous turbulence. A continuous model that predicts the correct

long time behavior for turbulence decay problem has

been introduced by Hsu and Chen [7]. But as shown

in Ref. [i], for practical combustion problems where

the long time statistical behavior is not crucial, the

modified Curl model gives acceptable results. In the

present study, both the modified Curl model and the

continuous model have been used.

Step 3: chemical reaction

Chemical reaction is represented in the Monte

Carlo simulation by the movement of sample particles

in the phase space due to reaction. A sample point

with a given composition {0 1 ,0 2 ,•• , 0 m 1 at time t

will acquire new composition at the next moment, and

the changing rate is

I dd, = wi(01,02,"

i = 1,2,- ,M,

where wi is the reaction rate for that specific sample

point. If we regard wi as the convection velocit y of

the particle in the phase space, the it is clear that

the above ODE represent the movement of a sample

particle in the phase space.

The various chemistry models considered in this

study are presented in the following section.

Chemistry Models for H, — F, Reaction

Finite rate chemistry_

The reaction of H2 -+ F2 = 2HF can be represented by the following two step chain reactions[4]:

H2 +F--'' HF+H

H+F, HF+F

with

610

k t = 2.6 x 10127.0.5exp

RT )

k 2 — 3 x 10 9 Tl sexp

RT

1680 )

where T is in K, R in cal cool -1 K-1 , and k in cm-3

mol' sec'.

The above expressions only provide the forward

reaction rate; the backward reaction rate can be calculated using the e q uilibrium coefficient:

k b = kflkeq,

and the equilibrium coefficient can be evaluated using

Gibbs free energy. Polynomials for Gibbs free energy

for various species can be found in Ref. [8]. A numerical experiment shows that the backward reaction

rates are much smaller than that of the forward reaction for the two reaction steps described above and

can be neglected.

Let C denote the mole concentration; for each

sample particle, we need to solve the following set of

ODES:

d

=— k i CH, CF = w t

dt

t

d CF2 = —k, CH CF, = w2

dt

d CH

dt — —w t + w 2 = ws

d CF

di --L —W2 =W4

d CH F

dt

_

— 

-wl — w2 = ws

d 

_ —

s

h j, Wi

dt

i-1

where h is the enthalpy and h' is the heat of formation.

80
Using n to denote the time level. we found that

the following semi-implicit scheme seems to be most

stable:

a few hundred degrees Kelvin. Under these circumstances, the variation of specific heat C, is negligible,

and the temperature can be calculated as

C" +1 H' = CHI

1+Atk1CF

Cn+1

F'

n

= F,

1 + At kzCH

C" +1

H

— C"

+ At k l C',CF

1 + At kzCH

Ci}1

F

_ CF + ^t k- CF, H

1+At k1CH,

CHF1 = CHF + At (k 1 CK^ CF + k, C;1 CF= )

The energy equation is solved explicitly. To start the

chain reaction 0.00002 mole fraction of F was released

into the flowfield initially.

Using the above equations to calculate the reaction in a larninar premixed flame, we found that the

time required for this reaction is the order of one microsecond (Figs. 1 and 2). The turbulence time scale

in the calculation. r = k/e, is the order of one second. Considering this large difference of time scales.

it is impractical to use finite rate chemistry in this

calculation.

Fast reaction

Since the chemical reaction is very fast compared

to the flowfield development and the reverse reactions

are negligible, we chose to use the following complete

irreversible reaction mechanism. We assume that after molecular mixing, complete reaction is achieved

within each particle during the time interval of one

marching step in the calculation. The reaction is calculated using the following equations.

For Cy7 > C"

C" + 1 = C

H n — Cn 7

H2 F7

CF+ 1 =0

HF1 = ^HF + 2CF,

For CHI < CFA :

CH71 = 0

C F+ 1 - CFA - cn

CHFI = CH F + 2CH,

The energy equation can be written as

Ah = —h f, ACHF (Mp F) ,

where MHF if the molecular weight of HF.

Since the concentrations of reactants are low, the

heat release is low, and the temperature rise is within

7'n+1 =7"+

Ah

C;'

The error involved in this approximation is less than

5%.

Results and Discussions

Code validation.

Before applied to the hydrogen-fluorine diffusion

flames, the computer algorithm is first validated using experimental data for a heated turbulent jet. The

non-reacting flow data serves as a check for the convection and molecular mixing process in the pdf solver

as well as the k — c model in the N.-S. solver.

Extensive experimental results for turbulent

plane jet have been reported by many authors. Measurements for mean velocity field in a turbulent jet

were first reported in the 1930's[9]; turbulent shear

stress measurements had been reported more recentiv[10,11]. To determine the effect of turbulence

on mixing, the temperature field of a heated turbulent jet had been studied by several authors. The turbulent jet has a slightly higher temperature than the

ambient. Measurements of both the the mean temperature and the rms of the temperature fluctuations

were given[12-17].

In the present study of this non-reacting flow

case, the temperature field is treated as a conserved

scalar and is simulated by the pdf of the temperature; the velocity field and turbulent shear stress are

obtained by solving the N.-S. equation and a k — e

turbulence model.

In the finite difference solution of the flowfield,

41 grid points are used across half of the jet width. A

symmetry boundary condition is used a-. the jet centerline. Fig. 3 shows the comparison of the present

solution of the mean velocity field with the experimental data, and Fig. 4 presents the numerical solution

of the turbulent shear stress as compared to the experimental data. Good agreements between numerical results and experimental data are observed for

both the mean velocity field and the turbulent shear

stress. The turbulent shear stress is calculated from

— < uV >= v,au/8y, where v_ = Ck2 /e. A good

prediction of the turbulent shear stress ensures that

the turbulent time scale, r = k/c, supplied to the

Monte Carlo simulation is correct.

For the Monte Carlo simulation of the temperature field, two sample sizes of 1000 and 1500 sample particles per cell are used. The predicted mean

temperature and the root mean square (rms) of the

temperature variation are presented in Figs. 5 and

6. The results show that both calculations produce

fairly good comparisons with the experimental data,

which means a sample size of 1000 particles per cell

is large enough.'

81
The above validation lends credibility to the N.S.

solver as well as the !Monte Carlo solver developed in

the present study.

H, — F, diffusion flames

The flow conditions for the H, — F, diffusion

flames are set according to an experiment performed

b y Hermanson and Dimotakis (1989). The flame consists two streams. The upper stream contains [V_, and

F2 , the flow velocity is UI = 22 m/s: the lower stream

contains of IV-2 and of H2i with velocity U2 = 8.8 m/s.

In the present study, 6 cases involving various percentage of H2 and F2 in the upper and lower stream

are considered; the conditions are listed in Table 1,

where AT,,v is the difference between the adiabatic

flame temperature and the free stream temperature.

tax

No.

eq uivalence

ram

lower stream,

mole fraction

of H,

upper stream,

mole fraction

of F,

OT.,f

1 1 0.02 0.02 186

2 1 0.04 0.04 368

3 1 0.06 0.06 554

4 1!4 0.01 0.04 151

5 1/4 0.02 0.08 302

61 1/4 0.04 0.16 600

Table 1. Initial conditions of the flames calculated in the present study.

Fig. 71 shows the calculated temperature rises due

to combustion for cases 2 and 5 of Table 1 and the corresponding experimental data. In the figure. ST is the

shear laver thickness determined by 1% of the temperature rise. AT is the actual temperature rise due

to combustion. ( the two streams have the same temperatures initially,) and ATadf is the adiabatic flame

temperature assuming complete reaction. The solution for case 2 agrees fairly well with the experimental

data for the same case. while the mean temperature

is slightly over predicted for case 5. The computation

shows the flame center shifts as a result of change in

equivalence ratio, which is consistent with the experimental results.

The rms of the temperature variance for cases 2

and 5 are given in Fig. 8. The results show that in

a diffusion flame, the temperature variance has two

peaks, and the highest values do not coincide with

the maximum values of the temperature distribution.

The shifting of the flame due to the change of equivalence ratio can also be observed from this figure.

The mass fractions of H2 , F2 and HF for cases

2 and 5 are presented in FIgs. 9 and 10. respectively.

One can see that although a complete reaction chemistry model was used. with the pdf method, results

simular to that of a finite rate computation are produced, which is one of the many advantages of the

pdf method.

One experimentally established fact is that for

the same equivalence ratio, with various mole concentrations of fuel and oxidizer in the flowfield, the

normalized temperature stays the same [4,51. The

present calculation confirmed this. Figs. 11 and 12

are the mean temperatures and rms's of the temperature variance for an equivalence ratio of one; three

mole concentrations of fuel were considered. One can

see that the curves coincide with each other. The

same agreements are found for equivalence ratio 1/4

(Figs. 13 and 14).

The pdf distributions for flame temperature at

the center and at the outer edge of the flame (for case

2 of Table 1) are plotted in Fig. 15, where the xaxis denotes (T— < T >)/Q, with < T >= T being

the mean temperature and v the rms. The pdf distributions show that at the center of the flame, the

temperature with the highest probability is not far

from the adiabatic flame temperature, while at the

outer edge of the flame, most of the time one would

find a temperature close to that of the free stream

temperature. Nonetheless, as a result of external and

internal intermittency, low temperature fluid does exist at the center and high temperature fluid the outer

edge. Pdf distributions for various species concentration can also be obtained from the solution. but will

not be presented here.

Concludinff Remarks

A Monte Carlo soi'ution algorithm for the pdf

evolution equation has been developed and successfully combined with a finite difference flow solver in

the study of turbulent combustion. The algorithm

was validated using turbulent mixing data from nonreacting flows. Turbulent diffusion flames of H2 -

F2 were computed using the pdf method, and good

agreements between numerical solution and experimental data were observed. The computation identified the change of equivalence ratio as the cause of

flame shift, and demonstrated that similarit y solution exists for flames with the same equivaience ratio. The present work showed that a grid dependent Monte Carlo scheme can be readily applied to

a general curvilinear coordinate system: therefore. it

is suitable for realistic reacting flow computations.

Acknowledeement

The author would like to thank Prof. S.B. Pope,

Drs. J.-Y. Chen and J.S. Shuen for their helpful discussions. This work is supported by by NASA Lewis

Research Center under Contract NAS-25266 with Dr.

L.A. Povinelli as program monitor.

References

1. Pope, S.B.. " PDF Methods for Turbulent Reactive

Flows." Prg. Energy Combust. Sci., 1985. 11, 119-

192.

2. Kollmann. W.. "The PDF Approach to Turbulent

Flow." Theoret. Comput. Fluid Dynamics. 1990. 1,

249-285.

3. Kuznetsov, V.R.. and Sabel'nikov, V.A., Turbulence and Combustion, Hemisphere Publishing Corp.,

1990.

82
4. Mungal, i\f.G. and Dimotakis. P.E.. "Mixing and

Combustion with Low Heat Release in a Turbulent

Shear Layer" J. Fluid Mech. 1984, 148, 349-382.

5. Hermanson, J.C. and Dimotakis, P.E., "Effects of

Heat Release in a Turbulent, Reacting Shear Laver."

J. Fluid Mech. 1989, 199, 333-375.

6. Janicka, J., Kolbe, W., and Kollmann, W., "Closure of the Transport Equation for the Probability

Density Function Scalar Field." J Non-Equilib. Thermodyn. 1979, 4, 47.

7. Hsu, A.T. and Chen, J.Y., "A Continuous Mixing

Model for PDF Simulations and its Applications to

Combusting Shear Flows." 8th Symposium on Turbulent Shear Flows, Munich, Germany, Sept. 9-11,

1991.

8. Gordon, S. and McBride, B.J., "Computer Program for Calculation of Complex Chemical Equilibrium Compositions, Rocket Performance, Incident

and Reflected Shocks, and Chap man-Jouguet Detonations." NASA SP-273, March. 1976.

9. Forthmann. E., Uber Turbulente Strvhlausbreitung,

Diss. Gottingen 1933; Schlischting, H., p 747.

10. Bradbury, L.S., "The Structure of Self-Preserving

Turbulent Plane Jet," J. Fluid Mech., 1965, 23, 31-64.

11. Gutmark. E. and Wvgnanski. I., "The Plane Jet,"

J. Fluid Nlech.. 1976, 73, 465-495.

12. Antonia, R.A., Browne, L.W.B., Chambers, A.J..

and Rajagopalan. S., "Budget of the Temperature

Variance in a Turbulent Plane Jet." Int. J. Heat Mass

Transfer, 1983, 26-1, 41-48.

13. Ashir, J. and Uberoi, M.S., "Experiments on Turbulent Structure and Heat Transfer in a Two Dimensional Jet." Physics of FLuids, 1975, 18-4, 405-410.

14. Browne, L.W.B., Antonia, R.A., and Cliembers,

A.J., "The Interaction Region of a Turbulent Plane

Jet." J. Fluid Mech., 1984, 149, 355-373.

15. Jenkins. P.E. and Goldschmidt, V.W., "Mean

Temperature and Velocity in a Plane Turbulent. Jet,"

ASME J. Fluids Eng., 1973, 95, 581-584.

16. Uberoi, M.S. and Singh, P.L. "Turbulent Mixing

in a Two-Dimensional jet," Physics of Fluids. 1975.

18-7, 764-769.

17. Bashir, J. and Uberoi, M.S.. "Experiments on

Turbulent Structure and Heat Transfer in a Two Dimensional Jet." Physics of Fluids, 1975, 14.4, 405-410.

83
ODU

700

Y 600

01

v

500

4001-

300 / V//

0.0 0.5 1.0 1.5 2.0 2.5

time (10 -6 sec)

Figure 1. Temperature rise in a premixed laminar H2 -F2 flame.

0.030

0.025

-.0.020

,A 0.015

V 0.010

0.005

0.000

O

0.0 0.5 1.0 1.5 2.0 2.5

Y/Y.5Ua

Figure 4. Turbulent shear stress in a 2D turbulent jet. o Ref [101, A Ref [11], - Present calculation.

0.05

0.04

C

0

v 0.03

0

0.02

0

E

0.0 0.5 i.0 1.5 2.0 2.5

;ime (10 -6 sec)

Figure 2. Changes of mass fractions in a premixed laminar H7 -F? flame. - - - H2 , - - - F2,

HF.

1.0

0.8

0 0.6

\ O

0.4

0.2

0.0

0.0 0.5 1.0 1.5 2.0 2.5

Y/Y.5Ua

Figure 3. Mean velocity distribution in a 2D turbulent jet. o Ref (91, - Present calculation.

1.01-

0.8

o "

a 0.6

r0.4

0.2

00

0.0 0.5 1.0 1.5 2.0 2.5

Y/Y.5To

Figure 5. Mean temperature in heated plane jet.

A Ref [141, n Ref [171, o Ref [161 o Ref (151. - pdf

solution. 1500 particles per cell. - pdf solution.

1000 particles per cell.

0.40

0.35

0.30 d a°

0.25

n^ 00 no 000Q>0,

n n^ D o O

0.20 0 0 0 c0

V 0.15 00 \

0.10 0

o^

0.05 O

O

0.001 O

0.0 0.5 1.0 1.5 2.0 2.5

Y/Y 5Ta

Figure 6. RMS of temperature variance in heated

plane jet. o Ref. [121, o Ref [171, D Ref [161, - pdf

solution. 1500 particles per cell; - - - pdf solution.

1000 particles per cell.

0.01

0 00

84
1.0

0.8

00.6

0.4

A(n-)-\\

0.2

0.0 -0.8 -0.4 0.0 0.4

Y16-

Figure 7. Mean temperature rises in H2-F2

diffusion flames. Case 2 (refer to 'fable 1): o Ref [5],

- present; Case 5: 0 Ref [5], - - - present.

U.J

0.4

o -

a 0.3

a ^ ^

0.2

u r

0.1 /

0.0-0.8 -0.4 0.0 0.4

Y/dFigure 8. RMS of Temperature variance in H 2 -

F, diffusion flames. - equivaience ratio 1 (case 2,

Table 1). - - - equivalence ratio 1/4 (case 5).

U. UO

- 0.07

0.06

0.05

Ln 0.04 6

N 0.03

0.02

0.01

0.00-0.8 -0.4 0.0 0.4

Y16T

Figure 10. Mass fraction distributions in an H 2 -

F2 flame, case 5, equivalence ratio 1/4. - H2 , - - -

F2 , - - HF.

l.0

0.8

Ti 0.6 o

a

a 0.4

0.2

0.

-0.8 -0.4 0.0 0.4

Y/6T

Figure 11. Mean temperature rises in H2-F,

flames, equivalence ratio 1. Referring to Table 1: -

case 1, o case 2, 0 case 3.

U.UJ U.

0.04

:V 0.03

0.02 ti

^ 0.01 1 / \

0.00-0.8 -0.4 0 0.4

Y/6T

Figure 9. Mass fraction distributions in an H,-

F2 flame, case 2, equivalence ratio 1. - H2 , - - - F2,

- - HF.

Jr

0.

9

0

a 0..

a

0_.

E

0.

0.-0.8 -0.4 0.0 0.4

Y/6T

Figure 12. RMS of temperature variance in H2 -

F2 flames, equivalence ratio 1. Referring to Table 1:

- case 1, o case ' 2, 0 case 3.

85
^.3

^.5

a

F- -

-.2

2.0

-0.8 -0.4 0.0 0.4

y/6T

Figure 13. Mean temperature rises in H2-F2

flames, equivalence ratio 1/4. Referring to Table 1:

— case 4, o case 5, 0 case 6.

-.5

-.4

0

a

3.2

E

--0 

-0.8 -0.4 0.0 0.4

y/6T

Figure 14. RMS of temperature variance in H2 -

F2 flames, equivalence ratio 1/4. Referring to Table

1: — case 4, o case 5, D case 6.

:.7

^.5

1 \

1 1

I \

I \

Z_3 1

i \

C1 I ^

-3 -2 -1 0 t 2 3

(T-<T> )/o

Figure 15. Probability density function distributions in an H2 -F2 flame (case 2). — pdf at the center

of the flame, - - - pdf at the out edge of the flame.

86
An Improved k-e Model for Near Wall Turbulence

T.H. Shih* and A.T. Hsu**

Center for Modeling of Turbulence and Transition

NASA Lewis Research Center

Cleveland, Ohio 44135

Abstract

This paper presents an improved k-e model for low Reynolds number turbulence

near a wall. The work is twofold: In the first part, the near-wall asymptotic behavior of

the eddy viscosity and the pressure transport term in the turbulent kinetic energy equation

are analyzed. Based on these analyses, a modified eddy viscosity model with the correct

near-wall behavior is suggested, and a model for the pressure transport term in the kequation is proposed. In addition, a modeled dissipation rate equation is reformulated,

and a boundary condition far the dissipation rate is suggested. In the second part of the

work, one of the deficiencies of the existing k — e models, namely, the wall distance (e.g.,

Y + ) dependency of the equations and the damping functions, is examined. An improved

modal th^.t does not depend on any v,, dl distance is introduced. Fully developed turbulent

channel flows and turbulent boundary layers over a flat plate are studied as validations for

the proposed new models. Numerical results obtained from the present and other previous

k-e models are compared with data from direct numerical simulation. The results show

that the present k-e model, with added robustness, performs as well as or better than other

existing models in predicting the behavior of near-wall turbulence.

1. Introduction

The k-e model is one of the most widely used_ turbulence models in enginee.iug

applications. Patel et al. [ '] recently reviewed existing two-equation models that can be

integrated directly to the wall. One of their conclusions was that the damping functions

used in turbulence models, especially the one for the eddy viscosity, need to be further

modified in order to improve model performance. In fact, as we shall see later, many

existing k-e models do not provide the correct near-wall behavior of the eddy viscosity.

Shihl'l recently proposed a new near-wall k — e model based on asymptotic analysis.

The present paper is a direct extension of that work.

In the present paper, we will first analyze, in section 2, the near-wall asymptotic

behavior of the eddy viscosity and the pressure transport term in the k-equation, and

in sections 3 and 4, propose models according to their near-wall behaviors. The model

equation for the dissipation rate is reformulated following an argument similar to that of

Lumley,lsl and a boundary condition for e is suggested.

An asymptotic analysis shows that, in the near wall region, while the pressure

transport term in the turbulent kinetic energy equation is small compared to the dissipation

and molecular diffusion terms, it is much larger than the turbulent transport term, and

* Institute for Computational Mechanics in Propulsion.

** Supervisor, Comput. Phys. Section, Sverdrup Technology, Inc., Member AIAA.

This paper is declared a work of the U.S. Government and

is not subject to copyright protection in the United States.

87
it plays an important role in the balance between the dissipation and molecular diffusion

terms. This near-wall behavior is also observed in direct numerical simulation of fully

developed channel flows (Mansour et a1. 13 ^) , Kim et al. [ ' ] ). However, in existing k-e models,

this pressure transport term is either ignored or lumped into the turbulent, transport model.

The present work introduces a model for the pressure transport term explicitly.

Most of the existing k-e models for near-waL turbulence use y+ (defined as ury/v,

where u r the friction velocity) as a parameter in constructing damping functions, with the

Jones-Launder model being the only exception. While the use of y + is perfectly fine for

simple attached boundary layer flows, it is inconvenient in more general applications such

as separated flows and flows with corners, where y + is not well defined. The Jones-Launder

model has the advantage of avoiding y + ; however, the model is known to perform poorly

in predicting near wall turbulent quantities, especially the turbulence kinetic energy. In

the present work, a new damping function is derived based on asymptotic analysis (Section

5). The new function is constructed upon a non-dimensional quantity that is independent

of the coordinate system.

The new models proposed in this paper were validated using direct numerical simulation data for fully developed turbulent channel flows and turbulent boundary layers over

a flat plate. These numerical results are reported in Section 6. Comparisons are also made

with other popular k — e models implemented in the same computer code. The numerical

results show that the present model, in general, performs better than the existing models

while providing added robustness.

2. Asymptotic Analysis

To analyze the near-wall asymptotic behavior of the eddy viscosity and other turbulent quantities, we expand the fluctuating velocities and pressure in Taylor series about

the wall distance as follows:

u i — b i y + c l y2 + dl y3 + ...

u 2 = C2 y2 + d2y3 + ...

U 3 = b3 y + C3y2 + d3 y3 + ...

(1)

p = ap + bpy + cpy2 + dpy3 + ...

where the coefficients ap , bl , C2 , ... are functions of x, z and t. Using the continuity and

momentum equations, Mansour et al J31 derived the following relations between the coefficients,

2C2 = — ( b i,i + b3,3)

ap,1 = 2vc 1

(2)

ap , 3 = 2vc3

The eddy viscosity is usually defined as

2

— ( u = ui) = vT(U` >i + Ui, i ) — 2kbii (3)

88
where ( ) stands for ellscinble average and k = (v rdi) /2 is the turbulent kinetic cnrrgy.

For plane shear flows, we can write from Eq. (3)

VT = aulay

()

and using Eq. (1), we obtain the near-wall asymptotic behavior of the eddy viscosity:

VT _ —( -

/

1

2) y 3 + — ( bl d2 + C1 

2

2 ( b , C2)( C 1) y4

+ 0(y 5) (rJ)

1

That is, near the wall vT is 0(y3 ). A correct eddy viscosity model should have this nearwall behavior. We shall see later that many existing models do not have this near-wall

behavior. For later use, let us examine also the near-wall asymptotic behavior of the

turbulent kinetic energy k and its dissipation rate E = v(u i,j u i j ). Using Eq. (1), we

obtain the following relations for the k and e:

( b ]) (b3 )yz 2

+((b,Cl) + ( b 3 l:3))y3 +OW) l6) 2

E — (b 1

2)

+ ( b3) + 4((b1c1) + ( b3 C3)) y + O ( y2 ) ( 7)

V

In addition, the pressure transport term in the k-equation, 11 = — v (u ip,i ), becomes (using

Eq.s (1) and (2))

II = —2v((blcl) + ( b3C3 ) ) y + O( y2 ) ( S)

The turbulent transport term in the k- equation, —(ku i) i , can be estimated as 0(y3).

Therefore, the pressure transport term is much larger than the turbulent transport term

near the wall.

3. Eddy viscosity model

In this section, we will propose a model for the eddy viscosity using its near-wall

behavior described in the previous section. In general, the eddy viscosity model can be

written as

11T = C UT (9)

where u' and t' are the turbulent characteristic velocity and length scale, respectively.

Depending on how the velocity and length scales are specified, the edd y viscosity model

can be a mixing length model, a one-equation (k) model or a tyro-equation (e.g. k-e)

model. For example, in plane shear flows, Prandtl's mixing length model specifies the

characteristic velocity with 69Ulay. For near wall turbulence, the Van Driest mixing

length model further damps the length scale to y[1 — exp(—y'/A)]. For more advanced

mixing length models, see Baldwin and Lomax [ '] , and King[ '] . One-equation (k) models

use k 1l2 as the characteristic velocity, which is determined by the turbulent kinetic energy

89
equation. In two-equation models, e.g. k-E models, the length scale is usually specified

by 012 1E, where E is determined by a dissipation rate equation. In this paper we will

concentrate on two-equation models, wherein the eddy viscosity is usually written as

k2

VT = Cμ f,, k (10)

E

where C,, = 0.09, and fμ is a damping function. The form of the damping function

is critical in such formulations, since the prediction of the mean velocity field depends

primarily on the eddy viscosity model. Thus it is important for an eddy viscosity model

to have the proper near-wall behavior. We have examined the near-wall behavior of eddy

viscosity models based on various k -E model equations. The results are listed in Table 1,

which shows that some of the k-E models do not have the correct near-wall behavior of the

eddy viscosity, namely, vt = 0(y3).

The quantity k3/ ' /E is usually considered a characteristic length scale, P', (or the

size) of the energy containing eddies. One expects that near the wall the size of these

eddies to be of the order of the distance from the wall, O(y). However, Eq.s (6) and (7)

show that k3 ^ 2 /E is O(y 3 . Hence, 0 12 /E is not an appropriate quantity to represent the

length scale of the large eddies near the wall. We therefore introduce a new variable E:

ak/ax i ak/axi

E =E -v

2k

which has the following property: E approaches E away from the wall and is 0(y 2) near the

wall. Therefore, k. 3/2 /E is a proper quantity to characterize the length scale of the large

eddies. With this length scale, the eddy viscosity should be written as

k2

LIT = C",

fμ 

e

(12)

Now in order for vT to have the correct near-wall behavior, the damping function fμ

must be O(y) near the wall and approaches 1 away from the wall. The damping functions

used in various k-E models are listed in Table 2. If we consider the presence of the gall

as the main effect on the eddy viscosity, then we may assume fμ is mainly a function of

Y+ . The form of fu can be determined quite accurately if we know vT , k and E from, for

example, the direct numerical simulation. One may optimize the following simple form by

numerical experiments:

f, = 1 - exp(-a l y+ - a2 y+

2

- a3 y+ 3

 - a4y+

4

 ) (13)

The optimal values for channel flows are a 1 = 6 x 10 -3 , a2 = 4 x 10-4 , a3 = -2.5 x

10 -6 , a4 = 4 x 10 -9 . It can be shown that this form of damping function does provide

the required near-wall behavior. As will be shown later, the above constants are valid for

general boundary layer flows.

(11)

90
4. Modeled k-e equation

To complete the eddy viscosity model, we need the modeled equations for the turbulent kinetic energy and its dissipation rate. In this section we will analyze the near-wall

behavior of the k-equation and propose a model for the pressure transport term Nvith a

proper near-wall behavior. The equation for the dissipation rate is also reformulated with

a formal invariant analvsis.

4.1 k-equation

We start with the equation for the turbulent kinetic energy,

k,:+U;k,;=D„+T+H+P— e (14)

where D,,, T and H represent the transport of the turbulent kinetic energy due to the

viscosity, turbulent velocity and pressure, respectively. P and e are the rate of production

and dissipation of the turbulent kinetic energy. The terms on the right hand side of Eq.

(14) are defined as follows:

D„ = vk,ii

T = —(kuj),i

P _ —(utuj)U;,j

E = v(u=,iu=,i)

tT Sing Eq.s (1) and (2), we obtain the budget of the k-equation near the wall,

Dk = O(y3)

Dt

D, = "(( bi) + (b3)) + 6v ( ( b i ci ) + ( b3 c3)) y + O(y2)

T = 0(y3) (16)

H = —2v((b l c l ) + (b3C3 ))y+O(y2)

P = O(y3)

E = v (( b2 ) + ( b3)) + 4v(( b l cl) + ( b3 c3)) y + O(y2)

This budget shows that the term 11 is much larger than the term T, and that H cannot

be neglected if we want the k-equation be balanced in the near-wall region. However, the

existing models either do not consider this term or simply combine it with the term T and

model them as

VT

( 17)

ak l'i

In this paper, we propose a model for 11 which has a similar form to that of the standard

turbulent transport model, but with a coefficient to ensure its correct near -wall behavior,

Eq. (8). The proposed model form of II is

Co VT— k,, ( 18 )

II = fμ [1 — exp(—y+)] Qk ^'j

91
where Co = 0.05 is a model constant. In some existing k -E models, it is assumed that E = 0

at the wall. In that case, in order to balance the term D,,, a nonzero artificial term D

must be added to the k-equation. The form of D for various k -E models is listed in Table

3. Finally, the modeled k-equation becomes

ak ak_ a vT ak avi aU)

avi

at +U^ax; ax^^ (v+ a k ) ax;^ +II+vT( ax; ax,)axj - E +D (19)

In the present model, D = 0, since E is nonzero at the wall. (The boundary condition for

E will be discussed later.)

4.2 E-equation.

The exact dissipation rate equation is

E ,t + U; E , ; = D' + T E + II` + PD (20)

where D', T` and 11` represent the diffusion rate of the dissipation rate due to the viscosity,

turbulent velocity and pressure, respectively, and PD stands for the entire mechanism of

the production and destruction of the dissipation rate E. The terms on the right hand

side of the above equation are as follows:

D;,=vE;;

T = f — v(ui,kui,kuj),i

f--2v

R P (P,ku.i,k),i (21)

PD = - 2v (( u

i,k uj,k) + (vk,iuk,i))Ui j - 2v(u;ui,k)Ui,k;

- 2v('ui,ku;,ku=,I) - 2v 2 

(ui,klui,kJ)

The term II` is usually neglected and the term T f is modeled as

To model PD, we define %P by

VT l

T` = T E , ;

J

(22)

PD = _ 

EE

k

At the level of the k-E model, we assume T is a function of v, VT, k, E, e, Uij and Ui,jkSince T is an invariant, it must be a function of the invariants that can be constructed

from these quantities. Therefore we can write

vT Ui,; Ui,; k

IP = T (R,, 1/VTUi,JkUi,Jk—)

E EE

where Rt is the turbulent Reynolds number k2 1VE. We expand IP in a Taylor series about

VTUijUijlI and vvTUi jkUi , ; k klEe, and take only the linear terms. We obtain

' I• V vT Ui,i Ui ,i .I, k

= WO + y l E + 4'2vvTUi,JkUi,Jk 

EE

(23)

92
where the coefficients VVo, 01 and 02 are in general functions of Rt . Finally, the modeled

dissipation rate equation becomes

aE aE a LT aE E EE —4- U;

= [( v +—) ^ 4-Crf]- VTUi,l Ui ,.7 —C2 f2— + E (24) at ax; ax; Gr F ax; k k

where, C, and C2 are the model constants, and fl and f2 are functions of Rt . The term

E in the present model comes from the last term in Eq. (23):

E = VIIT[Ti,Ik Ui,jk (25)

where we have taken 02 = —1. The form of E and Cl , C2i f1 and f2 for various k-E models

are listed in Tables 3 and 4.

4.3 Boundary Condition for E.

Many of the earlier k — E models use E = 0 as the boundary condition for the

dissipation rate. It is now generally agreed that this is not the physically correct boundary

condition. However, controversy still exists in what boundary condition should be used for

the dissipation rate. In some calculations, ac/ay = 0 is used. which clearly has no physical

background. Most models use the second derivative of the turbulent kinetic energy at the

boundary as the boundary condition for E, as listed in Table 1. This condition cornes

directly from the k-equation and is physically correct; however, it makes the problem very

stiff and thus put very stringent restrictions on the choice of initial profiles for k and E.

If the initial profile for k is not given correctly, the second derivative of k can become

negative and cause the solution procedure to diverge.

We propose the following boundary condition based on the asymptotic analysis. At

y = 0, it is obvious from eqs. (6) and (7) that

2

E = 2v 

a 

a

(26)

y

This expression is exact at the gall, and it does not add stiffness to the solution procedure.

5. Deficiencies and Improvements of Existing Models

5.1 Damping functions

One of the deficiencies of the existing near wall models is that most of the wall

damping coefficients are functions of a wall coordinate, such as the y+ . This types of

damping function works well only in the cases of attached boundary layer flows with simple

geometries where y+ is well defined. For practical engineering problems with corner flows

or separated flows, some ad hoc treatment to the damping function must be made. The

same is also true for the length scale in an algebraic model. The only exception to the above

is the Jones-Launder model in which the damping function is a function of Rt = k2/vE.

Although the Jones-Launder model has the advantage of independent of y + , it

is known that this model does not predict correctly the near-all turbulence, especially,

93
the near-wall turbulent kinetic energy is underpredicted. One natural cation could be to

modify the J-L model such that it would predict the correct near-wall turbulent quantities.

However, one basic difficulty met in such attempts is that the k-equation is very stiff. The

fact that we are putting k2 back into the equation by using the parameter Rt aggravates

the situation.

To avoid the above difficulty and yet still achieve the same end, we introduce the

following parameter:

U4

R,, 

= ve ,

(27)

where U is the total velocity. (Note: U is the total velocity in a coordinate frame fixed to

the solid boundary.) From the results given in Section 3, since e approaches a finite value

at the wall, it is obvious from eq. (1) that R,, = 0(y4 ) near the wall. Similar to eq. (13)

of Section 3, we write the damping function for the eddy viscosity as:

f^^ = 1 - exp(-a 1 Ru/4 - a2 R,i,/2 - a3 R,,) (28)

with a l = 5 x 10-3 , a 2 = 7 x 10-5 , and a3 = 8 x 10 -7 . One can easily verify that with this

damping function, the eddy viscosity has the correct near wall behavior, i.e., VT = O(y3).

One point worth mentioning is the wide applicability of the above damping function. Though developed with Shih model (Section 3, 4) in mind, it can be used with any

existing k - e model that uses a non-zero boundary condition for the dissipation. This new

parameter R,,, unlike Rt , by no means affect the stability of the solution procedure.

5.2 Pressure transport term.

In order to remove the coordinate dependency of the k-equation, we replace the

pressure transport term given in eq. (18) by the following expression:

L

II Co VT kti J = (29) 2 QA

11 J

,1

where the model constant is readjusted to Co = 0.01.

5.3 The formulation of Z.

In order to obtain the correct wall behavior for the eddy viscosity, we have introduced e` in eq. (11), Section 3. Theoretical analysis shows that e is always positive.

However, in numerical calculations, the value of E may become negative or even oscillatory

due to round of errors. ( Depending on the accuracy of the numerical procedure, this may

or may not be the case.) Here, an alternative definition of e is suggested:

E = (1 - exp(-R,1/2 ))e (30)

This expression has the same near wall behavior as eq. (11) but is less likely to cause

numerical instability.

94
%Vith the modifications suggested above, the model constants need slight adjustment. The constants used in the present modified model are as shown in Table 4.

6. Numerical Testing

Flows with self-similar solutions are particularly useful for accurate model testing,

because their solutions are independent of initial conditions, and one does not need to

choose carefully the initial conditions for the k and E. In this paper, we use a fully

developed channel flow and a flat plate boundary layer for model testing. These flows are

the simplest wall bounded turbulent shear flows with self-similar solutions. However, the

complex features of the turbulence, for example, the effect of the wall on shear turbulence,

are present. In the case of the channel flow, the k-E equations form a one-dimensional

problem, numerical calculations are easy and accurate. Recently, the measurementsl8l

confirmed the accuracy of the direct numerical simulation data. l41 These data are used for

model validations.

In legends of the figures presented at the end of the paper, .wliich will be discussed

in detail in the following paragraphs, the word "present" refers to results obtained using

the model suggested in Section 4 together with the new damping function and e given in

Section 5, while the label "Shill" refers to results obtained using strictly formulations given

in Section 4.

6.1 Fully developed turbulent channel flow

Let h be the half width of the channel, u, the friction velocity and Re, the Reynolds

number defined as u T h/v. Let U, k, E, VT and y be the non-dimensional quantities, normalized by u,, u

2

, uT/h, v and h, respectively. The modeled equations for the channel flow

becuine

AT 1—y

dy Re, 1 + v (31)

d { 1 [,+(1+C)vTJdk}+ VT(

Q

dU)2 1 -

E -0 (J2) dy Re, k dy dy Re,

d 1 vT dE E dU 2 1 EE d 2 2 1 —{ (1+—)—} +C, — vT(—) –C2 f2 + VT( d 2) Reg 0 (33) dy Re, a, dy k dy Re, - y T

where

k2

VT = Cμfμ Re, —E

( dk )2

l y

2k Re,

fμ = equation(13) (34)

0.4 ( Re, k2

f2 = 1 -

1.4 expl-( 6E )2]

_ Co

C fa[1 — exp(—y+)]

95
The boundary conditions are simple. At. the wall,

U=k=0

(dy ) z

( 35)

e 2k Re,

and at the center of the channel,

dk de

= 0 (36) dy - dy

The main results from different k-e models for a channel flow with Re, = 180

are plotted in figures 1 - 4. All the calculations are compared with the direct numerical

simulation data. Figure 1 shows mean velocity profile, Figure 2 shows the turbulent kinetic

energy, Figure 3 shows the turbulent shear stress distribution, and Figure 4 shows the

Dissipation rate. From these numerical results we reach the following conclusions: The

model of Jones and Launderl "I (JL) underpredicts the mean velocity as well as the peak

value of the turbulent kinetic energy. Chien's model [ "] performs better than the JL model,

but it overpredicts the mean velocity near the center of the channel as well as the turbulent

kinetic energy. In these two models, e = 0 at the wall is used as the boundary condition,

so the dissipation rate-near the wall cannot be correctly predicted. Lam and Bremhorst1`1

use a nonzero boundary condition for e and have made some improvement for the mean

velocity and turbulent kinetic energy compared with the results of the JL model. However,

the shear is much overpredicted, and the dissipation rate near the wall is not correct. The

model of Nagano and Hishidal1 I presents a very good prediction for the mean velocity and

shear stress, while the peak value of k is under-predicted. Their main modification to the

JL model is a change in the damping function fμ and the form of E. A zero dissipation rate

at the wall is used. The numerical results from the Shih model and the present modified

Shih model show improvements in the prediction of all quantities, including the dissipation

rate.

6.2 Boundary layer flows

The boundary layer equations and the corresponding k- and e-equations are solved

using a conventional semi-implicit finite difference scheme. In this scheme, the coefficients

for the convection terms are lagged one step in the x-direction, and the source terms in

the k- and e-equation are linearized in such a way that stability is ensured.

In the present study, a grid of 100 points in the y-direction is used. The grid is

stretched linearly with Dye+ 1 /Dye = 1.05. The grads expands in the y-direction according

the the boundary layer growth rate.

The results of this calculation are presented in Figure 5 through Figure 10.

Figure 5 and 6 show the comparison of the wall shear stress from various models to

experimental data and some DNS data. As shown in Figure 5, at low Reynolds number

(based on momentum thickness), J-L model overpredicts the shear stress while Chien

model and NH model under-predict the shear stress. One common character of these three

models is that they all used zero boundary condition for the dissipation, and thus unable

96
to predict the correct turbulence near wall behavior. At high Reynolds number (Figure

6) the JL model still overpredicts the wall shear stress, while the others seem to do a fare

job.

Figure 7 and 8 are the results for Rce = 1410, and Figure 9 and 10 are the results

for Ree = 7700. From these results one can reach similar conclusions as we did from the

channel flow case: The JL model generally underpredicts the peak value of the turbulent

kinetic energy while overpredicts it in the inertia layer; the model also underpredicts the

mean velocity profiles. The LB model and the NH model also under-predict the peak value

of the turbulent kinetic energy near the wall. The three models mentioned above either

have zero boundary- condition for e or do not have the correct order of magnitude for eddy

viscosity. The Figures show that, in general, the present model performs better than the

existing models.

Conclusions

Froin the model testing, we conclude that the present k -e model has made considerable improvement over previous k-e models according to the comparison with the direct

numerical simulation data. We find that the improvement is mainly due to the modified

eddy viscosity model and the model of the pressure transport term in the k quaticn. The

proposed dissipation rate equation also shows a better near-wall behavior than the previous ones. The correct boundary condition for E also seem to play an important role in the

accurate prediction of the turbulence near wall behavior.

Acknowledgement

The first author is grateful to Drs. N.N. Mansour, D. Driver and Professor P. Moin

for their many useful comments. The first author was supported in part by Center- for

Turbulence Research at Stanford; in part by Institutc for Computational Mechanics in

Propulsion, NASA Lewis Research Center. The second author is supported by NASA

Lewis under Contract NAS3-25266.

97
References

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

low-Reynolds-number flows: A Review," AIAA Journal, 23, 1985, pp. 1308-1319.

2. Shih, T.H., "An Improved k - e Model for Near Wall Turbulence and Comparison

with Direct Numerical Simulation," NASA TM-103221, August, 1990.

3. Mansour, N.N., Kim,J. and Moin. P., "Reynolds-Stress and Dissipation Rate

Budgets in a Turbulent Channel Flow," J. Fluid Mech., Vol. 194, 1988, pp. 15-44.

4. Kim, J., Moin, P. and Moser, Robert., "Turbulent statistics in full), developed

channel flow at low Reynolds number," J. Fluid Mech., 177, 1987, pp. 133-166.

5. Lumley, J.L., "Computational modeling of turbulent flows," in Advances in

Applied Mechanics, Vol. 18, 1978, Academic Press, New York, pp.123-176..

6. Baldwin, B.S. and Lomax H., "Thin Layer Approximation and Algebraic Model

for Separated Turbulent Flows," AIAA Paper no. 78-257, Huntsville, Ala., 1978.

7. King, L.S., "A comparison of turbulence closure models for transonic flows about

airfoils," AIAA Paper no. 87-0418, Reno, Nevada, 1987.

8. Nishino, K. and Kasagi, N., "Turbulent statistics measurement in a two dimensional channel flow using a three-dimensional particle tracking velocimeter," Seventh

Symposium on Turbulent Shear Flows, Stanford University, 1989.

9. Spalding, D.B., GENMIX: "A General Computer Program for Two-dimensional

Parabolic Phenomena," 1977, Pergamon Press.

10. Jones, W.P. and Launder, B.E., "The Calculation of Low -Reynolds Number

Phenomena with a Two-Equation model of Turbulence," International Journal of Heat

and Mass Transfer, Vol. 16, 1973, pp. 1119-1130.

11. Chien, K.-Y., "Predictions of Channel and Boundary -Layer Flow with a LowReynolds-Number Turbulence Model," AIAA Journal, Vol. 20, Jan. 1982, pp. 33-38.

12. Lam, C.K.G. and Bremhorst, K., "A Modified Form of the K - e Model for

Predicting Wall Turbulence, "ASME Transactions, Journal of Fluids Engineering Vol. 103,

Sept. 1981, pp.456-460.

13. Nagano, Y. and Hishida, M., "Improved Form of the K - e Model for Wall

Turbulent Shear Flows, "ASME Transaction, Journal of Fluids Engineering, Vol. 109,

June 1987.

98
Table 1 Eddy viscosity and boundary condition for E in various k-E models

Model vT BC: E,,

JL O(y3) 0

Reynolds O(y5) v aayz

LB O(y4) vayz

Chien O(y3) 0

NH O(y4) 0

z

Shih O(y3) 2v ( Wyk/

2

Present O(y3) 2v (aa k)

y

Table 2 Damping functions used in various k-e models

Model fl, f i f2

JL exp(1 +R=/50) 1.0 1 - .3 exp(-Ri )

Reynolds 1 - exp(-.0198Rk ) 1.0 [1 - .3 exp(-R= /9)] fμ

LB [1 - exp(-.0165Rk )] z 1 } (fY )3 1 - exp(- R

2 )

X 

(1 

+ 2R.5 }

Chien 1 - exp(-.0115y + ) 1.0 1 - .22 exp(-R2 /36)

NH [1 - exp(-y+ /26.5)] 2 1.0 1 - .3 exp(-R,2 )

Shih Eq. (13) 1.0 1 - .22exp(-R2 /36)

Present Eq. (28) 1.0 1 - .22 exp(-Rt /36)

Table 3 Model terms in various k-E models

Model II D E

JL

Reynolds

0

0

-2v(aay ^2

0

2vvT(aY2

0

LB 0 0 0

Chien

NH

Shih

0

0

Eq. (18)

- yz' 2

-2v(aay /2

0

- 2y2E exp(-.5y+)

vVT

2

(1 - fμ)(ay2 )

( a2z )z vvT 

ay

Present Eq. (29) 0

a?u 2 vvT( aye

Table 4 Model constants in various k-E models

Model

JL

C,,

.09

C,

1.45

C2

2.0

ak

1.0

QE

1.3

Reynolds .084 1.0 1.83 1.69 1.3

LB .09 1.44 1.92 1.0 1.3

Chien .09 1.35 1.8 1.0 1.3

NH .09 1.45 1.9 1.0 1.3

Shih .09 1.45 2.0 1.3 1.3

Present .09 1.5 2.0 1.3 1.3

Rt = K2 /vE, Rk = I^yl y , y+ = wryly.

99
20

15

10

5

0 L ----

100

4

3

2

1

i i i

Jones & launder

Lam & Bremhorst

Chien

. - - . . - - Nagano & Hishida

-------- I....... Shih

Present

O O O DNS Data

10' 102

y+

Figure 1. Mean velocity profile for channel flow, Re, = 180.

— — _ Jones & launder

Lam & Bremhorst

______ Chien

. . . - Nagano & Hishida

.................. Shih

^ 1- - •. ^ Present

O O O DNS Data

0.

0- \^

O

O

O

O _

t

O ^^_---

O O -

^- ---- -----

O O O O

5

0`

0

y+

Figure 2. Turbulent kinetic energy for channel flow, Re, = 180.

100
0.2

0.1

0.1

0.1

of

0.1

Q3'

0.1

I

0.0

0.0 r

0.0

L

0.0

no.

0

Jones & launder

- - - Lam & - Bre-ry;horsi

- - - - - - Chien

...... Nagano & Hishida

........ I ......... Shih

Present

o o o DNS Data

I.

0.

0.

0.

0.1

0.00

0 zo 4C 60 80 100 120 140 160 180 200

y+

Figure 3. Turbulent shear stress for channel flow, Re, = 180.

— — — Jones & launder

- - - Lam & Bremhorst

- - - - - - Chien

........ Magano & Hishida

.................. Shih

— Present

o o o DNS Data

y+

Figure 4. Dissipation rate for channel flow, Re, 180.

101
r

0.008

0.007

0.006

0.005

C^~

0.004

0.003

0.002

Jones & Launder

Lam & Bremhorst

Chien

- - . . - - . - . Nagano & Nishida

------------- --

Shih.

Present

+ +

+ Smith & Walker

O O O Wieghardt & Ti.11•nann

0 n o Kar-man & Schoenherr

0.ONO,

Jones & -Launder

Lam & Bremhorst

0.007 _ _ _ _ Chien

. . . Nagano & Nishida

...... Shih

Present

0.006 O O O Experiments

.t. + .r Numerical Simutaiion

1 O

0.005 a ^

O ^

i

0.004 ^^^^ — + — — — —

0.003

0.002

0.001' 1 1'1 1 1 1

1'

0 400 800 1200 1600 2000

Re-,

Figure 5. Wall shear stress coefficient, flat plate boundary layer flows, low Reynolds

number data.

0.001

0 2000 4000 6000 8000 10000

Re-,

Figure 6. Wall shear stress coefficient, flat plate boundary layer flows, high Reynolds

number data.

102
5

Jones & Launder

_ Lam & Prem.horst

b ^^ -- - - - - - Chien

. . . Nagano & Hishida

4

Q

. . .

Shih

Present

\ O O o Numerical Simulation

O ^\

3

+ G \

I ^

2 O

dI O \\\ 0 \ \

O

O \

o

0

0.0 0.2 0.4 0.6 0.8 1.0 1.2 y/b

Figure 7. Turbulent kinetic energy for flat plate boundary layer, Reg = 1410.

25

f

I

b- O O

20

^o

15 'O

i

10

/ — — — Jones & Launder

- Lam & Bremhorst

--_--- Chien

_ - - - _ _ - Nagano & Hishvda

5

..... _...--- Shih

Present

O O O Numerical Simulation

0

10 1

02 103

y+

Figure 8. Mean velocity profile for flat plate boundary layer, Reg = 1410.

103
Jones & Launder

Lam & Bremhorst

___--__ Chien .

... . . Nagano & Hishida

................... Shih

Present

O O O Klebnoff

O

O \

--O \

O

O

O

0 000 1_

0.0 0.2 0.4 0.6 0.8 1.0 1.2

y/b

Figure 9. Turbulent kinetic energy for flat plate boundary layer, Rea = 7700.

Mean velocity Profile (Re 19 = 7700

30

25

0.010

0.008

0.006

0.004

0.002

i

20

+ 15

10

5

O

C,

o

/

Jones & Launder

_ _ _ _ Lam & Bremhorst

Chien

... ... .. . Nagane & Hishida

----

......... Shih

Present

O O O Klebanoff

0

101 102 103 104

y+

Figure 10. Mean velocity profile for flat plate boundary layer, Rea = 7700.

104
SOURCE: NASA TM 104368, ICOMP-91-06, CMOTT-91-01

Low Reynolds Number Two-Equation

Modeling of Turbulent Flows

V. Michelassi; T.-H. Shih

Center for Modeling of Turbulence and Transition

NASA Lewis Research Center

March 2S, 1991

Abstract

A new k — e turbulence model that accounts for viscous and wall effects is

presented. The proposed formulation does not contain the local wall distance

thereby making very simple the application to complex geometries. The formulation is based on an existing k — c model that proved to fit very well with

the results of direct numerical simulation. The new form is compared with nine

different two-equation models and with direct numerical simulation for a fully

developed channel flow at Re = 3300. The simple flow configuration allows a

comparison free from numerical inaccuracies. The computed results prove that

few of the considered forms exhibit a satisfactory agreement with the channel

flow data. The new model shows an improvement with respect to the existing

formulations.

105
SOURCE: NASA TM 103222, ICOMP-90-0017

MODELING OF NEAR-WALL TURBULENCE

T.H. Shih*

Institute for Computational Mechanics in Propulsion

Lewis Research Center

Cleveland, Ohio 44135

N.N. Mansour

National Aeronautics and Space Administration

Ames Research Center

Moffett Field, California 94035

ABSTRACT

This paper presents an improved k-c model and a second order closure model for lowReynolds number turbulence near a wall. For the k-E model, a modified form of the eddy

viscosity having correct asymptotic near-wall behavior is suggested, and a model for the pressure diffusion term in the turbulent kinetic energy equation is proposed. For the second order

closure model, we modify the existing models for the Reynolds-stress equations to have proper

near-wall behavior. A dissipation rate equation for the turbulent kinetic energy is also reformulated. The proposed models satisfy realizability and will not produce unphysical behavior.

Fully developed channel flows are used for model testing. The calculations are compared with

direct numerical simulations. It is shown that the present models, both the k-e model and the

second order closure model, perform well in predicting the behavior of the near wall turbulence.

Significant improvements over previous models are obtained.

106
SOURCE: NASA TM 104413, ICOMP-91-09, CMOTT-91-03

ADVANCES IN MODELING

THE PRESSURE CORRELATION TERMS

IN THE SECOND MOMENT EQUATIONS

Tsan-Hsing Shih and Aamir Shabbir

Institute for Computational Mechanics in Propulsion

and Center for Modeling of Turbulence and Transition

Lewis Research Center

Cleveland, Ohio 44135

and

John L. Lumley

Cornell University

Ithaca, New York 14853

ABSTRACT

In developing turbulence models, different authors have proposed various model constraints in an attempt to make the model equations more general (or universal). The most

recent of these are the realizability principle (Lumley 1978, Schumann 1977), the linearity

principle (Pope 1983), the rapid distortion theory (Reynolds 1987) and the material indifference principle (Speziale 1983). In this paper we will discuss several issues concerning

these principles and will pay special attention to the realizability principle raised by Lumley (1978). R.ealizability (defined as the requirement of non-negative energy and Schwarz'

inequality between any fluctuating quantities) is the basic physical and mathematical principle that any modeled equation should obey. Hence, it is the most universal, important

and also the minimal requirement for a model equation to prevent it from producing unphysical results. In this paper we will describe in detail the principle of realizability, derive

the realizability conditions for various turbulence models, and propose the model forms

for the pressure correlation terms in the second moment equations. Detailed comparisons

of various turbulence models (Launder et al. 1975, Craft et al. 1989, Zeman and Lumley

1976, Shih and Lumley 1985 and one proposed here) with experiments and direct numerical simulations will be presented. As a special case of turbulence, we will also discuss the

two-dimensional two-component turbulence modeling.

107
Advancements in Engineering Turbulence Modeling

T.-H. Shih

Center for Modeling of Turbulence and Transition

ICOMP/NASA Lewis Research Center

9th National Aero-Space Plane

Technology Symposium

November 1-2, 1990

Paper Number 105

108
ABSTRACT

Some new developments in two-equation models and second order closure models will

be presented. Two-equation models (e.g., k-e model) have been widely used in CFD for

engineering problems. Most of low-Revnolds number two-equation models contain some

wall-distance damping functions to acount for the effect of wall on turbulence. However,

this often causes the confusions and difficulties in computing flows with complex geometry

and also needs an ad hoc treatment near the separation and reattachment points. In

this paper, a set of modified two-equation models is proposed to remove abovementioned

shortcomings. The calculations using various two-equation models are compared with

direct numerical simulations of channel flows and flat boundary layers.

Development of second order closure model will be also discussed with emphasis on the

modeling of pressure related correlation terms and dissipation rates in the second moment

equations. All the existing models poorly predict the normal stresses near the wall and fail

to predict the 3 dimensional effect of mean flow on the turbulence (e.g., decrease in the

shear stress caused by the cross flow in the boundary layer. The newly developed second

order near-wall turbulence model to be described in this paper is capable of capturing the

near-wall behavior of turbulence as well as the effect of three dimension mean flow on the

turbulence.

109
1. k-E model

The two-equation model, especially k -E model, is still the most widely used model

for computing engineering flows. We first list some of the commonly used two-equation

models,[1],[2],[a),[a],[sl,[s] and their predictions on the fully developed channel flows and

boundary layer flows compared with the corresponding direct numerical simulations. Then,

we propose a modified k -E model which does not contain any wall distance. The proposed

k-E model has been also tested using direct numerical simulation data.

The eddy viscosity VT is assumed in two-equation models as follows:

k2

VT = Cμ fμE

or

VT = C,, fμ kT

where T = k2/e

The general k-E (or k-T) model equations are of the following forms:

k,t + Ujk,j K VT + Uk v/ k1J J

E ,t + Uj E,J - VT

^ \ QE + v/ E^JJ ,7

T + UJT,J - VT

K U r2 + v^ T'JJ ,j

+ Il + VT (Ui , j + Uj, i) Ui ,j - E + D

+ CI k VT ( Ui,j + Uj,i) Ui,j - C~2 f2 k + E

+ 2 ( VT + v) k,iT,i - 2 C vT + v I T,iT,i

k Uri T Uri

+ (1 - CE1) 7kVT(Ui,j + Uj,i)Ui,j + ( C ,2f2 - 1)

110
This table lists the model terms, damping functions and model constants appeared in

various two-equation models

Model B E Cil Cie fl {

J2

JL 2vvi(—;)2

E 1.45 2.0 1.0 1 - .3e -

LB o E 1.44 2.0 1+( f ) 2 1-e-^

Chien -21E E 1.35 1.8 1.0 _ 4 

1 - .22e 36

NH vvt(1 - ( 31 -IUI ) 2 E 1.45 1.9 1.0 1 - .3e

Shih 1/.j(")'

071 

'1 " R' 2 1 . 45 2.0 1.0 1 — .22e — 36

Model H D f^ O'k

JL o —2v (aa )2 e^ 1T 5o)

l.o 

LB O O (1 — e -.ol65R12( 1 + 20.,')

1 R

1.0

Chien o —21/k /y2 e(--0155y 1 — ') 1.0

NH O —2v(^ ) 2 [1 - e(-y+/ 26.5)12 1.0

SAA o — O (1 +"- 15

le

 )tank( o)

v

1.36

Shih .05vik,i

CrkfN[ —exP(' y+)j

---]

O

1 — e(-6xlo-3y+ -4xlo-4y + {

2.5 x 10 -6y +3 -4 x to-9y+4)

1. 3

T/t = CItf/, k J+ = u-y ^L = P = 

Vk

v !v E v
o us

— — X.

- 111

--- Chien

----- wl

Pr es en C

0

N

0

r

e

0

0

0

The following figures show the predictions on fully develop channel flow using various

two-equation models compared with the direct numerical simulation data. [7) Ploted quantities include the mean velocity U, turbulent shear stress (uv), turbulent kinetic energy k

and dissipation rate EPS c. The open circle represents direct numerical simulation, and

the solid line represents the model prediction.

0

4

N

0

0

0

1e 0.0 25.0 50.0 7S0 1000 1250 1500 15.o 2000

Y+ y+

0 x

a

o 11

I °.

N

0

0

0

f. Is. /°. W W tM M. 1^°. IN. 11°

0.0 25.0 50.0 75.0 100D 125.0 150.0 173.0 200.0

Y+ y*

112
A" t A---

- - - -- - - Gve

SAN

Y

0 10 10 .o 10 100

y'

z

/

.ii

Sy., . A" . A—.--

J

Qi.w

Xyw k Xl,Aida

SAN

.. ♦ . ♦ J/wvri^al S:+nWntS...

10' r P

y ♦

The figures below show the model predictions on flat plate boundary layer flow. A

direct numerical simulation of boundary layer flow [8] is used for comparison. The skin

friction coefficient C f is also included in the comparisons.

Mean velocity profile (P.ey - 141C Turbulent kinetic .energy (Rey = 1410)

0.e \

♦ \ -

F

`e \

`A \

.. _ ... Msew . XI,AtN

JAN

♦ ♦ ♦ Xvw.nal SiwWtvw

0.0

ao os a. az az 1 .0 1s

y/a

I. ve• 1

SuJan., . (aw{.r

^^

r

__`

_ _ L++n . Ir.Mw,1

^

_ Gvn

I^ .... Map— . X4ti4. SAN

LODe I r1, O O O N 4^...4r4n

^'•,l O

7.003 I O

I _

7.003

I

1 I

1 I

a^J I UNCLASSIFIED

\ I

\ I

....a

.00 zao rzao r.oa zoos

Rey

113
Overall, Shih's k-E model gives better prediction in both fully developed channel flows

and flat plate boundary layer flows according to the comparisons with corresponding direct

numerical simulations. However, this model has the same problem as the others, that is it

contains the wall distance parameter y + , which is defined as

+ ury

y = —

v

where u T is the friction velocity. The difficulty would occure in some situations. For

example, near the seperation point uT aproaches zero and hence vt (through fμ (y+ )) will

approach zero everywhere when this u T is used. Another example is the flow with complex

geometry that the wall distance is not well defined. In the both cases, the ad hoc treatment

is needed in the model implimentation. We notice that Jonse-Launder's model [ ', does not

contain y+ . However, its present form does not perform very well in the simple testing

flows. Here, we based on the Shih's k-E model modify the parts which is related the y+

with another parameter R = k E^z L . R is a ratio of turbulence length scale to viscous

length scale. The modifications[91 made here are:

fμ = 1 — exp {C3 I1 — exp(C6 R1/4 )] }

II=

 ^

Co VT k

f μ ^^ ,J

r 1/2

E = E I 1 — exp(—Rt )

where Rt = k2 /vE, Co = .004, C3 = .0004, C6 = 1.2

114
The following figures show the predictions from the present model on fully developed

channel flow compared with other models (including Jonse and Launder's modelN ) and direct numerical simulation dataP ] The open symbols represent direct numerical simulation,

and the lines represent the model prediction.

5.0

4.0

3.0

x

2.0

1.0

0.0 Z

0.0

40.0

35.0

30.0

r^

25.0

m

"5

.^ 20.0

a^

p 15.0

10.0

5.0

0.0

0.0

50.0 100.D 150.0 0.0 50.0 100.0 150.0 200.

Y+ Y+

50.0 10010 150.0 00 50.0 100.0 150.0 200.0

Y+ Y+

115
50.0 100.0 150.0 C.0 50.0 100.0 150.0 200.

Y+ Y+

0.0

10-1.0

0.80

0.70

0.60

0.50

I^ 0.40

0.30

0.20

0.10

0.00

0.0

20.0

16.0

12.0

8.0

4.0

,o

'o - O

q DNS q DNS

ch --- lb

-------- jl -------- sh

- nh present

100.0 101"° 102.0 1 a-1.0 100.0 101.0 102.0

Y+ Y+

116
2. Second order modeling of near-wall turbulence

Using the near-wall asymptotic behavior of turbulence l10l as model constraints, we

formed a set of modeled transport equations for the Reynolds-stress tensor and the dissipation rate of turbulent kinetic energy. The main emphasis was on developing a near-wall

model for the pressure correlation and dissipation terms in the Reynolds-stress equation.

A modeled dissipation rate equation is derived more rationally. Asymptotic analysis shows

that near the wall, the viscous diffusion term in the Reynolds-stress equations becomes

the leading term and is balanced by the pressure correlation and dissipation terms. We

use this as a model constraint in the model development. The proposed models satisfy

realizibility and ensure no unphysical behavior will occur. Here, we briefly describe and

list the proposed models.

Reynolds stress equation

The exact equation for the Reynolds stress tensor is:

Dt ( u = uj) = Pij + Tij + Dij + HIj — Eij

where ( ) stands for an ensemble average, D/Dt = a/at + Uka/a2k. the terms Pij , Tij,

D(v) , rlij and e ij represent the production, turbulent diffusion, viscous diffusion, velocity

pressure-gradient correlation and dissipation tensor, and are identified as follows:

Pij = — ( u i u k) Uj, k — (ujUk)Ui,k

Tij = —(uiUjUk),k

D(^) = v(uiuj),kk

17 ij = —P (u iP,j + ujP,i)

Eij = 2v(ui,kuj,k)

The proposed near-wall model for n ij — E ij is:

Ilij — E ij = —fw 

(42) 

[2(u i uj) + 4(( u i u k) njn k + (ujuk) nink ) + 2(uk uj)nknrninj]

where ni is a unit vector normal to the surface, and fw = exp(—(Rt/C1 ) 2 ), Rt = 9vf ,

Cl = 1.358RO . 44 , R,! ,, = u r b/v. uT is the friction velocity, b is the thickness of the boundary

layer or the half width of the channel. -

Away from the wall, the velocity pressure-gradient correlation N ij is split into the

rapid part W' ) and the slow part 11(2):

llij = rl(1) + rl(2)

117
The proposed model for IIi,- Eij is:

( 2 {

—E (Q b

ij + 3 6ij )( 1 — JiU)

where

F 72

/j = 2 + 9 1 2 + 80.1 ln[1 + 62.4(—II + 2.3111)]} exp(— 7.77

)

Rt R 

1/2

t

F = 1 + 27111 + 911

II = —1bijbji

III = 13bijbjkbki

bij = (uiuj ) / ( g2 ) — bij/3

The rapid part of velocity pressure-gradient, H(1 ) is modeled as follows(Shih and

Lumley [11,12) ):

II^^ ) _ ( 5 + 2a5)( g2 )( Ui,j + Uj i ) — 3(1 — a 5 )(Pij — 2Pbij)

2 2 2 6

+ (3 + 316 

a5)(Dij — 2Pbij) + 

2 (Pij — Dij) + 6bijP

+ 5(g2) [(( U i Uk ) UJ,4 + ( UJ Uk ) Ui,9)( UkU 9) — ( U = Up ) (Ul U 4)( Up ,9 + U4,p)]

where,

Pij = — ( UiUk ) Uj , k — (UjUk)Ui,k

Dij = —(uiUk)Uk ,j — (Ujuk)Uk i

P = 

1

2

Pi i

a5 = — 10(1+C2F1/2)

C2 = 0.8[1 — exp(-(Rt/40)2)]

Finally the model for the third moments is modeled as:

(uiuj uk) = —.07(22) [(UkUp )( UiUJ),p + (U) UP)(UiUk ) ,p + (UiUp)(UJUk),p]

Dissipation rate equation

The modeled dissipation rate equation derived in this work is:

•/, E ,t + UiE, i = ( 1/E,i — (EUi)),i

EE

— `^0 (g2)

— 01 E ( (UiUj ) Ui,j — 02 v (q2) 

(Ukul)(Ui,jl — Ul,ij)Ui,jk

q 

2)

118
where

00 = 5 + 0.98[1 — 0.331n(1 — 55II)] exp(- 2.83R, 112)

Y'1 2.1

02 = —.15(1 — F)

y(42),i(g2),i

4(42)

The turbulent flux term (EU k ) is modeled as:

(Euk) = -.07 (q22) 

E (ukup)E,p

These figures show some existing Reynolds stress models (for example, Launder and

ShimP 31 , Lai and Sof "I ) and present model compared with the direct numerical simulations.^71

Mean velocity U/u'

70.o

]u 0

IoA

00

i0 ia•' l0 l0"

7•

Turbulent kinetic energy k/u'/u'

s.o

•.. - 160

.0

]A

S

+ 7A

IA

0.0

Turbulent stress —<uv>/u'/u'

0.7 p- 100

.. o.rs

^ — . lw../.r k Any

0.1

L

3

01

0.0 1 O 0.0 60.0 100.0 ,30.0 200.0

UNCLASSIFIED

Intensities

70

.. ons

7A

i`

IA

0.0

00 SOD IOOA 150.0 }00.0 QO 5010 100.0 150.0 200 0

Y.

119
o

1d

N

0

kleon velocity U/u' Turbulent sbass --<uv>/u'/u'

ta.a

m.o

^ loo

ao

o.0

0-u 10" 10" 10" 10" _ 0'0 OOA 100.0 IJp.O ]00.0

1 I•

UNCLASSIFIED

Turbulent kinetic energy k/u'/u' Intensities

5o

J.o

• ].o

00

0.0 50.0 IOOA 150.0 200.0 0.0 00.0 100 0 150.0 ]00.0

U - MEAN VELOCITY K - KINETIC ENERGY

0 0

present

t,! i

r ^/

10= 0.0 25.0 50.0 75.0 100.0 12-50 150.0 175.0 200.0

Y+ y+

u, v, w -- RMS OF FLUCTUATING VELOCITY <UV> PROFILE

0 0

0

N

0

0

a7

a^

6T

q

0

J aJ

q

0

0

0

J

0.0 25.0 50.0 Tao 100.0 125D 130.0 17'5.0 200.0 0.0 MD 50.0 7.0 100.0 123.0 150.0 175.0 200.0

ii.I- y I.

La

os

S a.^

S

t

of

o.o

JA

lA

IA

120
3. Second order modeling of a three-dimensional boundary laver

A study[ "] of three-dimensional effects on turbulent boundary layer were achieved by

direct numerical simulation of a fully developed turbulent channel flow subjected to transverse pressure gradient. The results show that, in agreement with experimental data [ "] , the

Reynolds stresses are reduced with increasing three-dimensionality and that, near the wall,

a lag develops between the stress and the strain rate. To model these three-dimensional

effects on the turbulence, we have tried various two equation models and second order

closure models. None of the current models can predict the reductions in the shear stress

observed using direct numerical simulations. However, we found that the newly proposed

second order closure model listed in the previous section do at least qualitatively capture

these three-dimensional effects.

121
O

O

OO

O UNCLASSIFIED

CO I

0

0

+ cn

^o

0

I^ o

I^

I

P o

^3 o

The following figures show the direct numerical simulation of the three dimensional

boundary layer flow and the model predictions from Launder and Shima, Lai and So and

the present models.

`^ dA

d^ ap Y

a2 2

Schematic of three-dimensional channel flow.

0

. I UNCLASSIFIED

C\2

se.

W'

0

ri

V'

0

0 i

0.0 0.2 0.4 0.6 0.8 1.0

Y16

Turbulent intensities in fixed coordinates at At =

0.15 intervals.

vw

0 uv

I

0.0 0.2 0.4 0.6 0.8 1.0

y/s

Reynolds stress component in fixed coordinate at

At = 0.15 intervals.

0

0

—20.0 0.0 20.0 40.0 60.0

13, -ys , -y-r , (degrees)

Distribution across the channel of the mean flow

angle (--- ), the Reynolds shear stress angle (----) ,

and the mean velocity gradient angle ( ).

122
OFF

LL

>, -Z

cd

G E

E

16 o

O

O

O O

cq

O

rQ

O

O

O

0

O

0 

+

V^l

O

O

O

(D

va

0, 0 zo— 1"0— go— 1 3 0— O'l-

<'LT"L>

0 p O"C O*Z . O'l 00 

O

(n 'a 'n

123
A

d

,d

.. S

u°

------ .

.

o e ...................w ....... _.fro__.._.._......

.................

A: ? o --'--.... gip....;..._...... __._..___....._..._...

r° 90

---.—•—' — — -- =----- — —'--

S > sN

....... _ ......—._....Aa.___.^...__—_—......._........ c

0

0

o

^^l

EF

O

pO

b

Q

O

LL

- - Un

ICI

U

Z

i

i'

y

4

4.

^R

0

124
References

1. Jones, W.P. and Launder, B.E., "The Calculation of Low-Reynolds Number Phenomena with a Two-Equation model of T`urbulence," International Journal of Heat and

Mass Transfer, Vol. 16, 1973, pp. 1119-1130.

2. Lam, C.K.G. and Bremhorst, K., "A Modified Form of the K-e Model for Predicting

Wall Turbulence, "ASME Transactions, Journal of Fluids Engineering, Vol. 103, Sept.

1981, pp.456-460.

3. Chien, K.-Y., "Predictions of Channel and Boundary-Layer Flow with a Low-ReynoldsNumber Turbulence Model," AIAA Journal, Vol. 20, Jan. 1982, pp. 33-38.

4. Nagano, Y. and Hishida, M., "Improved Form of the K - e Model for Wall Turbulent

Shear Flows,"ASl11E, Transaction, Journal of Fluids Engineering, Vol. 109, June 1987.

5. C.G. Speziale, R. Abid and E. Clay Anderson, "A Critical Evaluation of Two-Equation

Models for Near Wall Turbulence," AIAA 21st Fluid Dynamics, Plasma Dynamics and

Laser Conference, June 18-20, 1990/Seattle, WA.

6. T.-H. Shih, " An Improved k-c Model for Near-Wall Turbulence and Comparison with

Direct Numerical Simulation," NASA TM 103221, ICOMP-90-16.

7. Kim, J., Moin, P. and Moser, Robert., "Turbulent statistics in fully developed channel

flow at low Reynolds number," J. Fluid Mech., 177, 1987, pp. 133-166.

8. P.R. Spalart,"Direct simulation of a turbulent boundary layer up to R B = 1410," bf

J. Fluid Mech., vol. 187, 1988, pp. 61-98.

9. V. Michelassi, personal communication, 1990.

10. Mansour, N.N., Kim,J. and Moin. P., "Reynolds-Stress and Dissipation Rate Budgets

in a Turbulent Channel Flow," J. Fluid Mech., Vol. 194, 1988, pp. 15-44.

11. Shih, T.H. and Lumley, J.L., " Modeling of pressure correlation terms in Reynoldsstress and scalar flux equations," Rept. FDA-85-3, Sibley School of Mech. and

Aerospace Eng., Cornell University.

12. Shih, T.H. and Lumley J.L., " Second-order modeling of near-wall turbulence," Phys.

Fluids, 29, 1986, pp. 971-975.

13. Launder, B.E. and Shima, N., "Second-Moment Closure for the Near-Wall Sublayer:

Development and Application," AIAA Journal, Vol. 27, Oct. 1989, pp. 1319-1325.

14. Lai, Y.G. and So, R.M.C., Int. J. Heat Mass Transfer, to appear(1990).

15. Moin, P., Shih, T.H., Driver, D. and Mansour, N.N., "Numerical Simulation of a -

Three-Dimensional Turbulent Boundary Layer," AIAA 89-0373.

16. Bradshaw, P. and Pontikos, N., "Measurements in the Turbulent Boundary Layer on

an Infinite Swept-Wing," J. Fluid Mech, 159, 1985, pp. 105-130.

125
SOURCE: NASA TM 104369, ICOMP-91-07, CMOTT-91-02

Second Order Modeling of

Boundary-Free Turbulent Shear Flows

T.-H. Shih

Institute for Computational Mechanics in Propulsion

and Center for Modeling of Turbulence and Transition

Lewis Research Center

Cleveland, Ohio

J.-Y. Chen

Sandia National Laboratories

Livermore, California

J.L. Lumley

Cornell University

Ithaca, New York

Abstract

This paper presents a set of realizable second order models for boundary free turbulent

flows. The constraints on second order models based on the realizability principle are reexamined. The rapid terms in the pressure correlations for both the Reynolds stress and

the passive scalar flux equations are constructed to exactly satisfy the joint realizability.

All other model terms (return-to-isotropy, third moments and terms in the dissipation

equations) already satisfy realizability (Lumley 1978, Shih and Lumley 1986). To correct

the spreading rate of the axisymmetric jet, an extra term is added to the dissipation

equation which accounts for the effect of mean vortex stretching on dissipation. The test

flows used in this study are the mixing slicar layer, plane jet, axisvmmetric jet and plane

wake. The numerical solutions show that the new unified model equations (with unchanged

model constants) predict all these flows reasonably as the results compare «ell with the

measurements. lVe expect that these model equations would be suitable for more complex

and critical flows.

126
A. Shabbir Evaluation of Turbulence Models

Visiting Assistant Professor. for Predicting Buoyant Flows

D. B. Taulbee

Professor. Experimental data for the buoyant axirymmetric plume are used to validate certain

Mem. ASME closure hypotheses employed in turbulence model equations for calculating buoyant

flows. Closure formulations for the turbulent transport of momentum, thermal

Department of Mechanical energy, kinetic energy, and squared temperature used in the k-f and algebraic stress

and Aerospace Engineering, models are investigated. Experimental dotafor the mean velocity, mean temperature.

University at Buffalo. SUNY and kinetic energy are used in the closure formulation to obtain Reynolds stresses,

Amherst, NY 14260 heat fluxes, etc., which are then compared with their measured values.

1 Introduction

Various turbulence models have been formulated for prcdicting buoyancy-driven flows. Some of the parameters in these

models have been determined by keying the solution of the

model equations to experimental data for certain basic flows

Such as decay of grid turbulence. Other parameters have been

determined by calibrating closure formulations directly with

experimental data- However, this approach may be somewhat

inaccurate due to the lack of quality experimental data for

certain correlations, especially dissipation. Finally, certain

model parameters have been fine tuned or determined by requiring that the computed solution agree with experimental

data for more complex flows, such as shear flows. In addition

there have been instances where model parameters have been

adjusted or empirical corrective terms added so that agreement

with experimental data is accomplished for a particular flow.

When model parameters are adjusted to get agreement, say

for the mean velocity and temperature fields for a particular

flow, little regard is given for the internal integrity of the model.

In other words, are the various processes such as diffusional

transport, pressure-strain interactions, etc-, predicted correctly? Or arc there compensating assumptions where one process is overpredictcd at the expense of another and yet the end

predicted result for the mean flow agrees with experiment?

The lack of complete sets of data for higher moments, dissipation, and pressure-velocity correlations for various flows

has prevented detailed verification of closure models for the

various processes that have to be modeled.

The objective of this paper is to use the recently obtained

and comprehensive experimental data of Shabbir and George

(1987) and Shabbir (1987) on the axisymmetric buoyant plume

to assess the various closure relations proposed for the kineticcncrgy/dissipation and the algebraic stress models for buoyancy-dominated flows. The usual approach is to solve the

modeled differential equations numerically, and then compare

the computations with the experiment. However, this method

does not help pinpoint the drawbacks in the various terms of

the models. In this paper, instead of the usual approach, correlations obtained from measured velocity and temperature

are used directly to verify the closure hypotheses for the turbulcnt transport of momentum, thermal energy, and turbulent

kinetic energy.

2 Experimental Data

The data used were taken in an axisymmctric buoyant plume

by Shabbir and George (1987) and Shabbir (1987), who mcasContriSutcd by the Hut Tnnsfcr Division and presented at the National ttut

Tnnsfcr Confrrenee, Pittsburgh, Pennsylvania. August 9-12, 1987. Manuscript

received by the Hcat Transfcr Division March J, 1988; rrision received November 10, 1989. Keywords: Modeling and Scaling. Plumes, Turbulence.

ured velocity and temperature fields at several vertical levels

above a heated source of air. Here we .briefly summarize their

experimental technique and results.

The three-wire probe used consisted of a cross-μire and a

temperature wire. Thus the instantaneous values of the two

velocity components (vertical and radial) and temperature were

measured. The axisymmetry of the flow was established by

using an array of 16 thermocouples and also by rotating the

cross-wire by 90 deg. Profile for the correlations between the

velocity components and velocity components with temperaturc through the fourth order were determined from the instantaneous measurements.

Source conditions were continuously monitored in order to

calculate the rate at which buoyancy was added at the source.

The source Grashof number was 5.5. By integrating the mean

energy equation, an integral constraint can be obtained for a

buoyant plume. For a neutral environment this constraint implies that the rate at which buoyancy crosses each horizontal

section is constant and must equal the rate at which buoyancy

is added at the source, i.e., the ratio

F 1 — ll

F F, 

[2, ^o g,6(UAT+ur)rdr

1

(1)

must be unity (F, is the source buoyancy)- This integral constraint was satisfied within 7 percent.

The correlation profiles at various heights were found to be

similar in the coordinate n = r/z (z accounted for the virtual

origin) when the vcloaty is scaled by U, F."z" and the

temperature is scaled by T, = Fo az

-sn/gs. The measurements

agreed well with the earlier study by George et al. (1977), who

measured only the temperature and the vertical component of

velocity. The scatter in the measurements of higher moments

is typical for such flows and is also present in previous experiments, such as those of George et al. (1977). The primary

reason for the scatter is that slow time scales of the flow require

much longer averaging time for the higher moments in order

to obtain the same statistical convergence as for the mean

quantities. Other errors in the measurements arise from the

flow reversal on the hot wire—a phenomenon most likely to

occur toward the outer edges of the flow where local turbulent

intensities arc considerably higher. These are discussed in Shabbir and George (1987).

The various correlations in similarity variables were fitted

with curves using a least-squares fitting procedure. This rcpresentation allows easy evaluation of the terms in the governing

equations and closure formulations when they arc cast in similarity variables. Using these profiles the balances for the mean

momentum and energy differential equations were carried out

to check whether the flow satisfied the equations of motion it

is supposed to represent. Within the thin shear layer and the

Boussincsq assumption the mean momentum and energy equa127

Journal of Heat Transfer NOVEMBER 1990, Vol. 112 1945
-41C

e

6

z

0

J 

•40

T r(9Avt)

.2e

o

} ERROR

r

dr 9A( T - TJ5 ¢

J _2O

f

N

-w

rs0az 9A(T- T_)

0.04 r/z 0

.12 e.2e

«z ar(ruv)

+a

W

K •a

r

C

N

-a

r/z .

Fig- 1(t) Mechanical and thermal dissipation profiles! Full lines are

experimental; the chained line Is from model equation (10), which was

solved for a with all other quantities taken from experiment.

I ERROR

- ak -_ U –q2u;+–pull +P+G-e (4)

f ax

l

ax

J

2 P

where P = –uTi au,/ax, is the mechanical production and

G = –0gY,7 is the production due to buoyancy, Each term

except for the dissipation e and the pressure transport pu; is

determined from the experimentally determined correlations.

0.04 r/z o.rz o.zo

The pressure transport was evaluated from pu//p = –gzuj/

y= 5, a formula given by Lumley (1978). Although this closure

Fig- 1(a) Balances of mean energy and momentum equations (taken relation has not been verified experimentally, it was felt that

from Shabbir and George, 198n since the pressure transport is significant, some correction

should be included rather than simply neglecting it, as is often

tions can be respectively written as done. The dissipation determined from the balance of the

au aU 1 a turbulent kinetic energy equation is shown in Fig - 1(a) as a

U az + V ar r ar (r"")–gSAT

(2) solid line.

By a similar procedure the dissipation of the mean-square

aAT aAT _ 1 a temperature t' is determined from

U az + V ar r ar (rte) (3) arr a oT U — _ – u^ –2 i77 — –2e (5)

Since all the quantities appearing in these equations are meas- J axe axe ' ax`

ured, their profiles were substituted to see whether the mealuremcnts balance the equations. Figure 1(a), taken from All terms are evaluated from e xperi me ntal data and the reShabbir and George (1987), shows that the experiment satisfies suiting thermal dissipation is shown 

in

n Fig. ig. 1(b).

this nontrivial test within 10 percent. An error of such mag- The time scales qz /e and r2 /c, for the relaxation of the menitude is typical of turbulent shear flows. chanical and thermal dissipation, respectively, are shown in

The dissipation of mechanical energy was determined by Fig. 2, along with their ratio

balancing the turbulent energy equation R= (rz/e,)/(q^/e) (6)

Nomenclature

F. = buoyancy flux, equation (1)

g = acceleration due to gravity

G = turbulence production from

buoyancy

k = turbulent kinetic energy

p = fluctuating pressure

P = turbulence production by mean

flow

Pr T = turbulent Prandtl number

r = radial coordinate

R = time scale ratio, equation (4)

f = fluctuating temperature

T = mean temperature

u = fluctuating axial velocity component

U = mean axial velocity component

u = fluctuating radial velocity component

z = vertical coordinate

,6 = coefficient of thermal expansion

F = dissipation of mechanical energy

E, = dissipation of mean-square

temperature

P T = turbulent eddy viscosity

P = mean density

946 1 Vol. 112, NOVEMBER 1990 128 Transactions of the ASME
o ur(Rin

. , EL

—.— 1ILru.,C "I"S —tL

< ,0 0

/ \ o

\ o

i \ 0

0.04 r/z 

0.12 0.20

Fig. 4 Axial heat flux

2_'

rho JA 2.0

o.^

u

I =- 

0.1

I?

t0

0.0

0.3 0S

r= ^.2

r

a.l

0.w 0.08 0.12 0.16

r/z

Fig. 2 Mechanical and thermal time scale ratios and variation of R as

.calculated from experimental data

i

O uruin(xr

w1LL Ic STiCSs _CL

)-3-r^—^\\

0 ^, o

r ^\ 0

r

Oi D.Z % r o

r

r

r

r

J 1 r ^^

0 - 04 0.12 0.20

r1z

Fig. 3 Shear stress

These time scales appear extensively throughout the model

formulations and will be further discussed in the next sections.

3 Assessment of Closure Hypotheses of k-c Model

The form of the k-e model, which is considered to be the

standard one, is that used by Launder and Spalding (1974).

In this model the Reynolds stress is given by

(a Qj aU\ - 2

- —uu u^ = v T + J - kbu (7)

axe ax; 3

and the heat flux by

vT aT

-^ PrT ax;

(8)

where PT = C kt /c and C„ = 0.09 (see Launder and Spalding,

1974).

O unalRrxr

K-- —EL

—.— ALL .0 MUS —EL

1.5

0 A

0 \ ^_

0.04 0.12 0.20

r^z

Fig. 5 Radial heal flux

By invoking the thin shear layer assumption for a buoyant

plume, the above relations reduce to

au

uu= - v T -

ar

aT

U7 _ - {v TIP rT) 

aZ

aT

H = - (vT/PrT) ar

Taking PrT = 1.0, the right-hand sides of the above equations

were evaluated experimentally. These are compared with measurcd values of Yu-, i7, and W in Figs. 3-5. The points are

experimental values and the chain lines are from the model.

The modeled values of uu and 0 compare reasonably with

the experimental profiles except in the outer portion of the

curves. On the other hand the modeled profile of vertical heat

flux 7 is much smaller than the experimental one. It is well

known that the simple gradient models given by equations (7)

and (8) with an isotropic eddy viscosity arc inadequate for

determining streamwise turbulent momentum and heat fluxes.

Usually these quantities do not influence the prediction for

shear layers since only the radial fluxes arc important in these

flows. However, in the case of the buoyant plume the flux

lD is a dominant production term in the turbulent kinetic energy

Journal of Heat Transfer 129 NOVEMBER 1990, Vol. 112 / 947
r c - lo

,.r

---- "orgy

/

1-1

o

.o. _w .0 .11 .2u

Fig. 6 Transport of kinetic energy

equation and its correct calculation is very important for accurate prediction of k.

The diffusional transport in the kinetic energy equation (4)

for a thin shear layer is modeled as

I — 1 _ ak

2 uq2± pu -1 ` ar (9)

P

Using pu/p = - uq2 15 from Lumley (1978) gives uq2 12 =

- (513)v 7-ak/ar from which the result, with the right side evaluated from experimental results; is showrl in Fig. 6. It is seen

that the predicted and experimental data peak at different

radial locations; however, the predicted magnitude is more

accurate, which indicates that the pressure diffusion needed

to be taken into account.

In the k-e model the dissipation is calculated from

_ z U, aE _ - d 

e'u/+C,, 

e (P+G) -C, 2 E (10)

ax; ax; k k

where e'u; _ -(vr-/a,)aE/ax; and a, = 1.3, C,, = 1.44, and

C, 3 = 1.92 as given by Launder and Spalding (1974). In order

to get an indication of the validity of equation (10), it was

numerically solved for the dissipation e with all other quantities

needed to evaluate the coefficients determined from the experimental cori:clations. The result is shown in Fig. 1(b) and

it is seen that it compares reasonably well with the curve obtained from balancing the turbulent kinetic energy equation

with experimental data.

Launder et al. (1972) showed that the standard k-c model

yields a solution for the axisymmetric jet that overpredicts the

spreading rate by about 30 percent. The standard k-c model

also does not correctly predict the axisymmetric buoyant plume

(Hossain and Rodi, 1982). Proposals have been made (Pope,

1978; Hanjalic and Launder, 1980) for modifying the dissipation equation, based on arguments concerned with vortex

or eddy structures characteristic of axisymmetric flows. The

modified equation produces more dissipation, thus decreasing

the turbulent eddy viscosity, which results in a smaller spreading rate of the flow. Here we use the empirical correction given

by Rodi (1972) where C, 2 = 1.92 (1-0.03511) with H = 1 (yE/

U,,,)dU, Idxl° -2 and where U,,, is the maximum velocity and

YE is the distance from the centerline to the edge of the shear

layer. This correction decreases the destruction term in the

dissipation equation, hence producing an increased dissipation.

However, when this correction is used in equation (10) there

is very little change in the solution for f when experimental

data are used for the other quantities in the equation. This is

probably due to the approach taken here, which does not allow

for the nonlinear interactions between the various terms in the

closure_ If the kinetic energy and dissipation equations arc

solved simultaneously, then the axisymmetric correction will

produce a significant change in the solution of the k-< model.

4 Assessment of Closure Hypotheses for Algebraic

Stress itiodel

Chen and Rodi (1975), Tamanini (1978), Chen and Chen

(1979), and Hossain and Rodi (1982) have made predictions

for the buoyant jet using algebraic stress models. Many of the

ideas used in these models for calculating buoyant flows originatcd with Launder (1975, 1978). Algebraic stress models are

obtained by simplifying the convective transport equations for

Reynolds stresses and heat fluxes so they are no longer differential equations. The dynamic equation for the Reynolds

stress tensor is

2

(C-D)-r;=PJ+GJ- 3 E6J

rk

2 1

- C,e - 3 6 C2 ii 1 - (pu 3 P6J 1

30C2 -22 2 aU; 5 DJ- 3 P6;;1 - (8C2-6)k

(V

"J, + a,,,

-C3 GJ 2- 3 G6J) (t 1)

where P,•; _ —u;uk aUi

/aXk — u;Lk aU,laXk is the mechanical

production and GJ = -0g; Y17 - Og; lu^ is the buoyancy

production. The left side represents convection minus diffusional transport, the dissipation is assumed isotropic, and the

last three lines represent the closure formulation Jor

p(au;/ax;+au;/ax;) /p given by Launder et al. (1975) and

Launder (1975, 1978). Launder assumes: (1) an equilibrium

situation where convection is balanced by diffusion (C - D

= 0) and production is balanced by dissipation (P + G - e

= 0); (2) the second and third terms (third line) in the rapid

part of the pressure-velocity correlation are negligible; the

coefficient C2 is adjusted so that the first term approximates

the entire rapid part; (3) the parameter C3 is taken equal to

C2 . After applying all the assumptions

2 C,+C2-I 1-C2 k

uu;= 3C kbJ- C (PJ +G,,) (12)

t t e

where c, = 2.2 and C2 = 0.6. It should be pointed out that

in free shear flows the equilibrium condition (C - D = 0 and

P + G - E = 0) only applies in the outer portion of the flow_

Also, Zeman and Lumley (1976) found C3 = 0.3, after applying

all the constraints applicable to determining the contribution

of buoyancy to the pressure-strain correlation.

The dynamic equation for the heat flux is

a aU; (C-D) - u,u - - -ice— -ogf - - ' 

ax; ax;

- C,, ^+CL -.-,7' - 5 ia^ 'U. + C3,6gt2 (13)

axi

where the first line on the rig ht side is the production and the

second line is the closure for pat/ax;/p. Neglecting convection

and diffusion (C - D = 0) and the third term in the second

line, equation (13) becomes

1 k aT au,

u7 = Ct , e [ - u u' —ax, - 

(1 - C2,)V ax; - (1 - C31) fig

• (14)

where C„ = 3.0, C2i = 0.5, and Car = 0.5. Zeman and Lumley

(1976) show that C2, = 0.8 and C3 , = 0.2 from theoretical

considerations.

Neglecting the convection and diffusional transport in equation (5) and eliminating c, with equation (6) gives

r2 =-R E U17aX {ii)

948 1 Vol. 112, NOVEMBER 1990 130 Transactions of the ASME
O ,anr...,rr

--- ou

_ _-. — acur.^c :nns woc^

o- o

oo

o

rl

which was given by Launder (1975, 1978) and used by Hossain

and Rodi (1982).

Chen and Rodi (1975) and Chen and Chen (19 7 9) used the

differential equation (5), with e, eliminated by using equation

(6) to determine t' rather than using equation (15). In either

case R is taken to be a constant equal to 0.8 (Hossain and

Rodi, 1982; Chen and Rodi, 1975; Chen and Chen, 1979;

Launder, 1975, 1978). It is seen in Fig. 2 that the experimentally

determined value of R is much lower with an average value

across the profile of roughly 0.25. Launder (1978) cites cxperimcntal evidence for R being in the range of 0.5 to 0.8.

However, he found that the algebraic stress relations agreed

but with an experiment for a stably stratified homogeneous

shear flow with R = O.S. Hence, that value has been adopted

in the algebraic stress models- The experimental results of

Shabbir and George (1987) indicate R is much lower for strongly

buoyant flows_ When the algebraic stress model is applied to

this experiment with R = 0.8, the results arc very poor for

t2 . Therefore, in the following evaluation of the algebraic stress

model, the experimentally determined proiele for R (Fig- 2) is

used.

Equations (10), (12), and (13) represent a system of algebraic

equations that calk be solved for u,--u-j, i77, and t2 . Employing

the t::in shear layer approximation, where only gradients in

the radial direction are retained, Hossain and Rodi (1982) give

2C1 +C2 -1 1—C2k au

3 C k+ C E (

-2uu ar +2 I (16)

1 /

0.04 r/r 0.12 0.20

Fig. 7 Mean square temperature

— ecru r^rxr

—E—IC S-ESS woes

I^

o.s

.O• -08 .11 .Ib -l0

r/r

Fig. 8 Radial Reynolds stress vl

1—Cz k au

uu= 

C E (- err 

+8^ (17)

1

- 2 C2 +C2 - 1 k (18) 3 CI

u7= E ^-uU err—( 1 —C2,)^ aU +(1—C3r)FT,-] Cl,

(19)

_ -1 k—

Cl ,

aT

E .ar

(20)

k aT 2R 

E 

T7 

- err

(21.)

from which v, C„k'/E where

2(i-C2)(-I C;+C ..

r= 1+ % 3. C2

2 ) 1 k (

Cl, E d aU/er

aT/arr)

^

C

(22)

The system of equations (16)-(21) was solved to determine the

Reynolds stress and heat flux components with U, T, k, E, and

R given by the experiment. The following values for the constants were used:

C1 =2.2, 'C2 =0.6, C3=0.6

C1 ,=3.0, C2,=0.5, C3,=0.5

The value of C,,, which appears in the eddy viscosity relation

r, C„ k2/E and is given by equation (22), is roughly 0.125

and is reasonably constant across the flow. This value is considerably larger than the, value of C„ = 0.09 in the standard

k-E model. Rodi (1972), to correct for the discrepanau in the

prediction for the axisymmetric jet, developed an empirical

correction to C,,. Thc parameter C„ is replaced by (I-0.465H) C,,

whereH= 1(yE/U.)dU./dt'IQ ' ', U. is the maximum vclocity,

and yE is the distance from the centerline of the edge of the

jet. Hossain and Rodi (1982), Chen and Rodi (1975), and Chen

and Chen (1979) used this correction in their predictions for

turbulent buoyant jets. When the correction is applied, we get

approximately 0.09, which is the value of C„ for the standard

k-e model.

The question we are asking is, "given the turbulent energy,

dissipation, velocity, and temperature, does the proposed algebraic stress expression correctly predict the Reynolds stress

and heat flux components?" Figure 8 shows the radial Reynolds stress determined from equation (18). It is seen that the

predicted value ^'/k = 0.53 is a little smaller than the experimental curve in the center portion of the plume, but agrees

quite well with experiment in the outer portion. The predicted

shear stress uu, the axial heat flux i7, the radial heat flux uf,

and the mean squared temperature t2 are shown in Figs. 3, 4,

5, and . 7, respectively. Again the points arc experimental data

and the broken lines arc from the model- It is seen that the

shear stress uu and radial heat flux i7 are predicted reasonably.

However, the vertical heat flux 0 and temperature fluctuations

t2 are - predicted .poorly and-, have incorrect shapes; unlike the

experimental values they go to zero near the origin." .

Equation (21) gives t2 proportional to the radial temperature

gradient, which is zero at the ccntulinc. Then since 1'2 = 0 at

r 0, equation (19) gives u? = 0 at r = 0. In order to obtain

nonzero value: for i7 and t2 at. the centerline from the model

equations (12), (14), and (15); terms containing the axial gradicnt, i.c., aW3z and aT/az, were retained. These terms were

added to equations (19) and (21) and the system of equations

was solved again. Although the centerline -values of ut and

t2 were found to be nonzero, the predictions still decreased to

relatively small values new the centerline.

Another possibility for this behavior is.the neglect of advcaion and diffusion terms in the model. Gibson and Launder

(1976) have proposed the.following model for these-terms:

(C-D)—,= k .(P+G-E) (2-3)

(C-D)^;-^= 

2r' 

(P,-E,)+ 

2--k 

(P+G-c) (24)

where P, is the production term in the r2 equation. These were

Journal of Heat Transfer 13 1

NOVEMBER 1990, Vol. 112 1949
0 Ui .ur

r.^

— •tcsuac tnus wood

7.0 /

\\

/ O\

r /

f'o

I

o

.oe .12 .16

rig _

Fig. 9 Transport of temperature variance 12

incorporated in equations (16)-(21) and the resulting set of

nonlinear coupled algebraic . equations was solved simultancously. The results did improve the prediction for the vertical

heat flux 0 and temperature variance l but the comparison

for radial heat flux and shear stress became worse. As noted

by Gibson and Launder (1976), the above model is not good

near an axis of symmetry. This is why, by incorporating them

in the original model, no overall improvement in the prediction

is achieved.

Chen and Rodi (1975), Tamanini (1978), and Chen and Chen

(1979) use th e differential convective-transport equation (5) to

determine t in their predictions of buoyant jets- Equation (6)

was used to eliminate c,. Thus, the final form of the r z equation

becomes

ateat, 1 a k2 ail aT 1 E — U— +.V — _ -- C,---- —217— = — — ? (25)

az ar . r ar e ar ar R k

This equation was numerically solved for the temperature variance tZ with all other quantities needed to evaluate the coef-

-ficichts determined from experiments. The value of C, was

taken as 0.13. The bat agreement,. as shown in Fig. 7, was

achieved with R = 0.35. When the average experimental value

of R = 0.25 is used, the prediction peaks at about 6.0 (Ti =

0.04) as compared to the experimental value of t 2 of about 8.0

(71 0.04). When the standard value of R = 0.8 is used, tZ

is overpredicted by a factor of four.

When tt is calculated from the convective-transport equation

(5), the diffusive transport is given-by the simple gradient

closure

a^ (26)

with C, = 0.13 as given by Chen and Rodi (1980). The prediction for v1r, using experimental information to evaluate the

right-hand side of equation (26), is shown in Fig. 9. The predicted curve peaks somewhat above and toward the centerline

as compared to the data.

Finally, we ask that if the models do not depict the axial

heat flux 0 and the temperature variance P correctly, then

why do the predictions such as made by Hossain and Rodi

(1982), Chen and Rodi (1975), Tamanini (1978), and Chen and

Chen (1978) show reasonable agreement with the experiment

for the mean velocity and buoyancy? The answer to this is

that with R = 0.8 the temperature variance t from equation

(21) or (25) is too large. This makes (tic vertical heat flux

from equation (19) large enough so that the mean velocity and

buoyancy are reasonably predicted.

5 Summary and Conclusions

The experimental data on buoyant plumes were used to

evaluate various closure relations for turbulence transport. The

objective was not to propose new models, but to evaluate the

closure schema proposed by other workers for buoyancydominated flows. The closures evaluated were those used in

the k-E and algebraic stress models. The results arc summarized

below.

1 The closure relations of the k-c model compare rcasonably with experimental data, except for the axial turbulent

transport, which is drastically underpredicted. The axial heat

flux governs the production due to buoyancy in the kinetic

energy and dissipation equations and its correct prediction is

very important. This is a probable reason why the results of

Hossain and Rodi (1982) from the k-c model undci predict the

spreading rate for the plume by 10 percent even when w6syrnmetric jet corrections arc included.

2 The ratio R of the time scales, which is used to determine

the dissipation of the mean squared temperature in the algebraic stress model, was found to be considerably different from

the accepted value of R = O.S. Apparently R is not a universal

constant, but can vary from flow to flow and is influenced by

the strength of the buoyancy present. From the experimental

data on a plume it appears that R = 0.25 for strongly buoyant

flows3 The closure equations for the shear stress and radial heat

flux of the algebraic stress models also compared well with

experiment but are not better than the simple gradient closures

used in the k-e model. The axial heat flux and mean squared

temperature are predicted poorly in the central core of the flow

and had incorrect trends. This drawback could be attributed

to the assumption of local equilibrium, which resulted in the

neglect of convection and diffusion terms in the transport

equations for Reynolds stress and heat flux. However, no

substantial improvement was.achieved by keeping the secondary derivatives or by incorporating the model for the convection and diffusion terms. Therefore, the full dvnainic equations

for Reynolds stress and heat .flux with convection and diffusion

are required to predict the axial heat flux and temperature

variance properly.

Acknowledgments

This work was partially supported by the National Science

Foundation under Grants ATM-8023699 and MSM 8316833.

References

Chen. C J., and Chen. C. H., 1979 • "On Prediction and Unified Correlation

ror Decay of Vertical. Buoyant Jas," ASME JOURNAL or HEAT TaANsrEa, Vol.

101, pp. 532-437.

Chcn, J. C. and Rodi. W., 1975, "A Mathematical Model for Stntificd

Turbulent Flows and Its Application to Buoyant Jets." Proc. l61h Congress.

1AHR, Sao Paulo. Brazil, pp. 31-37.

George. W. K-, Alpert, R. L., and Tamanim, F., 1977, "Turbulence Mcasurcments in as Auisymmaric Buoyant Plume." /er. J. Hear Mass Transfer,

Vol. 20, pp. 1145-1151. .

Gibson. M. M., and Laundcr, B. 1_, 1976, "On the Calculation of Horiwntal

Non-cquilibrium Turbulent Shear Flows Under Gravitational Inrtucnee • " ASM E

JOURNAL Or HEAT Tr.ANsrstt, Vol. 98, pp. 81-87.

Hanialic, K., and Launder, B. E.. 1990, "Sensitizing the Dissipation Equations

to lrrotatiotsal Strains," ASME Journal of Fluids Engineering, Vol. 102, pp.

34-40.

Hossain, M. S., and Rodi, W., 1982, "A Turbulence Model for Buoyant

Flows and Its Application to Buoyant Jets," Turbulent Duoyant Jets and Plumes.

W. Rodi, cd., Pcrgamon. NY, pp. 121-178.

950 1 Vol. 112, NOVEMBER 1990 132 Transactions of the ASME
Launder, B. E., 1973, "On the Effects of a Gravitational Field on the Turbulent Transport of Hear and Momentum.- J. Fluid Mech., VCI. (.7 , pp, 569-

590.

Launder, B. E., 1978, '• Heat and Mass Transport.- Topics in Physics. Vol.

12: Turbulence, Springer-Vcrtag. Ncw York, pp. 231-287.

Launder, B. E-, Morse, A. P., Rodi, W., and Spalding, D. B., 1972, "The

Prediction of Free-Shear Flows—A Comparison of the Performance of Six

Turbulence Models," Proc. Lnntley Free Shear Flows Conf.., Vol. 1, NASA

SP 320.

Launder, B. E., Reese, G. J., and Rodi, W., 1975, "Progress in the Dcvclopmcnt of a Reynolds-Stress Turbulence Closure," J. Fluid Mech., Vol. 68,

pp. 537-566.

Launder, B. E., and Spalding, D. B., 1974, "The Numerical Computation

of Turbulent Flow, ­ Comp. Melh. in Appl. Mech. and Engr.. Vol. 3, p. 269.

Lumley, J. L., 1978, "Computational Modeling of Turbulent Flows," Advances in Applied Mechanics, Vol. 18, pp. 123-176.

Pope. S. B., 1972, "An Explanation of the Turbulcn; Round-Jct/Plane-Jet

Anomuy," AlAA J., Vol. 16, pp. 279-281.

Rodi, W., 1972, "Thc Prediaic•n of Free Turbulent Boundary Layert by Use

of a Two-Equation Model of Tutbulencc." Ph.D. Thesis, Univ. of London,

United Kingdom.

Shabbir, A., 1987, "An Experimcntal Study of an Axisymmaric Turbulcnt

Buoyant Plume and Evaluation of Closure Hypotheses,- Ph.D. Diucrtation,

Univ. at Buffalo, SUNY.

Shabbir, A., and George, W. K., 1987, "Energy Balartce Mcasurcmenu in

an A itymmetric Buoyant Plume. • Sixth Symposium on Turbulent Shear Flows.

Toulouse, pp. 9-3-1 to 9-3-6 (also submitted to J. Fluid Meeh.).

Tamanini, F., 1978, "The Effect of Buoyancy on the Turbulence Structure

of Vertical Round Jcts." ASME J otnt.NAL or HEAT Tx_txsyLR, Vol. 100; pp.

659-664.

Zcman, O., and Lumley, J. L., 1976, "Modeling Buoyancy Driven Mixed

Layers, ­ J. Arm. Sci., Vol. 33, p. 1988.

RPa dPr; of

The e ournal - 'of Heat Transfer 

Will. Be- Interested IA-. 

HTD-Vol.-138 -... _......._,. -..

Heat-Transfer in Turbulent Flow.

Editors: R.S. Amano, M.E. Crawford N.K: Anand..

;:fTopics covered include. fundamental researcli^o' •turbu{ence,iti,heat transfer. pror:esses,. boundary, ;_;n:

.layer flows, temperature _turbulence spectrum, turbulence-modeling; and applicationsao:"_ =,< c• :^:: ;:;

heat exchangers; and gas turbines.^-

_ - _-. - _ . _.: ....:. 1. -.. : JI: : al_-^:l^,lf^aT ^Tt.Z _ .iri,•• •- •-+ ..'..-t f;._.t_

_ • T_. .. :. . ^r5`53L'. fYr ;^%fIiT1Z? ;47:• ':ji ' ii; :>o^fa=^'^•jT'.i if ".rl --it .-

1990. - . -- .0rder No. H00591.' . = ISBN No. 0-7918-0483 L 6 .,2 G pp.-.

 ^{ -:^ .w=' ? =

$30 List/$15 ASME Members ' ...- . _ _........_

_ .- _ .. __ `—^^s _,.'J1 J•=i::L)^I?t' i2• ,t'c- cam 3. :% - -. -.-^

a' : _ ' •J - r. .T7 •-- ..F1^t •u:li..sr::^'Jl:i

To order, write Order Department, 22 Law Dnve, Box-2300; Fairfield,^NJ 07007-2300 _''?Z„ ::^ ^;;^; j;X_,:

or call 1--800-THE-ASME (843 ,2763) or FAX 1-201-882-1717 %. "'`F

"c • '-..,,'^:a,^F.^' _ _._.. - ;.i, 1 .. .. .... ^. 1 .^. `ii ^ r ^ti:.^.v-:'^cYLr-.<^ ^N-+=L:7',i^i.4^.1C^i11r^ ,l'..^^ [s'^^.'_G'iift. ^ti

133

Journa l of Heat Transfer NOVEMBER 1990, Vol. 112 / 951
SOURCE: CMOTT-91-06

On the Basic Equations for the Second-order Modeling of

Compressible Turbulence

W. W. Liou and TAL Shih

Center for Modeling of Turbulence and Transition

ICOMP/NASA Lewis Research, Center

Cleveland, OH 44135

Abstract

Equations for the mean and the turbulence quantities of compressible turbulent

flows are derived in this report. Both the conventional Reynolds average and the

mass-weighted Favre average were employed to decompose the flow variable into a

mean and a turbulent quantities. These equations are to be used later in developing

second-order Reynolds stress models for high-speed compressible flows. A few recent

advances in modeling some of the terms in the equation due to compressibility effects

are also summarized.

134
AIAA-91-1801-CP

UNSTABLE VISCOUS WALL MODES IN ROTATING PIPE FLOW

Z. Yang'

Center for Modeling of Turbulence and Transition

ICOMP, NASA Lewis Research Center

Cleveland, OH 44135

S. Leibovicht

Sibley School of Mechanical and Aerospace Engineering

Cornell University

Ithaca, NY 14853

Abstract

Linear stability of flow in rotating pipe is studied.

These flows depend on two parameters, which can be

taken as the anal Reynolds number Re and the rotating rate, q. In the region of Re » 1 and q = 0(1), the

most unstable modes are concentrated near the pipe

wall, the so-called "wall" modes. These wall modes

are found to satisfy a simpler set of equations containing two parameters rather than four parameters as in

the full linear stability problem. The set of equations

is solved numerically and asymptotically over a wide

range of the parameters. In the limit of Re —+ oo,

the eigenvalue goes to the inviscid limit. The eigenfunction shows a two layer structure. It reaches the

inviscid limit over the main part of the domain, while

near the wall of the pipe, the eigenfunction is represented by a viscous solution of boundary layer type.

1. Introduction

Swirling flow is common ir, nature and technology.

Fully-developed flow in a rotating pipe is an exact

solution of the Navier-Stokes equations, and is the

simplest available model of swirling flows. Swirling

flows are known to be subject to instability, and the

question of stability of flow in rotating pipe has consequently attracted a reasonable amount of attention.

Pedleyr showed that these flows are unstable to inviscid non-axisymmetric perturbations when the rotation is fast (in a sense which will made definite).

Stability to inviscid perturbations in the finite rotation rate region was studied numerically by Marlowe'-.

Later, Maslowe and Stewartson 3 extended this work

in a significant way and established by asymptotic

methods that the dominant unstable modes are wall

'Research Associate

t Profce or

modes (that is, modes of motion concentrated asymptotically close to the wall) in the limit of large azimuthal wave number.

For viscous perturbations, Pedley 4 , and simultaneously Jesoph and Carr is found that the critical

Reynolds number at which the perturbation is neutral in the limit of fast rotation. Cotton and Salwene

carried out comprehensive computations in search of

neutral stability curves and they discovered that the

neutral modes are center modes when Reynolds number is large. Center modes in rotating pipe flow were

later analyzed asymptotically by Stewartson, Ng and

Brown7 , and these authors speculated that center

modes dominate for large Reynolds number.

In this study, we investigate the effect of viscosity on the inviscid wall modes found by Maslowe and

Stewartson3 when the Reynolds number is large but

finite. We find the proper scaling for Reynolds number in order for viscous wall mode to exist and derive

simplified governing equations for them. These equations contain only two parameters instead of four parameters in the full linear stability problem. The wall

mode equations are then solved both numerically and

asymptotically.

The plan of this study is as follow: Linear stability analysis is formulated in section 2, where it is

shown numerically that the most unstable modes are

wall modes. The governing equations for viscous wall

modes are derived in section 3_ Numerical solutions

of the viscous wall mode equations are presented in

section 4. An asymptotic analysis for the viscous wall

modes equations is carried out in section 5, and section 6 concludes the paper.

2. Linear Stabilitv Formulation

If length is nondimensionalized by the radius of the

pipe L and velocity by the the axial velocity at the

axis U, the laminar base flow in a pipe rotating with

Copyright©1991 by the authors. Published by

the American Institute of Aeronautics and

Astronautics, Inc. with permission. 135
angular velocity Q is then described in a cylindrical where

(x, B, r) coordinate system by,

-t= k(1—r2)+mq—w

U = (1 — r 2 , q r, 0)

where the inverse Rossby number

_ nL

q U

(1)

The above equations for the perturbations are supplemented by the following boundary conditions. On

the wall of the pipe, the no slip boundary conditions

(2) require

measures the relative strength of rotation. This

nondimensionalizaLion also defines a Reynolds number

Re = UL (3)

I/

where v is the kinematic viscosity of the fluid.

Linear stability analysis concerns the stability of

the motion subject to infinitesimal perturbations.

Since the linear stability problem with base flow given

by (1) is then separable in the x and 0 directions, the

perturbation field may be written in the normal mode

form. i.e.

A (u(r), v(r), w(r), p(r)]exp [i(kx + mB — wt)]

where A is an arbitrary constant, k is the axial wave

number of the perturbation and m is the azimuthal

wave number of the perturbation. w is the complex

frequency, with its real part being the frequency and

the imaginary part being the growth rate. Without

loss of generality, m is taken as positive, k is any real

number, and w is to be found. If Im(w) is positive,

the flow is linearly unstable.

The Navier-Stokes equations linearized about the

base flow, and the equation of continuity are of the

following form in the cylindrical coordinate system

used,

i-fu — 2rw + ikp =

1 a, u 1 au m2 u — k 2

Re (ar e + r r2 a 

u

i-ry + 2qw + imp —r

1 1 av m2 + 1 2iw

2Re 

(a2v

ar 2 + r ar r2 v +m — k r2 v

i-yw — 2qv + ap

a =

1 (a'w 1 aw m2 + 1 2imv

art + r aT r2 w — 

k-w r2 

im aw w iku+ —v - —ar +— =0 (4)

r r

U(I) = v(1) = w(1) = 0 (5)

At the center of the pipe, the perturbation must

satisfy the following conditions to ensure that to be

single-valued8

for m=0

u'(0) = v(0) = w(0) = p'(0) = 0,

for IM! =1

u(0) = v(0) + im w(0) = p(0) = 0,

for m otherwise

u(0) = v(0) = w(0) = p(0) = 0. (6)

The above ordinary differential equations (4) and

the boundary conditions (5) and (6) form an eigenvalue problem with w as the eigenvalue. Nontrivial

solution exists only when w takes some specific values

given by

w = w(Re, q; m, k)

If Im(w) is less than zero, the flow is stable; if Im(w)

is positive, the flow is said to be linearly unstable. We

are interested in the regions where the instability occurs. Earlier studies show, and our results confirm,

that instability occurs only for mk < 0. In the following, we will take m positive, and k negative.

The linear stability problem was studied extensively by numerical means by Cotton and Salwene

and by Yang'. Cotton and Salwen searched in the

parameter space (Re, q; m, k) for the neutral modes

(modes with zero wi ), and found that the neutral

modes are center modes, with most of the uontrivial

activities confined near the center of the pipe. Yang's

study concentrated on the most unstable modes, and

we shall briefly describe the most unstable modes that

he found in the region Re » 1 and q = 0(1).

Table 1 gives the information for the most unstable modes at Re = 10000 and some different values

of q's. As q is increased, the the most unstable mode

has larger values of both the azimuthal and axial wave

numbers. In addition, we find for q = 0(1), although

the real part and the imaginary part of the eigenvalue have quite different size, the difference of the

real part of the eigenvalue from mq is of the same

order as that of the imaginary part. These findings

136
Table 1: The most unstable modes for different q's

with Re = 10000.

m: the azimuthal wave number,

k: the axial wave number,

w: the eigenvalue of the most unstable mode

q m k W

0.5 3 -0.50 (1.25, 0.177)

1.0 5 -0.81 (4.67, 0.296)

1.5 7 -1.04 (10.13, 0.379)

2.0 9 -1.20 (17.61, 0.439)

3.0 11 -1.22 (32.63, 0.518)

suggest the viscous wall mode scalings studied in the

next section. For q = 3, the most unstable mode has

the wave numbers m = 11, k = —1.22. The eigenfunction corresponding to this most unstable mode is

shown in Fig 1. The eigenfunction is normalized such

that the maximum axial velocity is 1. The real part

of the eigenfunction is drawn in solid lines, and the

imaginary part is drawn in dotted lines. The nontrivial behavior of this eigenfunction takes place in a

thin region near the wall of the pipe. This behavior

justifies the "wall mode" terminology used, and will

be the subject of further study in next section.

3. The viscous wall mode equation

From the numerical computations in the last section, it is clear that in the region of Re » 1 and

q = 0(1), the dominant modes are given by the

asymptotic wall modes, a type of modes with large

azimuthal wave number and with nontrivial behavior

concentrated near the wall.

When azimuthal wave number m is large, there

could be another type of mode for general swirling

flows, the ring mode, as demonstrated by Leibovich

and Stewartson 10 . In the case of rotating pipe flow,

ring modes do not exist. Because of the existence of

the pipe wall, wall modes characterize the behavior

of perturbation with large azimuthal wave number.

Maslowe and Stewartson 3 studied the stability of

inviscid pipe flow to perturbations of very large azimuthal wave number, and established that the prevailing modes are the wall modes in the limit of

M oo. Stewartson 11 analyzed the effect of the

viscosity on ring modes and found that viscous effects come into play in higher order terms because

the inviscid ring mode solution can satisfy the exact

boundary conditions of the problem.

We study the wall mode when the viscosity is taken

into account. Since the inviscid wall mode is close to

the wall, where the no-slip condition is required but

,..

u

(C)

Figure 1: Eigenfunction for Re = 10000, q = 3, m =

11,k = —1.22. (a) axial velocity, (b) azimuthal velocity, (c) radial velocity.

is not satisfied by the inviscid solution, it is expected

that the viscous effect might be important in this

case, at least in certain regions .

We are going to study the linear stability problem

for large azimuthal wave number and large Reynolds

number, i.e. both m and Re are large. In this case,

the proper wall mode scalings are

r = 1— 77

(.m2 + k2)1/2

W = iw

2p

P =

m

(.m2 + k2)3/2 Re = —Re

2k

_ 2k /

(.J mq (m2 +k2)1/2 l7)

The scaling for the radial variable means that the behavior of the wall mode is confined to a small distance

to the wall comparable with the wave length of the

perturbation. The scaling for the Reynolds number

giv , the balance between the. viscous term and the

inertia term. The form of the scaling for the complex

137
frequency is suggested by the numerical computations

of the full linear stability problem. In above, a minus sign is introduced in places where k appears, for

we know from the full linear stability problem that

instability occurs only when m k < 0.

Upon substituting the above expressions into the

linearized Navier-Stokes equations and the equation

of continuity for the perturbation, and dropping the

terms of order m- 1 and smaller, we find that w should

satisfy the following single equation, after u, v, p are

eliminated.

LDLDw - LLw + LDw - Quo = 0 (8)

D = d

T77

L Re(D2-1)-(,7

Q = - 4(m + qk)

(9) k

which emerges as one of the independent parameters

for the viscous wall mode.

The boundary conditions are

w=Dw=D2 Lw-Lw=0 at r7=0 (10)

and

nonlinear eigenvalue problem. This nonlinear eigenvalue problem is changed to a system of linear eigenvalue equations by letting

Y = Lw (12)

In Y, w, the governing equations are

L(D2 - 1)Y + 2DY - (Q - 2)w = 0

Y - Lw = 0 (13)

and the boundary conditions are

w=Dv,=D2 Y-Y=0 atr7=0 (14)

w= Div =D2Y-Y=0 asr7-»oo (15)

We employed a spectral method with Chebyshev

polynomials as the basis functions to solve equations

(13) - (15). Because the domain of definition extends

to infinity, while Chebyshev polynomials are only defined over [-1, 1], the method of domain trur:cation

was used to numerically truncate the domain of definition from [0, oo) to [0, b] and the boundary conditions at infinity are replaced by

w=Dw=D2 Y-Y=Oat 77=b (16)

In the spectral method used, we write

N

;-r

where

and

w=Dw=D2 Lw-Lv,= 0as r7 -+oo (11)

Thus, we have established the formulation for the

viscous wall modes. The governing equation takes

a simpler form compared with the full problem, and

the number of the independent parameters is reduced

from four to two. Of these two parameters, Re and Q,

Q measures the effect of rotation on the wall mode

and Re measures the effect of viscosity on the wall

mode. The solutions of above wall mode problem

will give:

(D = CD (Re ' Q)

If Im(LD) is less than zero, the flow is linearly stable;

if Im(C,)) is positive, the flow is said to be linearly

unstable.

4. Numerical Solution of Wall Mode Equations

In general, solutions of the wall mode equation have

to be found numerically. Because the eigenvalue enters quadratically, this wall mode equation gives a

N

Y = E Y,"T-i(y) (17)

;-r

where y is related to r1 by

y = 2r7/b - 1.

Thus the domain of definition for y is [-1,11. The

reduction from the differential equations to a set of

algebraic equations is made by a Galerkin-Tau projection —i.e. the Galerkin method is used to project the

equations while the Tau method is used to enforce the

boundary conditions on the spectral representations

of the perturbation field. The Galerkin-Tau projection results in an set of algebraic equations, which are

of the form of generalized eigenvalue problem with

complex matrices.

Two methods are used in this study to find the

eigenvalues and the eigenvectors of this complex generalized eigenvalue problem. One method uses the

IMSL subroutine EIGZC which uses the QZ transformation to find all the eigenvalues and, optionally, all

138
the eigenvectors. The other method used is an inverse

power iteration for the generalized eigenvalue problem developed by Kribus", which finds the eigenvalue closest to the initial guess and its corresponding

eigeavector. The inverse power iteration is faster, so

it is used whenever a good guess is available. The QZ

is used to provide the starting values for the Inverse

Power Iteration. It is also used when the phenomenon

of mode jumping is suspected to occur.

Since the eigenfunctions decay exponentially as

n — oo, b = 10 was found sufficient for most of the

calculations. The number of terms in the Chebyshev

representation, N, varies with Re. For order one Re,

we must take N = 60, and b = 10 to achieve three

digit accuracy in eigenvalue, and 0(10 -4 ) accuracy

for the eigenfunction. But when Re is large, the solution shows a behavior of boundary type for 71 near

zero, i.e. near the wall of the pipe, and a large value

of N is needed to resolved this region. The largest N

needed for the parameter range covered here is 115.

There is a symmetry in the eigenvalue problem due

to the replacement of the boundary conditions at infinity by the conditions at 77 = b. The solution is

invariant under

77 b — n

w .--. w'

Y Y'

where the star means the complex conjugate. This

symmetry signifies that for a given solution, its image about r) = b/2 is also a solution of this equation

with the same growth rate. Apparently, this solution

is a spurious mode in the sense that it is a solution

of the differential equation after enforcing the boundary condition at finite b rather than the solution of

the original differential equation defined over the infinite domain. This symmetry property can be used

to check the resolution of the numerical solution of

the algebraic problem, which should also exhibit this

symmetry.

Extensive computations were carried out in the

(Re, Q) plane. Fig 2 shows the imaginary part of the

eigenvalue C), which is proportional to the growth rate

of the perturbation, vs. Re for some different values

of Q. The eigenvalues shown are for Q = 5, 10, 20, 30,

40, 50 respectively, although numerical computations

were carried out for a larger range of Q. The real

part of the eigenvalue is shown in Fig 3. The growth

rate increases as Q is increased, which r ,!ans that rotation helps perturbations to extract energy from the

,.0

... r / ii

14

Figure 2: The imaginary part of the wall mode eigenvalues vs. Re for Q = 5, 10, 20, 30, 40, 50.

base flow. This result is in agreement with that of

Pedleyt , who found that the maximum growth rate

is 2.0 and is reached in the fast rotation limit. This is

exactly the same as the upper bound for the growth

rate for flow in the rotating pipe, as shown by Joseph

and Carmis.

The growth rate increases with Re. As Re gets

large, the growth rate will increase. The eigenvalues

seem to change smoothly as ' Re — oo and approach

to the results found by Maslowe and Stewartson from

their inviscid analysis. The limit of Re — co will be

analyzed asymptotically in the next section, and a

comparison of the results with the inviscid case will

be made.

Fig 4 shows the eigenfunction for Re = 2, and Q =

5. The eigenfunction plotted is normalized such that

its maximum modulus is 1 and its phase at the position of maximum modulus is 0. As in the eigenfunction plotted in Fig 1, the real part of the eigenfunction

is drawn in solid line while the imaginary part of the

eigenvalue is drawn in dotted line. To see the effect

of increasing Q on the eigenfunction, we show in Fig

5 the eigenfunction for Re = 2 and Q = 20. As Q

increases, the eigenfunction is pushed outward, but

this response is not very sensitive to Q. To see the

effect of changing Re, Fig 6 shows the eigenfunction

for Re = 100 and Q = 5. As Re increases, the eigenfunctions are pushed towards n = 0, i.e. toward the

wall. The limit of Re oo will be the further studied

in the next section.

5. The Limit of Re oo

139
I p x,^^

7

7

.., ]-•

/.I ].l 11,1 11.1 L.l /., ]., ..• ..l 11 1/./

Figure 3: The real part of the wall mode eigenvalues Figure 5: WaII mode eigenfunction, Re = 2, Q = 20

vs. Re for Q = 5, 10, 20, 30, 40, 50.

Figure 4: Wall mode eigenfunction, Re = 2, Q = 5 Figure 6: Wall mode eigenfunction, Re = 100, Q = 5.

140
Maslowe and Stewartson 3 carried out the wall

mode analysis for m » 1 for the inviscid case. One

of the purposes of our viscous wall mode analysis is

to see if the results of their inviscid analysis is the

limit of the viscous analysis for large Reynolds number. We use a perturbation technique to address this

question in this section.

For Re » 1, to can be expanded formally by taking

Re the small parameter, we write:

w = Woo ( 17) + lie- 1Wi(11) + ...

To leading order, the equation for Wo is:

D2 Wo - (1 + Q )Wo = 0 (18)

( 11 +1;1)2

with the following boundary conditions:

Wo(0) = 0 (19)

In addition, the solution is required to match the

outer solution in the matched asymptotic sense.

The general solutions of the inner problem in the

leading order can be found, and they are

Wo

o

 = Ql exp(A1C) +,3.2 C eZp (A 1C) +03 exP(A2()

+Q4 C eZp(A2t) +As +Qa t (23)

where

A1,2 = x/ )1/2

and 131, 02, 03, 04, )3s, /3a

1

are constants to be determined by the boundary conditions and the matchinf

conditions. CD is given by the outer solution.

The boundary conditions at C = 0 require thatt.

01 +)33 + 05 = 0

W0'('7 - -) - 0 (20) 31 A , +Q2 +03A3 +Q4 +Ot = 0

This poses the same eigenvalue problem that was 3 3

studied by Maslowe and Stewartson, thus their solu- 2a1 Q2 + 2^ '2

= 0

tions (both the eigenvalue and the eigenfunction) may As — oo the inner solution must also match the

be viewed, as might have been expected, as the first outer solution. Near 77 = 0, the outer solution is, to

term in a formal outer expansion in inverse powers of leading order,

Re.

The outer solution thus found can satisfy all the Wo = a 77

boundary conditions at infinity since the outer solution decays exponentially as 71 oo. But not all the where a = DWo(0) ,f 0.

boundary conditions at n = 0 can be satisfied. For To carry out the matching, we need to know the

example, DWo (0) # 0. Thus, another (inner) solu- behavior of Wo for large C. The behaviors of the

tion of boundary type near 77 = 0 is needed. exponential terms are determined by the sign of the

The inner variable is found to be: real parts of a l , A 2 . Since a1i A 2 are the square roots

1/2 of (i,), a l is the negative of a 2 . In this study, we

= 711e take

The inner expansion is assumed to be:

to = Red ( Wo(() + Re -1/2 Wi (C) + ...)

where d is to be determined.

After substituted into the governing equation, the

equation to leading order is found to be:

( D2 +^) 2 D2 W0 = 0 (21)

where

d

D = d(_

The boundary conditions for the inner solution are:

Wa(0) = 0

DWQ(0) = 0

(D2 + 4))" Wo(i0 ) = 0 (22)

Re(a1) < 0

Re(A 2 ) > 0

As -+ oo, the matching of the exponentially growing terms gives:

Q3 = 0

Q4 = 0

The exponentially small terms are immaterial, and

the matching of the algebraic terms gives

Red Qa = a rl

which gives

1

d=--

2

(24)

141
Table 2: Eigenvalue of viscous wall mode.

Q eigenvalue (Re = 100)

5.0 (-0.14917D+01, 0.72752D+00)

10.0 (-0.18163D+01, 0.15096D+01)

20.0 (-0.21617D+01, 0.26668D+01)

30.0 (-0.23766D+01, 0.35755D+01)

40.0 (-0.25358D+01, 0.43505D+01)

50.0 (-0.26633D+01, 0.503841)+01)

Table 3: Eigenvalue of inviscid wall mode.

Q eigenvalue (Re = oo)

5.0 (-0.14142D+01, 0.68305D+00)

10.0 (-0.17476D+01, 0.148751)+01)

20.0 (-0.21038D+01, 0.26574D+01)

30.0 (-0.23247D+01, 0.35710D+01)

40.0 (-0.24877D+01, 0.43487D+01)

50.0 (-0.26181D+01, 0.50383D+01)

and

Qe=o

All the other constants can be determined, yielding

Q1 = -c,/A1

'3, = 0

Qs = a/A1

The inner solution, to the leading order, is

Wo(() +

a1 

(1 - exp(A1())) (25)

Thus, we have the following composite solution to the

leading order:

Wo = Wo ('1) + WoM - -'I (26)

Thus, one sees that indeed, the solutions for the

most unstable modes found by Maslowe and Stewartson are the correct limit when Re oo, except in

a thin layer near the wall of the pipe where the flow

field is described by a viscous layer of the boundary

layer type. The eigenvalues found are the same as for

the inviscid case.

The above asymptotic analysis is confirmed by the

numerical calculation. In Table 2, we present the

eigenvalues for Re = 100 and in Table 3, we present

the eigenvalues from the inviscid calculation fora few

values of Q. It is seen that differences of the eigenvalues in these two cases are small. In Fig 7, we show

the eigenfunction for Q = 5.0 from the inviscid calculation. It is seen that the general shape agrees with

the viscous calculation presented in Fig 6. To see

the existence of a thin viscous layer near the wall, we

present in Fig 8a a blow-up of Fig 6, and in Fig 8b

a blow up of Fig 7 near the wall. It is seen that flow

fields near the wall are different, the viscous solution

has a zero slope while the solution from the inviscid

equation has a non-zero slope.

,..

n

..

6. Discussion

We have examined the linear stability of rotating I

pipe flow to perturbations of large azimuthal wave

number. It is found that when the azimuthal wave

number is large, the nontrivial behaviors are concentrated near the wall of the pipe, so the prevailing

variations are manifested as wall modes. The equa- Figure 7: WaIl mode eigenfunction for the inviscid

tions governing wall modes are found to contain two case, Q = 5.

parameters Re and Q, which measure the Reynolds

number and swirl, respectively.

142
expected that a local viscous critical layer based on

( a i our viscous wall mode formulation would get rid of

this singularity.

Acknowledgement

q

This work was supported by the Air Force Office

of Scientific Research under contract AFOSR-89-0346

monitored by Dr. L. Sakell.

References

L.

(b)

,J

7

Figure 8: Comparison of the wall mode eigenfunctions near the wall. (a) the viscous solution at Re =

100, (b) the inviscid solution.

For large azimuthal wave number, the wall mode

is a distance 0(m-1 ) away from the wall. Viscous effects are confined to a layer with thickness

of O(Re -113 ) near the wall. The case treated here

is when those two layers are of the same order, i.e.

Re = O(m2 ). When Re/m2 > 1, we would expect the solution to be mainly inviscid, as found by

Maslowe and Stewartson, except in a thin viscous

layer near the wall. This is indeed the case as shown

both numerically and asymptotically.

The region for the rotation rate considered is q

0(1). For large value of q, which corresponds to the

fast rotation case, Pedley4 shows that as Reynolds

number is increased, the dominant mode is taken by

perturbations of larger and larger values of m, the

azimuthal wave number. In the limit of Re oo, the

azimuthal wave number for the dominant modes will

also go to infinity, and these modes could be viewed

as wall modes. Thus, our work here can be viewed

as an extension of Pedley's work" to the region of

q = O(1).

Left unexamined is the case when the inviscid solution has a singularity, in which case the form of

the expansion we assumed would break down. In the

study of Maslowe and Stewartson, the leading mode

is a neutral mode when Q = 2. For Q > 2, the

leading mode is unstable, but there are still neutral

modes. The neutral modes have a 1/(r - ro) singularity, where ro is the location of the singularity. It is

1 T.J. Pedley. "On the instability of rapidly rotating shear flows to non-axisym metric disturbances."

J. Fluid Mech., 31:603-607, 1968.

' S.A. Maslowe. "Instability of rigidly rotating

flows to non-ausymmetric disturbances." J. Fluid

Mech., 74:303-317, 1974.

3 S.A. Maslowe and K_ Stewartson. "On the linear

inviscid stability of rotating Poiseuille flow." Phys.

Fluids, 25:1517-1523, 1982.

4 T.J. Pedley. "On the instability of viscous flow in

a rapidly rotating pipe." J. Fluid Mech., 35:97-115,

1969.

s D.D. Joseph and S. Carmi. "Stability of

Poiseuille flow in pipes, annuli, and channels." Quart.

Appl. Math., 26:575-599, 1969.

e F.W. Cotton and H. Salwen. "Linear stability

of rotating Hagen-Poiseuille flow." J. Fluid Mech.,

108:101-125, 1981.

T K. Stewartson, T.W. Ng, and S.N. Brown.

"Viscous center modes in the stability of swirling

Poiseuille flow." Phil. Trans- Roy., 324:473-512,

1988.

a G.K. Batchelor and A.E. Gill. "Analysis of

the stability of ansymmetric jets" J. Fluid Mech.,

14:529-551, 1962.

9 Z. Yang. Two theoretical studies of stability of

swirling flows. PhD Thesis, Cornell University, 1990.

ro S. Leibovich and K. Stewartson. "A sufficient

condition for the instability of columnar vortices." J.

Fluid Mech., 126:335-356, 1983.

rr K. Stewartson. "The stability of swirling flow

at large Reynolds number when subject to disturbances of large azimuthal wavenumber." Phys. Fluids, 25:1953-1957, 1982.

" A. Kribus. "Computation of leading eigenspaces

for generalized eigenvalue problems." Trans. of the

7th Army Conf. on Appl. Math and Comp., 1990.

143
Development of A Nei: Flux Splitting Scheme

Meng-Sing Liou' and Christopher J. Steffen, Jr.t

NASA Lewis Research Center

Cleveland, Ohio

Maximizing both accuracy and efficiency has been the

primary objective in designing a numerical algorithm

for computational fluid dynamics (CFD). This is especially important for solution of complex 3D problems

which often involve Navier-Stokes equations with turbulence modeling and chemical species equations. Upwind

schemes have been well received for both their capability of resolving discontinuities and their sound theoretical basis in characteristic theory for hyperbolic systems. Several flux splitting schemes, notably by StegerWarming, Van Leer, Osber, and Roe, have been tested

and discussed extensively in the past decade.

However, several inherent shortcomings exist in each

of these schemes. For example, while the Van Leer

scheme is simple, taking only O(n) operations, n being

the number of equations, and yields accurate solutions

for inviscid problems, it suffers inaccuracy for predicting

velocity and temperature fields in the viscous problems.

The Roe scheme, commonly accepted as the most accurate scheme available currently, however is a great deal

more complex and costly due to the matrix operation,

requiring O(n') operations. Moreover, the extension to

the chemically reacting flows renders no unique way of

defining the `Roe-averaged' states. The Osher scheme

has the smooth property and has recently been generalized by Suresh and Liou to deal with chemically reacting

flows. But the determination of the intermediate states

also requires O(n') operations. Thus it is logical to ask

whether there is room in the universe of upwind schemes

for improvement to arrive at a simple (O(n) operations)

and accurate scheme for a wide range of problems.

In this paper, we summarize recent successes of a

new splitting scheme for some model aerodynamic problems where Van Leer and Roe schemes failed. The new

scheme is based on a rather different idea of splitting in

which the convective and pressure terms are separated

and treated differently in accordance with the underlying physical intuitions. We propose an appropriately

defined cell-face advection Mach number using values

from the two straddling cells via associated characteristic speeds. This interface Mach number is then used to

determine the upwind extrapolation for the convective

quantities. Next the pressure splitting is weighted using polynomial expansions of the characteristic speeds.

Thus, the name of the present scheme is properly coined

as Advection Upstream Splitting Method (AUSM). The

scheme is remarkably simple and yet its accuracy in the

present study rivals and in some cases surpasses the Roe

scheme in the Euler and Navier-Stokes solutions at considerably reduced computational effort. The detailed

formulation of the scheme is net sbown here due to space

limitation. However it will appear elsewhere.

The calculation of the hypersonic conical flow demonstrates the accuracy of the splittings in resolving the

flow in the presence of strong gradients. The temperature and pressure profiles of the first order results are

shown in Figs. 1 (a) and (b). The Van Leer splitting is

seen to produce a thicker boundary layer which in turn

further displaces the shock wave. Both AUSM and Roe

solutions are in excellent agreement.

The second series of tests involve the 2D inviscid flow

over a NACA 0012 airfoil. The results, not included

here, demonstrate that the level of entropy generation

at the stagnation point is about three times smaller than

the Roe solution.

In the third case we calculate a series of supersonic

flows over a circular cylinder. The Roe splitting in all

conditions and grids tested yields anomalous solutions

(sometimes referred to as the carbuncle phenomenon),

which may appear as nonsymmetric, protuberant, or indented contours, see Fig. 2. The mode of these nonphysical solutions appears to be sensitive to changes in

Mach number or grid. The AUSM, however, gives expected solutions in all calculations.

The fourth test deals with a 2D shock wave/laminar

boundary-layer interaction. In Fig. 3, the AUSM is seen

to give excellent agreement with the data especially in

the reattachment region, which has been in defiance of

many previous calculations in the literature. Also the

oblique shock appears to be more tightly captured by

the AUSM.

In summary, it appears that the new splitting scheme,

AUSM, has delivered the promise by improving the accuracy as well as efficiency significantly. As usual, the

final judgment will be decided via many more tests and

further modifications of the scheme.

Senior Scientist, Internal Fluid Mechanics Division. Member AL4A

Aerospace Engineer, Computational Fluid Dynamics Branch. Member AIAA

144
15

fig. 1(b) 4

1

<n

Q) 13

4D

-d 12

°P

11

10L

- Val; Leer

Roe

o AUSM

M = 7.95

Pr = 1.0

O

CO

4

o_

x 2

U 0

2

U

4

0

r-4 2

X

U 0

-2

15

fig. 1(a) 14

Q' 1.3 a^

s`.

b 12

11

10

----------

Va Ir'

Ro

o AU I

M =

Pr I

I

i

0 1 2 3 4 5 0 4 8 12 1 J6

(Twall Tm)th6ory

Pf P- T/T^

Fig. 1 Conic flow with M. = 7.95, Re = 4.2 x 10 5 , and half cone angle = 10°: (a) pressure, and (b) temperature

distributions. ,00 ------ -

'00 -

1 0 0000

5

/ S 2 0.000

0.71 0.71 3 0 8000

16 / 1 (

/ a 1

^

2D00

fig. 2 a) \4 fig. 2(b) 0.3 - , 1 6 OODO

b' ( 0.43 3 15 3^ ] 2.000

14=?", i.

2

I /^- - // 8 28000

0 t. 0.,4 I 2 9 32000

y .^ 21 10 36000

16 _ C '^ n •0000

-0.14 -0.14 = ^^^ ,2 ..000

3 4 8000

14 \ 13 1 4 5 2000

-013 ^ ^/ 043 15 `\ ^^ ^ 15 °6000

\ / 5 16 6 1000

17 E 4000

-0 71 J.71 u

t8 fi800^

. ^\ t972000

I \ _. 20 76000

-100 +-T---r-a--- T--•- ''00

-1.00 -0.71 -0.43 -0.14 014 0.43 0.71 1.00 -1.00 r-0 71 -0.43 -0 1. 0.14 0.43 0.71 t 00

+ aza o-a.^s

Fig. 2 Mach contours for a Mach 6 inviscid flow over a circular cylinder: (a) the AUSM solution, (b) the Roe

solution displaying a protuberrant, two-shock contours.

1.300 -

_-- 1.300

1.114 1 Ill

fi g. 3(a) 

0.929 4 1414 fig 3(b) 0929

0.743 5. 13 07-3

Y-'u'n ^^^ o y-ans

12

X19

0.557 14 \

12 0.557

20

13

0.371 14 0.371

20 11 ,0

13

0.186 0.186

O 4

19 9 5/ 9

C

n ry,n 0.000

0.00 031 0.63 041 126 1.57 1.89 2.20 0.00 0.31 0.63 0.94 1.26 1.57 1.89 2.20

.aos ..

n

. a=K

1 0.5000

2 1.0000

3 12500

41.5000

5 1 7000

6 1 7200

7 1.7400

8 ,7600

g 17600

10 '8000

11 1.8200

12 1.8400

13 1.8600

14 1.8800

15 1.9000

16 1.9200

17 1 9400

18 1.9600

19 1 9800

20 2.0000

0.0 1.0 

X

2.0 0.0 1.0 X 2.0

Fig. 3 Mach contours and comparison of skin friction with data for M. = 2.0, Re = 2.96 x 10 5 [Hakkinen et al]:

(a) the AUSM solution, and (b) the Roe solution. 145
High-Order Polynomial Expansions(HOPE) for Flux- Vector Splitting

Meng-Sing Liout and Chris J. Steffen, Jr.*

NASA Lewis Research Center

Cleveland, OH 44135, U.S.A.

Summary

The Van Leer flux splitting is known to produce excessive numerical dissipation for

Navier-Stokes calculations. One example is the incorrect prediction of boundary-layer

profiles. We attempt in this paper to remedy this deficiency by introducing a higherorder polynomial expansion(HOPL for short) for the mass flux. In addition to Van

Leer's splitting, a term is introduced so that the mass diffusion error vanishes at M = 0.

Several splittings for pressure are proposed and examined. The effectiveness of the

HOPE scheme is illustrated for 1-D hypersonic conical viscous flow and 2-D supersonic

shock-wave/boundary-layer interactions. Also, we give the weakness and suggest areas

for further investigation of the scheme.

Introduction

In the past decade, upwind differencing schemes have gained considerable attention

for their accuracy and robustness in Euler flows with discontinuieties, shock waves in

particular. Naturally, significant research effort in the CFD community has been focused

on maxmizing the accuracy and efficiency, among other objectives. Four popular but

conceptually different flux splitting ideas have been utilized for nearly 10 years: Steger

and Warming, Van Leer, Roe, and Osher. However, each scheme has an associated

weakness when numerical accuracy and efficiency are considered.

In this paper, we deal specifically with the improvement of Van Leer's flux vector splitting[1]. Besides its simplicity, Van Leer's splitting has the following properties: (1) it can

be interpreted as a special member of a family of second-order polynomial expansions[2],

and (2) the associated flux Jacobian and eigenvalues are continuous at the sonic points.

Van Leer's choice allows one vanishing eigenvalue in the case of an ideal gas, thereby

resulting in a crisp shock representation. Furthermore, the continuous differentiability

is helpful for convergence acceleration, e.g., in multigrid schemes.

However, failing to recognize the contact discontinuity, the Van Leer splitting[1] produces excessive numerical diffusion and thus requires a huge number of cells to correctly

resolve the boundary-layer flow. Some improvements have been demonstrated recently

by Hanel et al[2] and Van Leer[3] for 1-D conical, hypersonic viscous flow, but a pressure

glitch arises. A new scheme by the present authors[4] has been proposed that not only

corrects this pressure difficulty, but also is remarkably simple to implement. Nevertheless, the above schemes[2-4] have already departed from the ideas of flux vector splitting

t Senior Scientist, Internal Fluid Mechanics Division

* Aerospace Engineer, Computational Fluid Dynamics Branch

146
and in fact become more like the flux difference splitting. Since the differentiability and

simplicity are desirable properties, one would still wish to search for a better splitting

scheme that is strictly based on the flux vector splitting.

In this paper, we propose a family of higher-order polynomial expansions for the mass

flux that diminishes the diffusion error as M —► 0. We give a detailed study of the

accuracy of the scheme for 1-D conical flow and 2-D shock wave/boundary-layer interactions. The weakness of the scheme is also pointed out and possibile improvements

suggested.

Analysis

To exemplify the concept, let us consider the quasi two-dimensional system of equations

for conical flows:

8U OF

+

= S at 8 rl

where UT = ( P, P u , P v , PE ), FT = (PV , P vu , P V ' + p , pvH), E = e + 1 /2(u' + 0), and

H = E + plp. The flow considered consists of a very thin shear layer at the wall

and a shock wave away from the wall. An algorithm must be capable of minimizing

the numerical smearing(diffusion) at the locations where an eigenvalue changes sign or

approaches zero. For example, Van Leer's splitting[l] can represent shock profile well,

while greatly diffusing the boundary layer. The Van Leer split mass fluxes are:

Fl = Fl + Fi ; F} = fpa/4(M f 1)'.

The net difference from the curve it approximates is largest at M = 0; its value equals

pa/2. This error, viz numerical diffusion, significantly broadens the boundary layer,

leading to incorrect velocity and temperature profiles. A simple way to remove the

diffusion at M = 0 is by adding an extra higher-order term that allows the split mass

fluxes to pass through the origin(Fig. 1), i.e.,

F± = fpa/4[(M f 1 )' + m i(M )(M ' — 1)'],

where the higher-order term has a coefficient m l , in general function of M. It should

have the following properties:

(1) m l —4 —1 as M --► 0;

(2) m i (M) = ml(—M);

(3) m l --+ 0 as M —+ fl.

A formula satisfying those properties is chosen as:

m l = ( M' — 1 )/(M' + 1)S,

where the exponent S is a free parameter; also shown in Fig. 1 is m l vs M with S = 4.

147
1.5

1.0

H 0.5

0.0

Co

—0.5

1.0

1.5

(F! +F,-)/pa

-- F:/pa

- - - - F,-/pa

^^ /^E^ mr(^)

1.0

0.8

^y 0.6

Q 0.4

0.2

0.0

(Z)

(3)

(4)

^j

/

—1.0 0.0 1.0 —0.6 0.0 0.6

Mach No. Mach No.

Fig. 1. Mass Flux Splitting in HOPE. Fig. 2. Pressure Splitting.

In the conical flow calculations, the accuracy and convergence appear to be insensitive

to the specified values of S = 2, 4, 6. Now, regarding the flux as a sum of convective and

pressure terms, we can write the splitting formula for the flux vector:

Fl } PV 1 0

Fz pvu _ 

Ff

u 0

F3pvv + p ' v + p}

F4 pvH H 0

With the realization in [5] that the pressure splitting could be considered separately in

the Van Leer formula[l], a whole host of freedoms for the pressure splitting becomes

possible. Following is the list of formulas tested:

(pl) P} = T-1/4(M f 1)'(M :F 2)P,

(p2) : p} = (pl) :F 3/4M (M' — 1)2p,

(p3) : p} = (pl) f 3/4m 1 M(M' — 1)2p,

and

(P4) : P} = 1/2(1 f 7M)p.

Figure 2 displays the distribution of the split pressure vs M. The first formula is that

used by Van Leer[1]. The second and third splits, (p2) and (p3), yield vanishing

slope at M = 0, thus corresponding to central differencing. However, no instability

was encountered in the conical flow problem with the (p2) or (p3) split used in an

implicit code. The fourth split (p4) is obtained from an approximate integration along

characterics. As will be seen later, the four formulas give essentially the same results

for the conical flow calculated.

sults And Discussion

In this paper, two cases were tested to check the accuracy and convergence of the

HOPE scheme. The first case is the 1-D self-similar conical flow over a 10-degree

half cone at hypersonic speed, for which a detailed comparison study was conducted.

148
14.0

13.5

0(e) HOPE

o O(A-) HOPE

M = 7.95

Pr = 1.0

0(e) HOPE

0 0(62) HOPE

M = 7.95

Pr = 1.0

4 8 12 1 16

TIT_ (T,.,a,3T_). ,

2 3 4 5

PAP

14.0

13.5

13.0

12.5

0 12.0

11.5 P

11.0

10.5

10.0

0

13.0

12.5 o^

o, 12.0

11.5

11.0

10.5

10.0

0

The flow conditions are: M,,. = 7.95, and Rem = 4.2 x 10 5 . Since P* = 1.0, exact

solution gives adiabatic wall temperature, 13.64T., . The second case is the 2-D shock

wave/laminar boundary-layer interactions, for which experimental measurements were

available[6]. The conditions are: M., = 2.0, Rem = 2.96 x 10 6 , and oblique shock

angle 0 = 32.585 degrees. In both cases, the results from the Roe splitting are also

included for comparison. An implicit Newton iteration procedure was used to achieve

steady-state solution with L,,,, residual dropped by five orders of magnitude.

Figures 3 and 4 show the pressure and temperature distributions from the first- and

second-order solution on a 65-grid; little difference is seen. A monotone solution across

the shock is obtained with the first-order scheme while oscillation appears in the secondorder scheme, which can be eliminated by a TVD procedure. It is noted that the firstorder pressure is smooth at the edge of the boundary layer, unlike the Roe solution which

shows a slight discontinuity(not shown here). Although the boundary layer exhibits a

steep temperature gradient, the HOPE scheme predicts the wall temperature correctly,

indicating removal of the numerical diffusion associated with the original Van Leer

splitting.

Fig. 3. Pressure Profile of Conic Flow. Fig. 4. Temperature Profile of Conic Flow.

Figure 5 displays the results using various pressure splittings; they are practically identical except the Van Leer pressure split (pl) shows some minor oscillation near the wall.

However, the pressure splittings show significant effect on the convergence rate. The

(p3) and (p4) splits are the best, comparable to the Roe splitting, while the other two

are roughly two to three times slower. These may indicate possible instability in a more

complex case.

Finally, for the 2-D case, the surface pressure and friction coefficient are plotted in

Figs. 7 and 8. The first-order HOPE results compare fairly with Roe's splitting and

experimental data. However, the second-order calculation experienced difficulty in convergence in which the residual was reduced by only two orders of magnitude and the

result is not presented here.

149
-- ROE Splitting

HOPE with P•,

HOPE with Pit

---- HOPE with P-3

----- HOPE with P•4

101

100

10-1

10-2

8

10-3

10-4

i0-5

i

oT

f

e

f

f

• ROE Splitting =

f

o HOPE with P-,

o HOPE with P=t

HOPE with Psj f

• HOPE with Pt4

M = 7.95

Pr = 1.0 i

t

_

_s

f

1 4.6

1 3.5

1 3.0

y 1 2.5

N

o, 12.0

11.5

11.0

10.5

10.0

0 1 2 3 4 5 0 2000 4000

p/pm

iterations

Fig. 5. Comparison of Pressure Splits. Fig. 6. Convergence History.

We suspect that a further investigation on other pressure splittings may lead to success

in stability and convergence. Nevertheless, a systematic study of the eigenvalues of

the split fluxes and the complete discretized system will prove to be a useful endeavor.

Above all, the present research suggests that there are still possibilities in flux-vector

splitting after Van Leer's appeared nearly 10 years ago. The possibilities may very well

still lie in the mass-flux splitting.

2.0 10

8

1.6

6

Cn

1.2

o_ 4

X 2

0.8

C^~ 0

0.4 -2

- 4

0 0

0.0 1.0 x 2.0 0.0 1.0 x 2.0

Fig. 7. Pressure at the Wall. Fig. 8. Friction Coefficient at the Wall.

References

[1] B. van Leer, Lecture Note8 in Phy8ic8, Vol. 170, 1982.

[2] D. Hanel and R. Schwane, AIAA Paper 89-0274, 1989.

(3] B. van Leer, CFD Symposium on Aeropropulsion, NASA CP-10045, 1990.

[4] M.-S. Liou and C. J. Steffen, Jr., (in preparation).

[5] M.-S. Liou, B. van Leer, and J. S. Shuen, J. Comput. Physics, Vol. 87, 1990.

(6] R. J. Hakkinen, I. Greber, L. Trilling, and S. S. Abarbanel, NASA Memo 2-18-

59 ], 1959.

150
Three-Dimensional Calculations of Supersonic Reacting Flows

Using an LU Scheme

Sheng-Tao Yu, Y.-L. Peter Tsai, and Jian-Shun Shuen

Sverdrup Technology, Inc.,

NASA Lewis Research Center, Cleveland, Ohio

ABSTRACT

A new three-dimensional numerical program that incorporates comprehensive real

gas property models has been developed to simulate supersonic reacting flows. The code

employs an implicit, finite volume, Lower-Upper (LU), time-marching method to solve

the complete Navier-Stokes and species equations in a fully-coupled and very efficient

manner. A chemistry model with nine species and eighteen reaction steps is adopted

in the program to represent the chemical reactions of H2 and air. To demonstrate the

capability of the program, flow fields of underexpanded hydrogen jets transversely injected

into the supersonic airstream inside the combustors of scramjets are calculated. Results

clearly depict the flow characteristics, including the shock structure, the separated flow

regions around the injector, and the distribution of the combustion products.

151
Simnlations of Free Shear Layers Using

a Compressible k — c Model

S.T. Yu'

Sverdrup Technology, Inc.,

NASA Lewis Research Group

Brook Park, OH 44142

C.T. Changt and C.J. Marek$

NASA Lewis Research Center

Cleveland, Ohio 44135

ABSTRACT

A two-dimensional, compressible Navier-Stokes

equations with a k — c turbulence model are solved

numerically to simulate the flows of compressible free

shear layers. The appropriate form of k and c equations

for compressible flows are discussed. Sarkar's modelling

is adopted to simulate the compressibility effects in the

k and c equations. The numerical results show that the

the spreading rate of the shear layers decreases with

increasing convective Mach number. In addition, favorable comparison was found between the calculated

results and Goebel and Dutton's experimental data.

INTRODUCTION

Recent national interest in trans-atmospheric vehicle has rekindled the hypersonic research. For this

vehicle, a supersonic combustion ramjet (scramjet) engine was proposed to provide the power. Inside this

scramjet engine, compressible mixing layers are important phenomena. The performance of the engine will

depend on the supersonic mixing and the flame holding

of shear layers.

The behavior of incompressible mixing layers has

been studied extensively. However, additional study

is required to understand the compressibility effects of

free shear layers at high speeds. Figure 1 shows a sketch

of a typical free shear layer. Two streams at different

temperatures, densities, and Mach numbers merge together to form a free shear layer. Various combinations

of flow conditions of high-speed streams and low-speed

streams allow for the systematic study

Research Engineer, Member AIAA.

t Aerospace Engineer, Member AIAA.

$ Senior Research Scientist.

of compressible shear layers. In this paper, we report

the incorporation of a k — c model with compressibility

effects to a two-dimensional flow equations solver for

the simulation of compressible shear layers. First. we

point out the derivation procedure of the compressible

k and c equations. Unlike the procedure in the incompress equations, both Favre and Reynolds averaging

procedures) are performed in deriving the flow and turbulence equations. Particularly, the additional terns in

the Navier-Stokes equations due to the averaging procedure are illustrated. These terms are often omitted

in CFD practices. The k and c equations are presented

in vector form for the convenience of illustrating the

numerical method.

The lower-upper (L U) scheme developed by loon

and Jameson' is adopted in this work. This method has

proven very efficient for large systems of equations. 3 For

completeness, a brief account of the numerical method

is presented in this paper. The newly developed solver

then is applied to simulate compressible free shear layers with five different convective Mach numbers from

0.05 to 1.48. One of five test conditions is a replica

of that reported by Goebel and Dutton.' The most

important feature of the compressible free shear layers one want to demonstrate in the calculations is the

decrease of the shear layer thickness with increasing

Mach number.' Favorable comparison were found between the experimental data and the calculated results.

THEORETICAL MODEL

In deriving the compressible flow equations of fully

turbulent flows, all the flow properties are Favre averaged (mass weighted averaged) except the densit y. p,

and pressure, p. The definition of the Favre average is

^ _ PIP• (1)

152
where 0 is any flow property. 'Thus flow variables are

decomposed in the following fashion,

0 = ^ + O" - (2)

On the other hand, conventional Reynolds' average are

used for the pressure and the density. According to

the definition of Favre averaging, the following relations

exist:

P'O'

P (3)

Pon = 0,

T = 0.

In doing so, all the terms associated with the density

fluctuation, e.g., p'u', in the Reynolds' averaged equations were eliminated in the Favre's averaged equations.

The resulting flow equations are much simpler compared to the equations derived by the Reynolds averaging proceduresWritten in a strong conservative form, the turbulent, compressible, flow equations. can be expressed as

follows:

Q

a

+ a (E—Ev )+ a (F—Fv )=H. (4)

Here x and y are Cartesian coordinates, Q is the dependent variable, E and F are the convective flux vectors

and E„ and F„ are the viscous flux vectors. The equations are similar to the laminar equations. However,

all variables in the equations are the averaged variables, and the transport properties are the effective,

i.e., laminar plus turbulent, properties. Here, we want

to point out that the viscosity multiplied by the dilatational terms in the normal stresses is still the laminar

viscosity as illustrated in the following relations:

au 2 au av Txx = 2 (^ + /at)ax — 3μ ax + ay)'

N2 au av 1 (5) TYY — 2(ta + pt ) ay 3 μ (ax + a Jy

The vector H on the right hand side of the flow

equations, Eq. 4, represents the additional terms introduced by the averaging procedure. The vector H can

be expressed as

0

z 8pk

— a ax

2 ank

H = —3 ay (6)

Tr_ [(^ + 0", ax a [(^ + at ay J

+2ua (Pk)+2v-L(pk)+G ay

where k is the turbulent kinetic energy and is defined

as

k

 = 2

p(u„ u„ +v"v„)

(7)

P

G is the generation of the turbulent kinetic energy and

is defined in the following section.

The equations of turbulent kinetic energy, k, and

energy dissipation rate, c, are derived by the manipulation of the flow equations and the averaging procedure.

The derived equations of k and c can not be solved directly due to the closure problem. Modelling of certain

terms in the equations is necessary to make the governing equations well posed. The details of the modelling

is beyond the scope of this work. Here, only the final

form of the equations are presented. The k and c equations in the Cartesian coordinate system can be cast

into the vector form: .

aQke + aEk^ + aFk , _ aEv ke + aFvke + S (8)

at ax ay ax ay

where

Tk 1

Qke = PC J ,

Gii(

uk Eke = '

pvk

£' ke = PVE '

^ ak (({L ak) as 1 ,

Evke

+

I\ ( 

ax

/ (μ ,k ) + aY Fvke 1`/

l{L+ o )aY

S—C G—Pc(1+aM=)1

, (9) (C1 G — C2pf) k

where G is the generation term of the turbulent kinetic

energy and can be expressed as:

1 au; aui 2 auk2— auk 

(10) G = μt [2 (axj + axi ) — 3 ax k — 3pkaxk 

where the subscripts i, j, and k follows the convention

of Cartesian tensor.

The eddy viscosity ta t is derived in terms of appropriate length and velocity scales. For the k — c turbulence model, the length scale and the velocity scale of

turbulent fluctuations are taken as

_ k3/2

P

c (11)

u»

r kt/2

153
These relations allow the eddy viscosity, Eft, to be modelled as:

_kZ

μt = CN P

e 

(12)

The constants used in the k - E model are the standard

Jones and Launder's values: 6 C,, = 0.09, Cl = 1.44,

C2 = 1.92, vk = 1.0, and Q, = 1.3. These constants

were never altered during the course of this work.

In order to accommodate the compressibility effect, the dissipation term pE in the source term of the k

equation is multiplied by a correction factor (1+aMt ).

Here M= is the local turbulence Mach number defined

as Mt = v/-k-/a where a is the local speed of sound. The

constant a in the term is taken as unity. This model

is developed by Sarkar et al.' The physical meaning of

the term is that for high turbulence Mach number (Mt)

flows, the dissipation of the turbulence kinetic energy

is enhanced by a factor of aMt . For free shear layers at

high convective Mach number, the turbulence intensity

is greatly reduced due to compressibility.

In calculating the turbulent free shear layers, the

inlet boundary conditions for mean velocities and temperature are specified based on the hyperbolic tangent

profile with specified initial shear laver thickness. The

hyperbolic tangent profile is an approximation of the

self-similar solution for fully developed turbulent free

shear layers. The inlet transverse velocities are set to

be zero. The turbulent kinetic energy and dissipation

are specified according to the local equilibrium assumption and a algebraic turbulence model :8

{^c = Pln, ay (13)

where 1,,, = 0.1256 and the shear layer thickness, 6, is

based on the distance between the two transverse locations where u = u l - OAAii and u = u2 + OAAu.

The dissipation can be related to the local length scale

which is specified based on the local shear layer thickness:

E 

_ 

C, 

k5

6

(14)

where C, = 1.23. Using Eqs. (13) and (14), and the

equation for the eddy viscosity, i.e, Eq. (12), k and

E can be readily obtained for the upstream boundary

conditions.

Numerical Method

The numerical solution of Eqs. (4) and (8) is

performed in a general, body-fitted coordinate system,

(f, p). For the purpose of discussion, we will concentrate on Eq. (4). However, the procedure is equally

applicable to Eq. (8). Coordinate transformation of

Eq. (4) results in:

aQ + ^^ (E - Ev ) + j^ (F - FH(15) v ) =

where

Q=hQ

E = h({SE+^,F)

F = h ( rl. E + r)v F) (16)

Ev = h (^.Ev +^,F,,)

Fv = h ( t7zEv + rlyFv)

H=hH

in which It is the cell volume.

The transformed equaticn, Eq. (15), is solved using a time-marching, LU scheme. The LU scheme can

be obtained by approximately factorizing the left-handside (LHS) of the equation. In time-marching form, the

implicit upwind difference scheme of Eq. (15) can be

written as

[I+ At(D{ A - + D+ B- - D+

D{ A + + D, h + )]AQ = AtRHS

(17)

In Eq. (17), At is the time-step. Backward-difference

operators are denoted by D- and D, , and forwarddifference operators by D{ and D,+. The flux Jacobians, A+ , B+ , A - , and B- are constructed such that

the eigenvalues of `+' matrices are nonnegative and

those of `-' matrices are nonpositive. The matrix D

is the Jacobian matrix of the source term.

The LHS matrix of Eq. (21) is usually too large for

direct inversion. An approximate-factorization procedure is implemented which results in the following LU

scheme :

I + At (D-,k + + D,^ B + - Q^ - 3 - D I J

+ At(B+-B-)-D I i

Atj

// A+ B+

/J

1111

I + At I D{ A- + D+ h - + Q^ + Oq - D I J AQ

= AtR \

ff

(18)

where the grid spacing in the general coordinate, A^

and An are usually taken to be one. R represents the

residual of each LU time marching step. In calculating

the residual, R, both the inviscid and viscous terms are

discretized using the central-difference approach:

R = D^(Ev - E) + D,(Fv - F) + II (19)

154
where Df and D, are the central difference operators.

Equation (18) is the generic form for the LU

scheme. Its derivation can be found in Ref. 3 and

will not be repeated here. This LU scheme requires

inversion of the matrix,

[I + At (A+ — A - + B+ — I3 - — D) ] (20)

for the L operator and

[I+At(A,+— A- +II+— B- -bjD)]

(21)

for the U operator.

Up to this point, no definition has been made to

the exact form of the split flux Jacobians. Yoon and

Jameson2 proposed that the split flux Jacobians are

defined as

A+ = 0.5(A + yAI)

A - = 0.5(A — yAI)

(22)

B+ = 0.5(B + yBI)

B - = 0.5(B — yBI)

where 7p and yII are greater than the spectral radii of

the associated flux Jacobians :

yA > max(IAAI)

(23)

yn > max(JA,3J)

The purpose of constructing split flux-Jacobians by Eq.

(22) is to make the matrices in Eqs. (20) and (21) diagonal for efficient inversion. Apparently, the eigenvalues

of the split flux-Jacobians are not the characteristics

speeds of the flow.

In solving the k and c equations, the aforementioned numerical method, i.e., the LU scheme on the

left hand side and central differencing on the right hand

side, is used. The solution procedure of the whole equation set is decoupled into flow solver and turbulence

solver. Thus, the turbulence solver stands alone and

can be easily turned on or off. This arrangement does

not affect the overall numerical stability due to the fact

that the feedback from k and c equations to the flow

equations depends on the turbulent transport properties only. Thus, it is more efficient and convenient to

separate the solution procedure into two parts.

The source terms of the k and c equations demand

special treatment. In linearizing the source terms for

the numerical method, the Jacobian matrix is obtained

through the derivative of the source terms with respect

to the dependent variables, i.e., pk and pc. Following

the usual practice, the form of the source terms guarantee a 2 x 2 full matrix for the Jacobian matrix. However, special treatment in deriving the Jacobian matrix

is applied in this work to enhance the numerical stability. In the k equation, c has been replaced by k(c/k)

where c/k is treated as a constant. A similar method is

also applied to the source term of the c equation. The

Jacobian matrix obtained is:

k(1 +aMi) CO', D — )

0 —CZ k

27

Note that off-diagonal terms are eliminated and the diagonal terms are always negative. Thus, the implicit

part of the source terms of k and c equations behaves

like a sink which always stabilize the numerical scheme.

Results and Discussions

Bogdanoff' introduced the convective Mach number as a parameter that collapses the growth rate date

of plane shear layers. The convective Mach number is

defined as

MC = 

Ul — U` 

= 

U` — 

U2 (28) 

cl c2

where Ul and U2 are the freestream velocities, and cl

and c2 are the freestream sound speeds. Uc is the convective velocity of large structures and is defined as

UC _ 

U1C2 + U2C1 (29)

C1 + C2

Table 1 shows the five test conditions of the simulated

compressible free shear layers. The range of the convective Mach numbers, M, is from 0.05 to 1.48. The last

two rows show the spreading rate and the ratio of the

compressible spreading rate to incompressible spreading rate in which b is the shear layer thickness and the

superscript ' represents the derivative of b with respect

to the streamwise distance. As indicated in the last

row of the table, the ratios of the spreading rate of

the free shear layers decreases as the convective Mach

number increases. The test conditions of Case 3 in the

table are the same as in the experiments reported by

Goebel and Dutton.' In the rest of the section, we will

first show the direct comparison between the simulated

results and the experimental data for Case 3. Then,

detailed numerical solutions of Case 3 in terms of velocity, turbulent kinetic energy, turbulence dissipation

rate, Reynolds stress, and eddy viscosity are presented

in a coalesced fashion. Finally, the solutions of five

cases are compared to each other in the figures of kinetic energy, Reynolds stress, and ratio of spreading

rates.

The solutions of Case 3 are examined in detail.

Figure 2 shows the development of the free shear layer.

155
Note that x and y axes are not on a 1:1 ratio for the

convenience of illustration. The definition of the shear

layer boundaries is the same that of the shear layer

thickness. III 2, the boundaries of the shear layer

corresponding to 10 — 90% are drawn. Circles are the

experimental data of Dutton et al. The calculated results agree well with the experimental data. After the

developing region, the boundary of the shear layer is

almost linear. Incidentally, Figs. 3a and 3b show the

numerical convergence trends of the flow and k — c

equations. In about 2500 iterations, the residuals drop

about 12 orders of magnitude and reach the machine

accuracy.

Figure 4 shows the Mach number profile at various

axial locations. The velocity gradient in the transverse

direction decreases as the flow goes downstream. Figure

5 is the coalesced version of Fig. 4. The nondimensionalized y coordinate defined as (y— y,)/b is used, where

y, is the transverse location at the center of the shear

layer, and b is the shear layer thickness. Note that the

upstream boundary condition of velocities is prescribed

according to a hyperbolic tangent curve which is an approximation of self-similar solution of free shear layers.

According to Fig. 5, this self similarity of velocity profiles never fail as the flow goes downstream.

Figure 6 shows coalesced turbulence kinetic energies at various locations. The turbulence kinetic energy is nondimensionalized by AU'. This figure clearly

shows that after the first three stations, i.e., about 150

mm, the turbulence kinetic energy retains self similarity. Thus, the developing region for the turbulence is

about 150 mm. Figure 7 shows the turbulence dissipation profiles. As flow goes downstream, the peak values

of turbulence dissipation at each stations decrease, and

the turbulence dissipation never reach a fully developed

condition. However, if the c value is nondimensionalized by AU3 16 the coalesced profiles appear as shown

in Fig. 7. A similar behavior is observed for the eddy

viscosity profiles (Fig. 8). Figure 9 shows the nondimensionalized Reynolds stress profiles. Again, the turbulence Reynolds stress becomes the fully developed at

about 150 mm downstream of the splitter plate.

Figure 10 shows the comparison between the predicted fully developed Reynolds stress and the experimental data reported by Dutton et al.' The predicted

solution underestimated the peak value of the Reynolds

stress profile by 6 — 8%; however, the overall trend

of the predicted results is correct. Many factors contribute to the discrepancy between the predicted result

and experimental data. Among them, the upstream

boundary conditions are simplified in the solution procedure, i.e., no effort was made to simulate two boundary layers merging at the tip of the splitter plate. This

could offset the solution in the developing region and

shift the fully developed solutions.

Figure 11 shows the comparison of the turbulent

kinetic energy between the five cases. Again, the turbulent kinetic energy is normalized by the square of the velocity difference of the two streams. It is clear that the

normalized turbulence intensity decreases with increasing convective Mach number. A similar situation is

observed in Fig. 12, the normalized Reynolds stress decreases with increasing convective Mach number. Figure 13 shows the distribution of the ratios of spreading rates for the five cases compared to experimental

data. The ratio of the spreading rates decreases from

about unity to 0.45 as the convective mach number increases from about 0. to 1.45. Both the experimental

data and the simulated results show the spreading rate

ratio reaches an asymptotic value after the convective

mach number exceeds unity. Simular phenomenon can

be seen in Figs. 11 and 12. Both figures show that

the normalized turbulence kinetic energy and Reynolds

stress reach asymptotic values as the convective Mach

number increases.

CONCLUDING REMARKS

In this paper, we report the incorporation of a compressible k, — c model with a two-dimensional Navier

Stokes solver to study compressible free shear layers.

Sarkar's modelling is adopted to simulate the compressibility effects of the k and e equations. The model enhances the turbulence dissipation rate of flows at high

speeds. In deriving the governing equations, the Favre

averaging procedure for the fully trubulent flow equations is elaborated. The equation sets are presented in

the vector form for the convenience of the discussion of

the numerical method. Yoon and Jameson's LU scheme

is used to solve the equation sets. Details of boundary

conditions and the treatment of the source terms of the

k and c equations are discussed. Then the program is

applied to simulate compressible free shear layers with

five different convective Mach numbers. The decrease

of the spreading rate of the shear layers with increasing

Mach number is observed in the calculated results. Results also show favorable comparison with Goebel and

Dutton's experimental data.

REFERENCES

1. Kuo, K.K., 1986 "Principles of Combustion," John

Wiely & Sons, p. 412.

2. Yoon, S. and Jameson, A., 1989 "Lower-Upper Implicit Schemes with Multiple Grids for the Euler

Equations," AIAA J. Vol. 25, No. 7, p. 929.

156
3. Yu, S.T., Tsai, Y.P., and Shuen, J.S., 1989 "ThreeDimensional Calculation of Supersonic Reacting

Flows Using an LU Scheme," AIAA paper 89-0391

(accepted for publications in J. of Compt. Phys.)

4. Goebel, S.G., and Dutton, J.C. 1991 "Experimental Study of Compressible Turbulent Mixing Layers," AIAA J. Vol. 29, No. 4, p.538.

5. Bogdanoff, D.W., 1983 "Compressibility Effects in

Turbulent Shear Layers," AIAA J. Vol. 21, No. 6,

p.926.

6. Jones, W.P., and Launder, B.E., 1972 "The Prediction of Laminarization with a Two-Equation

Model of Turbulence," Int. J. Heat Mass Transfer,

Vol. 15, p. 301.

u,

U,

_ _C ~'c

uz

U,

Fig. 1: Schematic of free shear layer.

7. Sarkar, S., Erleba.cher, G., Hussaini, M.Y., and

Kreiss, H.O., 1989 "The Analysis and Modeling of

Dilatational Terms in Compressible Turbulence,"

ICASE Report No. 89-79, NASA Langley Research Center.

8. Burr, R.F. and Dutton, J_C., 1990 "Numerical Modeling of Compressible Reacting Turbulent

Shear Layers," AIAA Paper 90-1463.

9. Papamoschou, D. and Roshko, A. 1988 "The Compressible Turbulent Mixing Layer: An Experimental Study," J. Fluid Mech., Vol. 197, p.453.

Table 1: Test conditions of the calculated shear layer.

CASE 1 2 3 4 5

M I ,M2 1.2, 1.1 2.0, 1.4 1.9 1, 1.37 2.5, 1.1 6.1, 3.15

Mc 0.0.5 0.31 0.45 0.70 1.49

U I ,U 2

[m]

419, 394 676,465 702,404 965, 390 3461, 1786

T 1,T2

[K]

300, 300 275,275 334,215 295,295 800,900

PI , p^

(Kg/-m 3)

0.64,0.6 0.7,0.7 0.57,0.89 0.64,0.64 0.24,0.24

P [atml 0.55 0.55 0.55 0.55 0.55

S' 0.007 0.021 0.027 0.035 0.024

1.01 0.7 0.54 0.47 0.40

157
15®

0

0 Dutton's Data

R

0

O

m

u 0

4

d

-15 1 Fig. 4: The velocity profiles of the free shear layer at dif- 0 100 200 300 400 500 mm

fercnt axial locations.

Axial Location

Fig. 2: The boundaries of the free shear layer (10--90%).

1.s

1

0.5

b

V

a

8

I

v

-1

-1.5

488

A

( a )

-s

d

0

v -SB

J

- ts ^

A

iteration nunoer

-6

(b) x 48.333 + 298.33

. 98.333 O 348.33

O 148.33 • 398.33

0 198.33 • 448.33

0 248.33 498.33

i.

0

J -III Fig. 5 Coalesced velocity profiles of different axial locations.

-21

8 1888 2888 3868 lBAB 5888

Iteration number

Fig. 3: Convergence trends of the flow and turbulence

equations.

158
-1

L 

1

^Xr ^

0.3

b

o x

a e

X 

X

x x x

8 0.81 0.62 6.03 0.91

K'/AU--2, ((^/.7^27

Fig. G: Coalesced turbulent kinetic energy profiles of different axial locations.

1.3

1

e.s+

. ,

b ` x

• X x

X

'

e

X

+ x

e 8.001 e.e00 0.812 0.816

<.*-'>/.W2 N. "..d R.L..Id 8—...

Fig. 9: Coalesced Reynolds stress profiles of different axial

locations.

1.3

1

e.s

b

a e

• +• •

x

x X X

X

—1 —.XX

—1 3

6 6.882 0.684 0.006 0.966

Fig. 7: Coalesced turbulence dissipation profiles of different axial locations.

Dutton's Data

^.s

(

7

( ^.w r. r( ^. ul ^.w

W'•'.tr.rt ..+.iwa ^•.ru. t.....

Fig. 10: Comparison of Reynolds stresses between experimental data and predicted results.

e e.eez e.ea4 8.886 0. ass 8.91

μ j(p.6nu)

Fig. 8: Coalesced eddy viscosity profiles of different axial

locations. 159
1.2

1

e,s

0,6

O

b

8.4

0.2

0

I

I.s

1

^. s

1.5

-i —

1 IM5 1.01 8.915 1.12 1.8S5

e l'1.W1 $—]I." TY WIW lC

Fig. 11: Comparison of turbulent kinetic energy with five

different convective Mach numbers.

I.5

1.5

L

-1.5

_S

-i."4

I<^'a•a l/.Y7 Yrw M" S—I

Fig. 12: Comparison of Reynolds stresses with five different

convective Mach numbers.

q — —_

El

o

q —

o q

A

. o

0,5

N^

1,5 2

A Papamoschou and R.oshko [9]

[{ Langley's data [7]

-^- With Compressibility Model

- - - Without Compressibility Model

Fig. 13: Comparison of ratios of spreading rate between experimental data and predicted results.

160
Form Approved

REPORT DOCUMENTATION PAGE OMB No. 0704-0188

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

gathering and maintaining the data needed, and completing and reviewing the collection of information. Send comments regarding this burden estimate or any other aspect of this

collection of information, including suggestions for reducing this burden, to Washington Headquarters Services, Directorate for information Operations and Reports, 1215 Jefferson

Davis Highway, Suite 1204, Arlington, VA 22202-4302, and to the Office of Management and Budget, Paperwork Reduction Project (0704-0188), Washington, DC 20503,

1. AGENCY USE ONLY (Leave blank) 2. REPORT DATE 3. REPORT TYPE AND DATES COVERED

October 1991 Technical Memorandum

4. TITLE AND SUBTITLE 5. FUNDING NUMBERS

Center for Modeling of Turbulence and Transition (CMOTT)

Research Briefs-1990

WU-505-62-21

6. AUTHOR(S)

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

REPORT NUMBER

National Aeronautics and Space Administration

Lewis Research Center E-6570

Cleveland, Ohio 44135-3191

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

AGENCY REPORT NUMBER

National Aeronautics and Space Administration NASA TM-105243

Washington, D.C. 20546-0001 ICOMP-91-17;

CMOTT-91-07

11. SUPPLEMENTARY NOTES

Report compiled by Louis A. Povinelli, Meng-Sing Liou, and Tsan-Hsing Shih, NASA Lewis Research Center (work

funded under Space Act Agreement C-99066-6). ICOMP Program Director: Louis A. Povinelli, (216) 433-5818.

12a. DISTRIBUTION/AVAILABILITY STATEMENT 12b. DISTRIBUTION CODE

Unclassified - Unlimited

Subject Category 01

13. ABSTRACT (Maximum 200 words)

This research brief contains the progress reports of the CMOTT Research staff from May 1990 to May 1991. It is

intended to be an informational report of the CMOTT activities as well as an annual report to ICOMP and NASA. The

Center for Modeling of Turbulence and Transition (CMOTT), a cooperative turbulence research team, was formally

established in May of 1990 as a focal group within the Institute of Computational Mechanics for Propulsion (ICOMP).

The location of CMOTT is shown in the organization structure chart in Appendix A. The objectives of the CMOTT are

to develop, validate and implement the models for turbulence and boundary-layer transition in the practical engineering

flows. The flows of interest are three-dimensional, incompressible and compressible flows with chemistry. The

schemes being studied include the two-equation (e.g., k-e) and algebraic Rcynolds-stress models, the full Reynoldsstress (or second moment closure) models, the probability density function (pdf) models, the Renormalization Group

Theory (RNG) and Direct Interaction Approximation (DIA), the Large Eddy Simulation (LES) and Direct Numerical

Simulation (DNS).

14. SUBJECT TERMS 15. NUMBER OF PAGES

Turbulence models 170

16. PRICE CODE

A08

17. SECURITY CLASSIFICATION 18. SECURITY CLASSIFICATION 19. SECURITY CLASSIFICATION 20. LIMITATION OF ABSTRACT

OF REPORT OF THIS PAGE OF ABSTRACT

Unclassified Unclassified Unclassified

NSN 1b40-U1-28U-5500 Standard Form 298 (Rev. 2-89)

Prescribed by ANSI Std. Z39-18

298-102