Mean velocity equation for fluctuating flow

Accepted Manuscript Mean velocity equation for fluctuating flow Juergen Piest

PII: DOI: Reference:

S0997-7546(16)30085-1 http://dx.doi.org/10.1016/j.euromechflu.2017.08.009 EJMFLU 3198

To appear in:

European Journal of Mechanics / B Fluids

Received date : 9 March 2016 Revised date : 21 July 2017 Accepted date : 15 August 2017 Please cite this article as: J. Piest, Mean velocity equation for fluctuating flow, European Journal of Mechanics / B Fluids (2017), http://dx.doi.org/10.1016/j.euromechflu.2017.08.009 This is a PDF file of an unedited manuscript that has been accepted for publication. As a service to our customers we are providing this early version of the manuscript. The manuscript will undergo copyediting, typesetting, and review of the resulting proof before it is published in its final form. Please note that during the production process errors may be discovered which could affect the content, and all legal disclaimers that apply to the journal pertain.

Mean velocity equation for fluctuating flow Juergen Piest July 21, 2017 Abstract Arguments are presented that the momentum equation derived by Zwanzig-Mori technique of N-particle statistical mechanics is identical to the Reynolds equation for the mean velocity field in fluctuating hydrodynamics; thus, calculating the non-linear stress of the momentum equation is perhaps a means to avoid the closure problem of turbulence theory. In this paper, the third-order term of the friction force is calculated. Multilinear mode coupling theory is applied in order to obtain formulas for the third- and fourth-order equilibrium time correlation functions appearing in the expression. The force term is obtained as a multilinear convolution integral containing higher order gradients of the velocity field. As an application, circular jet flow for low Reynolds numbers is investigated. Numerical solutions have been obtained from an iteration starting with the corresponding solution of the Navier-Stokes equation. Deviations become distinct for Re > 30. The velocity field is compared with the result of an existing experimental investigation in the same Reynolds number range.

1

Introduction

This paper deals with an attempt to derive an equation for the mean velocity field in fluctuating flow. We call a flow fluctuating if a repetition under equal experimental conditions as a rule results in a somewhat different velocity field. The mean velocity field is defined as the arithmetic mean over a large number of repetitions (ensemble mean). Turbulent flow is the most important example of fluctuating flow; but for small Reynolds numbers there are forms which, though fluctuating, are not called turbulent. In the past, searching for the mean velocity equation was identical to investigating the closure problem of turbulence theory. This approach has a long tradition in hydrodynamics (We do not include entries to the references for this paragraph, since there exist excellent collections; see, e. g., [1]). Since its definition by O. Reynolds in 1894, several generations of scientists worked on this topic, some of them for large parts of their lives. Two lines of approach can be distinguished: Reasonable additional formulations have been invented in order to arrive at a practicable system for describing the turbulent process; this starts with the eddy viscosity concept of Boussinesq, already in 1877, but includes also the large eddy simulation (Deardorff 1970), which combines direct numerical simulation with a statistical averaging concept for small eddies and is one of the most used practical computation methods presently. - In the other line, one uses special mathematical techniques in the intention to obtain essential results in spite of the closure problem. To this belong the long list of papers of R. Kraichnan starting 1958 which employ the concept of homogeneous isotropic turbulence; but also direct numerical simulation if an averaging procedure is used. This technique is especially interesting since it displays turbulence properties from data which stem directly from solutions of the Navier-Stokes equations. The mean velocity field u in fluctuating flow obeys the Reynolds equation (see [2] eq. (1.12)), which reads, with some denotations more suitable for the discussion in this paper: ∂u + u · ∇u) = −∇P + ∇ · R (1.1) ∂t P is the mean pressure field and ρ the mass density. R is the dissipative stress which consists of the Stokes part and the Reynolds stress: ρ(

R = ρ (ν∇u − h(ˆ u − u)(ˆ u − u)i)

(1.2)

ˆ denotes the velocity field of a single outcome of the experiment (realization), and the pointed u brackets describe the ensemble mean (expectation) of the enclosed quantity. ν is the kinematic viscosity. 1

The Reynolds equation differs from the Navier-Stokes equation by the appearance of the Reynolds stress term. This is an additional dependent variable in the equation; thus, the equation is not closed and therefore is not sufficient to determine the mean velocity field. This constitutes the closure problem of turbulence theory in its simplest form. The problem is directly connected to using the N.-S. equation as the basis of the theory. In the derivation of the equation from statistical mechanics it is seen that the fluid velocity is an expectation. In turbulence theory, the quantity is redefined as a stochastic variable; but in this approach it is not possible to obtain a probability distribution, which would be necessary in order to derive a mean velocity equation. In this paper, (classsical) N-particle statistical mechanics is employed in order to obtain a closed equation. It is obvious that for investigation of turbulent processes it is not necessary to deal with single molecules of the fluid. It is the probabilistic structure of statistical mechanics which enables the analysis. Statistical mechanics is the basis for describing classical processes including fluid dynamics; it would be a surprise if turbulent motion would be an exception. Thus, it is assumed here that it is principally possible to derive the mean velocity equation by this approach. The complication is that in a statistical mechanics analysis it is necessary to distinguish between macroscopic and microscopic parts of the motion, and to handle the latter in a suitable way. - The Navier-Stokes equation can be derived from statistical mechanics via the Boltzmann equation by Chapman-Enskog method; see e. g. [3]. But in the course of the derivation of the Boltzmann equation the molecular chaos assumption is introduced, which is correct for many applications including laminar fluid motion, but seems problematic for an investigation of fluctuating flow. In this paper, the Zwanzig-Mori projection operator technique of statistical mechanics (ZMT) is utilized, which avoids this problem. The technique starts with the basic rules of statistical mechanics. Still, it is not a completely deductive method. In order to separate slow from fast parts of the process, it is necessary to define a suitable projection space of slow variables. In the ZMT, the conserved variables are used. In general, physical experience is necessary in order to design such a technique. An essential step is defining the mean velocity in the frame of statistical mechanics. In general, the fluid velocity is defined: 1 (1.3) u = hpi ρ p the microscopic momentum density. In this paper, we consider incompressible constant density and temperature processes. Then u is, up to a constant factor, equal to the expectation of p. This definition is independent from the structure of the flow; if a flow exhibits macroscopic fluctuations, the arithmetic mean will account for these. Thus, we must conclude that the statistical mechanics u is originally the mean velocity field of a fluctuating flow. We employ relevant parts of the ZMT in the presentation of Grabert [4]. As we will describe in more detail in the next section, the final form of the derived equation has exactly the form of the Reynolds equation (1.1), but now with a fomula for the Reynolds stress which is closed: R = ρ ν∇u + RN L

(1.4)

The derivation of the explicit formula for R is the first topic of this paper. Because of the relations described above, we may equate the mean velocity field in the Reynolds equation and the statistical mechanics definition of this quantity, and correspondingly the Reynolds stress in (1.2) and the nonlinear part RN L of the dissipative stress derived by the ZMT. The identity cannot be proved on purely theoretical grounds and should therefore be considered the basic hypothesis of this paper; it is necessary to gain evidence of its correctness by comparing numerical solutions of the equation with experimental results. In sec. 2 the result of the ZMT is presented. The formula contains a local equilibrium correlation function which is a nonlinear functional of the mean velocity field. Presently, a technique for evaluating this quantity in closed form does not exist. Therefore, a functional power series expansion in the velocity is used; the coefficients then contain total equilibrium time correlation functions which can be calculated by applying the multilinear mode-coupling theory (MCT) of Schofield and co-workers [5], [6]. The lowest order nonlinear term of the expansion, which is of 3rd order, is used as an approximation for the friction force. An algebraical calculation of this quantity is performed. This is a lengthy procedure. It contains the following steps: - calculation of the 2nd functional derivative of the original formula; - transformation of the total equilibrium correlations into standard form; 2

- calculation of the standard correlations; - evaluation for the Reynolds stress formula. These are presented in secs. 3 (for the first two steps), 4 and 5. In sec. 3, after calculation of the third-order term, the space integrals are transformed to Fourier space, and the correlation functions are brought into standard form so that they can be processed by the multilinear mode-coupling theory. In sec. 4 (in appendix D) a short description of the MCT is given, and the results of the application are presented. In sec. 5 the description is confined to shear flow (sound and heat flow effects neglected), and the appearing static correlations are calculated. - In sec. 6 the final form of the third order term is presented. Detailed calculations are postponed to the appendices. In secs. 7 to 11, an application of the equation to circular jet flow is considered. The stream function equation is used in order to calculate the velocity fields. As a first step, solutions of the Navier-Stokes equation are calculated for a set of Reynolds numbers. Solutions of the mean velocity equation are found as the lowest order step of an iteration procedure. Deviations from the N.-S. solutions become distinct for Re > 30. In sec. 11, the flow structure is compared to the result of an experimental investigation [13].

2

Mean velocity equation

The ZMT permits to derive, from the basic equations of motion of classical statistical mechanics, a closed system of equations for the expectations of the conserved mechanical quantities mass, energy and momentum. In [4], these are called the Generalized Transport Equations, which are derived in sec. 2.2 of this book and applied for simple fluids in sec. 8 where they take the form (8.1.13). The ’organized fluxes’ are developed further and result in (8.4.15). In this paper, the momentum component is taken which has the form (1.1), with a formula for the dissipative stress (the ’disorganized fluxes’): Z t Z 0 R(x, t) = dt dx0 S(x, x0 , t, t0 )∇0 u(x0 , t0 ) (2.1) 0

0

0

S(x, x , t, t ) = βh[G(t0 , t)(1 − P(t)s(x)](1 − P(t0 ))s(x0 )iL,t0

(2.2)

In (2.2), s(x) are the fluxes assoziated with the conserved quantities; P(t) is the projection operator which is the central operator quantity of the technique, and G(t0 , t) is an exponential operator; hiL,t denotes the local equilibrium expectation which is generally a function of time. These quantities and the essential steps of the derivation are explained in appendix A. - In Addition to (1.1), the mass component of the equation system yields the continuity equation: ∇·u=0

(2.3)

The kernel function S is a time correlation function in local equilibrium which is a nonlinear functional of u. As far as the author knows, presently it is possible to calculate correlation functions for total equilibrium only. Therefore, it has been necessary to expand S into a functional power series in u. The expansion can be done with the set b(t) of conjugate parameters of the local equilibrium formula (see app. A (A.8b)). The expansion is performed at the point b = b0 , which corresponds to u = 0: b0 = {−β

µ , β , 0} m

1 b − b0 = {β u2 , 0 , −βu } 2

(2.4) (2.5)

The power expansion of S reads:

S

δS 1 δ2 S |b0 ∗ (b − b0 ) + |b ∗ ∗{(b − b0 ), (b − b0 )} + · · · δb 2! δb δb 0 S (0) + S (1) + S (2) + · · ·

= S| =

b0

+

(2.6)

The denotation of terms in the second row is for later reference. The ∗ indicates multiplication, summation over five elements and integration over space and time. For the present purpose, the expansion is cut after S (2) . When these terms are inserted into (2.1), one obtains corresponding parts R(1) , R(2) , R(3) of R; the upper index describing the order in u. The elements of b0 are the conjugate parameters of the total equilibrium ensemble. Thus, when S and its derivatives are taken at b = b0 , the quantities 3

in the integrand resemble total equilibrium space-time correlation functions, which have been calculated using the multilinear mode-coupling theory (MCT) of Schofield and co-workers [5], [6].- In order to switch from the expansion in b to that in u, certain corrections have to be made which will be mentioned in due course. The expansion in u is actually an expansion in the Reynolds number Re of the flow. The linear part R(1) of the stress tensor has been calculated by several authors including Grabert [4], sec. 4; the Stokes form of the stress tensor is obtained, with a microscopic definition of the friction matrix. Thus, the stress tensor has the form (1.4). For small Re, the equation formally reduces to the Navier-Stokes equation. By detailed calculation, it can be shown that R(2) = 0 [7], so that the lowest nonlinear part of R is of third order. Determination of R(3) will be the topic of the next sections.

3

Third-order term of the dissipative force

Derivation of the third-order term provides considerable effort; a computer algebra system (Mathematica) has been used for most of the calculations. In order to describe the third-order term, we switch, in (1.1), from the stress tensor R to the friction force D = ∇ · R. By introducing the corresponding term of (2.6) into (2.1), we obtain a formula for the third-order term D 3 of D. For the first part of the investigation, where functional derivatives are calculated, we need a form of the formula with component indices and space variables explicitly shown:

(D3 )a (x, t)

=

β2 ∇c 2

Z

t

0

dt

0

Z

dx

0

Z

∞

00

dt

0

Z

dx

00

Z

∞

0

000

dt

Z

dx000

δ 2 Sabcd (x, x0 , t, t0 ) |b × δbe (x00 , t00 )δbf (x000 , t000 ) 0

× ∇0d ub (x0 , t0 )ue (x00 , t00 )uf (x000 , t000 )

(3.1)

Latin letters denote indices which run over 3 elements; by contrast, greek indices (where they appear) run from 1 to 5. - In order to continue, the second-order functional derivative of S is to be calculated. This is described in appendix B. The resulting formula consists of 11 terms ((B.9) to (B.22) including those which are zero). Taking these at b = b0 ; changing from the conserved quantities and their fluxes a, s to the corresponding ortho-normalized quantities h, r, which alters the pre-factor of the formula; and performing some substitutions on the exponential operators yields (B.27) to (B.37). Gradient operators appearing in (3.1) are transferred to the kernel function; by partial integration, this alters the sign of the expression. We change to Fourier space to prepare for the incorporation of the results of the MCT. Here, we can switch back to a concise form of the formulas where the wave numbers are not explicitly shown; this is necessary to keep the presentation in appendix C reasonably short. The expression now reads: (D3 )0 (t) = −

1 2 ρ β 2

Z

t

0

dt0

Z

0

∞

dt00

Z

∞

dt000 K0123 (t, t0 , t00 , t000 )u1 (t0 )u2 (t00 )u3 (t000 )

(3.2)

0

We have introduced number indices. A number index is a combined index which contains a component index and a wave number variable. Moreover, double indices contain a summation over the indexed 1 elements and an integration over the corresponding wave number space, together with the factor (2π) 3. The kernel function K consists of the 11 Terms K1 to K11 (B.43). The δ- and θ-factors are temporarily omitted; when the time integrations start they will be added again. The correlations in these formulas are to be substituted by those which can be evaluated be the MCT; these are the Cn (C.12). The relations are elaborated in appendix C. Now the kernel function in (3.2) consists of 5 quantities M1 to M5 ; the δ- and θ-factors are reintroduced (C.39), (C.41). When (C.39) is inserted into (3.2), some of the time integrals can be executed. Moreover, we can apply an approximation since the time scale of the u-factors is much larger than that of the equilibrium time correlations. The upper limits of integrals over correlation functions can then be extended to infinity. The formula now reads:

(D3 )0 (t)

=

1 2 ρ βu1 (t)u2 (t)u3 (t) × 2 Z ∞ Z t Z ×( dt0 Φ10123 (t0 ) + lim dt0

−

t→∞

0

4

0

0

t−t0

dt00 Φ20123 (t0 , t00 ))

(3.3)

with Φ1 = M1 + M3 and Φ2 = M2 + M4 + M5 . The detailed expressions are given in the appendix, (C.42). - This is the part of the third-order term which stems from the term S (2) of the kernel function in expression (2.6) for the stress tensor. As has been mentioned earlier, still to be added is the part of the first-order term in (2.6) (the second-order term of D) with the index value  = 1 for b − (b0 ) . When one treats the second-order term in the same way which led to (3.2), one obtains:  − 1 Z t Z ∞ 1 3 12 ∂p 2 2 dt00 (K2 )012 (t, t0 , t00 ) u1 (t0 )uq(t00 ) dt0 (D2 )0 (t) = − ρ β 2 ∂ρ 0 0

(3.4)

In the rest of this section, the number index 2 contains a component index which has the fixed value 1, if not stated otherwise. The pre-factor contains a part which results from switching to normalized variables in the correlation functions. uq is the Fourier transform of u2 for the wave number variable k00 . The kernel function K2 has been calculated in appendix C, (C.43) to (C.45). In (3.4), integration over t00 and an approximation because of different time scales is performed: Z ∞ 1 3 1 −1 dt0 (K2 )012 (t0 ) (3.5) (D2 )0 (t) = − ρ 2 β 2 (∂ρ P |β ) 2 u1 (t)uq(t) 2 0 K2 , now a function of t0 only, is given in (C.46). - (3.5) has to be added to (3.3) in order to obtain the full third-order term D3 of the friction force:

(D3 )0 (t)

=

1 2 ρ β(u1 (t)u2 (t)u3 (t) × 2 Z ∞ Z t Z t−t0 ×( dt0 Φ10123 (t0 ) + lim dt0 dt00 Φ20123 (t0 , t00 )) t→∞ 0 0 0 Z ∞ − 21 0 +(ρβ∂ρ P |β ) u1 (t)uq(t) dt (K2 )012 (t0 )) −

(3.6)

0

To complete the derivation, the equilibrium correlation functions C3 and C4 appearing in the coefficients have to be calculated. This will be described in next section.

4

Multilinear mode coupling theory

In order to calculate the correlation functions, the MCT has been applied [5], [6]. An expression for C3 is explicitly given in [6], while C4 has been calculated with the aid of the rules given in these papers. Both formulas are given in terms of C2 , together with a differential equation for the latter quantity; at a later stage, the solution for C2 is inserted. - Only the final expressions are given in this section. A short exposition of the technique and some more detailed formulas can be found in appendix D. For C3 , the formulas derived in [6] are (30), (31); with the denotations of the present paper and in terms of number indices, they read: (C3 )123 (t) = (C2 )14 (t)(j3 )423 + (G12 )123 (t) (G12 )123 (t) = −

Z

0

t

dt0 (C2 )14 (t − t0 )(s3 )456 (C2 )25 (t0 )(C2 )36 (t0 ) (s3 )123 = hˆ r1 h∗2 h∗3 i

(4.1) (4.2) (4.3)

j3 is described in (C.29). Projected flux densities rˆ are defined in (B.26). s3 is somewhat different from the quantity defined in (C.32); but the difference affects only the δ-factor of the correlation function explained in appendix B. We will use the same notation. In the beginning of appendix D, a short introduction into the multilinear mode coupling theory is presented. The formulas used for the calculations in this paper are (D.5), (D.6). The formulas for C4 obtained from these for β = 3 read, after some manipulations: (C4 )1234 (t) = ((C3 )156 (t) − (j3 )156 )Φ23456 + (C2 )17 (t)(j4 )7234 + (G13 )1234 (t) (G13 )1234 (t) = −

Z

0

t

dt0 (C2 )15 (t − t0 )(σ 4 )5678 (C2 )26 (t0 )(C2 )37 (t0 )(C2 )48 (t0 ) 5

(4.4)

(4.5)

The constant matrices j3 , j4 and σ 4 are described in appendix D. C2 (t) is, in [6], the quantity G11 called the propagator. The equation for C2 is described in appendix D. Its exact form is needed if time derivatives of are to be calculated; otherwise, a simplified form is valid: ∂t (C2 )12 (t) = −κ13 (C2 )32 (t) For the coefficient matrix κ, again see appendix D. The solution of C2 matrix exponential function: (C2 )12 (t) = e−κ12 t

5

(4.6) can be expressed by the (4.7)

Evaluation of the correlation formulas

In order to keep the formulation reasonably short, in the past sections of this paper, a rather formal description has been used employing number indices; double indizes indicate a component summation and a wave number integration, together with a factor (2π)−3 . While this is reasonable for formal calculations, for the rest of this paper we will switch to an easily readable form where wave numbers are explicitly shown. The remaining indices are then component indices. In this section we will continue to denote them by numbers, instead of letters. The reason is that in the formulas to be described some indices would run from 1 to 3 (so they would be latin), while others initially would be greek, and during the analysis would change to latin. This would lead to a cumbersome notation. The explicit formulas are presented in appendix E. For all correlation functions, we will extract the δ-factor described in (B.39), denoting the remaining part with the same symbol as before. It then becomes apparent that the additional wave number integrals introduced via the definition of the projection operators (B.40), (B.41), all can be executed. The explicit formula for D 3 is (E.1) to (E.4), with (E.5) to (E.7) for the intermediate quantities Φ1, Φ1, K2 . The explicit expressions for the correlations in these quantities are (E.8) to (E.13). The formula for ∆3 contains time derivatives of the correlations C3 , C4 ; these can be substituted in such a way that they do contain the quantities G12 , G13 but not their time derivatives. This is shown in appendix E, (E.12) to (E.23). The results are inserted into Φ1, Φ2. It is seen that terms still containing G12 , G13 cancel each other. The remaining terms are products of C2 of different order (in Φ2 with different time variables) multiplied by certain coefficients. Now we are able to perform the time integrations in (E.2). Some definitions and intermediate results are presented in (E.24) to (E.28). (E.29) is the expression obtained for the kernel function ∆3 . After combining some similar quantities, it consists of 10 terms. A problem arises with terms no. 7 and 9 where wave number integrals appear which contain static correlations. Here it is not possible to localize the integrand factors totally since then the integral will diverge. We will discuss this in appendix E in the text after (E.30). I sometimes call a quantity greek or latin, depending on which type its indices are. In the present state, in ∆3 indices 0 to 3 are latin, all others are greek. To proceed further, we make use of the mode analysis of hydrodynamics. The greek matrix κ is identical to the matrix of linearized hydrodynamics; see, e. g., [10]. It is symmetric in its indices and therefore can be decomposed into hydrodynamic modes κ ˆ: κ12 = χ13 χ23 κ ˆ3 (5.1) In this formula, index 3 which appears three times indicates a summation. χ is the modal matrix which facilitates the decomposition. The modes are: two sound modes, a heat mode and two identical shear modes. Since in this investigation we are not interested in sound and heat effects, we neglect all parts of κ connected to the first three modes. Then κ simplifies considerably: It now is latin, and it reads: κ12 (k) = −ε12 (k) ν ε12 (k) = k 2 δ 12 − k1 k2

(5.2)

ν is the kinematic shear viscosity. - Now that κ is latin, this also applies for matrices ζ 2 , ζ 3 , ζˇ3 , ζ 4 and their inverses. One then finds that now all component indices of (E.29) are latin, with the exception of indices 5, 6 in the wave number integrals in terms 7 and 9. Static correlations can be calculated with a method presented in [11] which is described in appendix E. This is then applied to the static correlations appearing in (E.29); and an argument is given that the two terms with wave number integrals can be neglected. The remaining formula for ∆3 is (E.31). Finally it is shown in appendix E that the second-order term (the last row of (E.1)) just cancels the first row of (E.31). This completes the evaluation of the results of the MCT.

6

6

Final form of the 3rd order term

For this final formulation, the third order friction term D 3 is just denoted D. The mean velocity equation (1.1) then reads: ρ(

∂u + u · ∇u) = −∇P + D ∂t

(6.1)

Since the remaining component indices are all latin, we now denote them by latin letters. The final form of the 3 rd order term is (F.3), (F.4). In order to transform these expressions back to geometric space, is seems reasonable to switch to a different notation:

Da (k, t) =

ρ ν

Z

dk0 dk00 dk000 δ(k − k0 − k00 − k000 )

M1

=

M2

=

2 X

Mµ (k0 , k00 , k000 )(Qµ )a (k0 , k00 , k000 , t)

(6.2)

µ=1

1 (k0 + k000 )2 (k 002 + (k0 + k000 )2 ) 1 (k0 + k000 )2 (k 002 + (k0 + k000 )2 )(k 02 + k 0002 )

(6.3)

(Q1 )a = −i(kg0 + kg00 + kg000 )(s3 )aceg (ikh0 )(s3 )bdf h εef (k0 + k000 )ub (k0 , t)uc (k00 , t)ud (k000 , t)

(6.4)

(Q2 )0 = i(kg0 + kg00 + kg000 )(s3 )aceg i(kh0 + kh000 )(s3 )f bdh k 0002 εef (k0 + k000 )ub (k0 , t)uc (k00 , t)ud (k000 , t)

(6.5)

In geometric space, we have: Da (x, t)

ρ ν

Z

dx0 dx00 dx000

2 X

µ=1

Mµ (x − x0 , x − x00 , x − x000 )(Qµ )a (x0 , x00 , x000 , t)

(6.6)

The functions in this formula are the (backwards) Fourier Transforms of the equally denoted functions in (6.2). The transform of the kernel functions (6.3) can be found in closed form: M1 (x0 , x00 , x000 )

=

M2 (x0 , x00 , x000 )

=

˜ x x∗ x ˆ

x0 1 arctan( )δ(x0 − x000 ) 24 π 3 x0 x00 x00 1 x ˜ arctan( ) 8 4 00 ∗ 2 π x x ˜x x ˆ 1 0 (x − x000 ) 2 1 0 = (x + x000 ) 2 = x00 + x∗

(6.7)

=

(6.8)

As usual, vectors are in bold letters, while their magnitude is in normal letters. The velocity functions (6.4), (6.5) are transformed straightforwardly: 0 0 000 0 00 000 (Q1 )a = (s3 )aceg (s3 )bdf h (∇0g + ∇00g + ∇000 g )∇h εef (∇ + ∇ )ub (x , t)uc (x , t)ud (x , t)

(6.9)

0 000 0002 (Q2 )a = (s3 )ace (s3 )f bdh (∇0g + ∇00g + ∇000 εef (∇0 + ∇000 )ub (x0 , t)uc (x00 , t)ud (x000 , t) (6.10) g )(∇h + ∇h )∇

εbc (∇) = ∇2 δ bc − ∇b ∇c Expansion of (6.9), (6.10) leads to lenghty expressions which will not be presented here. 7

(6.11)

7

Application to circular jet flow

It is desirable to apply the formulae presented in the preceding section in order to either compare with the results of experiments, or at least have some numerical results for testing the theory. For such an application, several restrictions have to be observed. Important is the Reynolds number range. There is a range where the nonlinear dissipation term is small against the linear one. This is where the mean velocity equation takes approximately the form of the Navier-Stokes equation. Because of the power series approximation of the friction term, the theory should be applied for Reynolds numbers which do not excede too far this range. There is another restrictive feature of the theory in its present form: The statistical mechanics formalism is, as usual, for a fluid with no boundaries; the same is true for the correlation function calculations. In Navier-Stokes calculations, boundaries are accounted for by boundary conditions, but do not change the form of the equation. Here the situation is different: The existence of a boundary should influence the formulae for correlation functions, and via the dissipation term alter the form of the equation itself. The author has attempted to calculate the stationary flow in a tube with the formulae in their present form, and obtained the result that the third-order friction term is generally zero for this configuration. This apparently means that the theory in its present form cannot be applied to tube flow. The possibility to calculate correlation functions in the presence of boundaries is one of the major necessities for making the theory applicable. The theoretical formulae have been tested for the circular jet as a representative flow, since it is thought that sufficiently far from the nozzle boundary effects may be relatively small. In addition, the circular jet is well studied in models and experiment. It seemed possible to use an existing empirical formula for the jet velocity field as a test formula in order to investigate the various terms of the mean velocity equation, and perhaps find that the test formula is a solution. For small Reynolds numbers there are essentially the formulas by Schlichting [12] and by Schneider [14]. The Schlichting formulas are considered valid, with suitable coefficients, for laminar as well as for turbulent flow. For laminar flow, they are a solution of the Navier-Stokes equation under the condition of constant momentum flux. Schlichting used a point source at the orifice of the jet; this can be altered by introducing an ”effective point source” as has been installed by Andrade et al. [16], in order to secure convergence of the integrals of the momentum equation in the neighbourhood of the orifice. But with the Schlichting formulae, the integrals definitely diverge in the far region of the field. Thus, the present theoretical equations are in contradiction with the Schlichting model. - Schneider’s formulae are for laminar flow, i. e. for rather small Reynolds numbers. They are the result of a theoretical investigation which leads to a slowly varying momentum flux. The formulae show good agreement with the experimental results of Zauner [13] which confirm that for these Reynolds numbers the round jet has essentially a finite length and is part of a recirculation pattern. With these formulae, the integrals converge in the far field. But here another problem arose: The Schneider model again uses a point source at the orifice. The author did not succeed in substituting this by an effective point source in such a way that the velocity field behaves reasonably in other parts of the flow space. The solution of the mean velocity equation for D = 0 can be used as a test field in order to investigate the properties of the nonlinear dissipation force for small Reynolds numbers. The equation then has the form of the Navier-Stokes equation, and the solutions are the velocity fields of laminar fluid motion. In sec. 8, cylindrical coordinates are introduced. In sec. 9, numerical calculations of velocity fields for a set of Reynolds numbers are reported. Calculations of jet flow have already been performed 1985 by W. Schneider and co-workers [17]; many details, especially concerning boundary conditions, have been adopted from this paper. In the present paper, the stream function equation is used for the calculation, and the velocity field is calculated from the stream function. Some properties of the resulting flow fields are discussed and compared with results of publications mentioned before. In sec. 10, using these results, the nonlinear friction term is numerically calculated as a function of axial and radial coordinates. The computation is impeded by the large scatter of the data, the reason for which is perhaps a consequence of using the stream function equation instead of the N.-S. equation but is not well understood presently. In sec. 11, an approximation for the integro-differential equation (6.1) is formulated, which enables a numerical solution; the differences to the solution of the Navier-Stokes equation are considered. The computed velocity field is compared to the result of an experimental investigation in the same range of Reynolds numbers [13]. Some qualitative aspects mentioned in the paper of A. J. Reynolds [18] are discussed. 8

8

Evaluation of the nonlinear dissipation term

Generally, formulas (6.1), (6.6) are 3 nonlinear partial integro-differential equations of the Fredholm type; together with the continuity equation, they form a system of 4 equations for the three components of the mean velocity and the mean pressure. The formula for the nonlinear friction force contains, within the static correlation c3 , the quantity λ (E.34) which is a physical constant of the fluid: α γ cV

λ=

(8.1)

Here, α is the thermal expansion coefficient, γ the isothermal compressibility, cV the specific heat for constant volume. For standard conditions, we have λ = 0.40 for air and λ = 0.091 for water. By this quantity, the results depend to a certain extent on the physical properties of the fluid. The calculations presented in this paper are for a hypothetical fluid with λ = 0. For the application in mind, all formulas have been transformed to cylindrical coordinates and specified to cylinder-symmetric conditions. For most circular jet investigations, boundary layer approximations are considered to be valid; that is, radial components are small compared to axial components, and axial gradients are small to radial ones. In a jet, this will be true for regions not too far from the nozzle. In the present investigaton, the space integrals appearing in (6.6) extend over the whole flow area; but major contributions to the integral will be from the region in the neighbourhood of the axis and not too far from the nozzle. Thus, boundary layer conditions may still be a reasonable approximation. The formulas then considerably simplify. An orifice of diameter d exhibiting a parabolic velocity distribution with center velocity u0 is introduced. In the equations, all quantities are made dimensionless by using as parameters d and u0 . Then, for any choice of the parameters, in dimensionless coordinates the orifice width is 1. The third order term, instead of νρ , is now proportional to the Reynolds number Re = uν0 d . - It will be explained in sec. 10 that for detailed investigations we will only need the axial components of the two terms of (6.6), now called D1, D2 which are now functions of the axial and radial coordinates x, r. D1 simplifies from the Dirac function in (6.7):

Re D1(x, r) = 16π 3

Z∞ 0

dx

0

Z

0

∞

dr

0

Z

∞

00

dx

0

Z

0

∞

dr

00

Z

2π

0

dϕ

0

Z

2π

0

× z0 = z 00 =

p

p

dϕ00 ×

 0 z r0 r00 arctan( )Q1(x0 , r0 , x00 , r00 ) 0 00 zz z 00

(8.2)

(x − x0 )2 + r2 + r02 − 2rr0 cos(ϕ0 )

(8.3)

(x − x00 )2 + r2 + r002 − 2rr00 cos(ϕ00 )

(8.4)

Q1 is the axial component of the function (6.9) in cylindrical coordinates; it will be further specified in sec. 10. The formula for D2 has been elaborated in the same way but will not be given here.

9

Numerical calculation of laminar jet flow.

As is stated in the introduction, the author performed numerical calculations of the velocity field of circular jet flow as stationary solutions of eq. (6.1) for D = 0 which has then the form of the NavierStokes equation. These data are needed for the calculation of the nonlinear friction force (6.6) in sec. 10 which in turn will be used for an approximate solution of (6.1) in sec. 11. We select, among others, the Reynolds numbers chosen by W. Schneider and co-workers [17] so that we can compare the results. The computations in [17] are for momentum Reynolds numbers Rm = 3.5, 5.5: Rm =

√

m ν

(9.1)

m is the momentum flux at the orifice divided by 2π. The relation to Re is: √ Re = Rm 24 9

(9.2)

So, the numbers from above correspond to Re = 17, 27. For the present purpose, it is necessary to repeat the calculations for some other Reynolds numbers. - The author used the personal computer algebra system Mathematica. He intended to perform the calculations with the aid of the NDSolve command. There is now a choice between the so-called line method which uses a rectangular net with constant mesh size, and the finite element method. But in NDSolve, the finite element method is presently restricted to linear equations of at most second order. Thus it was necessary to use the line method. This method cannot process equations of the elliptic type; so it was not possible to use the stationary Navier-Stokes equation. It was decided to solve the instationary equation up to a time where further variations with time could be neglected, so that the final solution could be declared to be the stationary one. We switch to non-dimensional quantities, as is already done in the nonlinear dissipation term. Especially, the time is made non-dimensional with ud0 . For circular symmetric flow, in cylinder coordinates there exists the Stokes stream function ψ: u=

1 ∂r ψ r

(9.3)

1 v = − ∂x ψ (9.4) r To the Navier-Stokes equation, there corresponds the stream function equation (Milne-Thomson, [19], p. 647): ∂t ζ −

∂(ψ, 1r ζ) 1 1 2 − E (rζ) = 0 ∂(x, r) Re r 1 ζ = − E2ψ r

1 E 2 f = ∂x,x f + r∂r ( ∂r f ) r with ζ being the third component of the curl: ζ = (∇ × u)3

(9.5) (9.6) (9.7)

(9.8)

The different sign of the second term in (9.5) compared to [19] is due to the different sign in the definition of ψ by Milne-Thomson. Corresponding to equation (6.1), we have the stream function equation: ∂t ζ −

∂(ψ, 1r ζ) 1 1 2 − E (rζ) = ∆ ∂(x, r) Re r ∆ = (∇ × D)3

(9.9) (9.10)

Instead of determining u, v directly, the stream function has been computed by solving equation (9.5). The calculation is performed for a square of dimensionless length a in (x, r)-space. The boundary conditions have been copied from [17]. In dimensionless two-dimensional x-r-coordinates, the orifice has the width 21 ; thus for the wall we have the condition:  1 − 4r2 , 0 ≤ r ≤ 0.5 u(0, r) = (9.11) 0 , 0.5 < r On the jet axis, the conditions read:

∂r u(x, r)|r=0 = 0 v(x, 0) = 0

(9.12)

For the stream function ψ(t, x, r), the calculation starts from rest: ψ(0, x, r) = 0 Conditions (9.11), (9.12) are transformed for ψ:  r2 2 −r ψ(t, 0, r) = 1 16

10

, 0 ≤ r ≤ 0.5 , 0.5 < r

(9.13)

(9.14)

(r∂r,r ψ(t, x, r) − ∂r ψ(t, x, r))|r=0 = 0

(9.15)

ψ(t, x, 0) = 0

(9.16)

Unfortunately, it is necessary to prescribe conditions for the outer boundaries of the square also. This is a consequence of the calculation being limited to a finite area while the true motion extends over the total positive x, r-space. Some estimate for the stream function has to be found for the outer boundaries, for a suitable chosen square width a, so that the solution can be trusted for, say, x, r < a2 . Schneider and co-authors used the analytical stationary solution for the linearized Navier-Stokes equation (creeping flow), which they took from [20], p. 153; if the coordinates x, r are replaced by polar coordinates σ, θ, the stream function is independent of the radius σ: ψ(σ, θ) = c(1 − cos(θ)3 ) The constant c is determined from the value of ψ for θ = conditions for the outer boundaries are inferred:

ψ(t, a, r) =

from (9.14): c =

1 a3 (1 − 3 ) 16 (a2 + r2 ) 2

∂x ψ(t, x, r)|x=a = −

ψ(t, x, a) =

π 2

(9.17)

3r2 a2 5

16(a2 + r2 ) 2

1 x3 (1 − 3 ) 16 (x2 + a2 ) 2

∂r ψ(t, x, r)|r=a =

3ax3 5

16(x2 + a2 ) 2

1 16

. - From this, the

(9.18)

(9.19)

(9.20)

(9.21)

The boundary conditions have been somewhat adapted for t < 1 in order to conform them to the initial condition (9.13). - Calculations were performed for a set of Reynolds numbers. Values for a started with the values chosen in [17]; a suitable value for the end time te of the instationary calculation was found by determining the coordinates of the center of the toroidal eddy for a set of end times and observe for which t-values they become steady. Since from now we only consider results obtained for the final time, we cease writing the time as an argument. - A disturbing feature of the results was that the value u(0, 0) = 1 prescribed by the boundary conditions is not returned by the calculations; for most Reynolds numbers, it was considerably smaller. It became apparent that the computed value depends strongly on the width a chosen for the calculation, but nearly not on the Reynolds number. Smaller a led to a larger center velocity. There is a optimal value a0 = 56.4 where the boundary condition is met exactly, again nearly independent of the Reynolds number. - In a second run, the procedure was somewhat altered: After the calculation with individually chosen a, a > a0 (the ‘outer’ field), a second run with a0 was performed, where for the outer boundaries the results from the first run were taken (’inner’ field); finally, the ψ field was composed from the two results by taking it from the inner field for x, r < a0 and from the outer field for the rest of the area. - As an example, the result for Re = 27 is shown, for comparison with the result in [17]; the data are a = 200 and te = 100000 . Fig. 1 is the contour plot of the stream function, which shows the stream lines of the flow: 11

70 60 50 40 30 20 10 0 0

10 20 30 40 50 60 70

Fig. 1: Contour plot of the stream function ψ(x, r), Re = 27; dimensionless units The similarity to [17] Fig. 5 is visible. For this rather small Re, the jet spreads widely and produces a dominant toroidal eddy. The distance of the maximum from the orifice, rm , is somewhat larger here.The spacial range is chosen such that it can be seen how the streamlines are composed at the outer boundaries of the inner region. The irregularities of the stream lines in the neighborhood of the orifice have been produced by the program for preparing the contour plot. - It is interesting to compare the values rm obtained for several Re with a theoretical formula given by W. Schneider [14]: Re2 ) (9.22) 192C 2 Fig. 2 shows rm values from calculated velocity fields for a set of Reynolds numbers against the theoretical graph with fitted constants b, C 2 : rm = b exp(

rm 25 20 15 10 5 0

5

10

15

20

25

30

Re

Fig. 2: Distance of eddy center from orifice rm in dimensionless units, against Reynolds number Re The constants of best fit for the calculated values here are b = 4.02, C 2 = 2.31. b remains undetermined in Schneider’s theory, C 2 is a value between analytical and experimental valus in [17] Fig. 6. Thus, for this Reynolds number range, the calculated values fit reasonably in Schneider’s theoretical work. For larger Reynolds numbers, the method of composing the stream function field cannot be used any more: In the inner region, the eddy center starts to ’feel’ the artificial right boundary and gets distorted. Then it is presumably necessary to employ a solver which permits variable net width. In experimental jet investigations, often the axial velocity at the jet axis is considered. Fig. 3 shows the behavior due to several authors, in nondimensional variables, for Re = 27: u 1.0 0.8 0.6 0.4 0.2 0.0

0

5

10

15

20

x

Fig. 3: Axial velocity at the jet axis u against distance from the orifice x; dimensionless units. Re = 27 The dashed line is from Schlichting’s formula, with a suitably introduced effective point source in order to account for an orifice with finite width. The full lines are from Schneider’s theoretical investigations which interpret the experimental results in [13], and the result of the present calculation. It is not possible to evaluate the latter data directly from the stream function at r = 0, because of the factor 1r in (9.3), (9.4). From visual interpretation, it became apparent that the values u(x, ra ) with ra = 0.0001 can be used as the axis velocity. - The full lines are only marginally different. 12

10

Numerical calculation of third order term using the results for laminar flow

We are now in the position to perform a calculation of the third order term using the numerical results of the preceding section. This involves the calculation of the space integral in eq. (8.2). The first attempt led to a numerical divergence of the integral in the neighborhood of the jet axis. A preliminary investigation has been performed in order to explore in which way this problem can be avoided. From the velocity function Q1 in (8.2), a single example term Q1x has been chosen: Q1x = u(x00 , r00 )∂x u(x0 , r0 )∂x,r,r u(x0 , r0 )

(10.1)

We used the results for Re = 20; for this preliminary calculation, it was sufficient to use the data of the outer field. The gradients appearing in (10.1) were computed as Mathematica InterpolatingFunctions, and then it was attempted to calculate the corresponding partial term of D1 (8.2) with the aid of the Mathematica command NIntegrate. The result clearly showed that the integral diverges already for this trial term. - It was presumed that because of the 1r factors in (9.3), (9.4) some numerical errors have blown up in the neighborhood of the jet axis. As an example, we consider the Stokes friction term in the mean velocity equation, which reads in dimensionless variables and in boundary layer approximation: la(x, r) =

1 1 ∂r (r∂r u(x, r)) Re r

(10.2)

la has been calculated for several values of x as a function of r. Fig. 4 shows two curves for x = 5.51 and x = 5.52: la 0.001 -0.001

0.2

0.4

0.6

0.8

1.0

r

-0.002 -0.003 -0.004 -0.005

Fig. 4: Stokes force for x = 5.51, x = 5.52 as function of r; dimensionless variables. Re = 20 It is seen that the curves exhibit irregular singularities at the jet axis. For r = 0, la should have the finite value −0.0036. This behavior is the reason for the divergences. The prescriptions finally chosen for smoothing the velocity components in the neighborhood of the axis are:  u(x, r) , r ≥ 0.0001 us (x, r) = (10.3) u(x, 0.0001) , r < 0.0001  v(x, r) , r ≥ 0.0001 vs (x, r) = (10.4) 0 , r < 0.0001 The smoothed components are called u, v again. After this manipulations, the computed integrals converged. - As a first application, the first part D1 of the third order term (formula (8.2)) has been calculated. The velocity function originally consisted of two parts: Q1 = q0 + q1 cos(ϕ0 − ϕ00 )

(10.5)

Both parts have been computed separately, for x = 30, r = 5. The second term amounted to about 0.5% of the first one. For further calculations, this second part has been neglected. This has already been allowed for in the form of (8.2). For q0 a lengthy formula is obtained which is presented in appendix G. - Then, the two parts of the third order term D1, D2 mentioned in sec. 8 were calculated separately, on the jet axis for a set of x-values, for Re = 35 as an example. Larger values for D1 are of the order 10−3 to 10−4 , while D2-values generally are of 10−7 . Thus, D2 will be neglected in the numerical calculations. In order to solve the stream function equation for the mean velocity field (9.9), we need to calculate the quantity ∆ (9.10): 13

∆ = ∂x D2 − ∂r D1

(10.6)

D1 , D2 are the axial and radial components of the friction force. In boundary layer approximation, the first term will be two orders lower than the second, so that we approximate:

∆ = −∂r D1

(10.7)

∆ = −∂r D1

(10.8)

Since D2 is neglected, the final formula is:

where D1 is given by (8.2). ∆ has been calculated by differentiating analytically the kernel function in (8.2), obtaining:

∆(x, r) = −

DM 1 =

Re 16π 3

Z

0

∞

dx0

Z

∞

dr0

0

Z

∞

dx00

0

Z

∞

dr00

0

Z

2π

0

dϕ0

Z

2π

dϕ00 DM 1 q0

(10.9)

0

 0  00 z z0 z r0 r00 (− arctan( )( 0 (r − r0 cos(ϕ0 )) + 00 (r − r00 cos(ϕ00 )))+ 02 002 00 z z z z z z 002 (r − r0 cos(ϕ0 )) − z 02 (r − r00 cos(ϕ00 )) + z 02 + z 002

(10.10)

with z 0 , z 00 by ( 8.3 ,8.4) and q0 given in appendix G. - As an example, the calculation is described for Re = 30. The parameters are a = 300, te = 50000. ∆ has been calculated on a horizontal cross section at r = 10, and two vertical sections at x = 12 and x = 23.5 . Fig. 5 shows the data points for the horizontal cross section:

Δ 0.0010

0.0005

0.0000

10

20

30

40

50

60

70

x

-0.0005

Fig. 5: Calculated values for ∆ at horizontal cross section for r = 10 Immediately seen is the scatter of the points. The author presumes that the numerical calculation errors in the neighborhood of the jet axis have not been eliminated sufficiently by the prescriptions (10.3), (10.4). - In order to cope with these errors, it would perhaps be preferable to base the calculation on the Navier-Stokes equation itself instead the stream function equation. - Moving averages (9 points each) have been calculated, and a smooth line has been drawn to interpolate the points, as seen in Fig. 6:

14

Fig. 6: Moving averages of fig. 5; values of third-order term ∆ against horizontal variable x The cross point with the vertical cross section at x = 10 has been marked. The slightly positive values for x ≥ 40 have been ignored. The line has then been scanned with the Mathematica coordinate tool, to obtain (in this case) a set of about 50 points. In the same way, the vertical cross sections have been processed. These data are the basis for an estimate, by the Mathematica Interpolation command, for the function ∆(x, r) in the range 0 ≤ x, r ≤ 300. On the 4 boundaries, ∆ = 0 has been prescribed. The behavior of the obtained function is somewhat ragged, but is still a well-defined set of points in the prescribed square.

11

Approximate solution for the stream function of the mean velocity equation

In the stream function equation (9.9), on the right side we introduce, for ∆, the function determined in the preceding section. The equation can then be viewed as the first equation of an iteration procedure for calculating the stream function. It has the same structure as the equation for laminar flow, with an additional inhomogeneous term on the right. This equation can be solved with the NDSolve command. First, the calculation has been performed with the parameters for Re = 30, in the same way as described in sec. 9. The contour plots of the solutions look very similar. For the velocity at the axis, as a function of x, there is practically no difference between the two calculations. - The calculations have been repeated for Re = 32. Since here some problems of composing the two fields described in sec. 9 are encountered, only the stream functions of the original (outer) fields are compared. Figs. 7 and 8 show the solutions which belong to the Navier-Stokes and to the mean velocity equation: 350 300 250 200 150 100 50 0 0

50

100 150 200 250 300 350

Fig. 7: Stream function for Re = 32. Solution of eq. (9.5)

15

350 300 250 200 150 100 50 0 0

50

100 150 200 250 300 350

Fig. 8: Stream function for Re = 32. Solution of eq.(9.9) It is seen that the third order friction term now provides a distinct distortion of the flow field far from the nozzle. The jet is of finite length; for x > 100, the velocity on the axis is slightly negative, i. e. the flow is oriented towards the nozzle. u 0.00015 0.00010 0.00005

80

100

120

140

x

Fig. 9: Velocity at the jet axis u as a function of distance x. Upper curve: solution of (9.5) Fig. 9 shows the velocity at the jet axis as a function of distance, calculated for the outer field. Far from the nozzle the deviation from the solution of eq. (9.5) (upper curve) is distinct. As already mentioned, for x > 100 the velocity becomes negative. - Finally, a calculation has been performed for Re = 35. Here the stream function pattern shows some details which are considered to be not physically relevant. This demonstrates that for higher Reynolds numbers it is inevitable to obtain a numerical solution of eq. (6.1). To the knowledge of the author, there is only a single experimental jet investigation for this Reynolds number range by E. Zauner ([13]). Fig 10 is the fig. 3 of this paper which shows a flow in water for Re = 34. The visualization is with tiny air bubbles and illlumination.Fig. 10: Visualization of jet flow; Re = 34 .

Fig. 10: Visualization of jet flow; Re = 34 . From: E. Zauner [13] Fig. 3 The toroidal eddy and the finite length of the jet are clearly seen. From the data presented, the dimensionless jet length is calculated to be 86. This is somewhat less than the value shown in fig. 9 for Re = 32. A direct numerical solution of the integro-differential equation should be performed in order to see whether there is a better coincidence between theory and experiment. Usually, jet flow measurements are performed for the axial velocity at the jet axis. The investigations by Hamel and Richter [21] started at Re = 433. Rankin and co-workers [22] found a formula for the axial velocity which is proportional to Re; but the actual measurements started for Re = 1000, so it seems risky to extrapolate it for Re = 30. - A. J. Reynolds [18] describes observations of jet behavior for a range 16

of Reynolds numbers starting at about Re = 30. For Re < 150, the jet is steady for a certain length, roughly proportional to Re, above which the jet ”breaks up”. This may be interpreted as a fluctuating jet motion for which the mean velocity field may be calculated by solving eq. (6.1) described in this paper; but for this Reynolds number range it would be necessary to solve the integro-differential equation (9.9) together with (10.9).

12

Summary, and suggested further investigatons

A closed equation for the mean velocity field in fluctuating fluid flow has been derived employing ZwanzigMori technique of N-particle statistical mechanics. Arguments are presented that, for fluctuating flow, the statistical mechanics definition of fluid velocity is identical to the mean velocity of hydrodynamic turbulence theory; then, the Reynolds stress term is identical to the nonlinear part of the stress term in the mean velocity equation presented in the present paper. An analysis has been performed in order to bring this quantity into a form ready for numerical evaluation. It has been necessary to employ a power series approximation which requires that the theory in its present form is to be applied for small Reynolds numbers only. - The formula contains equilibrium time correlation functions up to fourth order. Multilinear mode coupling theory has been applied to calculate these quantities. Moreover, the formalism is restricted to shear motion (no sound and heat flow effects). The final result are the formulas (6.1), (6.6) to (6.11). Another major restriction of the investigation is that the N-particle theory in its present form is with no boundaries. This is not just a theoretical argument. The formulae have been applied to circular tube flow with the result that the calculated friction term is zero. This probably means that the theory in this form is unable to describe turbulent tube flow. As a test of the theory, circular jet flow is investigated. It became apparent that the theoretical formula is in contradiction with the jet flow formulae by Schlichting [12]. The theoretical formulae elaborated by Schneider and co-workers [14], [15], though in good agreement with the experiments of Zauner [13], could not be used in this test because of the point source construction for the inflow employed there. Thus, finally, an iteration procedure for the integro-differential equation derived in the present paper has been designed in order to discover for which Reynolds number the solutions start to deviate from the corresponding solutions of the Navier-Stokes equation. As is seen, noticeable deviations start for Re = 32; they are most distinct in the neighborhood of the axis and in some distance from the orifice. These details are compared with the results of an experimental investigation in the same Reynolds number range [13]. Some additional investigations would improve the understanding of the mean velocity equation and its solutions. An immediately arising problem is to calculate the laminar jet flow from the Navier-Stokes equation directly so that the solution returns the boundary condition at the orifice and the velocity profile at the jet axis correctly. That would perhaps also solve the second problem: Eliminate the scatter of the calculated values for the third-order term. - Since the iterative solution method employed in sec. 11 is for low Reynolds numbers only, it is another topic to develop a procedure to solve numerically the full integro-differential equation. Before the theoretical results described in this paper may possibly be applied to problems with higher Reynolds numbers, it is necessary to explore the analytical form of local equilibrium space-time correlations. It might be possible to proceed along the lines of the papers concerning total equilibrium correlations ([5],[6]). In addition, one needs a method to find the analytical form of total equilibrium correlations for flow configurations with boundaries. This is an important topic since there are many turbulent flow profiles well studied by experiment which belong to bounded flow configurations.

Appendices A

Derivation of mean velocity equation

In this appendix, we explain the definitions and show the main steps of the derivation of equation ( 1.1) with the formulas (2.1), (2.2) for the dissipative stress. The fluid is considered to be a system of N particles of mass m with positions y j and velocities v j which are combined to the phase space matrix z. Vector components are described by Latin indices, e. g. y j = {yja } . The system is enclosed in a box of 17

Volume V . Later on, the thermodynamic limit is performed. - A function g(z) is called a phase space function, or microscopic variable. Especially, we need the space densities of the conserved quantities mass, energy and momentum n, e, p ( [4] (8.2.4-6)) which are collected to a 5-element linear matrix a. They are functions of an additional space variable x: a=

N X j=1

e aj δ(x − y j )

(A.1)

ej = mv j , while the energy function contains For the particle functions e aj we have n ej = m , p the interparticle potential. m ist the particle mass. Quantities (like a) which are lists of 5 elements are denoted by normal letters while three-component vectors are bold. - The quantities a obey the conservation relations: · a = −∇ · s (A.2) The fluxes s have the same general structure as the a (2.1); especially, we have s1 = p . The time evolution of any phase space function A is described by the Liouville equation: ·

A = LA

(A.3)

L (often defined as iL) is the Liouville operator. From (A.3), the formal solution for A(t) given the initial value A is: A(t) = eLt A (A.4) eLt is an exponential operator. In the statistical model, z and N are considered random variables; that is, the probability density f (z, N ) is of grand canonical type. The ensemble mean value (expectation) of a phase space function A is defined in the ‘Heisenberg picture’: ∞ Z X hAi(t) = dz A(z, N, t)f (z, N ) (A.5) N =1

In this formula, f (z, N ) is the initial probability distribution, and A(z, N, t) is the value of A at time t if the initial positions and velocities of the particles are described by z . The operation (integration + Summation) is sometimes indicated by the symbol ‘tr’: tr{Ω} =

∞ Z X

dzΩ(z, N )

(A.6)

N =1

Certain probability densities (also called distributions here) are frequently used in the analysis. One of them is the (total) equilibrium distribution which corresponds to macroscopic rest: f0 = ψ(N ) exp(Φ0 + β(µN − H(z))) 1 m 3N ( ) ψ(N ) = N! h

(A.7a) (A.7b)

Here, h is Planck’s constant, β = 1/(kB T ) , kB being Boltzmann’s constant and T the temperature, µ is the chemical potential which is a function of mass density ρ = hni and temperature, and H(z) is Hamilton’s function which describes the total energy of the fluid. For the normalization constant, we have Φ0 = −βP V , P being the equilibrium pressure. Expectations with respect to the equilibrium distribution are denoted by hi0 . - In case of a simple fluid, the ’relevant probability distribution’ of Grabert’s formalism (see [4], sec. 2.2 and 8.3) is the local equilibrium distribution: fL (t) = ψ(N ) exp(Φ(t) − a(z) ◦ b(t)), 1 µ b = {β( u2 − ), β, −βu}, 2 m Φ(t) = − log(tr{ψ exp(−a ◦ b(t))}).

(A.8a) (A.8b) (A.8c)

Here the symbol ◦ is introduced for the operation: Multiplication, plus Summation over the 5 elements of the linear matrices a , b, plus Integration over geometrical space. The elements of b are called the conjugate parameters; they are functions of the quantities β , µ and u which we will sometimes call the thermodynamic parameters, and which will be considered to be slowly varying functions of space and 18

time. The b are defined such that the expectations of the a are identical to their expectations in local equilibrium: hai = haiL (A.9) The ZMT is a means for separating macroscopic and microscopic parts of the random variables. This is performed (for simple fluids) with the aid of a projection operator which projects out of any microscopic variable A the part which is proportional to the conserved quantities a. It reads: PA = hAiL + hA δaiL ◦ hδa δai−1 L ◦ δa

(A.10)

a−haiL ; hi−1 L denotes

Here, δa = the inverse of the expectation matrix in the formula. For general nonstationary flow, the local equilibrium distribution (A.8a), and therefore P, are time-dependent. (A.10) is the form of the projection operator which has been used in the present investigation. The definition in [4] uses an equivalent form (2.3.1). For eLt , Grabert uses the decomposition: eLt = eLt P(t) +

Z

0

t

0

˙ 0 ))(1 − P(t0 ))G(t0 , t) + (1 − P(t))G(0, t) dt0 eLt P(t0 )(L − P(t G(t0 , t) = exp− (

Z

t

t0

dt00 L(1 − P(t00 )))

(A.11) (A.12)

G(t0 , t) is called a time-ordered exponential Operator. The analysis in [4] consists in applying the Liouville equation (A.3) to a and inserting (A.4) and (A.11). Then, expectations in the Heisenberg picture (corresponding to (A.5)) of the equation are taken. It is postulated that the initial probability density is of the form of the ’relevant probability density’, which for simple fluids is the local equilibrium density (A.8a). Grabert states that this should not be considered a general restriction of the method but a means to form the general particle system into the type specially chosen (the simple fluid here); see [4], sec. 2.2 . It is shown that in this case the last term in (A.11) vanishes after averaging, while the first and the second term lead to the ’organized drift’ and ’disorganized drift’. Grabert’s final equation, the ’generalized transport equation’, can then be written, by using (A.2): ∂hai = −∇ · hsiL,t + D ∂t

(A.13)

This is a short form of [4], (8.1.13). As a last step, the momentum component of the equation is taken; the organized drift can be calculated further and results in the convolution and the pressure term of (1.1), while the disorganized drift is the dissipaative stress with formula (2.1).

B

Third-order term of the stress tensor

We need to calculate the secornd-order derivatives of the kernel function S. The detailed formula for the second derivative reads: δ 2 Sabcd (x0 , t0 ) δbe (x00 , t00 )δbf (x000 , t000 ) The first-order derivative consists of two parts called the parametric and the functional part:   δSabcd = −δ(t00 − t0 )β h[G(t0 , t)Q(t)sac (x)]Q(t0 )δa (x00 , t0 )Q(t0 )sbd (x0 )iL,t0 δb (x00 , t00 ) p

(B.1)

(B.2)

Q(t) = 1 − P(t)

(B.3)

δa (x, t) = a (x) − ha (x)iL,t

(B.4)



 δSabcd = βΘ(t00 − t0 )Θ(t − t00 )h[G(t0 , t00 )LP(t00 )δa (x00 , t00 )Q(t00 )G(t00 , t)Q(t)sac (x)]Q(t0 )sbd (x0 )iL,t0 δb (x00 , t00 ) f (B.5) 19

Θ is the unit step function: Θ(t) =



0, t < 0 1, t ≥ 0

(B.6)

For later reference, the first-order derivatives are formulated with respect to all five conjugated parameters (the greek index  is used). - These formulas contain operator chains which are, as usual, written like products. The question arises whether I can, when differentiating a chain, use the product rule; in other words, if A and B are operators, whether it is correct to write: δAB δA δB = B+A δb δb δb

(B.7)

Since I could not find a general proof in the literature, I checked (B.7) for any two consecutive operators appearing in (B.2), (B.5), and found it correct in all cases. The calculations are partly long and tedious and will not be shown here. - From (B.2), (B.5), the second-order functional derivatives are to be calculated. The formulas depend on b by 6 and 8 ’factors’, respectively; including in each formula the dependence of the local equilibrium probability density. The corresponding terms are given below. In the course of the calculation, for G(t0 , t) an identity formula is needed which can be checked by discretizing it: Q(t0 )G(t0 , t) = Q(t0 )G(t0 , t)Q(t) (B.8) Three of the terms are found to vanish; thus, the second-order derivative of S, and therefore the third-order term of R consist of 11 terms. For shortness, on the left hand side, indices and variables are omitted. Second order derivatives with respect to the momentum part of the conjugated parameters are needed only (that is, the greek index  is switched to e). For each term, the first step is just the definition of the quantity. δ(t) is the Dirac δ-function. (p,1) δfL (t0 ) δ2 S = −δ(t00 − t0 )β tr{ [G(t0 , t)Q(t)sac (x)][Q(t0 )sbd (x0 )]δae (x00 , t0 )} δb δb δbf (x000 , t000 )



= δ(t000 − t0 )δ(t00 − t0 )β h[G(t0 , t)ˆ sac (x, t)][Q(t0 )[Q(t0 )sbd (x0 )]δae (x00 , t0 )]δaf (x000 , t0 )iL,t0

(B.9)

(p,2) δ(δae (x00 , t0 )) δ2 S = −δ(t00 − t0 )β h[G(t0 , t)Q(t)sac (x)]Q(t0 )[Q(t0 )sbd (x0 )] iL,t0 δb δb δbf (x000 , t000 ) = −δ(t000 − t0 )δ(t00 − t0 )β h[ Q(t0 )sbd (x0 )Q(t)sac (x)(t0 , t)]iL,t0 hae (x00 )δaf (x000 , t0 )iL,t0 (B.10) 



(p,3) δ2 S δQ(t0 ) sbd (x0 )]δae (x00 , t0 )iL,t0 = −δ(t00 − t0 )β h[G(t0 , t)Q(t)sac (x)]Q(t0 )[ δb δb δbf (x000 , t000 )

= −δ(t000 − t0 )δ(t00 − t0 )×

× βh[G(t0 , t)Q(t)sac (x)]Q(t0 )[P(t0 )δaf (x000 , t0 )Q(t0 )sbd (x0 )]δae (x00 , t0 )iL,t0

(B.11)



(p,4) δ2 S δQ(t0 ) = −δ(t00 − t0 )β h[G(t0 , t)Q(t)sac (x)] [Q(t0 )sbd (x0 )]δae (x00 , t0 )iL,t0 δb δb δbf (x000 , t000 )

− δ(t000 − t0 )δ(t00 − t0 )×

× β h[G(t0 , t)Q(t)sac (x)]P(t0 )δaf (x000 , t0 )Q(t0 )[Q(t0 )sbd (x0 )]δae (x00 , t0 )iL,t0



(B.12)

(p,5) δ2 S δQ(t) = −δ(t00 − t0 )β h[G(t0 , t) sac (x)]Q(t0 )[Q(t0 )sbd (x0 )]δae (x00 , t0 )iL,t0 δb δb δbf (x000 , t000 ) = −δ(t000 − t)δ(t00 − t0 )×

β h[G(t0 , t)P(t)δaf (x000 , t)Q(t)sac (x)]Q(t0 )[Q(t0 )sbd (x0 )]δae (x00 , t0 )iL,t0 =0

20

(B.13)

The last step is correct because of (B.8) and QP = 0. 

(p,6) δ2 S δG(t0 , t) = −δ(t00 − t0 )β [ )sac (x)]Q(t0 )[Q(t0 )sbd (x0 )]δae (x00 , t0 )iL,t0 δb δb δbf (x000 , t000 )

= −δ(t00 − t0 )Θ(t000 − t0 )Θ(t − t000 )×

× β h[G(t0 , t000 )LP(t000 )δaf (x000 , t000 )Q(t000 )G(t000 , t)Q(t)sac (x)]Q(t0 )[Q(t0 )sbd (x0 )]δae (x00 , t0 )iL,t0 

δ2 S δb δb

(f,1)

(B.14)

= Θ(t00 − t0 )Θ(t − t00 ) ×

δfL (t0 ) [G(t0 , t00 )LP(t00 )δae (x00 , t00 )Q(t00 )G(t00 , t)Q(t)sac (x)]Q(t0 )sbd (x0 )} δbf (x000 , t000 ) = −δ(t000 − t0 )Θ(t00 − t0 )Θ(t − t00 ) × ×β tr{

×βh[G(t0 , t00 )LP(t00 )δae (x00 , t00 )Q(t00 )G(t00 , t)Q(t)sac (x)][Q(t0 )sbd (x0 )]δaf (x000 , t0 )iL,t0 (B.15)

(f,2) δ2 S = Θ(t00 − t0 )Θ(t − t00 ) × δb δb



×βh[G(t0 , t00 )LP(t00 )δae (x00 , t00 )Q(t00 )G(t00 , t)Q(t)sac (x)] =

δQ(t0 ) sbd (x0 )iL,t0 δbf (x000 , t000 )

δ(t000 − t0 )Θ(t00 − t0 )Θ(t − t00 ) ×

×βhG[(t0 , t00 )LP(t00 )δae (x00 , t00 )Q(t00 )G(t00 , t)Q(t)sac (x)]P(t0 )δaf (x000 , t0 )Q(t0 )sbd (x0 )iL,t(0 B.16) 

(f,3) δ2 S = Θ(t00 − t0 )Θ(t − t00 ) × δb δb

×βh[G(t0 , t00 )LP(t00 )δae (x00 , t00 )Q(t00 )G(t00 , t) =

δQ(t) sac (x)]Q(t0 )sbd (x0 )iL,t0 δbf (x000 , t000 )

δ(t000 − t)Θ(t00 − t0 )Θ(t − t00 ) ×

×βh[G(t0 , t00 )LP(t00 )δae (x00 , t00 )Q(t00 )G(t00 , t)P(t)δaf (x000 , t)Q(t)sac (x)]Q(t0 )sbd (x0 )iL,t0

= 0

(B.17)



(f,4) δ2 S = Θ(t00 − t0 )Θ(t − t00 ) × δb δb δQ(t00 ) ×βh[G(t0 , t00 )LP(t00 )δae (x00 , t00 ) G(t00 , t)Q(t)sac (x)]Q(t0 )sbd (x0 )iL,t0 δbf (x000 , t000 ) δ(t000 − t00 )Θ(t00 − t0 )Θ(t − t00 ) ×

=

0 ×βhG[(t0 , t00 )LP(t00 )δae (x00 , t00 )P(t00 )δaf (x000 , t00 )Q(t00 )G(t00 , t)Q(t)sac (x)]Q(t0 )sbd (x0 )iL,t(B.18)



δ2 S δb δb

(f,5)

= Θ(t00 − t0 )Θ(t − t00 ) ×

δ(δae (x00 , t00 )) Q(t00 )G(t00 , t)Q(t)sac (x)]Q(t0 )sbd (x0 )iL,t0 δbf (x000 , t000 ) = −Θ(t00 − t0 )Θ(t − t00 ) × δhae (x00 )iL,t00 ×βh[G(t0 , t00 )LP(t00 ) Q(t00 )G(t00 , t)Q(t)sac (x)]Q(t0 )sbd (x0 )iL,t0 δbf (x000 , t000 ) = δ(t000 − t00 )Θ(t00 − t0 )Θ(t − t00 ) × ×βh[G(t0 , t00 )LP(t00 )

=

×βh[G(t0 , t00 )LP(t00 )Q(t00 )G(t00 , t)Q(t)sac (x)]Q(t0 )sbd (x0 )iL,t0 hae (x00 )δaf (x000 , t00 )iL,t00

0

21

(B.19)

(f,6) δ2 S = Θ(t00 − t0 )Θ(t − t00 ) × δb δb δP(t00 ) ×βh[G(t0 , t00 )L δae (x00 , t00 )Q(t00 )G(t00 , t)Q(t)ˆ sac (x)]Q(t0 )ˆ sbd (x0 )iL,t0 δbf (x000 , t000 ) = −δ(t000 − t00 )Θ(t00 − t0 )Θ(t − t00 ) × 

×βh[G(t0 , t00 )LP(t00 )δaf (x000 , t00 )Q(t00 )δae (x00 , t00 )Q(t00 )G(t00 , t)Q(t)ˆ sac (x)] ×

× Q(t0 )ˆ sbd (x0 )iL,t0

(B.20)

(f,7) δ2 S = Θ(t00 − t0 )Θ(t − t00 ) × δb δb



δG(t00 , t) Q(t)sac (x)]Q(t0 )sbd (x0 )iL,t0 δbf (x000 , t000 ) Θ(t000 − t00 )Θ(t − t000 )Θ(t00 − t0 )Θ(t − t00 ) ×

×βh[G(t0 , t00 )LP(t00 )δae (x00 , t00 )Q(t00 ) =

×βh[G(t0 , t00 )LP(t00 )δae (x00 , t00 )Q(t00 ) ×

× G(t00 , t000 )LP(t000 )δaf (x000 , t000 )Q(t000 ))G(t000 , t)Q(t)sac (x)]Q(t0 )sbd (x0 )iL,t0(B.21)

(f,8) δ2 S = Θ(t00 − t0 )Θ(t − t00 ) × δb δb δG(t0 , t00 ) ×βh[ LP(t00 )δae (x00 , t00 )Q(t00 )G(t00 , t)Q(t)sac (x)]Q(t0 )sbd (x0 )iL,t0 δbf (x000 , t000 ) = Θ(t000 − t0 )Θ(t00 − t000 )Θ(t00 − t0 )Θ(t − t00 ) × 

×βh[G(t0 , t000 )LP(t000 )δaf (x000 , t000 )Q(t000 ))G(t000 , t00 ) × × LP(t00 )δae (x00 , t00 )Q(t00 )G(t00 , t)Q(t)sac (x)]Q(t0 )sbd (x0 )iL,t0

(B.22)

Results (B.17), (B.19) are found the same way as (B.13). - For application in (3.1), all terms have to be taken at b = b0 . This has effects listed below: - expectations in local equilibrium change to those in total equilibrium, denoted by h i0 ; - projection operators P(t), Q(t) change to time-independent operators P, Q; - δae (x, t) changes to ae (x); - G(t0 , t00 ) changes to the exponential operator g˜(t) with t = t00 − t0 : g˜(t) = eLQt

(B.23)

This operator has to be substituted by the similar operator g(t): g(t) = eQLt

(B.24)

The reason is that only g(t) can be substituted further to the non-projected exponential operator f (t) = eLt . - The substitution succeeds for all 11 terms with the aid of one of the following identities: Q g˜(t) g˜(t) L

= g(t) Q

= L g(t)

g˜(t)Q L = −PL + Lg(t)

(B.25a) (B.25b) (B.25c)

Next, we change from the microscopic variables a, s to the corresponding ortho-normalized quantities h, r; this extracts from each term the factor ( βρ )2 . Together with the factor β which appears explicitly in each factor, the pre-factor in (3.1) changes to 21 ρ2 β. The gradient operators appearing in (3.1) are transferred to the kernel function; by partial integration, this alters the sign of the expression. In the formula, we then have the expressions ∇c rac (x), ∇0d rbd (x0 ) which we denote ra (x), rb (x0 ) respectively; the distinction to the original r’s is by the number of indices appended. In addition, I used the abbreviation: rˆ = Qr 22

(B.26)

We obtain: (p,1) δ2 R = −δ(t000 − t0 )δ(t00 − t0 )h[Pr0 + Lg(t − t0 )h0 ][Qˆ r1∗ h∗2 ]h∗3 i0 δb δb 0  2 (p,2) δ R = −δ(t000 − t0 )δ(t00 − t0 ) h[g(t − t0 )ˆ r0 ]ˆ r1∗ i0 hh∗2 h∗3 i0 δb δb 0  2 (p,3) δ R = −δ(t000 − t0 )δ(t00 − t0 ) h[g(t − t0 )ˆ r0 ][P rˆ1∗ h∗3 ]h∗2 i0 δb δb 0





(p,4) δ2 R = δ(t000 − t0 )δ(t00 − t0 ) h[Pr0 + Lg(t − t0 )h0 ]Ph∗3 Qˆ r1∗ h∗2 i0 δb δb 0



(p,6) δ2 R = −δ(t00 − t0 )Θ(t000 − t0 )Θ(t − t000 )hg[(t000 − t0 )QLPh∗3 g(t − t000 )ˆ r0 ]ˆ r1∗ h∗2 i0 δb δb 0

(f,1) δ2 R = −δ(t000 − t0 )Θ(t00 − t0 )Θ(t − t00 )h[Lg(t00 − t0 )Ph∗2 g(t − t00 )ˆ r0 ]ˆ r1∗ h∗3 i0 δb δb 0  2 (f,2) δ R = δ(t000 − t0 )Θ(t00 − t0 )Θ(t − t00 )h[Lg(t00 − t0 )Ph∗2 g(t − t00 )ˆ r0 ]P rˆ1∗ h∗3 i0 δb δb 0  2 (f,4) δ R = δ(t000 − t00 )Θ(t00 − t0 )Θ(t − t00 )h[g(t00 − t0 )QLPh∗2 Ph∗3 g(t − t00 )ˆ r0 ]ˆ r1∗ i0 δb δb 0





δ2 R δb δb

(f,6) 0

= −δ(t000 − t00 )Θ(t00 − t0 )Θ(t − t00 )h[g(t00 − t0 )QLPh∗3 Qh∗2 g(t − t00 )ˆ r0 ]ˆ r1∗ i0

(f,7) δ2 R = Θ(t000 − t00 )Θ(t − t000 )Θ(t00 − t0 )Θ(t − t00 ) × δb δb 0





× h[g(t00 − t0 )QLPh∗2 g(t000 − t00 )QLPh∗3 g(t − t000 )ˆ r0 ]ˆ r1∗ i0

(f,8) δ2 R = Θ(t000 − t0 )Θ(t00 − t000 )Θ(t00 − t0 )Θ(t − t00 ) × δb δb 0

× h[g(t000 − t0 )QLPh∗3 g(t00 − t000 )QLPh∗2 g(t − t00 )ˆ r0 ]ˆ r1∗ i0

(B.27)

(B.28)

(B.29)

(B.30)

(B.31)

(B.32)

(B.33)

(B.34)

(B.35)

(B.36)

(B.37)

From now on all expectations are with respect to total equilibrium. We can therefore omit the 0 at the expectation symbol. - In wave number space, equilibrium correlation functions have a well-known property. The general form of such a function is: (1)

(2)

hB1 B2 · · · Bn(n) i

(B.38)

The upper indices just number the microscopic functions B (i) , while the lower ones are the number indices defined in sec. 3. The quantity contains a factor: (2π)3 δ(k1 + k2 + · · · + kn )

(B.39)

This fact simplifies some formulas considerably and will be used several times in the next sections. The kernel function of (3.1) has to be brought into a form so that it contains time correlation functions which can be evaluated by MCT. These are of the general form h[f (t)A]B ∗ i0 with f (t) = eLt . As a first step, the remaining projection operators P and Q are to be substituted. As can be seen in the formulas above, the terms contain 0 to 4 of these operators (the operator contained in rˆ remains unaltered). With the introduction of the normalized variables, the operators in wave number space read: P A = hAi + hA h∗4 ih4 23

(B.40)

Q A = A − hAi − hAh∗4 ih4

(B.41)

Generally, the component indices contained in the number indices in these formulas must be greek (running over the full list of 5 elements); it then would follow that I have to insert δh4 instead of h4 . But as will be seen later, finally all indices reduce to be ”latin”, i. e. will run over the 3 momentum components. Then, hh4 i = 0 , so that (B.40), (B.41) are correct. - The formulas are applied in several consecutive steps. After each step, certain auxiliary formulas are applied in order to simplify the expressions: hh1 rˆ2 i = h(Qh1 )r2 i = 0 hLAi = 0 hLg(t)hi = 0 hri = −hLhi = 0 Lh1 = −ˆ r1 − hr1 h4 ih4

(B.42)

The result for the 11 terms, now denoted K1 , · · · , K11 , is as follows. The δ- and θ-factors are temporarily omitted; when the time integrations start they will be added again: =

(K2 )0

=

(K3 )0

= −h[g(t − t1 )ˆ r0 ]h∗2 h∗4 ihˆ r1∗ h∗3 h4 i

(K4 )0 (K5 )0 (K6 )0 (K7 )0 (K8 )0

C

h[Lg(t − t1 )h0 ]h∗3 h∗4 ihˆ r1∗ h∗2 h4 i − h[Lg(t − t1 )h0 ]ˆ r1∗ h∗2 h∗3 i

(K1 )0

−h[g(t − t1 )ˆ r0 ]ˆ r1∗ ihh∗2 h∗3 i

= h[Lg(t − t1 )h0 ]h∗4 ihˆ r1∗ h∗2 h∗3 h4 i − h[Lg(t − t1 )h0 ]h∗4 ihˆ r1∗ h∗2 h5 ihh∗3 h4 h∗5 i = h[g(t − t3 )ˆ r0 ]h∗3 h∗5 ih[g(t3 − t1 )ˆ r5 ]ˆ r1∗ h∗2 i

= −h[g(t − t2 )ˆ r0 ]h∗2 h∗4 ih[Lg(t2 − t1 )h4 ]ˆ r1∗ h∗3 i

= h[g(t − t2 )ˆ r0 ]h∗2 h∗5 ih[Lg(t2 − t1 )h5 ]h∗4 ihˆ r1∗ h∗3 h4 i =

(K9 )0

=

(K10 )0

=

(K11 )0

=

(B.43)

−h[g(t − t2 )ˆ r0 ]h∗3 ih[g(t2 − t1 )ˆ r5 ]ˆ r1∗ ihh∗2 h∗5 i +h[g(t − t2 )ˆ r0 ]h∗3 h∗6 ih[g(t2 − t1 )ˆ r5 ]ˆ r1∗ ihh∗2 h∗5 h6 i −h[Lg(t − t2 )h0 ]h∗2 ih[g(t2 − t1 )ˆ r5 ]ˆ r1∗ ihh∗3 h∗5 i ∗ ∗ ∗ +h[g(t − t2 )ˆ r0 ]h2 h3 h5 ih[g(t2 − t1 )ˆ r5 ]ˆ r1∗ i +h[g(t − t2 )ˆ r0 ]h∗2 h∗6 ih[g(t2 − t1 )ˆ r5 ]ˆ r1∗ ihh∗3 h∗5 h6 i h[g(t − t3 )ˆ r0 ]h∗3 h∗7 ih[g(t2 − t1 )ˆ r5 ]ˆ r1∗ ih[g(t3 − t2 )ˆ r7 ]h∗2 h∗5 i ∗ ∗ ∗ h[g(t − t2 )ˆ r0 ]h2 h7 ih[g(t3 − t1 )ˆ r5 ]ˆ r1 ih[g(t2 − t3 )ˆ r7 ]h∗3 h∗5 i

Preparation of the kernel function

It is seen that the time correlations in the formulas (B.43) belong to one of two groups described by: h[g(t)ˆ r0 ]B ∗ i

∗

h[Lg(t)h0 ]B i

(C.1a) (C.1b)

B is a certain microscopic function. The formulas which are still built with the operator g(t) are to be transformed into those with f (t). This will be done in a number of steps. We start with (C.1a). First, we have: h[g(t)ˆ r0 ]h∗1 i = 0

(C.2)

(B.25a) and Qh1 = 0 have been used. Next, we introduce a symbol for the memory function of the process: Γ01 (t) + h[g(t)ˆ r0 ]ˆ r1∗ i

(C.3)

The memory function is, as usual, ”localized” in time, i. e. approximated by: Γ01 (t) = δ(t)(2π)3 δ(k0 − k1 )γ 01 24

(C.4)

γ 01 , when written with component indices a, b, for incompressible fluid reduces to: γ ab = νk02 δ ab

(C.5)

ν is the kinematic viscosity, and δ ab is the Kronecker symbol. - Next, we need an operator identity described in [8], Appendix B2; in the time-dependent form, and applied to the operators defined in this paper, it reads: Z t g(t) = f (t) − dt0 f (t0 )PLg(t − t0 ) (C.6) 0

We obtain the correlation function identity, for any microscopic function B: Z t h[g(t)ˆ r0 ]B ∗ i = h[f (t)ˆ r0 ]B ∗ i − dt0 Γ04 (t − t0 )h[f (t0 )h4 ]B ∗ i

(C.7)

0

After localization of the memory function, we get:

h[g(t)ˆ r0 ]B ∗ i = h[f (t)ˆ r0 ]B ∗ i − γ 04 h[f (t)h4 ]B ∗ i

(C.8)

The first term on the rhs of (C.7) is still to be transformed. The fifth formula (B.42) together with (A.3) and ω 04 = hr0 h∗4 i yields: rˆ0 = −∂t h0 − ω 04 h4 (C.9) We obtain: h[g(t)ˆ r0 ]B ∗ i = −∂t h[f (t)h0 ]B ∗ i − κ04 h[f (t)h4 ]B ∗ i

(C.10)

κ04 = ω 04 + γ 04

(C.11)

We denote with Cn die time correlation functions which can be calculated by MCT: (Cn )01···n−1 (t) + h[f (t)h0 ]h∗1 · · · h∗n−1 i

(C.12)

Then, we have the following applications of (C.7): h[g(t)ˆ r0 ]h∗1 h∗2 i = −∂t (C3 )012 (t) − κ08 (C3 )812 (t)

(C.13)

h[g(t)ˆ r0 ]h∗1 h∗2 h∗3 i = −∂t (C4 )0123 (t) − κ08 (C4 )8123 (t)

(C.14)

Somewhat more complicated is the transformation if B contains rˆ as a factor. Again, (C.9) is applied, together with ω ∗04 = −ω 04 : h[g(t)ˆ r0 ]ˆ r1∗ h∗2 i = −∂t h[f (t)h0 ]ˆ r1∗ h∗2 i − κ04 h[f (t)h4 ]ˆ r1∗ h∗2 i

= ∂t h[f (t)h0 ](∂t h∗1 )h∗2 i − ω 17 ∂t h[f (t)h0 ]h∗7 h∗2 i + κ04 (h[f (t)h4 ](∂t h∗1 )h∗2 i − ω 17 h[f (t)h4 ]h∗7 h∗2 i)

(C.15)

The first and the third term need further transformation. Because of the symmetries of (3.2), the following formula is correct if it is used only within this expression: (∂t h∗1 )h∗2 = 21 ((∂t h∗1 )h∗2 + h∗1 ∂t h∗2 ) = 12 ∂t (h∗1 h∗2 )

(C.16)

We obtain, for the correlation function (C.15): h[g(t)ˆ r0 ]ˆ r1∗ h∗2 i =

= − 12 ∂tt (C3 )012 (t) − ω 17 ∂t (C3 )072 (t) − κ04 ( 12 ∂t (C3 )412 (t) + ω 17 (C3 )472 (t))

(C.17)

We now switch to time correlations belonging to group (C.1b). We need an operator identity similar to (C.6) which again can be found in [8], Appendix 2B. Applied to this paper, it reads: Z t g(t) = f (t) − dt0 g(t0 )PLf (t − t0 ) (C.18) 0

25

We obtain a formula similar to (C.7): h[Lg(t)h0 ]B ∗ i = h[Lf (t)h0 ]B ∗ i −

Z

t

0

dt0 h[Lf (t − t0 )h0 ]h∗4 ih[Lg(t0 )h4 ]B ∗ i

(C.19)

For shortness, we introduce some new denotations: Ψ0 (t) = h[Lg(t)h0 ]B ∗ i ∗

(C.20) ∗

Φ0 (t) = h[Lf (t)h0 ]B i = ∂t h[f (t)h0 ]B i

(C.21)

Then, applying (C.9), (C.19) reads: Ψ0 (t) = Φ0 (t) −

Z

0

t

dt0 ∂t (C2 )04 (t − t0 )Ψ4 (t0 )

(C.22)

This is an integral equation for Ψ which can be formally solved by temporarily switching to Laplace space. The result is: Z t Ψ0 (t) = Φ0 (t) + κ04 dt0 Φ4 (t0 ) (C.23) 0

Taking (C.21) into account, we can integrate (C.19):

h[Lg(t)h0 ]B ∗ i = ∂t h[f (t)h0 ]B ∗ i + κ04 (h[f (t)h4 ]B ∗ i − hh4 B ∗ i)

(C.24)

This is the basic formula for time correlations of form (C.1b). Taking B = h1 , we have: h[Lg(t)h0 ]h∗1 i = ∂t (C2 )01 (t) + κ07 ((C2 )71 (t) − δ 71 )

(C.25)

The lowest-order non-projected time correlation function C2 obeys a simple differential equation ((D.12), with κ constant in t): ∂t (C2 )01 (t) = −κ07 (C2 )71 (t)

(C.26)

h[Lg(t)h0 ]h∗1 i = −κ01

(C.27)

This is introduced into (C.25):

The lowest-order time correlation of type (C.1b) is a constant. - For the next-higher order, we find from (C.24): h[Lg(t)h0 ]h∗1 h∗2 i = ∂t (C3 )012 (t) + κ07 ((C3 )712 (t) − (j3 )712 ) (j3 )012 = (C3 )012 (t = 0) = hh0 h∗1 h∗2 i

(C.28) (C.29)

Again the situation is somewhat more complicated if B contains the projected flux rˆ. For the lowest order, we obtain, by applying (C.24): h[Lg(t)h0 ]ˆ r1∗ h∗2 i = ∂t h[f (t)h0 ]ˆ r1∗ h∗2 i + κ04 (h[f (t)h4 ]ˆ r1∗ h∗2 i − hh4 rˆ1∗ h∗2 i)

(C.30)

We apply (C.9): h[Lg(t)h0 ]ˆ r1∗ h∗2 i =

−∂t h[f (t)h0 ](∂t h∗1 )h∗2 i + ω 17 ∂t h[f (t)h0 ]h∗7 h∗2 i +

+κ04 (−h[f (t)h4 ](∂t h∗1 )h∗2 i + ω 17 h[f (t)h4 ]h∗7 h∗2 i − hh4 rˆ1∗ h∗2 i)

(C.31)

We introduce a notation for the static correlation appearing in the formula: (s3 )124 = hˆ r1 h2 h∗4 i = −hˆ r1 h2 h∗4 i∗ = −hh4 rˆ1∗ h∗2 i

(C.32)

Finally, we can apply a symmetry consideration as in (C.15); by using (C.16) we obtain: h[Lg(t)h0 ]ˆ r1∗ h∗2 i =

1 2 ∂tt (C3 )012 (t) + ω 17 ∂t (C3 )072 (t) + +κ04 ( 12 ∂t (C3 )412 (t) + ω 17 (C3 )472 (t)

26

+ (s3 )124 )

(C.33)

For the next-higher order correlation we have, by (C.24): h[Lg(t)h0 ]ˆ r1∗ h∗2 h∗3 i = ∂t h[f (t)h0 ]ˆ r1∗ h∗2 h∗3 i + κ04 (h[f (t)h4 ]ˆ r1∗ h∗2 h∗3 i − hh4 rˆ1∗ h∗2 h∗3 i)

(C.34)

Introduction of (C.9) yields: h[Lg(t)h0 ]ˆ r1∗ h∗2 h∗3 i =

−∂t h[f (t)h0 ](∂t h∗1 )h∗2 h∗3 i + ω 17 ∂t h[f (t)h0 ]h∗7 h∗2 h∗3 i)

+κ04 (−h[f (t)h4 ](∂t h∗1 )h∗2 h∗3 i + ω 17 h[f (t)h4 ]h∗7 h∗2 h∗3 i − hh4 rˆ1∗ h∗2 h∗3 i) (C.35)

The symmetry consideration which is correct if it is applied within the integral formula (3.2) gives: (∂t h∗1 )h∗2 h∗3 = 31 ((∂t h∗1 )h∗2 h∗3 + h∗1 (∂t h∗2 )h∗3 + h∗1 h∗2 ∂t h∗3 ) = 31 ∂t (h∗1 h∗2 h∗3 )

(C.36)

Again a denotation for the static correlation is introduced: (s4 )1234 = hˆ r1 h2 h3 h∗4 i = −hˆ r1 h2 h3 h∗4 i∗ = −hh4 rˆ1∗ h∗2 h∗3 i

(C.37)

The final form for (C.33) is: h[Lg(t)h0 ]ˆ r1∗ h∗2 h∗3 i =

1 3 ∂tt (C4 )0123 (t) + ω 17 ∂t (C4 )0723 (t) + +κ04 ( 31 ∂t (C4 )4123 (t) + ω 17 (C4 )4723 (t)

+ (s4 )1234 )

(C.38)

The formulas obtained in this section are introduced into the terms (B.43). The Dirac and Heaviside functions omitted are added again. The kernel function K in (3.2) now consists of 5 main parts: K

=

δ(t000 − t0 )δ(t00 − t0 )M1

+δ(t00 − t0 )Θ(t000 − t0 )Θ(t − t000 )M2

+δ(t000 − t00 )Θ(t00 − t0 )Θ(t − t00 )δ(t00 − t0 )M3 000

00

000

00

0

00

(C.39) 00

0

+Θ(t − t )Θ(t − t )Θ(t − t )Θ(t − t )δ(t − t )M4

+Θ(t000 − t0 )Θ(t00 − t000 )Θ(t00 − t0 )Θ(t − t00 )δ(t000 − t0 )M5

We denote: (j2 )12 = hh∗1 h∗2 i

(C.40)

The coefficients M1 , · · · , M5 are: (M1 )0123

=

−(j2 )23 γ 01 δ(t − t0 )

−2(s3 )124 (κ05 (C3 )534 (t − t0 ) + ∂t (C3 )034 (t − t0 ))

(M2 )0123

=

−ω 18 ((C4 )7823 (t − t0 )κ07 + ∂t (C4 )0823 (t − t0 )) 1 − (κ07 ∂t (C4 )7123 (t − t0 ) + ∂tt (C4 )0123 (t − t0 )) 3 ((C3 )824 (t − t00 )κ08 + ∂t (C3 )024 (t − t00 ))

×(κ45 ((s3 )513 − (s3 )135 )

(M3 )0123

+2ω 16 ((C3 )763 (t00 − t0 )κ47 + ∂t (C3 )463 (t00 − t0 ) 1 +κ47 ∂t (C3 )713 (t00 − t0 ) + ∂tt (C3 )413 (t00 − t0 )) 2 = γ 51 (κ08 ((C3 )836 (t − t00 )(j3 )625 + (C3 )826 (t − t00 )(j3 )635 ) +∂t (C3 )026 (t − t00 )(j3 )635 + ∂t (C3 )036 (t − t00 )(j3 )625

(M4 )0123 (M5 )0123

−(C4 )8235 (t − t00 )κ08 − ∂t (C4 )0235 (t − t00 ))

= γ 51 (κ08 (C3 )837 (t − t000 ) + ∂t (C3 )037 (t − t000 )) ×(κ79 (C3 )925 (t000 − t00 ) + ∂t (C3 )725 (t000 − t00 )

= γ 51 (κ08 (C3 )827 (t − t00 ) + ∂t (C3 )027 (t − t00 ))

×(κ79 (C3 )935 (t00 − t000 ) + ∂t (C3 )735 (t00 − t000 ) 27

(C.41)

When the M -quantities are inserted into (3.2) and the approximations described in sec. 3 are performed, finally two parts of the kernel function remain: Φ10123 (t0 )

= −(j2 )23 γ 01 δ(t0 )

−2(s3 )124 (κ05 (C3 )534 (t0 ) + ∂t (C3 )034 (t0 ))

+2(j3 )625 γ 51 ((κ08 ((C3 )836 (t0 ) + ∂t (C3 )036 (t0 ))

Φ20123 (t0 , t00 )

−κ51 ((C4 )8235 (t0 )κ08 + ∂t (C4 )0235 (t0 )) 1 − (κ07 ∂t (C4 )7123 (t0 ) + ∂tt (C4 )0123 (t0 )) 3 = 2γ 51 (κ08 (C3 )837 (t0 ) + ∂t (C3 )037 (t0 )) 00

(C.42)

00

×(κ79 (C3 )925 (t ) + ∂t (C3 )725 (t ))

+((C3 )824 (t0 )κ08 + ∂t (C3 )024 (t0 )) ×(κ45 ((s3 )513 − (s3 )135 )

+2ω 16 ((C3 )763 (t00 )κ47 + ∂t (C3 )463 (t00 )) 1 +κ47 ∂t (C3 )713 (t00 ) + ∂tt (C3 )413 (t00 )) 2 In the relevant part of the second-order term of D (3.4), the kernel function K2 is the sum of (B.2) and (B.5), taken at b = b0 for  = 1; (B.25a) has been used to substitute the operator g˜: (K2 )012 (t, t0 , t00 ) = −δ(t00 − t0 ) h[g(t − t0 )ˆ r0 ]ˆ r1∗ δh∗2 i+

+ Θ(t00 − t0 )Θ(t − t00 )h[g(t00 − t0 )QLPδh∗2 g(t − t00 )ˆ r0 ]ˆ r1∗ i

(C.43)

For the first term, we have: h[g(t − t0 )ˆ r0 ]ˆ r1∗ δh∗2 i = h[g(t − t0 )ˆ r0 ]ˆ r1∗ h∗2 i − h[g(t − t0 )ˆ r0 ]ˆ r1∗ ihh∗2 i

(C.44)

The operator P in the second term in (C.43) is evaluated: h[g(t00 − t0 )QLPδh∗2 g(t − t00 )ˆ r0 ]ˆ r1∗ i = −h[g(t − t00 )ˆ r0 ]h∗2 h∗3 ih[g(t00 − t0 )ˆ r3 ]ˆ r1∗ i

(C.45)

The identity QLδh3 = −Qr3 = −ˆ r3 has been used. Originally, the full form (B.40) with δh instead of h had to be used in the first factor; but it is found that the terms with the constant parts hh2 i, hh3 i vanish, so that the final form is (C.45). For further preparation the formulas (C.3) and (C.4), (C.13), (C.17) are used. In (3.4) integration over t00 and an approximation because of different time scales is performed. This leads to formula (3.5) for the second-order term, with the kernel function K2 now depending only on the time variable t0 : (K2 )012 (t0 )

= γ 01 δ(t0 )hh2 i

+γ 31 (∂t0 (C3 )032 (t0 ) + κ04 (C3 )432 (t0 )) +∂t0 t0 (C3 )012 (t0 ) + i kc00 ∂t0 (C3 )012 | 0

+κ03 (∂t0 (C3 )312 (t ) + i

0 0 0 20 =c (t ) + ω 14 ∂t0 (C3 )042 (t ) 00 0 0 kc (C3 )312 | 0 20 =c (t ) + ω 14 (C3 )342 (t ))

(C.46)

The index 0 20 =c means that the component index of the number index 2 runs over the three values c = 3, 4, 5; in contrast to the general prescription for the second-order term where it has the fixed value 1.

D

Multilinear mode-coupling theory

We refer to the paper by van Zon and Schofield [6]. We continue to use the ortho-normalized microscopic variables h defined in sec. 3. In [6], multiple variables Aα are employed which are products of the single variables: Aα = h1 h2 · · · hα (D.1) 28

Greek indices, in addition to describing a component index running over 5 elements, are used in this appendix in connection with the multiple variables. Depending on the context, they denote either the order of the multiple variable or the full set of number indices: α =⇒ {1, 2, · · · , α}

(D.2)

Orthogonal multiple variables Qα and projection operators Pα are successively defined (do not mix up these symbols with the (non-indexed) P and Q used in earlier sections):  Q0 = 1   Pα−1  Qα = (1 − ν=0 Pν )Aα (D.3) P0 X = hXi    Pα X = hX Qα iQα

A greek symbol appearing twice includes summation over each of the component indices and integration over each of the wave numbers; for each wave number including a factor (2π)−3 . The last of the formulas (D.3) is somewhat simpler than [6] (6), owing to the fact that our basic variables h are normalized. - An important property of the theory is that it is formulated before performing the thermodynamic limit; that is, there is a finite number of particles, N . Finally, the limit N → ∞ is taken. In the formulas appear terms which are proportional to different powers of N ; in the limit, only those with the highest power of N ’survive’. In this paper, we will provide only the final formulas obtained after performing the thermodynamic limit. - The object of the theory are time correlations of the Q: Gαβ = hQα (t)Q∗β i

(D.4)

Again, formula (D.4) is somewhat simpler than [6] (15). In the correlations appearing in D3 and D2 , we have α = 1; for β = 1, (D.4) is directly the two-point correlation C2 , while the higher-order correlations C3 and C4 are obtained from (D.4) for β = 2 and 3, respectively. In the latter cases, formula [6] (26) (given in Laplace space) has to be applied. Since β 6= α, the first term vanishes. The case β = 2 is analyzed in [6] as an example. The N -ordering analysis shows that only the second term of the formula survives. For β = 3, the second and certain parts of the third term survive. But the detailed analysis shows that the third term contains a dimensionless combination of thermodynamic parameters which is  1, so that finally both cases are described by the same formula: Gαβ (ζ) = Gαα0 (ζ)Mα0 β 0 (ζ)Gβ 0 β (ζ)

(D.5)

Mαβ (ζ) = hQ˙ α Q∗β i − Γαβ (ζ)

(D.6)

Certain simplifications (compared with [6] (26)) stem from the fact that for the quantities considered β 6= α. ζ is the independent variable of Laplace space. α0 describes an index set of the order α with the same wave numbers as α but different component indices. Γαβ is the time correlation of the fluctuating forces defined in [6] (11). These quantities describe the friction effects of the fluid. Generally, Γαβ (z) is replaced by its value at z = 0 (’localization’ in time). The lowest order quantity, α = β = 1, for small wave numbers resembles the Green-Kubo expressions. In [6], there is no prescription for the higher order Γαβ . We refer to the older correlation function theory by Kawasaki [9], where these quantities are neglected. - In addition to (D.5, (D.6), it is shown ( [6] (25)) that, after transforming back to the time domain, for any α, in highest N -order Gαα0 (t) factorizes into the product of the corresponding two-point correlations. - Essentially, (D.5), (D.6) express the higher order correlations in terms of C2 . The constant matrices in the formulas for C4 (4.4), (4.5) read: Φ12345 = (j5 )12345 − (j3 )123 (j2 )45 − (j4 )1237 (j3 )745

(D.7)

(σ 4 )1234 = (s4 )1234 − Φ23456 (s3 )156  (j2 )12 = hh1 h2 i     (j3 )123 = hh1 h2 h3 i  ∗ ∗ ∗ (j4 )1234 = hh1 h2 h3 h4 i  (j5 )12345 = hh1 h2 h∗3 h∗4 h∗5 i     (s4 )1234 = hˆ r1 h∗2 h∗3 h∗4 i

(D.8)

29

(D.9)

The choice of conjugate complex variables in (D.9) looks somewhat arbitrary; the quantities are defined as they appear in (4.4) and (D.7), (D.8). j3 is different from (C.29) by the δ-factor, but again we use the same denotation. - For the applications in this paper, it is necessary to investigate the equation for C2 (t) (G11 in [6], called the propagator) in some detail. Formula [6] (34) is approximated to lowest order, which in the time domain leads to the equation: ∂t (C2 )12 (t) = −ω 13 (C2 )32 (t) −

Z

0

t

dt0 (Γ11 )13 (t − t0 )(C2 )32 (t0 )

ω 12 = −hh˙ 1 h∗2 i

(D.10)

(D.11)

Γ11 is localized in time, as in (C.4). But then ∂t C2 is discontinuous at t = 0 since ∂t (C2 )12 (0) = −ω 12 ; therefore, the formula resulting from (D.10 ) must be written: ∂t (C2 )12 (t) = −κ13 (t)(C2 )32 (t) κ13 (t) =  ω 13 + θ(t)γ 13 1, t > 0 θ(t) = 0, t ≤ 0

   

(D.12)

  

The time dependence of κ13 is important when one needs to calculate the second derivative of C2 ; otherwise, it can be neglected. We have: ∂t κ13 (t) = δ(t)γ 13

(D.13)

∂tt (C2 )12 (t) = −δ(t)γ 13 (j2 )32 + κ13 κ34 (C2 )42 (t)

(D.14)

In (D.14) the denotation (C2 )12 (0) = hh1 h∗2 i = (j2 )12 has been used.

E

Evaluation of the correlation formulas

The formula for D 3 , (3.6), is in explicit notation:

1 1 (D3 )0 (k, t) = − ρ2 β( 2 (2π)6

Z

dk0 dk00 dk000 u1 (k0 , t)u2 (k00 , t)u3 (k000 , t) δ(k − k0 − k00 − k000 ) × . × (∆3 )0123 (k0 , k00 , k000 ) 1

+ ξ− 2

(∆3 )0123 (k0 , k00 , k000 ) =

Z

∞

1 (2π)3

Z

 dk0 dk00 u1 (k0 ,t)uq(k00 t) δ(k − k0 − k00 )(∆2 )01 (k0 , k00 )

dt0 Φ10123 (k0 , k00 , k000 , t0 ) + lim

t→∞

0

(∆2 )01 (k0 , k00 ) =

Z

∞

Z

0

t

dt0

Z

t−t0

0

dt0 (K2 )01 k0 , k00 ,t

0

ξ = ρβ∂ρ P |β

(E.1)

dt00 Φ20123 (k0 , k00 , k000 , t0 , t00 ) (E.2)


0

(E.3) (E.4)

The factor (2π)3 δ(k − k0 − k00 − k000 ) has been extracted from the quantities Φ1, Φ2 in (C.42); also, (2π)3 δ(k − k0 − k00 ) has been extracted from K2 in (C.46). uq(k) is the Fourier transform of (u(x))2 . 30

These are the resulting formulas for the integrands in (E.2), (E.3): Φ10123 (k0 , k00 , k000 , t0 ) = −(2π)3 δ 23 δ(k − k0 )γ 01 (k)δ(t0 ) − 2s124 (k0 )(κ05 (k)(C3 )534 (k, k000 , k0 + k00 , t0 ) + ∂t (C3 )034 ((k, k000 , k0 + k00 , t0 )) + 2j625 γ 51 (k0 )(κ05 (k)(C3 )536 (k, k000 , k0 + k00 , t0 ) + ∂t (C3 )036 (k, k000 , k0 + k00 , t0 )) − κ51 (k0 )(κ06 (k)(C4 )6235 (k, k00 , k000 , k0 ,t0 ) + ∂t (C4 )0235 (k, k00 , k000 , k0 ,t0 )) 1 − (κ05 (k)∂t (C4 )5123 (k, k0 , k00 , k000 ,t0 ) + ∂tt (C4 )0123 ((k, k0 , k00 , k000 ,t0 )) 3

(E.5)

Φ20123 (k0 , k00 , k000 , t0 , t00 ) = = 2γ 51 (k0 )(κ06 (k)(C3 )637 (k, k000 , k0 + k00 , t0 ) + ∂t (C3 )037 (k, k000 , k0 + k00 , t0 ))× × (κ78 (k0 + k00 )(C3 )825 (k0 + k00 , k00 , k0 ,t00 ) + ∂t (C3 )725 (k0 + k00 , k00 , k0 ,t00 ) + 2(κ05 (k)(C3 )524 (k, k00 , k0 + k000 , t0 ) + ∂t (C3 )024 (k, k00 , k0 + k000 , t0 ))×

× (ω 16 (k0 )(κ47 (k0 + k000 )(C3 )763 (k0 + k000 , k0 , k000 ,t00 ) + ∂t (C3 )463 (k0 + k000 , k0 , k000 ,t00 )) 1 + κ47 (k0 + k000 )∂t (C3 )713 (k0 + k000 , k0 , k000 , t00 ) + ∂tt (C3 )413 (k0 + k000 , k0 , k000 ,t00 )) 2

(E.6)


(K2 )01 k0 , k00 , t0 = γ 01 (k)(2π)3 δ(k00 )hh2 iδ(t0 )

+ γ 31 (k0 )(∂t0 (C3 )032 (k, k0 , k00 , t0 ) + κ04 (k)(C3 )432 (k, k0 , k00 , t0 )) + ∂t0 t0 (C3 )012 (k, k0 , k00 , t0 ) + i kc00 ∂t0 (C3 )01c (k, k0 , k00 , t0 ) + ω 14 (k0 )∂t0 (C3 )042 (k, k0 , k00 , t0 ) + κ03 (k)(∂t0 (C3 )312 (k, k0 , k00 , t0 ) + i kc00 (C3 )31c (k, k0 , k00 , t0 ) + ω 14 (k0 )(C3 )342 (k, k0 , k00 , t0 ))

(E.7)

The wave number k is, in (E.5) and (E.6), to be substituted by k0 + k00 + k000 , and in (E.7) by k0 + k00 . - The equation for C2 (D.10 ), (D.12) now reads: ∂t (C2 )12 (k,t) = −κ13 (k,t)(C2 )32 (k,t)

(E.8)

We have, after factor extraction, (j2 )12 = δ 12 . Therefore, for the second derivative, one obtains from (D.14): ∂tt (C2 )12 (k,t) = −δ(t)γ 12 (k) + κ13 (k)κ34 (k)(C2 )42 (k,t)

(E.9)

The detailed form of (4.1), (4.2) is: (C3 )123 (k, k0 , k00 ,t) = (C2 )14 (k,t)(j3 )423 + (G12 )123 (k, k0 , k00 ,t) Z t 0 00 (G12 )123 (k, k , k ,t) = − dt0 (C2 )14 (k,t − t0 )(s3 )456 (k)(C2 )25 (k0 ,t0 )(C2 )36 (k00 ,t0 )

(E.10) (E.11)

0

We write C3 and G12 as functions of three wave numbers which correspond to the wave numbers of the three h-factors in the definition formula of C3 ; though, because of k00 = k − k0 , they are actually functions of k, k0 only. - (4.4) , (4.5) turn to: (C4 )1234 (k, k0 , k00 , k000 ,t) = ((C3 )156 (k, k0 , k00 ,t) − (j3 )156 )Φ23456 +

+ (C2 )17 (k,t)(j4 )7234 + (G13 )1234 (k, k0 , k00 , k000 ,t)

(G13 )1234 ((k, k0 , k00 , k000 ,t) = Z t =− dt0 (C2 )15 (k,t − t0 )(σ 4 )5678 (k)(C2 )26 (k0 ,t0 )(C2 )37 (k00 ,t0 )(C2 )48 ((k000 ,t0 ) 0

31

(E.12)

(E.13)

with k = k0 + k00 + k000 . Time derivatives of C3 , C4 will be reformulated. From (4.1) we have: ∂t (C3 )123 (k, k0 , k00 ,t) = −κ15 (k,t)(C2 )54 (k,t)(j3 )423 + ∂t (G12 )123 (k, k0 , k00 ,t)

(E.14)

∂tt (C3 )123 (k, k0 , k00 ,t) = = (−δ(t)γ 15 (k) + κ16 (k)κ65 (k)(C2 )54 (k,t))(j3 )423 + ∂tt (G12 )123 (k, k0 , k00 ,t)

(E.15)

From (4.2): ∂t (G12 )123 (k, k0 , k00 ,t) = −(s3 )156 (k)(C2 )25 (k0 ,t)(C2 )36 (k00 ,t)+ Z t dt0 κ17 (k,t − t0 )(C2 )74 (k,t − t0 )(s3 )456 (k)(C2 )25 (k0 ,t0 )(C2 )36 (k00 ,t0 ) + 0

= −(s3 )156 (k)(C2 )25 (k0 ,t)(C2 )36 (k00 ,t) − κ14 (k)(G12 )423 (k, k0 , k00 ,t)

(E.16)

We remember that k00 = k − k0 . For ∂tt G12 we obtain, after some manipulations: ∂tt (G12 )123 (k, k0 , k00 ,t)

=

(ζ 3 )145678 (k, k0 , k00 )(s3 )645 (k)(C2 )27 (k0 ,t)(C2 )38 (k00 ,t) + + κ14 (k)κ45 (k)(G12 )523 (k, k0 , k00 ,t)

(ζ 3 )123456 (k, k0 , k00 ) = κ14 (k)δ 25 δ 36 + δ 14 κ25 (k0 )δ 36 + δ 14 δ 25 κ36 (k00 )

(E.17) (E.18)

By a similar calculation, we obtain for C4 , from (4.4): ∂t (C4 )1234 (k, k0 , k00 ,t) = ∂t (C3 )156 (k, k0 ,t)Φ23456 + + ∂t (C2 )17 (k,t)(j4 )7234 + ∂t (G13 )1234 (k, k0 , k00 ,t)

(E.19)

∂tt (C4 )1234 (k, k0 , k00 ,t) = ∂tt (C3 )156 (k, k0 ,t)Φ23456 + + ∂tt (C2 )17 (k,t)(j4 )7234 + ∂tt (G13 )1234 (k, k0 , k00 ,t)

(E.20)

For the time derivatives of C2 and C3 , the foregoing formulas can be used for substitution. From (E.13), we find: ∂t (G13 )1234 ((k, k0 , k00 , k000 , t) = = −(σ 4 )1567 (k)(C2 )25 (k0 ,t)(C2 )36 (k00 ,t)(C2 )47 ((k000 ,t) − κ15 (k)(G13 )5234 ((k, k0 , k00 , k000 , t)

(E.21)

∂tt (G13 )1234 ((k, k0 , k00 , k000 , t) = = (ζ 4 )15678¯2¯3¯4 (k, k0 , k00 , k000 )(σ 4 )8567 (k)(C2 )2¯2 (k0 ,t)(C2 )3¯3 (k00 ,t)(C2 )4¯4 ((k000 ,t)+ + κ15 (k)κ56 (k)(G13 )6234 ((k, k0 , k00 , k000 , t)

(E.22)

(ζ 4 )12345678 (k, k0 , k00 , k000 ) = = κ15 (k)δ 26 δ 37 δ 48 + δ 15 κ26 (k0 )δ 37 δ 48 + δ 15 δ 26 κ37 (k00 )δ 48 + δ 15 δ 26 δ 37 κ48 (k000 )

(E.23)

We introduce these results into (E.2). For C2 , the matrix exponential function (4.7) explicitly reads:. (C2 )12 (k,t) = e−κ12 (k)t

(E.24)

The time integrations in (E.2) result in inverses of ζ 3 , ζ 4 and of: (ζ 2 )1234 (k, k0 ) = κ13 (k)δ 24 + δ 13 κ24 (k0 )

(E.25)

(ζˇ3 )123456 (k, k0 , k00 ) = −κ14 (k)δ 25 δ 36 + δ 14 κ25 (k0 )δ 36 + δ 14 δ 25 κ36 (k00 )

(E.26)

32

For instance, we have: Z

0

∞

0 dt(C2 )12 (k,t)(C2 )34 (k0 ,t) = (ζ 2 )−1 1324 (k, k )

(E.27)

After integration, the parts from Φ1 and Φ2 can be united. Some terms cancel, others can be simplified. For instance, the relation applies: 0 00 0 00 −1 (ζ 3 )123478 (k, k0 , k00 )(ζ 2 )−1 7856 (k , k ) = κ14 (k)(ζ 2 )2356 (k , k ) + δ 14 δ 25 δ 36

(E.28)

The calculation of the kernel function ∆3 (E.2) yields the expression (again k = k0 + k00 + k000 ): ∆3 (k, k0 , k00 , k000 ) = − 31 (s4 )0123 (k)

+ γ 04 (k)( 13 (j4 )4123 − (2π)3 δ(k − k0 )δ 14 δ 23 )

00 0 000 + 2(s3 )056 (k)(s134 (k0 ) − (s3 )413 (k0 + k000 ))(ζ 2 )−1 5624 (k , k + k )

00 0 000 − 2(s3 )056 (k)(j734 γ 41 (k0 ) − 41 (j3 )413 γ 47 (k0 + k000 ))(ζ 2 )−1 5627 (k , k + k )

000 0 00 00 0 0 −1 + 2(s3 )047 (k)(s3 )856 (k0 + k00 )(ζ 2 )−1 4738 (k , k + k )(ζ 2 )5629 (k , k )κ19 (k ) 1 000 0 00 00 0 0 00 −1 + (s3 )047 (k)(s3 )589 (k0 + k00 )(ζ 2 )−1 4736 (k , k + k )(ζ 2 )8921 (k , k )κ56 (k + k ) 2   Z 1 ˘ ˘ ˘ ˘ ˘ + − dk Φ47856 (k, −k)(s3 )056 (k, −k) + (s4 )0478 (k) × (2π)3

00 000 0 0 × (ζ 3 )−1 478239 (k , k , k )κ19 (k )

(k00 , k0 + k000 )× − 21 (s3 )047 (k)(s3 )568 (k0 + k000 )(ζ 2 )−1 472¯ 9 0 000 0 000 0 000 × (ζˇ3 )−1 99 (k + k ) 913568 (k + k , k , k )κ¯  Z 1 −1 ˘ Φ23956 (k, ˘ −k)(ζ ˘ ˘ ˘ + (s3 )047 (k) κ19 (k0 ) dk 2 )4756 (k, −k)+ (2π)3

00 0 000 0 000 + 21 (j3 )813 (ζ 2 )−1 4729 (k , k + k )κ98 (k + k )




0 000 −1 + 21 (s3 )047 (k)(s3 )¯5¯6¯7 (k0 + k000 )(ζ 2 )−1 (k00 , k0 + k000 )× 1389 (k , k )(ζ 2 )472¯ 8 0 000 0 000 0 000 0 000 × (ζˇ )−1 9¯ 4 (k + k ) 8¯ 9 (k + k )κ¯ ¯ ¯¯¯ (k + k , k , k )κ¯ 3 489567

(E.29)

Static correlations can be calculated with a method described in [11]. A modified equilibrium probability density is used which describes a fluid which moves with a uniform velocity u. Expectations of the conserved quantity densities and their fluxes are formulated; derivatives of several order with respect to the the conjugate parameters are taken; it is shown that the resulting formulas for u = 0 are the static correlations. Many elements of the static correlation matrices are zero. As a first step of application of the described technique, one needs a means for finding these zero elements. The rule is this: Let n be the order of the static correlation. Let η be the number of (generally greek) indices whose values are not latin. Build n + η.Write e for the case that this number is even and o that it is odd. Write, for short, jn and sn for the two types of static correlations which appear in the formulas. Then: n+η jn sn

e o 6= 0 0 0 6= 0

(E.30)

6= 0 means ’generally not zero’. For the majority of coefficients in (E.29) all indices are latin, so η = 0. We find that the latin elements of s4 and j3 are zero; so, in (E.29), the terms 1 an 4 as well as the second parts of terms 7 and 9 are zero. The wave number integrals in terms 7 and 9 remain being a problem. We will assume here that they are somehow related to the localized value of the integrand (i. e. the integrand with all wave numbers equal to zero). In term 9, all indices of Φ are latin. With j5 and j3 equal to zero, Φ is zero too. In term 7, indices 4, 7, 8 are latin while 5, 6 are either both latin or both non-latin (1 oder 2), since otherwise the second factor of the integrand will be zero. In the first case, Φ = 0. In the second case, in (D.7), the second term is zero; in the first term, for j5 , η = 2 , and for the third term, for j3 , η = 2 also; then, from (E.30) both coefficients are zero. Therefore, in the wave number integrals 33

in terms 7 and 9, the localized value of the integrand is zero. This does not necessarily include that the integrals vanish. But, with this argument, we will neglect them. - The formula for ∆3 then simplifies: ∆3 (k, k0 , k00 , k000 ) = γ 04 (k)( 31 (j4 )4123 − (2π)3 δ(k − k0 )δ 14 δ 23 )

00 0 000 + 2(s3 )056 (k)((s3 )134 (k0 ) − (s3 )413 (k0 + k000 ))(ζ 2 )−1 5624 (k , k + k )

000 0 00 00 0 0 −1 + 2(s3 )047 (k)(s3 )856 (k0 + k00 )(ζ 2 )−1 4738 (k , k + k )(ζ 2 )5629 (k , k )κ19 (k ) 1 000 0 00 00 0 0 00 −1 + (s3 )047 (k)(s3 )589 (k0 + k00 )(ζ 2 )−1 4736 (k , k + k )(ζ 2 )8921 (k , k )κ56 (k + k ) 2 − 21 (s3 )047 (k)(s3 )568 (k0 + k000 )(ζ 2 )−1 (k00 , k0 + k000 )× 472¯ 9 × (ζˇ )−1 (k0 + k000 , k0 , k000 )κ¯99 (k0 + k000 ) 3 913568

0 000 −1 + 21 (s3 )047 (k)(s3 )¯5¯6¯7 (k0 + k000 )(ζ 2 )−1 (k00 , k0 + k000 )× 1389 (k , k )(ζ 2 )472¯ 8 0 000 0 000 0 000 0 000 × (ζˇ )−1 8¯ 9 (k + k )κ¯ 9¯ 4 (k + k ) ¯ ¯¯¯ (k + k , k , k )κ¯ 3 489567

(E.31)

The non-zero static correlation components are calculated with the technique outlined above. One obtains (all indices latin): (j4 )1234 = ξ −1 (δ 12 δ 34 + δ 13 δ 24 + δ 14 δ 23 ) (E.32) 1 (s3 )123 (k) = ik4 √ (δ 12 δ 34 + δ 13 δ 24 − λδ 14 δ 23 ) ρβ λ=

∂β p|ρ ∂β ε|ρ

(E.33) (E.34)

p is the static pressure and ε the energy density, both as functions of mass density ρ and inverse kinetic temperature β. Via λ, the formulas are dependent on the physical properties of the fluid. We will show that the second-order term (the last row of (E.1)) just cancels the first row of (E.31). The part of D3 which contains this term simplifies somewhat when (E.32) is introduced: 1 1 − ρ2 β 2 (2π)6

Z

dk0 dk00 dk000 u1 (k0 , t)u4 (k00 , t)u4 (k000 , t) δ(k − k0 − k00 − k000 )× × γ 01 (k)(ξ −1 − (2π)3 δ(k − k0 ))

(E.35)

The quantity K2 (E.7) has been processed by introducing the formulas for C3 (6.8), (6.9) and its derivatives (E.14) to (E.18). Then, G12 does not explicitly show up any more, and K2 now reads: (K2 )01 (k0 , k00 , t) = (γ 01 (k)(2π)3 δ(k00 )hh2 i − γ 04 (k) j412 )δ(t) − i k600 (C2 )56 (k00 , t)(C2 )41 (k0 , t)s045 (k)
+ (C2 )52 (k00 , t) −(C2 )43 (k0 , t)s045 (k)γ 31 (k0 ) + (C2 )41 (k0 , t)s067 (k)(ζ 2 )6745 ((k0 , k00 )

(E.36)

Remember that the component index 2 has the fixed value 1. The formula now consists of 3 terms. 0 00 The second term, after time integration, contains the factor k600 (ζ 2 )−1 4156 (k , k ) which is zero. The third term shows a C2 factor one index of which has the fixed value 1. When the analysis is restricted to shear modes, as is introduced in the previous section, this term does not contribute. Thus, in (E.36), the first term remains. - We can write h2 = h1 , where on the left the index variable 2 appears and on the right 1 1 its fixed value 1; we have hh1 i = ξ 2 . For the coefficient of the second partial term we find j412 = ξ − 2 δ 41 . The last row of (E.1) now reads: Z 1 2 1 ρ β dk0 dk00 u1 (k0 ,t)uq(k00 t) δ(k − k0 − k00 )γ 01 (k)(ξ −1 − (2π)3 δ(k00 ) ) (E.37) 2 (2π)3 The main difference to (E.35), apart from the order of integration, is the appearance of uq. The relation to the u’s is: 34

uq(k, t)

= =

Z

dxe−ik·x u4 (x, t)u4 (x, t) Z Z 1 0 dk00 δ(k − k0 − k00 )u4 (k0 , t)u4 (k00 , t) dk (2π)3

(E.38)

When this is introduced into (E.37) (and some wave numbers are renamed), it is seen that (E.37) actually cancels (E.35), as was stated in the beginning of the paragraph.

F

Final form of the 3rd order term

After the reduction to shear modes, for the inverse matrices (ζ 2 )−1 , (ζˇ3 )−1 we obtain: 0 00 (ζ 2 )−1 bcde (k , k ) =

0 00 000 (ζ˜3 )−1 bcdef g (k , k , k ) =

εbd (k0 )εce (k00 ) k 02 k 002 ν(k 02 + k 002 )

(F.1)

εbe (k0 )εcf (k00 )εdg (k000 ) k 02 k 002 k 0002 ν(−k 02 + k 002 + k 0002 )

(F.2)

When this is inserted into (E.31) and then in (E.1), after expansion some terms which stem from the second term of (5.2) vanish because of the continuity condition. Then, it becomes apparent that the last three terms of (E.31) cancel each other, while the two remaining terms can be somewhat simplified. The final form in wave number space is: ρ 1 Da (k, t) = ν (2π)6

Z

dk0 dk00 dk000 ub (k0 , t)uc (k00 , t)ud (k000 , t) δ(k − k0 − k00 − k000 )× × ∆abcd (k0 , k00 , k000 )

(F.3)



(F.4)

∆abcd (k0 , k00 , k000 ) = (s3 )ace (k0 + k00 + k000 )εef (k0 + k000 )×  × −(s3 )bdf (k0 )

1 + (k0 + k000 )2 (k 002 + (k0 + k000 )2 )

+(s3 )f bd (k0 + k000 )

k 0002 0 000 2 002 (k + k ) (k + (k0 + k000 )2 )(k 02 + k 0002 )

We extracted the thermodynamic factors from the s3 coefficients in order to collect them in front of (F.3). We now have: (s3 )bcd (k) (s3 )bcde

= ike (s3 )bcde = δ bc δ de + δ bd δ ce − λδ be δ cd

(F.5)

The quantites on the left are distinguished by the number of indices.

G

Formula for the velocity function

The formula for the velocity function q0 (the part without a trigonometric factor) of the part D1 of the third order dissipation term in 10.5 has been calculated with Mathematica, brought into TraditionalForm (partial derivatives in superscript representation), copied as a graphic and shown here as Fig. 11:

35

Fig. 11: Formula for q0

References [1] J. M. McDonough, Introductory lectures on turbulence. Dpt. of mechanical engineering an mathematics, Univ. of Kentucky [2] W. D. McComb, The physics of fluid turbulence, Oxford Science Pub., Oxford 1990 [3] K. Huang, Statistical Mechanics. Wiley, N. Y. 1963 [4] H. Grabert, Projection operator techniques in nonequilibrium statistical mechanics. Springer, Berlin, Heidelberg, New York (1982) [5] J. Schofield, R. Lim, I. Oppenheim, Physica A181, p. 89, 1992; [6] R. van Zon, J. Schofield, Phys. Rev. E65, 011106, 2001 [7] J. Piest, arXiv 0803.3972, 2008 [8] Zubarev, D.; Mozorov, V.; R¨opke, G.: Statistical mechanics of nonequilibrium processes. Vol. 1, Akademie Verlag Berlin 1996. [9] K. Kawasaki, Ann. Phys. 61, (1970), 1 [10] M. H. Ernst, E. H. Hauge, J. M. J. van Leeuwen, Phys. Rev. A 4, 5, (1971), 2055 [11] A. M¨ unster: Statistische Thermodynamik. Springer, Berlin 1956 [12] H. Schlichting, Boundary layer theory, 7

th

ed., McGraw Hill 1979

[13] E. Zauner: J. Fluid Mech. 154, 111-119, 1985 [14] W. Schneider: J. Fluid Mech. 154, 91-110, 1985 [15] W. Schneider, E. Zauner, H. B¨ohm: J. Fluids Eng. 109, 237-241, 1987 36

[16] E. N. da C. Andrade, L. C. Tsien Proc. Phys. Soc. London 49, 381-391, 1937 [17] H. B¨ohm, E. Zauner, W. Schneider: Numerical Investigation of viscous flow induced by an axisymmetric laminar jet issuing form a plane wall. Num. Meth. in laminar and turbulent flow. Swansea, 1985, 96-106E [18] A. J. Reynolds, Observations of a liquid-into-liquid jet, 1962 [19] L. M. Milne-Thomson: Theoretical hydrodynamics. Reprint, Dover Publications, N. Y. 1996 [20] J. Happel, H. Brenner: Low Reynolds number hydrodynamics. 2nd ed., Noordhoff, Leydon 1973 [21] B. Hamel, E. Richter, Luft- und K¨altetechnik 1979 Heft 1 [22] G. W. Rankin, K. Sridhar, M. Arulraja, K. R. Kumar, J. Fluid Mech. 133, 21t7-231, 1983

37