THEORETICAL AND NUMERICAL STUDY OF THE PROBABILITY DISTRIBUTION OF THE CONCENTRATION IN FREE TURBULENT FLOWS V. A. Sabel'nikov

UDC 532.517.4:536.46

In the current literature there are two conventional approaches to the determination of the probability distribution of the concentration in turbulent flows. In the first approach the functional form of the probability density function is specified a priori, usually on the basis of its first two moments, i.e., the average concentration and the variance, for which semiempirical equations are formulated and solved [1-3] (see also [4]*). The second approach is universal and holds greater promise; here the probability density function is determined from the solution of semiempirical equations. The latter are derived with the application of exact analytical versions of the equation of motion and the diffusion equation [5-7] as well as model (for example, Langevin's) forms of the equations [8-10] (see [ii, 12] for the details). Specific examples of the analysis of turbulent combustion on the basis of the equations for the probability distribution of the concentration are virtually nonexistent in the literature (for model situations see, e.g., [8, 9, 12-14]). This situation is attributable not only to technical difficulties associated with the solution of the equations for the probability distribution function, but also to the inadequate confirmation of those equations for the description of chemically inert jets, for which, unlike flame jets, an abundance of experimental material is available. Clearly, the analysis of such flows and detailed comparisons with the experimental data constitute a first and necessary step toward the widespread application of the concentration distribution functions for the description of turbulent combustion. Moreover, such analyses provide a means for better delineating the limits of applicability of the hypotheses underlying the closure of the probability distribution equations. Bearing all of the foregoing in mind, in this study we use the semiempirical equation obtained in [5] for the probability distribution of the concentration of a passive impurity in order to determine the probability distribution and the intermittency factor in self-similar free flows. We show that in this case, because of the alternating coefficient of the derivative with respect to the self-similar coordinate, the statement of the boundary-value problem has much in common with the statement of the boundary-value problem for elliptic-type equations. We solve the stated boundary-value problem by two different methods: integral and finite-difference. We obtain good correspondence with the experimental data for a plane wake of a slightly heated circular cylinder [15, 16] as well as for heated plane [17] and axisymmetrical [18, 19] jets. The fields of the mean velocity, the mass flow, and the scalar dissipation, which determines the mixing rate of the substance to molecular scales, are assumed to be known in the given theory. i. Fundamental Equations. We consider the diffusion of a passive impurity in free turbulent flows such as jets and wakes. For large Reynolds numbers Re the region of turbulent motion in free flows, according to [20], is separated at every instant from the nonturbulent region by a sharp boundary. Velocity fluctuations are also present in the nonturbulent region due to the action-at-a-distance mechanism operating through the pressure fluctuations, but they have a potential character. Inasmuch as the paths of the fluid particles do not go beyond the boundary of the turbulent region [21] and as the concentration field does not have a mechanism analogous to the action at a distance in the case of the velocity field, it is fully justified to assume that the concentration of the impurity in the nonturbulent region can take only constant values if concentration fluctuations are absent in the entrant cross section [19]. Because the diffusion equation is linear, these values can be set equal to zero and unity. As *V. I. Golovichev has also used this approach in his candidate's dissertation: "Numerical stimulation of nonequilibrium processes in turbulent flows of reacting gases and continuous chemical laser systems," Novosibirsk (1977). Zhukovskii. Translated from Fizika Goreniya i Vzryva, Vol. 18, No. 2, pp. 77-88, March-April, 1982. Original article submitted April 22, 1981.

0010-5082/82/1802-0195507.5.0

9 1982 Plenum Publishing Corporation

195

a typical example we cite the diffusion problem in a turbulent jet propagating in a submerged space. Here the concentration is equal to zero in the submerged space and to unity in the potential (unsteady) core of the jet, while all its intermediate values are observed in the turbulent fluid. In this case the probability density function of the concentration P(c) is written in the form

P(c) =

~o~(c) + ~,8(c - I) +

~P~(c), ~lo+ ~, + ~ =I.

Here Pt(c) is the conditional density function of the concentration in the turbulent fluid; c(0 _< c ~ i), concentration of the passive impurity; Yo, YI, probabilities of the concentration values c = 0 and c = i, respectively; ~, intermittency factor; and 6, Dirac delta function. For Pt(c) we have the semiempirical equation and boundary conditions

O~lP,/Ot + &tP,/Ox~ + Ov=(c)~P,/Ox~ = - ~, 9O~P,/Oc~, vJc) = q,~-~(c " ),

P,(0) = P , ( i ) = 0.

in c [5]

(1.1)

(1.2)

Here t i s t h e t i m e ; Xk, C a r t e s i a n c o o r d i n a t e s (k = 1, 2, 3 ) ; N = D ( ~ c / ~ x a . ~ C / S X a ) , i n s t a n taneous scalar dissipation (D i s t h e m o l e c u l a r d i f f u s i v i t y ) ; Uk, h y d r o d y n a m i c v e l o c i t y ; o2 = <(c )2> , variance of the concentration fluctuations; qk -- <(Uk - ) (c -- )>, mass flux; < >, < >t, unconditional averaging and averaging over the turbulent fluid respectively. Iterated Greek-letter subscripts indicate summation from 1 to 3. -

-

We now enumerate the hypothesis underlying the derivation of Eq. (i.i) in [5]. The first hypothesis is qualitative and does not introduce any additional unknown functions. It postulates the existence of intermittency and is based on the experimental data of [20]. The second hypothesis states that the random variables N and c, which are governed by radically different space scales,are statistically independent in accordance with Kolmogorov's theory of locally homogeneous and isotropic turbulence [22], so that the conditional average of the scalar dissipation in the turbulent fluid can be regarded as independent of the concentration for a given value of the latter. The experiments of [23] are found to be satisfactorily consistent with this hypothesis. Thus, the second hypothesis introduces into Eq. (I.i) an unknown function t that is not expressed in terms of Pt" The third hypothesis replaces the conditional average of the fluctuation velocity Vk(C) for a given concentration by a linear function of the concentration* v~(c)

=

qko-~(c - )

[if the velocity and concentration have a normal distribution function P (Uk, c), the latter relation is exact]. The experimental data of [24] are in good agreement with this hypothesis. Here additional unknown functions are also introduced: the mass-flux vector [the variance a 2 is expressed in terms of Pt and ~ (1.5)]. Thus, along with the unknown characteristics, i.e., the probability distribution Pt and the intermittency factor y, Eq. (i.i) also contains the additional functions , t, and qi (governing the fields of the average velocity, the scalar dissipation, and the mass flux respectively), which must be specified by means of experimental data or other considerations (for example, from semiempirical equations). It must be recognized here that = yt, because N = 0 in the nonturbulent fluid. In a downstream region sufficiently far from the impurity source the probability of the equality c = i in the nonturbulent fluid is negligible, i.e., it can be assumed that YI = 0. Thenfor P(c) we obtain from the above-indicated relation P(c) = (i -

We confine the ensuing

discussion

"08(c) + ~PAc).

(1.3)

to the designated region.

*It is clear that the second and third hypotheses must be modified in the case of chemically active impurities. Experimental data are needed for a more definite conclusion. Also, in this case it appears more promising to use the joint distribution function of the velocity, concentration, and temperature.

196

The unconditional

and conditional averages of any function G(c) are determined according

to the well-known rule =-yGPdc, t ~_ ~ ~GPtdc of (i. 3) :

=

. The following relation is valid by virtue

( t - ~/)G(0) +

~~.

(1.4)

In the special case G = c 2 we obtain from (1.4) a relation between the total variance 2 the variance (7t2 with respect to the turbulent fluid:

o~/<~>~ (~ -

o't =

<(C --

v +

t)2>t

o,1<~>,)/~, = [. (r. - - t)zPtdc,

and

(1.5)

t = ~ cPtdc = /~?. Inasmuch as the relation (1.5) for the variance 2 , entering into Eq. (i.i), includes an integral of the unknown function Pt, Eq. (i.i) is therefore a nonlinear integodifferential equation in partial derivatives. After the integration of Eq. (i.i) with respect to all c, the conditions for the normalization of Pt, the definition of in terms of P, and the relation (1.5) between t and yield an equation for the intermittency factor: a?lat + 9a~lax:: + a.~v=( ,.)lax: =. ~, 9 OP,(O)/ac,

(1.6)

where vi(t) = qi ~ -2( t --) = qi ~ -2 < c> _i(i -- y) 9 On the other hand it follows from Eq. (1.6) that the Pt normalization condition is satisfied. Thus, if we integrate (I.i) with respect to all c and make use of (1.6), we readily infer from the resulting expression that if the normalization condition is satisfied at the initial time, it will also be satisfied at all subsequent times. 2. Equations for the Self-Similar Problem. The Cauchy problem for the inverse parabolic equation (i.i) is ill-posed; the smaller the scale of the perturbations i n the initial conditions the more rapidly they will grow, causing the solution to quickly become totally distorted. We note that one possigle regularization procedure in the steady-state case entails the integration of Eq. (i.i) in the backward direction (i.e., in the direction of decreasing longitudinal coordinate x), beginning with a certain cross section x = x~, in which an "initial" condition for Pt and y can be framed on the basis of some suitable considerations. According to the experimental data of [25], at sufficient distances from the "source," which can be a nozzle or a flow obstacle, the statistical characteristics, i.e., the hydrodynamic fields, are described with high accuracy by self-similar functions. It has been shown experimentally [26] that this conclusion carries over to the concentration distribution function. These considerations imply that a self-similar solution of Eq. (I.i) provides logical "initial" conditions for x~/d >> 1 (d is a characteristic linear dimension of the "source"). The analysis of the self-similar solutions is the objective of the ensuing discussion. For self-similar free flows with axial or planar symmetry (submerged circular and plane jets, the wake of a cylinder), Eqs. (i.i) and (1.6) can be written as follows in the boundary-layer approximation after some rather cumbersome manipulations: A

9

(2.1)

a~/laz + B 9 a~/lan = a ~ ] I O n ~ + c,#,

A ( / ~ , - - i ) d~,

af (o)l

d/~\ =

+

(2.2)

where

/ (~1, z) = tPt (c, x, g), ~ = dt, z = y / L (x), L = d (x/d) ~, 13 = (i + i)/2, Ik(z) = ~ ~lk/d~l, 0

A = A o + AI~I, B = B d l + B d i ~, C = Co + Ct~l,

A o = ~ - F O~h l I-2~"?, A I = - - A o / I 2 , B 2 = - - A 1- d i n O / d z ,

B 1 =

-- Bsls

0 = r/y, ~

(2.3)

al,

a 1 = ~u. O2/h, u = l / z i 9dF/dz, Co : - - B X . (al -t- A o ' d y / d z q: A ( ? . d l s / d z ) ( I s - - y)-l,. C1 = - - 2B~ H- [a~ -I- A t (dI~/dz - - d?/dz)] (I~ - - y)-l~ 197

d is the diameter of the cylinder for a wake, the width of the nozzle for a circular jet, or the width of the nozzle for a plane jet, and the parameter i is equal to zero for plane flows and to unity for axisymmetrical flows. For a wake the function F(z) = z, while for jets it is related to the velocity field as follows:

f = z ~ C u z - v)~ -~, where Uo is the nozzle

<~>

exit velocity

=

U~Cz),

=

U , v ( z ) , U, = Uo(d/x)% .

for jet flows or the freestream velocity

for a wake.

The functions r and h characterize the fields of the average concentration and the scalar dissipation in the turbulent fluid: = R(x)r(z), t = UsRlL-lh(z), R = Us/Uo. In the derivation of Eqs. (2.1) and (2.2) we have used the relation between the mass flux qy in the transverse direction and the averageconcentration : q. = B/ziFrUs R, which follows from the averaged diffusion equation (written in the boundary-laye~ approximation) after a single integration, as well as the relation (1.5) for 02 The solution of (2.1) is sought in the domain 0 ~ z < ~, 0 ~ < ~. It is assumed here that the maximum value of the coordinate n is equal to infinity, because in the self-similar zone i/ t >> i. The solution must satisfy the normalization condition Io = 1 and the condition II = 1 implied by the definition of t in terms of PtThe equivalence of the normalization condition and relation (2.2) is verified as in Sec. i. The principal singularities of the equation for the probability distribution of the concentration in self-similar flows are associated with the properties of the coefficient A of the derivative with respect to the coordinate z in (2.1). It follows from (2.3) that

It will be apparent from what and z = ~ are singular.

A>0,

0~<1,:

A=0,

~=12,

follows

"

A<0,

II
(2.4)

z=O.

that A = 0 also for z = ~.

Therefore,

the points

z = 0

Of considerable significance is the change of sign of the coefficient A for ~ = II(z). Consequently, despite the fact that Eq. (2.1) is parabolic, both directions along the z axis are of equal significance (the theory of equations of this type has been intensely developed in recent years; see, e.g., [27]). To understand some of the properties of Eq. (2.1) we analyze the behavior of the solutions at the singular points. 3.

Solution

on the Symmetry Axis.

At the point

1~ + al~f:

z = 0 Eq.

(2.1) has the form

+ (cl + d,~) fc = 0,

(3. i )

w h e r e a l = n l ( m T c ) - 1 , c l = a l ( 1 -- n l ) , d l = a~m, n = rc/~r m = y c h c / ( B U c ~ c ) , ~ = oZ/R2 = r2(%2 __ y ) y _ l . , t h e s u b s c r i p t c c o r r e s p o n d s t o t h e v a l u e o f t h e f u n c t i o n s a t z = 0; and t h e prime signifies differentiation w i t h r e s p e c t t o ~. Equation (3.1) has been investigated in [5] on t h e a s s u m p t i o n t h a t i n t e r m i t t e n c y can be neglected (i.e., Yc = 1 ) . The p a r a m e t e r n i s i n v e r s e l y p r o p o r t i o n a l to the intensity of the concentration fluctua t i o n s and i n t h e i n v e s t i g a t e d c a s e s a c q u i r e s f a i r l y l a r g e v a l u e s (n ~ 4 - 5 ) . It follows from t h e b a l a n c e e q u a t i o n f o r ~2 t h a t t h e v a l u e o f m i s e q u a l t o t h e r a t i o o f t w i c e t h e scalar dissipation 2 t o t h e n e g a t i v e o f t h e c o n v e c t i o n . One o f t h e b o u n d a r y c o n d i t i o n s i s g i v e n b y e x p r e s s i o n from t h e a s s u m p t i o n o f t h e e x i s t e n c e o f t h e moments o f a l l

( 1 . 2 ) , and t h e s e c o n d i s positive orders. Then

f~(0) = 0, lim~h/c = 0, where k is an arbitrary

positive

number.

From

/r Equation (3.3) is also obtained from applying the conditions Io = i and 11 = i by multiplying (3.1) by ~ and integrating equivalent to the normalization condition Since fc(0) = 0, we infer from (3.3) that Yc < i.

198

)

obtained

(3.2)

(2.2) we obtain

dl(l--Vc).

(3.3)

Eq. (3.1) by integrating with respect to all ~ and [the latter holds identically, as can be verified with respect to all ~]. Thus, relation (3.3) is Io = i and can be used in place of the latter. a nontrivial solution of Eq. (3.1) exists only for

Equation (3.1) includes the unknown value of the intermittency factor Yc. The boundary conditions (3.2), the n o r m a l i z a t i o n condition, and the d e f i n i t i o n of I2 cannot be used to find Tc, rather they only establish a relationship b e t w e e n I2 (or n) and Yc. Thus, the general solution of Eq. (3.1) is expressed in terms of the parabolic cylinder function Dv

[28]: f~ = exp (-- t~/4) [R,Dv (t~) +R~Dv (

-

-

t~)]. t., _~ ],z~l, (3.4)

t 1 = t~.

-

2 ]/~m,

-

v = n ~"(m

-

-

~,~) . ~ j l ,

where RI and R2 are arbitrary constants. Using the asymptotic expansion of the functions Dv for large values of the argument, we readily verify that the second b o u n d a r y condition in (3.2) is satisfied for all RI and R2. The other conditions, fc(0) = 0 and Io = i, as well as the d e f i n i t i o n of I2 in terms of fc e s t a b l i s h three relations between the constants RI, R2, Yc, and 12, i.e., give the relation b e t w e e n Tc and 12. place

The simplest p r o c e d u r e for d e t e r m i n i n g transform of Eq. (3.1):

d---p= Since Io = i, we have quired relation

F(0)

-

= i.

rn ~

-

F

aip -~

Then,

(l--~)~-~-/

this r e l a t i o n

alP -~- d1'

integrating

is the following.

F (p) =

the La-

we o b t a i n

the re-

e-~nlfl ~. o

the equations

exp --~n~p + m p - - ( v T i ) l n

We take

for F(p),

p+--~

(3.5)

dp=l.

0

Below,

Eq.

(3.5) will

be simplified

considerably

for i -- Tc << 1 and n >> i.

From the known properties of the functions D~ we infer that fc > 0 only for v < 0, i.e., for m < Vc (this c o n d i t i o n is indicated in [5] on the assumption that Tc = i). According to the e x p e r i m e n t a l data, in a plane w a k e [15] m = 2.6, in a plane jet [17] m = 1.84, in an a x i s y m m e t r i c a l jet [18, 19] m = 0.75, and in every case 1 -- Vc << i. Consequently, the condition of strict p o s i t i v i t y of fc fails in a w a k e and in a plane jet, apparently because of the error of the assumed linear (in c) a p p r o x i m a t i o n for the c o n d i t i o n a l - a v e r a g e v e l o c i t y v(c) for (c -- )o -i >> i. However, even for m > Tc the solution (3.4) can be regarded as fully s a t i s f a c t o r y if the roots fc are situated outside its principal d o m a i n of variation.* A s y m p t o t i c estimates and n u m e r i c a l calculations (to be discussed below) show that this situation obtains in the most important case 1 -- Tc << i, n >> 1 and for the aboveindicated values of the p a r a m e t e r m, in which case the roots fc lie in the interval (c -) -I > 3 (~ > i + B/n). As n § ~ Eq. (3.5) acquires method to evaluate the integral:

the following

--h

This relation is in good agreement 1 -- Tc 5 I0 -~ (see below).

with

form after

exp

the results

application

of the s t e e p e s t - d e s c e n t

(3.6)

v

of a n u m e r i c a l

solution

of Eq.

(3.1)

for

The b o u n d a r y - v a l u e problem for Eq. (3.1) with a specified value of Tc has b e e n solved n u m e r i c a l l y by the m e t h o d of Numerov [29] under the conditions fc(0) = 0 and (3.3) in the interval 0 ~ ~ ~ 5 w i t h a step A~ = 10 -3 [the order of a p p r o x i m a t i o n of the method is O(Ans)] by ranging w i t h respect to the second moment 12 up to c o i n c i d e n c e of the specified second moment with the value calculated a c c o r d i n g to fc" The n e c e s s a r y integrals were computed by Simpson's rule. The calculations show that the e x p e r i m e n t a l range of the values of n for the above-indicated values of m is attained fo~ 1 -- Tc = 10-2-10-4The functions fc in this case turn out to p r a c t i c a l l y coincide w i t h a G a u s s i a n d i s t r i b u t i o n function.

*The first-order s e m i e m p i r i c a l e q u a t i o n p r o p o s e d in [6] for P(c) to have a bounded solution in the general case for m > i.

can be readily

shown not

199

The calculated dependence of the intensity of the fluctuations g = ~t/t = /12 -- I, the momental skewness At, and the kurtosis E t on the intermittency factor Yc for wake flow after a cylinder (m = 2.6, r c = 0.268 [15]) is shown in Fig. i. It is difficult to compare the analytical and experimental data directly, because the intermittency factor has not been measured with the necessary accuracy in the experiments. For example, the experimental data of [15] m = 2.6, n = 4.75, according to Eq. (3.6), corresponds to 1 -- u = 10-4" The value of the momental skewness A t is in good agreement with the experimental results in [16], where A t = 0.4. There are discrepancies in the values of the kurtosis: E t = 0.I in [16]. This fact can be attributed both to errors of the assumptions underlying the derivation of the equation for the probability distribution P(c) and to the difficulties of measuring the higher-order moments. 4. Solution in the Neighborhood of z = ~, It is clearly impossible to specify f for z + ~. We therefore proceed as follows. We assume that for z + ~ the quantity f can be represented by an asymptotic series, the principal term of which is described by a selfsimilar function. Inasmuch as Io = Ii = i, this term does not depend on z, i.e., f(~, z) + f (~) as z § ~. The final result is consistent with the concepts developed in [25] and with the experimental data of [19], from which it follows that the intensity of the concentration fluctuations in the turbulent fluid vary only very slightly with increasing distance from the jet axis. For f~ we obtain from (2.1) an ordinary integrodifferential equation whose coefficients depend on the asymptotic behavior of the functions r, h, and F as z § ~. Let the principal terms in the expansions of these functions have the form

r = r0(z)exp[ -- pi(z)], h = h0(z)exp[ - p2(z)],

(4. l)

F = Fo(z).

Here pl > 0 and p2 > 0. It is assumed in regard to the functions grow (or decrease) slower than a certain power of z. The principal

term of the expansion

of the intermittency

ro, ho, and Fo that they

factor y~(z)

is sought

in the

form

?~ =

-~z hoz---(

exp

--~p(z)

(4.2)

Here p = 2p~ -- p= > 0, and a > 0 is for the time being an unknown expansions (4.1) and (4.2) into (2.1) and (2.2), we obtain

constant.

Substituting

the

/ ~ Jr +n (bB - - a ) / ~ + [b (2+ - - t ) N + a ( t - - ~ ) l / ~ = 0,

(4.3)

12 (0) = a - - b.

(4.4)

H e r e b=a/Iz; + = l i m p J p > O . An a n a l y s i s of the experimental d a t a o f [ 1 5 ] s h o w s t h a t i n Eq. (2.1) A § 0 as z ~ ~. This restriction on t h e b e h a v i o r o f t h e c o e f f i c i e n t A is also a natural consequence of physical considerations as to the existence of the limit f~. Since A ~ (dp/ dz) -~, the function p(z) grows more rapidly than the first power of z. It is established by a standard procedure that one of the linearly independent solutions o f Eq. ( 4 . 3 ) h a s f o r n + ~ t h e a s y m p t o t i c representation exp (--1/3b~n3), while the other has the asymptotic form ~(~-~-2). Inasmuch as the existence of the moments of all positive orders is postulated, in the case under consideration, by contrast with the solution for z = 0 analyzed in Sec. 3, the condition tim ~ = ~ 0 is nontrivial. We therefore have an eigenvalue problem for the constant a entering into expressions (4.2)-(4.4). By the condition of nonnegativity of f~ the constant a must be equal to the first eigenvalue. The function f= is completely determined, thus departing from the solution for z = 0 discussed in Sec. 3, where fc is determined up to a single parameter Yc' whose value, as will be apparent below, is determined solely from the global solution subject to the symmetry boundary condition at z = 0. It is convenient to solve the boundary-value independent variable ~ = ~b ~/3. Then

problem

for

(4.3) by making

l ~ + +~ (~ - - ~x) l g q- [~ (2+ - - t) q- ~ ( t - - +)] f~ = 0. Here ~

200

= ab -=/3,

and the prime signifies

differentiation

with respect

the change of

(4.5) to ~.

,4t,Et

40

0,5

I0

0

5,

A~ 6t

0,~

0

4 -l,g 'tf-Zc) Fig.

Fig. ness the wake

-0,5

i

1,0 Fig. 2

1

i. Calculated curves of the intensity e, momental skewAt, and kurtosis E t of concentration fluctuations versus intermittency factor Y c in the symmetry plane of the far of a circular cylinder.

Fig. 2. Calculated curves of the first eigenvalue a and the intensity 6, momental skewness At, and kurtosis E t of the concentration fluctuations versus the parameter ~. We determined the first eigenvalue D~ for q0 in the interval 0 ~ 0 with a step equal to 0.I by the ranging method. We initially used Numerov's procedure in the interval 0 < ~ ~< 7 with a step A~ = 10 -3 for (4..5) to solve the Cauchy problem f~(0) = 0, f'(0) = i. The value of ~i was selected so as to satisfy the condition f~ > 0 in the investigated interval. Then from the numerical results we tested the rate of convergence of the function f to zero in the neighborhood of the right endpoint of the interval of integration. In every case the rate turned out to be clearly greater than power-law. After finding the first eigenvalue ~I we determined the function f~ by ranging with respect to the second moment I2, solving the Cauchy problem f~(0) = 0, ~f~(0)/$q = a - b for Eq. (4.3) by Numerov's procedure in the interval 0 ~ n ~< 6 with a step A~ = 10 -3 . Figure 2 shows the results of the calculations of the boundary-value problem for the first eigenvalue a, the fluctuation intensity e, the momental skewness At, and the kurtosis E t as a function of ~. According to the experimental data of [15, 16, 19, 30], c = 0.580.74, A t = 0.8-0.9, and E t =--0.i. Figure 3 shows the results of calculations of the function f~ for ~ = i, 0.5, 0. Also shown in the same graph are the experimental data of [31, 32]. It is seen that the best correspondence with the experimental data is obtained for ~ = 0. The solution corresponding to = 0 has been found previously in [5]. It is expressed in terms of an Airy function [28]: __

f= The v a l u e 5. problem

o f ~z i n t h i s

case

b 1/3

b~,3

1,274 ] / - - - - - ~ Ai (y), Y = coincides

with

[q - - / 2 (oo)-I = ~ - - 9~.

the largest

Formulation and Numerical Solution of the f o r Eq. ( 2 . 1 ) i s f o r m u l a t e d a s f o l l o w s : lim~]=0

l i m ] = f| ~-,~

for k > 0 ,

root

of the

Boundary-Value

(4.6) equation

Problem.

Ai(y)

= 0.

The b o u n d a r y - v a l u e

](0, z)----0,

~ Oz2Z+l 0) = 0 (l = O, t, 2, . . - ) .

(5.1)

Nonrigorous arguments can be advanced in support of the solvability of (5.1). They are based on a separate inspection of Eq. (2-.1) in the intervals 0 ~<-n < I2 and ~ > le. In these intervals Eq. (2.1) is exactly equivalent to the equation describing the velocity distribution in a boundary layer. In the interval 0 ~ ~ < I2 the initial conditions must be prescribed at z = 0, and in the interval ~ > I2 at z = ~o. In each of the indicated intervals the solution for f is a functional of the a priori unknown functions 7, I2 and the value of f for ~ = I2 (say, fo(z) = f[I2(z), z]). These functions are then determined from the normalization condi-

201

12

i,OC~--~ ,~oo

//~,

|I

x [51]

f~

r

+ [s2]

i ,

o

z

2

Fig.

q

0

0~,4

Fig.

3

z 0~8

4

Fig. 3. Calculated and experimental results for the probability density function of the concentration in the turbulent fluid. i) ~ = i; 2) 0.5; 3) 0. Fig. 4. Calculated and experimental results for the intermittency factor 7 and the rms value of the concentration fluctuations in the far wake of a circular cylinder, i) Calculation of y by the integral method; 2) calculation of y by the finitedifference method; 3) calculation of o/o c by the integral method; 4) calculation of o/o c by the finite-difference method. tion I2 = i, the definition of I2, and the smootness of f for ~ = I2. (5.1) (symmetry condition) is required in order to determine Yc"

The last relation

in

We have solved the boundary-value problem (5.1) for Eq. (2.1) by two different methods. In both cases, bearing in mind the results of comparing the analytical and experimental data for f for z = ~ (see Fig. 3), for f~ we used the solution (4.6) corresponding to ~ = 0. The first is an integral method. Its fundamental i d e a r e s t s on the results obtained in Secs. 3 and 4, where solutions were obtained at z = 0 and z = ~. This method is similar to the integral method used to calculate the characteristics of a boundary layer. We specify f in the form f = ~fc + (i -- ~)f~, a = [12(z) -- I2(~)][12(0) -- I2(~)] -l. It is required to determine the second moment I2 and the intermittency factor y. The differential equation for I2 obtained by multiplying (2.1) by q2 and then integrating with respect to all q combines with the equation (2.2) for y to form a system of two first-order ordinary differential equations. The boundary conditions have the form lim I s - I2(oo), lira? ~ 0. Since I2(~)

is known,

from these equations we find I2 and y.

It is more practical to solve the inverse problem, where the experimental data of [1518, 30] and the theoretical results for z § ~ in Sec. 4 are used to specify the second moment I 2 ( z ) , a n d the unknown functions are considered to be the dissipation h in the turbulent fluid and the intermittency factor y. The solution thus obtained agrees satisfactorily with the experimental data on the large number of parameters measured in [15-18], such as the variance o 2 , the intermittency factor y, and all aspects of the balance in the equation for 02 , the momental skewness, and the kurtosis. As an example, Fig. 4 shows the results of the calculations of y and the rms value o of the concentration fluctuations, referred to its value o c at z = 0 (curves 1 and 3 respectively). Also shown in the graph are the experimental data on the ratio o/o c from [15] and on y from

[16]. The second method of integration of the boundary-value problem in Sec. 2 is a relaxation method similar to the one used for the calculation of separation flows in boundary layers [33]. We approximated the first derivatives with respect to z with allowance for the direction of information transmission: for A > 0 by means of left-sided differences and for A < 0 by means of right-sided differences. We approximated the first derivative w i t h respect to n analogously in accordance with the sign of B. We approximated the second derivative with respect to q in

202

the usual way. The difference scheme was first-order with respect to Az, An (Az, A~ are the integration steps). For quasilinearization of the nonlinear system of difference equations we used the solution obtained in the preceding iteration. The resulting linear equations were solved in each layer z = const by a double-sweep (modified Gaussian elimination) procedure with successive transition from one layer to the next in the direction of increasing z. The iterations were terminated w h e n the sum of the moduli of the differences between the values of f at all computing points in the current and preceding iterations became smaller than a predetermined small number. For convergence of the iteration process a relaxation coefficient of 10 -2 was used. We now discuss the results of the calculations by this method for a plane wake. It was assumed in this case that f = 0 for N = qm = 3 and that the boundary condition corresponding to z * ~ was prescribed at z = zm = 1 (some of the calculations were carried out for ~m = 4 and zm = 0.8). The main calculations w e r e carried out for An = qm/90 and Az = 1/40. The functions r and h were specified from the experimental results of [15, 16]. The calculations show that in the interval z > 0.25 the distributions of the intermittency factor and the second moment agree satisfactorily with the experimental data. However, in the interval 0 ~ z ~ 0.25 the distributions of both quantities oscillate. For example, y is 15-.20% greater than unity at several computing points. This shortcoming is not removed by refinement of the computing grid. These oscillations are clearly attributable to errors of the finite-difference approximation of Eq. (2.1) in the domain z ~ 0, ~ ~ 0. Thus, if we introduce the "natural" Concentration scale s = (c -- t)/o t the point c = 0 will correspond to the variable s = --t/o t. This value is s = 4-5 at z = 0, i.e., is very large, considering that the solution f varies approximately as exp (--i/2s 2) (see Sec. 3). Hence it is clear that for z ~ 0 and q ~ 0 the solution varies quite a b r u p t l y . Since y is related to ~f(0)/3n, the values of y are calculated with error, which clearly accounts for the oscillations. The indicated difficulty can be circumvented by disregarding the solution near the symmetry plane. It is necessary to specify f at a certain point far from the symmetry plane. It was remarked in Sec. 3 that f is almost normal at z = 0. The experiments of [16] show that in the wake of a cylinder for z < 0.25, i.e., in the region where the intermittency is inconsequential (y ~ i) the momental skewness and kurtosis are relatively small, i.e., f is also close to normal and is well approximated by the first few terms of the Gram--Charlier series

f : ng o (s) [ l q- At~3! (s 3 - - 3s) + Et/4! (s 4 - - 6s 2 -6 3) + -q- A ~ / 7 2 (s6.-- 15s 4 q- 45s ~ q- 15)],:

go = ( 2 g ) - ~

exp ( - - 1/2s2),~ s = n (~ - -

1).

In the calculations the quantities n, A t, and E t were taken from the experimental data of [15, 16]. The values of the coordinate z at which f is given by the above-indicated equation were set equal to z = zo = 0.075, 0.i, 0.125, and 0.15. The oscillations diminished with an increase in zo, and for zo = 0.15 the solution is quite satisfactory. This is the solution analyzed below. The calculations show that the transition of a fairly narrow interval: 0.3 ~ z ~_ 0.45. In Fig. intermittency factor y and the ratio a/~c (curves experimental results of [15, 16]. It is seen that satisfactorily with the experimental data.

f to the asymptotic form f~ takes place in 4 the results of the calculations of the 2 and 4 respectively) are compared with the the results of the calculations agree

We note in conclusion that the approach suggested in [S] and adopted as the basis of the present study for describing the mixing process in turbulent flows entails the derivation of exact open equations for the probability distributions by means of phenomenological (and sometimes speculative) hypotheses involving semiempirical constraints and relations between the various statistical properties of the turbulence, including the concentration distribution function. This approach, therefore, ignores the problem of closing the equations and formulating a system of equations without empirical constants. It can essentially only be used to express some functions in terms of others, and its usefulness will be dictated strictly by what kind of functions must be specified (and this depends in large measure on the nature of the particular hypotheses) in order to find the unknowns. In the rigorous setting of the investigated problem the unknown characteristics, i.e., the intermittency factor and the probability distribution of the concentration, are related to no less simple functions: the

203

conditional-average velocity for a given value of the concentration, and the scalar dissipation. In this sense the exact value of the concentration , the scalar dissipation. In this sense the exact open equation for the concentration distribution has little to offer. Only when the hypotheses stated in Sec. 1 are invoked is it possible to determine such subtle characteristics as the concentration distribution and the intermittency factor from such rather "gross" characteristics of the turbulence as the fields of the average velocity, average concentration, and scalar dissipation (or in the solution of the inverse problem from the field of the variance of the concentration). The author is grateful to V. R. Kuznetsov for interest in the study, several discussions of the results, access to the processed experimental data of [31, 32] used in Fig. 3, and the idea for the method of deriving Eq. (3.6). The author also thanks G. I. Barenblatt for interest in the study and critical comments, as well as V. L. Zimont for suggestions incorporated during preparation of the manuscript for publication. LITERATURE CITED i. 2. 3. 4.

5. 6. 7. 8. 9. i0. ii. 12. 13. 14. 15. 16. 17. 18. 19. 20. 21. 22. 23. 24. 25. 26. 27. 28. 29. 30. 31. 32. 33.

204

F. G. Lockwood and A. S. Naguib, Combust. Flame, 24, No. i, 109 (1975). V. R. Kuznetsov,and A. B. Lebedev, Izv. Akad. Nauk SSSR, Mekh. Zhidk. Gaza, No. 1 (1977). V. L. Zimont, E. A. Meshcheryakov, and V. A. Sabel'nikov, Fiz. Goreniya Vzryva, 14, No. 3 (1978). V. L. Zimont, E. A. Meshcheryakov, and V. A. Sabel'nikov, in: Proceedings of the Fourth Tsanderov Lectures, Section on Engine and Aircraft Design [in Russian], IIET Akad. Nauk SSSR, Moscow (1978). V. R. Kuznetsov, Izv. Akad. Nauk SSSR, Mekh. Zhidk. Gaza, No. 5 (1972). C. Dopazo, Phys. Fluids, 18, No. 4 (1975). S. B. Pope, Combust. Flame, 27, 299 (1976). V. A. Frost, Izv. Akad. Nauk SSSR, Energ. Transport, No. 6 (1973). P . M . Chung, Phys. Fluids, 15, No. i0 (1972). J. Janicka, W. Kolbe, and W. Kollman, J. Non-Equilib. Thermodyn., 4, 47 (1979). V. Z. Kompaniets, A. A. Ovsyannikov, and L. A. Polak, Chemical Reactions in Turbulent Gas and Plasma Flows [in Russian], Nauka, Moscow (1979). E. E. O'Brien, AIAA Paper No. 137 (1980). V. A. Sabel'nikov, in: Proc. Sixth All-Union Symp. Combustion: Combustion of Gases and Natural Fuels [in Russian], Chernogolovka (1980). R. J. Bywater, AIAA Paper No. 139 (1980). P. Freimuth and ~i. S. Uberoi, Phys. Fluids, 14, No. 12 (1971). J. C. LaRue and P. A. Libby, Phys. Fluids, 17, No. ii (1974). M. S. Uberoi and J. Bashir, Phys. Fluids, 18, No. 4 (1975). S. Corrsin and M. S. Uberoi, NACA Report 1040 (1951). H. A. Backer, H. C. Hottel, and G. S. Williams, J. Fluid Mech., 30, Part 2 (1967). S. Corrsin and A. L. Kistler, NACA Report 1244 (1955). L. D. Landau and E. M. Lifshits, Continuum Mechanics [in Russian], GITTL, Moscow (1953). G. K. Batchelor, Theory of Homogeneous Turbulence, Cambridge Univ. Press, Oxford (1959). V. R. Kuznetsov and V. I. Rasshchupkin, Izv. Akad. Nauk SSSR, Mekh. Zhidk. Gaza, No. 6 (1977). V. A. Bezuglov, in: Proc. Twentieth Sci. Conf. Moscow Physicotechnical Institute [in Russian], Dolgoprudnyi (1975). A. A. Townsend, The Structure of Turbulent Shear Flow, Cambridge Univ. Press, Oxford (1956). M. A. Meshkov and Yu. A. Shcherbina, in: Abstr. Fifth All-Union Congress on Theoretical and Applied Mechanics [in Russian], Alma-Ata (1981). N. V. Kislov, Dokl. Akad. Nauk SSSR, 255, No. 1 (1980). H. Bateman and A. Erd41yi, Higher Transcendental Functions, McGraw-Hill, New York (1953). R. W. Hamming, Numerical Methods for Scientists and Engineers (2nd ed.), McGraw-Hill, New York (1973). P. A. Antonia, A. Prabhy, and S. E. Stephenson, J. Fluid Mech., 72, Part 3 (1975). A. D. Birch, D. R. Brown, et al., J. Fluid Mech., 88, Part 3 (1978). I. Ebrahimi, R. Gdnther, and F. Haberda, Fortsch. Flugwes., 43, No. 2 (1977). J. E. Carter, AIAA Paper No. 583 (1974).