- Email: [email protected]
Variational method in the statistical theory of turbulence
PHYSICA ELSEVIER
Physica D 105 (1997) 185-202
Variational method in the statistical theory of turbulence C.K. Zoltani a, S. Kovesi-Domokos b,,, G. Domokos b a Advanced Computational and Information Services Directorate, Army Research Laboratory, Aberdeen Proving Ground, MD 21005, USA b The Henry A. Rowland Department of Physics and Astronomy, The Johns Hopkins University, Baltimore, MD 21218, USA
Received 21 December 1995; revised 22 October 1996; accepted 13 December 1996 Communicated by A.M. Albano
Abstract Variational principles are introduced in the statistical theory of turbulence. With the help of those, the various correlation functions are determined by extremizing certain functionals. As an illustration of the method, some sample calculations of mean flow profiles and two-point correlation functions of cylindrically symmetric free jets are carried out, both in a single-phase flow and a non-reactive two-phase one. A Rayleigh-Ritz method is used for the determination of the correlation functions. A judicious choice of the trial functions leads to the determination of the correlation functions with a modest amount of computation. The results are in reasonable agreement with experimental data. Keywords: Turbulence; Variational principles
1. Introduction The study of fully developed turbulence in a fluid flow at high Reynolds numbers is of great practical importance. It is also a challenging problem from the theoretical point of view, due to the very large number of coupled degrees of freedom involved. Hence the description of such fluid flows cannot proceed via usual techniques of approximating the system by a linear one and using perturbation techniques around the linearized limit. Despite the existence of an enormous body of theoretical literature on the subject, there is still no "standard" theoretical treatment of fully developed turbulence. In our previous works we took a statistical approach to this problem [1] following the pioneering works of Martin et al. [2] and DeDominicis and Peliti [3]. In this approach one considers the classical equations of motion of a fluid, viz. the Navier-Stokes equations, perturbed by a random force. From the physical point of view, the random force represents the fluctuations in the fluid; its presence is necessary in order to avoid unstable, non-turbulent solutions of the equations of motion. The statistics of the random force, in turn, generates a statistical distribution of the components of the flow velocity. The correlation functions computed from the latter distribution can be compared directly with the results of measurements. An application of the technique developed * Corresponding author: E-mail: [email protected]. 0167-2789/97/$17.00 Copyright © 1997 Elsevier Science B.V. All rights reserved PH SO 167-2789(96)00009-2
186
C.K. Zoltani et al./Physica D 105 (1997) 185-202
in [ 1] to the case of turbulent channel flow lead to a reasonable agreement with the available experimental data, cf. Burgett [4]. The approach adopted in [4] consisted of a combination of analytical and numerical methods: an analytic approach was followed and numerical computations were used only to the computation of Fourier transforms, integration, etc. While this approach leads to a considerable reduction of computational complexity, it became clear that a straightforward application of those techniques to problems like fluid flow in jets would lead to an unreasonably large amount of computing time. There are basically two features of the approach described in [1-4] which make the computations difficult. First, there was no easy way available in order to satisfy the constraint imposed by the equation of continuity on the various velocity correlation functions; second, the approximation methods used are essentially perturbative in nature. As a result, any attempt at achieving some reasonable accuracy in the computations required a great deal of computational effort. In this work we describe some further developments in the statistical theory of fully developed turbulence which are designed to avoid the difficulties outlined above. First, we describe a technique by means of which the equation of continuity can be automatically satisfied. (We explicitly describe the procedure for the case when the fluid may be regarded an incompressible one, although the method can be generalized to arbitrary, compressible flows.) Second, we develop a variational approach to the computation of correlation functions. This enables us to avoid perturbative approaches altogether. The paper is organized as follows. In Section 2 we briefly review the statistical approach to turbulent flows as described in [1 ]. In Section 3 we introduce vector potentials which enable us to satisfy the equation of continuity identically. We also discuss the method of imposing various symmetry requirements on the correlation functions. We work out explicitly the constraints imposed by the requirement of cylindrical symmetry. The variational approach is described in Section 4. We present the result of a sample calculation in Section 5: some correlation functions of an axisymmetric jet are calculated by means of the technique developed here. Two subsequent sections are then devoted to a generalization of the formalism to the case of non-reactive multiphase flows.
2. Review of the statistical approach Throughout this paper, we frequently use a condensed notation. We consider a class of systems described by a classical equation of motion of the general form:
OtX + F[X] = f ( x ) .
(1)
Here X stands for an element of the vector space o f dynamical variables. F[X] is an autonomous map of the vector space upon itself; x stands for both spatial coordinates and the time variable, x = (x, t). Finally, f ( x ) represents a Gaussian random force acting upon the system; it plays the role of a disordering field, driving the system described by (1) away from a non-chaotic behavior, cf. [ 1]. Specifically, if the system is an incompressible fluid described by the Navier-Stokes equations, X is identified with the velocity field, u. The map F[X] reads: F [ u ] i = (uSOs)U i + Oi p - v V 2 u i
(2)
with the velocity field satisfying the constraint, V • u = O. In that case, the quantity p (the pressure divided by the density) is not an independent dynamical variable: it can be explicitly expressed in terms of u. Due to the constraint upon u, we are free to assume that the perturbing random force is solenoidal. The quantity playing the central role in theories of this type is the generating functional of the correlation functions. It can be expressed in terms of the probability distribution of the random force as follows. Let K stand
C.K. Zoltani et al./Physica D 105 (1997) 185-202
187
for the correlation operator of the random force. Then the generating functional of the correlation functions of X is given by the functional integral: Z[jl =
f
D X e x p ( - [ ( ( O t X - F), K ( O t X - F ) ) + (j, X)]),
(3)
where (., -) stands for a scalar product over the vector space, involving integration over space-time variables and summation over tensor indices, whereas j is an arbitrary function: the functional derivatives of Z with respect to j give the correlation functions. The cumulants are generated by the functional W = - in Z. Letting f to be a white noise and of zero range ("generalized white noise"), viz. K (Xl, x2) = k 8 (tl - t2)83 (Xl - x2) (where k is a constant) is often a satisfactory choice. It is to be emphasized, however, that this choice is by no means intuitively obvious. Rather, one expects the perturbations representing the interaction of the fluid with its surroundings, thermal fluctuations, etc. to be dominated by large wavelengths. This is borne out to some extent by model calculations carried out by Burgett (loc. cit.) for simple geometries. Burgett used a variety of correlations, their form being inspired by the Omstein-Zemike process and obtained a somewhat better agreement with the available data (in particular at short wavelengths) than with a generalized white noise. The functional weight of integration, DX, is proportional to the (infinite) determinant, Det(0t - 8 F / ~ X ) .
(4)
It was shown in [1] that the latter can be expressed in terms of a functional integral over Fadeev-Popov ghosts, [5]. The central questions are: a satisfactory implementation of the constraints imposed upon X (in the case of an incompressible fluid flow, the constraint V. u = 0) as well as the development of suitable approximation techniques for the computation of Z (or its functional derivatives).
3. Incompressible flow: Vector potentials and symmetries An incompressible flow is characterized by the fact that the velocity field is solenoidal. Any solenoidal vector field is obtained as the curl of a vector potential. Hence, we write u = V x A, where A is the vector potential. It is also well known that, given the velocity field, the vector potential is determined only up to the gradient of an arbitrary scalar function: A and A + V H give rise to the the same velocity field. (This is called the gauge freedom in electromagnetic theory: we adopt the same terminology here.) The scalar function H can be chosen so as to simplify the problem of determining the vector potential. The equation satisfied by the vector potential is obtained by substituting the relationship, u = V × A into the Navier-Stokes equation. This is straightforward and we do not reproduce the result here. Correspondingly, in the statistical theory discussed in the preceding section, one first determines the correlation functions of the vector potential; the velocity correlation functions are then obtained by taking the curl with respect to every argument. It is worth remarking that this way of satisfying the incompressibility constraint is not the only possibility: it would be possible to take the constraint into account by means of the Fadeev-Popov method, see [5]. However, the method described here is simpler and leads to results more quickly than any other approach we know of. This is due to the gauge freedom just described. In what follows, we use a gauge such that the component of the vector potential along the mean flow vanishes. ("Axial gauge".) On denoting the mean flow by U, the condition to be satisfied is: U • A = 0. This can always be achieved. In fact, let B be a vector potential which reproduces the velocity field, but its projection onto U is not
188
C.K. Zoltani et al./Physica D 105 (1997) 185-202
necessarily zero. Then, in order to satisfy U • A = 0, with A = B + V H , the scalar H has to satisfy the differential equation U . B + U . V H = 0.
(5)
There remains a residual gauge freedom, viz. a function h satisfying U . V h = 0 can always be added to any solution of 5. Any symmetry requirement should be now imposed upon the correlation functions of the vector potentials instead of the correlation functions of the velocity field itself. A symmetry of a flow means that the correlation functions are invariant under a subgroup of the Euclidean group. For instance, homogeneous turbulence means invariance under translations, hence, the n-point correlation function of the vector potential depends on the n - 1 coordinate differences only, not on the coordinates themselves. Likewise, in the case of isotropic turbulence, the correlation functions of order 2n are proportional to the n-fold direct product of the metric tensor with itself, whereas correlation functions of odd order vanish. (One has to recall that in the case of an isotropic turbulence, the mean flow vanishes.) Let us now concentrate on the two-point correlation function; the construction of the general n-point correlation function proceeds in a similar fashion. For the sake of definiteness, we choose a coordinate system such that its third axis coincides with the direction of the mean flow U. The relationship between the vector potential and the velocity field is a linear one. Hence, a Reynolds decomposition of the vector potential leads to that of the velocity field. Let us write, A:
(.4)+a,
U--V
× (A),
(a):0,
u'--V
×a,
(6)
where u ~ stands for the fluctuating part of the velocity field and (. • -) denotes, as usual, the expectation value of a quantity. In an axial gauge we have A3 = a 3 = 0 ,
(7)
so that only the components of the vector potential transverse to the mean flow are non-vanishing. (These components are denoted by subscripts and superscripts in capital letters, A, B, etc.) The two-point correlation function of the vector potential a (x, t) can be decomposed in a basis of second rank tensors in the plane perpendicular to the mean flow. A basis of such tensors consists of three elements. Correspondingly, the two-point correlation function can be written as follows:
ZAB ~ (aa(Xl, tl)aB(X2,
t2)) :
~ABZI + ~ABfRsxRxSZ2 ~- I (XIAXZB -~ X2AX1B)Z3,
(8)
where ~A8 and GAB are the Kronecker and Levi-Civita tensors, respectively. The functions Zl, Z2 and Z3 are invariant under rotations around the third axis and they are symmetric functions of their arguments, (xi, tl) and (x2, t2). (The antisymmetric product of the two coordinates has been inserted in front of Z2 so as to satisfy the permutation symmetry of the correlation function with all three invariant functions being symmetric.) Thus we found that the two-point correlation tensor of an incompressible fluid has three independent components only. Cylindrical symmetry of the flow further reduces the number of independent elements to 2. In fact, cylindrical symmetry means that the correlation function is invariant under rotations around the third axis; the tensors E and 3 are the only invariant ones under such rotations; hence Z3 = 0. We now give the explicit expressions of the velocity correlation functions for a flow of cylindrical symmetry. GAB : 6AB013023 Z1,
G33
=
GA3
:
--OI3[XIAZ2 q- OzAZ1],
6RS oIROzsZI q- [2 + X~ OIR
+ X~O2R]Z2.
(9)
C.K. Zoltani et al./Physica D 105 (1997) 185-202
189
The notation used in this and subsequent equations is the following. The first subscript denotes the argument in the correlation function; the second subscript or superscript refers to the vector component. The summation convention is used throughout. Cylindrically symmetric turbulent flows of incompressible fluids have been, of course, discussed previously, cf. [6-8]. For such flows, there is some overlap between our results and those obtained in the works just quoted; the results are in agreement with each other. The present approach is, however, more general, since it can be applied to flows of arbitrary symmetry and it is more explicit than the treatment of [8]. It is worth mentioning that the present construction can be generalized in a straightforward fashion to compressible flows where both the velocity field and the density field are dynamical variables. In that case one has to introduce two vector potentials with a correspondingly enlarged group of gauge transformations; this problem will be discussed elsewhere.
4. Variational principles for the correlation functions Variational principles play an important role in every branch of mechanics, cf. [9], in fact, in every branch of physics. From a practical point of view, variational principles are useful in developing non-perturbative approximation techniques. Briefly, the procedure is the following. Given an equation of motion, one finds an action functional from which the equation of motion is obtained as the Lagrangian derivative of the action. The solution of the equation of motion (with the appropriate boundary conditions) then can be approximated by means of the Rayleigh-Ritz method. One devises a form of the solution containing a certain number of arbitrary parameters (the "trial function"); for any value of the parameters, however, the boundary conditions are satisfied. (Finding the appropriate trial function is a matter of intuition; there is no recipe for constructing it.) Next, one inserts the trial function into the expression of the action, which then becomes a function of the parameters. The approximate solution is obtained by extremizing the action as a function of the parameters. Obviously, extremizing the action as a function of a finite number of parameters is a simpler problem than extremizing it on the class of all, sufficiently smooth, functions; the latter being equivalent to solving the equation of motion. The strength of the variational method lies in the fact that, given a class of functions (the trial functions for all physically admissible values of the parameters), the one found by extremizing the action gives an "optimal" solution (typically, in the sense of an £2 norm: the distance between the approximate and exact solutions being minimal). This is, in general, not the case if, given a trial function, the values of the parameters are determined by other means, e.g. using integral invariants. The procedure can be generalized to systems where fluctuations are present. The key to the generalization is the existence of an effective action. (One obvious property of the effective action should be that if the fluctuations are absent, it should coincide with an action described above.) We now proceed to review briefly the effective action formalism and then construct an effective action for the problem at hand. We temporarily return to the condensed notation used in Section 2. Consider the expression of the generating functional of the cumulants, as quoted there. We have
w[j] =
-,n
f
DX exp (t) - [((0t - F), K(Ot - F)) + i(j, X)].
(10)
The averages of the various quantities are given by functional derivatives of W, ~W G(1)=(X(xl))--
3j(xl)
~2W G(1 2 ) = ( X ( x l ) X ( x 2 ) )
6j(xl)gj(x2) +G(1)G(2),
(11)
C.K.Zoltanietal./PhysicaD105(1997)185-202
190
etc. Instead of j, we can now introduce G(1) as a functional argument, by means of a Legendre transformation in function space. Let us define a new functional,
S1 ~- W - (j, G1).
(12)
One readily verifies with the help of (10) that $1 is a functional of G(1), its first functional derivative being given by the expression aSl 6G(1) - - -
j(1).
(13)
The functional Sj can be determined from a functional differential equation obtained by combining Eqs. (10)-(13). This gives
•SI SI-(GI,ujaSI~ =-lnf DXexp(-I{(OtX-F)'K(OtX-F))-\~'X)I ).
(14)
The differential equation can be solved by iteration; one usually takes as a zeroth approximation the functional & o = S[G1] = ((OtG1 - F[G1]), K (OtGl - FIG1])).
(15)
The principal advantage of the functional $1 is that it becomes stationary if the external "source", j, is put equal to zero, see (13). Hence, one may determine $1 up to a certain accuracy, but thereafter a calculation of Gl does not have to rely upon any perturbative approximation scheme. This procedure can be generalized if one is interested in determining the higher-order correlation functions as well as the average of the dynamical variable itself. We outline the procedure for a scheme aimed at determining the functions GI and G2. It is necessary to introduce a bilinear source, h(1, 2), in addition to j, viz. W [ j , h] = - In f DX exp(-[((OtX - F), K (OtX - F ) ) + ( j , X ) + (X, hX)]).
(16)
Similarly, to the procedure outlined above, one now performs a double Legendre transform in order to obtain a functional, $2 [Gl, G2] : $2 = W [ j , h] - ((G1, j ) - Tr(hG2)),
(17)
where the trace is understood in the operator sense. One readily verifies the relations
•$2
- - -
6G2
h,
aS2
- - -
~G1
j.
(18)
Thus, $2 is stationary in both G1 and G2 in the limit of vanishing sources. The equation obeyed by the functional $2 can be determined in the same way as Eq. (14). We merely quote the result
$2-(GI, ~-1//~$2"]- Wr(G2, ~~2j~$2"]
(~$2 X)- {X, ~$G2 ~$$2X~]~. =-lnf DXexp(-[((OtX--F) K (OtX-F))-\~1' J_]J
(19)
Just as for the equation obeyed by $1, the only known method of solving Eq. (19) is by iteration. We give here the first approximation to $2 : $2 ~ = S [ G I ] + T r \ 6 G 1 6 G
1G2
-Tr(lnG2).
(20)
C.K. Zoltani et al./Physica D 105 (1997) 185-202
191
(This form of $21 is obtained by changing to an integration variable, Y = X - G j in Eq. (19); thereafter S [Y + G 1 ] is expanded in powers of Y up to quadratic terms. The functional differential equation is then solved by quadrature.) Functionals of the type (19) were first used in [ 10,11]; the differential equations obeyed by such functionals was obtained by Cornwall et al. [ 12]. In the papers just quoted, it is also shown that higher-order correlation functions can be computed by the application of either one of the variational principles described above. In this work, however, we concentrate upon the average flow (corresponding to G 1) and the two-point correlation function. For this purpose, the use of the variational principle based on the functional $2 is best suited.
5. Calculation of a cylindrically symmetric jet As an illustration of the technique developed above, we now perform a computation of some properties of a free jet, with a cylindrically symmetric mean flow. We compute the functional $2 in the first approximation, cf. (20). The power of the variational principle lies in the fact that one can use a Rayleigh-Ritz method in order to minimize $2. This leads to a considerable reduction in the complexity of computation: instead of solving coupled partial differential equations for the correlation functions, one tries to guess the form of the vector potential for the mean flow and of the functions ZI and Z2 in (9), leaving a few free parameters. On substituting this form into the expression of $2, the problem reduces to the computation of integrals and the minimization of the expression so obtained as a function of the parameters. The success of such an approach depends on finding a reasonable functional form of (A) and of ZI and Z2. This can be accomplished only by a judicious use of one's physical intuition and by trial and erroL
We want to simplify the calculation as much as possible: we guess some simple functional forms which could reproduce the main features of the experimental data. Let us examine now the Navier-Stokes operator entering the expression of the functional $2: N i [ u ] -~ Ot ui - v • 2 u i 4- (uSOs)Lt i 4- Oi p.
(21)
We are interested in fully developed turbulence: that means that it is described by a stationary ensemble. Hence the term containing the time derivative can be omitted from (21). Next we notice that if one introduces dimensionless quantities (as we shall do in what follows), the coefficient of the viscous term, V2u, becomes proportional to 1/Re, where Re stands for the Reynolds number. Hence, for flows of large Reynolds numbers - Re ~ 104, say - it is safe to omit the viscous term too, as long as we are interested in the behavior of the fluid on scales substantially larger than the dissipation scale. Thus, we are going to work with the truncated Navier-Stokes operator, ni[u] 4- Oi p,
(22)
n i [u]
(23)
with =
us
Os ui .
The next task is to eliminate the pressure from the generating functional. To this end, we regard the pressure one of the fluctuating dynamical variables and carry out the functional integration over it explicitly. Using an operator notation as before, we have to compute a functional integral of the form:
f
1
D p e x p - ~ [ n + V p l K [ n + Vpl.
(24)
C.K. Zoltani et al./Physica D 105 (1997) 185-202
192
This is a standard Gaussian integral. The result of the computation, to be inserted into the expression of the functional $2, is the following:
f
1
Dpexp-~[n
+ VplK[n
+ V p ] ----exp
-
eKTn
,
(25)
with K T = K - (VK)t[VKV + ff2]-IVK,
(26)
where 't' denotes the transpose of an operator. In the last equation we inserted a constant ot2 in order to make the inverse of the operator V K V well defined on large length scales (equivalently, at small wave numbers). Evidently, K T is the transverse part of K, V. K T = 0. We now have to make a physical assumption about the correlation operator, K. We argue that the stirring force should point in the direction of the mean flow: this is both intuitively plausible and it is the simplest way of satisfying the requirement of cylindrical symmetry of the problem. Otherwise, we assume the correlation to be of short range. Hence we take (cf. Section 2):
g i j (X, X t) = k ~i3~j363 (x - x').
(27)
In order to obtain a finite result even in the limit ot --+ 0, we have to assume that k is proportional to 1/or, With this, the transverse part of K becomes
KTij (x, x') = a e-C~lZ-Z'l(~2(x a -- x'A)~i3~j3,
(28)
where a is a constant independent of or. At the end of the calculation one may take o~ -+ 0, although there is some evidence that a finite ~e (just as in an Ornstein-Zernike process) gives slightly better results, cf. Burgett [4, Chapter 6]. We are now ready to compute the functional $2 l, Eq. (20). Clearly, the first term reads f
S[U]
d3x d3x ' n i [U(x)] KTij (x, x')n j [U(x')].
(29)
On expressing the average velocity as the curl of its vector potential, U = V x A and using (26), we obtain in the axial gauge: S = a
1
1
d2x dzl dz2 e-etlzJ-z21[--fASO3AseRMOAORAM] 4- ~[03(eRSORAs)2lzl
[ ' ' ']z2,
(30)
where the second factor in square brackets has the same structure as the first one, but it is evaluated at z2. We now notice that the second term in both square brackets in the last equation is a total derivative. Consequently, upon integration by parts, its contribution becomes proportional to ~ and hence it is small for small values of that parameter. Therefore, that term can be omitted without substantially affecting the results and we shall do so in this paper. We next compute the functional derivatives of (30) in order to generate the second term of $21 in Eq. (20). This is a straightforward procedure and we merely quote the result: 32S
= 2a e-~Jz~-z21{2ot ~(zl
--
Z2)~-NArRBORON62(Xl
--
X2)
~AA(Xl, Zl)~AB(X2, Z2) %-or2 ENA~RMONORAM~PB 6QSOpOQAs~Z(xI -- X2)}.
(31)
C.K. Zoltani et al./Physica D 105 (1997) 185-202
193
In this equation and from now on, x~, etc. stand for the transverse components of the position vector, whereas zl, etc. denote the longitudinal component of the same vector. The symbol ~ (z) stands for the sign function, ~ (z) = z~ Izl. We now notice that the quantity ~ can be scaled out of the expression $21. In fact, upon the rescaling, z --+ (1/~)z, A --+ (1/~1/2)A, and Z A 8 ~ (1/Ot3)ZAB, the expression of $21 is multiplied by an overall factor 1/or 4. This, however, can be absorbed into the constant entering the expression of the force correlation function; at any rate, the location of the stationary point of the functional $21 is independent of the value of or. (We remark, however, that higher-order approximations to $2 do not have this scaling property.) The reader certainly notices that so far our approach has been quite general and the results obtained are applicable to almost any flow geometry. We now specialize to the case of a cylindrically symmetric free jet. We choose as characteristic units of length and velocity the diameter of the injection pipe, d, and the mean velocity at the centerline of the jet exit, Uom, respectively. From now on, as it is customary, we measure all velocities and distances in these units. However, in order to keep the notation simple, this is not indicated explicitly in the formulae. (For instance, a velocity component denoted by U is understood to mean U/UOm, or an axial distance, z, is to be interpreted as z / d . ) A cylindrically symmetric mean flow can be described in the axial gauge by means of a vector potential which has a tangential component only. Specifically, in cylindrical coordinates (x 1 = r cos ~0, x 2 = r sin ~p), we choose a trial function for the vector potential of the form A~o = I U (O, z) R(z)2[1 - e x p ( - r 2 / R(z)2) ]
(32)
with the other two components vanishing. Here U(0, z) stands for the mean velocity along the jet axis; the quantity R ( z ) characterizes the radius of the jet at distance z from the exit. The choice of such a functional form is motivated by its simplicity and by various fits and approximate calculations of free jets; see, e.g. the classic text of Hinze [13]: This vector potential gives rise to an axial and a radial component of the mean velocity of the form: 2
.2
Uz(r, z) = U(0, z)e -r /R(~) , Ur(r,z)--
d l n UZ( )0d, z
21R(z)2(1-e-rZ/R(z)2)+R(z)-~z
1
1
e -r2/R(z)2
.
(33)
Next we impose the requirement of cylindrical symmetry on the two-point correlation function, ZaB cf. Eq. (8): (34)
ZAB = 8ABZ 1 4- eABZ 2.
In order to simplify the computation, we arbitrarily set Z2 = 0. There is no strong physical motivation for this choice; it simplifies the calculations to a considerable extent, albeit at the cost of some loss of accuracy. We assume the following functional form for ZI: Z1 = A e x p [ - f ( z l
-- z2) 2 --
g ( X l a -- x 2 A ) ( x I A -- x 2 A ) ] [ U ( O , z I ) U ( 0 ,
Zz)]d(1 + B M ) e -~M.
(35)
In this equation, the parameters A, B, d, S, f , and g are treated as variational parameters, although later we found that the results are rather insensitive to the choice of d and we ended up using a plausible value, d = ½. The expression M is defined as follows: 1
M=
(
XCX¢ +
kR(zj) 2
xzxa) .
(36)
The radius of the jet was assumed to be a linear function of z, R ( z ) = a 4- bz, with the parameters a and b being also varied. We did not, however, vary the functional form of the mean velocity on the jet axis: although this is
C.K. Zoltani et al./Physica D 105 (1997) 185-202
194
possible in principle, we used a fixed form in order to simplify the calculational task. We used the following simple modification of Spalding's formula (cf. [13]), which is quite accurate for z > 2: U(0, z) --
1.35 1 + 0.03z 3/2"
(37)
The next task is to insert the expressions of the vector potential of the average flow and of the two-point correlation function into the functional $21 and to minimize that expression with respect to the parameters. The operations are elementary, but extremely tedious to carry out by hand, although this is not impossible. They are, however, easily carried out by a symbolic manipulation program. The search for an extremum of the functional was carried out numerically. In order to carry out the search for the extremum, one needs starting values of the parameters. These were obtained by comparing Uz(r, z), Eq. (33), to a few data points from [14]; we also chose, arbitrarily, f = g = 3 = 1 as starting values. (The values of f and g cannot be read off directly from the data.) The search resulted in the following values of the parameters: a = 0.25,
b = 0.076,
A = 0.0013,
B = 2.3,
3 = 1.79,
f ~ 5.6,
g ~ 3.7.
(38)
In the approximation used, the value of the functional $2 changes rather slowly with variations of the parameters f and g. This is due to the fact that in this approximation, the interaction between the fluctuations and the mean flow is taken into account, but the mutual interaction of the fluctuations is neglected. (The functional $21 is almost local: it contains 6 functions and derivatives of 6 functions of finite order.) Hence, the values of f and g cannot be determined efficiently. This is not a serious deficiency, however: most measurements determine correlation functions at coincident arguments, which are rather insensitive to these parameters. The correlation functions of velocity fluctuations are determined by inserting the expression of Z into Eq. (9). We summarize the results by listing the expressions of the axial component of the mean flow and of the correlation functions at coincident arguments: Gll = G22 = 2 f A Uz(O, z)e-XZ{1 + B X 2 q- (1/2 f ) [ d ( z ) 2 + X 2 B d ( z ) ( d ( z ) - (B - 1 ) / B b / R ( z ) ) + X4 B b / R ( d ( z ) + [(1 - 2 B ) / 4 B ] b / R ( z ) ) + x S B / 4 ( b / R ( z ) ) 2 ] } , 2
G33 = 4gA Uz(O, z)e - x {1 + xZ[B - (2B - 1)/4R(z)Zg] + x 4 B / 4 R ( z ) 2 g } , 2
G13 = G23 = A x e - x U(O, z){(l - B) d(z) + xZ[Bd(z) + (1 - 2 B ) b / R ] + X 4 b B / R ( z ) } ,
Uz(r, z) = U (O, z ) e -
X2
.
Here, X 2 = Sr2/R(z) 2 and d(z) = (d/dz)(U(O, z)) 1/2.
6. Generalization to multiphase flows The generalization of the techniques and results described in the previous sections to non-reactive multiphase flows is a straightforward matter. Due to the practical importance of such flows, we outline the procedure and carry out a sample calculation. For the sake of simplicity, we restrict ourselves to a two-phase flow: a fluid carries along a dispersed particulate phase. In what follows, variables characterizing the carrier fluid will be given a subscript f, those characterizing a dispersed particulate phase a subscript p. As before, we work in terms of dimensionless variables, by dividing velocities with some characteristic speed, coordinates by a characteristic size of the system under consideration. Both phases are assumed to obey the equations of hydrodynamics (leading approximation in a Chapman-Enskog expansion). The volume fractions of the fluid and particulate phase are denoted by ef and ~p, respectively. Of course,
C.K. Zoltani et al./Physica D 105 (1997) 185-202
¢f -F- 6p •
195
(39)
1.
We now have the equations of continuity, 0
OtCf -t- OX-----[ (6fU~) ~- O,
(40)
o+
(41)
Otep +
(~pUp)= O.
The equations of motion read:
• 3 i j ~f 3p C • 3t (~fu~) + 0~7 (~fufuf) + Pf 3x i -- ~f~p--pf(Up - U~) + f[, •
0
i
j
6p Op OX i
Ot('p";) q- ~ x j ('p"p/'tp) q- PP
C -- ,f,p--(U~pp -- @) + f;.
(42)
(43)
Here ft and fp stand for the perturbing Gaussian random forces acting on the fluid and particulate phases, respectively. In Eqs. (42) and (43) p stands for the external pressure. The quantity C, in general, is a function of I uf - Up I. However, with the exception of some unusual cases, such as highly viscous carriers, very heavy loading, etc. the velocity difference between both phases is not too big: in that case, C can be, to a good approximation, replaced by its value given by Stokes' law. This results in a considerable simplification of the computations. The reader will notice that a viscous term has been omitted from Eq. (42). Despite the fact that C is proportional to the viscosity of the carrier fluid, this approximation is a permissible one unless one is interested in very small scales (large wave numbers). The coupling term is proportional to a velocity difference (in the Stokesian approximation), whereas the viscous term is proportional to the scalar curvature of a velocity field: on large and moderate length scales the latter is less important than the former. Let us now take a look at Eqs. (39), (40)-(41). One immediately realizes that if the flow of one of the phases is approximated by an incompressible one, the flow of the other phase becomes incompressible too. In a large number of practically important situations, the approximation of an incompressible flow is a rather good one. In what follows, we make that approximation. Hence, the pressure can be integrated out and the velocity fields can be derived as curls of appropriate vector potentials. Finally, the expression of the stirring force (or, more precisely, its correlation operator) has to be discussed. In this paper we are mainly concerned with problems where the mean flow is cylindrically symmetrical, with the mean velocity having a large component along the axis of symmetry and a rather small radial component. In such cases one gets satisfactory results by taking a stirring force of the same symmetry and, in fact, neglecting the radial component of the stirring force altogether. For the sake of definiteness, expressions identical in form to the stirring force used for a single-phase flow. We also take for both phases an identical form of the matrix elements of the correlation operator. The calculation has been carried out to the first approximation to the solution of the Cornwall-Jackiw-Tomboulis functional differential equation. Unfortunately, the result of the calculation is neither transparent, nor revealing. The reader will be spared the sight of the result. Instead, in Section 7, we present the results of a calculation of the correlation functions of a steady two-phase jet, sufficiently far from the plane of injection. (In this way, one can assume that the turbulence is fully developed: transients died away and the correlation functions are stationary.)
196
C.K. Zoltani et al./Physica D 105 (1997) 185-202
uf/Uo 0
0.2
0.4
0.6
0.8
1.0
0 i~lni~tlll,,,,I,,,,I,,, o
theoretical
•
measured
I
5
qmo
o II)
z./d 10
o O0 0 O0 0
15
•o o go o
2O
Fig. 1. Mean flow: streamwise component. Single phase and/or carrier fluid.
7. Correlation functions of a cylindrically symmetric two-phase jet The calculation of a cylindrically symmetric two-phase jet proceeds along lines as described in the previous sections for the single-phase case. All velocities are measured in units of the carrier fluid velocity, U0, on the axis at the plane of injection. The unit of length is the diameter of the injection pipe. All vector potentials and correlation functions are assumed to be time independent. We chose essentially the same trial functions for both phases as in the case of a single-phase flow. (Of course, the optimal values of the parameters are different.) We used the 3 `/component of the mean velocity on the axis as an input. The data were taken from [14,15]. The following forms give a good fit to the data: 1.35
Uf(z, r = 0) -- 1 + 0.035Z 3/2
- 0.35 e x p ( - 0 , lz2),
UP(z, r = 0) = 0.78 ( 11. 1+3 0.007Z 3/2
0,13 exp(_0.03z2) )
(44) (45)
The functional form of the average of the vector potential was assumed to have a Gaussian radial dependence for both phases. In cylindrical coordinates we have
(A~) = ~R(z) exp
R-~) 2
(46)
with R(z) = a + bz. The optimized values of the parameters a and b are: af = 0.25,
bf = 0.076;
ap = 1.38,
bp = 0.013.
(47)
Likewise, we use the same functional form for the quantity Z] as for the case of a single-phase flow for both phases:
C.K. Zoltani et al./Physica D 105 (1997) 185-202
197
Z
r/d 13.o. . . .
,~.o. . . .
,~.o . . . .
0,0
olo. . . .
,o
2.0
I
3.0
x ~ x ~
I
0.2
[
;~
t
UflUo
0.4
L
2 o
t
0.6-
zlo~9
.~ f
0.8
o
~,odlae~ ~ured
I I J
1.0
Z
r/d 3 'o'
" ,o ' ' ' " , o
....
oo' ....
,'o ....
0.0
,'o
.x
0.2
UflUo
~o' ....
o
0.4
0.6 z/o. la 0.8
o
n~un~
1.0
Fig. 2. Mean flow: radial component at various distances from the jet exit. Single phase and/or carrier fluid.
ZI ~---A e x p [ - f ( z l - z2) 2 - g(xT, ] -- XT,2) 2] exp ( - S N ) (1 + B N ) [U3(Zl, 0)U3(z2, 0)] 1/2.
(48)
Here,
N
=
~
z,7 + -R- -(z2) -Y
(Just as in the single-phase case, we set Z2 = 0.) The optimal values of the parameters are: A f'f ~ 0.0019,
A p'p ~ 0.0013,
B f'f ~ 2.3,
B p'p ~ 3.5,
~f'f ~ 1.79,
~P'P ~ 2.63.
(49)
The values of the parameters f, g are practically identical in both phases, fP'P ~ f f ' f ~ 5.6 and gP.P ~ gf'f ~ 3.7. These values have been computed for the particle loading extracted from [ 14,15]. The cross correlation functions, Z1p'f have basically the same shape; the parameters B, 3, f, g are practically identical. However, the overall scale is much smaller (A p'f << Af'f). This is intuitively obvious: with the loading fraction used in [14,15], the coupling
198
C.K. Zoltani et al./Physica D 105 (1997) 185-202
between both phases is a rather weak one. As a consequence, it is rather hard to obtain a reliable estimate of A p'f. Several iterations fluctuate around a value of the cross correlation scale at least an order of magnitude smaller than (the comparable) fluid and particle correlation scales; however, a stable maximum of the generalized entropy is hard to achieve. Once we have these parameter values, one can take the curl of the mean vector potential and the double curl (with respect to both arguments) of Z,~,# in order to obtain the mean velocity and the correlation functions of the fluctuations, respectively. Experimental data on the fluctuation correlations are often taken at coincident arguments; this was the case in [14,15] too.
8. Discussion The results of the calculation are illustrated in the Figs. 1-6. We plot the calculated and measured values of various mean values and fluctuations in both single-phase and two phase-flows. The results listed above are evaluated using the parameter values given in the preceding sections and they are compared with data taken from [ 14,15]. The experimental circumstances in these references were not exactly identical: in particular, the Reynolds numbers differ by approximately 30% at the jet exit. However, at large Reynolds numbers the profile of the mean flow and the correlation functions should be insensitive to the precise value of Re. The data bear out this conclusion. Moreover, one notices that the optimized parameter values for a single-phase flow and for the fluid component of a two-phase flow are identical within the calculational errors. For this reason, the data and theoretical results of single-phase flows and of the fluid component of two-phase flows have not been plotted in separate figures. (We conjecture that this situation is characteristic of light loads and it changes in the case of heavily loaded fluids. However, we have not carried out a detailed study of this problem so far.
UP/U o 0.2
0.4
0.6
0.8
1.0
,I
0 IIllill IIO ©
theoretical measured o 413 o 113 o
z/d 10
o
o co o co
15
o
o
20 Fig. 3. Mean flow: streamwise component at various distances from the jet exit. Particulate phase.
C.K. Zoltani et al./Physica D 105 (1997) 185-202
~
199
.r
z dd
,o
oo
i!o' ' ' '2!o ' ' ' ~!o
0.0 x 0.2.
UP/uo
0.4 ztD= 9.0
0.6
[] o~
o
nredicled m~sumd
~a x
0.8
r z r/d
10J i i 1210J i i 11101 i L 1O.°
1.0
20
3.0
o.oo o.1 0.2 UP/U o 0.3 0.4'
, a
predicted ~asuroa
0.5
~t~ ~ 0.6
F i g . 4. M e a n flow: radial c o m p o n e n t at various distances from the jet exit. Particulate phase.
In the figures we use conventional notation. The correspondence between the notation used in this paper (better suited for theoretical calculations) and the conventional one, as generally used in the literature and also in [14,15] is the following. (All arguments are suppressed, so, for instance, Gll stands for Gll(X, x). This does not lead to any ambiguity.) (U 2 )
(0 2 )
(tO 2 )
(UV)
G33 - - U0m2, G l l - U0m2, G22- U0m2, G13- U0m2,
Uz=U.
Overall, the simple trial functions used in our variational calculation, with the parameters determined by extremizing $2, are in a reasonable agreement with the data. (The agreement with the data is better at large values of z. This is understandable on physical grounds: at larger values of z, the statistics of the flow is closer to the stationary ensemble we have been working with.) In particular, our assumption of cylindrically symmetric correlation functions appears to be well justified. This is expected on physical grounds: to a very good approximation, the mean flow is cylindrically symmetric. Although this does not necessary imply that the correlation functions themselves should possess the symmetry, one expects violations of this symmetry to occur on relatively short time scales.
C.K. Zoltani et al./Physica D 105 (1997) 185-202
200
Z
r/d
3Jo....
~o~....
2ol . . . .
ool ....
11.o '
F
O. 0 0 0
,
,
,21oJ
i
i
,
31o
X
I
0.005 I I
0.010
-
0.015
-
I
[]
%
0.020
z.,c, = 90
[]
ID I
, =
z X
pte,aict~ ~=umd
F I
0.025
Z
rid 3b 0 i ~ i i 2I0 i i ~ i II.oL ~ J bJ.Oh = I =1]0 b I ~ 121.or J I 13] 0 0.000
.
.
.
.
.
.
.
.I
.
~×
I
0.005
P
~oD
o olo
× ×o
~
clx x
0.015 a
~asur~
0.020
0.025
Fig. 5. Fluctuation of the streamwise component at various distances from the jet exit. Single phase and/or carrier fluid.
Likewise, the assumption about "transverse scaling", namely that the radial dependence of both the mean flow and the correlation functions occurs only in the scale invariant combination, X 2 --- 3r2/R(z) 2, appears to be reasonably well satisfied. The correlation functions G 13 (not plotted) are different from the others for all flows considered: they are predicted to be about two orders of magnitude smaller than the measured values. Various modifications of the form of the trial functions have been tried, with a larger number of parameters. In particular, the assumption Z2 -- 0 was relaxed. As usual in a variational calculation, such a procedure leads to some improvement. However, in this case, the improvement is not dramatic (at least as long as one does not introduce drastic modifications into the functional form of the trial functions). Hence, this feature appears to indicate the limits of the approximations made in this paper. In extremizing the first approximation to the functional $2, one neglects the mutual interaction of the various components of the fluctuations. While the mean flow provides sufficient driving to the diagonal components of the correlation tensor, apparently the off-diagonal ones are predominantly driven by the other components of the fluctuations themselves. One will be able to take into account this effect by computing higher-order approximations to $2 and, at the same time, using improved functional forms in the trial functions.
C.K. Zoltani et al./Physica D 105 (1997) 185-202
201
r ~ Z
r/d ~o, . . . ~o . . . . . ,'o' . . . . . oo' o.oo
' ,'o'"';o~"'2o
~- ,:~ .×.x
~ . . . .
0.05 •
2
2
/Uo
0.10 za3=9
•
pcoa~toa
m~u.nd
0.15
0.20
~_~r z
rid
~.o,,,, ~o,,, ,~o
o,o,,,, ,o,
......
0,000 0.005
0.010 2 2/U°
%~
0.o15 rio
= 18
0.020 malLu~
0,025 0.030
Fig. 6. Fluctuation of the streamwise component at various distances form the jet exit. Particulate phase.
In summary, the use of a variational method appears to be a useful technique in computing statistical properties of fully developed turbulent flows. In particular, by means of a judicious choice of the trial functions one is able to avoid any reference to linearized approximations. In addition, the application of these techniques is likely to lead to a considerable reduction of computing efforts needed for such calculations.
Acknowledgements This research has been supported in part by the Army Research Office, Project No. I L l 6 1 1 0 2 A H 4 3 . We thank K.R. Sreenivasan for reading an earlier version of the manuscript and for useful comments.
References [1] [2] [3] [4]
G. Domokos, S. Kovesi-Domokos and C.K. Zoltani, Random systems, turbulence and disordering fields, Physica A 153 (1988) 84. P.C. Martin, E.D. Siggia and H.A. Rose, Statistical dynamics of classical systems, Phys. Rev. A 8 (1973) 423. C. De Dominicis and L. Peliti, Field theory, renormalizationand critical dynamics above Tc, Phys. Rev. B 18 (1978) 353. W.S. Burgett, The calculation of correlation functions in the statistical theory of turbulent flows, Dissertation, Johns Hopkins University, TIPAC Technical Report No. 8914 (1989) (unpublished).
202 [5] [6] [7] [8] [9] [10] [11] [12] [13] [14] [15]
C.K. Zoltani et al./Physica D 105 (1997) 185-202
C. Itzykson and J.B. Zuber, Quantum Field Theory (McGraw-Hill, New York, 1980). G. K. Batchelor, The theory of axisymmetric turbulence, Proc. Roy. Soc. London Ser. A 186 (1946) 480. S. Chandrasekhar, The theory of axisymmetric turbulence, Philos. Trans. Roy. Soc. London Ser. A 242 (1950) 557. G. Trevino, An introduction to the theory of nonhomogenous turbulence, Texas J. Sci. 34 (1982) 35. C. Lanczos, Variational Principles of Mechanics, 4th Ed. (University of Toronto Press, Toronto, Canada, 1970). C. De Dominicis and EC. Martin, Stationary entropy principle and renormalization in normal and superfluid systems, J. Math. Phys. 5 (1964) 14. G. Domokos and P. Suranyi, Spontaneous symmetry breaking in quantum field theory, Soy. J. Nuclear Phys. 2 (1964) 361. J.M. Cornwall, R. Jackiw and E. Tomboulis, Effective action for composite operators, Phys. Rev. D 10 (1974) 2428. J.O. Hinze, Turbulence (McGraw-Hill, New York, 1975) Ch. 6. C.K. Zoltani and A.E Bicen, A laser Doppler anemometry investigation of single phase and particle jet flows, Experiments in Fluids, 9 (1990) 295. C.K. Zoltani and A.E Bicen, BRL Report No. 3100 (1990).
Physica D 105 (1997) 185-202
Variational method in the statistical theory of turbulence C.K. Zoltani a, S. Kovesi-Domokos b,,, G. Domokos b a Advanced Computational and Information Services Directorate, Army Research Laboratory, Aberdeen Proving Ground, MD 21005, USA b The Henry A. Rowland Department of Physics and Astronomy, The Johns Hopkins University, Baltimore, MD 21218, USA
Received 21 December 1995; revised 22 October 1996; accepted 13 December 1996 Communicated by A.M. Albano
Abstract Variational principles are introduced in the statistical theory of turbulence. With the help of those, the various correlation functions are determined by extremizing certain functionals. As an illustration of the method, some sample calculations of mean flow profiles and two-point correlation functions of cylindrically symmetric free jets are carried out, both in a single-phase flow and a non-reactive two-phase one. A Rayleigh-Ritz method is used for the determination of the correlation functions. A judicious choice of the trial functions leads to the determination of the correlation functions with a modest amount of computation. The results are in reasonable agreement with experimental data. Keywords: Turbulence; Variational principles
1. Introduction The study of fully developed turbulence in a fluid flow at high Reynolds numbers is of great practical importance. It is also a challenging problem from the theoretical point of view, due to the very large number of coupled degrees of freedom involved. Hence the description of such fluid flows cannot proceed via usual techniques of approximating the system by a linear one and using perturbation techniques around the linearized limit. Despite the existence of an enormous body of theoretical literature on the subject, there is still no "standard" theoretical treatment of fully developed turbulence. In our previous works we took a statistical approach to this problem [1] following the pioneering works of Martin et al. [2] and DeDominicis and Peliti [3]. In this approach one considers the classical equations of motion of a fluid, viz. the Navier-Stokes equations, perturbed by a random force. From the physical point of view, the random force represents the fluctuations in the fluid; its presence is necessary in order to avoid unstable, non-turbulent solutions of the equations of motion. The statistics of the random force, in turn, generates a statistical distribution of the components of the flow velocity. The correlation functions computed from the latter distribution can be compared directly with the results of measurements. An application of the technique developed * Corresponding author: E-mail: [email protected]. 0167-2789/97/$17.00 Copyright © 1997 Elsevier Science B.V. All rights reserved PH SO 167-2789(96)00009-2
186
C.K. Zoltani et al./Physica D 105 (1997) 185-202
in [ 1] to the case of turbulent channel flow lead to a reasonable agreement with the available experimental data, cf. Burgett [4]. The approach adopted in [4] consisted of a combination of analytical and numerical methods: an analytic approach was followed and numerical computations were used only to the computation of Fourier transforms, integration, etc. While this approach leads to a considerable reduction of computational complexity, it became clear that a straightforward application of those techniques to problems like fluid flow in jets would lead to an unreasonably large amount of computing time. There are basically two features of the approach described in [1-4] which make the computations difficult. First, there was no easy way available in order to satisfy the constraint imposed by the equation of continuity on the various velocity correlation functions; second, the approximation methods used are essentially perturbative in nature. As a result, any attempt at achieving some reasonable accuracy in the computations required a great deal of computational effort. In this work we describe some further developments in the statistical theory of fully developed turbulence which are designed to avoid the difficulties outlined above. First, we describe a technique by means of which the equation of continuity can be automatically satisfied. (We explicitly describe the procedure for the case when the fluid may be regarded an incompressible one, although the method can be generalized to arbitrary, compressible flows.) Second, we develop a variational approach to the computation of correlation functions. This enables us to avoid perturbative approaches altogether. The paper is organized as follows. In Section 2 we briefly review the statistical approach to turbulent flows as described in [1 ]. In Section 3 we introduce vector potentials which enable us to satisfy the equation of continuity identically. We also discuss the method of imposing various symmetry requirements on the correlation functions. We work out explicitly the constraints imposed by the requirement of cylindrical symmetry. The variational approach is described in Section 4. We present the result of a sample calculation in Section 5: some correlation functions of an axisymmetric jet are calculated by means of the technique developed here. Two subsequent sections are then devoted to a generalization of the formalism to the case of non-reactive multiphase flows.
2. Review of the statistical approach Throughout this paper, we frequently use a condensed notation. We consider a class of systems described by a classical equation of motion of the general form:
OtX + F[X] = f ( x ) .
(1)
Here X stands for an element of the vector space o f dynamical variables. F[X] is an autonomous map of the vector space upon itself; x stands for both spatial coordinates and the time variable, x = (x, t). Finally, f ( x ) represents a Gaussian random force acting upon the system; it plays the role of a disordering field, driving the system described by (1) away from a non-chaotic behavior, cf. [ 1]. Specifically, if the system is an incompressible fluid described by the Navier-Stokes equations, X is identified with the velocity field, u. The map F[X] reads: F [ u ] i = (uSOs)U i + Oi p - v V 2 u i
(2)
with the velocity field satisfying the constraint, V • u = O. In that case, the quantity p (the pressure divided by the density) is not an independent dynamical variable: it can be explicitly expressed in terms of u. Due to the constraint upon u, we are free to assume that the perturbing random force is solenoidal. The quantity playing the central role in theories of this type is the generating functional of the correlation functions. It can be expressed in terms of the probability distribution of the random force as follows. Let K stand
C.K. Zoltani et al./Physica D 105 (1997) 185-202
187
for the correlation operator of the random force. Then the generating functional of the correlation functions of X is given by the functional integral: Z[jl =
f
D X e x p ( - [ ( ( O t X - F), K ( O t X - F ) ) + (j, X)]),
(3)
where (., -) stands for a scalar product over the vector space, involving integration over space-time variables and summation over tensor indices, whereas j is an arbitrary function: the functional derivatives of Z with respect to j give the correlation functions. The cumulants are generated by the functional W = - in Z. Letting f to be a white noise and of zero range ("generalized white noise"), viz. K (Xl, x2) = k 8 (tl - t2)83 (Xl - x2) (where k is a constant) is often a satisfactory choice. It is to be emphasized, however, that this choice is by no means intuitively obvious. Rather, one expects the perturbations representing the interaction of the fluid with its surroundings, thermal fluctuations, etc. to be dominated by large wavelengths. This is borne out to some extent by model calculations carried out by Burgett (loc. cit.) for simple geometries. Burgett used a variety of correlations, their form being inspired by the Omstein-Zemike process and obtained a somewhat better agreement with the available data (in particular at short wavelengths) than with a generalized white noise. The functional weight of integration, DX, is proportional to the (infinite) determinant, Det(0t - 8 F / ~ X ) .
(4)
It was shown in [1] that the latter can be expressed in terms of a functional integral over Fadeev-Popov ghosts, [5]. The central questions are: a satisfactory implementation of the constraints imposed upon X (in the case of an incompressible fluid flow, the constraint V. u = 0) as well as the development of suitable approximation techniques for the computation of Z (or its functional derivatives).
3. Incompressible flow: Vector potentials and symmetries An incompressible flow is characterized by the fact that the velocity field is solenoidal. Any solenoidal vector field is obtained as the curl of a vector potential. Hence, we write u = V x A, where A is the vector potential. It is also well known that, given the velocity field, the vector potential is determined only up to the gradient of an arbitrary scalar function: A and A + V H give rise to the the same velocity field. (This is called the gauge freedom in electromagnetic theory: we adopt the same terminology here.) The scalar function H can be chosen so as to simplify the problem of determining the vector potential. The equation satisfied by the vector potential is obtained by substituting the relationship, u = V × A into the Navier-Stokes equation. This is straightforward and we do not reproduce the result here. Correspondingly, in the statistical theory discussed in the preceding section, one first determines the correlation functions of the vector potential; the velocity correlation functions are then obtained by taking the curl with respect to every argument. It is worth remarking that this way of satisfying the incompressibility constraint is not the only possibility: it would be possible to take the constraint into account by means of the Fadeev-Popov method, see [5]. However, the method described here is simpler and leads to results more quickly than any other approach we know of. This is due to the gauge freedom just described. In what follows, we use a gauge such that the component of the vector potential along the mean flow vanishes. ("Axial gauge".) On denoting the mean flow by U, the condition to be satisfied is: U • A = 0. This can always be achieved. In fact, let B be a vector potential which reproduces the velocity field, but its projection onto U is not
188
C.K. Zoltani et al./Physica D 105 (1997) 185-202
necessarily zero. Then, in order to satisfy U • A = 0, with A = B + V H , the scalar H has to satisfy the differential equation U . B + U . V H = 0.
(5)
There remains a residual gauge freedom, viz. a function h satisfying U . V h = 0 can always be added to any solution of 5. Any symmetry requirement should be now imposed upon the correlation functions of the vector potentials instead of the correlation functions of the velocity field itself. A symmetry of a flow means that the correlation functions are invariant under a subgroup of the Euclidean group. For instance, homogeneous turbulence means invariance under translations, hence, the n-point correlation function of the vector potential depends on the n - 1 coordinate differences only, not on the coordinates themselves. Likewise, in the case of isotropic turbulence, the correlation functions of order 2n are proportional to the n-fold direct product of the metric tensor with itself, whereas correlation functions of odd order vanish. (One has to recall that in the case of an isotropic turbulence, the mean flow vanishes.) Let us now concentrate on the two-point correlation function; the construction of the general n-point correlation function proceeds in a similar fashion. For the sake of definiteness, we choose a coordinate system such that its third axis coincides with the direction of the mean flow U. The relationship between the vector potential and the velocity field is a linear one. Hence, a Reynolds decomposition of the vector potential leads to that of the velocity field. Let us write, A:
(.4)+a,
U--V
× (A),
(a):0,
u'--V
×a,
(6)
where u ~ stands for the fluctuating part of the velocity field and (. • -) denotes, as usual, the expectation value of a quantity. In an axial gauge we have A3 = a 3 = 0 ,
(7)
so that only the components of the vector potential transverse to the mean flow are non-vanishing. (These components are denoted by subscripts and superscripts in capital letters, A, B, etc.) The two-point correlation function of the vector potential a (x, t) can be decomposed in a basis of second rank tensors in the plane perpendicular to the mean flow. A basis of such tensors consists of three elements. Correspondingly, the two-point correlation function can be written as follows:
ZAB ~ (aa(Xl, tl)aB(X2,
t2)) :
~ABZI + ~ABfRsxRxSZ2 ~- I (XIAXZB -~ X2AX1B)Z3,
(8)
where ~A8 and GAB are the Kronecker and Levi-Civita tensors, respectively. The functions Zl, Z2 and Z3 are invariant under rotations around the third axis and they are symmetric functions of their arguments, (xi, tl) and (x2, t2). (The antisymmetric product of the two coordinates has been inserted in front of Z2 so as to satisfy the permutation symmetry of the correlation function with all three invariant functions being symmetric.) Thus we found that the two-point correlation tensor of an incompressible fluid has three independent components only. Cylindrical symmetry of the flow further reduces the number of independent elements to 2. In fact, cylindrical symmetry means that the correlation function is invariant under rotations around the third axis; the tensors E and 3 are the only invariant ones under such rotations; hence Z3 = 0. We now give the explicit expressions of the velocity correlation functions for a flow of cylindrical symmetry. GAB : 6AB013023 Z1,
G33
=
GA3
:
--OI3[XIAZ2 q- OzAZ1],
6RS oIROzsZI q- [2 + X~ OIR
+ X~O2R]Z2.
(9)
C.K. Zoltani et al./Physica D 105 (1997) 185-202
189
The notation used in this and subsequent equations is the following. The first subscript denotes the argument in the correlation function; the second subscript or superscript refers to the vector component. The summation convention is used throughout. Cylindrically symmetric turbulent flows of incompressible fluids have been, of course, discussed previously, cf. [6-8]. For such flows, there is some overlap between our results and those obtained in the works just quoted; the results are in agreement with each other. The present approach is, however, more general, since it can be applied to flows of arbitrary symmetry and it is more explicit than the treatment of [8]. It is worth mentioning that the present construction can be generalized in a straightforward fashion to compressible flows where both the velocity field and the density field are dynamical variables. In that case one has to introduce two vector potentials with a correspondingly enlarged group of gauge transformations; this problem will be discussed elsewhere.
4. Variational principles for the correlation functions Variational principles play an important role in every branch of mechanics, cf. [9], in fact, in every branch of physics. From a practical point of view, variational principles are useful in developing non-perturbative approximation techniques. Briefly, the procedure is the following. Given an equation of motion, one finds an action functional from which the equation of motion is obtained as the Lagrangian derivative of the action. The solution of the equation of motion (with the appropriate boundary conditions) then can be approximated by means of the Rayleigh-Ritz method. One devises a form of the solution containing a certain number of arbitrary parameters (the "trial function"); for any value of the parameters, however, the boundary conditions are satisfied. (Finding the appropriate trial function is a matter of intuition; there is no recipe for constructing it.) Next, one inserts the trial function into the expression of the action, which then becomes a function of the parameters. The approximate solution is obtained by extremizing the action as a function of the parameters. Obviously, extremizing the action as a function of a finite number of parameters is a simpler problem than extremizing it on the class of all, sufficiently smooth, functions; the latter being equivalent to solving the equation of motion. The strength of the variational method lies in the fact that, given a class of functions (the trial functions for all physically admissible values of the parameters), the one found by extremizing the action gives an "optimal" solution (typically, in the sense of an £2 norm: the distance between the approximate and exact solutions being minimal). This is, in general, not the case if, given a trial function, the values of the parameters are determined by other means, e.g. using integral invariants. The procedure can be generalized to systems where fluctuations are present. The key to the generalization is the existence of an effective action. (One obvious property of the effective action should be that if the fluctuations are absent, it should coincide with an action described above.) We now proceed to review briefly the effective action formalism and then construct an effective action for the problem at hand. We temporarily return to the condensed notation used in Section 2. Consider the expression of the generating functional of the cumulants, as quoted there. We have
w[j] =
-,n
f
DX exp (t) - [((0t - F), K(Ot - F)) + i(j, X)].
(10)
The averages of the various quantities are given by functional derivatives of W, ~W G(1)=(X(xl))--
3j(xl)
~2W G(1 2 ) = ( X ( x l ) X ( x 2 ) )
6j(xl)gj(x2) +G(1)G(2),
(11)
C.K.Zoltanietal./PhysicaD105(1997)185-202
190
etc. Instead of j, we can now introduce G(1) as a functional argument, by means of a Legendre transformation in function space. Let us define a new functional,
S1 ~- W - (j, G1).
(12)
One readily verifies with the help of (10) that $1 is a functional of G(1), its first functional derivative being given by the expression aSl 6G(1) - - -
j(1).
(13)
The functional Sj can be determined from a functional differential equation obtained by combining Eqs. (10)-(13). This gives
•SI SI-(GI,ujaSI~ =-lnf DXexp(-I{(OtX-F)'K(OtX-F))-\~'X)I ).
(14)
The differential equation can be solved by iteration; one usually takes as a zeroth approximation the functional & o = S[G1] = ((OtG1 - F[G1]), K (OtGl - FIG1])).
(15)
The principal advantage of the functional $1 is that it becomes stationary if the external "source", j, is put equal to zero, see (13). Hence, one may determine $1 up to a certain accuracy, but thereafter a calculation of Gl does not have to rely upon any perturbative approximation scheme. This procedure can be generalized if one is interested in determining the higher-order correlation functions as well as the average of the dynamical variable itself. We outline the procedure for a scheme aimed at determining the functions GI and G2. It is necessary to introduce a bilinear source, h(1, 2), in addition to j, viz. W [ j , h] = - In f DX exp(-[((OtX - F), K (OtX - F ) ) + ( j , X ) + (X, hX)]).
(16)
Similarly, to the procedure outlined above, one now performs a double Legendre transform in order to obtain a functional, $2 [Gl, G2] : $2 = W [ j , h] - ((G1, j ) - Tr(hG2)),
(17)
where the trace is understood in the operator sense. One readily verifies the relations
•$2
- - -
6G2
h,
aS2
- - -
~G1
j.
(18)
Thus, $2 is stationary in both G1 and G2 in the limit of vanishing sources. The equation obeyed by the functional $2 can be determined in the same way as Eq. (14). We merely quote the result
$2-(GI, ~-1//~$2"]- Wr(G2, ~~2j~$2"]
(~$2 X)- {X, ~$G2 ~$$2X~]~. =-lnf DXexp(-[((OtX--F) K (OtX-F))-\~1' J_]J
(19)
Just as for the equation obeyed by $1, the only known method of solving Eq. (19) is by iteration. We give here the first approximation to $2 : $2 ~ = S [ G I ] + T r \ 6 G 1 6 G
1G2
-Tr(lnG2).
(20)
C.K. Zoltani et al./Physica D 105 (1997) 185-202
191
(This form of $21 is obtained by changing to an integration variable, Y = X - G j in Eq. (19); thereafter S [Y + G 1 ] is expanded in powers of Y up to quadratic terms. The functional differential equation is then solved by quadrature.) Functionals of the type (19) were first used in [ 10,11]; the differential equations obeyed by such functionals was obtained by Cornwall et al. [ 12]. In the papers just quoted, it is also shown that higher-order correlation functions can be computed by the application of either one of the variational principles described above. In this work, however, we concentrate upon the average flow (corresponding to G 1) and the two-point correlation function. For this purpose, the use of the variational principle based on the functional $2 is best suited.
5. Calculation of a cylindrically symmetric jet As an illustration of the technique developed above, we now perform a computation of some properties of a free jet, with a cylindrically symmetric mean flow. We compute the functional $2 in the first approximation, cf. (20). The power of the variational principle lies in the fact that one can use a Rayleigh-Ritz method in order to minimize $2. This leads to a considerable reduction in the complexity of computation: instead of solving coupled partial differential equations for the correlation functions, one tries to guess the form of the vector potential for the mean flow and of the functions ZI and Z2 in (9), leaving a few free parameters. On substituting this form into the expression of $2, the problem reduces to the computation of integrals and the minimization of the expression so obtained as a function of the parameters. The success of such an approach depends on finding a reasonable functional form of (A) and of ZI and Z2. This can be accomplished only by a judicious use of one's physical intuition and by trial and erroL
We want to simplify the calculation as much as possible: we guess some simple functional forms which could reproduce the main features of the experimental data. Let us examine now the Navier-Stokes operator entering the expression of the functional $2: N i [ u ] -~ Ot ui - v • 2 u i 4- (uSOs)Lt i 4- Oi p.
(21)
We are interested in fully developed turbulence: that means that it is described by a stationary ensemble. Hence the term containing the time derivative can be omitted from (21). Next we notice that if one introduces dimensionless quantities (as we shall do in what follows), the coefficient of the viscous term, V2u, becomes proportional to 1/Re, where Re stands for the Reynolds number. Hence, for flows of large Reynolds numbers - Re ~ 104, say - it is safe to omit the viscous term too, as long as we are interested in the behavior of the fluid on scales substantially larger than the dissipation scale. Thus, we are going to work with the truncated Navier-Stokes operator, ni[u] 4- Oi p,
(22)
n i [u]
(23)
with =
us
Os ui .
The next task is to eliminate the pressure from the generating functional. To this end, we regard the pressure one of the fluctuating dynamical variables and carry out the functional integration over it explicitly. Using an operator notation as before, we have to compute a functional integral of the form:
f
1
D p e x p - ~ [ n + V p l K [ n + Vpl.
(24)
C.K. Zoltani et al./Physica D 105 (1997) 185-202
192
This is a standard Gaussian integral. The result of the computation, to be inserted into the expression of the functional $2, is the following:
f
1
Dpexp-~[n
+ VplK[n
+ V p ] ----exp
-
eKTn
,
(25)
with K T = K - (VK)t[VKV + ff2]-IVK,
(26)
where 't' denotes the transpose of an operator. In the last equation we inserted a constant ot2 in order to make the inverse of the operator V K V well defined on large length scales (equivalently, at small wave numbers). Evidently, K T is the transverse part of K, V. K T = 0. We now have to make a physical assumption about the correlation operator, K. We argue that the stirring force should point in the direction of the mean flow: this is both intuitively plausible and it is the simplest way of satisfying the requirement of cylindrical symmetry of the problem. Otherwise, we assume the correlation to be of short range. Hence we take (cf. Section 2):
g i j (X, X t) = k ~i3~j363 (x - x').
(27)
In order to obtain a finite result even in the limit ot --+ 0, we have to assume that k is proportional to 1/or, With this, the transverse part of K becomes
KTij (x, x') = a e-C~lZ-Z'l(~2(x a -- x'A)~i3~j3,
(28)
where a is a constant independent of or. At the end of the calculation one may take o~ -+ 0, although there is some evidence that a finite ~e (just as in an Ornstein-Zernike process) gives slightly better results, cf. Burgett [4, Chapter 6]. We are now ready to compute the functional $2 l, Eq. (20). Clearly, the first term reads f
S[U]
d3x d3x ' n i [U(x)] KTij (x, x')n j [U(x')].
(29)
On expressing the average velocity as the curl of its vector potential, U = V x A and using (26), we obtain in the axial gauge: S = a
1
1
d2x dzl dz2 e-etlzJ-z21[--fASO3AseRMOAORAM] 4- ~[03(eRSORAs)2lzl
[ ' ' ']z2,
(30)
where the second factor in square brackets has the same structure as the first one, but it is evaluated at z2. We now notice that the second term in both square brackets in the last equation is a total derivative. Consequently, upon integration by parts, its contribution becomes proportional to ~ and hence it is small for small values of that parameter. Therefore, that term can be omitted without substantially affecting the results and we shall do so in this paper. We next compute the functional derivatives of (30) in order to generate the second term of $21 in Eq. (20). This is a straightforward procedure and we merely quote the result: 32S
= 2a e-~Jz~-z21{2ot ~(zl
--
Z2)~-NArRBORON62(Xl
--
X2)
~AA(Xl, Zl)~AB(X2, Z2) %-or2 ENA~RMONORAM~PB 6QSOpOQAs~Z(xI -- X2)}.
(31)
C.K. Zoltani et al./Physica D 105 (1997) 185-202
193
In this equation and from now on, x~, etc. stand for the transverse components of the position vector, whereas zl, etc. denote the longitudinal component of the same vector. The symbol ~ (z) stands for the sign function, ~ (z) = z~ Izl. We now notice that the quantity ~ can be scaled out of the expression $21. In fact, upon the rescaling, z --+ (1/~)z, A --+ (1/~1/2)A, and Z A 8 ~ (1/Ot3)ZAB, the expression of $21 is multiplied by an overall factor 1/or 4. This, however, can be absorbed into the constant entering the expression of the force correlation function; at any rate, the location of the stationary point of the functional $21 is independent of the value of or. (We remark, however, that higher-order approximations to $2 do not have this scaling property.) The reader certainly notices that so far our approach has been quite general and the results obtained are applicable to almost any flow geometry. We now specialize to the case of a cylindrically symmetric free jet. We choose as characteristic units of length and velocity the diameter of the injection pipe, d, and the mean velocity at the centerline of the jet exit, Uom, respectively. From now on, as it is customary, we measure all velocities and distances in these units. However, in order to keep the notation simple, this is not indicated explicitly in the formulae. (For instance, a velocity component denoted by U is understood to mean U/UOm, or an axial distance, z, is to be interpreted as z / d . ) A cylindrically symmetric mean flow can be described in the axial gauge by means of a vector potential which has a tangential component only. Specifically, in cylindrical coordinates (x 1 = r cos ~0, x 2 = r sin ~p), we choose a trial function for the vector potential of the form A~o = I U (O, z) R(z)2[1 - e x p ( - r 2 / R(z)2) ]
(32)
with the other two components vanishing. Here U(0, z) stands for the mean velocity along the jet axis; the quantity R ( z ) characterizes the radius of the jet at distance z from the exit. The choice of such a functional form is motivated by its simplicity and by various fits and approximate calculations of free jets; see, e.g. the classic text of Hinze [13]: This vector potential gives rise to an axial and a radial component of the mean velocity of the form: 2
.2
Uz(r, z) = U(0, z)e -r /R(~) , Ur(r,z)--
d l n UZ( )0d, z
21R(z)2(1-e-rZ/R(z)2)+R(z)-~z
1
1
e -r2/R(z)2
.
(33)
Next we impose the requirement of cylindrical symmetry on the two-point correlation function, ZaB cf. Eq. (8): (34)
ZAB = 8ABZ 1 4- eABZ 2.
In order to simplify the computation, we arbitrarily set Z2 = 0. There is no strong physical motivation for this choice; it simplifies the calculations to a considerable extent, albeit at the cost of some loss of accuracy. We assume the following functional form for ZI: Z1 = A e x p [ - f ( z l
-- z2) 2 --
g ( X l a -- x 2 A ) ( x I A -- x 2 A ) ] [ U ( O , z I ) U ( 0 ,
Zz)]d(1 + B M ) e -~M.
(35)
In this equation, the parameters A, B, d, S, f , and g are treated as variational parameters, although later we found that the results are rather insensitive to the choice of d and we ended up using a plausible value, d = ½. The expression M is defined as follows: 1
M=
(
XCX¢ +
kR(zj) 2
xzxa) .
(36)
The radius of the jet was assumed to be a linear function of z, R ( z ) = a 4- bz, with the parameters a and b being also varied. We did not, however, vary the functional form of the mean velocity on the jet axis: although this is
C.K. Zoltani et al./Physica D 105 (1997) 185-202
194
possible in principle, we used a fixed form in order to simplify the calculational task. We used the following simple modification of Spalding's formula (cf. [13]), which is quite accurate for z > 2: U(0, z) --
1.35 1 + 0.03z 3/2"
(37)
The next task is to insert the expressions of the vector potential of the average flow and of the two-point correlation function into the functional $21 and to minimize that expression with respect to the parameters. The operations are elementary, but extremely tedious to carry out by hand, although this is not impossible. They are, however, easily carried out by a symbolic manipulation program. The search for an extremum of the functional was carried out numerically. In order to carry out the search for the extremum, one needs starting values of the parameters. These were obtained by comparing Uz(r, z), Eq. (33), to a few data points from [14]; we also chose, arbitrarily, f = g = 3 = 1 as starting values. (The values of f and g cannot be read off directly from the data.) The search resulted in the following values of the parameters: a = 0.25,
b = 0.076,
A = 0.0013,
B = 2.3,
3 = 1.79,
f ~ 5.6,
g ~ 3.7.
(38)
In the approximation used, the value of the functional $2 changes rather slowly with variations of the parameters f and g. This is due to the fact that in this approximation, the interaction between the fluctuations and the mean flow is taken into account, but the mutual interaction of the fluctuations is neglected. (The functional $21 is almost local: it contains 6 functions and derivatives of 6 functions of finite order.) Hence, the values of f and g cannot be determined efficiently. This is not a serious deficiency, however: most measurements determine correlation functions at coincident arguments, which are rather insensitive to these parameters. The correlation functions of velocity fluctuations are determined by inserting the expression of Z into Eq. (9). We summarize the results by listing the expressions of the axial component of the mean flow and of the correlation functions at coincident arguments: Gll = G22 = 2 f A Uz(O, z)e-XZ{1 + B X 2 q- (1/2 f ) [ d ( z ) 2 + X 2 B d ( z ) ( d ( z ) - (B - 1 ) / B b / R ( z ) ) + X4 B b / R ( d ( z ) + [(1 - 2 B ) / 4 B ] b / R ( z ) ) + x S B / 4 ( b / R ( z ) ) 2 ] } , 2
G33 = 4gA Uz(O, z)e - x {1 + xZ[B - (2B - 1)/4R(z)Zg] + x 4 B / 4 R ( z ) 2 g } , 2
G13 = G23 = A x e - x U(O, z){(l - B) d(z) + xZ[Bd(z) + (1 - 2 B ) b / R ] + X 4 b B / R ( z ) } ,
Uz(r, z) = U (O, z ) e -
X2
.
Here, X 2 = Sr2/R(z) 2 and d(z) = (d/dz)(U(O, z)) 1/2.
6. Generalization to multiphase flows The generalization of the techniques and results described in the previous sections to non-reactive multiphase flows is a straightforward matter. Due to the practical importance of such flows, we outline the procedure and carry out a sample calculation. For the sake of simplicity, we restrict ourselves to a two-phase flow: a fluid carries along a dispersed particulate phase. In what follows, variables characterizing the carrier fluid will be given a subscript f, those characterizing a dispersed particulate phase a subscript p. As before, we work in terms of dimensionless variables, by dividing velocities with some characteristic speed, coordinates by a characteristic size of the system under consideration. Both phases are assumed to obey the equations of hydrodynamics (leading approximation in a Chapman-Enskog expansion). The volume fractions of the fluid and particulate phase are denoted by ef and ~p, respectively. Of course,
C.K. Zoltani et al./Physica D 105 (1997) 185-202
¢f -F- 6p •
195
(39)
1.
We now have the equations of continuity, 0
OtCf -t- OX-----[ (6fU~) ~- O,
(40)
o+
(41)
Otep +
(~pUp)= O.
The equations of motion read:
• 3 i j ~f 3p C • 3t (~fu~) + 0~7 (~fufuf) + Pf 3x i -- ~f~p--pf(Up - U~) + f[, •
0
i
j
6p Op OX i
Ot('p";) q- ~ x j ('p"p/'tp) q- PP
C -- ,f,p--(U~pp -- @) + f;.
(42)
(43)
Here ft and fp stand for the perturbing Gaussian random forces acting on the fluid and particulate phases, respectively. In Eqs. (42) and (43) p stands for the external pressure. The quantity C, in general, is a function of I uf - Up I. However, with the exception of some unusual cases, such as highly viscous carriers, very heavy loading, etc. the velocity difference between both phases is not too big: in that case, C can be, to a good approximation, replaced by its value given by Stokes' law. This results in a considerable simplification of the computations. The reader will notice that a viscous term has been omitted from Eq. (42). Despite the fact that C is proportional to the viscosity of the carrier fluid, this approximation is a permissible one unless one is interested in very small scales (large wave numbers). The coupling term is proportional to a velocity difference (in the Stokesian approximation), whereas the viscous term is proportional to the scalar curvature of a velocity field: on large and moderate length scales the latter is less important than the former. Let us now take a look at Eqs. (39), (40)-(41). One immediately realizes that if the flow of one of the phases is approximated by an incompressible one, the flow of the other phase becomes incompressible too. In a large number of practically important situations, the approximation of an incompressible flow is a rather good one. In what follows, we make that approximation. Hence, the pressure can be integrated out and the velocity fields can be derived as curls of appropriate vector potentials. Finally, the expression of the stirring force (or, more precisely, its correlation operator) has to be discussed. In this paper we are mainly concerned with problems where the mean flow is cylindrically symmetrical, with the mean velocity having a large component along the axis of symmetry and a rather small radial component. In such cases one gets satisfactory results by taking a stirring force of the same symmetry and, in fact, neglecting the radial component of the stirring force altogether. For the sake of definiteness, expressions identical in form to the stirring force used for a single-phase flow. We also take for both phases an identical form of the matrix elements of the correlation operator. The calculation has been carried out to the first approximation to the solution of the Cornwall-Jackiw-Tomboulis functional differential equation. Unfortunately, the result of the calculation is neither transparent, nor revealing. The reader will be spared the sight of the result. Instead, in Section 7, we present the results of a calculation of the correlation functions of a steady two-phase jet, sufficiently far from the plane of injection. (In this way, one can assume that the turbulence is fully developed: transients died away and the correlation functions are stationary.)
196
C.K. Zoltani et al./Physica D 105 (1997) 185-202
uf/Uo 0
0.2
0.4
0.6
0.8
1.0
0 i~lni~tlll,,,,I,,,,I,,, o
theoretical
•
measured
I
5
qmo
o II)
z./d 10
o O0 0 O0 0
15
•o o go o
2O
Fig. 1. Mean flow: streamwise component. Single phase and/or carrier fluid.
7. Correlation functions of a cylindrically symmetric two-phase jet The calculation of a cylindrically symmetric two-phase jet proceeds along lines as described in the previous sections for the single-phase case. All velocities are measured in units of the carrier fluid velocity, U0, on the axis at the plane of injection. The unit of length is the diameter of the injection pipe. All vector potentials and correlation functions are assumed to be time independent. We chose essentially the same trial functions for both phases as in the case of a single-phase flow. (Of course, the optimal values of the parameters are different.) We used the 3 `/component of the mean velocity on the axis as an input. The data were taken from [14,15]. The following forms give a good fit to the data: 1.35
Uf(z, r = 0) -- 1 + 0.035Z 3/2
- 0.35 e x p ( - 0 , lz2),
UP(z, r = 0) = 0.78 ( 11. 1+3 0.007Z 3/2
0,13 exp(_0.03z2) )
(44) (45)
The functional form of the average of the vector potential was assumed to have a Gaussian radial dependence for both phases. In cylindrical coordinates we have
(A~) = ~R(z) exp
R-~) 2
(46)
with R(z) = a + bz. The optimized values of the parameters a and b are: af = 0.25,
bf = 0.076;
ap = 1.38,
bp = 0.013.
(47)
Likewise, we use the same functional form for the quantity Z] as for the case of a single-phase flow for both phases:
C.K. Zoltani et al./Physica D 105 (1997) 185-202
197
Z
r/d 13.o. . . .
,~.o. . . .
,~.o . . . .
0,0
olo. . . .
,o
2.0
I
3.0
x ~ x ~
I
0.2
[
;~
t
UflUo
0.4
L
2 o
t
0.6-
zlo~9
.~ f
0.8
o
~,odlae~ ~ured
I I J
1.0
Z
r/d 3 'o'
" ,o ' ' ' " , o
....
oo' ....
,'o ....
0.0
,'o
.x
0.2
UflUo
~o' ....
o
0.4
0.6 z/o. la 0.8
o
n~un~
1.0
Fig. 2. Mean flow: radial component at various distances from the jet exit. Single phase and/or carrier fluid.
ZI ~---A e x p [ - f ( z l - z2) 2 - g(xT, ] -- XT,2) 2] exp ( - S N ) (1 + B N ) [U3(Zl, 0)U3(z2, 0)] 1/2.
(48)
Here,
N
=
~
z,7 + -R- -(z2) -Y
(Just as in the single-phase case, we set Z2 = 0.) The optimal values of the parameters are: A f'f ~ 0.0019,
A p'p ~ 0.0013,
B f'f ~ 2.3,
B p'p ~ 3.5,
~f'f ~ 1.79,
~P'P ~ 2.63.
(49)
The values of the parameters f, g are practically identical in both phases, fP'P ~ f f ' f ~ 5.6 and gP.P ~ gf'f ~ 3.7. These values have been computed for the particle loading extracted from [ 14,15]. The cross correlation functions, Z1p'f have basically the same shape; the parameters B, 3, f, g are practically identical. However, the overall scale is much smaller (A p'f << Af'f). This is intuitively obvious: with the loading fraction used in [14,15], the coupling
198
C.K. Zoltani et al./Physica D 105 (1997) 185-202
between both phases is a rather weak one. As a consequence, it is rather hard to obtain a reliable estimate of A p'f. Several iterations fluctuate around a value of the cross correlation scale at least an order of magnitude smaller than (the comparable) fluid and particle correlation scales; however, a stable maximum of the generalized entropy is hard to achieve. Once we have these parameter values, one can take the curl of the mean vector potential and the double curl (with respect to both arguments) of Z,~,# in order to obtain the mean velocity and the correlation functions of the fluctuations, respectively. Experimental data on the fluctuation correlations are often taken at coincident arguments; this was the case in [14,15] too.
8. Discussion The results of the calculation are illustrated in the Figs. 1-6. We plot the calculated and measured values of various mean values and fluctuations in both single-phase and two phase-flows. The results listed above are evaluated using the parameter values given in the preceding sections and they are compared with data taken from [ 14,15]. The experimental circumstances in these references were not exactly identical: in particular, the Reynolds numbers differ by approximately 30% at the jet exit. However, at large Reynolds numbers the profile of the mean flow and the correlation functions should be insensitive to the precise value of Re. The data bear out this conclusion. Moreover, one notices that the optimized parameter values for a single-phase flow and for the fluid component of a two-phase flow are identical within the calculational errors. For this reason, the data and theoretical results of single-phase flows and of the fluid component of two-phase flows have not been plotted in separate figures. (We conjecture that this situation is characteristic of light loads and it changes in the case of heavily loaded fluids. However, we have not carried out a detailed study of this problem so far.
UP/U o 0.2
0.4
0.6
0.8
1.0
,I
0 IIllill IIO ©
theoretical measured o 413 o 113 o
z/d 10
o
o co o co
15
o
o
20 Fig. 3. Mean flow: streamwise component at various distances from the jet exit. Particulate phase.
C.K. Zoltani et al./Physica D 105 (1997) 185-202
~
199
.r
z dd
,o
oo
i!o' ' ' '2!o ' ' ' ~!o
0.0 x 0.2.
UP/uo
0.4 ztD= 9.0
0.6
[] o~
o
nredicled m~sumd
~a x
0.8
r z r/d
10J i i 1210J i i 11101 i L 1O.°
1.0
20
3.0
o.oo o.1 0.2 UP/U o 0.3 0.4'
, a
predicted ~asuroa
0.5
~t~ ~ 0.6
F i g . 4. M e a n flow: radial c o m p o n e n t at various distances from the jet exit. Particulate phase.
In the figures we use conventional notation. The correspondence between the notation used in this paper (better suited for theoretical calculations) and the conventional one, as generally used in the literature and also in [14,15] is the following. (All arguments are suppressed, so, for instance, Gll stands for Gll(X, x). This does not lead to any ambiguity.) (U 2 )
(0 2 )
(tO 2 )
(UV)
G33 - - U0m2, G l l - U0m2, G22- U0m2, G13- U0m2,
Uz=U.
Overall, the simple trial functions used in our variational calculation, with the parameters determined by extremizing $2, are in a reasonable agreement with the data. (The agreement with the data is better at large values of z. This is understandable on physical grounds: at larger values of z, the statistics of the flow is closer to the stationary ensemble we have been working with.) In particular, our assumption of cylindrically symmetric correlation functions appears to be well justified. This is expected on physical grounds: to a very good approximation, the mean flow is cylindrically symmetric. Although this does not necessary imply that the correlation functions themselves should possess the symmetry, one expects violations of this symmetry to occur on relatively short time scales.
C.K. Zoltani et al./Physica D 105 (1997) 185-202
200
Z
r/d
3Jo....
~o~....
2ol . . . .
ool ....
11.o '
F
O. 0 0 0
,
,
,21oJ
i
i
,
31o
X
I
0.005 I I
0.010
-
0.015
-
I
[]
%
0.020
z.,c, = 90
[]
ID I
, =
z X
pte,aict~ ~=umd
F I
0.025
Z
rid 3b 0 i ~ i i 2I0 i i ~ i II.oL ~ J bJ.Oh = I =1]0 b I ~ 121.or J I 13] 0 0.000
.
.
.
.
.
.
.
.I
.
~×
I
0.005
P
~oD
o olo
× ×o
~
clx x
0.015 a
~asur~
0.020
0.025
Fig. 5. Fluctuation of the streamwise component at various distances from the jet exit. Single phase and/or carrier fluid.
Likewise, the assumption about "transverse scaling", namely that the radial dependence of both the mean flow and the correlation functions occurs only in the scale invariant combination, X 2 --- 3r2/R(z) 2, appears to be reasonably well satisfied. The correlation functions G 13 (not plotted) are different from the others for all flows considered: they are predicted to be about two orders of magnitude smaller than the measured values. Various modifications of the form of the trial functions have been tried, with a larger number of parameters. In particular, the assumption Z2 -- 0 was relaxed. As usual in a variational calculation, such a procedure leads to some improvement. However, in this case, the improvement is not dramatic (at least as long as one does not introduce drastic modifications into the functional form of the trial functions). Hence, this feature appears to indicate the limits of the approximations made in this paper. In extremizing the first approximation to the functional $2, one neglects the mutual interaction of the various components of the fluctuations. While the mean flow provides sufficient driving to the diagonal components of the correlation tensor, apparently the off-diagonal ones are predominantly driven by the other components of the fluctuations themselves. One will be able to take into account this effect by computing higher-order approximations to $2 and, at the same time, using improved functional forms in the trial functions.
C.K. Zoltani et al./Physica D 105 (1997) 185-202
201
r ~ Z
r/d ~o, . . . ~o . . . . . ,'o' . . . . . oo' o.oo
' ,'o'"';o~"'2o
~- ,:~ .×.x
~ . . . .
0.05 •
2
2
0.10 za3=9
•
pcoa~toa
m~u.nd
0.15
0.20
~_~r z
rid
~.o,,,, ~o,,, ,~o
o,o,,,, ,o,
......
0,000 0.005
0.010 2 2
%~
0.o15 rio
= 18
0.020 malLu~
0,025 0.030
Fig. 6. Fluctuation of the streamwise component at various distances form the jet exit. Particulate phase.
In summary, the use of a variational method appears to be a useful technique in computing statistical properties of fully developed turbulent flows. In particular, by means of a judicious choice of the trial functions one is able to avoid any reference to linearized approximations. In addition, the application of these techniques is likely to lead to a considerable reduction of computing efforts needed for such calculations.
Acknowledgements This research has been supported in part by the Army Research Office, Project No. I L l 6 1 1 0 2 A H 4 3 . We thank K.R. Sreenivasan for reading an earlier version of the manuscript and for useful comments.
References [1] [2] [3] [4]
G. Domokos, S. Kovesi-Domokos and C.K. Zoltani, Random systems, turbulence and disordering fields, Physica A 153 (1988) 84. P.C. Martin, E.D. Siggia and H.A. Rose, Statistical dynamics of classical systems, Phys. Rev. A 8 (1973) 423. C. De Dominicis and L. Peliti, Field theory, renormalizationand critical dynamics above Tc, Phys. Rev. B 18 (1978) 353. W.S. Burgett, The calculation of correlation functions in the statistical theory of turbulent flows, Dissertation, Johns Hopkins University, TIPAC Technical Report No. 8914 (1989) (unpublished).
202 [5] [6] [7] [8] [9] [10] [11] [12] [13] [14] [15]
C.K. Zoltani et al./Physica D 105 (1997) 185-202
C. Itzykson and J.B. Zuber, Quantum Field Theory (McGraw-Hill, New York, 1980). G. K. Batchelor, The theory of axisymmetric turbulence, Proc. Roy. Soc. London Ser. A 186 (1946) 480. S. Chandrasekhar, The theory of axisymmetric turbulence, Philos. Trans. Roy. Soc. London Ser. A 242 (1950) 557. G. Trevino, An introduction to the theory of nonhomogenous turbulence, Texas J. Sci. 34 (1982) 35. C. Lanczos, Variational Principles of Mechanics, 4th Ed. (University of Toronto Press, Toronto, Canada, 1970). C. De Dominicis and EC. Martin, Stationary entropy principle and renormalization in normal and superfluid systems, J. Math. Phys. 5 (1964) 14. G. Domokos and P. Suranyi, Spontaneous symmetry breaking in quantum field theory, Soy. J. Nuclear Phys. 2 (1964) 361. J.M. Cornwall, R. Jackiw and E. Tomboulis, Effective action for composite operators, Phys. Rev. D 10 (1974) 2428. J.O. Hinze, Turbulence (McGraw-Hill, New York, 1975) Ch. 6. C.K. Zoltani and A.E Bicen, A laser Doppler anemometry investigation of single phase and particle jet flows, Experiments in Fluids, 9 (1990) 295. C.K. Zoltani and A.E Bicen, BRL Report No. 3100 (1990).