Incompressibility of the Leray-? model for wall-bounded flows

Incompressibility of the Leray-

␣

model for wall-bounded flows

M. van Reeuwijk,a兲H. J. J. Jonker, and K. Hanjalić

Department of Multi-Scale Physics and J.M. Burgers Center for Fluid Dynamics, Delft University of Technology, Lorentzweg 1, 2628 CJ Delft, The Netherlands

共Received 12 October 2005; accepted 8 December 2005; published online 24 January 2006兲 This study shows that the Leray-␣model does not explicitly enforce a divergence-free field for the filtered velocity. While this condition is automatically satisfied in the absence of boundaries, bounded domains require extra attention. It is shown, both analytically and through simulations of Rayleigh–Bénard convection, that incompressibility of the filtered velocity field cannot be guaranteed in the current formulation. Several suggestions are made to restore the incompressibility of the filtered velocity, and it is shown that free-slip boundary conditions for the filtered velocity do guarantee incompressibility for the domain under consideration. © 2006 American Institute of

Physics. 关DOI:10.1063/1.2166459兴

The Leray-␣model1is a promising method for simulat-ing three-dimensional turbulent flows. The model is inspired by the Lagrangian averaged Navier-Stokes-␣ 共LANS-␣兲 model,2–8 which can be derived using variational principles from a Lagrangian that has been averaged along fluid particle trajectories3,6,7,9 and adding a suitable viscous term. The LANS-␣ model has been applied to isotropic homogeneous turbulence and channel flow.4,10 So far the Leray-␣ model has been applied to isotropic homogeneous turbulence.1The governing equations are given by

⳵tui+ u˜j⳵jui=␯⳵j 2 ui−⳵ip + fi, 共1兲 ⳵juj= 0, 共2兲 u ˜i−␣2⳵2ju˜i= ui. 共3兲

Here implicit summation is performed over repeated indices and ⳵i=⳵/⳵xi. The rationale is to introduce a second,

smoother velocity field u˜i, obtained by applying a smoothing

filter 共3兲 with a filter size ␣ to ui, that advects the fluid,

thereby reducing the nonlinearity of the Navier–Stokes equa-tions. This principle has recently been proposed as a regular-ization model for large-eddy simulation11that allows a sys-tematic derivation of the implied subgrid model. The explicit presence of both unfiltered and filtered velocity in 共1兲 sug-gests similarity with the approximate deconvolution method12 and the variational multiscale method 共see, e.g., Ref. 13兲. The distinguishing aspect of the regularization ap-proach is that it modifies the distribution of energy in spec-tral space, see, e.g., Ref. 5. In this respect it differs signifi-cantly from other subgrid closures, as these normally employ eddy-viscosity concepts, thus modifying the dissipative pro-cesses.

An attractive feature of the Leray-␣ model is that it can be easily implemented in existing computational fluid dy-namics共CFD兲 codes, as the filtered velocity u˜ican be treated

as a dependent variable. The main reason for this is that the continuity equation共2兲 is in terms of ui, instead of u˜ias is the

case for the LANS-␣model. However, as can be seen from 共3兲, ⳵j˜uj= 0 automatically implies that ⳵juj= 0, but not vice

versa. Consequently, extra attention is required to ensure that ⳵j˜uj= 0 for the Leray-␣model.

One may question the importance of u˜i not being

divergence-free. Indeed, no “real” mass is lost in terms of the

uifield, as u˜imerely represents an average velocity by which

the velocity field ui is advected. However,⳵j˜uj⫽0 has

seri-ous implications, for example, for the budget of kinetic en-ergy e =1₂uiui, as the important identity ui˜uj⳵jui=⳵j˜uje fails

to hold because the term e⳵j˜ujis no longer zero. Because of

this, the advection is no longer a pure redistributor of energy, but becomes an active source/sink of kinetic energy as well. This is an artifact that changes the dynamics of the system fundamentally and uncontrollably, and thus⳵j˜uj⫽0 is

physi-cally unacceptable. As shown below, naively imposing no-slip boundary conditions for both uiand u˜iresults in

signifi-cant compressibility effects in the near-wall region, especially for higher values of ␣. For many engineering problems this region is of crucial importance and violation of the incompressibility of u˜iwould render the Leray-␣ model

useless for these cases.

Our aim in this Brief Communication is to show that the Leray-␣model in its present formulation cannot guarantee a divergence-free u˜ifield, even for simple domains. First, this

is illustrated by applying the model to simulate Rayleigh– Bénard convection. Then it is shown analytically, for a do-main bounded by impermeable walls at z = −1₂H and z =1₂H,

how enforcing no-slip boundary conditions for both the fil-tered and unfilfil-tered velocity results in serious compressibil-ity effects, specifically in the near-wall region. In the last part several strategies are discussed to resolve this issue, and it is shown that at least for this simple geometry,⳵ju˜j= 0 can be

ensured by applying free-slip boundary conditions for u˜i.

To illustrate the importance of incompressibility for u˜i,

numerical simulations of Rayleigh–Bénard convection are carried out. Rayleigh–Bénard convection is generated when a layer of fluid enclosed by two flat plates perpendicular to the gravitational acceleration is subjected to a positive tem-perature difference between the top and bottom plate; see, a兲_{Electronic mail: m.van.reeuwijk@ws.tn.tudelft.nl}

PHYSICS OF FLUIDS 18, 018103共2006兲

1070-6631/2006/18共1兲/018103/4/$23.00 18, 018103-1 © 2006 American Institute of Physics

e.g., Ref. 14. The Helmholtz equation共3兲 has been added to our code for direct numerical simulation共DNS兲15 and a di-rect solver共taking advantage of the homogeneous directions兲 is used to invert it. Simulations are performed at Ra= 105and Pr= 1 for a⌫=L/H=4 aspect ratio domain with H=0.15 m, ⌬⌰=2 K, ␤= 1.74⫻10−4_K−1_{, and g = 9.81 ms}−2_{. Periodic} boundary conditions are applied for the sidewalls, and no-slip boundary conditions for ui with a fixed temperature on

the bottom

共

z = −1₂H

兲

and top

共

z =1₂H

兲

walls. A mesh of 64⫻64⫻64 cells is used for all simulations, which is of sufficient resolution for DNS. Thus, the reduction of active scales by the Leray-␣model is not exploited here in order to minimize numerical errors. The typical convective turnover time␶= 36 s and the nondimensional heat-flux, i.e., the Nus-selt number, is Nu= 4.7. The filter size is set to␣/ H = 0.1.

In the simulations, the unfiltered uifield is

divergence-free up to machine precision. Applying the standard decom-position into mean and fluctuating parts, ui= u¯i+ ui

⬘

, this

means that the mean divergence⳵j¯ujand the standard

devia-tion⳵ju

⬘

j⳵kuk

⬘

1/2 are zero too. The same should be expected

from the mean divergence and the standard deviation of the filtered velocity, u˜. However, when applying no-slip 共Nsl兲

boundary conditions for u˜i, significant compressibility effects

can be observed 关Fig. 1共a兲兴. On average, the u˜i field is

divergence-free, as␶ ⳵j˜uj= O共10−30兲, but there are strong

de-viations specifically near the wall, as can be concluded from the profile of the standard deviation ␶ ⳵j˜uj

⬘

⳵k˜uk

⬘

1/2, which is

O共1兲. For comparison, the normalized mean and standard

deviation of the divergence of the uifield are O共10−30兲.

The real problem is the loss of the purely redistributive character of advection. Shown in Fig. 1共b兲 is the transport of turbulent kinetic energy 共TKE兲 by fluctuations in the near-wall region. This term comes from the advection of fluctua-tions ui

⬘

by the fluctuations u˜i, multiplying it by uiand

aver-aging,

ui

⬘

⳵j共u˜

⬘

jui

⬘

兲 =⳵j˜u

⬘

je

⬘

+ e

⬘

⳵j˜u

⬘

j, 共4兲

with e

⬘

=1₂u_i

⬘

u_i

⬘

. Note that the advective term on the left-hand side is written in divergence form, in accordance with its implementation in the code. The first term on the right-hand side is in divergence form so is purely redistributive. The second term is a source/sink, which normally vanishes be-cause the fluctuating field is divergence-free. The results are normalized by the volume-averaged dissipation rate 具⑀典V=␯3H−4Ra共Nu−1兲Pr−2, which is one of the important

exact results for Rayleigh–Bénard convection.14 The

z-coordinate is normalized by the thermal boundary layer

thickness, which can be estimated by ␭_⌰= H /共2 Nu兲. With the no-slip共Nsl兲 boundary conditions for u˜i, it can be seen

关Fig. 1共b兲兴 that there is a significant production of turbulent kinetic energy due to the term −e

⬘

⳵j˜uj

⬘

, which is up to 25% of

the average production of TKE well inside the thermal boundary layer. This artificial energy injection is an unac-ceptable side effect, and below follows an analysis of the cause of this behavior.

In the absence of boundaries it is straightforward that the incompressibility of ui is a sufficient condition to ensure

⳵j˜uj= 0. Applying a three-dimensional Fourier transform to

共3兲 gives that u ˜ˆi= 1 1 +␣2k2j uˆi, 共5兲

which shows that the effect of the filter is roughly to damp wave numbers k⬎␣−1. Taking the divergence of 共3兲, intro-ducing g =⳵juj and f =⳵j˜uj, and applying the Fourier

trans-form results in

fˆ = 1

1 +␣2kj

2gˆ, 共6兲

indicating that the divergences are algebraically dependent. From this it can be concluded that in the absence of bound-aries, condition共2兲 is sufficient to ensure the incompressibil-ity of u˜i.

Care has to be taken when using the Leray-␣framework in 共partially兲 bounded domains. A domain bounded by im-permeable walls at z = −1₂H and z =1₂H will be considered.

After performing a Fourier transform in the two unbounded horizontal directions x and y,共3兲 turns into

␣2_共K2_{− d} z 2_兲u˜ˆ i= uˆi, 共7兲 with K defined as K =

冑

␣−2+ kx 2 + ky 2 . 共8兲

Taking the divergence in共7兲 and using 共2兲 results in

FIG. 1. Effect of boundary conditions for u˜ifor a simulation of Rayleigh–

Bénard convection at Ra= 105_{, Pr= 1, and ␣/ H = 0.1. Nsl and Fsl denote}

results with no- and free-slip boundary conditions for u˜i, respectively.共a兲

Mean square divergence of the fluctuating u˜ifield;共b兲 transport of TKE by

turbulent fluctuations.

018103-2 van Reeuwijk, Jonker, and Hanjalić Phys. Fluids 18, 018103共2006兲

␣2_共K2_{− d}

z

2_{兲fˆ = 0, 共9兲}

where fˆ = ıkx˜ˆ + ıku y˜ˆ + dv zw˜ˆ . From 共9兲 it is clear that

incom-pressibility of u˜i is not guaranteed, as the divergence fˆ is

governed by an ordinary differential equation in z, which allows a homogeneous solution of the form

fˆ = AeKz+ Be−Kz. 共10兲

This elucidates the importance of boundary conditions, as

fˆ共kx, ky, z兲 is directly determined by A and B. Assuming that

identical boundary conditions are applied on both sides, it is easily shown that fˆ共kx, ky, z兲=0 can be obtained by setting fˆ

or any z derivative of fˆ to zero at the boundaries. Conse-quently, the velocity field u˜i is divergence-free when the

boundary conditions, either directly or indirectly, ensure that

dz n fˆ = ıkxdz n u ˜ˆ + ıkydz n v ˜ˆ + d_zn+1_{w˜ˆ = 0 共11兲}

at z = ± H / 2 for at least one n苸N.

When no-slip boundary conditions uˆ =vˆ = wˆ = 0 and u

˜ˆ = v˜ˆ = w˜ˆ = 0 are applied for ui and u˜i, respectively, a

divergence-free field for u˜i is obtained only under special

circumstances. At the wall, using共2兲 and 共7兲, the following terms are zero: dzwˆ = dz

2 u

˜ˆ = d_z2v˜ˆ = d_z2w˜ˆ = 0. Using these

rela-tions, the first three terms of共11兲 are given by

fˆ = dzw˜ˆ ; dzfˆ = ıkxdz˜ˆ + ıku ydz˜ˆ ;v dz

2 fˆ = dz

3 w˜ˆ .

No terms vanish based on the boundary conditions a priori for n⬎2. Thus the velocity field is only conditionally divergence-free; only under special circumstances, e.g., when dzw˜ˆ = 0 will the field u˜ibe divergence-free. This occurs,

for example, in laminar flows where w = 0, such as Pois-seuille flow. However, the Leray-␣ model is intended as a subgrid model for high Re turbulent flows, so, in general, the filtered velocity field will not be divergence-free under these circumstances.

The filtered velocity field is not automatically divergence-free due to the elliptic nature of共3兲, which im-plies nonlocal interactions, as the filtered velocity at the wall is influenced by the entire unfiltered velocity field. This can be made explicit by the solution of 共7兲 in terms of the Green’s functions: u ˜ˆi共z兲 =

冕

−H/2 H/2 G共z,␨兲uˆi共␨兲d␨,

where the dependence of uˆ and u˜ˆ on the wave numbers kx

and kyhave been omitted for convenience. The Green’s

func-tion for w˜ˆ , applying the Dirichlet boundary conditions w˜ˆ = 0

at z = ± H / 2, is given by GD共z,␨兲 =

冦

sinh关K共H/2 −␨兲兴sinh关K共H/2 + z兲兴 K␣2sinh共KH兲 , z艋␨, sinh关K共H/2 +␨兲兴sinh关K共H/2 − z兲兴 K␣2sinh共KH兲 , z艌␨.

冧

Differentiating with respect to z and integrating by parts, the derivative dzw˜ˆ at the bottom wall is related to the unfiltered

derivatives dzwˆ共z兲 by

dzw˜ˆ共− H/2兲 =

冕

−H/2

H/2 _{cosh关K共␨− H/2兲兴}

K␣2sinh共KH兲 d␨wˆ共␨兲d␨.

This relation shows that dzw˜ˆ 共and thus the divergence兲 at the

wall depends on the entire field of dzwˆ , and it is clear that the

constraint dzwˆ = 0 at the wall关following from substitution of

the boundary conditions into共2兲兴 is not sufficient to ensure that dzw˜ˆ = 0. When␣ⰆH, the mutual influence of the walls is

negligible and KHⰇ1. In this case, the integration kernel approaches K−1_␣−2_{exp共−Kd兲 with d=␨+ H / 2 as the distance} from the wall. Hence the typical region that influences dzw˜ˆ is

proportional to K−1_{⬇␣共for␣sufficiently small兲. In the limit} of ␣→0, the integration kernel converges to a Dirac delta function and the field of u˜ˆi will be divergence-free as

dzw˜ˆ→dzwˆ = 0 at the wall. However, the added value of the

Leray-␣ model is for nonzero values of ␣, and it can be concluded that the configuration of no-slip boundary condi-tions for uiand u˜ido not suffice for turbulent flows.

Several strategies could be used to prevent the undesir-able behavior of the divergence, u˜i. First, one could opt to

enforce⳵ju˜j= 0 instead of共2兲, which would avoid the

diver-gence issues a priori. Doing this, ui can be eliminated by

substituting共3兲 in 共1兲, which after some calculation can be written in a LES-type formulation, but with an asymmetric residual stress tensor mij given by mij= u˜ u
−uj i ˜ ui˜ .j11 This is

the most generic and intuitive approach, but from a compu-tational perspective it seems rather unfavorable, as it in-volves solving the Helmholtz equation for all components of

mij.

Second, ␣ could be made to vanish when approaching the wall. By doing this, u˜i would be automatically

divergence-free, as dzw˜ˆ→dzwˆ when␣→0 at the wall. This

would mean that no-slip boundary conditions could be used for u˜i. This is in fact what is proposed for the more general

Lagrangian averaging models:8the wall–normal components of the displacement-fluctuation covariance vanishes at the wall. However, the results obtained above, particularly共11兲, are not directly transferable to ␣=␣共x,y,z兲 because the Leray filter no longer commutes with differentiation.

An alternative that would permit no-slip boundary con-ditions for u˜ could be to add a pressure-correction procedure

to u˜i to make it divergence-free. Here u˜iwould be an

inter-mediate result, and the divergence-free velocity u˜_i*_{could be} obtained from u˜_i*_{= u˜ −⳵}

i␾, with ␾ a pressure-like scalar.

This procedure would involve solving the Poisson equation ⳵j

2_␾

=⳵j˜uj, and u˜i

*

would be used in the momentum equations

018103-3 Incompressibility of the Leray-␣model Phys. Fluids 18, 018103共2006兲

共1兲 instead of u˜i. The challenge of this approach would be to

attribute physical relevance to␾, other than that it enforces a divergence-free field for u˜*_.

Another option would be to use free-slip boundary con-ditions for u˜i, while applying no-slip boundary conditions for

ui. Using 共2兲 and 共7兲, only two terms are zero at the wall,

namely dzwˆ = dz

2

w˜ˆ = 0. Substituting d_z2w˜ˆ = 0 and the free-slip

boundary conditions dz˜ˆ = du zv˜ˆ = w˜ˆ = 0 into 共11兲 for n=1

im-mediately gives dzfˆ⬅0, from which it can be concluded that

this configuration results in an unconditionally divergence-free field for u˜i. To check this result, simulations have been

carried out with free-slip共Fsl兲 boundary conditions for u˜iand

the results are shown in Fig. 1. As expected, the divergence ⳵ju˜jis zero up to machine precision and the purely

redistribu-tive character of advection is restored.

One can arrive at free-slip boundary conditions more rigorously by using the boundary conditions in the wall– normal direction only, and have the wall-parallel boundary conditions follow from the restriction of incompressibility. In the wall–normal direction, the relevant condition is that the wall is impermeable, so that wˆ = 0 and also w˜ˆ = 0. Substituting wˆ = 0 and w˜ˆ = 0 into共7兲 yields d_z2_{w˜ˆ = 0, from which an} incom-pressible field can be constructed when dzfˆ = kxdz˜ˆ + ku ydz˜ˆv

= 0. The only solution yielding boundary conditions indepen-dent of wave number and other velocity components, is to set both dzu˜ˆ = 0 and dz˜ˆ = 0. Transforming back, it follows thatv

free-slip boundary conditions⳵z˜ =u ⳵z˜ = wv ˜ = 0 are required to

warrant a divergence-free u˜ifield under all circumstances. As

a physical interpretation, one could argue that the nonlocality of the filter causes the filtered velocity at the wall to be influenced by the flow in the bulk resulting in nonzero slip velocities. It is emphasized that the divergence condition共11兲 relies on the domain. Therefore, further analysis is required to verify whether the free-slip conditions are applicable for general domains as well.

It should be noted that if u˜iwere the Lagrangian average

of ui and uiwere to vanish at the wall, then its Lagrangian

average u˜i would also vanish at the wall.8 Thus, u˜i in the

Leray-␣ model deviates from the Lagrangian average veloc-ity, especially near the wall, if it needs a different boundary condition from ui. However, it can be expected that the

re-gion that will be influenced by the difference in boundary conditions is quite local and will vanish as␣→0. In some applications,共partial兲 free-slip boundary conditions have ac-tually been shown to enhance the global accuracy of the solution16on coarse grids.

As before, the Green’s functions will be used to study the influence of different boundary conditions for u共no slip兲 and u˜ 共free slip兲. The free-slip boundary conditions dz˜ˆ = 0 atu

z = ± H / 2 are of the Neumann type and the Green’s function

is GN共z,␨兲 =

冦

cosh关K共H/2 −␨兲兴cosh关K共H/2 + z兲兴 K␣2sinh共KH兲 , z艋␨, cosh关K共H/2 +␨兲兴cosh关K共H/2 − z兲兴 K␣2sinh共KH兲 , z艌␨.

冧

It is convenient to analyze the departure from no-slip bound-ary conditions for u˜ by the decomposition GN= GD+ G_⌬.

Here G_⌬represents the difference between the Green’s func-tion for free-slip共GN兲 and no-slip 共GD兲 boundary conditions,

yielding

G_⌬共z,␨兲 =cosh关K共z +␨兲兴

␣2_{K sinh共KH兲}. 共12兲

When␣ⰆH and thus KHⰇ1, the integration kernel G_⌬can be estimated by K−1␣−2exp共−Kd兲, and the influence region is proportional to K−1⬇␣ for ␣ sufficiently small. In the limit of␣→0, G_⌬converges to a Dirac delta function, and the wall slip velocity u˜→u=0. The influence region goes to

zero as well, from which it follows that this configuration globally converges to the Navier–Stokes equations.

This work is part of the research program of the Stichting voor Fundamenteel Onderzoek der Materie共FOM兲, which is financially supported by the Nederlandse Organi-satie voor Wetenschappelijk Onderzoek共NWO兲. The compu-tations were sponsored by the Stichting Nationale Comput-erfaciliteiten共NCF兲.

1_{A. Cheskidov, D. Holm, E. Olson, and E. Titi, “On a Leray-alpha model of}

turbulence,” Proc. R. Soc. London, Ser. A 461, 629共2005兲.

2_{S. Chen, C. Foias, D. Holm, E. Olson, E. Titi, and S. Wynne, “The}

Camassa-Holm equations and turbulence,” Physica D 133, 49共1999兲.

3_{S. Chen, C. Foias, D. Holm, E. Olson, E. Titi, and S. Wynne, “A}

connec-tion between the Camassa-Holm equaconnec-tions and turbulent flows in channels and pipes,” Phys. Fluids 11, 2343共1999兲.

4_{S. Chen, D. Holm, L. Margolin, and R. Zhang, “Direct numerical}

simula-tions of the Navier-Stokes alpha model,” Physica D 133, 66共1999兲.

5_{C. Foias, D. Holm, and E. Titi, “The Navier-Stokes-alpha model of fluid}

turbulence,” Physica D 152–153, 505共2001兲.

6_{D. Holm, “Variational principles for Lagrangian-averaged fluid}

dynam-ics,” J. Phys. A 35, 679共2002兲.

7_{D. Holm, J. Marsden, and T. Ratiu, “Euler-Poincaré models of ideal fluids}

with nonlinear dispersion,” Phys. Rev. Lett. 80, 4173共1998兲.

8_{D. Holm, “Fluctuation effects on 3D Lagrangian mean and Eulerian mean}

fluid motion,” Physica D 133, 215共1999兲.

9_{J. Marsden and S. Shkoller, “The anisotropic Lagrangian averaged Euler}

and Navier-Stokes equations,” Arch. Ration. Mech. Anal. 166, 27共2003兲.

10_{K. Mohseni, B. Kosovic, S. Shkoller, and J. Marsden, “Numerical}

simu-lations of the Lagrangian averaged Navier-Stokes equations for homoge-neous isotropic turbulence,” Phys. Fluids 15, 524共2003兲.

11_{B. Geurts and D. Holm, “Regularization modeling for large-eddy}

simula-tion,” Phys. Fluids 15, L13共2003兲.

12_{S. Stolz and N. A. Adams, “An approximate deconvolution procedure for}

large-eddy simulation,” Phys. Fluids 11, 1699共1999兲.

13_{S. S. Collis, “Monitoring unresolved scales in multiscale turbulence}

mod-eling,” Phys. Fluids 13, 1800共2001兲.

14_{E. D. Siggia, “High Rayleigh number convection,” Annu. Rev. Fluid}

Mech. 26, 136共1994兲.

15_{M. van Reeuwijk, H. J. J. Jonker, and K. Hanjalić, “Identification of the}

wind in Rayleigh-Bénard convection,” Phys. Fluids 17, 051704共2005兲.

16_{W. Layton, “Weak imposition of “no-slip” conditions in finite element}

methods,” Comput. Math. Appl. 38, 129共1999兲.

018103-4 van Reeuwijk, Jonker, and Hanjalić Phys. Fluids 18, 018103共2006兲