Universal behavior of crossover scaling functions for continuous phase transitions

a r X i v :0708.4157v 1 [m a t h .C V ] 30 A u g 2007ON THE SCALING LIMIT OF A SINGULARINTEGRAL OPERATOR 1I.M.Krichever and D.H.PhongDepartment of MathematicsColumbia University,New York,NY 100271IntroductionRecently an alternative approach to Birkhoff’s theory of difference equations [1]has been proposed in [4].This approach leads naturally to local monodromies of difference equa-tions,which should converge in principle to monodromy matrices of differential equations,thus providing a missing link in the theory of isomonodromic transformations of systems of linear difference equations (see e.g.[2,3,7]and references therein).The key to the convergence process in [4]is the scaling limit of a certain singular integral operator I ,arising from a Riemann-Hilbert problem.The operator I acts on functions φdefined on the vertical line with fixed abscisse at a ,and its kernel k (z,ξ)is given explicitly byk (z,ξ)=e πi (z −a )+e −πi (z −a )e λ(y −η)−e −λ(y −η)e −λy +e λy1Research supported in part by National Science Foundation grants DMS-02-45371and DMS-04-055192Schauder estimates with exponential growthWe introduce the following norms of Schauder type for functions on R.Fixκ∈R,m∈Z, 0<α<1,and letΛα(m,κ)be the space of functionsφon R satisfying the conditions|φ(x)|≤C(1+|x|)m eκ|x|,|φ(x)−φ(y)|≤C|x−y|α{(1+|x|)m eκ|x|+(1+|y|)m eκ|y|},(2.1)for all x,y∈R.We define||φ||Λα(m,κ)to be the infimum of the constants C for which theseinequalities hold.We also require the spaceΛα(log,κ)and the corresponding norm||φ||Λα(log,κ)defined by the conditions|φ(x)|≤C log(1+|x|)eκ|x|,|φ(x)−φ(y)|≤C|x−y|α{log(1+|x|)eκ|x|+log(1+|y|)eκ|y|},(2.2) The singular integral operator I can be expressed asI(φ)(y)=(e−y+e y)H(1e y−η−e−(y−η)ψ(η)dη≡limǫ→0 |y−η|>ǫ1e z−e−z.(2.5) Then the kernel K(z)is C∞(R\0),odd,and satisfies|K(z)|≤C |z|−1,if|z|≤1;e−|z|,if|z|>1.|∂z K(z)|≤C |z|−2,if|z|≤1;e−|z|,if|z|>1.(2.6)In particular,these are better estimates than for the standard Hilbert transform kernel K0(z)=z−1,and it follows at once that the operator H is bounded on the standard Schauder spaces(see e.g.[6,5]).To obtain estimates for the operator I,we need the boundedness of H on the above spacesΛα(m,κ),and this is provided by the following theorem: Theorem1Fix0<α<1,m∈Z.The operator H is bounded on the following Schauder spaces,||Hψ||Λα(m,κ)≤C m,α,k||ψ||Λα(m,κ),−1<κ<1.(2.7)Forκ=−1,we have the following bounds,for m∈Z,m≥−1,||Hψ||Λα(m+1,−1)≤C m,α||ψ||Λα(m,−1),m≥0,||Hψ||Λα(log,−1)≤Cα||ψ||Λα(m,−1),m=−1.(2.8) 2Proof.The method of proof is the standard method for Schauder estimates for singular integral operators.The only new feature here is the control of Hψ(x )for x large.In view of the fact that K (z )is odd and exponentially decreasing,we can writeHψ(x )=∞−∞K (x −y )(ψ(y )−ψ(x ))dy,(2.9)where the integrals on the right hand side are now convergent for ψ∈Λα(m,κ)with 0<α<1,κ<1.In particular,|Hψ(x )|≤|x −y |<1|x −y |−1+α((1+|x |)m e κ|x |+(1+|y |)m e κ|y |)dy +|x −y |≥1e −|x −y |((1+|x |)m e κ|x |+(1+|y |)m e κ|y |)dy.(2.10)These are clearly bounded for |x |bounded,so we may assume that |x |≥3.In this case,12|x |e−(1−κ)|z |(1+|x −z |)m dz +eκ|x ||x −z |<12|x |to write (1+|x −z |)m ≤2−m (1+|x |)m ,and the inequality still holds.Sinceκ<1,the desired bound follows for the first integral.Next,in the second integral,wehave 12|x |,and we can write|x −z |<12|x ||z |≤3for arbitrary N.This proves the lemma when0≤κ<1.When−1<κ<0,we write insteade−|z|eκ|x−z|=e−(1+κ)|z|eκ(|z|+|x−z|)≤e−(1+κ)|z|eκ|x|(2.16) and bound the integral on the left hand side of the lemma byeκ|x| |x−z|≥12||x|e−(1+κ)|z|(1+|x−z|)m dz.(2.17) The bounds for these integrals are now the same as in the previous case.This establishes the estimate(2.11).Finally,consider the caseκ=−1.In the region of integration |x−z|>4|x|,we have the integrand can be crudely bounded by e−2|x|(1+|x−z|)m e−1Since we have|K(x−y)−K(x′−y)|≤|x−y|·|∂x K(z)|(2.23) for some z in the segment between x−y and x′−y,and hence|z|∼|x−y|when |y−x|>3|x−x′|,we can write,in view of the bounds for the|∂z K(z)|,| |y−x|>3δ(K(x−y)−K(x′−y))(ψ(y)−ψ(x))dy|≤δ||ψ||Λαm,κ 3δ<|x−y|<1|x−y|−2+α{(1+|x|)m eκ|x|+(1+|y|)m eκ|y|}dy+δ||ψ||Λαm,κ |x−y|>1e−|x−y|{(1+|x|)m eκ|x|+(1+|y|)m eκ|y|}dy.(2.24) Thefirst integral on the right hand side is bounded byδ||ψ||Λαm,κ 3δ<|x−y|<1|x−y|−2+α{(1+|x|)m eκ|x|+(1+|y|)m eκ|y|}dy≤δ||ψ||Λαm,κ(1+|x|)m eκ|x| 3δ<|x−y|<1|x−y|−2+α≤Cδα(1+|x|)m eκ|x|.(2.25) Applying Lemma(1),we obtain similar bounds for the second integral.Altogether,we have shown that|Hψ(x)−Hψ(x′)|≤C||ψ||Λα(m,κ)|x−x′|α(1+|x|)m eκ|x|(2.26) for|x−x′|small,and the theorem is proved when−1<κ<1.The caseκ=−1is established exactly in the same way,using the corresponding estimates in Lemma1for κ=−1and m≥−1.The proof of Theorem1is complete.Q.E.D.Theorem2For0<κ<2,m∈Z,and0<α<1,the operator I is bounded on the following Schauder spaces,||I(φ)||Λα(m,κ)≤C m,κ,α||φ||Λα(m,κ),0<κ<2.(2.27)Forκ=0,m∈Z,m≥−1,the operator I satisfies the following bounds,||I(φ)||Λα(m+1,κ)≤C m,α||φ||Λα(m,κ),m≥0||I(φ)||Λα(log,κ)≤Cα||φ||Λα(m,κ),m=−1.(2.28)Proof.This is an easy consequence of Theorem1,the fact that the mapφ→ψ(y)= 13The scaling limit of IλWe come now to the operators Iλ.The estimates for I in the previous section show that Iλcannot be treated by simple scaling arguments from I.Instead,we shall study the bounds and limits for Iλasλ→+∞directly.It is simplest to carry this out for functions φsatisfying conditions of the form,|∂lφ(x)|≤C k(1+|x|)m,0≤l≤k,(3.1)forfixed m∈N,k∈N,and norms||φ||Λk(m)defined to be the best constant C k for whichthe above condition holds.The following theorem describes the limit of Iλin these spaces, although it should be clear from the proof and from the previous section that other more precise versions can be formulated as well:Theorem3Fix m∈N.Then we have the following bounds,uniform inλand inφ∈Λ1(m),||Iλ(φ)(y)− y0φ(ξ)dξ||Λ0(m+1)≤C mλ−1e−λη+eλη1eλt−e−λt(ψ(y−t)−ψ(y+t))dt+ |t|>y eλy+e−λyeλη+e−ληφ(η).(3.7)6To estimate(A),we apply Taylor’s formulaψ(y−t)−ψ(y+t)=t 1−1ψ′(y−ρt)dρ(3.8) which gives in this particular case,ψ(y−t)−ψ(y+t)=t 1−1dρ(1(eλ(y−ρt)+e−λ(y−ρt))2φ(y−ρt))(3.9) Thus(A)can be rewritten as(A)=1eλ(y−ρt)+e−λ(y−ρt)φ(y−ρt)≡A1+A0(3.10) where the functionχλ(ρ,t)is defined byχλ(ρ,t)=λteλ(y−ρt)+e−λ(y−ρt).(3.11)The following sharp estimates forχλ(ρ,t)play an essential role in the sequel: Lemma2For all0<t<y,the functionsχλ(ρ,t)satisfy the following properties(a)11−e−2λte−λt(1−ρ)<χλ(ρ,t)≤2λt2≤ 1−1χλ(ρ,t)dρ≤2.(3.12)Proof.In the region0<t<y,we have y−ρt>0for all|ρ|<1,and thus1eλ(y−ρt)+e−λ(y−ρt)≤2eλρt.(3.13)The upper bound implies(a),while the lower bound implies(b),when combined with thefollowing explicit formula1−1eλρt dρ=1Lemma3The term involvingφ′above tends to0at the following rate,|A1|≤C m1λ||φ||Λ1(m) y0dt(1+|y|m+|t|m)1−1dρχλ(ρ,t)(3.16)and the desired estimate follows from the statement(b)of Lemma2.Q.E.D.The estimates in Lemma2show thatχλ(ρ,t)provide an approximation of the Dirac measure concentrated atρ=1,1e−λ(y−ρt)+eλ(y−ρt)φ(y−ρt)−φ(y−t)(3.18) = 1−1dρχλ(ρ,t)(e−λ(y−ρt)−eλ(y−ρt)eλ(y−ρt)+e−λ(y−ρt)−1 |φ(y−ρt)|≤C e−2λ(y−t)||φ||Λ0(m)(1+y)m.(3.21)Proof.To prove(a),we writee−λy+eλy1+e−2λ(y−ρt) (3.22)≤eλρt(e−2λy+e−2λ(y−ρt))≤eλρt(e−2λy+e−2λ(y−t)).In particular,1−1dρχλ(ρ,t)−λt eλt−e−λt 1−1dρeλρt(e−2λy+e−2λ(y−t)).(3.23)8Since|ρ|<1dρe λρt =(λt )−1(e λt −e −λt ),the statement (a)follows.To establish the statement (c),we begin by noting thate−2λ(y −ρt )≤1−e λ(y −ρt )−e −λ(y −ρt )1+e −2λ(y −ρt )≤2e −2λ(y −ρt ).(3.24)Using the estimate for χλin Lemma 2and carrying out explicitly the integral in ρgives1−1dρχλ(ρ,t )eλ(y −ρt )−e −λ(y −ρt )1−e −2λte−λ(t +2y )1−1dρe 3λρt=11−e −2λt≤C e −2λ(y −t ),which implies immediately (c).To establish (b),let δ>0be any number sufficiently small and to be chosen suitably later.Write1−1dρχλ(ρ,t )(φ(y −ρt )−φ(y −t ))=1−δ−1+11−δ≡I δ+II δ.(3.25)The second term on the right hand side can be estimated by,|II δ|≤δt ||φ||Λ1(m )(1+y )m11−δdρχλ(ρ,t )≤δt ||φ||Λ1(m )(1+y )m,(3.26)while the first term can be estimated using Lemma 2,|I δ|≤2||φ||Λ1(m )(1+y )mλt1−e −2λt≤C sup [0,2y ]|φ|e −λδt .(3.27)The proof of Lemma 4is complete.Q.E.D.We can now carry out the integral in t .The precise estimates are given in the nextlemma:Lemma 5For any 0<y ,we have the following estimates,ydt1−1dρλt(e λ(y −ρt )+e −λ(y −ρt ))2φ(y −ρt )−ydt φ(y −t )≤C ||φ||Λ1(m )(1+y )m(yλ1Proof.In view of the defining formula(3.11)for the functionχλ(ρ,t)and the break up (3.18),the left hand side of the desired inequality is bounded by the integral in t of the three inequalities in Lemma4.This gives the following upper bound,||φ||Λ1(m)(1+y)m y0dt(e−λy+2e−2λ(y−t)+δt+e−δλt).(3.29) The integral can be evaluated explicitly,and wefindye−λy+12δy2+1λ(1−e−2λy)∼yδλ(3.32)which attains its lowest size yλ−12.This gives the estimate stated in the lemma.Q.E.D.We return now to the estimate of the contribution to Iλ(φ)(y)of the integral in t from the region|t|>y.Lemma6For any0<y,we have the following estimate|t|>y1eλ(y−t)+e−λ(y−t)φ(y−t)dt ≤C m1λy) Proof.Considerfirst the contribution from the region t>y.In this region,we have1eλ(y−t)+e−λ(y−t)≤2eλy1−e−2λt|φ(y−t)|dt= ∞0e−2λs1−e−2λ(s+y)(1+|s|m)ds.(3.34)10We claim that for all m∈N,we have∞0e−2λsλ(1+log(1+12λ 10du2λlog12and12 0du2λlog12(1+1u)m≤2u)m du+2m+1 1 2du2λlog12duλy).This is the desired estimate.Next,consider the contribution of the region t<−y.In this region,we have instead1eλ(y−t)+e−λ(y−t)≤2eλy1−e2λt (1+|t|)m dt≤C m1λy),(3.40)as was to be shown.Q.E.D.The bound provided by Lemma6involves a log(λy)−1term,and is not adequate for y close to0.This is because the integral is only a principal value integral when|t−y|is small,and the estimates we have just derived for the contribution of the region t>y do not take into account the cancellations inherent to principal value integrals.This issue is addressed in the next lemma:Lemma 7Assume that 0<λy <1.Theny<|t |<11e λ(y −t )+e−λ(y −t )φ(y −t )dt≤C ||φ||C 0[0,2](y +1e λt −e −λte λy +e −λyλ.(3.42)Proof.Since we can assume that λis large,the condition that λy <1implies that y <1,say.We can exploit the cancellation by writing the integral over the region y <|t |<1in the form,y<|t |<1=1y1e −λ(y −t )+e λ(y −t )φ(y −t )−e −λy +e λye −λ(y −t )+e λ(y −t )φ(y −t )−e −λy +e λye −λ(y −t )+e λ(y −t )(φ(y −t )−φ(y +t ))+φ(y +t ){e −λy +e λye −λ(y +t )+e λ(y +t )}(3.44)The contribution of the first term on the right hand side can be estimated as follows,1y1e −λ(y −t )+e λ(y −t )(φ(y −t )−φ(y +t ))≤||φ||C 1[−2,2]1yte λt −e−λt e λ(2y −t )dt ≤e 2e 2λt −1dt ≤Ce −λ(y −t )+e λ(y −t )−11+e−2λ(t +y )(1e λ(t +y ))+11+e −2λ(t −y )−11+e −2λ(t +y )(1eλ(t +y ))≤e −λt (e λy −e −λy )≤C λye −λt ,(3.46)and the corresponding integral in turn by,1ydt11+e −2λ(t +y )(1eλ(t +y ))≤||φ||C 0[0,2]λy1ydte 2λt −1=1e 2λt −1+11e 2λt −1≤C (1λt+11e 2λt)≤C1λy+1),(3.48)and hence,since λy <1,1ydt11+e −2λ(t +y )(1eλ(t +y ))≤||φ||C 0[0,2](1eλ(t +y )11+e −2λ(y +t )≤(e λy −e −λy )e −λ(3t +y )≤C λye −3λt ,(3.50)and hence1ydt1eλ(t +y )11+e−2λ(y +t )·|φ(y+t )|≤||φ||C 0[0,2]λy1ye −3λte λt −e −λte λy +e −λye λt1References[1]Birkhoff,G.D.,“The generalized Riemann problem for linear differential equationsand allied problems for linear difference and q-difference equations”,Proc.Amer.Acad.of Arts and Sciences49(1913)521-568.[2]Borodin,A.,“Isomonodromy transformations of linear systems of difference equa-tions”,Ann.of Math.160(2004)1141-1182.[3]Harnad,J.and A.R.Its,“Integrable Fredholm operators and dual isomonodromicdeformations”,Comm.Math.Phys.226(2002)497-530.[4]Krichever,I.M.,“Analytic theory of difference equations with rational and ellipticcoefficients and the Riemann-Hilbert problem”,Russian Math.Surveys59(2004) 1117-1154.[5]Ladyzhenskaya,O.and N.Uraltseva,“Linear and quasilinear elliptic equations”,Aca-demic Press,1968.[6]Stein,E.M.,“Harmonic Analysis:real-variable methods,orthogonality,and oscilla-tory integrals”,Princeton University Press,1993.[7]Tracy,C.A.and H.Widom,“Fredholm determinants,differential equations,and ma-trix models”,Comm.Math.Phys.163(1994)33-72.。

Constitutive model for the triaxial behaviour of concreteAuthor(en):William, K.J. / Warnke, E.P.Objekttyp:ArticleZeitschrift:IABSE reports of the working commissions = Rapports des commissions de travail AIPC = IVBH Berichte der ArbeitskommissionenBand(Jahr):19(1974)Persistenter Link:/10.5169/seals-17526Erstellt am:22.08.2011NutzungsbedingungenMit dem Zugriff auf den vorliegenden Inhalt gelten die Nutzungsbedingungen als akzeptiert. Dieangebotenen Dokumente stehen für nicht-kommerzielle Zwecke in Lehre, Forschung und für dieprivate Nutzung frei zur Verfügung. Einzelne Dateien oder Ausdrucke aus diesem Angebot könnenzusammen mit diesen Nutzungsbedingungen und unter deren Einhaltung weitergegeben werden.Die Speicherung von Teilen des elektronischen Angebots auf anderen Servern ist nur mit vorherigerschriftlicher Genehmigung des Konsortiums der Schweizer Hochschulbibliotheken möglich. DieRechte für diese und andere Nutzungsarten der Inhalte liegen beim Herausgeber bzw. beim Verlag.SEALSEin Dienst des Konsortiums der Schweizer Hochschulbibliothekenc/o ETH-Bibliothek, Rämistrasse 101, 8092 Zürich, Schweizretro@seals.chhttp://retro.seals.chIABSE AIPC IVBHSEMINAR on:«CONCRETE STRUCTURES SUBJECTED TO TRIAXIAL STRESSES»17th-19th MAY,1974-ISMES-BERGAMO(ITALY)III-lConstitutive Model for the Triaxial Behaviour of Concrete Stoffmodell für das mehrachsiale Verhalten von BetonModlle de Constitution pour le Comportement Triaxial du BitonK.J.WILLAMPh. D.,Project LcaderInstitut für Statik und Dynamik der Luft-und Raumfahrtkonstruktionen University of Stuttgart E.P.WARNKEDipl.-Ing.,Research Associate Institut für Statik und Dynamik der Luft-und Raumfahrtkonstruktionen University of StuttgartSUMMARYThis paper describes different modeis for the failure surface and the constitutive behaviourof concrete under triaxial conditions.The study serves two objectives,the working stressdesign and the ultimate load analysis of three-dimensional concrete components.In the first part a three parameter failure surface is developed for concrete subjected to triaxial loading in the tension and low compression regime.This model is subsequently refined by adding two additional parameters for describing curved meridians,thus extend ing the ränge of appli¬cation to the high compression zone.In the second part two constitutive modeis are formulated for elastic perfectly plastic be¬haviour in compression and elastic perfectly brittle behaviour in tension.Based on the normality principle,explicit ex pressions are developed for the inelastic deformation rate and the correspon¬ding incremental stress-strain relation Thus these modeis can be readily applied to ultimate load analysis using the initial load technique or the tangential stiffness method.Dedicated to the60th birthday of Professor Dr.Drs.h.c.J.H.Argyris.1.INTRODUCTIONOver the last two decades a profound change has taken place with the appearance of digital Computers and recent advances in structural analysis[l],[2],[3].The close symbiosis between Computers and structural theories was instrumental for the development of large scale finite element Software packages[4]which found a wide ränge of application in many fields of eng i nee ring sciences.The high degree of sophistication in structural analysis has clearly left behind many other disciplines,one of them being the field of material science.The proper description of the relevant constitutive phenomena has posed a major limitation on the analysis when applied to complex ope¬rating conditions.In the following a constitutive model is presented for the over load and ultimate load ana¬lysis of three-dimensional concrete structures, e.g.Prestressed Concrete Reactor Vessels and Con¬crete Dams.Considering the size of finite elements in a typical idealization one is clearly deal ing with material behaviour on the continuum level,in which the micro structure of piain and rein¬forced concrete components can be neglected.This scale effect of the analysis allows a macro-scopic point of view according to which material phenomena such as cracking can be simulated bythe behaviour of an equivalent continuum.The objective of this study is twofold:First a mathematical model is developed for the description of initial concrete failure under triaxial conditions.Subsequently,this formulation is applied to construct a constitutive model for the over load and ultimate load analysis of three-dimensional concrete Structures.Alternative ly,the failure surface can be applied to working stress design using relevant safety philosophies.In the first part a three parameter model is developed which defines a conical failure sur¬face with non-circular base section in the principal stress space,thus the strength depends on the hydrostatic as well as deviatoric stress state.The proposed failure surface is convex,continuousand has continuous gradient directions furnishing a close fit of test data in the low compression ränge. In the tension regime the model may be augmented by a tension cut-off criterion.This basic formu¬lation is refined in AppendixII by a five parameter model with curved meridians which provides a close fit of test data also in the high compression regime.Subsequently,a material model is constructed based on an elastic perfectly plastic formu¬lation which is augmented by a brittle failure condition in the tensile regime.In this context equivalent constitutive constraint conditions are developed,based on the"normality"principle, which can be readily applied to the finite element analysis via the concept of initial loads.In the past considerable experimental evidence has been gathered which could be used for the construction of a triaxial failure envelope of concrete.However,most of the data were ob-tained from tests with proportional loading and uniform stress or strain conditions which were distorted by unknown boundary layer effects.For the ultimate load analysis via finite elements these two assumptions are clearly invalid.The non-linearity is responsible for local unloading even if the structure is subjected to monotonically increasing stresses.Moreover,the action of a curved thick-walled structure is control led by non-uniform stress distributions,even if global bending effects and local stress concentrations are neglected for the time being.However,for obvious reasons it is customary to assume that test results from uniform stress-or strain experiments can be used to predict the failure behaviour of structural components subjected to non-uniform stress or strain conditions. One should be aware that this fundamental hypothesis has little justification,except that it is at present the only realistic approach for construct ing a phenomenological constitutive law.The actual mechanism of crack initiation and crack propagation could in fact differ fundamental ly between uniform and non-uniform stress distributions.Considerable test data has accumulated on the multiaxial failure behaviour of mortar and concrete specimens subjected to short term loads.The experimental results can be classified into tests in which either two or three stress components are varied independently.To the first category belong the classical triaxial compression tests on cylindrical specimens(triaxial cell experiments) [5],[6],[7],[8],[9]and the biaxial tension-compression tests on hollow cylinders[lO],[l1J In addition,there is the class of biaxial compression and tension-compression tests on slabs [12]/D3L t14L L15L L16l E17"L t18L[191-The second category contains experi¬ments in which cubic specimens are subjected to arbitrary load combinations[20J,[21J Some of these types of tests are present ly still being processed[22],[23],[24JSo far few attempts have been made to utilize this experimental evidence for construct inga mathematical model of the triaxial failure behaviour of concrete.A comprehensive study of this problem was undertaken in[25],for which similar conclusions were reached in[26J,[27].All three modeis fall into the class of pyramidal failure envelopes which have been examined extensive ly within the context of brittle material modeis as general izations of the Mohr-Coulomb criterion[28].In the same publication different modifications of the Griffith criterion are discussed,which have also been applied in[20]to model the failue surface of cubic mortar specimens in the tension-compression regime.None of these previous studies on failure envelopes was directed towards the non-linear analysis of concrete structures.To this end a number of rather simple material formulations were reviewed in[29],[30],[31]and applied to the ultimate load analysis of different concrete structures.TRIAXIAL FAILURE SURFACEIn the following a mathematical model is developed for the triaxial failure surface of con¬crete type materials.Assuming isotropic behaviour the initial failure envelope is fully describedin the principal stress space.Figure1shows the triaxial envelope of concrete type materials.The failure surface is basically a cone with curved meridians and a non-circular base section.The limited tension capacity is responsible for the tetrahedral shape in the tensile regime,while in compression a cylindrical form is ultimate ly reached.For the mathematical model only a sextant of the principal stress space has to be considered, if the stress components are ordered according to S,>Ct>Gs The surface is conveniently represented by hydrostatic and deviatoric sections where the first one forms a meridianal plane which contains the equisectrix S.«6,.»Gh as an axis of revolution The deviatoric section lies in a plane normal to the equisectrix,the deviatoric trace being described by the polar coordinatesr,ösee Fig. 2.Basically,there are four aspects to the mathematical model of the failure surface:1Close fit of experimental data in the operating ränge.2.Simple identification of model parameters from Standard test data.3.Smoothness-continuous surface with continuously varying tangent planes.4.Convexity-monotonically curved surface without inflection points.Close approximation of concrete data is reached if the failure surface depends on the hydrostatic as well as the deviatoric state,whereby the latter should distinguish different strength values according to the direction of deviatoric stress.Therefore,the failure envelope must be basically a conical surface with curved meridians and a non-circular base section.In addition,in the tensile regime the failure suface could be augmented by a tension cut-off criterion in the form of a pyramid with triangulär section in the deviatoric plane.Simple identification means that the mathematical model of the failure surface is definedby a very small number of parameters which can be determined from Standard test data, e.g. uniaxial tension,uniaxial compression,biaxial compression tests,etc.The description of the failure surface should also encompass simple failure envelopes for specific model parameters.In other words,the cylindrical von Mises and the conical Drucker-Prager model should be special cases of the sophisticated failure formula tion.Continuity is an important property for two reasons:From a computational point of view,it is very convenient if a single description of the failure surface is valid within the stress space under consideration.From the theoretical point of view the proposed failure surface should havea unique gradient for defining the direction of the inelastic deformations according to the1normolity principle1.The actual nature of concrete failure mechanisms also supports the conceptof a gradual change of strength for small variations in loading.Geometrical ly,the smoothness condition implies that the failure surface is continuous and has continuous derivatives.Therefore,the deviatoric trace of the failure surface must pass through r,and ru with the tangents4:,and tr at9-O*and0-CO°,see Fig. 2.Recall thatfor Isotropie conditions only a sextant of the stress space has to be considered,O^O*60*.Convexity is an important property since it assures stable matenal behaviour according to the postulate of Drucker[32],if the"normolity"principle determi es the direction of inelastic deformations.Stability infers positive dissipation of inelastic work during a loading cycle according to the coneepts of thermodynamics Figure3indicates that convexity of the overall deviatoric trace can be assured only if there are no inflection points end if the position vector satisfies the basic convexity conditionJL>J_where r;r(0.:*,i2o",z4o#)Ti.*r(ö*£3,iSo^soo)W Continuity infers compatibility of the position vectors and the slopes et0«O*and0¦GOD. Consequently,there are at least four conditions for curve-fitting the deviatoric trace withinO**0*GOt In addition,the convexity condition implies that the curve should have no inflection points in this interval,thus the approximation can not be based on trigonometric functions[30]or Hermitian Interpolation.If the curve should also degenerate to a circle for r,-r^then an elliptic approximation has to be used for the functional Variation of the deviatoric trace.The ellipsoidal surface assures smoothness and convexity for all position vectors r satisfy ing|rl4r<rL(2)The geometric construction of the ellipse is shown in Fig.4,the details of the derivation are given in the Appendix I.The half axes of the ellipse a,b are defined in terms of the position vectorsErt-4.r,(3)a.r,1'-Sr.ty*Lr%¦4.r,-BrtThe elliptic trace is expressed in terms of the polar coordinates r,B by(r>-r.*)cos9*rt(ir.-r«)&(tf-tf)o»9?Sif-Ar.rj(4a) rf8x1ft,with the angle of similarity9^,^gt"2-g3>(4b) aIn the following the deviatoric trace is used as base section of a conical failure surface with the equisectrix as axis of revolution.A linear Variation with hydrostatic stress generates a cone with straight line meridians.In this case the failure surface is defined in principal stress space by a homogeneous expansion in the"average11stress components$>Ä«:«.and the angle of similarity9IdW^MJ.i^^J-i(5)The average stress components6«,t«.represent the mean distribution of normal and shear stresses on an infinitesimal spherical surface.These values are normolized in the failure condition eq.(5)by the uniaxial compressive strength f^The stress components are defined in terms of principal stresses by(6)^-jL[(».-*£-(*-*£>+t*»-**]*These scalar representations of the state of stress at a point are related to the stress components on the"octahedral"plane60^o by**s*o(7a)**«ff t0The average stress components also correspond to the first principal stress invariant I,and the second deviatoric stress invariant I%according to**-T*.(7b)».-nfi^-RW«^For material failure,-f(>)*0the following constraint condition must hold between the average normal stress and the average shear stresst-rt^D-Ht]i«.1(8)The free parameters of the failure surface model fc H and f«.are identified below from typical concrete test data,such as the uniaxial tension test«f t the uniaxial compression testf4a and the biaxial compression test^Introducing the strength ratios ocft,«c0*z*ft/fco(9) the three tests are characterized by6.TEST Wf«,«/?«.er.*,-ft+--(¥->6Ä»feo3K CO***xei*^"£fc"*««R-O*n(10)Substituting these strength values into the failure condition eq.8,the model parametersare readily obtainedn<*o Af«*»-**«Tu**£«r+ä**oätjt*«««?<*.-«z(11)The apex of the conical surface lies on the equisectrix at*»«*The opening angle<P of the cone varies betweenandtan<y,«-^«+0-GO°(12)(13) The proposed three parameter model is illustrated in Fig.5for the strength ratios0(O'and42*o.i The hydrostatic and deviatoric sections indicate the convexity and smoothness failure envelope.The proposed failure surface degenerates to the Drucker-Prager model of a circular cone ifIn this case the conical failure surface is described by the two parameters2and r0±_**^-^-1i.i of the(14)(15)The single parameter von Mises model is obtained/if in addition7.2.-t>oo(16)[21].In this case the Drucker-Prager cone degenerates into a circular cylinder whose radius is defined byr.fo (17)with the strength ratios**m **"¦(18)Figure 6shows a comparison between the failure surface and experimental data reported inClose agreement can be observed in the low pressure regime for the strength ratios *oc'-Äand *i«o.i?In the high compression regime there is considerable disagreement mainly along the compressive branch.Therefore,the three parameter model is refined in the Appendix II by two additional parameters,extend ing the ränge of application to the high compression regime.This five parameter model establishes a failure surface with curved meridians in which the generators are approximated by second order parabolas along 0s O*and 0*Go with a common apex at the equisectrix,see also Fig.11Figure 7shows the biaxial failure envelope of the three parameter model for three differentstrength ratios otu*l&*^»o-U &u~)o,**.*o.o%and *u-l-&j **-o.is A comparison with test data from [18]^2l]indicates that the shear strength is overestimated consi-derably because of the acute intersection with the biaxial stress plane.However,if we consider the dominant influence of the post-failure behaviour on the structural response [30J,there is little reason for further refinements of the initial failure surface model.3.CONSTITUTIVE MODELIn the following the previous model of the failure envelope is utilized for the developmentof an elastic perfectly plastic material formula tion in compression.The constitutive model is sub¬sequently augmented by a tension cut-off criterion to account for cracking in the tension regime.In both cases it is assumed that the normo lity principle determines the direction of the inelastic deformation rates for ductile as well as brittle post failure behaviour.3.1Elastic Plastic FormulationInviscid plasticity is the classical approach for describing inelastic behaviour via incremen-tal stress-strain relations.The constitutive model is based on two fundamental assumptions,an appropriate description of the material failure envelope and the definition of inelastic deformation rates e.g.via the normo lity principle.a.Yield ConditionThe yield surface serves two objectives,it distinguishes linear from non-linear andelastic from inelastic deformations.The failure envelope is defined by a scalar function of stress,$(J5»o indicating plastic flow if the stress path intersects the yield surface.For concrete type of materials the yield condition can be approximated by the three parameter model shown in Fig.5or more accurately by the five parameter model developed in the Appendix II.b.Flow RuleFor perfectly plastic behaviour the yield surface does not change its configurationduring plastic flow,hence the stress path describes a trajectory on the initial yield surface,whilethe inelastic strains increase continuously.In this case the inelastic deformations do not contri-bute to the elastic strain energy,thus the inner product of plastic strain and elastic stressrates must be zeron**-°09)-In other words,the plastic strain rate must be perpendicular to the yield surfaceV^(2°) where the normal n is the unit gradient vector of the yield surface9j/9*(21)mExplicit expressions of2f/d9are developed in Appendix III for different yield surfaces.The normal defines the direction of the plastic strain rate,the length of which determines the loading parameter a The normality condition follows from Drucker's stability postulate which assures non-negative work dissipation during a loading cycle,also infernng convexity of the yield surface.For perfectly plastic behaviour the material stability is"indifferent"in the small, corresponding to the"neutral"loading condition for which initial yield and subsequent flow is governed by*«>-0and fC*>-o(22) The consistency condition implies that?(«>-Tt*"°(23) This Statement is clearly equivalent to the normality principle stated in eq.(19).c.Incremental Stress-Strain RelationsIn the following an elastic perfectly plastic consitutive model is derived using the previous Statements and the kinematic decomposition of the total deformationsy-C+Ylf,and TF"*+V(24) The linear elastic material behaviour is given by the rate formulation of generalized Hooke's law*Ei-£(i-^(25) Substituting the stress rate into the consistency condition,eq.(23),we obtain£«-n4E(JHW(26) This expression yields for the undetermined loading parameter ahtE(T-n\)-o(27)and hereby±i'WTT"E*(28) x The dot indicates the rate of change.9. The plastic strain rate follows from eq.(20)(29) The incremental stress-strain relations are obtained by Substitut ing y*v into the expression ofthe stress rate,eq.(25)Note the linear relationship between the stress and deformation rates in eq.(30)&F Y od The tangential material law Tr is defined byFor perfectly plastic behaviour,T depends only on the elastic properties and the instantaneous stress state via II The second term of eq.(32)represents the degradation of the material Constitution due to plastic flow.3.2Elastic Cracking FormulationSmall tensile strength is the predominant feature of concrete-type materials.In the following a simple constitutive model is developed for perfectly brittle behaviour in the tensile regime.In analogy to the elastic plastic formulation the elastic cracking model is based on two fundamental assumptions,a tension cut-off criterion for the prediction of cracking and an appro-priate description of inelastic deformation rates e.g.via the normality principle.a.Crack ConditionThe tension cut-off criterion distinguishes elastic behaviour from brittle fracture,i.e. Separation of the material constituents due to excess tension.To this end it is assumed that the scale of Observation justifies a continuum approach.For concrete-type materials cracking may be predicted by the single one parameter model based on the major principal stress^t**)<5-.-^e wifh*>i*^i>^i(33) where C\corresponds in general to the uniaxial tensile strength ft The failure surface is shown in Fig.8,which indicates the pyramidal shape and the triangulär base section in the deviatoric plane.Alternative ly,the tension cut-off condition could also be expressed in terms of the three parameter model of the previous section or the five parameter model developed inthe Appendix II.b.Fracture RuleFor ductile behaviour in the post failure ränge the inelastic deformation rate due to cracking is derived exactly along the formulation of an elastic plastic solid.The ductilepost failure behaviour forms an upper bound of the actual soften ing behaviour[30J,which may develop in concrete components due to reinforcements,dowel action and aggregate interlock.In the following/the case of perfectly brittle post-failure behaviour is discussed,since it requires slight modifications of the previous constitutive model for an elastic perfectly plastic solid. In analogy to elasto-plasticity the inelastic deformations due to cracking tlc do not contribute to the elastic strain energytYI.5-°/iü\10.This normality principle corresponds to the flow rule of plasticity stating that the inelastic strain rate due to cracking is perpendicular to the plane of fractureT|c-nX(35) For the maximum stress tension cut-off criterion the normal vector f\is defined by the direction of the major principal stress;thus in the principal stress spacem »f/»««(36)n|9$/9«M^iwhere£,is the unit vector«,-[l,o,o,o,o.°i(37) For perfectly brittle behaviour the loading parameter X is determined from the sofrening condition|C«}-0and£C«^*-*t(360 In this case the consistency condition infers that|t**t(39)c.Inelastic Strain IncrementsIn the following an expression is derived for the inelastic deformation rotes due to cracking.Substituting the stress rate expression eq.(25)into the consistency condition eq.(39)&*-**(*-M w we obtain an expression for the undetermined loading parameter A«*6(*-«.X)--«e.(") and hereby*-«Fg-*,(*,*T+<)(«)Note the equivalence to the elastic plastic formulation in eq.(28)except for the release of^due to brittle softening.The resulting inelastic fracture strain rate follows from eq.(35)(43)The first portion of this expression can be used to construct incremental stress-strain rela-tions in analogy to the elastic plastic formulation,see eq.(30).This part would correspondexactly to a ductile cracking model in which the major principal stress is kept constant at thetensile strength S,^e The corresponding tangential material law would become transversely isotropic with zero stiffness along the major principal axis.Additional cracking in other directions can be considered according ly.The second portion of eq.(43)represents the sudden stress release due to brittle fracture,*\, which is projeeted onto the structural level by a single initial load step in the analysis.11.3.3Transition ProblemThe previous rate formulation for elastic plastic and brittle behaviour is valid in a diffe¬rential sense only.In a numericaj environment clearly finite increments prevail during numerical Integration of the rate equations[33J,[34J.This approximation problem is magnified by the sudden transition from elastic to plastic or elastic to brittle behaviour.In the latter case the dis-continuity of the process is further increased due to the immediate stress release if the failure condition has been reached.Clearly,the success of the numerical technique depends primärily on the proper treatment of the transition problem for finite increments.Consider the most general case of a finite load step shown in Fig.9.At the outset we assume that the stress path has reached point A for which^C«*V°indicates an elastic state. Due to the finite load increment a fully elastic stress path would reach point B penetrating the yield surface at C for proportional loading.The condition$C^O>°violates Tne constitutive constraint condition^«^°(45) and suggests two strategies for numerical implementation.a.^P^^^J.^6.^^*J£n MethodAssuming proportional loading the load increment is subdivided into two parts,an elastic portion for the path A-C and an inelastic portion governing the behaviour after the failure surface has been reached at C.The evaluation of the penetration point C reduces to the geometric problem of intersecting a surface with a line,a task which is non-linear for curved failure envelopes.The computation of the stress trajectory on the yield surface involves the numerical integration of5^1F if(46)since the tangential material law varies with the current state of stress.In addition we have to assume that the inelastic strains increase proportional ly from ycto<jf£.In numerical calculations additional corrections are required at each iteration step to place the stress path back onto the yield surface[37]b.Normal Penetration MethodIn this scheme we assume that the elastic path reaches the yield surface at the inter-section with the normal ns The evaluation of the foot point D reduces to the geometric-prob¬lem of minimizing the distance between B and the failure envelope,see Fig.9<A-(«^«^(«^«J"*Minimum(47)The extremum condition is used to determine the components of C^by solving the linear system of equations.subjected to the constraint conditionft«*)*ö(49)Note that the loading parameter X is proportional to the distance d,thus the length of the inelastic deformation increment is determined from。

Draft:Deep Learning in Neural Networks:An OverviewTechnical Report IDSIA-03-14/arXiv:1404.7828(v1.5)[cs.NE]J¨u rgen SchmidhuberThe Swiss AI Lab IDSIAIstituto Dalle Molle di Studi sull’Intelligenza ArtificialeUniversity of Lugano&SUPSIGalleria2,6928Manno-LuganoSwitzerland15May2014AbstractIn recent years,deep artificial neural networks(including recurrent ones)have won numerous con-tests in pattern recognition and machine learning.This historical survey compactly summarises relevantwork,much of it from the previous millennium.Shallow and deep learners are distinguished by thedepth of their credit assignment paths,which are chains of possibly learnable,causal links between ac-tions and effects.I review deep supervised learning(also recapitulating the history of backpropagation),unsupervised learning,reinforcement learning&evolutionary computation,and indirect search for shortprograms encoding deep and large networks.PDF of earlier draft(v1):http://www.idsia.ch/∼juergen/DeepLearning30April2014.pdfLATEX source:http://www.idsia.ch/∼juergen/DeepLearning30April2014.texComplete BIBTEXfile:http://www.idsia.ch/∼juergen/bib.bibPrefaceThis is the draft of an invited Deep Learning(DL)overview.One of its goals is to assign credit to those who contributed to the present state of the art.I acknowledge the limitations of attempting to achieve this goal.The DL research community itself may be viewed as a continually evolving,deep network of scientists who have influenced each other in complex ways.Starting from recent DL results,I tried to trace back the origins of relevant ideas through the past half century and beyond,sometimes using“local search”to follow citations of citations backwards in time.Since not all DL publications properly acknowledge earlier relevant work,additional global search strategies were employed,aided by consulting numerous neural network experts.As a result,the present draft mostly consists of references(about800entries so far).Nevertheless,through an expert selection bias I may have missed important work.A related bias was surely introduced by my special familiarity with the work of my own DL research group in the past quarter-century.For these reasons,the present draft should be viewed as merely a snapshot of an ongoing credit assignment process.To help improve it,please do not hesitate to send corrections and suggestions to juergen@idsia.ch.Contents1Introduction to Deep Learning(DL)in Neural Networks(NNs)3 2Event-Oriented Notation for Activation Spreading in FNNs/RNNs3 3Depth of Credit Assignment Paths(CAPs)and of Problems4 4Recurring Themes of Deep Learning54.1Dynamic Programming(DP)for DL (5)4.2Unsupervised Learning(UL)Facilitating Supervised Learning(SL)and RL (6)4.3Occam’s Razor:Compression and Minimum Description Length(MDL) (6)4.4Learning Hierarchical Representations Through Deep SL,UL,RL (6)4.5Fast Graphics Processing Units(GPUs)for DL in NNs (6)5Supervised NNs,Some Helped by Unsupervised NNs75.11940s and Earlier (7)5.2Around1960:More Neurobiological Inspiration for DL (7)5.31965:Deep Networks Based on the Group Method of Data Handling(GMDH) (8)5.41979:Convolution+Weight Replication+Winner-Take-All(WTA) (8)5.51960-1981and Beyond:Development of Backpropagation(BP)for NNs (8)5.5.1BP for Weight-Sharing Feedforward NNs(FNNs)and Recurrent NNs(RNNs)..95.6Late1980s-2000:Numerous Improvements of NNs (9)5.6.1Ideas for Dealing with Long Time Lags and Deep CAPs (10)5.6.2Better BP Through Advanced Gradient Descent (10)5.6.3Discovering Low-Complexity,Problem-Solving NNs (11)5.6.4Potential Benefits of UL for SL (11)5.71987:UL Through Autoencoder(AE)Hierarchies (12)5.81989:BP for Convolutional NNs(CNNs) (13)5.91991:Fundamental Deep Learning Problem of Gradient Descent (13)5.101991:UL-Based History Compression Through a Deep Hierarchy of RNNs (14)5.111992:Max-Pooling(MP):Towards MPCNNs (14)5.121994:Contest-Winning Not So Deep NNs (15)5.131995:Supervised Recurrent Very Deep Learner(LSTM RNN) (15)5.142003:More Contest-Winning/Record-Setting,Often Not So Deep NNs (16)5.152006/7:Deep Belief Networks(DBNs)&AE Stacks Fine-Tuned by BP (17)5.162006/7:Improved CNNs/GPU-CNNs/BP-Trained MPCNNs (17)5.172009:First Official Competitions Won by RNNs,and with MPCNNs (18)5.182010:Plain Backprop(+Distortions)on GPU Yields Excellent Results (18)5.192011:MPCNNs on GPU Achieve Superhuman Vision Performance (18)5.202011:Hessian-Free Optimization for RNNs (19)5.212012:First Contests Won on ImageNet&Object Detection&Segmentation (19)5.222013-:More Contests and Benchmark Records (20)5.22.1Currently Successful Supervised Techniques:LSTM RNNs/GPU-MPCNNs (21)5.23Recent Tricks for Improving SL Deep NNs(Compare Sec.5.6.2,5.6.3) (21)5.24Consequences for Neuroscience (22)5.25DL with Spiking Neurons? (22)6DL in FNNs and RNNs for Reinforcement Learning(RL)236.1RL Through NN World Models Yields RNNs With Deep CAPs (23)6.2Deep FNNs for Traditional RL and Markov Decision Processes(MDPs) (24)6.3Deep RL RNNs for Partially Observable MDPs(POMDPs) (24)6.4RL Facilitated by Deep UL in FNNs and RNNs (25)6.5Deep Hierarchical RL(HRL)and Subgoal Learning with FNNs and RNNs (25)6.6Deep RL by Direct NN Search/Policy Gradients/Evolution (25)6.7Deep RL by Indirect Policy Search/Compressed NN Search (26)6.8Universal RL (27)7Conclusion271Introduction to Deep Learning(DL)in Neural Networks(NNs) Which modifiable components of a learning system are responsible for its success or failure?What changes to them improve performance?This has been called the fundamental credit assignment problem(Minsky, 1963).There are general credit assignment methods for universal problem solvers that are time-optimal in various theoretical senses(Sec.6.8).The present survey,however,will focus on the narrower,but now commercially important,subfield of Deep Learning(DL)in Artificial Neural Networks(NNs).We are interested in accurate credit assignment across possibly many,often nonlinear,computational stages of NNs.Shallow NN-like models have been around for many decades if not centuries(Sec.5.1).Models with several successive nonlinear layers of neurons date back at least to the1960s(Sec.5.3)and1970s(Sec.5.5). An efficient gradient descent method for teacher-based Supervised Learning(SL)in discrete,differentiable networks of arbitrary depth called backpropagation(BP)was developed in the1960s and1970s,and ap-plied to NNs in1981(Sec.5.5).BP-based training of deep NNs with many layers,however,had been found to be difficult in practice by the late1980s(Sec.5.6),and had become an explicit research subject by the early1990s(Sec.5.9).DL became practically feasible to some extent through the help of Unsupervised Learning(UL)(e.g.,Sec.5.10,5.15).The1990s and2000s also saw many improvements of purely super-vised DL(Sec.5).In the new millennium,deep NNs havefinally attracted wide-spread attention,mainly by outperforming alternative machine learning methods such as kernel machines(Vapnik,1995;Sch¨o lkopf et al.,1998)in numerous important applications.In fact,supervised deep NNs have won numerous of-ficial international pattern recognition competitions(e.g.,Sec.5.17,5.19,5.21,5.22),achieving thefirst superhuman visual pattern recognition results in limited domains(Sec.5.19).Deep NNs also have become relevant for the more generalfield of Reinforcement Learning(RL)where there is no supervising teacher (Sec.6).Both feedforward(acyclic)NNs(FNNs)and recurrent(cyclic)NNs(RNNs)have won contests(Sec.5.12,5.14,5.17,5.19,5.21,5.22).In a sense,RNNs are the deepest of all NNs(Sec.3)—they are general computers more powerful than FNNs,and can in principle create and process memories of ar-bitrary sequences of input patterns(e.g.,Siegelmann and Sontag,1991;Schmidhuber,1990a).Unlike traditional methods for automatic sequential program synthesis(e.g.,Waldinger and Lee,1969;Balzer, 1985;Soloway,1986;Deville and Lau,1994),RNNs can learn programs that mix sequential and parallel information processing in a natural and efficient way,exploiting the massive parallelism viewed as crucial for sustaining the rapid decline of computation cost observed over the past75years.The rest of this paper is structured as follows.Sec.2introduces a compact,event-oriented notation that is simple yet general enough to accommodate both FNNs and RNNs.Sec.3introduces the concept of Credit Assignment Paths(CAPs)to measure whether learning in a given NN application is of the deep or shallow type.Sec.4lists recurring themes of DL in SL,UL,and RL.Sec.5focuses on SL and UL,and on how UL can facilitate SL,although pure SL has become dominant in recent competitions(Sec.5.17-5.22). Sec.5is arranged in a historical timeline format with subsections on important inspirations and technical contributions.Sec.6on deep RL discusses traditional Dynamic Programming(DP)-based RL combined with gradient-based search techniques for SL or UL in deep NNs,as well as general methods for direct and indirect search in the weight space of deep FNNs and RNNs,including successful policy gradient and evolutionary methods.2Event-Oriented Notation for Activation Spreading in FNNs/RNNs Throughout this paper,let i,j,k,t,p,q,r denote positive integer variables assuming ranges implicit in the given contexts.Let n,m,T denote positive integer constants.An NN’s topology may change over time(e.g.,Fahlman,1991;Ring,1991;Weng et al.,1992;Fritzke, 1994).At any given moment,it can be described as afinite subset of units(or nodes or neurons)N= {u1,u2,...,}and afinite set H⊆N×N of directed edges or connections between nodes.FNNs are acyclic graphs,RNNs cyclic.Thefirst(input)layer is the set of input units,a subset of N.In FNNs,the k-th layer(k>1)is the set of all nodes u∈N such that there is an edge path of length k−1(but no longer path)between some input unit and u.There may be shortcut connections between distant layers.The NN’s behavior or program is determined by a set of real-valued,possibly modifiable,parameters or weights w i(i=1,...,n).We now focus on a singlefinite episode or epoch of information processing and activation spreading,without learning through weight changes.The following slightly unconventional notation is designed to compactly describe what is happening during the runtime of the system.During an episode,there is a partially causal sequence x t(t=1,...,T)of real values that I call events.Each x t is either an input set by the environment,or the activation of a unit that may directly depend on other x k(k<t)through a current NN topology-dependent set in t of indices k representing incoming causal connections or links.Let the function v encode topology information and map such event index pairs(k,t)to weight indices.For example,in the non-input case we may have x t=f t(net t)with real-valued net t= k∈in t x k w v(k,t)(additive case)or net t= k∈in t x k w v(k,t)(multiplicative case), where f t is a typically nonlinear real-valued activation function such as tanh.In many recent competition-winning NNs(Sec.5.19,5.21,5.22)there also are events of the type x t=max k∈int (x k);some networktypes may also use complex polynomial activation functions(Sec.5.3).x t may directly affect certain x k(k>t)through outgoing connections or links represented through a current set out t of indices k with t∈in k.Some non-input events are called output events.Note that many of the x t may refer to different,time-varying activations of the same unit in sequence-processing RNNs(e.g.,Williams,1989,“unfolding in time”),or also in FNNs sequentially exposed to time-varying input patterns of a large training set encoded as input events.During an episode,the same weight may get reused over and over again in topology-dependent ways,e.g.,in RNNs,or in convolutional NNs(Sec.5.4,5.8).I call this weight sharing across space and/or time.Weight sharing may greatly reduce the NN’s descriptive complexity,which is the number of bits of information required to describe the NN (Sec.4.3).In Supervised Learning(SL),certain NN output events x t may be associated with teacher-given,real-valued labels or targets d t yielding errors e t,e.g.,e t=1/2(x t−d t)2.A typical goal of supervised NN training is tofind weights that yield episodes with small total error E,the sum of all such e t.The hope is that the NN will generalize well in later episodes,causing only small errors on previously unseen sequences of input events.Many alternative error functions for SL and UL are possible.SL assumes that input events are independent of earlier output events(which may affect the environ-ment through actions causing subsequent perceptions).This assumption does not hold in the broaderfields of Sequential Decision Making and Reinforcement Learning(RL)(Kaelbling et al.,1996;Sutton and Barto, 1998;Hutter,2005)(Sec.6).In RL,some of the input events may encode real-valued reward signals given by the environment,and a typical goal is tofind weights that yield episodes with a high sum of reward signals,through sequences of appropriate output actions.Sec.5.5will use the notation above to compactly describe a central algorithm of DL,namely,back-propagation(BP)for supervised weight-sharing FNNs and RNNs.(FNNs may be viewed as RNNs with certainfixed zero weights.)Sec.6will address the more general RL case.3Depth of Credit Assignment Paths(CAPs)and of ProblemsTo measure whether credit assignment in a given NN application is of the deep or shallow type,I introduce the concept of Credit Assignment Paths or CAPs,which are chains of possibly causal links between events.Let usfirst focus on SL.Consider two events x p and x q(1≤p<q≤T).Depending on the appli-cation,they may have a Potential Direct Causal Connection(PDCC)expressed by the Boolean predicate pdcc(p,q),which is true if and only if p∈in q.Then the2-element list(p,q)is defined to be a CAP from p to q(a minimal one).A learning algorithm may be allowed to change w v(p,q)to improve performance in future episodes.More general,possibly indirect,Potential Causal Connections(PCC)are expressed by the recursively defined Boolean predicate pcc(p,q),which in the SL case is true only if pdcc(p,q),or if pcc(p,k)for some k and pdcc(k,q).In the latter case,appending q to any CAP from p to k yields a CAP from p to q(this is a recursive definition,too).The set of such CAPs may be large but isfinite.Note that the same weight may affect many different PDCCs between successive events listed by a given CAP,e.g.,in the case of RNNs, or weight-sharing FNNs.Suppose a CAP has the form(...,k,t,...,q),where k and t(possibly t=q)are thefirst successive elements with modifiable w v(k,t).Then the length of the suffix list(t,...,q)is called the CAP’s depth (which is0if there are no modifiable links at all).This depth limits how far backwards credit assignment can move down the causal chain tofind a modifiable weight.1Suppose an episode and its event sequence x1,...,x T satisfy a computable criterion used to decide whether a given problem has been solved(e.g.,total error E below some threshold).Then the set of used weights is called a solution to the problem,and the depth of the deepest CAP within the sequence is called the solution’s depth.There may be other solutions(yielding different event sequences)with different depths.Given somefixed NN topology,the smallest depth of any solution is called the problem’s depth.Sometimes we also speak of the depth of an architecture:SL FNNs withfixed topology imply a problem-independent maximal problem depth bounded by the number of non-input layers.Certain SL RNNs withfixed weights for all connections except those to output units(Jaeger,2001;Maass et al.,2002; Jaeger,2004;Schrauwen et al.,2007)have a maximal problem depth of1,because only thefinal links in the corresponding CAPs are modifiable.In general,however,RNNs may learn to solve problems of potentially unlimited depth.Note that the definitions above are solely based on the depths of causal chains,and agnostic of the temporal distance between events.For example,shallow FNNs perceiving large“time windows”of in-put events may correctly classify long input sequences through appropriate output events,and thus solve shallow problems involving long time lags between relevant events.At which problem depth does Shallow Learning end,and Deep Learning begin?Discussions with DL experts have not yet yielded a conclusive response to this question.Instead of committing myself to a precise answer,let me just define for the purposes of this overview:problems of depth>10require Very Deep Learning.The difficulty of a problem may have little to do with its depth.Some NNs can quickly learn to solve certain deep problems,e.g.,through random weight guessing(Sec.5.9)or other types of direct search (Sec.6.6)or indirect search(Sec.6.7)in weight space,or through training an NNfirst on shallow problems whose solutions may then generalize to deep problems,or through collapsing sequences of(non)linear operations into a single(non)linear operation—but see an analysis of non-trivial aspects of deep linear networks(Baldi and Hornik,1994,Section B).In general,however,finding an NN that precisely models a given training set is an NP-complete problem(Judd,1990;Blum and Rivest,1992),also in the case of deep NNs(S´ıma,1994;de Souto et al.,1999;Windisch,2005);compare a survey of negative results(S´ıma, 2002,Section1).Above we have focused on SL.In the more general case of RL in unknown environments,pcc(p,q) is also true if x p is an output event and x q any later input event—any action may affect the environment and thus any later perception.(In the real world,the environment may even influence non-input events computed on a physical hardware entangled with the entire universe,but this is ignored here.)It is possible to model and replace such unmodifiable environmental PCCs through a part of the NN that has already learned to predict(through some of its units)input events(including reward signals)from former input events and actions(Sec.6.1).Its weights are frozen,but can help to assign credit to other,still modifiable weights used to compute actions(Sec.6.1).This approach may lead to very deep CAPs though.Some DL research is about automatically rephrasing problems such that their depth is reduced(Sec.4). In particular,sometimes UL is used to make SL problems less deep,e.g.,Sec.5.10.Often Dynamic Programming(Sec.4.1)is used to facilitate certain traditional RL problems,e.g.,Sec.6.2.Sec.5focuses on CAPs for SL,Sec.6on the more complex case of RL.4Recurring Themes of Deep Learning4.1Dynamic Programming(DP)for DLOne recurring theme of DL is Dynamic Programming(DP)(Bellman,1957),which can help to facili-tate credit assignment under certain assumptions.For example,in SL NNs,backpropagation itself can 1An alternative would be to count only modifiable links when measuring depth.In many typical NN applications this would not make a difference,but in some it would,e.g.,Sec.6.1.be viewed as a DP-derived method(Sec.5.5).In traditional RL based on strong Markovian assumptions, DP-derived methods can help to greatly reduce problem depth(Sec.6.2).DP algorithms are also essen-tial for systems that combine concepts of NNs and graphical models,such as Hidden Markov Models (HMMs)(Stratonovich,1960;Baum and Petrie,1966)and Expectation Maximization(EM)(Dempster et al.,1977),e.g.,(Bottou,1991;Bengio,1991;Bourlard and Morgan,1994;Baldi and Chauvin,1996; Jordan and Sejnowski,2001;Bishop,2006;Poon and Domingos,2011;Dahl et al.,2012;Hinton et al., 2012a).4.2Unsupervised Learning(UL)Facilitating Supervised Learning(SL)and RL Another recurring theme is how UL can facilitate both SL(Sec.5)and RL(Sec.6).UL(Sec.5.6.4) is normally used to encode raw incoming data such as video or speech streams in a form that is more convenient for subsequent goal-directed learning.In particular,codes that describe the original data in a less redundant or more compact way can be fed into SL(Sec.5.10,5.15)or RL machines(Sec.6.4),whose search spaces may thus become smaller(and whose CAPs shallower)than those necessary for dealing with the raw data.UL is closely connected to the topics of regularization and compression(Sec.4.3,5.6.3). 4.3Occam’s Razor:Compression and Minimum Description Length(MDL) Occam’s razor favors simple solutions over complex ones.Given some programming language,the prin-ciple of Minimum Description Length(MDL)can be used to measure the complexity of a solution candi-date by the length of the shortest program that computes it(e.g.,Solomonoff,1964;Kolmogorov,1965b; Chaitin,1966;Wallace and Boulton,1968;Levin,1973a;Rissanen,1986;Blumer et al.,1987;Li and Vit´a nyi,1997;Gr¨u nwald et al.,2005).Some methods explicitly take into account program runtime(Al-lender,1992;Watanabe,1992;Schmidhuber,2002,1995);many consider only programs with constant runtime,written in non-universal programming languages(e.g.,Rissanen,1986;Hinton and van Camp, 1993).In the NN case,the MDL principle suggests that low NN weight complexity corresponds to high NN probability in the Bayesian view(e.g.,MacKay,1992;Buntine and Weigend,1991;De Freitas,2003), and to high generalization performance(e.g.,Baum and Haussler,1989),without overfitting the training data.Many methods have been proposed for regularizing NNs,that is,searching for solution-computing, low-complexity SL NNs(Sec.5.6.3)and RL NNs(Sec.6.7).This is closely related to certain UL methods (Sec.4.2,5.6.4).4.4Learning Hierarchical Representations Through Deep SL,UL,RLMany methods of Good Old-Fashioned Artificial Intelligence(GOFAI)(Nilsson,1980)as well as more recent approaches to AI(Russell et al.,1995)and Machine Learning(Mitchell,1997)learn hierarchies of more and more abstract data representations.For example,certain methods of syntactic pattern recog-nition(Fu,1977)such as grammar induction discover hierarchies of formal rules to model observations. The partially(un)supervised Automated Mathematician/EURISKO(Lenat,1983;Lenat and Brown,1984) continually learns concepts by combining previously learnt concepts.Such hierarchical representation learning(Ring,1994;Bengio et al.,2013;Deng and Yu,2014)is also a recurring theme of DL NNs for SL (Sec.5),UL-aided SL(Sec.5.7,5.10,5.15),and hierarchical RL(Sec.6.5).Often,abstract hierarchical representations are natural by-products of data compression(Sec.4.3),e.g.,Sec.5.10.4.5Fast Graphics Processing Units(GPUs)for DL in NNsWhile the previous millennium saw several attempts at creating fast NN-specific hardware(e.g.,Jackel et al.,1990;Faggin,1992;Ramacher et al.,1993;Widrow et al.,1994;Heemskerk,1995;Korkin et al., 1997;Urlbe,1999),and at exploiting standard hardware(e.g.,Anguita et al.,1994;Muller et al.,1995; Anguita and Gomes,1996),the new millennium brought a DL breakthrough in form of cheap,multi-processor graphics cards or GPUs.GPUs are widely used for video games,a huge and competitive market that has driven down hardware prices.GPUs excel at fast matrix and vector multiplications required not only for convincing virtual realities but also for NN training,where they can speed up learning by a factorof50and more.Some of the GPU-based FNN implementations(Sec.5.16-5.19)have greatly contributed to recent successes in contests for pattern recognition(Sec.5.19-5.22),image segmentation(Sec.5.21), and object detection(Sec.5.21-5.22).5Supervised NNs,Some Helped by Unsupervised NNsThe main focus of current practical applications is on Supervised Learning(SL),which has dominated re-cent pattern recognition contests(Sec.5.17-5.22).Several methods,however,use additional Unsupervised Learning(UL)to facilitate SL(Sec.5.7,5.10,5.15).It does make sense to treat SL and UL in the same section:often gradient-based methods,such as BP(Sec.5.5.1),are used to optimize objective functions of both UL and SL,and the boundary between SL and UL may blur,for example,when it comes to time series prediction and sequence classification,e.g.,Sec.5.10,5.12.A historical timeline format will help to arrange subsections on important inspirations and techni-cal contributions(although such a subsection may span a time interval of many years).Sec.5.1briefly mentions early,shallow NN models since the1940s,Sec.5.2additional early neurobiological inspiration relevant for modern Deep Learning(DL).Sec.5.3is about GMDH networks(since1965),perhaps thefirst (feedforward)DL systems.Sec.5.4is about the relatively deep Neocognitron NN(1979)which is similar to certain modern deep FNN architectures,as it combines convolutional NNs(CNNs),weight pattern repli-cation,and winner-take-all(WTA)mechanisms.Sec.5.5uses the notation of Sec.2to compactly describe a central algorithm of DL,namely,backpropagation(BP)for supervised weight-sharing FNNs and RNNs. It also summarizes the history of BP1960-1981and beyond.Sec.5.6describes problems encountered in the late1980s with BP for deep NNs,and mentions several ideas from the previous millennium to overcome them.Sec.5.7discusses afirst hierarchical stack of coupled UL-based Autoencoders(AEs)—this concept resurfaced in the new millennium(Sec.5.15).Sec.5.8is about applying BP to CNNs,which is important for today’s DL applications.Sec.5.9explains BP’s Fundamental DL Problem(of vanishing/exploding gradients)discovered in1991.Sec.5.10explains how a deep RNN stack of1991(the History Compressor) pre-trained by UL helped to solve previously unlearnable DL benchmarks requiring Credit Assignment Paths(CAPs,Sec.3)of depth1000and more.Sec.5.11discusses a particular WTA method called Max-Pooling(MP)important in today’s DL FNNs.Sec.5.12mentions afirst important contest won by SL NNs in1994.Sec.5.13describes a purely supervised DL RNN(Long Short-Term Memory,LSTM)for problems of depth1000and more.Sec.5.14mentions an early contest of2003won by an ensemble of shallow NNs, as well as good pattern recognition results with CNNs and LSTM RNNs(2003).Sec.5.15is mostly about Deep Belief Networks(DBNs,2006)and related stacks of Autoencoders(AEs,Sec.5.7)pre-trained by UL to facilitate BP-based SL.Sec.5.16mentions thefirst BP-trained MPCNNs(2007)and GPU-CNNs(2006). Sec.5.17-5.22focus on official competitions with secret test sets won by(mostly purely supervised)DL NNs since2009,in sequence recognition,image classification,image segmentation,and object detection. Many RNN results depended on LSTM(Sec.5.13);many FNN results depended on GPU-based FNN code developed since2004(Sec.5.16,5.17,5.18,5.19),in particular,GPU-MPCNNs(Sec.5.19).5.11940s and EarlierNN research started in the1940s(e.g.,McCulloch and Pitts,1943;Hebb,1949);compare also later work on learning NNs(Rosenblatt,1958,1962;Widrow and Hoff,1962;Grossberg,1969;Kohonen,1972; von der Malsburg,1973;Narendra and Thathatchar,1974;Willshaw and von der Malsburg,1976;Palm, 1980;Hopfield,1982).In a sense NNs have been around even longer,since early supervised NNs were essentially variants of linear regression methods going back at least to the early1800s(e.g.,Legendre, 1805;Gauss,1809,1821).Early NNs had a maximal CAP depth of1(Sec.3).5.2Around1960:More Neurobiological Inspiration for DLSimple cells and complex cells were found in the cat’s visual cortex(e.g.,Hubel and Wiesel,1962;Wiesel and Hubel,1959).These cellsfire in response to certain properties of visual sensory inputs,such as theorientation of plex cells exhibit more spatial invariance than simple cells.This inspired later deep NN architectures(Sec.5.4)used in certain modern award-winning Deep Learners(Sec.5.19-5.22).5.31965:Deep Networks Based on the Group Method of Data Handling(GMDH) Networks trained by the Group Method of Data Handling(GMDH)(Ivakhnenko and Lapa,1965; Ivakhnenko et al.,1967;Ivakhnenko,1968,1971)were perhaps thefirst DL systems of the Feedforward Multilayer Perceptron type.The units of GMDH nets may have polynomial activation functions imple-menting Kolmogorov-Gabor polynomials(more general than traditional NN activation functions).Given a training set,layers are incrementally grown and trained by regression analysis,then pruned with the help of a separate validation set(using today’s terminology),where Decision Regularisation is used to weed out superfluous units.The numbers of layers and units per layer can be learned in problem-dependent fashion. This is a good example of hierarchical representation learning(Sec.4.4).There have been numerous ap-plications of GMDH-style networks,e.g.(Ikeda et al.,1976;Farlow,1984;Madala and Ivakhnenko,1994; Ivakhnenko,1995;Kondo,1998;Kord´ık et al.,2003;Witczak et al.,2006;Kondo and Ueno,2008).5.41979:Convolution+Weight Replication+Winner-Take-All(WTA)Apart from deep GMDH networks(Sec.5.3),the Neocognitron(Fukushima,1979,1980,2013a)was per-haps thefirst artificial NN that deserved the attribute deep,and thefirst to incorporate the neurophysiolog-ical insights of Sec.5.2.It introduced convolutional NNs(today often called CNNs or convnets),where the(typically rectangular)receptivefield of a convolutional unit with given weight vector is shifted step by step across a2-dimensional array of input values,such as the pixels of an image.The resulting2D array of subsequent activation events of this unit can then provide inputs to higher-level units,and so on.Due to massive weight replication(Sec.2),relatively few parameters may be necessary to describe the behavior of such a convolutional layer.Competition layers have WTA subsets whose maximally active units are the only ones to adopt non-zero activation values.They essentially“down-sample”the competition layer’s input.This helps to create units whose responses are insensitive to small image shifts(compare Sec.5.2).The Neocognitron is very similar to the architecture of modern,contest-winning,purely super-vised,feedforward,gradient-based Deep Learners with alternating convolutional and competition lay-ers(e.g.,Sec.5.19-5.22).Fukushima,however,did not set the weights by supervised backpropagation (Sec.5.5,5.8),but by local un supervised learning rules(e.g.,Fukushima,2013b),or by pre-wiring.In that sense he did not care for the DL problem(Sec.5.9),although his architecture was comparatively deep indeed.He also used Spatial Averaging(Fukushima,1980,2011)instead of Max-Pooling(MP,Sec.5.11), currently a particularly convenient and popular WTA mechanism.Today’s CNN-based DL machines profita lot from later CNN work(e.g.,LeCun et al.,1989;Ranzato et al.,2007)(Sec.5.8,5.16,5.19).5.51960-1981and Beyond:Development of Backpropagation(BP)for NNsThe minimisation of errors through gradient descent(Hadamard,1908)in the parameter space of com-plex,nonlinear,differentiable,multi-stage,NN-related systems has been discussed at least since the early 1960s(e.g.,Kelley,1960;Bryson,1961;Bryson and Denham,1961;Pontryagin et al.,1961;Dreyfus,1962; Wilkinson,1965;Amari,1967;Bryson and Ho,1969;Director and Rohrer,1969;Griewank,2012),ini-tially within the framework of Euler-LaGrange equations in the Calculus of Variations(e.g.,Euler,1744). Steepest descent in such systems can be performed(Bryson,1961;Kelley,1960;Bryson and Ho,1969)by iterating the ancient chain rule(Leibniz,1676;L’Hˆo pital,1696)in Dynamic Programming(DP)style(Bell-man,1957).A simplified derivation of the method uses the chain rule only(Dreyfus,1962).The methods of the1960s were already efficient in the DP sense.However,they backpropagated derivative information through standard Jacobian matrix calculations from one“layer”to the previous one, explicitly addressing neither direct links across several layers nor potential additional efficiency gains due to network sparsity(but perhaps such enhancements seemed obvious to the authors).。

Velocityfield statistics in homogeneous steady turbulence obtained using a high-resolution direct numerical simulationToshiyuki Gotoh a)Department of Systems Engineering,Nagoya Institute of Technology,Showa-ku,Nagoya466-8555,JapanDaigen FukayamaInformation and Mathematical Science Laboratory,Inc.,2-43-1,Ikebukuro,Toshima-ku,Tokyo171-0014,JapanTohru NakanoDepartment of Physics,Chuo University,Kasuga,Bunkyo-ku,Tokyo112-8551,JapanVelocityfield statistics in the inertial to dissipation range of three-dimensional homogeneous steadyturbulentflow are studied using a high-resolution DNS with up to Nϭ10243grid points.The rangeof the Taylor microscale Reynolds number is between38and460.Isotropy at the small scales ofmotion is well satisfied from half the integral scale͑L͒down to the Kolmogorov scale͑␩͒.TheKolmogorov constant is1.64Ϯ0.04,which is close to experimentally determined values.The thirdorder moment of the longitudinal velocity difference scales as the separation distance r,and itscoefficient is close to4/5.A clear inertial range is observed for moments of the velocity differenceup to the tenth order,between2␭Ϸ100␩and L/2Ϸ300␩,where␭is the Taylor microscale.Thescaling exponents are measured directly from the structure functions;the transverse scalingexponents are smaller than the longitudinal exponents when the order is greater than four.Thecrossover length of the longitudinal velocity structure function increases with the order andapproaches2␭,while that of the transverse function remains approximately constant at␭.Thecrossover length and importance of the Taylor microscale are discussed.©2002AmericanInstitute of Physics.͓DOI:10.1063/1.1448296͔I.INTRODUCTIONKolmogorov studied the statistical laws of a velocity field for small scales of turbulent motion at high Reynolds numbers.1,2Two hypotheses were introduced in his theory ͑hereafter K41for short͒:local isotropy and homogeneity exists;and there is an inertial range in the energy spectrum of theflow that is independent of viscosity and large-scale properties at sufficiently high Reynolds numbers.The most prominent conclusion of his theory is the presence of the Kolmogorov spectrum E(k)ϭK⑀¯2/3kϪ5/3in the inertial range,where⑀¯is the average rate of energy dissipation per unit mass and K is a universal constant.Since K41,there has been a considerable amount of ef-fort made to study the turbulent velocityfield statistics in the inertial range,and the energy spectrum has been a central quantity of interest.The Kolmogorov spectrum and constant have been measured infield and laboratory experiments.3–7 The exponent for the inertial range spectrum is now widely accepted asϪ5/3,with a small correction to account forflow intermittency.The Kolmogorov constant K is between1.5 and 2.After studying the results of many experiments, Sreenivasan stated that K is1.62Ϯ0.17.7The spectral theory of turbulence has also been used to predict the Kolmogorov constant.The value of K is1.77when the Lagrangian history direct interaction approximation͑LHDIA͒is used,8,9and 1.72when the Lagrangian renormalized approximation ͑LRA͒is used.10,11These are fully systematic theories that do not contain any ad hoc parameters.Direct numerical simulations͑DNSs͒of turbulentflows are now performed at higher Reynolds numbers,due to the recent dramatic increase in computational power.In the early 90’s,the resolution of DNS reached Nϭ5123grid points with a Taylor microscale Reynolds number R␭of 210ϳ240.12–20Most high-resolution DNSs have been per-formed for steady turbulence conditions to achieve high Rey-nolds numbers and obtain reliable statistics.Although results were reported with R␭greater than200,an inertial range spectrum was observed only for the lowest narrow wave number band at which forcing was applied.The Kolmogorov constant was inferred to be about1.5ϳ2in Ref.14and1.62 in Ref.16,but these results are not convincing,due to the insufficient width of the scaling range,anisotropy of theflow field,limited ensemble size,forcing techniques used,and numerical limitations of the simulations.Intermittency has also attracted the interest of research-ers.Since Kolmogorov’s intermittency theory͑hereafter K62͒,21many theoretical and statistical models of intermit-tency have been developed.4,22,23The scaling exponents of higher order structure functions for velocity differences in the inertial range were studied intensively.Intermittency in-creases with a decrease in the size of the scales of motion. The small-scale statistics gradually deviate from a Gaussian distribution,and the scaling exponents differ from those pre-dicted by K41.Experiments at very high Reynolds numbers have beena͒Electronic mail:gotoh@system.nitech.ac.jpPHYSICS OF FLUIDS VOLUME14,NUMBER3MARCH200210651070-6631/2002/14(3)/1065/17/$19.00©2002American Institute of Physicsperformed in the atmospheric boundary layer and in huge wind tunnels,and the measured scaling exponents were found to deviate from K41scaling.5,6,24,25,84However,there have been arguments made about the lack of small-scaleflow isotropy and homogeneity in these experiments,which might be affected by the large-scale shear.25,26For experiments at moderate Reynolds numbers under relatively well-controlled laboratory conditions,the width of the scaling range is usually not large enough to determine the scaling exponents precisely.Extended self-similarity͑ESS͒has been exploited to overcome this difficulty and applied to various turbulentflows in both experiments and DNSs.28–31 The idea is to measure the scaling exponents of the structure functions when they are plotted against the third order lon-gitudinal structure function,rather than to use the separation distance.The width of the scaling range is longer than that obtained with the usual method at low to moderate Reynolds numbers.The scaling exponents are anomalous,but do agree with those obtained from high Reynolds number experiments up to a certain order.24,28–31However,there is no consensus as to why the structure functions give a longer inertial range, or what is missing from theflow statistics as a result.Also, there is no unique way to determine the scaling exponents for the transverse and mixed velocity structure functions,be-cause those higher order structure functions can be plotted against other types of third order structure functions as well as the third order longitudinal structure function.There also have been arguments about whether the scal-ing exponents for the longitudinal and transverse structure functions at small scales are equal.25–27,32–38Many experi-ments and DNSs have reported that higher order longitudinal scaling exponents are larger than transverse ones.However, some researchers have argued that the difference is due to deviation from the assumed conditions,such as local homo-geneity,isotropy,and the independence of small scales from macroscale parameters.They have suggested that when the Reynolds number becomes large enough,the difference will vanish.36,37,39,40In many aspects of turbulence research,there have been questions posed about the extent to which the local homoge-neity and isotropy of the turbulent velocityfield are attained. This will affect the small-scale statistics significantly.Recent experimental studies have shown that local isotropy is par-tially satisfied for lower order moments.25,26,37However,it is not sufficient to examine only the conditions assumed in the above studies,and only a limited knowledge of the trueflow conditions is available so far.26,37,38A DNS with a sufficiently large grid size provides abetter opportunity to examine the points raised above.It hasthe advantage that any physical quantity can be measureddirectly without deforming theflowfield.In the presentstudy,a series of large scale DNSs have been performed at ahigh resolution of up to Nϭ10243and R␭ϭ460.41–45The inertial range of the turbulencefield has a considerablelength,and useful velocity statistics can be extracted such asthe Kolmogorov constant,the energy spectrum,velocitystructure functions up to the tenth order,their scaling expo-nents,and probability density functions for velocity differ-ences.To the authors’best knowledge,these are thefirstDNS data in the inertial range;the data provide new insightinto the inertial and dissipation ranges.The main purposes of the present paper are to describethe statistics of the velocityfield in an incompressible steadyturbulentflow obtained from the DNS,and to reexaminecurrent knowledge of turbulence,developed since K41.Thepaper is organized as follows.The numerical aspects of thepresent DNS are described in Sec.II,and the energy spec-trum is examined in Sec.III.The variation of single pointquantities and probability density functions͑PDFs͒with theReynolds number is discussed in Sec.IV.The isotropy of thesecond and third order moments of the velocity difference isexamined in Sec.V,and the energy budget is examined interms of the Ka´rma´n–Howarth–Kolmogorov equation inSec.VI.The structure functions and scaling exponents arediscussed in Sec.VII.Section VIII presents an analysis ofthe crossover lengths of the structure functions.Finally,asummary and conclusions are provided in Sec.IX.II.NUMERICAL SIMULATIONThe Navier–Stokes equations are integrated in Fourierspace for unit density:ͩץץtϩ␯k2ͪuϭP͑k͒•F͓uÃ␻͔kϩf,͑1͒͗f͑k,t͒f͑Ϫk,s͒͘ϭP͑k͒F͑k͒4␲k2␦͑tϪs͒,͑2͒where␻is the vorticity vector,P͑k͒is the projection opera-tor,F denotes a Fourier transform,and f is a solenoidal Gaussian random force that is white in time.The spectrum of the random force F(k)is constant over the low wave number band and zero otherwise;the force is normalized asTABLE I.DNS parameters and statistical quantities of the runs.T eddya v is the period used for the time average.R␭N k max␯c f Forcing range T eddya v E⑀¯L␭␩(ϫ10Ϫ2)K 38128360 1.50ϫ10Ϫ2 1.30ͱ3рkрͱ1222.6 1.99 1.190.8910.501 4.10¯5425631217.00ϫ10Ϫ30.70ͱ3рkрͱ1214.9 1.390.6270.8290.393 2.72¯702563121 4.00ϫ10Ϫ30.50ͱ3рkрͱ1249.7 1.160.4570.7850.318 1.93¯1255123241 1.35ϫ10Ϫ30.50ͱ3рkрͱ12 5.52 1.250.4920.7440.1850.841¯2845123241 6.00ϫ10Ϫ40.501рkрͱ6 3.03 1.960.530 1.2460.1490.449 1.64 38110243483 2.80ϫ10Ϫ40.511рkрͱ6 4.21 1.740.499 1.1390.09890.258 1.63 46010243483 2.00ϫ10Ϫ40.511рkрͱ6 2.14 1.790.506 1.1500.08410.199 1.64 1066Phys.Fluids,Vol.14,No.3,March2002Gotoh,Fukayama,and Nakano͵ϱF ͑k ͒dk ϭ⑀¯in ,͑3͒where ⑀¯in is the average rate of the energy input per unit mass.A pseudo-spectral code was used to compute the con-volution sums,and the aliasing error was effectively re-moved.The time integration was performed using the fourth order Runge–Kutta–Gill method.Physical quantities of turbulent flow include the total energyE ͑t ͒ϭ12͗u 2͘ϭ32u ¯2ϭ͵ϱE ͑k ͒dk ,͑4͒the average energy dissipation per unit mass⑀¯ϭ2␯͵ϱk 2E ͑k ͒dk ,͑5͒the integral scaleL ϭͩ3␲4͵ϱk Ϫ1E ͑k ͒dkͪͲE ,͑6͒the Taylor microscale␭ϭͩ5EͲ͵ϱk 2E ͑k ͒dkͪ1/2,͑7͒the Taylor microscale Reynolds numberR ␭ϭu ¯␭␯,͑8͒and the Kolmogorov scale␩ϭͩ␯3⑀¯ͪ1/4.͑9͒The range of the Taylor microscale Reynolds number was 38to 460.The characteristic parameters of the DNS are listed in Table I.43Most of these are identical to Gotoh and Fukayama,43but the averaging time for R ␭ϭ381was ex-tended to 4.21large eddy turnover times.A statistically steady state was confirmed by observing the time evolutionof the total energy,the total enstrophy,and the skewness of the longitudinal velocity derivative.The statistical averages were computed as time averages over tens of large eddy turnover times for the lower Reynolds number flows,and over a few large eddy turnover times for the higher Reynolds number flows.The resolution condition k max ␩Ͼ1was satis-fied for most runs,except for R ␭ϭ460in which k max ␩was slightly less than unity (k max ␩ϭ0.96).This does not ad-versely affect the results in the inertial range.The computational time required for runs at a N ϭ10243resolution varied,depending on the statistical data that was gathered.Typically,60h was required for one large eddy turnover time.The total time of the computations was more than 500h for the longest run (R ␭ϭ381).Data col-lected during the transition period to steady state ͑about six large eddy turnover times ͒were discarded.The relatively long time required to attain steady state was due to the low wave number band forcing.This imposes a severe computa-tional putations with R ␭р284were per-formed on a Fujitsu VPP700E parallel vector machine with 16processors at RIKEN.Simulations of higher R ␭were per-formed on a Fujitsu VPP5000/56with 32processors at the Nagoya University Computation Center.III.ENERGY SPECTRUMFigure 1shows the three-dimensional energy spectrum calculated for each run.All of the curves are scaled to the Kolmogorov units and multiplied by k 5/3.As the Reynolds number increases,the curves extend toward lower wave numbers.The curves of flows with Reynolds numbers larger than R ␭ϭ284contain a finite plateau,which indicates that E (k )ϰk Ϫ5/3.There is a bump when 0.04рk ␩р0.3at the high end of the inertial range,which is consistent with pre-vious experimental and numerical observations.6,16The nor-malized energy transfer flux,defined by1⑀¯⌸͑k ͒ϭ1⑀¯͵kϱT ͑k Ј͒dk Ј͑10͒is shown in Fig.2,where T (k )is a nonlinear energy transfer function in the energy spectrum equation.4,22Between0.007рk ␩р0.04,⌸(k )/⑀¯is approximately constantand FIG.1.Scaled energy spectra,⑀¯Ϫ1/4␯Ϫ5/4(k ␩)5/3E (k ).The inertial range is 0.007рk ␩р0.04and K ϭ1.64Ϯ0.04.A horizontal line indicates K ϭ1.64.FIG.2.Normalized energy transfer flux,⌸(k )/⑀¯for R ␭ϭ381and 460.1067Phys.Fluids,Vol.14,No.3,March 2002Velocity field statistics in homogeneousclose to unity;thus the flow is in an equilibrium state over the inertial range of the energy spectrum,corresponding to the plateaus in Fig.1.The Kolmogorov constant given in Table I is determined using a least square fit between 0.007рk ␩р0.04on the R ␭Ͼ284curves.In Ref.43,the Kolmog-orov constant was reported as K ϭ1.65Ϯ0.05.However,the averaging time has since been extended for the R ␭ϭ381run.The R ␭ϭ478run differs slightly from statisticalequilibrium,since ⌸(k )/⑀¯is not exactly one;for this reason,the R ␭ϭ478data were not used for this analysis.The Kol-mogorov constant,computed using the data only from the R ␭ϭ381and 460runs,isK ϭ1.64Ϯ0.04,͑11͒which is in good agreement with experimental values and recent DNS data.7,16There are many DNSs reporting the Kolmogorov constant higher than the value 1.64.However,the length of the inertial range in those DNSs is not long enough to clearly observe the k Ϫ5/3range,and the top of the bump of the compensated energy spectrum k 5/3E (k )is un-derstood as the inertial range,so that the Kolmogorov con-stant is read as about 2as seen in Fig.1.16The Kolmogorov constant 1.64is also close to the value obtained using the LHDIA ͑1.77͒,8,9the LRA ͑1.72͒.10,11These spectral theories of turbulence are consistent with Lagrangian dynamics,arederived systematically,and contain no ad hoc parameters.Figure 3shows the one-dimensional energy spectrum ob-tained from the present DNS with R ␭ϭ460,from experi-ments,and from the LRA.The agreement between the curves is satisfactory.Therefore we conclude that the present DNS has successfully calculated a homogeneous turbulent flow field in the inertial range of the energy spectrum.IV.ONE-POINT STATISTICS A.MomentsSome one-point moments of the velocity field areS 3͑u ͒ϵ͗u 3͗͘u 2͘3/2,S 3͑u x ͒ϵ͗u x 3͗͘x 2͘3/2,͑12͒K 4͑u ͒ϵ͗u 4͗͘u 2͘2,K 4͑u x ͒ϵ͗u x 4͗͘u x 2͘2,K 4͑u y ͒ϵ͗u y 4͗͘u y 2͘2,͑13͒where u is the velocity component in the x direction.The variation of these moments with the Reynolds number is shown in Fig.4and listed in Table II.The general behavior of the curves is consistent with previous DNS and experi-mental data.13,14,18,19,26,46,47There are small effects of rela-tively low resolution on S 3and K 4for the velocity deriva-tives for R ␭ϭ381and 460data.The skewness factor of the velocity u is very small for runs with the R ␭р125,and isofparison of one-dimensional energy spectra.Symbols:experi-ments,solid line:present DNS (R ␭ϭ460),dashed line:statistical theory ͑LRA and MLRA ͒.FIG.4.Variation of the moments of the velocity and velocity gradient withthe Reynolds number.Line:present DNS,circle:K 4(u y )͑Jime´nez et al.,Ref.13͒,solid square:K 4(u )͑Jime ´nez et al.,Ref.13͒,square:K 4(u x )͑Wang et al.,Ref.14͒,plus:K 4(u x )͑Vedula and Yeung,Ref.18͒,star:ϪS 3(u x )͑Wang et al.,Ref.14͒.TABLE II.Moments of the velocity and velocity derivatives.R ␭S 3(u )K 4(u )S 3(ץu /ץx )K 4(ץu /ץx )K 4(ץu /ץy )380.0227 2.89Ϫ0.520 4.14 5.16540.00563 2.86Ϫ0.517 4.47 6.00700.00473 2.93Ϫ0.519 4.81 6.621250.0820 2.94Ϫ0.529 5.658.192840.0231 2.77Ϫ0.531 6.6310.1381Ϫ0.246 2.98Ϫ0.5747.9012.2460Ϫ0.1682.89Ϫ0.5457.9111.71068Phys.Fluids,Vol.14,No.3,March 2002Gotoh,Fukayama,and Nakanothe order of 0.2for runs with the R ␭у284.The relatively large values of the velocity skewness are caused by the shorter averaging time used compared to the low Reynolds number runs.Since most of the energy resides in the lowest wave number band,there are persistent large fluctuations of the large scales of motion over longer time period.The longer time average or the forcing at larger wave numbers would yield smaller velocity skewness.The flatness factor of the velocity field is close to three,which is the Gaussian value.The skewness factor of the longitudinal velocity deriva-tives is very insensitive to the Reynolds number,S 3͑u x ͒ϰR ␭0.0370,͑14͒where the exponent is determined by a least square fit.Theaverage value is Ϫ0.53,which is consistent with experimen-tal observations over the range of Reynolds numbers studied in the present work.However,the exponent is smaller than indicated by the experimental data.26,46The flatness factors for the longitudinal and transverse velocity derivatives in-crease with the Reynolds number asK 4͑u x ͒ϰR ␭0.266,K 4͑u y ͒ϰR ␭0.335.͑15͒The exponent of K 4(u y )is larger than that of K 4(u x );thus,the PDF for the transverse velocity derivative has longer tails than those of the longitudinal velocity derivative.From ex-perimental observations,Shen and Warhaft reported thatK 4(u x )ϰR ␭0.37and K 4(u y )ϰR ␭0.25.26Since there is scatter in the experimetal data,the exponents in Eq.͑15͒by the present DNS are not inconsistent with the experimental data.Van Atta and Antonia studied the Reynolds number dependence of S 3(u x )and K 4(u x ),46and found thatS 3͑u x ͒ϰR ␭0.12,K 4͑u x ͒ϰR ␭0.32for ␮ϭ0.2,͑16͒S 3͑u x ͒ϰR ␭0.15,K 4͑u x ͒ϰR ␭0.41for ␮ϭ0.25,͑17͒where ␮is the exponent defined by ͗⑀r 2͘ϰr Ϫ␮for the locally averaged energy dissipation rate.4,21Generally,the Reynolds number dependency of S 3and K 4in our DNSs is weaker than observed in the experiments,irrespective of the type of forcing used.We believe this is because the range of Rey-nolds numbers in DNS is smaller than experimental flows,and there remain small-scale anisotropy effects in the experi-ments.B.Probability density functionsThe probability density function conveys information about single-point velocity statistics.It has been one of the central issues of turbulence research in the last decade.Single-point PDFs for the velocity and its derivatives are shown in Figs.5–7.A longer time period was necessary for the time average to obtain well-converged PDF for the ve-locity Q (u ).The distribution Q (u )is close to Gaussian,and its tail extends to very low values of the order of 10Ϫ10.Such values have not been reported in the literature.The Q (u )curve for R ␭ϭ381is skewed negatively,but this is attributed to the insufficient time-averaging period ͑four large eddy turnover times ͒that was used.The overall trend is that Q (u )decays faster than a Gaussian distribution at large ampli-tudes.This behavior was also observed in one-dimensional decaying and forced Burgers turbulence.48,49Jime´nez has shown that the PDF Q (u )is slightly sub-Gaussian as the energy spectrum decays faster than k Ϫ1.50FIG.5.Variation of velocity PDF with the Reynoldsnumber.FIG.6.Variation of the longitudinal velocity derivative PDF with the Rey-noldsnumber.FIG.7.Variation of the transverse velocity derivative PDF with the Rey-nolds number.1069Phys.Fluids,Vol.14,No.3,March 2002Velocity field statistics in homogeneousThis is consistent with the present DNS results.Studies of the Q (u )tail predict that Q (w )ϰexp(Ϫc ͉w ͉3)when the forc-ing has a short correlation time.51,52Here,w ϭu /͗u 2͘1/2is the normalized velocity amplitude and c is a nondimensional constant.The asymptotic form of Q (u )was examined by plotting ln ͓Ϫln(Q (w )͔against ln ͉w ͉;however,the Q (w )tails were too short to determine the true asymptotic form.The PDF for the longitudinal velocity derivative is slightly skewed,as expected from the finite negative value of the skewness factor.The tail becomes longer as the Reynolds number increases.Figure 7shows that the PDF of the trans-verse derivative is symmetric and has a longer tail than the longitudinal derivative.There are many theories for the PDF of the velocity derivative.The asymptotic tail of Q (ץu /ץy )is presented in Fig.8,in which both the positive and negative sides are plotted by assuming that the PDF is symmetric.The tails gradually become longer as the Reynolds number increases;therefore,Q (s )is Reynolds-number dependent,and cannot be represented in a single stretched exponential form as Q (s )ϰexp(Ϫb ͉s ͉h ),where s is the normalized amplitude of ץu /ץy and b is a nondimensional constant that is a function of the Reynolds number.53V.ISOTROPYThe hypothesis of isotropy of the flow field is one of the key components of K41.There are various methods to ex-amine the degree of isotropy.One measure of isotropy can be obtained from the relations between the second and third order longitudinal and transverse velocity structure func-tions.These areD LL ϵ͗͑␦u r ͒2͘,D TT ϵ͗͑␦v r ͒2͘,͑18͒D LLL ϵ͗͑␦u r ͒3͘,D LTT ϵ͗␦u r ͑␦v r ͒2͘,͑19͒where␦u r ϵ͑u ͑x ϩr ͒Ϫu ͑x ͒͒•r /r ,͑20͒␦v r ϵ͑u ͑x ϩr ͒Ϫu ͑x ͒͒•͑I Ϫrr /r 2͒•e Ќ,͑21͒and e Ќis the unit vector perpendiculer to r ,and I is the unit tensor.Then the isotropy and incompressibility relations areD TT ͑r ͒ϭD LL ͑r ͒ϩr 2dD LL ͑r ͒dr ,͑22͒D LTT ͑r ͒ϭ16ddrrD LLL ͑r ͒.͑23͒In DNS,the solenoidal property of the Fourier amplitude velocity vector u ͑k ͒is always satisfied to the level of nu-merical error,which is smaller than 10Ϫ15.Thus,the accu-racy of the above relations depends solely on the deviation from isotropy.The two sides of Eqs.͑22͒and ͑23͒are com-pared for R ␭ϭ125,381,and 460in Figs.9and 10.The curves in the figures are divided by r 2/3and r ,respectively,and the vertical axes of the plots are linear.The thick lines represent the left hand sides of Eqs.͑22͒and ͑23͒,and the thin lines correspond to the right-hand sides.The isotropy of the second and third order moments is excellent for scales less than L /2.The difference at larger separations is caused by the anisotropy due to the small number of energy-containing Fourier modes.The curves for R ␭ϭ381and460FIG.8.Variation of the asymptotic tail of the transverse velocity derivative PDF with the Reynolds number.Both positive and negative sides are plot-ted.The rightmost curve corresponds to R ␭ϭ460.FIG.9.Isotropy relation at the second order.Thin line:D TT (r )r Ϫ2/3,thick line:(D LL (r )ϩ(r /2)(dD LL (r )/dr ))r Ϫ2/3.L /␩and ␭/␩are shown for R ␭ϭ460.FIG.10.Isotropy relation at the third order.Thin line:D LTT (r )r Ϫ1,thick line:((1/6)(d /dr )rD LLL (r ))r Ϫ1.L /␩and ␭/␩are shown for R ␭ϭ460.1070Phys.Fluids,Vol.14,No.3,March 2002Gotoh,Fukayama,and Nakanoin Fig.9are not horizontal,suggesting that the second order structure function does not scale as r 2/3.The scaling expo-nents will be examined later in this paper.The isotropic re-lations,such as D 1122ϭD 1133and D 2222ϭ3D 2233ϭD 3333,and Hill’s higher order relations were not computed.54VI.KA´RMA ´N–HOWARTH–KOLMOGOROV EQUATION The energy budget for various scales is described by the Ka´rma ´n–Howarth–Kolmogorov ͑KHK ͒equation,45⑀¯r ϭϪD LLL ϩ6␯ץD LL ץrϩZ ͑24͒for steady turbulence,4,55,56where Z (r )denotes contributions due to the external force given by Z ͑r ,t ͒ϭ͵Ϫϱt͗␦f ͑r ,t ͒•␦f ͑r ,s ͒͘dsϭ12r͵0ϱͩ115ϩsin kr ͑kr ͒3ϩ3cos kr ͑kr ͒4Ϫ3sin kr ͑kr ͒5ͪF ͑k ͒dk .͑25͒Since the external force spectrum F (k )is localized in arange of low wave numbers,the asymptotic form of Z (r )for small separations is given asZ ͑r ͒ϭ235⑀¯in k f 2r 3,k f 2ϵ͐0ϱk 2F ͑k ͒dk͐0ϱF ͑k ͒dk.͑26͒A generalized Ka´rma ´n–Howarth–Kolmogorov equation has also been derived:57–6343⑀¯r ϭϪ͑D LLL ϩ2D LTT ͒ϩ2␯ץץr͑D LL ϩ2D TT ͒ϩW ,͑27͒where W ͑r ͒ϭ4r͵0ϱͩ13ϩcos kr ͑kr ͒2Ϫsin kr͑kr ͒3ͪF ͑k ͒dk ,Ϸ215⑀¯in r 3k f 2for ͉k f r ͉Ӷ1.͑28͒Equation ͑24͒is recovered by substituting Eqs.͑22͒and ͑23͒into Eq.͑27͒.Figure 11shows the results obtained when each term of Eq.͑24͒is divided by ⑀¯r for R ␭ϭ460.Curves in which r /␩is larger than r /␩ϭ1200are not shown,because the sign of D LLL changes.A thin horizontal line indicates the Kolmog-orov value 4/5.When the separation distance decreases,the effect of the large scale forcing used in the present DNS decreases quickly,while the viscous term grows gradually.The third order longitudinal structure function D LLL quickly rises to the Kolmogorov value,remains there over the iner-tial range ͑between r /␩Ϸ50and 300͒,and then decreases.In the inertial range,the force term decreases as r 3according to Eq.͑26͒,while the viscous term increases as r ␨2Ϫ1(␨2Ͻ1)when r decreases.͓Since each term in the figure is divided by (⑀¯r ),the slope of each curve is 2and ␨2Ϫ2,respectively.͔The sum of the three terms in the right hand side of Eq.͑24͒divided by ⑀¯r is close to 4/5,the Kolmogorov value.The deviation of the sum from the 4/5law at the smallest scales is due to the slightly lower resolution of the data at these scales ͑k max ␩is close to one ͒.At larger scales greater than r /␩ϭ700,the deviation is caused by the finiteness oftheFIG.11.Terms in the Ka´rma ´n–Howarth–Kolmogorov equation when R ␭ϭ460.Thin solid line:4/5.FIG.12.Kolmogorov’s 4/5law.L /␩and ␭/␩are shown for R ␭ϭ460.The maximum values of the curves are 0.665,0.771,0.781,and 0.757for R ␭ϭ125,284,381,and 460,respectively.FIG.13.Terms in the generalized Ka´rma ´n–Howarth–Kolmogorov equation for R ␭ϭ460.Thin solid line:4/3.1071Phys.Fluids,Vol.14,No.3,March 2002Velocity field statistics in homogeneousensemble,which indicates the persistent anisotropy of the larger scales.The above findings are consistent with the cur-rent knowledge of turbulence developed since Kolmogorov,although confirmation of some aspects of turbulence using actual data is new from both a numerical and experimental point of view.56,59–65It is interesting and important to observe when the Kol-mogorov 4/5law is satisfied as the Reynolds numberincreases.6,66–69Figure 12shows curves of ϪD LLL (r )/(⑀¯r )for various Reynolds numbers.In this figure,the 4/5law applies when the curves are horizontal.The portion of the curves in which r /␩Ͼ1200is not shown.Although there is a small but finite horizontal range when R ␭Ͼ284,the level of the plateau is still less than the Kolmogorov value.The maximum values of the curves are 0.665,0.771,0.781,and 0.757for R ␭ϭ125,284,381,and 460,respectively.The value 0.781for R ␭ϭ381is 2.5%less than 0.8.An asymptotic state is approached slowly,which is consistent with recent studies.However,the asymptote is approached faster than predicted by the theoretical estimate.66,69Theslow approach is due to the fact that D LLL (r )is the third order structure function and most positive contributions are canceled by negative ones.Thus only the slight asymmetry of the ␦u r PDF contributes to D LLL .The level of the plateau of the R ␭ϭ460curve is slightly less than the others.A higher value would be expected if the time average period used for the R ␭ϭ460run were longer.The generalized Ka´rma ´n–Howarth–Kolmogorov equa-tion Eq.͑27͒is also examined in a similar fashion.Figure 13shows each term of the equation divided by ⑀¯r ;a horizontal line indicates the 4/3law.The agreement between the present data and theory is satisfactory.The third order moment slowly approaches the Kolmogorov value 4/3,as shown in Fig.14.The maximum values of the curves of the 4/3law are 0.564,1.313,1.297,and 1.259for R ␭ϭ125,284,381,and 460,respectively.VII.STRUCTURE FUNCTIONS AND SCALING EXPONENTSThe velocity structure functions are defined asS p L ͑r ͒ϭ͉͗␦u r ͉p͘,S p T ͑r ͒ϭ͉͗␦v r ͉p͘,FIG.14.Kolmogorov’s 4/3law.L /␩and ␭/␩are shown for R ␭ϭ460.The maximum values of the curves for the 4/3law are 0.564,1.313,1.297,and 1.259for R ␭ϭ125,284,381,and 460,respectively.FIG.15.Variation of the ␦u r PDF with r for R ␭ϭ381.From the outermostcurve,r n /␩ϭ2n Ϫ1dx /␩ϭ2.38ϫ2n Ϫ1,n ϭ1,...,10,where dx ϭ2␲/1024.The inertial range corresponds to n ϭ6,7,8.Dotted line:Gaussian.FIG.16.Variation of PDF for ␦v r with r at R ␭ϭ381.The classification of curves is the same as in Fig.17.FIG.17.Convergence of the tenth order accumulated moments C 10(␦u r )at R ␭ϭ381for various separations r n /␩ϭ2.38ϫ2n Ϫ1,n ϭ1,...,10.Curves are for n ϭ1,...10from the uppermost,and the inertial range corresponds to n ϭ6,7,8.1072Phys.Fluids,Vol.14,No.3,March 2002Gotoh,Fukayama,and Nakano。