ELSEVIER
Dynamics of Atmospheres and Oceans 27 (1997) 233-241
Eddy viscosity and the statistical theory of turbulence J.R. Herring NCAR, Box 3000, Boulder, CO 80307, USA
Received 31 December 1995; revised 30 September 1996; accepted 14 October 1996
Abstract We examine, for the case of stationary turbulence, the correlation between fully resolved turbulence (in the sense of Fourier modes) and a system composed of the same equations on a reduced wavenumber span. This 'reduced system' is subjected to a wavenumber-dependent eddy viscosity contrived so that the full and reduced systems have the same energy spectrum, on the reduced span. The particular numerical problem studied is isotropic turbulence, and the model of turbulence is a Langevin representation of the test-field model (Kraichnan (1971, J. Fluid Mech., 47: 513-524). The correlation between the full and reduced systems is surprisingly small, a fact we attribute to be related to the unpredictability of the underlying Navier-Stokes equations. The (qualitative) form of the eddy viscosity introduced by such modeling is also discussed, and it is noted that at large Reynolds numbers, there exists an optimal cut-off wavenumber which permits an elimination of explicit empiricism from this form of large eddy simulation (LES). Our computations permit a spectral estimate of the magnitude of the backscatter term in LES, along the lines recently proposed by Schumann (1995, Proc. R. Soc. London, Ser. A, 451: 293-318). ©1997 Elsevier Science B.V. © 1997 Elsevier Science B.V.
1. Introduction Geophysical flows have orders of magnitude more degrees of freedom than can be accommodated by supercomputers. Thus, only limited statistical information can be obtained from their numerical simulations. One standard technique of obtaining such statistics is to average the Navier-Stokes equations over small, contiguous volumes, so that the number of equations pertaining to averaged degrees of freedom fits into the computer's memory. This is large eddy simulation (LES), by which is meant that only scales larger than the averaging volumes are retained in the computer. The effects of scales smaller than the averaging volumes on those retained are represented by an eddy viscosity. In addition to damping terms, small scales also perturb the larger ones in a 0377-0265/97/$17.00 © 1997 Elsevier Science B.V. All fights reserved. PII S 0 3 7 7 - 0 2 6 5 ( 9 7 ) 0 0 0 1 2 - 2
J.R. Herring/Dynamics of Atmospheres and Oceans 27 (1997) 233-241
234
way inaccessible from the perspective of the retained scales. Such effects, called 'backscatter', are usually treated as random perturbations of the large, retained scales. Another technique for obtaining statistical information is to form equations for ensemble averaged moments of Navier-Stokes. Moment closure is often called spectral closure, as the key ingredient is the two-point moment, R i j ( r , r ' , t ) - (ui(r,t)uj(r',t)), and Rgj is most economically characterized (in the homogeneous context) by its Fourier transform with respect to r - r'. As moments are average properties, they are smoother than the velocity field itself, and consequently, their simulation may be done with fewer numerical degrees of freedom. For example, an average of the velocity field ( u ( r , t ) ) and its moments Ro(r,r',t) has a smoother dependence on (r,r') and t than u(r,t) itself. Hence interpolation may be used to economize the representation of moments with impunity. Neither LES nor spectral closure is, at present, able to make a priori error estimates. An alternative approach to those described above is to replace the velocity field with simple stochastic equations, and constrain these to be consistent with certain known properties of moments, and other stochastically specified information. These constraints should include consistency with energy conservation, and the known scale-size distributions of the turbulence. Such stochastic models must also provide a representation of the cascade of energy from large to small scales, where molecular dissipation destroys the eddies. In such Langevin modeling, interactions among pairs of modes with similar scale randomly excite a given degree of freedom (or mode), whereas interactions of the mode in question with much smaller-scaled degrees of freedom tend to damp it. In the language of this model, the collision of a pair of modes is the essential origin of 'backscatter' mentioned above. (Pair interactions stem from the bilinear form of Navier-Stokes.) Langevin modeling differs from the traditional LES counterpart in that it contains the same number of degrees of freedom as Navier-Stokes, but the terms representing nonlinearities make its equations of evolution much easier to solve than the latter. In addition, Langevin models may be formulated so that their predictions for moments are entirely consistent with a given moment closure. Hence, they may be viewed as equivalent formulations of moment closures, which may offer computational advantages or insights into the underlying dynamics of moment closures. Recently, Schumann (1995) has proposed a method of computing the 'backscatter' term that enters LES. We shall discuss here an extension of these ideas to the more detailed spectral domain. To capture the gist of the problem, we consider first the simple Langevin model of turbulence, whose velocity field we denote as u(t). We model u by a Langevin equation,
d u / d t = - I~u + F( t) (1) Here, /z represents eddy damping, and F(t) is a random force, which represents both the forcing necessary to maintain a steady state, and the 'turbulence force' that arises from the nonlinear terms in Navier-Stokes. We suppose that u contains all existing scales of turbulence, resolved as well as unresolved. Now we let v denote a reduced system, designed to reproduce as much of the statistics of u as is practical. We shall here be content with variance statistics, and shall require that
(u2)(t)
= (v2)(t)
(2)
J.R. Herring/Dynamics of Atmospheresand Oceans27 (1997) 233-241 We know from predictability studies that R = ( u u ) / / ~ l ( ( u 2 ) ( u 2 ) ) an equation of the form
d v / d t = -I.dV + F ( t ) + F'( t)
235
<( 1. For v we write
(3)
where the augmented value of the damping and the existence of the new random force, F', represent the damping and backscatter attributable to the omission of the unresolved scales in Eq. (3). For the steady-state form of Eq. (1) and Eq. (3) we deduce
L~ds(F'(t)F'(s))
= 2(1/R-
1)f_~ds(F(t)F(s))
(4)
where we have used Eq. (2). Hence the decorrelation associated with the unpredictability of turbulence is, in the sense noted here, equivalent to an extra stirring force in a rough, Langevin model of turbulence. Moreover, the magnitude of the force F' is derivable from a knowledge of the decorrelation R. It is frequently remarked that R << 1 does not necessarily vitiate LES, as much of its error is associated with the divergence of eddies with increasing time. The structures and their intensity computed by LES may then be trusted, so that its statistics could be accurate. We shall not comment further on this point except to note that both phase errors and structural errors contribute to the reduction of R in Eq. (4). We shall return to our discussion of the backscatter term after presenting a spectral estimation of the spectral version of Eq. (4) via a statistical theory, the test-field model (TFM) (Kraichnan, 1971; Herring and Kraichnan, 1972). LES may be formulated in a variety of ways, but one simple and precise manner is in terms of a spectral closure, as was done by Kraichnan (1976) and Chollet and Lesieur (1981). In this case, one asks for an eddy viscosity that, if applied to Fourier modes u(k,t) on a reduced domain 0 < k < k c (called hereafter ~ ) , yields velocity amplitudes whose spectra on _~ are identical to a spectrum that would be found for an unreduced or full domain (0 < k < oo). (From here on, we call the reduced domain velocity amplitude v(k,t), and the full-domain velocity field u(k,t).) In such a formulation, all that is claimed about v(k,t) is that its covariance is related to that of the untruncated u(k,t) in the manner just outlined: ( u ( k , t ) u ( - k , t ) ) = ( v ( k , t ) v ( - k , t ) ) , Ikl <- kc, v(k,t) = O, Ikl > kc, ~ = 0 < Ikl < k~. Thus, no claim is made that the amplitude, v(k,t), contains details of the actual flow. In particular, the reduced system's skewness (which measures the net energy transfer out of the energy-containing range into the dissipating range) will be much smaller for the reduced system. Hence, this method avoids the predictability dilemma inherent in the filtered velocity amplitude approach. In this representation of LES, interactions among the reduced system's Fourier modes occur both in the eddy-viscosity and in the advective-pressure term. The latter is just u • Vu - Vp projected (in the sense of Galerkin) onto ~ . This is an identical form for the acceleration to that for the unreduced equations, provided, of course, that the modes exterior to .~ are excluded. It is not clear that this familiar form for the advective non-linearity is that which would emerge from a normal form analysis of the equations of motion. Nevertheless, it seems plausible that the dynamics of v resemble those of u; hence, the folklore that v is to be regarded as a 'typical' (low-pass-band-filtered) realization of u.
236
J.R. Herring~Dynamics of Atmospheres and Oceans 27 (1997) 233-241
The LES program just outlined assumes the solution to Navier-Stokes on the full domain, but such knowledge is as yet far beyond analytic or computationaI capacity. However, we can apply the program cleanly to a Langevin model, which is designed to capture certain variance-level information of Navier-Stokes. It is the application of the program to this surrogate system that constitutes the subject of the present paper. We note, however that similar studies of the Navier-Stokes dynamics (at small Reynolds numbers) have been made by Merilees (1975) and others. Here we use spectral closure (TFM) (Kraichnan, 1971; Herring and Kraichnan, 1972) to estimate ~ ( k ) = (u(k,t). v(-k,t))/((]u(k,tl2))l/2(Iv(k,t)[z)J/2. Results for typical large-Reynolds-number flows indicate a surprisingly small value, approximately 0.3 for ~.~(k) across a broad range of k. Section 2 develops equations for the normalized correlation coefficient ~ ( k ) , and the eddy viscosity for stationary turbulence. Such equations are exact for the Langevin model, and depend on the fidelity of the TFM variance to the actual dynamics for their verisimilitude.
2. Langevin equations To implement this program for flows with finite (but large) microscale Reynolds numbers, Ra, we employ the Langevin model realization of the TFM as discussed by Kraichnan (1971), Herring and Kraichnan (1972) and Leith and Kraichnan (1972). We use a compact notation in which both the vector and wavenumber degrees of freedom (k) are tagged by i. Variance spectra are denoted by U, = (u~). Langevin-model equations of motion for u i are
(a, +
=
(t)Ec,.,.,
0i.
(5)
and in the reduced domain 2 , (V, = (L,~)) the equivalent equation is +
+
= o(OEc,.,.,
" v'Oii, ~i~, + fi( t)
(6)
In Eq. (5) and Eq. (6), ~r is specified by (~(t)o'(t'))= 6 0 - t ' ) , and ~(i) is an independent white-noise forcing, which allows U~ and ~ to become steady in time. Eddy dissipation, v~ddy, in Eq. (6) is fixed by the condition that V~= Ui, i ¢ 2 . /z, in Eq. (5) and /z~ in Eq. (6) are prescribed so that ensemble mean energies ½(u~) and 7(vit2) are conserved, if (bare) viscosity and forcing are turned off (this requires that /zi = 4Y'vtCijicj, utoijl). The equation of motion for U~ follows from Eq. (5) as (8 t + 2 p,i) U, =
2EC~j,Oij,UiUI + 2F,
where /, - f~dt'(fi(t)fi(t')). domain ~ for V/ is (a t q- 2/.Zff + 2 peddy)~ =
(7)
An equivalent equation corresponding to Eq. (6) on the
2ECii, Oij,~Vi +
2/,
(8)
J.R. Herring~Dynamics of Atmospheres and Oceans 27 (1997) 233-241
237
Defining the modal energy transfer function ~y-v.v as
g-i v = 2EC, j, gUtO~jt- 2 ~ U /
(9)
and
Y i v = 2~Ci~,VjV I Oij~ - 2 p~/~V~
(10)
we note that the equations of motion for U~ and V~ may be more simply written as ot~ = y u + 2Fi
(11)
oyi = J i v + 2F,
(12)
Requiring V/= U/, i ~.@ then gives for v eaay
(13)
vyY = (g,~ - ~)/(2u,)
We may also use Eq. (5) and Eq. (6) to evaluate the covariance between u i and v i, which we call Wi = (uivi): _ 2 (O t -{" ~.Li + $.L/~ + /,,'TddY)w/ - 2 E f i j l ~
..~
W j Wl + 2 F i
(14)
Putting Ai-~- Ui - Wi,
~Li "~ [£i"~ + 12i : 2['Li~ -~- ~i
= Oij~(1 - Pij,)
(15a) (15b)
we may write the steady-state equation for Ai as
(21~+Ei)Ai=e,+ECi2t{Oij~tpij, U j U , + ~ ( 2 U j A t - A j A , ) }
(16)
where we have used the symmetry CO~ = Cil j. It shoule be noted that Eq. (15a) defines Ei •
3. Results Fig. 1 shows the energy spectrum for a stationarily forced system, in which the forcing is confined to a small range of wavenumbers around k = 0.5 (for details, see the figure caption). Inferior to the forcing scales, we note a thermal equilibrium k 2 range, and for k > 0.5, the spectrum develops an extensive inertial range, with a spectral 'bump' in the decade of scales inferior to the dissipation range. The 'bump' region is not inconsistent with E ( k ) = k -3/2. Here, the Lagrangian relaxation rates comprising O(k,p,q) are computed as is (/xi(k)) (see above). Several values of the eddy viscosity, Veddy, according to Eq. (13) are given in Fig. 2 for various k c. Fig. 3 compares T(k) for the truncated system with that for the full system. The value for Veddy(k) is small in the energy-containing range, and may become slightly negative before the development of the cusp near k c, if k c is near the 'bump' region. We note a close fidelity of T~(k) to T(k) over the energy-containing scales. An
238
J.R. Herring~Dynamics of Atrnospheres and Oceans 27 (1997) 233-241 10~
' ' '"'"1
' I I'IIIll
I
i i i IIrn I !
'''""1
lO-Z ! 10 ~ E(k)
10 ~ r 10 .5 10 -6 10 -7 ]
10 ~
!
10 .9
104o
~ ,~,,,,,I
10 1
, i liiliil
10 0
l ,,,,,,,I
101 k
i i i iiHiI
10 2
10 3
Fig. 1. E(k ) = 2 zr k 2U(k) for stationarily forced isotropic turbulence, according to TFM, at Ra = 1700. Here, the modal forcing, F ( k ) = k 4 e x p [ - ( k - k f ) / w ] 2, k r = 0.5, w = 0.4 is applied to U(k) as in Eq. (7), and u = 6 . 3 2 × 10 5.
i n t e r e s t i n g p o i n t is that i f w e c h o o s e k c j u s t a b o v e t h e w a v e n u m b e r
k 0 at w h i c h T(k)
c r o s s e s f r o m n e g a t i v e to p o s i t i v e v a l u e s , w e m a y e s t i m a t e 2Ueady(k)=- T c ( k ) / ( 2 k 2) o n t h e r a n g e k 0 < k _< kc, a n d z e r o o t h e r w i s e , as T(k) in this r e g i o n is n e g l i g i b l e at l a r g e R~. I f this r e m a r k h o l d s f o r D N S , w e m a y d i s p e n s e w i t h t h e t h e o r y , a n d rely e n t i r e l y o n D N S f o r t h e c o m p u t a t i o n o f Ueddr. I f h o w e v e r , w e m a k e k s >> k o, T~(k) is c o m p a r a b l e
30
' ' '"'"1
' ''"'"1
' ' '"'"1
' ' '"'"l
I i i
20
i
i
i i 'Ve~ly
10
;
i
:
I
: :
i , i
":
#~
2 !"
0
_~:z]
t
10 -1
,,
i,,,,ll
10 0
i,
,,,H,I
.-'
, , ...... I
101
10 2
. . . . kc=lO0 t . . . .
__ k c = 1 0 1 I l ililill
"1
10 3
k Fig. 2. Eddy viscosity Peddy(k), according to Eq. (13), for k c = 10, 50, 100.
J.R. Herring/Dynamics of Atmospheres and Oceans 27 (1997) 233-241 i
i iii1,1]
,
i
....
I i i,li~
l
, ,i;,,,[
,
239
, ,ill,,[
.-°
i T(k)
-2
......
-4
-
i
i
Reduced Full
-
i
lIHHi
10-1
101 k
100
t i,illl
102
10 3
Fig. 3. Comparison of energy transfer function T(k)=- 2 ~ ' k 2 ~ - ( k ) for the full system with that for the reduced system, Tc(k), for k c = 10.
with T(k) near k c, and the theory must be invoked to compute Veddy(k). S i m i l a r comparisons at low R~ cannot be used in this way, as the positive value of T(k) is comparable with the negative, energy-containing region, and this height is not known from the reduced system. Fig. 4 shows 1 -~9~(k) = 1 - A(k)/U(k) for several k c. We note larger ~ ( k ) near the energy peak, where the forcing f~ is identical for both systems. Otherwise, the
1
' ' '"'"1
' ''"'"1
' ''"'"1
' ''"'"1
0.8
0.6 1-R 0.4
. . . . kc= 10I 0.2
.... kc=50-" ~kc=lO
0
i
10 "1
i i,ii,,[
i
100
, i1,,1,[
,
10 ~ k
i iiiii1[
i
102
i
ll[lld
103
Fig. 4. Error function 1 - .9~'(k) according to Eq. (16) for k c = 10, 50, 100.
240
J.R. Herring/Dynamics of Atmospheres and Oceans 27 (1997) 233-241
normalized error is substantial everywhere. Moreover, other calculations at a variety of larger k c (not shown) indicate that a ( k ) is not significantly reduced, until k~ ~ k d -= ( ~ / v 3 ) j/~, ~ - f~p2dpE(p)"
4. Discussion and concluding comments The magnitude of the reduced system's error, and the fact that it is not significantly decreased until k, ~ k d, is somewhat surprising. However, our results here are similar to those of Merilees (1975), who found that increasing the resolution of an eddy viscosity formulation of a simple baroclinic flow resulted in little error reduction. The discrepancy noted here has nothing to do with the presence of structures in the unreduced system and their lack in the reduced system, as the theory used is essentially structureless. We may relate the fact that A(k) saturates at k > kf with the nature of its transfer function as given in Eq. (16). It is known, by asymptotic analysis of the linear input term, that approximately 2 ~ A / leads to a production of A ( k ) = k r (Herring, 1984; M6tais and Lesieur, 1986). As the energy spectrum is steeper than this, saturation must occur at large k, in which the nonlinear term Aj A/) prevents A(k) from exceeding U(k). We have noted that if the cut-off wavenumber k~ is just larger than k o, (k o ~ T(k o) = 0), the evaluation of veddY(k) is approximately T~(k)/(2k2), k 0 _< k _< k~,, and zero elsewhere. Hence, the computation of v edoy is totally contained in the reduced system, and if this qualitative behavior holds for DNS as well, we have a formula for eddy dissipation without empiricism. We may call such a procedure a 'boot-strap' eddyviscosity prescription, as no theoretical closure assumptions need enter the evaluation of v~ddY(k). (Such a boot-strap method is similar to the 'constrained Euler system' proposed by Jim6nez (1993) and She and Jackson (1993).) However, we should note that for decaying turbulence, k 0 is a fraction (about one-fifth) of the dissipation wavenumber, k~ =- ( ~ / v 2 ) 1/4, whereas for the forced problem, k o / k d is much smaller. Hence, the computational load of implementing the boot-strap eddy viscosity program for decaying turbulence may be higher than for the forced system described here. However, such comments apply for the closure-based eddy viscosity formula, not necessarily for LES in general. Finally, we note that what we wish to examine is ,~(k) for Navier-Stokes, yet our calculations here are posed in terms of a Langevin model with enhanced dissipation and random forcing. We have thus posed the question of how to arrive at a reduced system and evaluate its parameters in terms of a surrogate model, and then expect to apply the results to Navier-Stokes (both full, and reduced) expecting edification. Why is this a reasonable expectation? In this connection, we note that second-order statistics are the same for Langevin models as for results of renormalization perturbation theory (such as embedded in TFM) despite the fact that the former has strictly Gaussian statistics, and the latter non-Gaussian statistics. It is known that certain higher-order statistics of the latter resemble those of Navier-Stokes (Chen et al., 1989). This statement applies to the velocity derivative skewness, and to those fourth-order moments needed for the computation of the velocity acceleration.
J.R. Herring/Dynamics of Atmospheres and Oceans 27 (1997) 233-241
241
Acknowledgements Part o f this w o r k w a s c o m p l e t e d w h i l e the a u t h o r w a s a v i s i t i n g P r o f e s s o r at the L a b o r a t o i r e de M o d e l i s a t i o n e n M r c a n i q u e h l ' U n i v e r s i t 6 P i e r r e et M a d e Curie. I a m g r a t e f u l f o r s e v e r a l b e n e f i c i a l d i s c u s s i o n s w i t h P r o f e s s o r M. L a r c h e v ~ q u e o f t h a t University.
References Chen, H.D., Herring, J.R., Kerr, R.M., Kraichnan, R.H., 1989. NonGaussian statistics in isotropic turbulence. Phys. Fluids A 1 (11), 1844-1854. Chollet, J.P., Lesieur, M., 1981. Parameterization of small scales of three-dimensional isotropic turbulence utilizing spectral closure. J. Atmos. Sci. 38, 2747-2757. Herring, J.R., 1984. The predictability of quasi-geostrophic flows. AIP Conf. Proc. 106, 321-332. Herring, J.R. and Kraichnan, R.H., 1972~ Comparison of some approximations for isotropic turbulence. In: M. Rosenblatt and C. van Atta (Editors), Lecture Notes in Physics, 12. Springer, Berlin. Jimrnez, J., 1993. Energy transfer and constrained simulations in isotropic turbulence. Annual Research Briefs --1993, Center for Turbulence Research, Stanford, CA, pp. 171-186. Kraichnan, R.H., 1971. An almost markovian Galilean invariant turbulence model. J. Fluid Mech. 47, 513-524. Kraichnan, R.H., 1976. Eddy viscosity in two- and three-dimensional turbulence. J. Atmos. Sci. 33, 1521-1536. Leith, C.E., Kraichnan, R.H., 1972. Predictability of turbulent flows. J. Atmos. Sci. 29, 1041-1058. Merilees, P.E., 1975. The effect of grid resolution on the instability of a simple baroclinic motion. Mon. Weather Rev. 103, 101-104. Mrtais, O., Lesieur, M., 1986. Statistical predictability of decaying turbulence. J. Atmos. Sci. 43, 857-870. Schumann, U., 1995. Stochastic backscatter of turbulence energy and scalar variance by random subgrid-scale fluxes. Proc. R. Soc. London, Ser. A 451,293-318. She, Z.-S., Jackson, E., 1993. Constrained Euler system for Navier-Stokes turbulence. Phys. Rev. Lett. 70, 1255.