PHYSICS OF FLUIDS 17, 025105 s2005d A scale-dependent Lagrangian dynamic model for large eddy simulation of complex turbulent flows Elie Bou-Zeida! Department of Geography and Environmental Engineering and Center for Environmental and Applied Fluid Mechanics, Johns Hopkins University, 313 Ames Hall, 3400 North Charles Street, Baltimore, Maryland 21218 Charles Meneveaub! Department of Mechanical Engineering and Center for Environmental and Applied Fluid Mechanics, Johns Hopkins University, 127 Latrobe Hall, 3400 North Charles Street, Baltimore, Maryland 21218 Marc Parlangec! School of Architecture, Civil, and Environmental Engineering, Swiss Federal Institute of Technology at Lausanne, Building GC-Ecublens, CH-1015 Lausanne, Switzerland sReceived 25 August 2004; accepted 26 October 2004; published online 19 January 2005d A scale-dependent dynamic subgrid model based on Lagrangian time averaging is proposed and tested in large eddy simulations sLESd of high-Reynolds number boundary layer flows over homogeneous and heterogeneous rough surfaces. The model is based on the Lagrangian dynamic Smagorinsky model in which required averages are accumulated in time, following fluid trajectories of the resolved velocity field. The model allows for scale dependence of the coefficient by including a second test-filtering operation to determine how the coefficient changes as a function of scale. The model also uses the empirical observation that when scale dependence occurs ssuch as when the filter scale approaches the limits of the inertial ranged, the classic dynamic model yields the coefficient value appropriate for the test-filter scale. Validation tests in LES of high Reynolds number, rough wall, boundary layer flow are performed at various resolutions. Results are compared with other eddy-viscosity subgrid-scale models. Unlike the Smagorinsky–Lilly model with wall-damping swhich is overdissipatived or the scale-invariant dynamic model swhich is underdissipatived, the scale-dependent Lagrangian dynamic model is shown to have good dissipation characteristics. The model is also tested against detailed atmospheric boundary layer data that include measurements of the response of the flow to abrupt transitions in wall roughness. For such flows over variable surfaces, the plane-averaged version of the dynamic model is not appropriate and the Lagrangian averaging is desirable. The simulated wall stress overshoot and relaxation after a jump in surface roughness and the velocity profiles at several downstream distances from the jump are compared to the experimental data. Results show that the dynamic Smagorinsky coefficient close to the wall is very sensitive to the underlying local surface roughness, thus justifying the use of the Lagrangian formulation. In addition, the Lagrangian formulation reproduces experimental data more accurately than the planar-averaged formulation in simulations over heterogeneous rough walls. © 2005 American Institute of Physics. fDOI: 10.1063/1.1839152g sSGSd stress tensor sits traceless partd is modeled according to I. INTRODUCTION Large-eddy simulation sLESd has become an important tool for the study of high-Reynolds number environmental1–8 and engineering9–11 turbulent flows. LES resolves the flow at scales larger than a certain size D, while the smaller scales are parametrized. The classic, most often used parametrization sthe Smagorinsky model12d is based on the concepts of eddy-viscosity and mixing length, in which the subgrid-scale ad Telephone: 1-410-516-5031. Fax: 1-410-516-8996. Electronic mail: eliebz@jhu.edu bd Telephone: 1-410-516-7082. Fax: 1-410-516-7254. Electronic mail: meneveau@jhu.edu cd Also at Department of Geography and Environmental Engineering and Center for Environmental and Applied Fluid Mechanics, Johns Hopkins University, Baltimore, MD 21218. Telephone: 141-21-693-6391. Electronic mail: Marc.Parlange@epfl.ch 1070-6631/2005/17~2!/025105/18/$22.50 = − 2nTS̃ij = − 2scs,DDd2uS̃uS̃ij . tSMAG ij s1d Above, nT is the eddy viscosity, S̃ij = 0.5s] jũi + ]iũ jd is the resolved strain rate tensor swhere ũi is the resolved velocity fieldd, and the strain-rate magnitude is given by uS̃u Î = 2S̃ijS̃ij. The only undetermined parameter in the above expression is the Smagorinsky coefficient cs,D. Even though many SGS models that do not rely on the viscous analogy concept have been developed, this paper will focus on eddyviscosity models since they continue to be the most widely used in practice, either individually or in conjunction with the similarity type models in the so-called mixed models.13–18 17, 025105-1 © 2005 American Institute of Physics Downloaded 24 Jan 2005 to 128.220.2.42. Redistribution subject to AIP license or copyright, see http://pof.aip.org/pof/copyright.jsp 025105-2 Bou-Zeid, Meneveau, and Parlange In 1991, an important development took place in LES with the introduction of the dynamic model and the Germano identity.19,20 By relating stresses at different scales, the Germano identity allows unknown model coefficients, such as cs,D, to be determined from the smallest resolved scales between the grid-scale D and a test-filter scale aD sa . 1, usually a = 2d. A major assumption of the original dynamic approach is scale-similarity, i.e., that model coefficients are the same at different scales: cs,D = cs,aD ssee the discussion by Meneveau and Katz21d. Scale invariance is a reasonable assumption if D pertains to an idealized inertial range of turbulence, but it is not expected to hold if D falls near a transition scale that separates different physical processes occurring in distinct ranges of scales. One example where D falls near a transition scale occurs when the grid scale D approaches the integral scale, a limit that is of relevance when LES approaches the Reynolds-averaged Navier-Stokes formulation in certain parts of the flow. Of specific interest to applications to be examined in this paper, such a situation occurs in LES of high-Reynolds number wall-bounded flows where the integral scale is on the order of the distance to the wall. At the high-Reynolds numbers that occur in applications to atmospheric or oceanic boundary layers, LES cannot resolve the viscous sublayer sdue to computational power limitationsd. In such applications, the first few cells near the surface have a grid scale on the order of the local integral scale and inaccurate results are obtained from the traditional scale-invariant dynamic model ssee SGS comparison section in this paper and Refs. 22 and 23 for a discussion and illustration of this effectd. Moreover, in this situation the subgrid stresses carry a significant fraction of the total mean momentum fluxes, and hence the LES results are particularly sensitive to the SGS model. To address this shortcoming of the traditional dynamic model, Porté-Agel et al.22 proposed a scale-dependent version of the dynamic model in which a second test filter determines how the coefficient changes across scales, thus providing more accurate estimation of the coefficient at the grid scale. In the tests of atmospheric boundary layer flow over homogeneous surfaces performed in Porté-Agel et al.,22 the scale-dependent dynamic model was implemented using planar averaging, i.e., the averages required to enforce the Germano identity were evaluated over horizontal planes parallel to the ground. This was appropriate for the simple geometries envisioned in those tests, where horizontal planes correspond to directions of statistical homogeneity of the turbulence. An important question is how to treat complexgeometry flows that do not possess directions of statistical homogeneity and thus do not present obvious spatial domains over which to evaluate averages during LES. For the dynamic model, the issue of averaging the terms in the Germano identity has been the subject of considerable research19,24,25 salso see discussion in Pope26d. Especially in the context of eddy-viscosity closures, averaging is crucial to reduce the large amount of noise that is present when no averaging is performed srecall that in LES one must evaluate the divergence of the modeled SGS stress tensor; therefore, unphysical fluctuations in the coefficient can lead to significant errors, not to mention numerical instability if the coef- Phys. Fluids 17, 025105 ~2005! ficient becomes negatived. To reduce the noise for applications in complex geometry flow, time averaging was proposed by Meneveau et al.24 In order to comply with Galilean invariance, time must be considered by following fluid parcels of the flow. Thus the Lagrangian dynamic model was proposed and tested in a number of flows.15,16,24,27–30 In these applications, viscous sublayers were resolved and hence the assumption of scale invariance was justified; the scaleinvariant Lagrangian approach allowed determination of the coefficient in these complex-geometry situations. However, for applications where the viscous sublayer cannot be resolved, such as high-Reynolds number se.g., atmosphericd boundary layers over complex terrain, that model is not applicable due to the lack of scale invariance near the ground. Thus, an important issue remains, namely, the formulation of a scale-dependent Lagrangian dynamic model. This paper is devoted to this task. In Sec. II we review the basic elements of the scale-invariant dynamic model, the scaledependent dynamic formulation with planar averaging of Porté-Agel et al.,22 and the Lagrangian scale-invariant dynamic model of Meneveau et al.24 Section III of the paper presents the proposed Lagrangian scale-dependent dynamic model, enabling applications of the dynamic model to complex-geometry flows without assuming scale-invariance or spatial homogeneity. Section IV describes the numerical code used in this work to simulate high-Reynolds number boundary-layer flows and presents test results in horizontally homogeneous flows for which the performances of the various SGS models described in Secs. II and III can be compared in detail. Section V describes applications of the Lagrangian scale-dependent dynamic model to high-Reynolds number atmospheric flow over rough surfaces with abrupt changes in wall roughness si.e., horizontally nonhomogeneousd, and compares the results to existing field measurements data and to the results obtained with a planar-averaged version of the scale-dependent model. Conclusions are presented in Sec. VI. II. REVIEW OF THE DYNAMIC, SCALE-DEPENDENT, AND LAGRANGIAN DYNAMIC SGS MODELS The original nondynamic Smagorinsky–Lilly model sdenoted SMAG belowd has already been introduced in Sec. I. For isotropic homogeneous turbulence, with D falling in the inertial range, the analysis of Lilly31 yields cs , 0.16 sfor the spectral cutoff filterd, a value that provides good results in LES of idealized isotropic turbulence. It remains to point out that for applications to high-Reynolds number boundary layers in which the viscous sublayer is not resolved ssee Pope26d a wall-damping function needs to be included in the specification of the coefficient, otherwise turbulence generation is excessively damped and insufficient kinetic energy occurs in the resolved scales of the simulation. A classic wall-damping function was proposed by Mason and Thompson32 where the SGS mixing length, l = cs,DD is decreased close to the surface to merge smoothly with the l , z behavior expected there. The resulting damping function is Downloaded 24 Jan 2005 to 128.220.2.42. Redistribution subject to AIP license or copyright, see http://pof.aip.org/pof/copyright.jsp 025105-3 Phys. Fluids 17, 025105 ~2005! A scale-dependent Lagrangian dynamic model 1 1 1 = , + ln ln0 fksz + zodgn s2d where k is the von-Karman constant s<0.4d and l0 = cs,0D is the mixing length away from the wall sin a region of nearly homogeneous isotropic turbulenced. In simulations using the Lilly–Smagorinsky model in this study, cs,0 is taken as 0.16 and a value of 2 is assigned to n, which fixes the damping function shape salternative values of n = 1 and cs,0 = 0.1 were also tested in other studies, see Porté-Agel et al.22 for more detailsd. Despite the use of this wall-damping function swith varying values of cs,0 and nd, the Smagorinsky–Lilly model remains overdissipative22,32 and would require further detailed calibrations of the coefficient to yield more accurate results. To avoid the need for such case-by-case calibration of parameters, the dynamic model of Germano et al.19 was proposed. A. The dynamic model The dynamic model19 consists of using the smallest resolved scales to measure the model coefficient during the simulation. The model is based on a relation between SGS stresses at different scales sthe grid scale D and a test filter scale aD, where a is usually taken as 2d expressed by the following identity: s3d Lij = Tij − s̄ij = ũiũ j − uD iuD j . Here sij is the SGS stress tensor at scale D, Tij is the SGS stress tensor at the test-filter scale aD and Lij is the SGS stress tensor defined from scales intermediate between D and aD. Lij is the resolved stress tensor and can be computed exactly from the resolved velocity field using Eq. s3d. Throughout the paper, a tilde s˜d denotes the filtering operation at the grid-scale D and a bar s¯d denotes test-filtering at the test-filter scale aD, typically a = 2. Later on, a caret sˆd denotes the second test-filtering at a second test-filter scale a2D. Ensemble averaging will be denoted by brackets kl. Brackets followed by dimensions subscripts will denote averaging in all the indicated dimensions; for example, kulx,y is the velocity averaged over the x and y directions. Using the Smagorinsky model to express the deviatoric parts of SGS stresses at the scales D and aD and assuming that the coefficient cs,D does not fluctuate strongly in space to justify extracting it from the test-filtering operation24,25 results in the following expressions: 2 D2uS̃uS̃ij, t̄ij = s̄Dij = − 2cs,D 2 2 D D TD ij = − 2cs,aDsaDd uSuSij . s4d The superscript D denotes the deviatoric strace-freed part of the tensor. Replacing in Eq. s3d yields an error in that identity induced by the use of the Smagorinsky model. This error is 2 D D eij = LD ij − sTij − t̄ijd = Lij − cs,D M ij , s5d where M ij is given by M ij = 2D2fuS̃uS̃ij − a2buSD uSD ijg; 2 b = cs,2 aD / cs,D s6d is a parameter that accounts for possible scale dependence of cs,D. Usually, the use of this model makes the assumption of scale invariance, i.e., cs,D = cs,aD or b = 1. To obtain an optimal value of cs,D, the square error eijeij, is minimized contracted over all tensor terms.33 Nevertheless, the local determination of cs,D at every grid point yields a highly variable coefficient that is numerically unstable, mainly due to a high frequency of negative values. Some averaging is hence needed to stabilize the coefficient. Ghosal25 showed that averaging over homogeneous spatial directions yields a system that is consistent with the extraction of cs,D from the filter operation and that is equivalent to Lilly’s expression averaged over homogeneous directions. Interestingly, Pope also showed that if the coefficient is optimized to minimize sin a least-square sensed the dependence of relevant turbulence statistics on the grid scale D, the expression obtained is again exactly Lilly’s expression.34 In wall-bounded flows, horizontal planes are usually selected as the homogeneous directions for averaging. In the absence of homogeneous directions, the Lagrangian approach24 can be used to average the coefficient over time along fluid pathlines sthe approach is further explained in Sec. II C belowd. However, for any type of averaging, denoted by brackets kl, the Smagorinsky coefficient determined by a least-square error minimization of keijeijl can be written as 2 = cs,D kLij M ijl . kM ij M ijl s7d Note that the contraction of Lij with M ij eliminates the need to distinguish between Lij and LD ij , since M ij is a deviatoric stracelessd tensor in incompressible flows. With planar averaging, this scale-invariant version of dynamic model will henceforth be denoted PASI. B. The planar-averaged scale-dependent dynamic model To account for scale effects in the dynamic model, two approaches exist. If prior knowledge of the variation of cs,D with scale is available, the parameter b can be prescribed a priori. Such a “semi-dynamic” approach with an imposed parameter b was tested by Meneveau and Lund35 to capture scale dependence in the transition from LES to direct numerical simulation sDNSd sD → the Kolmogorov scale hd for finely resolved LES. The approach was also used in simulations by Bou-Zeid et al.36 for wall-bounded flows in the limit where D tends to the local integral scale near the wall. In both studies, the semi-dynamic scale-dependent model was shown to give better results than the scale-invariant model. The other option is to implement a fully dynamic formulation where the scale-dependence parameter b is measured through an additional filtering operation. The latter approach has been successfully implemented with the planar-averaged dynamic approach for atmospheric boundary-layer sABLd flows by Porté-Agel et al.22 The main assumption used in this model is that a power-law behavior describes the scale dependence of the coefficient, i.e., cs,D , Df or, in a dimensionally more appropriate form, cs,D = cs,aDsD / aDdf. As a consequence, b evaluated as the ratio of coefficients at scales Downloaded 24 Jan 2005 to 128.220.2.42. Redistribution subject to AIP license or copyright, see http://pof.aip.org/pof/copyright.jsp 025105-4 Phys. Fluids 17, 025105 ~2005! Bou-Zeid, Meneveau, and Parlange aD and D is equal to b evaluated between scales a2D and aD si.e., the power-law assumption is equivalent to the assumption that the scale-dependence parameter b = a2f is itself scale-invariant—see discussion in Porté-Agel et al.22d. A second Germano identity, written between scales D and a2D yields In the Lagrangian SGS model, the coefficient cs,D is obtained by minimizing the weighted time average of the local error contraction eijeij over pathlines; this weighted time average can be written as E= E t eijfzst8d,t8geijfzst8d,t8gWst − t8ddt8 , s10d −` 2 = cs,D kQijNijl , kNijNijl s8d where Qij sthe resolved stress tensor between D and a2Dd and Nij are given by w ˆ ˆ Qij = ũˆ iũ j − ũiũ j, ˆ ˆ ˆ Nij = 2D fuS̃uS̃ij − a4b2uS̃uS̃ijg. 2 where zst8d are the previous positions of the fluid elements and Wstd is a relaxation function that typically allocates larger weights to the more recent history of the coefficients si.e., Wstd is a decreasing function of td. By filtering between D and aD and using a scale invariant form of the Germano identity, the coefficient is obtained by setting the variation of 2 to zero: E with respect to cs,D ]E 2 = ]cs,D s9d By equating the right-hand sides of Eqs. s7d and s8d, Porté-Agel et al.22 obtained a fifth-order polynomial in b, including ten terms stensor contractions involving various filtered strain-rates and resolved stress tensorsd that need to be averaged over directions of statistical homogeneity.22 The polynomial is solved for b and the solution is then replaced in Eq. s7d to obtain a scale-dependent estimate of cs,D. This implementation of the scale-dependent model, used in PortéAgel et al.,22 worked well with the planar averaging approach. E t 2eij −` ]eij 2 Wst − t8ddt8 = 0, ]cs,D s11d 2 : which results in the following expression for cs,D 2 cs,D = JLM , s12d Lij M ijfzst8d,t8gWst − t8ddt8 s13d M ij M ijfzst8d,t8gWst − t8ddt8 . s14d JMM where JLM = E t −` and C. The Lagrangian-averaged scale-invariant SGS model The requirement for homogeneous directions in the flow field limits the use of the dynamic model to relatively simple flows, excluding many practical flows in complex geometries. Local formulations of the model have been developed.25,37 As outlined in Sec. I, the need to evaluate the divergence of the SGS stress makes the highly intermittent coefficient fields, which often need to be clipped at zero, resulting from purely local dynamic determinations undesirable. Moreover, conceptually some averaging is needed to recover the statistical basis of the eddy-viscosity model.21,24 An alternative approach, combining features from the local and averaged formulations, was developed by Meneveau et al.24 The model averages the Smagorinsky coefficient in time following fluid pathlines and hence it is called the Lagrangian-averaged scale-invariant model sLASId. The Lagrangian averaging enforces to some degree the statistical basis that supports the use of an eddy-diffusion model and is physically justifiable since turbulent eddies with sizes about the grid scale are likely to be convected along fluid pathlines. Meneveau and Lund38 also showed that the turbulent energy cascade is most apparent when viewed in a Lagrangian frame of reference. The model is very well suited for the applications with heterogeneous spatial conditions since it preserves local variability, preserves Galilean invariance, and does not require homogeneous directions. It has already been applied in LES of flows in complex domains such as flows in internal combustion engines,27 flow over wavy walls,28 flows in thrust reversers,29 and flow of impinging jets.30 J MM = E t −` For the weighting function, a choice of an exponential form, Wst − t8d = s1 / Tde−st−t8d/T, allows replacing cumbersome evaluations of backward time integrals with forward relaxation-transport equations. Based on DNS results and dimensional self-consistency,24 the time scale T is chosen as T = 1.5DsJLM JMM d−1/8. This choice of the time scale offers the practical advantage of allocating less weight to recent history si.e., increasing the model’s memoryd if the current values of Lij M ij are negative. The time scale is effectively infinite if JLM reaches zero, thus preventing negative values of cs,D. The relaxation transport equations thus obtained for JLM and JMM are DJLM Dt = ]JLM ]t + ũ · ¹ JLM = 1 sLij M ij − JLM d TD s15d and DJ MM Dt = ]JMM ]t + ũ · ¹ JMM = 1 sM ij M ij − JMM d. TD s16d Using first-order numerical approximations in space and time, these equations can be discretized and included rather economically in an LES code. The resultant formulation to update from time-step “n” to “n + 1” at a grid point located at x is Downloaded 24 Jan 2005 to 128.220.2.42. Redistribution subject to AIP license or copyright, see http://pof.aip.org/pof/copyright.jsp 025105-5 Phys. Fluids 17, 025105 ~2005! A scale-dependent Lagrangian dynamic model n+1 Jn+1 sxd + s1 − «dJnMM sx − ũnDtd, MM sxd = «fM ij M ijg n+1 n JLM sxd = Hh«fLij M ijgn+1sxd + s1 − «dJLM sx − ũnDtdj, these terms relating to the test-filter scale dominate the error term fEq. s5dg. Using this observation and using a = 2, we may write 2 cs,2D = where «= Dt/Tn , 1 + Dt/Tn n Tn = 1.5DsJLM JnMM d−1/8 and Hhxj = ramp function = U s17d Here Dt is the time step. Bilinear spatial interpolation is used to evaluate the previous values at position x − ũnDt, i.e., “upstream” of the grid point in question. The ramp function is needed to ensure that numerical inaccuracies smainly due to the discretization of the equationsd do not yield slightly negative cs,D values despite the choice of the time scale to avoid such occurrences, and to avoid infinities when evaluating Tn. For wall-bounded flow simulations ssuch as LES presented later in this studyd, periodic boundary conditions for JLM and J MM are used in the horizontal directions. At the lower and upper boundaries, zero-gradient shomogeneous Neumannd boundary conditions are imposed, i.e., the values at the boundary are set equal to the values at the closest node inside the domain. s18d . JQN , s19d QijNijfzst8d,t8gWst − t8ddt8 s20d NijNijfzst8d,t8gWst − t8ddt8 . s21d JNN where JQN = E t −` and JNN = E t −` Here Qij and Nij are computed using their definitions given in Eq. s9d and a value of b = 1; however, the weighting function now involves exponential decay with a time constant, n n n −1/8 = 1.5DsJQN JNN d . T4D s22d In the Lagrangian scale-dependent model the averages IQN and INN are evaluated from two additional relaxation transport equations: ]JQN III. SCALE-DEPENDENT LAGRANGIAN DYNAMIC MODEL Implementation of the Lagrangian averaging approach described in Sec. II C with the scale-dependent method described in Sec. II B would require accumulating ten different Lagrangian averages ssee Sec. II B and Porté-Agel et al.22d, as well as solving a fifth-order polynomial for b at every grid point in the domain. This is a prohibitively expensive procedure, and the choice of the proper polynomial root, if it exists, is difficult due to the more noisy characteristics of Lagrangian averaging compared to planar averaging. Hence, a simplified procedure that is better suited for Lagrangian averaging is sought here The procedure is based on the observation that the dynamic model, in its scale-invariant formulation fi.e., using b = 1 when evaluating M ij sEq. s6dg, yields a coefficient that corresponds to the test-filter scale aD rather the grid-filter scale D. That is to say, in practice one approximately obtains a value suitable for cs,aD instead of cs,D when evaluating the right-hand side of Eqs. s7d or s12d with M ij using b = 1. This effect has been observed to be true in numerical experiments of scale dependence both when the coefficient increases with increasing scale, such as when D approaches the Kolmogorov scale,35 and when the coefficient decreases with increasing scale, e.g., when D is close to the integral scale.22 This effect can be traced to the fact that, in the Germano identity fEq. s3dg, the tensor terms relating to the test-filter scale aD are significantly greater than the terms relating to the grid-filter scale D sabout 4 times greater if a = 2d. Hence JMM Similarly, for a = 4, we write 2 cs,4D = if x ù 0 x u. −32 10 otherwise JLM ]t + ũ · ¹ JQN = 1 sQijNij − JQNd T4D s23d + ũ · ¹ JNN = 1 sNijNij − JNNd. T4D s24d and ]JNN ]t When they are discretized according to the same procedures and approximations described in the preceding section, they read n+1 n sxd = «4DfNijNijgn+1sxd + s1 − «4DdJNN sx − ũnDtd, JNN n+1 JQN sxd = Hh«4DfQijNijgn+1sxd + s1 − «4Dd n 3JQN sx − ũnDtdj, where «4D = n Dt/T4D n 1 + Dt/T4D , n n n −1/8 T4D = 1.5DsJQN JNN d and Hhxj = ramp function = U s25d if x ù 0 x u. −32 10 otherwise Boundary conditions are similar to those used for JLM and JMM speriodic in the horizontal directions and zero gradient at the lower and upper boundariesd. In the source terms, both M ij and Nij are evaluated assuming that b = 1 swhich can be considered as a “first guess” in an iterative, explicit Downloaded 24 Jan 2005 to 128.220.2.42. Redistribution subject to AIP license or copyright, see http://pof.aip.org/pof/copyright.jsp 025105-6 Phys. Fluids 17, 025105 ~2005! Bou-Zeid, Meneveau, and Parlange procedure to take scale-dependence into accountd. From the locally known approximate values of cs,2D and cs,4D, the local value of b can be evaluated according to s26d 2 2 b = cs,4D /cs,2D . 22 Next, as in Porté-Agel et al., we use the assumption 2 2 2 2 that b is scale invariant, i.e., b = cs,4D / cs,2D = cs,2D / cs,D swhich amounts to postulating power-law dependence of the coefficient as a function of scale, with arbitrary exponentd. We may then solve for the unknown coefficient at scale D, 2 cs,D = where b is given by Eq. s26d. In principle, b obtained from Eq. s26d can vary between 0 and infinity, depending on the local values obtained for cs,4D and cs,2D. In practice, of course, some of these limits may cause numerical difficulties. We have observed that in the limit of b → ` si.e., when cs,2D tends to zero while cs,4D does notd, the local coefficient used in LES goes to zero ssmoothlyd and this limit does not pose any difficulties in LES. However, when cs,4D tends to zero while cs,2D does not, b → 0, leading to very large values of cs,D from Eq. s27d. In simulations, this may lead to numerical instabilities associated with viscous stability conditions being locally violated. Therefore, some clipping of b away from zero is required to allow simulations to proceed. We choose a lower limit of b ù 0.125. The clipping limit is significantly below physically expected limiting behaviors. As the ground is approached, the mixing length g = Dcs,D is expected to become proportional to the integral scale z. Therefore, close to the ground Dcs,D , z or cs,D , z / D. Therefore, as the wall is approached, the average value of b tends to z2/s2Dd2 = 1/4. z2/D2 s28d Away from the wall, b is expected to increase to a mean value of about 1. Hence, the clipping limit of 0.125 allows the coefficient to decrease locally to half its average minimum value. Note that this clipping limit already corresponds to very severe scale dependence, namely, the eddy-viscosity coefficient at scale D being eight times larger than at scale 2D. No limit is imposed for high b, thus the allowable range is 0.125ø b ø `. Tests with other clipping limits were conducted sb was clipped at 0.1 and 0.167d without significant impact on the results. This was expected since it is observed in simulations that the clipping at 0.125 is only needed about 15% of the time: the clipping limit cannot significantly affect the results or be used as a tuning parameter to adjust the results. Therefore, the model coefficient to be used in the simulation is obtained at every grid point and time step from the Lagrangian averaged quantities according to 2 = cs,D 2 cs,D = s27d 2 cs,2D /b , 2 2 b = cs,2D /cs,D , test-filtering operations to evaluate the local source terms Lij M ij, M ij M ij, QijNij, and NijNij, updating the four Lagrangian averages JLM , JMM , JQN, and JNN following Eqs. s17d and s25d, and evaluating the local coefficient value to be used in LES according to Eq. s29d. The same scale-dependent approach can obviously be used with planar averaging where the planar averaged coefficient at the grid-filter scale can be obtained using 2 cs,2D = maxsb,0.125d S max JLM /JMM JQNJMM JNNJLM D . s29d ,0.125 In summary, the Lagrangian-averaged scale-dependent dynamic model sLASDd involves test-filtering and second kLij M ijl/kM ij M ijl . kQijNijlkM ij M ijl kNijNijlkLij M ijl S D s30d The brackets denote planar averaging. This planar-averaged scale-dependent model will be referred to as PASD. No clipping is needed with the planar averaged approach since the heavy averaging associated with the planar formulation eliminates the fluctuations that led to the need for clipping. IV. NUMERICAL CODE AND TESTS IN HOMOGENEOUS, HIGH-REYNOLDS NUMBER, BOUNDARY LAYER FLOWS A. Numerical code The isothermal LES equations are solved in rotational form to ensure conservation of mass and kinetic energy of the inertial terms:39 ]ũi = 0, ]xi s31d 1 ] p̃* ]ũi ]ũi ]ũ j ] + ũ j =− − s2scs,DDd2uS̃uS̃ijd + F̃i . + ]t ]x j ]xi r ]xi ]x j S D Here Fi is the mean streamwise pressure forcing. Note that in the above equations, the molecular viscous term is neglected because the paper focuses on very high-Reynolds number flows where viscosity is negligible at the resolved scales and the wall layer is modeled sas opposed to resolving the viscous sublayer, see Pope26d. The modified pressure term p̃ * = p̃ + s1/3drskk + s1/2drũ jũ j s32d is computed as usual from a Poisson equation sdivergence of momentum equation set to zero due to continuity conditiond. The code uses a pseudospectral approach in the horizontal directions. A second-order accurate centered-differences scheme, requiring a staggered grid, is used in the vertical direction. This entails storing the variables at heights jdz or sj + 1 / 2ddz; where j goes from 0 to N sthe number of vertical grid pointsd. The fully explicit second-order accurate AdamsBashforth scheme is used for time advancement. Aliasing errors can be detrimental to the accuracy of the SGS parametrization since they affect the smallest resolved scales used to compute the dynamic Smagorinsky coefficient. To overcome this problem, the 3 / 2 rule40 is used to fully dealias the convective terms. More details about some numerical aspects of the code sunrelated to the SGS modeld can be found in Refs. 7 and 8. Downloaded 24 Jan 2005 to 128.220.2.42. Redistribution subject to AIP license or copyright, see http://pof.aip.org/pof/copyright.jsp 025105-7 Phys. Fluids 17, 025105 ~2005! A scale-dependent Lagrangian dynamic model TABLE I. Simulation parameters for homogeneous surface simulations. Domain height Domain length L and width W H L = W = 2pH Mesh spacings dx = dy = 2p dz Wall roughness z0 imposed in lower boundary condition z0 / H = 10−4 Number of grid points N3 = 643 Initial conditions Mean velocity: logarithmic profile near the surface merging smoothly with a zero-gradient profile at the top. Velocity fluctuations: imposed randomly sin space and among componentsd on the mean profile using a prescribed turbulent kinetic energy profile sfrom results by Andren et al.ad. Warm-up period Warm-up simulations are run until the normalized total stress profile adjusts to a straight line reaching 1 at the surface and the mean resolved kinetic energy is stable. Forcing Imposed pressure gradient: s1 / rd = p = u*2 / H. Simulation time step, dt dtu* / H = 0.000 25 snondimensional time unitsd. Lagrangian model time step Number of simulation time steps Dt = 5dt 200 000 Output sampling frequency Every ten time steps a Reference 43. The boundary conditions in the horizontal directions are periodic. A stress free condition is imposed at the top of the domain by setting s33d ]3ũ1,2 = ũ3 = 0, where 1, 2, and 3 sor x, y, z in other parts of the paperd refer to the streamwise, cross-stream, and vertical directions, respectively. At the bottom of the domain, the vertical velocity is set to 0 at the surface. As a consequence of the staggered grid formulation, no boundary conditions are needed for the horizontal velocities since they are only stored at a distance dz / 2 above the surface. Stresses at the surface are imposed through a local similarity theory formulation41 ssee Piomelli and Balaras42 for a review of wall modeling in LESd. However, velocities filtered at twice the grid scale are used to compute the surface stress; this is needed to ensure that the average stress over the wall is close to the stress predicted by the classic log law. The need for this formulation and its derivation are explained in the Appendix. The resulting lawof-the-wall formulation is twsx,yd = − F k lnfsdz/2d/zog G 2 sfuD 1sx,y,dz/2dg2 + fuD 2sx,y,dz/2dg2d. s34d Subsequently, the stress is partitioned into its streamwise and cross-stream components in the usual manner: wall sx,yd = twsx,yd ti,3 FÎ uD isx,y,dz/2d uD 21 + uD 22 G , i = 1,2. s35d A sharp spectral cutoff filter is used in the wall stress and SGS computations. The fully scale dependent dynamic model increases the computational cost by about 20% compared to the Smagorinsky–Lilly model with imposed coefficient. Half of this increase is related to the dynamic computations of the coefficient while the other half is related to the Lagrangian averaging operations. These computational re- quirements are assessed when the coefficient cs,D is updated every fifth time step; this update frequency of the coefficient requires a time step for Lagrangian model Dt fEq. s17d and s25dg equal to five times the time step of the LES code dt. Tests for this study confirm that results obtained when updating the coefficient every single time step and every fifth time step are similar. B. Tests in homogeneous boundary layer flow Simulations over a homogeneous rough surface swith constant roughness zo imposed through the stress boundary condition at the lower surfaced were performed using: the simple Smagorinsky or Smagorinsky–Lilly model with a wall damping function sSMAGd, the PASI, the LASI, the PASD, and the LASD models. Table I details the parameters of the simulations. The cs,D coefficients computed by the different models sor imposed for the SMAG modeld for the homogeneous surface simulations are compared in Fig. 1. Close to the surface, all the dynamic formulations predict a lower value for cs,D than the one assumed by the Mason–Thompson damping function, whereas in the core of the flow they tend towards the classic value between 0.1 and 0.22. Regardless of the averaging method, the scale-dependent formulation yields larger coefficients than the scale-invariant formulation, consistent with the fact that near the surface the coefficient increases with decreasing scale. In addition, the planar averaged formulations predict lower cs,D than the equivalent Lagrangian formulation away from the wall and a higher cs,D close to the wall. We have verified that the larger Lagrangian averaged mean coefficients occur due to infrequent large values of cs,D obtained locally that dominate the mean values plotted but do not necessarily increase the mean SGS flux or mean dissipation to the same degree. Specifically, results presented later will show that the average dissipation character- Downloaded 24 Jan 2005 to 128.220.2.42. Redistribution subject to AIP license or copyright, see http://pof.aip.org/pof/copyright.jsp 025105-8 Bou-Zeid, Meneveau, and Parlange Phys. Fluids 17, 025105 ~2005! FIG. 1. Vertical profiles of the Smagorinsky coefficient kcs,Dlx,y,t for the different SGS eddy-viscosity models. istics for the different models are not influenced significantly by the averaging method but are much more sensitive to the scale-dependence of the model. Next, we present mean velocity profiles resulting from the simulations. For the homogeneous surface simulations, the velocity profile is expected to be logarithmic close to the surface sin the bottom 10%–20% of the simulation domaind, following kũl = su * / kdlnsz / zod. Figure 2 depicts the mean velocity fFig. 2sadg and the nondimensional velocity gradient skz / u * d]kũl / ]z fFig. 2sbdg obtained using the different SGS models. The solid black line in Fig. 2sad is the log law prediction with zo / H = 10−4. The Smagorinsky model results in a high gradient near the wall. This is in agreement with the previous findings22,32 suggesting that the model overdissipates resolved kinetic energy close to solid boundaries and thus the total Reynolds stresses are too low there leading to excessive mean velocity towards the core. Inversely, the two scale-invariant formulations sPASI and LASId produce insufficient dissipation leading to low velocity gradients and low streamwise velocity near the wall. On the other hand, the scale-dependent formulations yield a value of skz / u * d]kũl / ]z close to 1 near the surface suggesting that they are dissipating energy at a more appropriate rate and hence reproducing the log-law region more successfully than the other models. Note also that the present plane averaged scale-dependent results are slightly inferior to those presented previously.22 The discrepancy could be due to the fact that, in Ref. 22, a fifth-order polynomial was solved for b, instead of the more approximate method employed here based on Eq. s30d. Reproducing the log-law region depends on the ability of the SGS model to provide the correct dissipation rate close to the wall. However, a more complete insight into the energy dissipation characteristics of the SGS closure can be obtained by examining the streamwise velocity spectra. In the inertial subrange sk1z . 1, where k1 is the wavenumber and z is the distance to the walld, the effects of viscosity, boundary conditions, and large scale structures are not important and the turbulence is essentially isotropic. The energy cascade in this subrange follows the Kolmogorov spectrum yielding a FIG. 2. sad Normalized streamwise mean velocity profiles; the solid black line is the log law profile with zo / H = 10−4, sbd nondimensional mean velocity gradients for different SGS models. slope of −5 / 3. In the production range sk1z , 1d, the energy cascade is affected by the flow configuration. In wallbounded flows with neutral stability, there is evidence that the longitudinal sin the streamwise directiond energy spectrum of streamwise velocity in the production range follows a slope of −1.44–46 We remark that there exists evidence ssee, e.g., Refs. 47–49d that the k−1 regime does not extend over significant ranges of wavenumbers, and thus the matter is not conclusively settled from the data. Nevertheless, at the relatively short range of scales afforded by the resolution level of our simulations, the expectation of an approximate k−1 region still provides a useful criterion to test the various models. Figure 3 depicts the longitudinal u spectra produced by the different SGS models. The conclusions are similar to what was discussed for the log-law prediction. With the Smagorinsky model, too much energy is dissipated and the spectra decay much too fast at high wavenumbers. The Lagrangian scale-invariant formulation sLASId gives spectra that are too flat at small distances from the surface, indicating insufficient energy dissipation and a buildup of energy at the smallest resolved scales. The spectra of the PASI model are not shown here but depict underdissipation problems near the wall very similar to the LASI model results. The spectra for Downloaded 24 Jan 2005 to 128.220.2.42. Redistribution subject to AIP license or copyright, see http://pof.aip.org/pof/copyright.jsp 025105-9 A scale-dependent Lagrangian dynamic model Phys. Fluids 17, 025105 ~2005! FIG. 3. Normalized streamwise spectra of streamwise velocity vs k1z sthe lines represent the spectra obtained at heights z / H = 0.008, 0.024, 0.04, 0.056, 0.087, 0.119, 0.151, 0.183, 0.214, 0.246, 0.31, 0.373, 0.437, and 0.5d. the SMAG, LASI, and PASI models are most obviously in error close to the wall, where the assumption of isotropic homogeneous turbulence and the assumption of scale invariance clearly do not hold. The spectra obtained with the LASD and PASD models follow the −1 and −5 / 3 slopes well in the two ranges. This confirms that a dynamic scaledependent formulation is important for non wall-resolving LES in the vicinity of walls. Notice also that the LASD model is slightly less dissipative than the PASD model. The same is true when comparing the scale-invariant versions of the planar-averaged and Lagrangian-averaged model sthis holds for the results in this study and in Ref. 24d. C. Lagrangian scale dependent model detailed results In the preceding section, basic results of the LASD model were presented and compared with other SGS models. More comprehensive results obtained with the LASD model are presented in this section along with an analysis of the results sensitivity to the resolution of the simulation. The LES simulations over a homogeneous surface were run at resolutions of 323, 643, 963, and 1283. Figures 4sad and 4sbd depict the resolved, subgrid-scale, and total stress profiles for the 643 and 1283 resolutions. As expected, close to the surface most of the stress is in the SGS part while away from the wall most of the stress is resolved. The results are in qualitative agreement with many LES studies ssee Ref. 22 or the comparative study in Ref. 43, for exampled. Quantitative comparison is difficult since the partitioning of the stress into resolved and SGS components depends on the resolution and the SGS model used. The effect of resolution can be easily observed in Fig. 4 where the 1283 resolution results in higher resolved stresses sand lower SGS stressesd compared to the 643 simulation. The total stresses remain the same but the higher resolution simulation can resolve a larger part of these stresses. The variances of the resolved velocities for the different resolutions are plotted in Fig. 5 versus z / D. The results of the different resolutions match reasonably well. The profiles collapse at small z / D whereas further up an increase in resolution leads to an increase in the resolved variances. The variances of the current LES results fall well within the range reported by Andren et al.43 in their comparative study of four LES codes. In addition, the results coincide well with the variances reported by Porté-Agel et al.22 in their LES of ABL flow with a planar scale-dependent dynamic model that uses a different approach to “measure” scale dependence ssee Sec. II Bd. We confirm that the mean coefficient determined from LASD is, to a good approximation in this case of neutrally buoyant homogeneous boundary layer flow, a universal function of height z when normalized with the filter scale D. The vertical profiles of kcs,Dl computed at different resolutions are plotted versus z / D in Fig. 6. It is clear that, except for Downloaded 24 Jan 2005 to 128.220.2.42. Redistribution subject to AIP license or copyright, see http://pof.aip.org/pof/copyright.jsp 025105-10 Bou-Zeid, Meneveau, and Parlange Phys. Fluids 17, 025105 ~2005! FIG. 6. Collapse of the vertical profiles of the mean dynamic coefficient obtained from LASD model for different resolutions when plotted against z / D. FIG. 4. Vertical profiles of the resolved stress −kũ8w̃8lt,y,x s…d, the subgridscale stress −ktxzlt,y,x s1d, and the total stress ssolid line,—d: sad 643 resolution, sbd 1283 resolution. small deviations for the 323 simulation sdue to the effect of the top boundary condition in that cased, the results collapse well near the wall. Next, the variability of the coefficient computed from LASD is documented. Figure 7 depicts the probability den- sity functions spdfsd of the dynamic coefficient at three different heights; for each height, the pdfs obtained from the four different resolutions are shown. The distribution has low standard deviation and mean at z = dz / 2, the mean as well as the spread of the data increases further up. Consistent with the pdfs of Lagrangian coefficients determined using the scale-invariant dynamic model in wall-resolving LES,24 the peaks of the pdfs are close to the mean values and the pdfs do not exhibit secondary peaks or unusual features. Also note the spike in the pdf at zero associated with the use of a time scale that becomes infinite as cs tends to 0; this is in agreement with the scale-invariant version of the model.24 For the 323 resolution, considerable difference can be noted between the pdfs at z = H / 8 and z = H / 4. For the higher resolution runs, this difference decreases; this is a direct consequence of the decrease in the ratio of the grid scale D to the integral length scale ,z, at a given height, as the resolution is increased. Therefore, as the resolution increases, we note a decrease in the height at which the turbulence near the gridscale approaches isotropy. Figure 8 depicts the probability density function of the scale-dependence parameter b and the threshold of 1 / 8 FIG. 5. Normalized variances of resolved velocity components obtained with the LASD model at different resolutions. Downloaded 24 Jan 2005 to 128.220.2.42. Redistribution subject to AIP license or copyright, see http://pof.aip.org/pof/copyright.jsp 025105-11 A scale-dependent Lagrangian dynamic model Phys. Fluids 17, 025105 ~2005! FIG. 7. Probability density functions for the dynamic coefficient cs,D computed using LASD, at different heights and different resolutions: sad 323, sbd 643, scd 963, and sdd 1283. above which it is clipped. Peaks of the pdfs are well above the clipping limit except at z = dz / 2 where the peak is close to the clipping limit; however, even at dz / 2, most of the b values occur above the clipping limit. A small increase in the pdf can be noticed at the value of zero, caused by the use of a Lagrangian model time scale to avoid negative values of cs,4D, which yields a slightly increased probability of cs,4D = 0 fi.e., b = 0 according to Eq. s26dg. model sabout 10% more than the equivalent planar averaged versiond. Moreover, one should consider whether the variations actually affect relevant parameters such as fluid velocities and stresses. These two questions will be addressed in this section where simulations over heterogeneous walls are performed using the LASD and the PASD models and compared to experimental data. The other models sSMAG, LASI, and PASId, having already performed poorly in simulations over homogeneous surfaces, will not be tested. V. SIMULATIONS OF BOUNDARY LAYER FLOW OVER A HETEROGENEOUS ROUGH SURFACE A. Bradley’s experimental setup The previous results focused on flow over homogeneous rough surfaces and showed that the scale-dependent approach yields better results than the ad hoc Smagorinsky model or the scale-invariant formulations near the surface. Nevertheless, the motivation for using a Lagrangian local model rather than a plane averaged model becomes apparent only when examining more complex, nonhomogeneous flows si.e., when homogeneous directions are not available for averagingd. In the simulations over heterogeneous walls presented below, while averaging over horizontal planes is possible, the streamwise direction is not homogeneous and hence the variability of the coefficient in that direction would be suppressed by planar averaging. In such flows, an important question is whether spatial variations in the coefficient are great enough to justify the extra cost of the Lagrangian The atmospheric boundary layer measurements by Bradley50 have often been used to validate theoretical and numerical models for flow over an abrupt change in surface roughness.51–54 Bradley measured the surface stress swith drag platesd and the velocity profiles upstream and downstream of a sudden jump in surface roughness. The measurements were performed over a tarmac ssurface roughness zo = 0.002 cmd. A patch with a higher roughness szo = 0.25 cmd was created inside the tarmac by laying artificial roughness mats consisting of vertical spikes with reinforcing mesh in between. For the low-to-high roughness transition measurements, the high-roughness patch was created at the downstream end of the tarmac and measured 26 m 3 20 m in the streamwise and cross-stream directions, respectively. The tarmac had an effective area upstream of the measurement of Downloaded 24 Jan 2005 to 128.220.2.42. Redistribution subject to AIP license or copyright, see http://pof.aip.org/pof/copyright.jsp 025105-12 Bou-Zeid, Meneveau, and Parlange Phys. Fluids 17, 025105 ~2005! FIG. 8. Probability density functions of the scale-dependence parameter b obtained dynamically from the LASD approach, at different heights and different resolutions: sad 323, sbd 643, scd 963, and sdd 1283. 134 m 3 92 m in the streamwise and cross-stream directions, respectively. For the high-to-low roughness measurements, the high-roughness patch was slightly shorter, measuring 22 m in the streamwise direction, and the patch was placed near the upstream end of the tarmac. In the LES, in the cross-stream direction, a domain width of 64 meters is simulated with a resolution dy = 1 m. Furthermore, the high-roughness patch is assumed to extend over the entire cross-stream direction of the simulation si.e., infinitely wide patchd. For the results examined in this work, the cross-stream dimensions of the tarmac and the highroughness patch have very little impact on the results; this was confirmed in LES tests not presented in this work that included rectangular patches of various widths. The computational domain height is H = 20 meters, i.e., only the nearsurface layer of the ABL is simulated to allow for a high vertical resolution sdz = 10 cmd that can capture internal boundary layers originating at the transition between patches of different roughness. Figure 9 depicts the computational domain used to simulate conditions similar to Bradley’s field experiments. Two simulations were performed to reproduce Bradley’s low-to-high and high-to-low roughness transitions. Note that, due to the pseudospectral approach used in the code, what is actually being modeled is an infinite sequence of high-roughness and low-roughness starmacd patches. This does not affect the comparison results since we look at stress and velocity at the roughness jump, very close to the wall where the upstream conditions should not be very important. In practice, since the upstream conditions outside of the tarmac are unknown for the Bradley experiment, no better alternative to define inflow conditions exist. FIG. 9. LES parameters and simulation domain for reproducing Bradley’s field experimental study of atmospheric surface layer flow over abrupt roughness transitions. Downloaded 24 Jan 2005 to 128.220.2.42. Redistribution subject to AIP license or copyright, see http://pof.aip.org/pof/copyright.jsp 025105-13 Phys. Fluids 17, 025105 ~2005! A scale-dependent Lagrangian dynamic model C. Mean velocity profiles FIG. 10. Wall shear stress downstream of a low-to-high roughness jump, normalized by the equilibrium stress of the upstream surface: comparison of LES results and Bradley’s field data. B. Mean surface shear stresses Bradley plotted the mean surface shear stress measured downstream of the jump in roughness divided by the surface stress directly upstream of the jump. This ratio eventually relaxes downstream of the jump to the equilibrium value corresponding to the downstream surface. The ratio mainly depends on the two surface roughness values. The LES has to correctly predict this equilibrium ratio as well as the relaxation rate or the distance downstream of the jump where the ratio reaches its equilibrium value. Since LES codes with wall modeling “measure” the surface stress from the law-ofthe-wall using the velocity at the first grid point from the wall, they cannot accurately predict the departure from equilibrium of the stress ratio immediately after the roughness jump. This is due to the fact that the first grid points away from the wall do not lie in the internal equilibrium layer sIELd of the downstream patch and hence are affected by the upstream patch. Since the first grid point is 5 cm sdz / 2d above the wall, one expects that the values of the stress ratio will be affected by this source of error up to 5 m downstream of the jump in roughness sabout 5 cm3 100, see Brutsaert55 for scaling approximations for the IELd. Figures 10 and 11 depict the stress ratios measured by Bradley and simulated by the LES for low-to-high and high-to-low roughness transitions, respectively. The LASD model results agree well with the experimental results, especially past the downstream distance of 5 m after which the flow is in equilibrium with the underlying surface and the LES can accurately measure surface stress. The agreement between PASD model results and experimental data is slightly less satisfactory. The streamwise velocity profiles produced by the LES are compared to the profiles measured by Bradley at various downstream distances. Following Bradley’s approach, the velocities are normalized by the value at some height outside the internal boundary layer of the downstream patch s220 cm for the low-to-high roughness transition and 112.5 cm for the high-to-low roughness transitiond. This is necessary since Bradley’s data at different downstream distances are obtained at different times and hence cannot be normalized by u* su* will not be constant at the various times the data are acquiredd. Figure 12 depicts the comparison results. For the low-to-high roughness transition, the two SGS models and the experimental data agree very well. For the high-to-low roughness transition, the two SGS models produce different results and the LASD data agrees better with the experimental data. This has been traced to a greater difference in SGS dissipation between the two models over the low-roughness surface scompared to the difference over the high-roughness surfaced as will be shown in the following section. D. Coefficient variability over different patches and its effects In Fig. 13sad, the value of the Smagorinsky coefficient averaged in the cross-stream direction and in time, kcs,Dly,t, is plotted for the simulation with the LASD model reproducing Bradley’s data for the high-to-low roughness transition. The coefficient is divided by its average over horizontal planes, kcs,Dlx,y,t, to remove the effect of vertical variations and focus on horizontal variability. Up to 150% sminimum is ,0.5 of mean and maximum ,1.3 of meand variations between the high-roughness and low-roughness patches can be observed near the ground, and the effect of surface heterogeneity extends well into the lowest 10% of the domain. This sensitivity of cs,D to surface heterogeneity will not be captured by the planar-averaged model resulting in a difference in the SGS dissipation produced by the two models. This difference in estimating the SGS dissipation impacts the stress and velocity results produced by the models as depicted in Figs. 11 and 12. The dissipation of the LASD model divided by the dissipation of the PASD model is depicted in Fig. 13sbd. It is clear that the difference is more significant and extends higher above the low roughness surface scompared to the difference above the high-roughness surfaced. E. Sensitivity to grid resolution FIG. 11. Wall shear stress downstream of a high-to-low roughness jump normalized by the equilibrium stress of the upstream surface: comparison of LES results to Bradley’s field data. Similar to the simulations over homogeneous surfaces, different grid resolutions were tested for the heterogeneous surfaces simulations reproducing Bradley’s experiment. The results presented above pertain to the highest resolution used s1603 643 200 nodesd. Two other resolutions were tested: a low resolution of 803 323 100 nodes and a medium resolution of 1203 483 150 nodes. The results showed a consistent improvement in the reproduction of Bradley’s experimental data as the grid resolution was increased; this applies for both the PASD and LASD subgrid scale models. In this section, we only present the results obtained with the LASD Downloaded 24 Jan 2005 to 128.220.2.42. Redistribution subject to AIP license or copyright, see http://pof.aip.org/pof/copyright.jsp 025105-14 Phys. Fluids 17, 025105 ~2005! Bou-Zeid, Meneveau, and Parlange FIG. 12. Evolution of velocity profiles after an abrupt change in surface roughness: comparison of LES results to Bradley’s field data. model. Figures 14 and 15 depict the wall shear stress downstream of a low-to-high and high-to-low roughness jumps, respectively. One observes that as the resolution decreases and the first grid point moves further away from the wall, the stress profile becomes flatter. Figure 16 shows the velocity profile adjustment at about four meters downstream of the roughness jump. For the lowto-high transition, all resolutions are able to capture the experimental velocity profile rather well; this is in agreement with the results of the sensitivity to the SGS model. For this type of transition, the LES seems to be able to capture the velocity profile regardless of the numerical details and the SGS model. On the other hand, for the high-to-low roughness transition, the results from the low resolution simulation do not match the experimental data well. The improvement for the medium resolution is significant and little further improvement is obtained by passing to high resolution; this suggests that the results are converging to the experimental data as the resolution is increased. This type of roughness jump shigh-to-lowd is much more sensitive to the numerical details and the SGS model used. VI. CONCLUSIONS A scale-dependent dynamic subgrid scale model has been formulated based on the Lagrangian time-averaging approach sLASDd. The model extrapolates the Smagorinsky coefficient measured dynamically at two test-filter scales towards the grid-scale where the coefficient is unknown and needed for LES. Results show that the model performs well Downloaded 24 Jan 2005 to 128.220.2.42. Redistribution subject to AIP license or copyright, see http://pof.aip.org/pof/copyright.jsp 025105-15 A scale-dependent Lagrangian dynamic model Phys. Fluids 17, 025105 ~2005! FIG. 13. sColord. Sensitivity of SGS model to surface roughness for a high-to-low roughness transition: sad plot of the dynamic coefficient for the LASD model, sbd SGS dissipation predicted by the LASD model sPLASDd divided by the SGS dissipation predicted by the PASD model sPPASDd. in reproducing the log-law region in high-Reynolds number boundary layers that do not resolve the viscous sublayer, where scale-dependence is important. The model performs better than other SGS models tested here, including the Smagorinsky–Lilly model with a prescribed wall damping function, and the scale-invariant dynamic model. Streamwise velocity spectra indicate that the LASD model yields more accurate mean SGS dissipation properties. The model produces only a moderate increase in computational cost on the order of 20% scompared to the Smagorinsky model with imposed coefficientd, since the dynamic coefficient need not be updated at every time step of the LES shere we updated it at every fifth time step, but this choice depends on Courant– Friedricks–Lévy stability constraints of particular simulation parametersd. To test the model in an inhomogeneous flow under controlled conditions, it has been applied to LES to compare with the experimental field results of Bradley50 for highReynolds number boundary layer flow over a roughness discontinuity. The experimental results consist of stress and ve- FIG. 14. Wall shear stress downstream of a low-to-high roughness jump, normalized by the equilibrium stress of the upstream surface: sensitivity to grid resolution for the LASD SGS model. HR is for high resolution s160 3 643 200d, MR is for medium resolution s1203 483 150d, and LR is for low resolution s803 323 100d. FIG. 15. Wall shear stress downstream of a high-to-low roughness jump normalized by the equilibrium stress of the upstream surface: sensitivity to grid resolution for the LASD SGS model. HR is for high resolution s160 3 643 200d, MR is for medium resolution s1203 483 150d, and LR is for low resolution s803 323 100d. Downloaded 24 Jan 2005 to 128.220.2.42. Redistribution subject to AIP license or copyright, see http://pof.aip.org/pof/copyright.jsp 025105-16 Phys. Fluids 17, 025105 ~2005! Bou-Zeid, Meneveau, and Parlange FIG. 17. Variance of uD 1 and ratio of kuD 1l2 / kũ¯21l as a function of filter size. APPENDIX: A LOCAL LAW-OF-THE-WALL FORMULATION In this appendix, an improved local law-of-the-wall is introduced. The common approach used in LES of highReynolds number boundary layer flow where the viscous sublayer is not resolved2,8,56,57 is to impose the law-of-thewall in a strictly local sense: tw = − FIG. 16. Evolution of velocity profiles after an abrupt change in surface roughness: sensitivity to grid resolution for the LASD SGS model. HR is for high resolution s1603 643 200d, MR is for medium resolution s1203 48 3 150d, and LR is for low resolution s803 323 100d. locity measurements downstream of the roughness jump. The LES was successful in predicting wall stress adjustment as a function of downstream distance. Similarly, the velocity profiles from LES data coincided well with the experimental data at several downstream distances included in the analysis. Results obtained with the Lagrangian model matched the experimental data better than results obtained with the equivalent planar-averaged model. This difference was traced back to the sensitivity of the Smagorinsky coefficient to the roughness height of the underlying surface. This sensitivity affects the SGS dissipation rate and cannot be captured by the planar-averaged formulation. This indicates that a local formulation, such as the Lagrangian one, is better suited for simulations of flows in complex geometries. ACKNOWLEDGMENTS This study was funded by NSF under Grant No. EAR9909679 and by Agreement No. R828771-0-01 from the U.S. Environmental Protection Agency’s Science to Achieve Results sSTARd program. The authors are thankful to the Scientific Computing Division of the National Center for Atmospheric Research sNCARd for the use of their computer clusters. F k lnsz/zod G 2 sũ21 + ũ22d. sA1d Here k is the von-Karman constant s<0.4d and z = dz / 2. tw is the kinematic stress t / r sthe squared friction velocityd at the wall. The use of this relation imposes an average stress obtained from LES: ktwLESl = − F k lnsz/zod G 2 skũ21l + kũ22ld. sA2d However, the log law was developed and validated to be used in an average sense, i.e., ktwlogl = − F k lnsz/zod G 2 kũ1l2 , sA3d where the mean cross-stream component kũ2l is zero. Since the velocity at z = dz / 2 fluctuates, kũ21l . kũ1l2 sSchwartz inequalityd and ktwLESl . ktwlogl. Therefore, imposing the wall stress in a local formulation leads to increased average stresses for a given near-wall velocity. In LES with prescribed pressure gradient and mean stress, this would yield a slower flow near the surface. A potential solution is to divide the stress into a mean contribution and a local contribution similar to the approach used to impose a velocity gradient at the surface. The local contribution should average to zero, yielding a proper estimate of the average stress. This formulation could be used in homogeneous terrain to impose a local stress. However, when complex or heterogeneous areas are to be simulated, defining the mean and the variation parts is not always feasible. Fortunately, it turns out that filtering the velocity sat z = dz / 2 and only to prescribe the wall stressd at a scale 2D already significantly reduces the small-scale fluctuations so Downloaded 24 Jan 2005 to 128.220.2.42. Redistribution subject to AIP license or copyright, see http://pof.aip.org/pof/copyright.jsp 025105-17 Phys. Fluids 17, 025105 ~2005! A scale-dependent Lagrangian dynamic model 7 FIG. 18. Nondimensional velocity gradient for different wall model formulations. that the velocity variance becomes quite small. A formulation that does not require averaging can be derived using this filtered local velocity to impose a stress as in Eqs. s34d and s35d of the main text. The filtering preserves the large scale variations slarger than 2Dd which are important. At the same time, use of the filtered velocities results in an average stress that is very close to the stress predicted by the average similarity formulation for homogeneous surfaces. Figure 17 depicts the variance of the streamwise velocity ũ and the ratio kuD 1l2 / kũ¯21l for the unfiltered data sdata at a scale Dd and for increasing filter size sobtained from our LASD simulationsd. A filtering at a scale of 2D is sufficient to increase uD 2 / ũ¯21 to about 0.985. Note that there is no guarantee that this will yield a correct stress for heterogeneous surfaces. While a universal formulation to compute the average local or average stress over heterogeneous surfaces is yet to be found, this formulation can be implemented for heterogeneous surfaces and will certainly give better results that the fully local formulation. Figure 18 presents the nondimensional velocity gradients for the different formulations, from simulations over homogeneous surfaces. The local formulation gives a very high gradient near the wall sand strong oscillationsd, this is related to the high stress and low velocity, predicted by this formulation, at the wall. On the other hand, the filtered formulation gives better results, as does the mean formulation which is known to be well suited for homogeneous cases. 1 J. W. Deardorff, “Three-dimensional numerical study of the height and mean structure of a heated planetary boundary layer,” Boundary-Layer Meteorol. 7, 81 s1974d. 2 C.-H. Moeng, “A large-eddy simulation for the study of planetary boundary layer turbulence,” J. Atmos. Sci. 41, 2052 s1984d. 3 C.-H. Moeng and J. C. Wyngaard, “Statistics of conservative scalars in the convective boundary layer,” J. Atmos. Sci. 41, 3161 s1984d. 4 R. H. Shaw and U. Schumann, “Large-eddy simulation of turbulent flow above and within a forest,” Boundary-Layer Meteorol. 61, 47 s1992d. 5 P. J. Mason, “Large-eddy simulation: A critical review of the technique,” Q. J. R. Meteorol. Soc. 120, 1 s1994d. 6 B. Kosović, “Subgrid-scale modeling for the large-eddy simulation of high-Reynolds-number boundary layers,” J. Fluid Mech. 336, 151 s1997d. J. D. Albertson and M. B. Parlange, “Natural integration of scalar fluxes from complex terrain,” Adv. Water Resour. 23, 239 s1999d. 8 J. D. Albertson and M. B. Parlange, “Surface length-scales and shear stress: implications for land-atmosphere interaction over complex terrain,” Water Resour. Res. 35, 2121 s1999d. 9 M. Lesieur and O. Métais, “New trends in large eddy simulations of turbulence,” Annu. Rev. Fluid Mech. 28, 45 s1996d. 10 U. Piomelli, “Large-eddy simulation: achievements and challenges,” Prog. Aerosp. Sci. 35, 335 s1999d. 11 P. Sagaut, Large Eddy Simulation for Incompressible Flows, 2nd ed. sSpringer, Berlin, 2003d. 12 J. S. Smagorinsky, “General circulation experiments with the primitive equations: I. The basic experiment,” Mon. Weather Rev. 91, 99 s1963d. 13 J. Bardina, J. H. Ferziger, and W. C. Reynolds, “Improved subgrid-scale models for large-eddy simulation,” AIAA Pap. 80, 1357 s1980d. 14 B. Vreman, B. Geurts, and H. Kuerten, “On the formulation of the dynamic mixed subgrid-scale model,” Phys. Fluids 6, 4057 s1994d. 15 F. Sarghini, U. Piomelli, and E. Balaras, “Scale-similar models for largeeddy simulations,” Phys. Fluids 11, 1596 s1999d. 16 R. Anderson and C. Meneveau, “Effects of the similarity model in finitedifference LES of isotropic turbulence using a Lagrangian dynamic mixed model,” Flow, Turbul. Combust. 62, 201 s1999d. 17 C. Higgins, M. Parlange, and C. Meneveau, “Alignment trends of velocity gradients and subgrid-scale fluxes in the turbulent atmospheric boundary layer,” Boundary-Layer Meteorol. 109, 59 s2003d. 18 H. S. Kang, S. Chester, and C. Meneveau, “Decaying turbulence in an active-grid-generated flow and comparisons with large-eddy simulation,” J. Fluid Mech. 480, 129 s2003d. 19 M. Germano, U. Piomelli, P. Moin, and W. Cabot, “A dynamic subgridscale eddy viscosity model,” Phys. Fluids A 3, 1760 s1991d. 20 M. Germano, “Turbulence: the filtering approach,” J. Fluid Mech. 238, 325 s1992d. 21 C. Meneveau and J. Katz, “Scale-invariance and turbulence models for large-eddy simulation,” Annu. Rev. Fluid Mech. 32, 1 s2000d. 22 F. Porté-Agel, C. Meneveau, and M. B. Parlange, “A scale-dependent dynamic model for large-eddy simulation: application to a neutral atmospheric boundary layer,” J. Fluid Mech. 415, 261 s2000d. 23 J. Kleissl, C. Meneveau, and M. B. Parlange, “On the magnitude and variability of subgrid-scale eddy-diffusion coefficients in the atmospheric surface layer,” J. Atmos. Sci. 60, 2372 s2003d. 24 C. Meneveau, T. Lund, and W. Cabot, “A Lagrangian dynamic subgridscale model of turbulence,” J. Fluid Mech. 319, 353 s1996d. 25 S. Ghosal, T. S. Lund, P. Moin, and K. Akselvoll, “A dynamic localization model for large eddy simulation of turbulent flows,” J. Fluid Mech. 86, 229 s1995d. 26 S. B. Pope, Turbulent Flows sCambridge University Press, Cambridge, UK, 2000d. 27 D. C. Haworth and K. Jansen, “Large-eddy simulation on unstructured deforming meshes: towards reciprocating IC engines,” Comput. Fluids 29, 493 s2000d. 28 V. Armenio and U. Piomelli, “A Lagrangian mixed subgrid-scale model in generalized coordinates,” Flow, Turbul. Combust. 65, 51 s2000d. 29 L. Blin, A. Hadjadj, and L. Vervisch, “Large eddy simulation of turbulent flows in reversing systems,” J. Turbul. 4, 1 s2003d. 30 M. Tsubokura, T. Kobayashi, N. Taniguchi, and W. P. Jones, “A numerical study on the eddy structures of impinging jets excited at the inlet,” Int. J. Heat Fluid Flow 24, 500 s2003d. 31 D. K. Lilly, “The representation of small scale turbulence in numerical simulation experiments,” Proceedings of the IBM Scientific Computing Symposium on Environmental Sciences sIBM, Yorktown Heights, NY, 1967d, p. 195. 32 P. J. Mason and D. J. Thomson, “Stochastic backscatter in large-eddy simulations of boundary layers,” J. Fluid Mech. 242, 51 s1992d. 33 D. K. Lilly, “A proposed modification of the Germano subgrid-scale closure method,” Phys. Fluids A 4, 633 s1992d. 34 S. B. Pope, “Ten questions concerning the large-eddy simulation of turbulent flows,” New J. Phys. 6, 35 s2004d. 35 C. Meneveau and T. Lund, “The dynamic Smagorinsky model and scaledependent coefficients in the viscous range of turbulence,” Phys. Fluids 9, 3932 s1997d. 36 E. Bou-Zeid, C. Meneveau, and M. B. Parlange, “Large-eddy simulation of neutral atmospheric boundary layer flow over heterogeneous surfaces: Blending height and effective surface roughness,” Water Resour. Res. 40, W02505 s2004d. Downloaded 24 Jan 2005 to 128.220.2.42. Redistribution subject to AIP license or copyright, see http://pof.aip.org/pof/copyright.jsp 025105-18 37 Phys. Fluids 17, 025105 ~2005! Bou-Zeid, Meneveau, and Parlange U. Piomelli and J. Liu, “Large-eddy simulation of rotating channel flows using a localized dynamic model,” Phys. Fluids 7, 839 s1995d. 38 C. Meneveau and T. Lund, “On the Lagrangian nature of the turbulence energy cascade,” Phys. Fluids 6, 2820 s1994d. 39 S. A. Orszag and Y. H. Pao, “Numerical computation of turbulent shear flows,” Adv. Geophys. 18, 224 s1974d. 40 S. A. Orszag, “Transform method for calculation of vector coupled sums: application to the spectral form of the vorticity equation,” J. Atmos. Sci. 27, 890 s1970d. 41 J. W. Deardorff, “A numerical study of three-dimensional turbulent channel flow at large Reynolds numbers,” J. Fluid Mech. 41, 453 s1970d. 42 U. Piomelli and E. Balaras, “Wall-layer models for large-eddy simulation,” Annu. Rev. Fluid Mech. 34, 349 s2002d. 43 A. Andren, A. R. Brown, J. Graf, P. J. Mason, C.-H. Moeng, F. T. M. Nieuwstadt, and U. Schumann, “Large-eddy simulation of a neutrally stratified boundary layer: a comparison of four computer codes,” Q. J. R. Meteorol. Soc. 120, 1457 s1994d. 44 B. A. Kader and A. M. Yaglom, “Spectra and correlation functions of surface layer atmospheric turbulence in unstable thermal stratification,” in Turbulence and Coherent Structures, edited by O. Métais and M. Lesieur sKluwer Academic, Dordrecht, 1991d pp. 387–412. 45 A. E. Perry, K. L. Lim, and S. M. Henbest, “An experimental study of the turbulence structure in smooth- and rough-wall boundary layers,” J. Fluid Mech. 177, 437 s1987d. 46 G. G. Katul, C. R. Chu, M. B. Parlange, J. D. Albertson, and T. A. Ortenburger, “The low-wavenumber spectral characteristics of velocity and temperature in the atmospheric surface layer,” J. Geophys. Res. 100, 14243 s1995d. 47 J. F. Morrison, W. Jiang, B. J. McKeon, and A. J. Smits, “Reynolds number dependence of streamwise velocity spectra in turbulent pipe flow,” Phys. Rev. Lett. 88, 214501 s2002d. 48 T. B. Nickels and I. Marusic, “On the different contributions of coherent structures to the spectra of a turbulent round jet and a turbulent boundary layer,” J. Fluid Mech. 448, 367 s2001d. 49 J. C. Del Alamo, J. Jimenez, P. Zandonade, and M. D. Moser, “Scaling of the energy spectra of turbulent channels,” J. Fluid Mech. 500, 135 s2004d. 50 E. F. Bradley, “A micrometeorological study of velocity profiles and surface drag in the region modified by a change in surface roughness,” Q. J. R. Meteorol. Soc. 94, 361 s1968d. 51 C. C. Shir, “A numerical computation of air flow over a sudden change of surface roughness,” J. Atmos. Sci. 29, 304 s1972d. 52 K. S. Rao, J. C. Wyngaard, and O. R. Coté, “The structure of the twodimensional internal boundary layer over a sudden change of surface roughness,” J. Atmos. Sci. 31, 738 s1974d. 53 W. H. Schofield, “Turbulent shear flows over a step-change in surface roughness,” J. Fluids Eng. 103, 344 s1981d. 54 H. P. Schmid and B. Bünzli, “The influence of surface texture on the effective roughness length,” Q. J. R. Meteorol. Soc. 121, 1 s1995d. 55 W. Brutsaert, “Land-surface water vapor and sensible heat flux: Spatial variability, homogeneity, and measurement scale,” Water Resour. Res. 34, 2433 s1998d. 56 P. J. Mason and N. S. Callen, “On the magnitude of subgrid-scale eddy coefficient in large-eddy simulations of turbulent channel flow,” J. Fluid Mech. 162, 439 s1992d. 57 H. Schmidt and U. Schumann, “Coherent structures of the convective boundary layer derived from large-eddy simulations,” J. Fluid Mech. 200, 511 s1989d. Downloaded 24 Jan 2005 to 128.220.2.42. Redistribution subject to AIP license or copyright, see http://pof.aip.org/pof/copyright.jsp