Jtam.dvi JOURNAL OF THEORETICAL AND APPLIED MECHANICS 45, 3, pp. 705-724, Warsaw 2007 MODELLING OF NEAR-WALL TURBULENCE WITH LARGE-EDDY VELOCITY MODES Marta Wacławczyk Jacek Pozorski Institute of Fluid Flow Machinery, Polish Academy of Sciences, Gdańsk, Poland e-mail: mw@imp.gda.pl; jp@imp.gda.pl In the paper, low-order modelling of the turbulent velocity field in the near-wall region is performed using the Proper Orthogonal Decomposi- tion (POD) approach. First, an empirical eigenfunction basis is compu- ted, basing on two-point velocity correlations.Next, theGalerkinprojec- tion of the Navier-Stokes equations on the truncated basis is performed. This results in a system of Ordinary Differential Equations (ODEs) for time-dependent coefficients. Evolution of the largest vortical structures in the near-wall zone is then obtained from the time dependent coeffi- cients and eigenfunctions. The systemapplied in the presentwork consi- sts of 20ODEs, the reconstructedvelocity field is two-dimensional in the plane perpendicular to the main flow direction. Moreover, the filtering procedure associatedwith the PODmethod is discussed, the PODfilter is derived and compared with LES filters. Key words: proper orthogonal decomposition, coherent structures, turbulent channel flow 1. Introduction In the near-wall buffer and logarithmic regions, generation of turbulent energy and transport of heat andmomentum are mainly controlled by the dynamics of large turbulent eddies. Analysis of experimental data performed by Bake- well and Lumley (1967) proved that these so-called coherent vortices have the characteristic shape of rolls elongated in the streamwise direction. The insight into the instantaneous phenomena in the near-wall zone is useful in many ap- plications, including the study of heat transfer (Kasagi et al., 1988), diffusion and dispersion problems (Joia et al., 1998; Picciotto et al., 2005; Łuniewski et al., 2006) as well as study of flow control (Gunes and Rist, 2004). However, the numerical cost associated with the Direct Numerical Simulations (DNS) 706 M. Wacławczyk, J. Pozorski as well as Large Eddy Simulations (LES) of the near-wall region is very high. Hence, there is a need for other simplified methods which would predict, at least qualitatively, the time evolution of the largest eddies. Adevelopment along these lines is thePODapproach (Aubry et al., 1988). The method is based on the Lumley (1967) proposal to expand velocity in a series of eigenfunctions computed from the eigenvalue problem with experi- mentally determined two-point correlations. The problem is defined in such a way that the function basis is optimal for the representation of kinetic tur- bulent energy for a given flow. Initially, the POD method was used to educt coherent structures from the turbulent field (Bakewell and Lumley, 1967; Mo- in and Moser, 1989; Deville and Bonnet, 2002). Another application of the POD is to construct a low-order system of ODEs and resolve the dynamics of the largest coherent structures. For this purpose, the momentum equation is subjected to the Galerkin projection on the eigenfunction basis and a set of ordinary differential equations are solved for time-dependent coefficients. The instantaneous velocity field can be then reconstructed from eigenfunc- tions and their time-dependent coefficients. The PODmethod was applied to simulations of various flow configurations like, e.g., boundary layers and chan- nel flows (Holmes et al., 1996; Podvin, 2001; Aubry et al., 1988; Ball et al., 1991;Webber et al., 1997), mixing layers and jets (Citriniti andGeorge, 2000; Deville et al., 2000), ventilated cavities (Allery et al., 2006) and other flow cases. The drawback of such simulations is that they can be applied only to a specific geometry for which the two-point statistics have been predetermined experimentally. Yet, this approach can still be very useful, e.g., for analysing dynamics of passive scalars. The present paper deals with POD simulations of the velocity field in the near-wall region of a turbulent channel flow. It can be regarded as continu- ation of the previous work (Wacławczyk and Pozorski, 2002) where two-point velocity correlations were computed from the PIV experiment and a typical picture of large eddies was reconstructed; however, simulations of dynamical system were not performed. In this work, the DNS results of the turbulent channel flow at Reτ = 140 (courtesy of Dr. Bérengère Podvin) are used to compute the necessary input data for POD: the two-point correlations, eigen- functions characteristic for the near-wall flow, and coefficients of the ODEs. Contrary to the work of Aubry et al. (1988) where the simulations in the near-wall region were performed with 10 POD modes, 20 active modes are employed here. Henceforth, the present approach is called 20D model. The dynamics of large eddies is reconstructed from the time-dependent solution of ODEs and the spatially-dependent eigenfunction basis. The simulations are two-dimensional; this decreases the numerical cost but, on the other hand, the resulting picture of turbulentmotion is only qualitative. Another contribution of the present work is the derivation of the filter function associated with the Modelling of near-wall turbulence... 707 PODmethod. This is done in order to discuss the relations between the POD and LESmethods. 2. Eigenfunction expansion The POD is a method of analysis of stochastic fields, which leads to the em- pirical eigenfunction basis, optimal for the energy representation (Berkooz et al., 1993; Holmes et al., 1996); this means that the series of basis functions converge more rapidly (in the quadratic mean) than any other representa- tion. The POD method relies on numerically or experimentally determined two-point velocity correlations, so the basis is not given analytically (unlike the Fourier basis or the Chebyshev polynomials used in numerical spectral methods for LES and DNS). If we recall that the coherent structures can be defined as high-energy regions, the choice of the POD basis, optimal in the energy representation, seems to be adequate to describe the dynamics of large eddies. In order to explain in more detail how the two-point correlations are in- volved into the considerations, we first express the fluctuating velocity of a turbulent field realisation ω as a linear combination of empirical orthogonal functions ψ(n)(x) with random and statistically independent coefficients a (n) ω uω(x)= ∞∑ n=1 a(n)ω ψ (n)(x) (2.1) where ( ψ(n)(x),ψ(l)(x) ) = ∫ Ω ψ(n)(x)ψ(l)∗(x) dx = δnl (2.2) the parenthesis (·, ·) denotes the scalar product; x could either be a position vector or any spatial or temporal coordinate; Ω is the domain of integration, * denotes the complex conjugate and δnl is the Kronecker delta. To find the optimal basis, we seek for functions ψ(n) which maximise the expression 〈(uω,ψ(n))(uω,ψ(n))∗〉 (ψ(n),ψ(n)) = 〈(a(n)ω )2〉= λ(n) (2.3) (the Einstein sum notation convention does not apply here). Next, the va- riational method (Lumley, 1970) is used in order to find the maximum of λ; this leads to the eigenfunction and eigenvalue problem with the two-point correlations ( 〈u(x)u(x′)〉,ψ∗(x′) ) = λψ(x) (2.4) 708 M. Wacławczyk, J. Pozorski The solution to (2.4) gives a class of functions ψ(n)(x) which form, for a given flow case, the optimal POD basis. The basis can be further truncated since only a few first eigenfunctions contain the majority of kinetic energy of the turbulent flow field and characterise its inhomogeneity. In the particular case of homogeneous turbulent fields, Eq. (2.4) is satisfied by the Fourier basis. The considerations can be readily extended to the case of 3D fully develo- ped turbulent channel flow, statistically homogeneous in the streamwise and spanwise directions. The fluctuating velocity ui = Ui−〈Ui〉 reads ui(x,y,z,t) = = 1√ LxLz ∞∑ n=1 ∞∑ kx=−∞ ∞∑ kz=−∞ ankxkz(t)ψ (n) ikxkz (y)exp ( ix 2πkx Lx +iz 2πkz Lz ) = = 1√ LxLz ∑ p ap(t)φ (p) i (y) (2.5) where Lx and Lz denote the size of the domain (in wall units), in x and z direction respectively, ∑ p stands for the triple sum ∑ n ∑ kx ∑ kz , and the symbol φ (p) i (y) is introduced for brevity.The intervals x−x′ and z−z′ denoted by ∆x and ∆z, respectively, are related to the wavenumbers kx, kz by the formulae ∆x = Lx/kx and ∆z = Lz/kz. In the fully developed turbulent channel flow, the eigenvalue problem (2.4) is solved in the y direction only, for each pair of kx and kz Ly∫ 0 Φij(kx,y,y ′,kz)ψjkxkz(y ′) dy′ = λψikxkz(y) (2.6) where Ly denotes the domain length in the wall-normal direction and Φij is the Fourier transform of the two-point velocity correlation tensor Qij Φij(kx,y,y ′,kz)=Fxz{Qij(∆x,y,y′,∆z)}= = 〈Fxz{ui(x,y,z)}F∗xz{uj(x′,y′,z′)}〉 The Fourier transform of a given quantity A(x,z) is defined as Fxz{A(x,z)}= 1√ LxLz Lx∫ 0 Lz∫ 0 A(x,z)exp ( −ix2πkx Lx − iz2πkz Lz ) dx dz (2.7) The complex coefficients ankxkz(t) in Eq. (2.5) determine the time evolution of the fluctuation field and are themselves governed by a system of ordinary differential equations. TheseODEs are obtained by introducing the eigenfunc- tion expansion into theNavier-Stokes equation for fluctuating velocity (Aubry et al., 1988) and subjecting the resulting formula to the Galerkin projection (details are presented in Section 3). Modelling of near-wall turbulence... 709 3. Galerkin projection 3.1. Governing equations As stated in Section 1, the second application of the POD method is the reduced-order modelling. In this case, the time evolution of complex coef- ficients ankxkz(t) is sought for and the instantaneous velocity field of large eddies is reconstructed fromEq. (2.5). In a general case, the knowledge of the complete two-point correlations Qij(x,x ′,y,y′,z,z′), i.e. a tensor function of 6 variables, is required to fullymodel the large eddydynamics.However, in the case of homogeneity in the streamwise and spanwise directions, the tensor Qij becomes a function of ∆x = x− x′ and ∆z = z − z′. Aubry et al. (1988), who considered the turbulent boundary-layer flow, assumed additionally the zero streamwise variation, i.e. analysed only Qij = Qij(0,y,y ′,∆z). In spite of this simplification, the resulting velocity field still reveals the major features of near-wall turbulent flow like the intermittent, bursting behaviour. Now, thederivationof theordinarydifferential equations for thecoefficients ankxkz(t) will be outlined (cf. Holmes et al., 1996). The new contribution of this theoretical part of the paper is a discussion of relations between the POD and the LES approach. In particular, a formula for the filtering function asso- ciated with the POD method is derived. Holmes et al. (1996) chose the fully developed turbulent boundary layer as a particular case of study. First, the averaging procedure needs to be defined. As the flow is assumed to be stati- stically homogeneous in the streamwise and spanwise directions, the mean of a physical quantity A(x,y,z,t) is written as 〈A(y,t)〉= 1 LxLz Lx∫ 0 Lz∫ 0 A(x,y,z,t) dx dz (3.1) and hence it is, in general, a function of time.When the averaged momentum equation is subtracted from the instantaneous one (the Navier-Stokes eq.), we obtain ∂ui ∂t +〈U1〉 ∂ui ∂x + ∂〈U1〉 ∂y u2δi1+ [∂uiuj ∂xj −∂〈uiuj〉 ∂xj ] =−1 ̺ ∂p ∂xi +ν ∂2ui ∂xk∂xk (3.2) Here we also assumed that 〈U2〉 = 〈U3〉 = 0 and that 〈U〉 = 〈U1〉 is slowly varying in time, hence its time derivative can be neglected. In the present derivation, following Aubry et al. (1988), the mean velocity profile will be described by the analytical formula 〈U〉(y) = 1 ν y∫ 0 〈uv〉 dy+ u2 ∗ ν ( y− y2 2H ) (3.3) where H is the channel half-width. 710 M. Wacławczyk, J. Pozorski The next step is to decompose the fluctuating velocity into the coherent part u<, and the unresolved, backgroundpart, u>. The coherent term u< will be represented by sum (2.5) after truncation at some n = N, kx =±Kx and kz =±Kz, while the influence of the unresolved scales u> has to bemodelled. In order to derive the evolution equation for u<, Eq. (3.2) will be filtered by analogy to the procedure in LES. 3.2. POD filter Wecanassumethat thePODmethodprovides a certain kindof filter G(x,x′), hence the coherent part of the flow is prescribed as ui<(x)= ∫ Ω Gij(x,x ′)uj(x ′) dx′ (3.4) If we replace u< by the truncated series of eigenfunctions, the above equation can be rewritten N∑ k=1 akφik(x)= ∞∑ p=1 ∫ Ω Gij(x,x ′)apφjp(x ′) dx′ (3.5) Below, it is proved that the filtering function is given by the formula Gij(x,x ′)= N∑ k=1 φik(x)φ ∗ jk(x ′) (3.6) After substituting (3.6) into (3.5) we obtain N∑ k=1 akφik(x)= ∞∑ p=1 N∑ k=1 apφik(x) ∫ Ω φjp(x ′) φ∗jk(x ′) dx′ (3.7) which, fromthenormalisation condition for eigenfunctions, Eq. (2.2), gives the identity. Hence, we have shown that the POD method provides a particular type of filter (3.6); the result of filtering of a given vector Ai(x, t) will be denoted by Ãi(x, t). In the case of turbulent channel flow, homogeneous in the streamwise and spanwise directions, because of the relation between ψ and φ, Eq. (2.5), the filter Gij is given by the formula Gij(∆x,y,y ′,∆z)= (3.8) = 1 LxLz N∑ n=1 Kx∑ kx=−Kx Kz∑ kz=−Kz ψ (n) ikxkz (y)ψ (n) jkxkz (y′)exp ( i∆x 2πkx Lx +i∆z 2πkz Lz ) Modelling of near-wall turbulence... 711 In this case, Gij is a tensor function of 4 variables; this formmight be confu- sing andnot easy to interpret at first sight ifwe refer to standardhomogeneous filters used in LESmethods (Pope, 2000). Hence, we have found it instructive to compare thePODexpansionwith the truncated series expansion of velocity in theChebyshev space in the y direction and theFourier space in the x and z directions. Such an expansion is used in spectral DNS andLESmethods (Pey- ret, 2002). As an example, for a non-homogeneous turbulent channel flowwith non-moving walls at y = H and y =−H the truncated velocity is expressed as a sum ũi(x,y,z) = (3.9) = 1√ LxLz N∑ n=1 Kx∑ kx=−Kx Kz∑ kz=−Kz aTinkxkz(t)fn(y)exp ( ix 2πkx Lx +iz 2πkz Lz ) where in the above, the time-dependent coefficients aTinkxkz(t) are vectors (i = 1,2,3), while fn(y) are functions constructed from the Chebyshev po- lynomials of order 2n, shifted in a way which assures that velocity at the wall is zero, i.e., fn(y) = T2n(y/H)− 1. The function fn is the same for all components of velocity and all kx, kz. On the contrary, in POD series, cf. Eq. (2.5), ankxkz(t) is a scalar, identical for all components of velocity, while the POD eigenfunctions ψ (n) ikxkz (y) are vectors. The eigenvalue problem is solved for each kx, kz pair separately which assures that the basis optimally repre- sents the energy contained in each mode and in each component of velocity. Consequently, the filter function in the POD method is a tensor, while the filter function for the truncated Chebyshev-Fourier expansion reads G(∆x,y,y′,∆z)= (3.10) = 1√ LxLz N∑ n=1 Kx∑ kx=−Kx Kz∑ kz=−Kz fn(y)fn(y ′)exp ( i∆x 2πkx Lx +i∆z 2πkz Lz ) In Section 4, we present selected cross-sections of the PODfilter based on the computed POD modes and compare it with the filter constructed from the Chebyshev polynomials. 3.3. Filtered equations When the POD filtering procedure is applied to Eq. (3.2), we obtain the go- verning equation for the resolved part of the fluctuating velocity field in the near-wall region 712 M. Wacławczyk, J. Pozorski ∂ũi ∂t + 〈U1〉 ∂ũi ∂x + ∂〈U1〉 ∂y ũ2δi1+ ∂ũiũj ∂xj − ∂〈ũiũj〉 ∂xj = (3.11) = [∂ũiuj ∂xj − ∂〈ũiuj〉 ∂xj − ∂ũiũj ∂xj + ∂〈ũiũj〉 ∂xj ] − 1 ̺ ∂p̃ ∂xi +ν ∂2ũi ∂xk∂xk Thebracketed termon theRHSwill be replacedherebyamodel (cf.Holmes et al., 1996), analogously to the subgrid viscosity closure used in LES (cf. Pope, 2000) [·] =−2α1νT ∂s̃ij ∂xj + 4 3 α2(l>) 2 ∂ ∂xi [s̃kls̃kl−〈s̃kls̃kl〉] =−α1νT ∂2ũi ∂xj∂xj + ∂pT ∂xi (3.12) where s̃ij is the resolved strain rate tensor: s̃ij = (ũi,j + ũj,i)/2, α1 and α2 aremodel constants, and l> is the length scale characteristic of the unresolved motion. The second term on theRHS of Eq. (3.12) involves pT which is called the pseudo-pressure; the term can be added to the pressure and the two form together the modified kinematic pressure P = p̃/̺+pT . In the PODapproach, the subgrid viscosity νT is modelled in terms of the first neglected mode νT = u>l> = Ly∫ 0 〈uk>uk>〉 dy · √√√√√Ly Ly∫ 0 〈∂ui< ∂xj ∂ui< ∂xj 〉 dy (3.13) where u> is the velocity scale of unresolved motion reconstructed from the first neglected mode, and l> is the length scale of that mode. A justification for such a choice is presented in Holmes et al. (1996). Filtered equation (3.2), including themodel (3.12) and formula for the mean velocity (3.3) now reads ∂ũi ∂t =− [1 ν y∫ 0 〈ũ1ũ2〉 dy+ u2 ∗ ν ( y− y2 2H )]∂ũi ∂x + −ũ2δi1 [1 ν 〈ũ1ũ2〉+ u2 ∗ ν ( 1− y H )] + (3.14) +(ν +α1νT) ∂2ũi ∂xk∂xk − [∂ũiũj ∂xj − ∂〈ũiũj〉 ∂xj ] − ∂P ∂xi The coherent part of fluctuations can be replaced by truncated series (2.5) ui<(x,y,z)= ũi(x,y,z)= (3.15) = 1√ LxLz N∑ n=1 Kx∑ kx=−Kx Kz∑ kz=−Kz ankxkz(t)ψ (n) ikxkz (y)exp ( ix 2πkx Lx +iz 2πkz Lz ) Modelling of near-wall turbulence... 713 From now on, the notation (̃·) is skipped for the sake of brevity. In order to derive the system of ordinary differential equations for the coefficients ankxkz(t), the Fourier transform is applied to Eq. (3.14). This is followed by the Galerkin projection procedure, i.e. the equation is multiplied by the eigenfunction ψ (m) ikxkz (y) and integrated over y keeping inmind norma- lisation condition (2.2). It is noted that in (2.2) the summation over index i is performed, hence, the resulting ODEs do not depend on i. This procedure leads to the set of ordinary differential equations for the coefficients ankxkz (Holmes et al., 1996) dankxkz dt = A (mn) 1 (kx,kz)amkxkz +(ν +α1νT)A (mn) 2 (kx,kz)amkxkz + +B (lmqn) k′ x k′ z (kx,kz)alkxkzamk′xk′za ∗ qk′ x k′ z +C(t)+ (3.16) +D (lmn) k′ x k′ z (kx,kz)alk′ x k′ z am(kx−k′x)(kz−k′z) The coefficients A (mn) 1 and B (lmqn) k′ x k′ z represent the effect of the mean velocity profile on eigenmodes, (ν +α1νT)A (mn) 2 account for the dissipation of kine- tic energy. The influence of the unresolved scales on the system is modelled via νT ; moreover, D (lmn) k′ x k′ z , and B (lmqn) k′ x k′ z represent interactions between the ba- sis functions. Finally, the term C(t) represents communication between the inner and the outer region via pressure interactions which can bemodelled as impulsive random perturbations (Holmes et al., 1996). By solving the set of equations (3.16), one obtains time evolution of the coefficients ankxkz. 4. POD eigenfunctions Two-point correlations, needed for the PODmethod are computed fromDNS data at Reτ =140 (Podvin, priv. comm.). The data contain the Fourier coef- ficients ûi(kx,y,kz) of 200 independent snapshots of the fluctuating velocity field. The velocity was Fourier-transformed in the streamwise and spanwise directions where the flow is assumed to be periodic. The first wavenumbers, kx = 1, kz = 1 correspond to separations Lx = 600 in the streamwise and Lz =300 in the spanwise direction. The data file includes spanwise wavenum- bers kz ranging from −5 to 5 and streamwise numbers kx = 0,1,2. There are 20 computational points in the wall-normal direction, up to y+ = 50. The current truncation of available coefficients to Kz =5and Kx =2 is quite drastic, therefore only the velocity field of the largestmodes can be computed. Nevertheless, the data are sufficient for the simple POD simulations, with a 714 M. Wacławczyk, J. Pozorski limited number of modes, like the 10D or 32D models of Aubry et al. (1988), Shangi and Aubry (1993). The Fourier transform of the two-point correlation tensor is computed as an ensemble average over N =200 realisations Φij(kx,y,y ′,kz)= 1 N N∑ n=1 û (n) i (kx,y,kz)û (n)∗ j (kx,y ′,kz) (4.1) Additionally, the number of data can be increased by the use of symmetries (cf. Omurtag and Sirovich, 1999). As an example, it follows from the spanwise reflection symmetry that if the functions u1(x,y,z), u2(x,y,z), u3(x,y,z) sa- tisfy theNavier-Stokes equation, then u1(x,y,−z), u2(x,y,−z),−u3(x,y,−z) are also a solution to this equation. The spanwise reflection symmetry results in the following relation for the two-point correlation tensor Φij(kx,y,y ′,kz)= ωijΦij(kx,y,y ′,−kz) (4.2) where ωij = 1 if i and j are both equal to 3 or both different from 3, otherwise ωij = −1. Important for further simulations is the property of the two-point correlation tensor in the case of zero streamwise separations (kx = 0). As the velocity field is real, its Fourier transform must satisfy the relation ûi(kx,y,kz) = û ∗ i(−kx,y,−kz) which, together with symmetry (4.2) written for kx =0, imply that components of the eigenfunction vector ψ1,0,kz and ψ2,0,kz are purely real and ψ3,0,kz is purely imaginary. Given the statistics Φij(kx,y,y ′,kz), eigenvalueproblem(2.6) canbesolved for eachpair of kx and kz.However, as current simulationsare restricted to2D, it is sufficient to compute the eigenvalues and eigenvectors of Φij(0,y,y ′,kz) for kz =1, . . . ,5. In practice, only the two components ψ (n) 1,0kz and ψ (n) 3,0kz are found from eigenvalue problem (2.6), whereas the second component ψ (n) 2,0kz is computed from the continuity equation (cf. Aubry et al., 1988) 2πikzψ (n) 3,0kz + dψ (n) 2,0kz dy =0 (4.3) For this purpose, for each kz the above formula is integrated over y. In each eigenvalue problem (for kz =1, . . . ,5), the first two eigenfunctions and eigenvectors (i.e. for the indices n = 1 and 2) are computed. Here, the work differs from the approaches of Aubry et al. (1988) and Podvin (2001), which were restricted to the first eigenfunction. In the approach of Aubry et al., 10 ordinary differential equations are solved for the real and imaginary parts of the coefficients a1kz, kz =1, . . . ,5. Hence, themodel is called 10D. In the current approach the simulations of 20Dmodel are performed and twenty Modelling of near-wall turbulence... 715 equations for the real and imaginary parts of a1kz and a2kz, kz =1, . . . ,5 are solved. As will be shown further, the second eigenfunction (n = 2) improves the reconstructed velocity statistics of the wall-normal and spanwise compo- nent. The resulting eigenfunctions are presented inFig.1.Their corresponding eigenvalues (which are the energy ofmodes) are presented inTable 1 as a frac- tion of the total kinetic energy of the fluctuating field. As could be expected, the eigenfunction for n =1 and kz = 2 is the most energetic one. The same observation was made by Aubry et al. (1988) and Podvin (2001). For n = 2 we find that the mode kz =1 is the most energetic. It is also seen in Table 1 that the eigenmodes n = 2 are one order of magnitude smaller than those for n =1. Fig. 1. Eigenfunctions for kx =0, five different spanwise wavenumbers (from top to bottom) kz =1, . . . ,5, n =1; (a) ℜ[ψ(1)1,kz], (b) ℜ[ψ (1) 2,kz ]: (——), ℑ[ψ(1)3,kz]: (– – –), and n =2; (c) ℜ[ψ(2)1,kz], (d) ℜ[ψ (2) 2,kz ]: (——), ℑ[ψ(2)3,kz]: (– – –) Thefilter Gij(0,y,y ′,∆z), cf. Eq. (3.8), can nowbe reconstructed fromthe computed eigenfunctions. The cross-sections of the trace Gii(0,y,y ′,∆z) = Gii(y,y ′,∆z) for different values of y = y′ are presented in Fig.2a. 716 M. Wacławczyk, J. Pozorski Table 1. The energy of selected modes as a fraction of the total kinetic energy of the fluctuating field kz 1 2 3 4 5 n =1 7.99% 15.07% 7.38% 0.74% 0.94% n =2 1.27% 0.35% 0.22% 0.1% 0.05% In the z direction, the eigenfunction basis becomes the Fourier basis and the POD method simply provides a sharp spectral cut-off filter. Next, the separation ∆z was set to 0 and the trace of the filter Gii(y,y0,0) was plot- ted for a few selected values of y0, see Fig.2b. The functions attain zero at y = 0. To better illustrate the shape of the function Gii(y,y ′,0) its isolines are plotted in Fig.3. There exists a similarity between this function and the filter constructed from theChebyshev polynomials, cf. Eq. (3.10). The isolines of such a filter, with the truncation at N =6 are presented in Fig.4. In order to simplify the comparisonwith the PODfilter, the variable of theChebyshev polynomials has been appropriately shifted, so that the wall is now placed at y =0. Hence, we have shown that although the form of the POD filter func- tion, Eq. (3.8), may seem complicated at the first sight, it is in fact not far fromthefilter that follows fromthe truncatedFourier-Chebyshev expansion of velocity. Fig. 2. Cross-sections of the trace Gii(y,y,∆z) for (a) y + =2.6: (——), 5.7: (−− · −− · −−), 10.7: (– – –), (b) cross-sections of the trace Gii(y,y0,0) for y+0 =2.6: (——), 10.7: (−− · −− · −−), 30.5: (– – –) From the eigenfunctions and eigenvectors, the velocity statistics can also be reconstructed. As only a few eigenmodes are used, the reconstructed ve- locity field contains less energy than the DNS field. The statistics of velocity components are computed from the formula 〈uiuj〉= N∑ n=1 5∑ kz=1 λnkzψ (n) i,kz ψ (n)∗ j,kz (4.4) Modelling of near-wall turbulence... 717 Fig. 3. Isolines of the trace Gii(y,y ′,0) Fig. 4. Isolines of a filter G(y,y′,0) constructed from the shifted Chebyshev polynomials where i =1,2,3 and N is the number of eigenmodes used in the summation (N =1or N =2). Figure 5a presents the r.m.s. of fluctuating velocity for the streamwise component and the Reynolds-stress component 〈uv〉, the r.m.s. of wall-normal and spanwise components are presented in Fig.5b. The statistics reconstructed for N = 1 and N = 2 are compared with the DNS data of Iwamoto (2002). It is seen that the second eigenmode n = 2 improves the wall-normal and spanwise statistics, while the r.m.s. of streamwise component remains almost unchanged. The profiles presented in Fig.5a compare reaso- nably with the DNS data. In the region very close to the wall (y/H < 0.1 or y+ < 15), the comparison of wall-normal and spanwise fluctuation variance with DNS, cf. Fig.5b, is also acceptable, if one takes into account the simpli- city of the model used in the present study. In the wall vicinity, the flow is dominated by large eddies which form streamwise rolls. Further from thewall, the range of length scales of eddies increases. Hence, the POD method with a few eigenmodes is not able to accurately reproduce the turbulence statistics further away from the wall. 718 M. Wacławczyk, J. Pozorski Fig. 5. Turbulent channel flow at Reτ =140. Velocity statistics reconstructed from the POD eigenfunctions with n =1: (– – –) and n =2: (——), compared with the DNS data of Iwamoto (2002): (a) u+r.m.s.: •, −〈u+v+〉: , (b) v+r.m.s.: •, w+r.m.s.: 5. Results of POD simulation TheNavier-Stokes equations can now be projected on the eigenfunction basis which results in a set of ordinary differential equations for the time coeffi- cients ankz of the 20D model (n = 1,2; kz = 1, . . . ,5). As the evaluation of coefficients requires knowledge of the first and second derivatives of eigenfunc- tions, the computed profiles of Ψ (n) ikz (y) must be curve-fitted, especially in the near-wall region. For this purpose, the least-squares algorithm is used. The streamwise and spanwise components Ψ (n) ikz (y) (i =1 or 3) are fitted to curves aiy+biy 2+ciy 3+diy 4 and thewall-normal component Ψ (n) 2kz (y) is fit to a curve b2y 2 + c2y 3 + d2y 4, so that the near-wall scalings u,w ∼ y and v ∼ y2 are satisfied.Wehave found in the numerical tests that 4th-order polynomials are a good compromise in the interpolation process. The next step is to estimate the influence of the unresolved modes by computing the subgrid viscosity. Here, the approach of Holmes et al. (1996) is followed and νT is estimated from the statistics of the first unresolved mode, cf. Eq. (3.13). The resulting value of νT equals approximately 8.6. The solution to the set of ODEs differs with the varying parameter α. For large values of α (α > 10), the coefficients attain constant non-zero values. In the language of the dynamical system theory, this kind of solution is cal- led stable fixed point. As noticed by Aubry et al. (1988), in such a case the reconstructed velocity field forms two steady counter-rotating rolls, similar to those obtained in the pioneering work of Bakewell and Lumley (1969). The time evolution of the real parts of selected coefficients for α =15 is presented in Fig.6a, where Xnkz =ℜ[ankz]. Modelling of near-wall turbulence... 719 Fig. 6. Selected coefficients ankz computed from the set of ODEs in the 20Dmodel for (a) α =15, (b) α =0.15 Ifwe comparedynamical behaviour of the10Dmodel ofAubryet al. (1988) with the current 20Dmodel we find that the former is less noisy. For example, the intermittent region or steady periodic motion (for details see Aubry et al., 1988) was not found in the present approach. Instead, all coefficients for values of α between 0 and 10 fall in a more complex, chaotic region. The more chaotic time behaviour of the coefficients was also observed in the three- dimensional model of 32 modes, which was considered in the paper of Sanghi and Aubry (1993). Figure 6b presents the time behaviour of selected POD coefficients for α =0.15. A doubt follows from the fact that the solution is in fact sensitive to small changes in values of the coefficients. The point is that in order to compute the coefficients, in particular Amn2 (kx,kz), it is necessary to differentiate twice the eigenfunctions. It is recalled that the two-point correlations are computed from a finite number of snapshots, then the eigenvalue problem is solved, the eigenvectors are curve-fitted, the second component ψ (n) 2,kz is computed from the continuity equation, then the first and second order derivatives of eigen- functions are evaluated. Thewhole procedure introduces a numerical error. So far, the parameter α was chosen to be α =1.5. The resulting time behaviour of the coefficients is presented in Fig.7. The real Xnkz and imaginary parts Ynkz of the coefficients are rescaled, so that |ankz|= X2nkz +Y 2 nkz =1. For the chosen α = 1.5, we obtain the lowest value of the time scale associated with the most energetic mode (kz = 2, n = 1), equal approximately 250 viscous units. A lower value of the time scale is associated with the mode (kz = 1, n =1) and even lower with themodes kz =4and kz =5.The decrease of the time scale with the increasing wavenumber kz is also observed in the analysis of coefficients obtained from the direct projection of the DNS velocity field onto the POD modes (cf. Podvin, 2001). In the present dynamical system, the only exception is the thirdmodewhich has the largest characteristic time 720 M. Wacławczyk, J. Pozorski scale, equal 1300 approximately. Results for modes kz = 1, 2, 4 and 5 are comparable with the experimental data (cf. discussion in the paper of Kasagi et al. (1989) which show a scatter of interburst durations in the range 50-150 (in viscous units). Fig. 7. Real parts of coefficients a1kz computed from the set of ODEs in the 20D model for α =1.5 Thespatial information contained in thePODeigenfunctions togetherwith the time behaviour of POD coefficients obtained from the solution of low- order dynamical system, give a fluctuating velocity field in the near-wall zone, cf. Eq. (3.15). The field reconstructed at a given instant t0 is presented in Fig.8. In the current two-dimensional model, only modes with no streamwise dependence (kx =0)are retained.Thesemodesare associatedwith the streaks elongated in the streamwise direction that are observed in thewall turbulence (Omurtag andSirovich, 1999). In the cross-section of the flow, those structures formpairs of counter-rotating rolls which can also be identified inFig.8. Their characteristic spanwise length scale is around 100 viscous units,whichhas also been observed in other POD studies (e.g., Aubry et al., 1988). Fig. 8. The reconstructed velocity field for a given instant t0, cf. Eq. (3.15) Modelling of near-wall turbulence... 721 6. Conclusions and perspectives In thepaper, thePODapproachwasused toperformsimulations of large-eddy dynamics in the near-wall region of a turbulent channel flow. The first step was to compute a two-point velocity correlation tensor from 200 independent snapshots of DNS data. Next, the eigenvalue and eigenfunction problem was solved to obtain an empirical basis that optimally represents kinetic energy in the considered flow domain. The contribution of the present work is also the derivation and discussion concerning the form of the filter function associated with the POD method. The POD filter constructed from the eigenfunctions is generally a tensor. Additionally, each of its components is a function of 8 variables; however, in the case of homogeneity the number of variables is reduced. The form of the filter might be difficult to analyse at first sight, if one compares it with the standard filters used in the LESmethod on uniform grids. However, as it was shown, the shape of the trace Gii is similar to the filter function associated with the Fourier-Chebyshev expansion of velocity which is used in spectralmethods (Peyret, 2001). The tensor formof thePOD filter follows from the fact that POD eigenfunctions are vectors, while in the Chebyshev expansion the same basis is used for each component of velocity. ThePODsimulations are based on the experimental data of two-point cor- relations, however somenew informationabout thevelocity field canbeextrac- ted, like, e.g., the time evolution of large eddies. To study the flow dynamics, theNavier-Stokes equationswereprojected on the empirical basis to forma set of ordinary differential equations for time-dependent coefficients. Contrary to theworkofAubryet al. (1988)whoconsidered thebasis of 5 eigenmodes in the z direction andonly one eigenmode in the y direction, in thepresentpaper one additionalmode in the y directionwas considered.Hence, the total number of ODEs solved for real and imaginary parts of the time-dependent POD coeffi- cients was 20. The inclusion of additional modes resulted in a randomisation of time behaviour of coefficients in comparison to the 10D model of Aubry et al. (1988). The flow-field resulting from the simulations was two-dimensional, in the plane perpendicular to the main flow direction (y-z plane). This is a result of assumption that the quantities describing the flow are slowly-varying in the x direction in comparison to the variations in the y and z directions. The additional mode improves wall-normal and spanwise components of re- constructed velocity, which is particularly important for further applications like studies on the passive contaminant: in the y-z plane the passive contami- nant is traced in the velocity field of those two components (Wacławczyk and Pozorski, 2004). Analysis of turbulent flows with a passive scalar and two-phase dispersed flows in the near-wall geometry are promising directions for future investiga- 722 M. Wacławczyk, J. Pozorski tions. The low numerical cost of a dynamical system based on POD modes allows, e.g. for a detailed study of the passive contaminant field at varying Schmidt/Prandtl numbers. In the case of two-phase flows, a perspective for a future work is examination of preferential particle concentration and their separation at walls. Acknowledgements The work reported in the present paper has been supported by the IMP PAN statutory funds (O1/Z2/T1)andbythe researchprojectEDF/4300013122(Electricité deFrance,Chatou) coordinatedbyJ.-P.Minier.Moreover, the authorswish to express their gratitude to Bérengère Podvin, D.Sc., for making her DNS results of velocity available. References 1. Allery C., Béghein C., Hamdouni A., Verdon N., 2006, Particle disper- sion in turbulent flowsbyPODlowordermodel usingLES snapshots, In:Direct and Large-Eddy Simulation VI, E. Lamballais, R. Friedrich, B.J. Geurts and O.Mtais (Edit.), 755-762, Springer, the Netherlands 2. Aubry N., Holmes P., Lumley J.L., Stone E., 1988, The dynamics of coherent structures in the wall region of turbulent boundary layer, J. Fluid Mech., 192, 115-173 3. Bakewell H.P., Lumley J.L., 1967, Viscous sublayer and adjacent wall region in turbulent pipe flow,Phys. Fluids, 10, 1880-1889 4. Ball K., Sirovich L., Keefe L., 1991, Dynamical eigenfunction decompo- sition of turbulent channel flow, Int. J. Num. Meth. Phys., 12, 585-604 5. Berkooz G., Holmes P., Lumley J.L., 1993, The proper orthogonal de- composition in the analysis of turbulent flows, Annu. Rev. Fluid Mech., 25, 539-571 6. Citriniti J., George W., 2000, Reconstruction of the global velocity field in the axisymmetric mixing layer utilizing the proper orthogonal decomposition, J. Fluid Mech., 418, 137-166 7. Deville J., Bonnet J.P., 2002, Two-point correlation in fluid dynamics: POD, LSE and relatedmethods, In: Lecture Series onPost-processing of expe- rimental and numerical data, vonKarmanInstitute forFluidDynamics,Rhode- Saint-Genèse, Belgium 8. Deville J., Lamballais E., Bonnet J.P., 2000, POD, LODS and LSE: their links to control and simulations of mixing layers, ERCOFTAC Bulletin, 46, 29-44 Modelling of near-wall turbulence... 723 9. Gunes H., Rist U., 2004, Proper orthogonal decomposition reconstruction of a transitional boundary layer with and without control, Phys. Fluids, 16, 2763-2784 10. Holmes P., Lumley J.L., Berkooz G., 1996,Turbulence, Coherent Struc- tures, Dynamical Systems and Symmetry, Cambridge University Press 11. Iwamoto K., 2002, Database of Fully Developed Channel Flow, THTLAB Internal Report, No. ILR-0201, Dept. of Mech. Eng., The University of Tokyo 12. Joia I.A., Gobeau N., Ushijima T., Perkins R.J., 1998, POD study of bubble and particlemotion in turbulent channel flow, International Conference on Multiphase Flow, June 8-12, Lyon, France 13. Kasagi N., Kuroda A., Hirata M., 1989, Numerical investigations of near- wall turbulent heat transfer taking into account the unsteady heat conduction in the solid wall, J. Heat Transfer, 111, 385-392 14. ŁuniewskiM., J. Pozorski,T.Wacławczyk, 2006,LESof turbulent chan- nel flowwith dispersed particles [in Polish], Systems, 11, 224-232 15. MoinP.,MoserR.D., 1989,Characteristic-eddydecompositionof turbulence in a channel, J. Fluid Mech., 200, 471-508 16. Omurtag A., Sirovich L., 1999, On low-dimensional modelling of channel turbulence,Theoret. Comput. Fluid Dynamics, 13, 115-127 17. PeyretR., 2002,SpectralMethods for Incompressible Viscous Flow, Springer- Verlag, NewYork 18. Picciotto M., Marchioli C., Soldati A., 2005, Characterization of near- wall accumulation regions for inertial particles in turbulent boundary layers, Phys. Fluids, 17, 098101 19. Podvin B., 2001, On the adequacy of the ten-dimensional model for the wall layer,Phys. Fluids, 13, 210-224 20. Pope S.B., 2000,Turbulent Flows, Cambridge University Press 21. Sanghi S.,AubryN., 1993,Mode interactionmodels for near-wall turbulence, J. Fluid Mech., 244, 455-488 22. Wacławczyk M., Pozorski J., 2002, Two-point velocity statistics and the POD analysis of the near-wall region in a turbulent channel flow, J. Theor. Appl. Mech., 40, 895-916 23. Wacławczyk M., Pozorski J., 2004, Conjugate heat transfer modelling using the FDF approach for near-wall scalar transport coupled with the POD method for flow dynamics, In: Advances in Turbulence X, H.I. Andersson and P.-A. Krogstad (Edit.), CIMNE, Barcelona 24. Webber G., Handler R., Sirovich L., 1997, The Karhunen-Loéve decom- position of minimal channel flow,Phys. Fluids, 9, 1054-1066 724 M. Wacławczyk, J. Pozorski Modelowanie turbulencji za pomocą wielkoskalowych modów prędkości Streszczenie Przedmiotem pracy jest modelowanie turbulentnego pola prędkości w obszarze przyściennym za pomocą niskowymiarowego systemu dynamicznego, opartego o de- kompozycję w bazie funkcji własnych POD (ang.Proper Orthogonal Decomposition). Empiryczna baza funkcyjnaPODzostaławyznaczona z rozwiązania zagadnieniawła- snego, w którym obecne są dwupunktowe korelacje prędkości. Następnie, w wyniku projekcji Galerkina równań pędu na podprzestrzeń rozpiętą na tej bazie funkcyjnej, otrzymanoukład równań różniczkowychzwyczajnychna zależne od czasuwspołczyn- niki. Na podstawie funkcji własnych oraz z wyznaczonych współczynników rozkładu uzyskano ewolucję w czasie charakterystycznych struktur wirowychw obszarze przy- ściennym. Systemdynamiczny rozpatrywanywpracy składa się z 20 równań różnicz- kowych zwyczajnych. Zrekonstruowane pole prędkości jest dwuwymiarowe (w płasz- czyźnie prostopadłej do głównego kierunku przepływu). Ponadto w pracy dyskuto- wana jest procedura filtrowania związana z metodą POD. Wyprowadzony filtr POD porównano z formułą używaną wmetodzie symulacji dużych wirów. Manuscript received January 22, 2007; accepted for print May 11, 2007