Computational Modeling of Turbulent Flows*

.

ADVANCES IN APPLIED MECHANICS VOLUME

18

Computational Modeling of Turbulent Flows?. JOHN L. LUMLEY Sibley School of Mechanical and Aerospace Engineering Cornell University Ithaca. New York

I . Introduction . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . A. History and Generalities . . . . . . . . . . . . . . . . . . . . . . . . . . . B. General Assumptions . . . . . . . . . . . . . . . . . . . . . . . . . . . . . I1 . Mathematical Preliminaries . . . . . . . . . . . . . . . . . . . . . . . . . . . A. Representations . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . B. Realizability . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 111. The Return to Isotropy . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . A . Introduction . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . B. The Reynolds Stress . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . C. The Heat Flux . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . IV . The Rapid Terms . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . A. Introduction . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . B. The Heat Flux Integral . . . . . . . . . . . . . . . . . . . . . . . . . . . . C. The Temperature Variance Integral . . . . . . . . . . . . . . . . . . . . . D. The Reynolds Stress Integral . . . . . . . . . . . . . . . . . . . . . . . . . V . The Dissipation Equations . . . . . . . . . . . . . . . . . . . . . . . . . . . . A. The Mechanical Dissipation . . . . . . . . . . . . . . . . . . . . . . . . . B. The Thermal Dissipation . . . . . . . . . . . . . . . . . . . . . . . . . . . C. The Transport Terms . . . . . . . . . . . . . . . . . . . . . . . . . . . . . VI . The Transport Terms . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . A. Introduction . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . B. A Gaussian Model . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . C. Order of Magnitude Analysis . . . . . . . . . . . . . . . . . . . . . . . . D . The Pressure Transport . . . . . . . . . . . . . . . . . . . . . . . . . . . . E. Zeroth-Order Transport Terms . . . . . . . . . . . . . . . . . . . . . . . F. Modifications . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . References . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .

124 124 127 128 128 131 133 133 137 142 143 143 147 147 150 152 152 156

159 160 160 162 165 168 170 171 174

t This work supported in part by the U S. National Science Foundation. Meteorology Program under Grant Number ATM77.22903. and in part by the U . S. Office of Naval Research. Fluid Dynamics Branch. It is a pleasure to acknowledge fruitful discussion with B. Brumley. and the computational help of D. Hatziavromidis . 123

Copyright @ 1778 by Academic Press. Inc All rights of reproduction in any form reserved. ISBN 0-12-002018-1

124

John L Lnmley

I. Introduction

A. HISTORY AND GENERALITIES At the 1968 AFOSR-IFP-Stanford Conference on Computation of Turbulent Boundary Layers (Kline et al., 1969),all the methods formally presented were oriented specifically toward the turbulent boundary layer, and most involved some variant of the mixing length formulation. During the Monday afternoon discussion, however, C. Donaldson made an informal presentation (reported on pp. 114-118) in which he strongly supported an idea he was then developing, which he called “invariant modeling.” He said “. ..it has been my feeling that one should keep track of the dynamics of all the second-order correlations of importance.. . . one seeks to express those terms which are unknown ... in terms of the second-order correlations themselves.. . . In making these choices one is constrained only by the requirements of symmetry and the general conservation laws.” And in another place, “Using general tensor notation . .. one.. . seeks.. .the simple invariant form that will reduce to one of the forms generally seen in traditional boundary-layer studies.. . .” Later in the conference, an ad hoc committee was formed to deal with the importance of invariance (consisting of Bradshaw, Donaldson, and Mellor). This committee’s brief report (p. 426) states in part: “Invariance of a particular formulation to coordinate transformation is important . .. particularly .. . if one is attempting to describe threedimensional flows.. . . If one chooses models which are of invariant form, it should be found that these models have the greatest generality; they should have the highest probability of describing turbulent flows which depart from the particular geometry for which the parameters in the model were adjusted to agree with experimental data.” This was the birth of the technique which has become known as “secondorder modeling” (as well as invariant modeling, and which the French school refers to as the one-point closure). During the last decade it has undergone very rapid development in the hands of numerous authors whose work will be mentioned explicitly in the appropriate sections below. The development at the present time is by no means complete, and the present work must be regarded as an interim report. One thing is already clear, however; in many situations of practical importance this technique makes possible computations which often agree with what data are available. Inevitably, the technique is also being applied in many situations in which data do not exist, which must be regarded as a dangerous practice, since the limitations of the technique are not known with any precision. It is primarily the possibility of practical computation which has been responsible for the great interest in this method.

Computational Modeling of Turbulent Flows

125

The method is not without historical antecedents. A good description of the early development of these ideas is contained in Monin and Yaglom (1971, p. 318). Specifically, Kolmogorov (1942) appears to have been the first to suggest characterizing turbulence entirely by its intensity and scale and using this to simplify the equations, an idea which is also used by all authors. Chou (1945a,b, 1947) suggested a number of closure schemes, and in particular was the first to use the equations for the third moments, eliminating the fourth moments by various hypotheses, something which we shall discuss later. Finally, the suggestions of Rotta (1951a,b) and Davidov (1958, 1959a,b, 1961) for the modeling of the pressure-strain correlations and related matters have had extensive influence on the development. of this technique. It is not an exaggeration to say that there is little in use at the present time that was not suggested by these authors. These authors all predated the easy availability of large-scale computers, so that their suggestions, for the most part, could not be explored extensively. Second-order modeling, even in its most stripped-down form, results in general in the simultaneous solution of four partial differential equations in the domain of interest; more elaborate models in a three-dimensional situation might require the simultaneous solution of as many as 36 partial differential equations to obtain the mechanical field only. Fortunately, this is within the capabilities of present computers at a reasonable price, which cannot be said of any other technique. Direct simulation (Orszag and Patterson, 1972) is limited to relatively low Reynolds numbers, but this is not serious, since the large scales (which are responsible for transport) are dominated by inertia, and thus are essentially independent of Reynolds number; however, if there is no homogeneous direction in which averages may be taken, several hundred realizations must be generated to obtain stable statistics, which is prohibitively expensive. There are also problems with initial conditions (Lumley and Newman, 1977): problems of differencing errors in current codes restrict the initial conditions on turbulent structure to fairly unrealistic ones. Almost the entire computational time is used in setting up a realistic turbulence, by which time the mean initial conditions have already changed. Direct simulation is thus not an alternative for practical computation. The various sophisticated closures (Leslie, 1973) suffer from essentially the same problems as the direct simulations and hence are also limited to homogeneous situations. Thus, the second-order modeling is at present the only possibility for practical computation. Second-order modeling may also be said to have as antecedent the work of the school of Rational Mechanics, which had its roots in the work of Stokes (and to a lesser extent Navier), and the modem development of which is associated with the name of Truesdell and his co-workers (with, of course, many others unnamed in between). In early work on non-

126

John L Lamley

Newtonian fluids, closures were developed for flows in particular geometries, relating the stress to the deformation state. Typifyingsuch closures is the power-law fluid, in which the shear stress is assumed proportional to a power of the local strain rate. This is not an invariant formulation and cannot be easily generalized to other geometries and flows. The school of Rational Mechanics, leading to the work of Coleman and No11 (1961), adopted the point of view that it was necessary to discover first the general form for a general geometry and flow, which the dependence of stress on deformation could take, considering the various mathematical and physical restrictions to which it was subject; this would permit the identification of a relatively small number of invariant functions, which could then be determined by experiment in particular geometries, the results being applicable to flows in all geometries. This is practically a statement of Donaldson’s point of view, quoted at the beginning. In fact, the techniques used for the computation of turbulent boundary layers at the Stanford Conference are reminiscent of the situation in continuum mechanics before the advent of rational mechanics; in just a few decades this approach has revolutionized nonNewtonian fluid mechanics, and we may hope that a similar approach will do the same for the computation of turbulent flows. Regarding rational mechanics, it is fair to say that there are those who feel that it hrts been philosophically productive but has not been greatly useful for computations; that in its general form it is too complex, and there are many more fluids in nature than can be easily encompassed by it; and that it is consequently necessary to fall back on the old empiricism for many practical computations. Just so, in turbulence computations, there are those, typified by Bradshaw (personal communication), who feel that the behavior of turbulence is so complex that the search for general closures is probably futile and that practical computations will require empirical techniques developed for the specific geometry. There are also those who feel that the general second-order modeling produces forms too complex to be of use in practical situations. My position is not diametrically opposed to these. Rather, I would say that I believe in the ultimate possibility of developing general computation procedures based on first principles; and under certain circumstances I believe that it is possible to do this rationally by the techniques of second-order modeling. While it may be necessary to fall back on empiricism for computations in complex situations, I believe that rational second-order modeling can at least provide a guide for the construction of the more empirical models and can certainly serve in general to indicate the range of applicability of these techniques. Quite apart from its utility as a practical computational tool, I have found that the attempt to devise models of this kind which behave like turbulence brings to light aspects of turbulent behavior that might never have been noticed. That is, one often finds that experiments do qualitative things which

Computational Modeling of Turbulent Flows

127

models cannot be made to do, no matter how the constants are adjusted. This indicates that a basic physical mechanism has been omitted from the model; if it can be identified and included (often difficult), something fundamental will have been learned about turbulence. This, incidentally, shows that Gauss was not quite correct when he remarked that, given seven constants, he could produce an elephant on a tightrope, and that with nine, he could make it dance.

B. GENERAL ASSUMPTIONS The development of second-order models has proceeded on a somewhat ad hoc basis, the degree depending on the predelictions of the particular author. Although nearly everyone recognizes the desirability of having all terms be tensors of the appropriate rank, with the correct symmetry and other properties (such as the vanishing of traces where appropriate), there has been little consideration of turbulence dynamics beyond this, and almost any convenient quantity with the proper tensor properties has been fair game, without consideration for its behavior at large or small Reynolds numbers, large or small anisotropy, etc. We will try to develop here fragments of a rational approach, from which we will see that the second-order models appear to be an orderly expansion about a homogeneous, stationary turbulence, the large scales of which have a Gaussian distribution. In fact, homogeneous turbulence is observed to be approximately Gaussian in the large scales (Frenkiel and Klebanoff, 1967a,b),even in the presence of homogeneous distortion (Marechal, 1972); what is envisioned here is a turbulence which is made non-Gaussian in the large scales by inhomogeneity, but which would, on the removal of the inhomogeneity, relax to a Gaussian state. The expansion about the homogeneous, stationary state suggests that the ratios of the turbulence length scale to the length scale of the mean flow inhomogeneities and of the turbulence time scale to the time scale of the mean flow evolution are both small. This is essentially a kinetic theory type of approximation. It is known experimentally that these ratios of scales are not small in real turbulence, being in general of order unity (since the turbulence structure and the mean flow inhomogeneity are generally produced by the same mechanism, unlike the artificial situation in a wind tunnel, where a homogeneous turbulence and the associated homogeneous mean shear may be carefully produced by different mechanisms). It is legitimate to ask why one might expect second-order modeling to resemble real turbulence; that is, why one has a right to expect the first term in an expansion in a small parameter to be applicable when the parameter is of order unity. There are many other examples of this phenomenon: for instance, the first two terms in an expansion in small Reynolds number for laminar flow around a cylinder work well up to Reynolds numbers of order 10 (Van

128

John L Lumley

Dyke, 1964). That, of course, is not an explanation. The explanation here probably is this: by following a rational procedure we have created a physically possible phenomenon, not quite real turbulence perhaps, but one which conserves momentum and energy; transports the right amount of everything budgeted (momentum, energy, Reynolds stress, heat flux, etc.), although not by quite the right mechanism; satisfies realizability (so that nonnegative quantities are never negative, Schwarz’ inequality is always satisfied, etc.); behaves correctly for both large and small Reynolds numbers; and reduces to real turbulence in one limit (weak inhomogeneity and unsteadiness). Probably any mechanism that satisfied all of these restrictions would behave about the same. The physical input from experiment is essentially used to fix the amount of transport. Of course, it is also possible, as Donaldson suggests (personal communication), that there is some other basis on which these equations can be derived, under which they have broader applicability. A given set of equations can often be derived from different sets of hypotheses of different degrees of generality. An example is the equations for global energy of a disturbance to Couette flow, which can be derived both from the exact equations and from the equations in the small disturbance approximation (Lumley, 1972).We only claim here to have found a consistent basis. We have already seen that we are making a specific assumption about the probability density of the turbulence; we are also making an assumption about the spectrum: it can be parametrized by the large scales. That is, if we know the characteristics of the large scales and the Reynolds number, then the shape of the rest of the spectrum follows; in fact, all of the statistics of the turbulence are determined by the large scales and the Reynolds number. This suggests two things: the turbulence has had time to come to equilibrium spectrally with the large scales and hence changes in the large scales are slow enough for the small scales to follow, and boundaries and initial conditions are far enough removed to have no direct influence on the present state. Otherwise, these would, in general, introduce another parameter (or several) which should be included.

11. Mathematical Preliminaries

A. REPRESENTATIONS In what follows, we will use a number of ideas familiar to workers in continuum mechanics, but perhaps not so familiar to workers in other areas. The principal concept that we will need is that of the form which an isotropic

Computational Modeling of Turbulent Flows

129

tensor function of other tensors may assume. Many references exist (e.g., Lumley, 1970a), but none covers the subject in quite the way we will need it. Thus, we will give a brief introduction here. Suppose that we have an isotropic symmetric second-rank tensor function of a symmetric second-rank tensor, both of zero trace c$ij(bpp), say. This form will occur when we consider the return to isotropy. We may consider that bij and 4ijare average second-order turbulence moments. First, by an isotropic function we mean that the functional relation is isotropic; that is, it is not dependent on any other suppressed variable, such as the direction of a magnetic field or the axis of a wind tunnel. The values of 4ijare not isotrop ic, but any anisotropy in 4ijis induced by anisotropy of bij;if bij is isotrop ic, 4i, is isotropic. A general technique for determining the dependence of 4ijon bij is to select two arbitrary vectors A , and Bj and form an invariant 4 i j A iBj. This is now a tensor of zeroth rank, having no free indices, and must consequently be the same in any coordinate system, hence, an invariant. As an invariant, it must be a function only of invariants, those that can be made from A , , Bj, and bij. We now must ask what are the independent invariants of this collection of quantities. From the Cayley-Hamilton theorem (Lumley, 1970a, Section A2.6) we know that in three dimensions only the Oth, lst, and 2nd powers of bij are linearly independent since bij must satisfy its own secular equation:

b;

- lb$

+ IIbij - I l l b ; = 0,

(2.1) where b$ = 6,, Kronecker’s delta, and b$ = bikbkj,etc. The quantities I, 11, and 111 are, respectively,

+

I1 = (biibjj - b,Z,)/2, 111 = (biibjjbkk - 3biib;j 2b;)/3! (2.2) Note that bif.is the trace of b$ and not biib,. These quantities (2.2) are the only independent invariants of bij. In addition tofhese, we have the invariants of A , and B,, which are the lengths and included angle of this vector pair, or equivalently A i A i ,BiBi,and A i B i . We have in addition, the invariants that can be constructed between the vectors and the tensor; what is invariant here is essentially the orientation of the vectors relative to the principal axes of the tensor, which should give no more than six quantities. These can conveniently be gotten from the collection A i b i j A j ,A,b$A,, Ai bijB,, A , b$ Bj, BibijBj, B, b: Bj. Higher powers of bijare not independent by the Cayley-Hamilton theorem. We are assuming, incidentally, invariance under the full rotation group; that is, we include invariance under improper rotations, or reflections, presuming that there is no preferred spin to the turbulence relations (not to the turbulence!). We are also working in Cartesian coordinates; the generalization to non-Cartesian systems is relatively straightforward but adds complexity. Note that if bij were not symmetric

I = bii,

130

John L Lumley

there would be other independent invariants, since b , b,, # b, b,, and A, bijBj # Ai bji Bj. Now, 4ijAi Bj is bilinear in Ai Bj, and thus the right-hand side must also be since these vectors are arbitrary. We may thus exclude all forms quadratic in Ai and Bj, and set 4 i j A iBj = a(I, II, III)b$ Ai B, + @(I, II, IIl)b, Ai Bj + y(I, II, III)bi Ai Bj. (2.31 We may now remove the A, and Bj and obtain

4ij= adij + @bij+ y b i ,

(2.4)

where a, @, and y are unknown scalar (invariant) functions of the invariants I, II, and III. If we now include the condition that b, has zero trace, I = 0, I1 = - b i / 2 , and I I I = b i / 3 ; if we take account of the fact that 4ii= 0, then a = 211~13so that

4ij= fibij + y(b$ + 211dij/3).

(2.5)

This form is essentially completely general, granting only that there are no suppressed variables. We have now reduced the determination of the dependence to the determination of two unknown scalar functions of scalar arguments, independent of flow geometry. These considerations can be extended to more complex situations. For example, suppose that chi, were a function of b , and a vector, ci. The invariants of the vectors with each other and with the tensor would now be AB, Ac, Bc, cc, AbB, Abc, Bbc, cbc, Ab’B, Ab’c, Bb’c, cb’c in an obvious notation. The general expression would not only have terms proportional to AB, AbB, Ab’B, but also proportional to AcBc, AcBbc, AcBb’c, BcAbc, BcAb’c, etc., making all possible pairs bilinear in A and B leading to 12 terms in all. The invariant functions would depend on cc, cbc, and cb’c in addition to the previous terms. It is clear that we are including in this way more invariants than we need (more than are independent) since we begin with 18 (two from the set A, B, and c, and one from the set bo, b, and b’), while only 15 appear to be independent (from a consideration of which ones are necessary to determine the projections of the three vectors on the principal axes of bij).It is possible to construct a reduced basis (Lumley, 1970b)which considerably simplifiesthe forms; however, in most cases we will consider if the turbulence is isotropic, b will be proportional to Kronecker’s delta, and the reduced basis will no longer span the space. Thus, it is better to be somewhat redundant, and include extra invariants, in order to have a basis which is valid in all cases.

Computational Modeling of Turbulent Flows

131

B. REALIZABILITY Recently, Schumann (1977) introduced the concept of realizability. He pointed out that no matter what equation is used to predict values o f w , it must have the property that does not allow correlations to exceed the limiting values dictated by Schwarz’ inequality and does not permit negative component energies. One might argue that this is a mathematical nicety so far as the component energies are concerned (although it is clearly important in regard to the correlations); that is, one might feel that the occurrence of vanishingly small component energies is so rare in nature that requiring correct behavior of model equations in that situation is unnecessary. However, anyone who has had experience computing with second-order models knows that a frequent cause of aborted calculations is the occurrence of negative component energies, often associated with recovery from poor initial conditions. There are, in addition, real physical situations in which one component is strongly suppressed, notably in stably stratified buoyant turbulence. Hence, it is physically important and computationally convenient, as well as mathematically nice, to require realizability. The most convenient way of satisfying the requirement is to transform the equation to principal axes of uiuj; now there are no off-diagonal terms; we may designate the component energies in this coordinate system (the eigenvalues) by We now require that if were to vanish, its time derivative too would vanish, so that it arrives at zero with zero slope and cannot cross to negative values. Just how this is implemented, we will see in Section I11 when we discuss the return to isotropy. Note that it was erroneously stated in Lumley and Newman (1977) that ,is foolproof only the method just presented, satisfying realizability for if the direction of the principal axes is not changing with time. That this is not true can be seen by differentiating the eigenvalue relation. That is, if we write at every instant

z.

v,

then differentiation of the second relation and substitution of the first gives: BijX$k)+ ~ X k f=)U:k,Xlk)+ %Xik).

(2.7) If (2.7) is now multiplied by XIk)and if the second equation of (2.6) is used, the second and fourth terms, containing the time derivative of the eigenvector, cancel, leaving -

(2.8) and if it is required that the right side of (2.8) is to vanish when the eigenvalue vanishes, then the condition is satisfied.

,i(k) - X !I k ’ B . . X y , IJ

132

John L h m l e y

This general concept of realizability has many other applications, which have not been explored before the present work. For example, the equations for intrinsically positive quantities such as the dissipation of energy F must also have the property that the time derivative will vanish if the quantity vanishes. Further, we look at the equation for ?,the total turbulent energy; if ;TI vanishes, the time derivative of must vanish. This is a slightly different requirement from Schumann’s realizability in practice, although included in it in principle. That is, we will implement Schumann’s realizability by requiring that the time derivatives of u& vanish if the component vanishes, even if the other components do not. The sum of the three time derivatives, of course, must be equal to production minus dissipation; if vanishes, the requirement that correlation coefficients not exceed unity will assure the vanishing of the production. However, we must separately assure that the dissipation vanishes when vanishes to guarantee that the time derivative vanishes. Roughly speaking, this means that, as vanishes and if cdoes not, the time derivative of B must go to - 00 to assure that it will also vanish at the same point. The same consideration applies to the temperature variance. Finally, we must consider correlations between unlike quantities, like velocity and temperature, which are also constrained by Schwarz’inequality. These must ke treated by a slightly different technique, because the two do not form a tensor. However, we may still examine the eigenvalues of the correlation matrix. These eigenvalues are nonnegative and vanish in only three cases: if the variance of either variable is zero, or if the correlation coefficient becomes unity in absolute value. If we write the secular equation for the matrix and require that the time derivative of the eigenvalue vanish if the eigenvalue vanishes, we obtain as a condition that the time derivative of the determinant will vanish if any of these conditions is met:

a

a

a

__

-2abub

a

+ > p + 23 = 0

(2.9) for two arbitrary correlated quantities u and b. If we arrange to have 2 vanish when 2 vanishes and the same for b, Eq. (2.9)will be satisfied if ab also vanishes; if the correlation coefficient cannot exceed unity (in absolute value), this will be assured. Hence, we must consider the case of the correlation coefficient approaching unity. If p is the correlation coefficient, we have A _

*

p / p = ab/ab - 2 / 2 2 - $/2p.

(2.10)

Substitution in (2.9) permits us to write

--

2pp = ( 1 - p2)(ue2/a2 + b2/b2). -L-

(2.11 )

Hence, if the variances are vanishing, but are not yet zero, and the correlation coefficient rises to unity (i.e., the correlation is not vanishing fast enough), the time derivative of the correlation coefficient must vanish, and

Computational Modeling of Turbulent Flows

133

the equation for the correlation must be constructed in such a way as to assure this. This type of realizability, with regard to velocity-temperature correlations, has been completely neglected until now. In the case of the correlation between a scalar and a vector, such as temperature and velocity, in addition to the constraint on the heat flux along each axis, we may obtain a more general condition:

__

( G j P- e u i e U j ) A i A2j o

(2.12)

for an arbitrary vector A i . The conditions on the individual components are included in this and can be obtained by taking Ai successively parallel to the three axes. However, (2.12) is more general. The easiest way to implement the condition is to transform to principal axes of the tensor in parentheses in (2.12); then a statement equivalent to (2.12) is that all eigenvalues of that tensor must be nonnegative. In these coordinates then, we may write equations like (2.9) in each coordinate direction and assure ourselves of satisfying the general condition (2.12), since each eigenvalue is in the form of the usual Schwarz’ inequality. Again, the development of (2.6-2.8) is applicable, and the fact that the directions of the principal axes are changing in time is irrelevant.

111. The Return to Isotropy A. INTRODUCTION It is part of the folklore in the field that anisotropic turbulence tends to become more isotropic; that is, in the absence of other influences the components interchange energy so as to become more nearly equal, and this exchange process proceeds faster than the decay. The experimental evidence for this is individually not very strong; Fig. I, from Lumley and Newman (1977) is typical. Nevertheless, the evidence taken together is incontrovert-

10

15

20

25

jn

35

4n

45

50

55

i M

FIG.1. The return to isotropy of homogeneous turbulence following a contraction. Measurements of Uberoi. (Lumley and Newman, 1977).

John L Lumley

134

ible: such a return to isotropy does exist, but it appears to be very slow for weak anisotropy. This return to isotropy is produced by the pressure terms. Let us consider a homogeneous flow, without mean velocity gradients. The equation for the Reynolds stress is:

(m+ m)/p- 2VUi.k

(34 Although it makes no difference in this flow, it is customary to remove the transport effects (Tennekes and Lumley, 1972,Section 3.2) from the pressure term, writing =-

-(m+ p,iui)/p = -[(puj).i + (pui),jl/~+

uj,k*

+ ui,j)/p*

(3.2) Although Lumley (1975b) has pointed out that this separation is not unique, and (3.2) may not be the most appropriate choice, there are reasons of convenience (which will become apparent in Section VI) for making the separation as in (3.2). We concern ourselves with this matter here because we intend to use the modeling which we devise for this pressure term also in inhomogeneous situations. It is evident from the presence of the strain rate ~ ) the trace vanishes; that is, the term serves only to interin ~ ( u+~u ,~ ~, that change energy among the components, not to create or destroy energy. The viscous term in (3.1), although its primary function is to dissipate uiuj stuff to heat, can also cause interchange of energy among the components at any Reynolds number. Although the term is observed (Monin and Yaglom, 1975, p. 453) to become more isotropic with increasing Reynolds number, in agreement with Kolmogorov’s (1941) hypothesis of local isotropy, if u1 = 0 (which can occur conceptually or computationally, and even approximately in reality) then U1.k U1.k = 0 also. We will add and subtract the trace, which is twice the dissipation of energy: J4Uj.i

2 v u x = [ 2 v w k - 2Edij/3] + 2E6ij/3.

(3.3) The deviatoric part (in square brackets) now acts to interchange energy among components, but neither creates nor destroys total energy. If we define -E4ij = p(ui,j uj,i)/p - 2~2Zdij/3 (3.4) then Eq. (3.1) may be written as u,ul = -E4ij - 2Edij/3 and 4ijis dimensionless, has zero trace, and is solely responsible and responsible only for the return to isotropy. The tendency toward equipartition which returns the turbulence to (or toward) isotropy can also be regarded as a decay of Reynolds stress, since it is only a question of being, or not being, in principal axes of the Reynolds stress tensor. It seems natural to consider also the decay of other correlations, notably the heat flux. There is very little evidence for this decay; Fig. 2 is the only known data, reproduced from Warhaft and Lumley (1978a).

+

+

Computational Modeling of Turbulent Flows

135

XlM

FIG.2. Decay of temperature-velocity autocorrelation coefficient in grid-produced turbulence after Warhaft and Lumley (1978a; see also Newman et a/., 1978a). Curves were produced by the analysis of Section 1II.C.

Fortunately, the data show incontrovertible decay of the correlation coefficient. In the same homogeneous flow, without mean velocity or temperature gradients, the equation for the heat flux may be written as

iiJ

= --/p

- (v

+ y)B,iui.i,

(3.5) where y is the thermometric diffusivity. Although it is not at all customary, we will separate the first term on the right-hand side just as we did (3.2):

-DIP = -@),i/p

+ KIP.

(3.6) This is not ordinarily done, since the term obtained (the first on the right) is not in the form of a transport term, although it is a divergence (a transport term should be of the form (q),i for the turbulent transport of a quantity A). However, realizability requires that this term be separately modeled. To see this, suppose that all other terms in the equations are realizable. Then, in

John L Lumley

136 --

--

principal axes of f12uiuj - flui flu,, applying Eq. (2.9), realizability requires -

--

when Pu,ui = flu, flui. This is clearly true for the unmodeled terms (in their natural form as given in (3.7)). More to the point, however, since we will be modeling the term on the left for inclusion in the Reynolds stress equation, we must include a model of the term on the right in the heat flux equation, and the model must satisfy (3.7). Let us define _-

-4!jflujElq

-

2 - 3 -

- P ,iIp - ( V

so that Eq. (3.5) becomes simply -

~ i =8

-

+ Y)fl,jui.j

-- &,flu, E/q2.

(3.8)

(3.9)

This is not a completely general form, but is sufficiently general for our purposes. The molecular transport part of (3.8) will become essentially zero for high Reynolds-Peclet numbers under local isotropy, but that need not concern us. We also need the equation for the temperature variance:

f12 = -2E0.

(3.10)

We can now apply realizability to these equations. First, if we transform to principal axes of ~(iu~, and take = 0, then the right-hand side of (3.4) must vanish, so that we require 411= -4 (since the dissipation and -other _ components _ are notzero). Next, -~ we transform to principal axes of f12uiu, flui fluj and take pu,uj = 8 ~ fluj. 1 Then, applying Eq. (2.9), we require

Note that, if in addition the one-heat flux vanishes, all components of the Reynolds stress containing u1 vanish; this is equivalent to the vanishing of the one-eigenvalue of the Reynolds stress tensor, and hence by the realizability requirement for the Reynolds stress, the second term on the left of (3.11) must also vanish. From this, if the one-heat flux vanishes, the right -side of (3.11) must vanish. So @ must , have the same principal axes as f12uiuj 0%. BU, . Thus, it must be a function only of that tensor, although it may be a function of invariants of other tensors, parameters, etc.; for, if it were a function of any other vectors or tensors, the principal axes would not be the same.

Computational Modeling of Turbulent Flows

137

B. THEREYNOLDS STRESS Now let us consider dij.If Eq. (3.4) is considered as an equation for q5ij, it is clear that it is determined by the history of the Reynolds stress and the present value of the dissipation. We could write in general

4ij = 4ij(uVj, c v},

(3.12)

that is, a functional of the past values. The kinematic viscosity has been included for generality, although it is not explicitly present in the equation. Time does not appear explicitly, since we anticipate that the relationship will not change with time (although the quantities involved will). The functional will involve integrals over time, however, and the most satisfactory normalization_is_to define a new time d.r = dte/?; this is equivalent to T - to= f In (q2/qo2).This new dimensionless time is monotone in the true time since the turbulence is decaying. The argument of (3.12) contains eight quantities involving two dimensions; hence it is possible to form six independent dimensionless groups. We choose to do this in such a way as to make the role of anisotropy most evident; we will define an anisotropy tensor (3.13) This vanishes identically if the turbulence is isotropic, it is dimensionless, symmetric, and has zero trace; hence, it has five independent components. For the sixth group, we take a Reynolds number R, = (?)'/~Ev.

(3.14)

The factor of nine has been included so that, if we take ? = 3u2, B = u3/l, (3.14) becomes the traditional Reynolds number based on this length and velocity. Hence, (3.12) can be written as a functional of (3.13) and (3.14). Now, we may presume that the turbulence has a fading memory and will ultimately forget its initial state; if changes in the mean state are sufficiently slow relative to the turbulence time scale, (3.12) can be expanded in a functional Taylor series about the present state (see Lumley and Newman, 1977, for details); keeping only the first term, (3.12) reduces to a function of the present values of (3.13) and (3.14). The ratio of the turbulence time scale to the time scale of change of the mean field is almost always of order unity in real turbulence, so such an expansion is not formally justified; however, we will find that the resulting expressions work remarkably well. Thus, (3.15)

John L h m l e y

138

This is now of the form discussed in Section II,A, and we can immediately write the form (2.5)

dij = B(11, Ill, R,)bij+ ~ ( 1 1I,l l , RJ(b$ + 211dij/3).

(3.16)

We must now determine the forms of /3 and y in their dependency on the invariants and on the Reynolds number. (This y is, of course, not to be confused with the thermometric diffusivity.) First, it is instructive to consider how turbulence can be characterized in terms of the invariants. A plot of all turbulence on axes of 11 and I l l is shown in Fig. 3. (Note that our definitions (1.2) are the classical ones and differ from those used in Lumley and Newman (1977) by simple factors: 11 = - 11'12, 111 = 111'13, where the prime indicates the expressions used in that paper.) In Fig. 4, the possible region of variation of the eigenvalues of bij, say bl and b, (since the trace is zero, b3 = -b, - b,) is shown. Each eigenvaluecan be no smaller than - 3, corresponding to the vanishing of that component, and can be no larger than +$, corresponding to the vanishing of the other two. Thus, b , and b, are constrained to lie within the triangle in

-II

(2/27.1/3) . 1D

AXISYHMETRIC

ISOTROPIC 2D

AXISYMMETRIC

I -0.05

0

0.05

m

01

FIG.3. Possible states of turbulence parametrized by the independent invariants of the anisotropy tensor after Lumley and Newman (1977). Turbulence must occur within the region (or on its boundaries) delimited by the axisymmetric and two-dimensional states.

Computational Modeling of Turbulent Flows

139

FIG.4. The possible range of variation of the two independent eigenvalues of the anisotropy tensor. Turbulence must occur within the triangular region, or on its edges, which correspond to the two-dimensional state. The maxima and minima of 111 on the curve of constant I1 are indicated, where this curve crosses the lines corresponding to axisymmetric states.

Fig. 4. A constant value of I1 corresponds to the ellipse; if I1 is small in absolute value, the ellipse lies entirely within the triangle; if 1 I1 I is larger, it lies partly outside, and the parts outside are excluded (that is, they do not correspond to possible values of bl and b2).At the largest possible value of -11 = 3, the ellipse just touches the corners, corresponding to the three possible one-dimensional turbulences. On the segments of the ellipse lying within the triangle it is possible to determine the total derivative of 111 with respect to b, or b2 for fixed 11; it is readily found that relative maxima occur where the ellipse crosses the lines joining the apexes to the origin, corresponding to axisymmetric turbulence. Relative minima occur where these lines cross the other sides of the ellipse, when the ellipse is small enough to lie entirely within the triangle, also corresponding to axisymmetric turbulence. The relative maxima correspond in each case to two components being equal and the third being larger than the average, while the minima are the opposite. The minima disappear when the component that is smaller than the average reaches zero. Moving away from the maxima, the value of I11 decreases monotonically until an edge of the triangle is reached, corresponding to two-dimensional turbulence. In the two-dimensional and axisymmetric cases, there are simple relations between I1 and 111 which may readily be found. In the axisymmetric case, if

bij = ( 0b -b/2 0 0

8)

-b/2

, (3.17)

John L Lumley

140

then I1 = -3b2/4 and 111 = b3J4; hence 111 = +2(-11/3)3/2. In the twodimensional case, if bij is placed in principal axes and bl = -f the expressions for I1 and 111 can be readily reduced to give

+ 3111 + I1 = 0.

(3.18)

These lines are shown on Fig. 3. By inspection it is clear that a stronger statement can be made relative to the quantity in Eq. (3.18):

6 + 3111 + I1 2 0

(3.19)

everywhere. Hence, the vanishing of (3.19) can be used as an indicator of two-dimensionality of the turbulence, and we will avail ourselves of that possibility. First, Lumley and Newman (1977) have shown from the data of ComteBellot and Corrsin (1966) that for axisymmetric turbulence in vanishingly small anisotropy,

41 Jb = 2.0 + 8.0/R,?12.

(3.20)

The value of 2 corresponds to no return to isotropy. This can be seen by forming the equation for bij: (3.21) dbijldr = - (4ij - 2bij) so that if q511 = 2b in the axisymmetric case, the value of bij remains unchanged. Hehce, (3.20) suggests that there is no linear return to isotropy at an infinite Reynolds number; the effect is thus either a viscous or a nonlinear effect. By forming the equations for I1 and I l l , and considering very small anisotropy, Lumley and Newman (1977) have shown that y must vanish as the anisotropy vanishes at infinite Reynolds number. We have an additional requirement. The so-called final period of decay (Batchelor, 1956) describes turbulence at a Reynolds number so low that the nonlinear terms may be neglected. Since they are responsible for the return to isotropy, there is consequently no longer any return to isotropy. Consequently, the right-hand side of (3.21) must be zero for a zero Reynolds number. We may now attempt to satisfy these various requirements in the simplest way possible. We will be able to satisfy the conditions with y = 0, giving a relatively simple expression. [Lumley and Newman (1977) give a complete expression, but the form is cumbersome.] If 4ij= /3bij, then by (3.21) we must have /3 2 2 always, since f l < 2 would correspond to the spontaneous increase of anisotropy in the absence of any external agency. There is no proof that this should not happen, but it seems unlikely. Since the expression (3.19) is always positive and vanishes if one of the eigenvalues of bij takes the (i.e., the turbulence becomes two-dimensional), it is tempting to value

-4

141

Computational Modeling of Turbulent Flows

take fl = 2 + F(R,, ZI, Ill)($+ 3111 + 11).This would automatically satisfy realizability for the Reynolds stress. It has the disadvantage that, if the turbulence became two-dimensional, there would be no return to isotropy of the remaining two components. However, it is not our intent here to make a model that behaves properly for two-dimensional turbulence, which is fundamentally different in so many respects. We are interested only in a simple, workable form which satisfies realizability, and since a form & j = C1bij is widely used in the literature (see, e.g., &man and Lumley, 1978), the form proposed is likely to be satisfactory. The unknown function F must be determined so that it vanishes at a zero Reynolds number, takes on the value 8.O/R:l2 for small anisotropy and a large Reynolds number, is always positive, and otherwise fits the existing data. The form

fl = 2 + e~p[-D/R;/~](72/Rj"~

+ A h[l + B(-ZZ + CZZZ)])($ + 3111 + ZZ),

A = 80.1,

B = 62.4,

C = 2.3,

(3.22)

D = 7.77

fits the data given in Lumley and Newman (1977) nearly as well as the more complicated expression given there (Fig. 5 ) ; the agreement with the form

1

2

3

4

5

6

7

8

-ux10*

FIG. 5. The return to isotropy + , , / b (the function from the analysis of Section II1,B) for an axisymmetric flow with negative 111. The light curves are the more complete description of Lumley and Newman (1977), while the heavy curves are the analysis of Section II1,B. The experimental points are those reproduced in Lumley and Newman (1977).

142

John L Lumley

given by Lumley and Newman (1977) is equally good for the axisymmetric case with 111 > 0, for which no data exist. The behavior in Fig. 5 displays the various characteristics we have discussed. Note particularly the fact that the curve arrives at the twodimensional state (- 211 = $) with finite slope which indicates that this state is unstable; if perturbed, there will be a finite rate of return to isotropy. C. THEHEATFLUX

We may now consider +!, beginning from the realizability condition (3.11). Let us introduce a tensor

__

(3.23) D,, = iqqa - eu, euj/P?. This has nonnegative eigenvalues that vanish only if the correlation is perfect in that component. From the realizability condition, we have already , be a function of Dij. It is straightforward but tedious to found that &must show that (in principal axes of D i j ) if D l l = 0, then 2

4 + 3111 + 11 = 3DzzD3,8U1 /@?.

(3.24)

Substituting this in (3.22), and that in (3.11), we obtain for realizability in principal axes of Dij with D l l = 0: =

1

+ r + (F/2)($ + 3111 + 11 -

D22D33),

(3.25)

where F is the expression in (3.22) which multiplies the last set of parentheses. We may identify D Z 2D3, as ll,, the second invariant of the tensor Dij under these conditions. r is the time scale ratio, r = (~/P)(42@). Thus, to satisfy realizability, we could assume = [I

+ Y + (F/2)(3 + 3111 + I I -

llD)]dij

+ gDij + hD$,

(3.26)

where g and h are functions of the invariants of D,,, I,, ll,, and 111,. If we consider decay of the correlation coefficient in an isotropic turbulence, we find, if p is the correlation coefficient,

I, = -p(f$&

We find 11, = 2( 1 - p2)/9 + f and

- r - l)E/?.

- 1 - r = (F/2)(3- 11,)

+ gD11 + hD:l.

(3.27) (3.28)

With D1 = ( 1 - p2)/3 finally we have d In p/dz = -(8/3 - F/9)(1 - p’) - (h/9)(1 - p2)’.

(3.29)

Excellent agreement with the data of Warhaft and Lumley (1978a) is obtained by taking h = 0, g/3 - F/9 = 1.6 (the curves in Fig. 3 were obtained

Computational Modeling of Turbulent Flows

143

using these values). Until more data are available indicating variation with the Reynolds number, etc., we will use these values. We have a choice; the F appearing in (3.29) is F with I I = I I I = 0, say F,. In defining g in general, we could use either the general F, or just F,. Lacking any information, we will take F because it produces the simplest equation. If we form the equation for the decay of the correlation coefficient of an axial heat flux in an axisymmetric turbulence, we obtain with this choice d In p/dr = -(B - 2)/2 - 3bll/(l

- (1

+ 3b11)

+ 3bll)(l - p2)(1.6+

b11 F/6)

(3.30)

so that the decay rate is increased if b l l is positive (and vice versa). This must await experimental confirmation but has a certain appeal. The axisymmetric form (3.30) is applicable to the conditions of &man and Lumley (1976), where a buoyancy-driven atmospheric surface mixed layer was considered. The heat transfer is entirely vertical, and the turbulence is axisymmetric about the vertical axis. Using Zeman and Lumley’s (1976)fixed value for B = 3.25, we find for the right-hand side of (3.30) a value of -2.44, whereas the value using their equations, which assume a fixed value for 4; = 7.0, is - 4.79; the value we obtain for 4; = 4.65 under their conditions. It seems very likely that the larger value was necessitated by the failure of their rapid terms to satisfy realizability conditions, driving the heat flux to impossibly high values, and requiring a larger return to isotropy term to kekp the heat flux under control. IV. The Rapid Terms A. INTRODUCTION

In what we have done so far, we have considered only homogeneous situations without mean velocity gradients or buoyancy. If we begin from the Navier-Stokes equations for the fluctuating velocity in an incompressible situation (in the Boussinesq approximation, if buoyancy is importantsee Lumley and Panofsky, 1964, Section 2.1),

+ ui,juj+ ui,juj + ui,juj -

+

+

- p , ~ p pie vui,jj and take the divergence, we obtain an expression for the pressure: tii

=

+

(4.1)

-v2pip = 2ui,juj,i ui,juj,i- pie,i. (4.2) This contains two types of terms on the right-hand side. The first and third are linear in the fluctuating quantities, while the second is quadratic. The second term is the nonlinear scrambling (Hanjalic and Launder, 1972) or the

144

John L Lumley

return-to-isotropy term; it is the pressure field associated with the mixing of the turbulence by itself generated by this term which is responsible for the gradual equalization of the energy in the various components of the turbulence. It is the other terms which we wish to examine here. We may conveniently split the pressure into two parts, defining

The correlations with p") are those which have already been modeled in Section 111. The correlations with p") and its gradients are those which we refer to as the "rapid" terms. The term "rapid" is used for two reasons. First, there is a classical problem in turbulence called the rapid distortion problem (Batchelor and Proudman, 1954), in which one imagines that the turbulence is subjected to a velocity gradient so intense that for some little time one may neglect the transport and distortion of the turbulence by itself, and hence neglect p(') and the nonlinear and viscous terms in Eq. (4.1). Since the pressure p") is the only pressure present during this rapid distortion, it is natural to refer to it as the rapid pressure. Second and more important for us, however, it is not difficult to show that, if an initially isotropic turbulence is subjected to a sudden distortion, not necessarily rapid, all correlations in the equations not involving p(l)begin from isotropy and gradually develop anisotropy in response to the distortion; the same is true for the sudden application of a gravitational field. The terms involving p(l), however, are instantly anisotropic. This also justifies the use of the term rapid, and in addition makes clear why it is necessary to model these terms with considerable care; in nearly every situation, they exert a very strong influence on the structure of the anisotropy of the velocity field (see, e.g., Townsend, 1970). Incidentally, the concept of rapid distortion has been extended to the case of the sudden application of a very strong gravitational field (Gence, 1977)within the same approximation. The first of Eq. (4.3) may readily be solved by standard techniques to give p") explicitly in terms of the right-hand side. Of the many solution techniques available, we will use Fourier transforms, which are appropriate to a homogeneous field. Although there are, in general, terms in the various correlations with p") arising from inhomogeneity, consideration of which would require a solution technique for p(') that did not suppose homogeneity, these terms are never considered, and there is no indication that the homogeneous forms are insufficient in any practical situation. For an introduction to the use of Fourier transforms for homogeneous stochastic fields, see Lumley (1970a). Indicating the Fourier transform off by [f], we may write

Computational Modeling of Turbulent Flows

145

To obtain the various correlations that appear in the equations, we must now multiply by the Fourier transform of the appropriate quantity, average, and integrate over the Fourier space. Taking the simplest case first, let us multiply by [el* (the complex conjugate), average, and integrate over Fourier space, designating S the spectrum of temperature variance, s, the spectrum of the heat flux, and s,,, the spectrum of the Reynolds stress for future reference. In all cases the integral of the appropriate spectrum over the entire Fourier space gives the quantity in question. We will indicate an integral over three-dimensional Fourier space by 1dK.

Taking the second part of the integral only, without the buoyancy parameter, we will consider

j” (KqKi/K2)SdK= Iqi,

say.

Now, Iqi must clearly be symmetric. If i is summed on q, the term in parentheses goes to unity and the integral is, by definition, the temperature variance, so that we have as a condition I,, = p.In addition, if the turbulence is isotropic, S is a function only of the magnitude of the wavenumber vector, and hence is spherically symmetric. If the integral is carried out in two steps, first over spherical shells, and second, over the radius, the term in parentheses can be integrated immediately to give

on a spherical shell. If the integral is now carried out over the radius, we have for the isotropic case We must develop a model for Iqi which satisfies these restrictions. Most authors use the isotropic value for Iqi for all situations. However, we must also satisfy the realizability conditions. Writing Iqij and Iqijk for quantities constructed like Zqi in (4.6),but using S j and S , and applying the realizability condition (2.9)to the equations for the heat flux, taking in account the fact that the realizability conditions must be separately satisfied by the terms multiplying the buoyancy parameter since it may be given any value and orientation (the same being true of the terms multiplying the mean velocity gradient), and using the fact that the spectra and hence Iqi,etc., may be

John L Lumley

146

(say, by unfelicitous initial_conditionsJ-from _ fl, and U i S jwe , obtain (in principal axes of the tensor p u i u j - 8ui 8 u j )

manipulated separately

&I,,

=

F Z , ~ , , GI,,= PI,,,,

G I =, ~ FPz,~~ (4.9)

as the conditions which must be satisfied if either the corresponding velocity component vanishes, or the corresponding correlation coefficient achieves the (absolute) value of unity, or the temperature variance vanishes. If we can assure ourselves that the corresponding condition (4.9) is satisfied when the correlation coefficient achieves unity in absolute value, then by (2.11) the correlation coefficient will be bounded. Hence, if either the temperature variance or the corresponding velocity variance vanishes, the corresponding heat flux will vanish. If we arrange that Iql vanishes when the one-heat flux vanishes, and similarly for the other components, then (4.9) will be satisfied under all conditions. Thus, we must (1) arrange for I , , , to vanish if vanishes, etc., and (2) must satisfy the corresponding member of (4.9) if the correlation coefficient in question becomes f 1. Wyngaard (1975) and &man and Lumley (1978) looked at the horizontally homogeneous case of a vertical temperature and velocity gradient and considered a strong stable stratification which annihilated the vertical velocity and heat flux. The equation for the vertical heat flux in these circumstances is -

ew =

-83

+ p3P.

(4.10)

Since it seems likely that the vertical heat flux should remain zero under these circumstances, they were led to require that 1 3 3 = F as the vertical velocity vanishes. As a general prescription this is, however, not consistent with (4.9), if we assume that the value of lqiis determined entirely by the temperature and velocity fields as we shall. For then, if the velocity and heat flux vanish in any direction, the diagonal component of Iqi in that direction must take on the value p,while the other two must vanish, since lii= p. The condition is immediately violated if two components of the velocity vanish. In fact, an exact condition can be obtained from (4.9). If the horizontal velocity and heat flux vanishes, then the first two of (4.9) are satisfied; if the correlation coefficient in the vertical is unity, then 1 3 3 must vanish, since Ipii= 0 by continuity, and hence 1,33 = 0 since I,, and I,,, are zero (since 8ul and = 0) (see,for example, Lumley, 1970a, Section 4.7). Hence, 1 and I , , cannot both be zero. This exact condition was also invoked by Wyngaard (1975) and &man and Lumley (1978). We may add one final condition: the diagonal components of lqimust be nonnegative, since the temperature spectrum is nonnegative. In the case when the correlation coefficient achieves the value of 1, the conditions (4.9) have a simple interpretation. If u , and 8 are perfectly cor-

,

,,

Computational Modeling of Turbulent Flows

147

related, they are proportional, and consequently the spectra S, S , , and S I 1 are all proportional. Thus, we see that, not only must we construct a model for Iqi consistent with all these conditions, but we must model I,, at the same time, since their limiting values must be coupled. The equivalent conditions on Iqij are I4V. . = I.W J..’ I44J. = 8u.; J I W.. = 0 ; and Inij -,0 as -,0. These have all been discussed before or are completely analogous to conditions on Iqi. We have, in addition, the conditions on Iqij which arise from realizability of the Reynolds stress. We may obtain these conditions by placing the equations for the Reynolds stress in principal axes; if the intensity in, say, the onedirection vanishes, then for the vanishing of the time derivative we require that I,, vanish also (since the one-heat flux vanishes) for arbitrary p . Thus this is not an independent condition.

6

B. THEHEATFLUX INTEGRAL In the past it has been customary to take for Iqij a simple linear combination of G.For example, &man and Lumley (1978) and nearly all other authors use

Iqij = (3)[Si,euj

-

(+)(SijK + Sqj&)].

(4.11 )

This satisfies the first, second, and third of the conditions in Iqij, but it does not satisfy the important condition that it vanish when the j-heat flux vanishes, both required by realizability and by the fact that S j vanishes when the j-heat flux vanishes. This almost surely causes computational difficulties (Zeman and Lumley, 1978). We could attempt to remedy the situation by constructing a very general formalism, expressing Iqij as a function of bijand the heat flux, forming invariant functions, etc. Fortunately, it is possible to construct by inspection a simple expression that at least satisfies all the necessary conditions : 1,ij

= (6,i

-- _ _ -8~,0~i/tl~,8~,)8~j/2.

(4.12)

C. THETEMPERATURE VARIANCE INTEGRAL

The construction of an expression for Iqi which satisfies all requirements of its own and is properly related to the expression (4.12) by realizability conditions is not quite so simple, and we have been unable to find a completely general expression. The following expression has nonnegative eigenvalues, has the proper trace and the proper value under isotropy, but

148

John L Lumley

satisfies realizability exactly only on the total heat flux, not on the components: I 1.4 = eZsi,/3

+ Phi, - (3/2a)(G&,

-

__

-.e~,e~,s,/3).

(4.13)

To see in what senserealizabilityis satisfied, suppose we are in principal axes of the tensor @uiuj - flui 8 u j , and suppose that the correlation is perfect in the one-direction. Then, because off-diagonal terms of the tensor and because -- the leading diagonal term also vanvanish in principal axes, ishes, we may write @u,u, = Bu, Ou, ;substituted in the expression for bqi, this may then be combined with the other terms in (4.13), producing

__

-- _ _

I,, = (eu, eu,/2q2)(s,1 - eu, eu,leu, eu,).

(4.14)

At the same time we have from (4.12) the expression

-- _ _

I,, = (8U,/2)(sq1- eu, eu, /eu,eu,).

(4.15)

These have the same dependence on q ; the realizability condition produces --

eu,eu,

=

(4.16)

which is true only if all components are perfectly correlated (in which case all three realizability conditions will be satisfied exactly). This can occur approximately in free convection, in which the horizontal velocity essentially vanishes; htnce, the expressions (4.12) and (4.13) may be useful in that situation. When 24, = 0, I , , = Ou, 8u,,/2q2; if the horizontal heat flux is small, this will be essentially zero, as is desirable (see Eq. (4.10)).No computations have been done yet with (4.12)-(4.13), and we must regard them as an interim solution since it seems likely that additional terms may be required. Lumley (1975a) suggested that perhaps these rapid terms should be functions of the mean field parameters, that is, the buoyancy parameter and the mean velocity gradient. However, Reynolds (1976) quite correctly pointed out that the Iqi,etc., are determined by the form of the spectra; if the spectra are set by initial conditions in a particular state, the Iqi, etc., will be completely determined independent of the buoyancy parameter and the mean velocity gradient. These will exert their effect as the field evolves, of course, but only by their influence on the form of the spectra. Hence, in determining forms for these terms, we exclude all quantities other than those characterizing the state of the turbulent field, the second-order correlations. Lumley (1975a) also suggested that the expressions for Iqi, etc., must be linear in the second-order quantities. The reasoning was as follows: the rapid terms are those which would appear in rapid distortion theory. Rapid distortion theory is a linear calculation, and therefore results are superposable. We may model rapid distortion theory by using our rapid terms, neglecting return to isotropy and dissipation; hence, our modeling results must also be

Computational Modeling of Turbulent Flows

149

superposable. Hence, the expressions for lqimust be linear in the secondorder quantities. This position seems to have found acceptance (Reynolds, 1976). However, it must be incorrect in some sense, since it is impossible to satisfy realizability with linear expressions. This is a perplexing problem, which deserves further thought. The answer could, of course, be that it is not possible to parametrize even approximately the terms lqiin terms of the second-order quantities; however, I feel that that is unlikely. A more appealing explanation is the following: in rapid distortion theory the solutions are dependent on the time multiplied by the mean velocity gradient due to the linearity; since geometrically different mean velocity fields produce different solutions, the solutions can be expressed as functions of dimensionless structural parameters, constructed from the mean velocity gradient (Lumley, 1975a); due to the way in which time appears, these structural parameters may, in each solution, be replaced post hoc by expressions constructed from the solution variables themselves, producing an apparent nonlinearity. Hence, a general expression, parametrized in terms of the second-order quantities, in the course of a rapid distortion problem, would have certain combinations of these quantities which remained constant, playing the role of structural parameters. It has also been suggested (Lumley, 1975a)that the rapid distortion problems could be used to calibrate the modeling of the rapid terms. While it is possible to construct rapid terms which model the rapid distortion problem quite successfully (Gence, 1977), it now seems likely that this is not a useful procedure. The problem is, that the spectra, for the same values of the large scale parameters (the second order quantities) will, during rapid distortion, have quite different forms from their usual ones, since the strain rate, or buoyancy, is supposed to be sufficiently strong to dominate the nonlinear effects. Thus, the strain rate, or buoyancy, will be felt directly in the small scales of the spectra, which will not have an equilibrium form. It is conceivable that the rapid terms could be modeled incorporating a parameter representing the ratio of the time scale of the distortion to the time scale of the dissipative eddies; during the distortion this is essentially zero, while for the cases we are interested in modeling, it is essentially infinite. One could envision calibrating the model against rapid distortion with the parameter set at one value, and afterwards setting it at the other value. However, this would be a complicated and delicate procedure and seems unlikely to be worth the effort. Since even the extremely crude models for the rapid terms currently in use (which do not satisfy realizability and sometimes not even incompressibility) give results which are in many ways satisfactory, it appears probable that we can, by devising models which do satisfy these conditions, produce results which are completely satisfactory without going to extremes.

John L. Lumley

150

D. THEREYNOLDSSTRESS INTEGRAL We may now consider the rapid terms which multiply the velocity gradient. From Eq. (4.5) we may evidently write the appropriate terms in the heat flux equation as h i

=***

-

- Ui,kUkd + 2Up,,1piq,

(4.17)

and in the Reynolds stress equation as 2up,q1piqj+ 2upsqIpjqi - ui,km - U j , k m ~ i . (4.18) Again, if we rotate to principal axes of the Reynolds stress, and consider the vanishing of, say, the one-component, then for the time derivative of the a dpq. We have used the fact that one-component to vanish we require lplql the velocity gradient U i , j is arbitrary relative to the principal axes of -the Reynolds __ stress and that Ui,i= 0. If we now rotate to principal axes of 02uiuj - 8ui Buj , and apply Eq. (2.9), supposing that the one-eigenvalue, say, vanishes, i.e., that the correlation is perfect in the one-direction, we obtain

uiuj =

* * *

-

-

~21p,q1 = 8% 1,1,

+ Adpq,

(4.19)

--

where A is arbitrary, and use was made of the fact that P u i u, = Bui Bu, ,the mean velocity gradient is arbitrary, and Ui,i = 0. We have in addition the various symmetry and other conditions on lPw:lpiei= lpijq; lpiqj = lipqj; lppqj = ii&(which we will term normalization); lpiij = 0 (incompressibility). We can add, finally, that if the turbulence is isotropic, we must have (Crow, 1968; Rotta, 1951a) ~

-

1P4V.. = (4dijdpq - dpid, - dpjdqi)q2/30

(4.20)

obtained by integrating the isotropic form of the spectrum (Lumley, 1970a). The realizability conditions above are new, but all other conditions have been known for some time. The conditions other than realizability can all be satisfied by a linear combination of terms in bij: Iijpq

+ (APijbpq - (i+)(aipbjq + diqbjp + djpbiq + djqbip)l + (?/11)[5dijbpq - (dipbjq + d i q b j p + d j p b i q + djqbip)I

= C[bijapq

+ (?/30)[46ijdpq - dipdjq- diqdjp).

(4.21)

This form was independently devised in different ways by a number of different workers (for a full discussion, see Launder, 1975). The constant C is not determined by the various conditions, and is usually determined empirically. However, different workers obtain different values, when different flows are used for calibration. This is probably because the form (4.21) does not satisfy realizability for any value of C, as can be easily checked. If we

Computational Modeling of Turbulent Flows

151

,

transform to principal axes of b,, suppose that b , = -4,and form Iilplrwe find that it is zero for i # p. However, for i = p # 1 we find ZPlp1 = -bp,(3C

+ a ) / l l + C/11 - ?/330

(no sum)

(4.22)

while for i = p = 1 we find

l l l I l= -C/11

+ 2?/55.

(4.23)

Thus, there is no value of C which will satisfy realizability. Depending on the value which is chosen, this term will either oppose the growth of the offdiagonal terms of the Reynolds stress, not permitting them to grow large enough in some geometries, or permit them to grow too large. In case the other two components of the velocity are equal, we could take C = 2?/5, which would then satisfy realizability, since the diagonal components would then be equal. Launder (1975) suggests a value of 0 . 2 7 5 a (our C = c2 2 / 2 in Launder’s notation). The situation with regard to realizability relative to the heat flux is even worse. There is essentially no relation between IP,,,and It is evident that the failure of realizability is due primarily to the terms of the form 6 , bip, etc., in (4.21), since on setting j = q = 1, we are left with bip,which cannot be made to behave properly. Hence, the entire group of four terms (required by symmetry) must be excluded. It is possible to satisfy all symmetry, incompressibility, normalization, and isotropy conditions with the following expression :

This does not satisfy realizability, however (relative to the Reynolds stress), since in principal axes of bij, if b , , = -$, we have Iilpl

= (?/30)(36i,6p,

- dip).

(4.25)

This, at least, has only delta functions, which could be canceled by the addition of another term. We must create a function to add to (4.24)which is symmetric in i - j and p - q, which has zero trace when i is summed on j and when j is summed on p, which vanishes under isotropy, and which, in principal axes of b,,, if b , = -4, will produce only terms like those of (4.25) so that a coefficient can be chosen to cancel them. Note that the d i p can remain (since the whole expression will be multiplied by U i , ,which is incompressible). It is necessary to go to third order in order to obtain such a term (note that, to be sure of having only terms like those of (4.25)we must never have bi, appearing):

,

bib,,

+ (11/3)bijhPq+ (211/3)6ijb,, + (2Z11/10)6ij6pq- (3111/10)(6ip6jq+ 6jp6iq).

(4.26)

152

John L Lumley

This vanishes under isotropy, has both traces zero, and in principal axes of bij, if b,, = -4, reduces to (taking the ilpl component) dil dP1(-& - 11/3 - 111/10) - (3111/10)dip.

(4.27)

We can satisfy realizability relative to the Reynolds stress by adding (4.26) to (4.24) with a coefficient 27?/( 10 + 9011 + 27111). If we examine the realizability relations relative to the heat flux, we again find that it is not possible the_ general relation, but it is possible to satisfy the relation for the to satisfy_ case of 8ui 8ui = pa,when all components are perfectly correlated (or only one component exists) by the combination of (4.24) and (4.26). Note that there are eight independent expressions which all satisfy the requirements of symmetry, incompressibility, etc., any of which can be multiplied by (4 + 3111 + 11) and added to the combination of (4.24) and (4.26) without condition _ on either heat flux (approximately) or affecting the realizability_ Reynolds stress, since if 8ui 8ui = then 4 + 3111 + 11 = 0. These expressions are forbiddingly complicated, but where is it written that turbulence must be simple? The principle is simple; the result is an odd tensor polynomial on a bounded domain. The complexity is primarily so that the expression will say the same thing regardless of the geometry. It should be emphasized that these expressions are by no means the last word, especially since they have not been tested by computation. They are predicated on the simple form for Zijk, which is, in some respects, the easiest to model. A more complex form could be constructed for Zijkr with constants which could be adjusted to conform to experiment, and which would then lead to more complex forms for Z i j and Zijkl. The failure to satisfy realizability exactly relative to the heat flux may be serious, and it may be necessary to construct these more complicated forms for this reason. These interim forms at least point the way toward improvements on the forms in current use, which do not satisfy realizability in any sense.

e'z

V. The Dissipation Equations A. THEMECHANICAL DISSIPATION The equations for the mechanical and thermal dissipation are in the sorriest state. Whereas realizability and considerable physical reasoning has given substantial form to the equations for the second and third moments, when the equations for the dissipations are reduced by high-ReynoldsPeclet number assumptions (Lumley and Khajeh-Nouri, 1974) what is left is

B + E , i ui +

io+ Eo,i ui+

= - (2/?)+, = - (EoZ/P)Jlg,

(5.1)

Computational Modeling of Turbulent Flows

153

where the dimensionless invariant functions on the right-hand side contain the entire mechanism for the production and destruction of the dissipations. The situation is particularly difficult, because while the dissipation occurs at small scales, it is fixed by large-scale mechanisms; what we are really trying to determine is the rate at which energy or temperature variance is passed along the spectral pipeline to scales at which it can be attacked by molecular transport (see Tennekes and Lumley, 1972, Section 3.2). If the molecular diffusivity is changed in a turbulent flow, it simply changes the scale at which the dissipation occurs, but does not change the level of the dissipation. Since the spectral transport is related to the skewness of the velocity differences at two points (Monin and Yaglom, 1975, Section 21.4) for intermediate separation distances, it would probably make sense to try to model the equation for this quantity, as suggested by D. C. Leslie (personal communication). In the meantime, we must do what we can with what we have. If the Reynolds stress equation is regarded as an equation for the dissipation, we may probably conclude that I(/ is a functional of the history of the Reynolds stress, the heat flux, and the dissipation, as well as the mean velocity gradient, the gravity vector, and possibly the viscosity. Making the same assumption regarding slow changes, we can write II/ as a function of these quantities:

*

*(w,G,E, ui,j,

(54 We cannot exclude the mean velocity gradient and buoyancy vector as we did previously (in other equations) because it is not at all clear that (5.2) is determined solely by the state of the turbulence at a given instant. It is certain, for example, that if the buoyancy vector is zero, a heat flux cannot affect the mechanical dissipation. It also seems clear that the sudden imposition of a mean velocity gradient cannot have any immediate effect on the rate of change of dissipation if the turbulence is isotropic; on the other hand, additional sources of production in the energy equation must cause changes in the level of dissipation. We can thus conclude that the heat flux and buoyancy vector must occur as a product and the mean velocity gradient probably cannot occur without the anisotropy tensor. We may select a collection of nondimensional groupings in (5.2), as =

pi,

-

v).

$ = $(blj, q2ui,jfGGPjh RI)* (5.3) This must now be a function of the invariants that can be constructed from these quantities. Unfortunately, they are numerous, particularly because the mean velocity gradient is not symmetric. Pope (1974; see also Spencer and Rivlin, 1959, 1960) has examined a related problem. A complete list includes 22 invariants, even with the restrictions we have placed on the problem. It is clearly impossible to determine the dependence of a function of 22 variables by any imaginable set of experiments. The list is usually restricted to include only quantities of first degree in the mean velocity gradient and the buoy-

John L Lumley

154

ancy vector on the principle that these may be regarded in the first instance as weak. Such a restricted list is: -

$ = $(ll, I l l , bijq2Ui,j/E,b;TUi,j/E, & f i i / E ,

Gbijfij/E, &bifij/E, R,). (5.4)

It is also consistent with the position that the mean velocity gradient and the buoyancy are weak to expand (5.4) in a series in these quantities, keeping only linear terms: $ = $0

+

+

-

$1

bijq2Ui,j/E $2b;TUi,j/E

+ 4b3&Bi/Z + $4&bijfij/E + t,bsGbifij/E,

(5.5)

where the coefficients are functions of I I , I l l , and the Reynolds number. Even these are not all used. For example, Zeman and Lumley (1978) use $0

$3

= 3.8

+ 6011/(1 + 3( -211)1'2),

=

+0.95,

=

-3.8,

= 20/( 1

+ 3( -211)'12).

(5.6) These forms are determined on a largely empirical basis; the forms in the denominator of $o and t+b4 are found necessary in computation of the buoyancy-driven atmospheric surface mixed layer, which becomes very anisotropic near the inversion base. The leading constant is the only one that can be determined cleanly: the data of Comte-Bellot and Corrsin (1966) give $4

= 3.78

- 2.77R;

'I2

+ 18.1811 +

(5.7) for large Reynolds number and very small anisotropy (Lumley and Newman, 1977). The correlation coefficient for the Reynolds number variation is0.8, while that for the variation with anisotropy is 0.33, and hence this coefficient should not be taken too seriously. If we consider turbulence without production from either mean velocity gradient or buoyancy, then there are two limiting cases for which values of the leading coefficient can be determined: in the final period of decay, for very small Reynolds number and all anisotropies, Lumley and Newman (1977) determine that $o = they also obtain the same result for one-dimensional turbulence and arbitrary Reynolds number, i.e., for If = -3. In the absence of any other information, they suggest a simple interpolation formula relating these values: $0

v;

+ 0.980 exp[ -2.83R; '/'I[

1 - 0.33

In( 1 - 5511)l.

(5.8) At observed values of anisotropy, this does not differ much numerically from the form used by Zeman and Lumley (1978). $o =

Computational Modeling of Turbulent Flows

155

It has been known for some time that computation schemes of the type discussed here, when calibrated against plane flows, will not predict the development of axisymmetric flows properly, as discussed by Pope (1977). Specifically, the spreading rate of a round jet is overestimated by about 40% (Pope, 1977). Pope suggests that this is due to the omission in (5.5) of a term of the form (Ui,j - Uj,i)(Uj,k- Uk,j)(Ui,k+ Uk,i). The spectral transfer is produced by the stretching of large-scale vorticity (Tennekes and Lumley, 1972, Section 8.2), and this term is intended to represent this process; in particular, the term vanishes in two-dimensional mean flow, when the mean vorticity is normal to the plane of the mean strain rate. In axisymmetric flows, however, this is not true, and the term is nonzero. The inclusion of such a term restores the correct spreading rate for the axisymmetricjet. I feel, however, that the term suggested cannot be quite right, because it is independent of the anisotropy of the turbulence. That is, if the turbulence is isotropic, it should make no initial difference what the geometry of the mean velocity field is. Another possibility presents itself immediately: in developed turbulence, to a zeroth-order approximation, one can set approximately bij oc (Ui,j+ U j , i )(Lumley, 1967), and this can be substituted in the expression proposed by Pope. Evaluating the expression in a plane shear flow, we fmd that it will vanish if b,, (normal to the plane of the shear) vanishes, which it approximately does in plane deformations. This suggests other possibilities that have the same property, in particular 111; it seems more rational that +o should depend on 111 than on Pope’s (1977) suggestion, but this must be determined from a comparison of the decay of isothermal anisotropic turbulence without mean velocity gradients in the axisymmetric case, when 111 # 0, and in the case of plane deformation when it is approximately zero. It is unfortunate that the data are so poor; the determination of t,b from experimental data is equivalent to determining the second derivative of the data, and very little data are sufficiently accurate to permit this. We can adduce one rigorous condition, which is unfortunately of no help in determining the form of our equation. For realizability, it is necessary that e vanish if 7 vanishes. Hence, the F equation must be arranged so that, if vanishes without F, will become infinite so as to drive B to zero. As it stands (5.9, or even (5.3), has this property. We may introduce one further possible criterion. At the edge of a wake or jet, or thermal plume, the production closely balances the dissipation as both vanish (Tennekes and Lumley, 1972, Section 4.2). We have already considered the necessity for assuring that the dissipation vanishes when the velocity vanishes; here, however, we must say something stronger to assure that they vanish together at the proper rate. If they are to vanish at such a rate that (?)”/E is bounded for some n, and production and dissipation are of

+

156

John L Lamley

the same order, then we should have

zz = nz?

= 2nE(P, - e),

(5.9) where we have written P , for the total production, and (5.9) should be valid as both energy and dissipation vanish and production is of the order of dissipation. All the terms in Eq. (5.5) will be of the same order; hence, this suggests that $?, $4, $5 should not be present, and that -4b3 = = $o = 2n. Since +~I! is at least 3.78, n = 1.89 or larger so that the dissipation will go to zero faster than the energy as observed (Tennekes and Lumley, 1972; Section 4.2). Elimination of 4b4 would have caused difficulty in the calculation of Zeman and Lumley (1978); however, what was needed was a reduction of the dissipation production in the case of nearly onedimensional turbulence; this could be achieved in other ways, say, by the introduction of a factor of 4 + 3111 + 11. The equality between and -$3 we will see again when we discuss the equation for the dissipation of temperature variance. B. THETHERMAL DISSIPATION We can now consider the second part of (5.1), the equation for the dissipation of temperature variance. Proceeding in the same way, we obtain the same set of variables, with the addition of Ee, p,and These allow the addition of several invariants and a ratio, r = ie;fz/esE, the time scale ratio, which can appear as a variable in all the coefficients. If we consider the case of decaying anisotropic turbulence without production of any kind, we have

_ -

---

I+P= $ o e ( ~ 111, ~ , r, R,, 8U,i&$FTTI euibijeuj/iPq2,euieujb;/iPz), (5.10)

Some authors avoid the entire problem by assigning a constant value to r (Spalding, 1971). However, it seems likely that r may vary depending on local conditions of transport and production. We may form an equation for r, obtaining d In r/dr = r(2 - i,hoe) - (2 - $o). (5.11) Extensive measurements by Warhaft and Lumley (1978b) indicated that, when the turbulence was isotropic, there was absolutely no tendency for r to change during the decay; it could be set at any value initially and would maintain that value until the end of the tunnel. This suggests taking, at least for the isotropic case, (5.12) i,ho0 = 2 - (2 - i,b0)/r. It also suggests that any tendency for the ratio r to tend toward an equilibrium value must be due to anisotropy and/or the production terms, which

Computational Modeling of Turbulent Flows

157

are producing temperature and velocity fluctuations from gradients by the same mechanism. We can write by analogy with (5.5) (5.13) $* = *oe + $18zi&3,i/Eg, where we have suppressed the other terms, taking the position (at least temporarily) that the temperature equation should not depend directly on the buoyant production or mechanical production (Newman et al., 1978a). If we form the equation for r, we find (designating the mechanical production by P,, the buoyant production by P b and the thermal production by P,) d In r/dz = i - ( t / ~ ~-* 2jP,/&

+ (2 - $l)Pm/E+ (2 + $3)Pb/E- r$$ +

(5.14) where we have designated by $: and the departures of $oe and $o from their isotropic values. Let us consider first an anisotropic decaying turbulence without production. If r has reached its equilibrium value re, which must be determined later, then we expect the time derivative to vanish, so that we must have re $$ =

(5.15)

If it is assumed that $$ is not a function of r, (5.15) fixes its value. Now consider an equilibrium situation with production, but without buoyancy so that P, = 0. We expect that P , = E, and P, = E ~ If. r has reached its equilibrium value re, we must again have the time derivative vanish, which gives (5.16) re($1*- 2 ) = $ l - 2. If we now introduce buoyancy, and let P , + P b = i, we find - $ 1 = $ 3 . Although their forms were somewhat different from those we have assumed, the above relations (with the exception of - $ 1 = t+b3) were approximately satisfied in the work of Zeman and Lumley (1978). These relations provide for an orderly relaxation of the time scale ratio whenever there is anisotropy or production. We obtain, making use of these various relations, d In r/dr = (2 - $,)[(P,

+ P&

-

+

(P,/EB)r/re] $o'(l - r/re).

(5.17)

The fact that the equation for dissipation of temperature variance does not contain, with these assumptions, terms quadratic in the heat flux, is not serious. Newman et al. (1978a) found from an analysis of the data of Alexopoulos and Keffer (1971) that the inclusion of such a term was not warranted by the data. We must now determine a value for re. For vanishing anisotropy and infinite Reynolds-Peclet number it is clear that the number should be near

158

John L. Lumley

unity. A test field model simulation (Newman et al., 1978a; Newman and Herring, 1978) suggests the value 1.0. However, the measurements of Warhaft (Warhaft and Lumley, 1978b) suggest a value of r = 1.34 when the three-dimensional velocity and temperature spectra peak at the same wavenumber, which one would expect in a real equilibrium situation. For a vanishing Reynolds number and all anisotropy, we have the final period of decay. Corrsin (1951)has shown that then re = 0.6. For large anisotropy we may consider one-dimensionalturbulence. Newman et al. (1978a)consider a somewhat special type of one-dimensional thermal turbulence which is isothermal in the direction of the unique velocity component, and take a limit from the not-quite-one-dimensionalcase (because the limit is singular); they find that re has the same value for all Reynolds numbers as it does in the final period of decay. Having only three limiting values of re,it is reasonable to construct an interpolation formula among them, the simplest possible. Newman et al. (1978a) take the variation to be of the same type as was assumed in Eq. (5.8). If the value of 1.34 is taken as the value for re at infinite Reynolds number and vanishing anisotropy, then we have (5.18)

i+bo - 2 = 4r,/3

for all Reynolds numbers and anisotropies, since the limiting values of the two quantities are then related. Putting these various relations together, we obtain

he = 2 - (2 - h 0 ) / 1 + 4 ( h - h 0 ) / 3 ( h - 2), = 2 + 4(Ic11 - 2)/3(h - 21,

(5.19)

where $oo is the isotropic value of @o. The following forms were used by Zeman and Lumley (1978): t,hoe = 3.0(1 =

__

__

+ r/4) - 3ohi&~/r2e2 q2,

+0.97.

(5.20)

Distinction among these various forms will require more experiments on isotropic decay with heat flux and anisotropic decay with heat flux, followed by the same flows with mean velocity and temperature gradients. From the work done so far, it is certain that terms involving the mean gradients of temperature and velocity must be present in the equations, as well as the terms involving anisotropy. Note that we have a realizability requirement on the dissipation of temperature variance: if the temperature variance vanishes without the dissipation, V ,I must become infinite to force the dissipation to zero. However, this condition is satisfied.

Computational Modeling of Turbulent Flows

159

C. THETRANSPORT TERMS In the following section we will discuss the transport terms in the equations for the Reynolds stress and the heat flux. However, the transport terms in the dissipation equations d o not fit in the same picture, because it is not possible to write exact equations for the transport of dissipation, as it is for the transport of velocity and temperature variance, etc. Hence, we will discuss these terms here. We may again apply the somewhat tentative condition that as the dissipation and the velocity variance vanish, they should d o so in such a way that the left member of (5.9) is satisfied. When we wrote (5.9) we had in mind a homogeneous situation, but now the transport terms must also be included. Thus, we must require

+ 2p/p)uk],k = ($b(q?u,),k

=~[(qz

?(Ek),k

(5.21)

as the velocity variance and the dissipation vanish, where we have made use of the homogeneous approximation for the pressure transport (see Section VI). We may satisfy this requirement by replacing throughout in the expression for &k (see Section VI) ( w j ) , k

(5.22)

=wjE,k/M

42

which supposes that b,, remains finite as vanishes; if we then write where q2ukis the expression with the replacement (5.22), the condition (5.21)will be automatically satisfied. This also has the advantage of being independent of the value of n. For example, in the nonbuoyant case we obtain from (6.54)

mk = (3%/5?)&k,

m k

= -(+)(?/&)[3/(10 -(:)Eipu,/u,/q2

+ 4C,)][&,p(u,p + 2 W k w / ? ) -(

~ ~ , ~ p u p / q z ] . (5.23)

Exactly the same consideration can be applied to the transport of Ee. If we and Ee vanish together so that (e')n/Ee remains bounded, presume that then one should require

n(82u,),kEe= @(m),k. This may be satisfied by replacing in the expression for

(eUJj =eu,EB,i/2nZe,

6

= e'ze,j/nEe

(5.24) (5.25)

and calling the resulting expression (p/rEe)EE. In the nonbuoyant case this results in ___ --= - (e'/~~)[j( 1 C ~ / ~ ) ] [ E ~ , ~ euj ( We u k / e z )- ( + ) ~ ~ . ~ e ~ ~ (5.26)

+

+

e~~/e~

John L Lamley

160 VI. Transport Terms

A. INTRODUC~ION In the equations for second-order quantities, in inhomogeneous flows transport terms appear which are of third order. For example, in the equafor the transport term is ),k; in the equation for it is tion (ui0uk),k. These represent the divergence of the flux of the quantity in question, produced by the fluctuating velocity. In many flows, these terms represent the principal source of energy, heat flux, etc.: for example, in the atmospheric surface layer driven primarily by surface heating, the erosion of the inversion base, and gradual thickening of the surface mixed layer, is due entirely to these terms. To close the equations for second-order quantities, we must somehow express these terms as functions of the second-order quantities. If we approach this problem from the point of view of invariant modeling (Lumley and Khajeh-Nouri, 1974), we would say that, in the nonbuoyant case, u,ulu, must be a functional of B and (ignoring the Reynolds number dependence). Since the third moment must vanish in the case of homogeneity, it must depend on gradients of the arguments of various orders. It is possible to write an expansion, supposing weak anisotropy and inhomogeneity; to first order this gives:

(w

w,

q,

where a and b are unknown constants. This is clearly unsatisfactory,since it makes no distinction among the components. To second order, six more constants are introduced. Each of the fluxes (say of temperature variance or heat flux) will contain as many more undetermined constants. From a practical point of view, there are not enough well-documented experiments to unambiguously determine all these constants. In addition, it is not clear that second order is high enough. We have already seen, in connection with realizability, that it is not generally sufficient to consider anisotropy small in any sense. What is needed is a physical model which will produce a form for the fluxes, consistent with our general form, but with all the constants determined. Classically (see, e.g., Tennekes and Lumley, 1972) a simple mixing length or gradient transport argument was used in a simple shear: (q2/2 + P/p)O

- vT a(?/2)/ay

(6.2) with no attempt made to consider tensor properties. The eddy viscosity used is the same as that obtained (not unambiguously) from the ratio =

Computational Modeling of Turbulent Flows

161

-iiB/(dU/dy). This type of model is currently used in invariant form (Lewellen, 1975):

6= -c(F),jqEjqz/E,

(6.31

where the constant still must be determined from comparison of predictions with some class of experiments. The ?/E may be replaced by some other scalar combination of the correct dimensions, depending on the author’s preferences. With regard to uiuiu,,it was recognized early on (Hanjalic and Launder, 1972) that there was no justification for writing -

- ( ~ ) , ~ q 2 / ~

uiujuk

(6.4)

since the three velocities have equal standing and should be treated equally. Hence, we should write something like uiujuk

(mu,),* +

-(T/c)[(wj),lm +

(vk),m]. (6.5)

The constant of proportionality was optimized by Launder e f al. (1975) at 0.055. Many authors find this forbiddingly complicated, particularly in more

complex flows. Since there is not much theoretical justification for it, and since the constant still must be determined from experiment, many prefer to use a simpler form such as (6.4). For example, Daly and Harlow (1970) use (6.4); Launder (1975) recommends a value of the coefficient of 0.125. Hanjalic and Launder (1972) tried to provide justification for the form (6.5). If the equation for ~ i ~ isj written, ~ k it consists of the substantial derivative; production terms of the form U i , , , T and its permutations; and its permutations; dissipressure gradient terms of the form - p . i / u , l p pation terms of the form - vui,,,uj,,, uk - vui,,,uj uk,,, and its permutations; and finally, on the left-hand side a collection of terms of the form:

(m),, -

-

wj(mU, ),l wk(m),I u j ( q I ) , I -

(6.6)

Hanjalic and Launder (1972) suggested the following set of assumptions: neglect all substantial derivatives, production, and dissipation terms; since any rapid terms arising from the pressure gradient correlations will be of the same form as the production terms, neglect them also; replace the pressure gradient correlation with a relaxation term proportional to -

- ui ujuk q2/c

(6.7)

with an undetermined coefficient; in the expression (6.6) introduce a quasiGaussian assumption (Lumley, 1970a): ___. -- -UiUjUkUl

= uiujukul ~~

+

uiukujul

+ ujukuiul.

This collection of assumptions leads to the form (6.5).

(6.8 1

162

John L Lumley

The same set of assumptions was applied to the buoyant case by Zeman and Lumley (1976) with excellent results. Forms for the fluxes were produced which can be identified with definite physical mechanisms (Lumley et al., 1978). The relative success of this technique suggests that the assumptions may be justifiable in some way; that is, there may be some single guiding principle from which they all follow naturally. We can ask a number of questions and make several comments which may help to shed light on this: homogeneous turbulence is observed to be Gaussian in the energy containing range (Frenkiel and Klebanoff 1967a,b) even in the presence of nonzero velocity gradients (Marechal, 1972); departure from Gaussian behavior, therefore, is associated with inhomogeneity (which is clear because nonzero values o f i p j i i k are non-Gaussian and are fluxes); hence, (6.8) is presumably correct only in a homogeneous situation and should carry a correction for inhomogeneity; it would be possible to justify neglecting fourth cumulants (recall that all cumulants of second order or higher vanish for a Gaussian) in an equation for third cumulants if cumulants of successive orders had relaxation times which were successively shorter; if inhomogeneity is responsible for nonzero third moments, a term like U i , , , m is of order &/El relative to (6.7), and if this quantity is small, we could justify neglecting these terms (we have taken Ui,,,to be of order u/l where u and 1 are the scales of the energy containing eddies); finally, how d o we know that (6.7) is the right form? Other linear combinations of third moments are equally attractive. It seems clear that what is needed is a model of turbulence that will relax to Gaussian behavior in the absence of inhomogeneity, successive cumulants relaxing faster, and from the equation for the density (or equivalent quantity) equations for the moments may be obtained. The model must be constructed in such a way that our second moment equations are reproduced unchanged, and the third moment equations are reproduced with as few changes as possible; that is, there are presumably constraints on possible forms, consistent with all moments being derivable from a single density which will relax to a Gaussian distribution in the absence of inhomogeneity. This is a type of realizability although somewhat more subtle. B. A GAUSSIAN MODEL In order to examine these possibilities further, it is simpler to work at first with a passive scalar quantity; the added complexity of a vector quantity such as ui is a nuisance. Consider a quantity 8:

B + eViui + oSiui+ e,iui- 6=

(6.9)

Computational Modeling of Turbulent Flows

163

Since both 8 and ui appear in this equation, we will have to generate an equation equivalent to that for the joint density of 8 and ui. It is convenient to work in terms of the moment generating function, since its derivatives at the origin are the cumulants and it is quadratic for a Gaussian distribution: F = In If the equation for ui iri

G,

G = exp[ik,u, + i ~ ] .

+ ui,juj+ ui,,uj+ ui,juj -a= - p , , / p + pie + vui,jj

(6.10) (6.11)

is multiplied by iki that for 8 is multiplied by il, the two are added, and if the resulting equation is multiplied through by G, averaged and the resulting equation divided by G, we obtain (with a little rearrangement)

F

+ aF/aXj uj + (ik, U,,, + ilO,,) aF/aik, - ik,& aF/ail

+ a2F/axj aik, + (aF/axj)(aF/aikj)- (u,In G),, ~= -ik,p,,G/pG + vik,Gu,,/G + yilG8,jj/G. -

-

--

(6.12)

We have used the fact that G-I

a2G/axjdikj = a2F/axj aik, + (dF/axj)(aF/aikj)

(6.13)

from the definition. Now, the terms on the right-hand side of (6.12) may be rearranged in suggestive forms, and physically appealing assumptions may be made directly regarding the terms. However, one has relatively little feeling for the moment generating function and how it may interact with other quantities. It is more productive to adopt the position that we know what equation we wish to obtain for second order (specifically, the terms on the right of (6.12) must give rise to the rapid terms, the return to isotropy, the pressure transport, and the dissipation) and select for the right-hand side of (6.12) the simplest form that will (a) relax to Gaussian in the absence of inhomogeneity, and (b) will produce the correct second-order terms. Let us consider a joint Gaussian distribution. The moment generating function is given by:

F = - & [ k i k j m + 2 k i l d + l2P].

(6.14)

Let us consider only the return to isotropy terms and the dissipation terms; we will deal with the rapid and pressure transport terms later. Using the simplest forms, we have (with C1 and C, as defined by Zeman and Lumley, 1978) u:u. 1 1 + ... = -C1(E/&(iZpj - 6,,7/3) - 6,,2-?/3, ui +

$+

* a *

=

- Ce(E/z)Q,

*.. = -2$.

(6.15)

John L Lumley

164

The ellipses indicate the omitted terms: production, transport, pressure transport, rapid, etc. We can form the time derivative of F :

p = ... = -i{ki k j [ - Cl(C/z)(ii& - Sij?/3) - 6i,Z/3]

+ 2ki I[ - C , ( E / ~ ) Q +] lz[-2(G/?P)F]},

(6.16)

where the ellipses indicate the other terms omitted on the left side of (6.21)as well as those arising from the rapid and pressure transport terms. Now, in a Gaussian distribution iipq = a2F/aikidikj = F,ij, -

uie = a2F/aikiail = F,il, -

(6.17)

e2 = a2F/aii ail = F , ~ ~ ,

and we propose as the simplest possible generalization of (6.16)

p + ...= - f { k i

k j [ - ~1 ( E / Z M F , i j

- dij ~ , p p / 3 ) S i j ~ . p p W @ I

+ 2ki I[ - C,j(E/z)F,iJ+ 12[ -2(F,j/P)F,11]}-

(6.18)

In the absence of inhomogeneity, buoyancy, and mean gradients, the terms not shown explicitly in (6.18)vanish. The equation then has solutions other than Gaussian, but which all relax to Gaussian. This is easiest to see in the scalar case, setting k = 0. Then

p = 12 ( E- e /7e

)F,ii*

(6.19)

Indicating the nth cumulant by C,, we have (differentiating n times with respect to i l ) (6.20) C, = -n(n - I)(C,j/eT)C,. If we normalize by Cyz, designating P, = C,/C;'2, we have P, = - n(n - ~)(E,/P)P,.

(6.21)

Hence, all normalized cumulants above the second decay to zero, no matter what the initial values, leaving a Gaussian distribution. The choice of right-hand side in (6.18)is not as arbitrary as it seems. To second order, the - lz?P in (6.16)can be replaced only by - 12F,11,ilFv1,or 2F. The second does not converge to Gaussian, producing constant P,; the third diverges, the coefficient in (6.21)being positive. Hence, the only further generalization would be the inclusion of third and higher derivatives. Equation (6.18)appears to be adequate, however, since it produces terms at third order which are also obtained by simple physical reasoning. For example, in the equation for e", the dissipation term may be written as -

e3/3

-~

+ -.. = y e z e , j j = y ( e 2 e , j ) , , - 2ye,je,je 1:

-22z -2(~,j/?P)8"

(6.22)

Computational Modeling of Turbulent Flows

165

(Zeman and Lumley 1976), where we have used the usual high-Reynolds number approximation, and the idea that fluctuations of O2 and of Eg should be well correlated, so that regions of high O2 dissipate most of it locally. From (6.18) we have (6.23)

which is the same as (6.22). The same is true of the dissipation terms produced at third order from the mechanical part of (6.18Fthey are identical with those obtained from physical reasoning [like that leading to (6.2211.

C. ORDEROF MAGNITUDE ANALYSIS Let us consider now the full equation for the nth cumulant of 8. This equation contains no terms arising from pressure transport or rapid terms, or buoyancy, so that we may obtain it directly from (6.12), using (6.18) as a right-hand side. Cumulants involving velocity also appear, and we will need a designation for them; let us say C . for the cumulant involving 8 n times and uj once. That is, C3j = - 38 uj8. We will consider n greater than 2, so that the last term on the left-hand side of (6.12), being linear in the independent variables, will not contribute. We obtain

& Y-

(6.24)

We wish now to carry out an order of magnitude analysis on this equation, in the case of weak inhomogeneity,in an attempt to identify the more important terms. Let us designate by 8 and u the rms values of 8 and ui, respectively. We will write C , = One,,

Cnj= Pue,,

(6.25)

and we will make the assumption that enj= O(en+l).That is, we presume that since both variables share a distribution they will approach Gaussian together. To fix ideas, let us assume that we are in a nearly parallel, weakly inhomogeneous, steady, two-dimensional flow such as a late wake. Then, from the wake scaling relations, if L and 1 are, respectively, the downstream and cross-stream length scales, U1/L= O(u/l) (Tennekes and Lumley, 1972, p. 109). We also have @ , j = 0(8/l).Such terms as CnjSj are dominated by the

John L..hmley

166

cross-stream gradients, and are hence of order Cnj/l.The orders of magnitude become

c , ,uj: ~ endnull= n(n - l)(Eg/P)dnen[uP/GIn(n- I)], nO,jCn-lj:ne,Pu/I = n(n - l)(Eg/P)Pen[uF/coIn(n- I)],

c , , ~en+ , ~pull : = n(n - l)(G/eZ)Pen[(en+l/en)uP/EgIn(n- 111,

n(n - l)(Eg/@)Cn:n(n

-

l)&/e’)Pe, = n(n - l)(~/@)O”en[l]. (6.26)

Let us designate by 4 = up/&.Of course, in a real wake, or in any naturally occurring turbulent flow, q = O(1). However, we must consider an artificial situation, in which q is small compared to unity, a weakly inhomogeneous situation which could be produced in the laboratory. As q goes to zero, some term in (6.26)must balance the last term. Clearly the first two cannot, since they both vanish. The third is a possibility; this would give qen+l/enn(n - 1) = O(1)

(6.27)

but this would give e3 = O(2/q), e4 = O(12/q2),etc. Here we have used the fact that ez = 1 by definition. Hence, as q goes to zero, all normalized cumulants go to infinity. This is clearly contrary to observation since (as we have already commented) homogeneous turbulence is observed to be Gaussian. The only other possibility is that the quadratic terms balance the last term. Thus, for example, qezen- l/en = O(1)

(6.28)

(ignoring for a moment the various factors of n). This gives en = O(qe,- 1) (again using e2 = 1) and e3 = O(q)so that all normalized cumulants vanish as they ought to, and higher orders vanish faster than the lower orders do. All these terms are of the same order, the general term being given by qek en- k + 1 /en*

(6.29)

Computational Modeling of Turbulent Flows

167

If (6.28) holds, en = O(q'-2), (6.29) is also of order 1 . Hence, the entire collection of quadratic terms must be of order unity. Since there are n - 2 of them, each must be approximately of order ( n - 2 ) - :

'

(J

If

qe,e,-,+,/e,n(n - 1 ) = O(n - 2 ) - ' .

(6.30)

e, = O(k!(q/2)k-2),

(6.31 )

+

then the left-hand side of (6.30)is given by 2(n - k l)/n(n- 1). If we sum this from 2 to n - 1, we have 1 - 2/n(n - I ) , which approaches unity for large n. Hence, (6.31)must be approximately true for large k. Note that (6.31)does not imply that ek+ /ek+ 0 as k + 00 for fixed q, but rather that e,, /ek 0 for fixed k as q + 0. In fact, ek/ek- = kq/2; by picking a small enough q, we can make this ratio as small as we like for the first few cumulants; eventually, however, for k > 2/q, the cumulants will begin to increase. This is to say that the distribution can be made to look Gaussian in as large a neighborhood of the origin as we like, but that far enough out in the tails it will never look Gaussian. This is not a problem, however, since our concern was simply to be able to neglect fourth cumulants in the equation for third cumulants, and this is justified if q is small enough. Using (6.31),the third term becomes --+

qen+ 1 /enn(n - 1 ) = (q2/2)(n+ l)/n(n- 1 )

(6.32)

so that this is a second-order term in q. Hence, to the zeroth order we can write

(4)

C,,jC,-

,j

+ ... +

in

C,- I,jClj= - n(n - l)(Eo/p)Cn.

(6.33)

The neglect of the second-order term (6.32)corresponds to neglecting cumulants of order n + 1 in the equation for cumulants of order n. The neglect of the first-order term corresponds to the neglect of the substantial derivative and the production terms. Note that, in the case of cumulants involving velocity, which can be handled in exactly the same way, the rapid velocity gradient terms are of the same form as the production terms, and hence will also be of first order. We will deal separately with the problem of the buoyant terms, which requires special treatment. Thus, to zerdh order we have

__

30: Buj = - 6(%/8')8"

(6.34)

John L Lumley

168

and

6 E K j + 4e3eUj = - 12(G/e’)(F

- 3822).

(6.35)

Equation (6.34) is the standard mixing length approximation; if we took

o2 a - q u , ( i P / ~ J

(6.36)

we would obtain (6.34) on multiplying by 0, but with an unknown coefficient. Equation (6.34) involves the neglect of the fourth cumulant. To first order in (6.34), we should keep the substantial derivative and production terms (presuming that the time derivative is of the same order as the advective term):

8 + q U j + 30,,-

+ 3cBuj= -6(FB/ p”)e .

(6.37)

To second order we should include a zeroth-order approximation for the neglected term C,,,,; this order approximation for C,, will presumably be of the form (6.35), but we must derive it explicitly, since there may be other complications. Deardorff (1978) has pointed out that, if an equation such as (6.37) is carried for the third-order quantities, then diffusion terms (which will be provided by C , , , ) are essential to stabilize the solutions and prevent the development of spurious peaks. It should be noted that when we made our assumption (6.18), we were essentially neglecting derivatives of order 3 and higher, which we now see will be of first order in q. Thus, while (6.34) and (6.35) are consistent approximations, it is possible that (6.37) neglects terms of the order of those retained. TRANSPORT D. THEPRESSURE We must now turn to consideration of cumulants involving velocity, dealing first with the nonbuoyant case. The rapid terms, as already noted, will cause no difficulty, since they are of the same form as the production terms. The only problem is pressure transport. Traditionally, p,ij is split into a deviatoric part and an isotropic part, which is a transport term (though the split is not unambiguous: Lumley, 1975b). There appears to be no reason to k in a similar manner, particularly since it leaves as remainder split m Pi,ik + which one might expect to be poorly correlated (although one has little intuition for such correlations; in fact ui,,is evidently well correlated with p , and perhaps the two have a common large-scale part which is well correlated with uk). On the other hand, the entire pressure gradient correlation in the second-order equation consists of the rapid and

m,

169

Computational Modeling of Turbulent Flows

the return-to-isotropy parts and the pressure transport; if the pressure gradient terms in higher order equations are obtained from an equation such as (6.18), all must be present, their combination representing the entire term; the term obtained by generalizing the return-to-isotropy part alone being only a piece of the total. The pressure term in Eq. (6.12) may be written as - -

+

~-

-ikip,,G/pG = -iki(pG),i/pG i k i z / p c . The second term on the right may be rewritten as - -

- -

-kikjpui3jG/pG- kilpO,iG/pG

-

(6.38) (6.39)

which is of the form of (6.18). In the equation for ui uj , this produces P(uj,i + ui.j)/P while the first term on the right of (6.38) gives

(6.4)

Thus, the separation (6.38) represents the traditional separation into pressure transport and pressure-strain correlation. Since we are carrying to higher orders the generalization of the second term on the right of (6.38) (from Eq. (6.18)), we must also carry to higher orders the generalization of the first term on the right, so that we will have the entire term. Beginning from Eq. (4.3) in a homogeneous flow we may write

- [p'2']/p = ( K i K j / K 2 ) [ U i U j ] .

(6.42)

Multiplying by [uk]* and averaging and integrating, we obtain (6.43) where Sijk is the spectrum of w k . Although it proved to be dangerous with regard to realizability in Section IV, we will attempt here to express the integral as a linear combination of this triple correlation. If we define

we have the following conditions: Iijpqr = Zjipqr; Iijpqr = Iijqpr; Iijpqj =0; IiiPq = the last but one resulting from incompressibility.Realizability only requires that the term vanish if uk vanishes. The most general linear form in w k contains five coefficients. The application of these conditions determines all the coefficients, resulting in the form :

w,, Iijpqr

= ( f ) G i j u p u q u , - (&)(GirJqq&

+

Gj,ii&&).

(6.45)

John L Lamley

170 Finally, we obtain Iijijr

=

, = (314 1 ur/P - P(2

ur2

(6.46)

and since this does vanish if ur vanishes, it satisfies realizability. Just as our viscous terms-could be obtained through third order by generally replacing E by q2c/q2, the form of (6.46) suggests replacing everywhere (6.47) -PIP = (4*- ?)/5 so that the first term on the right-hand side of (6.38) becomes

- ?)I, (6.48) and this term could be added to the right-hand side of the equation for F. (+)iki[(a/axi

+aF/axi)(F,pp + ~

, p ~ , p

This cannot be correct for all orders, however, since for fifth order and above the cubic term in (6.48)introduces terms of lower degree in q than those we have kept, that is, terms of order q-'. This seems unlikely to be true and suggests that there is more wrong with (6.48) than meets the eye. In fact, there is the same ambiguity regarding the correct form as there was relative to the right-hand side of (6.18). We know the form we wish to produce at second order, (6.46), but several possible general expressions will produce this form. In particular, we could eliminate from (6.48) any term which does not contribute to second order, such as the troublesome cubic term; this could be done by eliminating either the second term in the first parenthesis of (6.48), or the second term in the second, or both. The1 first possibility leaves quadratic terms which contribute to third-order cumulants to zeroth order; computation with these forms shows that the contributions from the pressure terms very much over-correct, leading to negative diffusivities. Hence, we must discard this alternative. The second and third possibilities give no zeroth-order contributions to third-order quantities, but the second will produce zeroth-order contributions to fourth-order quantities. This seems unlikely, and lacking better information we select simply (f)iki(a/aXi)(F,pp

- q2).

(6.49)

Note, incidentally, that the inclusion of this expression produces in the equation for the heat flux a term of the form

-( e p ) , k / P = ( 8 q z ) , k / 5

(6.50)

(see Eqs. (3.6) and (3.7)). This term does satisfy realizability as required. E. ZEROTH-ORDER TRANSPORT TERMS We may now proceed to consider the equation for cumulants including velocity. Let us take first C , j , since it appears in the equation for the temper-

Computational Modeling of Turbulent Flows

171

ature skewness. We know from our discussion above that we may ignore substantial derivatives, production and rapid terms, and higher order cumulants, since our reasoning is equally applicable to equations for cumulants involving velocity. Expression (6.49) may be neglected like the third term in (6.26). Using the notation c,jk for the cumulant involving uj, uk, and temperature n times, we obtain for C,,

qG + 3 ( G ) , k G + 330% + 3(G),k82U, =

+ C & / ? ) ( P / G ) ] ( G - 3eZeUj).

-3@,/F)[1

(6.51)

We have neglected C 3 p pIf. this were inserted in Eq. (6.37), we would have the diffusion terms which stability requires. Returning now to third-order quantities, let us consider Differentiating, we obtain to zeroth order

6.

This is close to the form assumed by previous authors, but here we have a definite value for the coefficient. The zeroth-order approximation for Bui uj is (8ui),kvj

=

Finally, for m

+

+

(%j),k%&

(mj),kak

-[c1(@)(G - 428dij/3) + (2E/3z)fiSij _+ 2c&//q2)8ui~j]. k

(6.53)

we obtain

(uij),pulrup + (u,u,),puiup + (uju,)*puiu, =

-3Cl(E/?)[wk

- ($)(dijfi

- (2E/32)(dij&

+

dikfi

+

+ 6ikfi + djkfi)]

djkfi).

(6.54)

Note that if i f j # k # i, then the right-hand side of (6.54) reduces to the form (6.5), but with an explicit coefficient. However, the coefficients are quite different in the different directions. F. MODIFICATIONS In deriving Eqs. (6.15) and the following, we used very simple return-toisotropy forms. If we wish to be more elegant and make use of the material in Section 111, we must consider what is the proper generalization of Eq. (6.18).

172

John L Lumley

In fact, there are many generalizations possible; probably the most convenient is to replace Cl by fi and Ce by &,, leaving the arguments of these functions unchanged; then, when the differentiations are carried out with respect to iki, etc., these quantities act just like C1 and C,, so that the expressions we have derived (6.51)-(6.54) will remain unchanged, except for the simple substitution. The next level of complexity, the substitution of expressions (6.17) in the arguments of /3 and 6fj,would result in the appearance of additional terms in (6.51)-(6.54) corresponding to the differentiation of these coefficients with respect to iki,etc. The additional complexity does not appear to be worth the effort, unless it proves to be necessary. That is, suppose that on computation the coefficients in (6.51)-(6.54)prove to be too large or small; they could be modified by including these terms and adjusted by adjusting the dependency of the coefficients in the Reynolds stress and heat flux equations. We must now consider the effect of buoyancy. If the equation for Cnjis written, a term of the form flj C,, appears on the ri t-hand side. Relative to the return to isotropy term, this term is of order POq /z&n(n - 1). We must make a decision relative to the magnitude of this quantity. One possibility is to say that the influence of buoyancy is weak, that is, the turbulence is primarily shear produced, and the buoyant influence is secondary. PO?/@ is jOlt/u2 where It is a turbulent integral length scale; in supposing this small we are s'aying that the buoyant acceleration PO is small relative to the turbulent acceleration u2/lt. If it is of order q, we can neglect it and the associated rapid terms just like the mean velocity gradient terms. On the other hand, if we are considering a flow that is produced by buoyancy, we must have jOlt/u2 = O(1) (Tennekes and Lumley, 1972, Section 4.6). The buoyancy terms will then be of the same order as the return to isotropy terms, and must be kept to zeroth order, that is, in the expressions (6.51)-(6.54). This, of course, increases the complexity considerably, since the rapid terms in buoyancy must also be kept. These rapid terms may be obtained by the substitution of the expressions (6.17) in the various rapid terms (see Section IV) followed by appropriate differentiation. A simplified version of this gives excellent predictions of the evolution of the buoyancydriven atmospheric surface mixed layer (Zeman and Lumley, 1976, 1978). In the order of magnitude analysis which we carried out on Eq. (6.24), other choices are possible, corresponding to other physical situations. For example, if we had said that &,/eZ = u/l, which is realistic, but that C,j,j = O(Cnj/L),while U i , j= O(u/l)and a,, = O(O/l),which would correspond to a near-equilibrium situation of a turbulence maintained by almost uniform gradients, with weak inhomogeneity, then our small parameter becomes 1/L and the order of magnitude analysis remains the same with the exception of the terms involving the mean velocity and mean temperature

e

Computational Modeling of Turbulent Flows

173

gradients, which become of order one; thus we would have to keep the mean velocity and mean temperature gradient terms, as well as the rapid terms in the mean velocity and mean temperature gradients in the expressions (6.51)-(6.54). Each situation should be considered on its merits. For example, at the inversion base over the atmospheric surface mixed layer, both mean velocity and mean temperature gradients are strong, while the turbulence has been produced primarily by buoyancy. In this situation we may expect that it will be necessary to keep both the buoyancy terms and the terms in the mean gradients in the expressions (6.51)-(6.54). An approximate form of this, keeping the buoyancy terms and the mean temperature gradients, has been worked out in Lumley et al. (1978). Finally, a word of caution is in order. Some situations are not in any sense almost Gaussian. The temperature fluctuations in a heated wake are a case in point. In the intermittent region, the probability density for the temperature fluctuations will have a spike, corresponding to the free-stream temperature, lower than the-local mean (Antonia and Sreenivasan, 1977). This will be combined with an almost Gaussian distribution corresponding to the temperature fluctuations in the turbulent portion of the wake with the average offset to higher temperatures. As one moves toward the axis of the wake, the magnitude of the spike decreases, but never completely disappears, even on the axis, since there is still a small but finite probability of finding oneself in a tongue of nonturbulent fluid engulfed by the wake. This combined density cannot be approximated in any sense by expansions about the Gaussian state. The part of the density corresponding to the fluctuations within the turbulent fluid probably can; the spike, however, must be treated separately. Moments of any order can be written explicitly as moments of the continuous distribution plus terms involving the displacement of the spike from the mean of the continuous distribution and the intermittency (the time spent in the turbulent fluid). A rational approach probably would involve dealing separately with the two parts, using our technique for the continuous part and handling the spike explicitly (Libby, 1975). These considerations apply to any intermittent situation, such as the edge of a wake or jet or boundary layer, the entrainment region at the inversion base in the atmospheric surface rnixed layer, etc. They apply, in addition, to the modeling of chemical reactions. Initially, the reactants are not mixed on a molecular level but may be macroscopically mixed. The density will consist of two spikes, if we imagine reactants such as acid and base that can be identified on a single numerical scale. As the reaction progresses, the spikes will disappear, and a continuous distribution of reactant will appear. Again, we may be able to apply our ideas to the continuous part, but the spikes should be dealt with separately.

174

John L Lumley REFERENCES

J. F. (1971). Turbulent wake in a passively stratified field. ALEXOPOULOS, C. C., and KEFFER, Phys. Fluids 14, 216224. ANMNIA, R. A., and SREENIVASAN, K. R. (1977). Conditional probability densities in a turbulent heated round jet. Australas. Hydraul. Fluid Mech. Con$, 6th, Inst. Eng. Aust. Nat. Conf. Publ. NO. 77/12, pp. 411-414. BATCHELOR, G. K. (1956). “The Theory of Homogeneous Turbulence.” Cambridge Univ. Press, London and New York. BATCHELOR, G. K.,and PROUDMAN, R. 1. (1954). The effect of rapid distortion of a fluid in turbulent motion. Q.J. Mech. Appl. Math. 7, 83-103. CHOU,P.-Y. (1945a). On velocity correlation and the solution of the equation of turbulent fluctuation. Q. Appl. Math. 3, 38-54. CHOU,P.-Y. (1945b). Pressure flow of a turbulent fluid between parallel planes. Q. Appl. Math. 3, 198-209. CHOU,P.-Y. (1947). Turbulent flow along a semi-infinite plane. Q. Appl. Math. 5, 346-353. COLEMAN, B. D., and NOLL,W.(1961). Recent results in the continuum theory of viscoelastic fluids. Ann. N.Y. Acad. Sci. 89,672-714. COMTE-BELLOT, G., and CORRSIN, S. (1966). The use of a contraction to improve the isotropy of grid-generated turbulence. J. Fluid Mech. 25, 657-682. CORRSIN, S. (1951). The decay of isotropic temperature fluctuations in an isotropic turbulence. J . Aeronaut. Sci. 18, 417-423. CROW,S. C. (1968). The viscoelastic properties of he-grained incompressible turbulence. J. Fluid Mech. 33, 1-200. DALY,B. J., and HARLOW, F. H. (1970). Transport equations in turbulence. Phys. Fluids 13, 2634-2649. DAVIDOV, B. I. (1958). Phenomenological equation of statistical dynamics of an incompressible turbulent fluid. Zh. Eksp. Teor. Fiz. 35, 527-529. DAVIDOV, B. I. (1959a). Statistical dynamics of an incompressible turbulent fluid. Dokl. Akad. Nauk SSSR 127, 768-771. DAVIDOV, B. I. (1959b). Statistical theory of turbulence. Dokl. Akad. Nauk SSSR 127,98&982. DAVWV,B. I. (1961). Statistical dynamics of an incompressible turbulent fluid. Dokl. Akad. Nauk SSSR 136,47-50. DEARDORFF, J. W.(1978). Closure of second- and third-moment rate equations for diffusion in homogeneous turbulence. Phys. Fluids 21, 525-530. FRENKIEL, F.N., and KLEBANOFF, P. S. (1967a). Higher order correlations in a turbulent field. Phys. Fluids 10, 507-520. FRENKIEL, F. N., and KLEBANOFF, P. S. (1967b).Correlation measurements in a turbulent flow using high speed computing methods. Phys. Fluids 10, 1737-1747. GENCE,J. N. (1977). Turbulence homogdne associee A des effets de gravite. Thdse, Docteur Ingenieur, Universite Claude Bernard, Lyon. HANJALIC, K.,and LAUNDER, B. E. (1972). A Reynolds stress model of turbulence and its application to thin shear flows. J . Fluid Mech. 52, 609638. KLME,S. J., MORKOVIN, M.J., SOVRAN, G., and COCKRELL, J. J. (4s.) (1969) Proc. Cornput. Turbulent Boundary Layers-1 968 AFOSR-IFP-Stanford. Thermosci. Div., Stanford University, Stanford, California. KOLMOOOROV, A. N. (1941). Local structure of turbulence in an incompressible fluid at very high Reynolds numbers. Dokl. Akad. Nauk SSSR 30,299-303. KOLMOOOROV, A. N. (1942) Equation of turbulent motion of an incompressiblefluid. Izu. Akad. Nauk SSSR, Ser. Fiz. 6, 56-58.

Computational Modeling of Turbulent Flows

175

LAUNDER, B. E. (1975). Progress in the modeling of turbulent transport. In “Lecture Series 76: Prediction Methods for Turbulent Flows.” von Karman Inst. Fluid Dyn., Rhode-St.Genese, Belgium. LAUNDER, B. E., REECE,G. J., and RODI, W. (1975). Progress in the development of a Reynolds stress turbulent closure. J. Fluid Mech. 68, 537-566. LESLIE,D. C. (1973).“Developments in the Theory of Turbulence.” Oxford Univ. Press, London and New York. LEWELLEN, W. S. (1975). “Use of Invariant Modeling,” ARAP Rep. No. 243. Princeton, New Jersey. LIBBY,P. A. (1975). On the prediction of intermittent turbulent flows. J. Fluid Mech. 68, 273-295. LUMLEY, J. L. (1967). The applicability of turbulence research to the solution of internal flow problems. In “Fluid Mechanics of Internal Flow” (G. Sovran, ed.), pp. 152-169. Elsevier, Amsterdam. LUMLEY, J. L. (1970a). “Stochastic Tools in Turbulence.” Academic Press, New York. LUMLEY, J. L. (1970b). Toward a turbulent constitutive relation. J. Fluid Mech. 41, 413-434. LUMLEY, J. L. (1972). Some comments on the energy method. In “Developments in Mechanics” (L. H. N. Lee and A. H. Szewczyk, eds.), Vol. 6, pp. 63-88. Notre Dame University Press, Notre Dame, Indiana. LUMLEY, J. L. (1975a). Introduction. In “Lecture Series 76: Prediction Methods for Turbulent Flows.” von Karman Inst. Fluid Dyn., Rhode-St.-Genese, Belgium. LUMLEY, J. L. (1975b). Pressure strain correlation. Phys. Fluids 18, 750. LUMLEY, J. L., and KHAJEH-NOURI, B. (1974).Computational modeling of turbulent transport. Adv. Geophys. 18, 169. LUMLEY, J. L., and NEWMAN, G. R. (1977). The return to isotropy of homogeneous turbulence. J. Fluid Mech. 82, 161-178. LUMLEY, J. L., and PANOFSKY, H. A. (1964). “The structure of Atmospheric Turbulence” (R. E. Marshak, ed.), Interscience Monographs and Texts in Physics and Astronomy, Vol. 12. Wiley (Interscience), New York. O., and SIES, J. (1978). The influence of buoyancy on turbulent transLUMLEY, J. L., ZEMAN, port. J. Fluid Mech. 84, 581-597. MARBCHAL,J. (1972). Etude experimentale de la deformation plane d’une turbulence homogene. J. Mec. 11, 263-294. MONIN,A. S., and YAGLOM, A. M. (1971).“Statistical Fluid Mechanics” (J. L. Lumley, ed.), Vol. 1. MIT Press, Cambridge, Massachusetts. MONIN,A. S., and YAGLOM, A. M. (1975).“Statistical Fluid Mechanics” (J. L. Lumley, ed.), Vol. 2. MIT Press, Cambridge, Massachusetts. NEWMAN, G. R., and HERRING, J. (1978). Second order modeling and statistical theory modeling of a homogeneous turbulence. J. Fluid Mech. (to be submitted). NEWMAN, G. R., LAUNDER, B., and LUMLEY, J. L. (1978a). Modeling the behavior of homogeneous scalar turbulence. J. Fluid Mech. (to be submitted). NEWMAN, G. R., WARHAFT, Z., and LUMLEY, J. L. (1978b). The decay of heat flux in gridgenerated turbulence. J. Fluid Mech. (to be submitted). ORSZAG,S. A., and PAITERSON, G. S. (1972). Numerical simulation of turbulence. In “Statistical Models and Turbulence,” Lecture Notes in Physics, Vol. 12, pp. 127-147. Springer-Verlag, Berlin and New York. POPE,S. B. (1974). A more general eNective viscosity hypothesis. J . Fluid Mech. 72, 331-340. POPE,S. B. (1977). “An Explanation of the Turbulent Round Jet/Plane Jet Anomaly,” Rep. No. FS/77/12. Imperial College, London. REYNOLDS,W.C. (1976). Computation of turbulent flows. Annu. Reo. Fluid Mech. 8, 183-208.

176

John L Lamley

ROITA,J. C. (1951a). Statistische Theorie Nichthomogener Turbulenz. 1.2. Phys. 129,547-572. ROITA, J. C. (1951b). Statistiche Theorie Nichthomogener Turbulenz. 2. 2.Phys. 132, 51-77. SCHUMANN, U. (1977). Realizability of Reynolds stress turbulence models. Phys. Fluids 20, 721-725. SPALDING,D. B. (1971). Concentration fluctuations in a round turbulent free jet. Chem. Eng. Sci. 26, 95-107. SPENCER, A. J. M., and RIVLIN,R. S. (1959). The theory of matrix polynomials and its application to the mechanics of isotropic continua. Arch. Ration. Mech. Anal. 2, 309-336. SPENCER, A. J. M., and RIVLIN,R. S.(1960). Further results in the theory of matrix polynomials. Arch. Ration. Mech. Anal. 4, 214-230. TENNEKES, H., and LUMLEY,J. L. (1972). “A First Course in Turbulence.” MIT Press, Cambridge, Massachusetts. TOWNSEND, A. A. (1970). Entrainment and the structure of turbulent flow. J. Fluid Mech. 41, 13-46. VAN DYKE,M. (1964). “Perturbation Methods in Fluid Mechanics.” Academic Press, New York. WARHAFT, Z., and LUMLEY, J. L. (1978a).The decay of temperature fluctuations and heat flux in grid-generated turbulence. In “Lecture Notes in Physics,” Vol. 76, pp. 113-123. SpringerVerlag, Berlin and New York. WARHAFT, Z., and LUMLEY, J. L. (1978b). An experimental study of the decay of temperature fluctuations in grid-generated turbulence. J. Fluid Mech. 88,659-684. J. C. (1975). Modeling the planetary boundary layer-extension to the stable case. WYNGAARD, Boundary-Layer Meteorol. 9, 441-460. ZEMAN, 0.. and LUMLEY, J. L. (1976). Modeling buoyancy driven mixed layers. J. A m o s . Sci. 33, 1974-1988. O., and LUMLEY, J. L. (1978).Buoyancy effects in entraining turbulent boundary layers. ZEMAN, A second order closure study. In “Lecture Notes in Physics.” Springer-Verlag, Berlin and New York. To be published.