Variable density formulation of the dynamic Smagorinsky model Randy McDermott April 15, 2004

1

Introduction

There are no new ideas presented in this document. The aim here is to clarify the modeling approach which is more or less standard in the LES literature, and to conform the nomenclature in the literature (which varies considerably, especially considering the sign convention difference between mechanical and chemical engineers for the stress) to that used in our (CRSIM) codes. Additionally, it never hurts to see things explained in a different way. Creation of this document was as much an education for the author as it hopefully will be for any new graduate students to our group. What follows is basically a mixture of two papers: the original paper by Moin et al. [6], and the review by Martin et al. [5]. The formulation presented in these papers is suitable for high Mach number flows. For the low Mach, variable density case the procedure is little changed. There simply exists a modeling choice, of whether to explicitly model the isotropic component of the subgrid stress, or to absorb this force into the pressure term. We will first derive the deviatoric part of the subgrid stress for variable density flows. The most significant departure from the incompressible case is the form of the Leonard stress. Then, we will show how to dynamically model the isotropic component based on the approximation of Yoshizawa [9].

2

Newton’s law of viscosity

The isotropic part of a stress is indistinguishable from pressure in a hydrodynamic sense. In contrast, the deviatoric stress is what is usually associated with Newton’s law of viscosity. Let us start from what Newton actually said about the stress and see what this implies. All Newton said was that the surface stress is a linear function of the velocity gradient. Let the total surface stress be σij . For isotropic fluids (the fluid has the same properties no matter which direction one is looking from a point in the fluid) with no angular acceleration (this implies the stress must be symmetric, else there would be no surface moment to counter balance the stress and a fluid element would accelerate to infinite angular velocity) it can be shown [7] that 1

the linearity assumption amounts to the following statement about the stress (throughout this document the summation convention over repeated indices is assumed), σij = a1 δij + a2 Sij + a3 Skk δij , where the strain rate, Sij (the symmetric part of the velocity gradient tensor), is given by, µ ¶ 1 ∂ui ∂uj , Sij ≡ + 2 ∂xj ∂xi

(1)

(2)

and the coefficients, a1 , a2 , and a3 , are to be defined. It is easy to set the first coefficient, a1 , because if the fluid is not moving the stress must be equal to the equilibrium thermodynamic pressure. Hence, a1 = p. The other two coefficients (respectively, the first and second coefficients of viscosity) are typically written as, a2 = −2µ, and a3 = −λ. Now, the mechanical pressure, pm is defined to be the average of the trace (the sum of the diagonal components) of the total stress tensor, pm ≡

1 σii . 3

(3)

Contracting on i and plugging this definition into Equation 1 gives, 3pm

=

pm

=

p − pm

=

3p − 2µSkk − 3λSkk , µ ¶ 2 p− µ + λ Skk , 3 µ ¶ 2 µ + λ Skk . 3

(4)

Assuming the thermodynamic and mechanical pressures are equal amounts to setting the second coefficient to, λ = −2/3 µ. This is known as Stoke’s assumption. Kinetic theory suggests this is exactly true for monatomic gases, and it is often argued that this assumption is valid for all Newtonian fluids, except for extremely high Mach cases. With all this, Newton’s law for the surface stresses becomes, σij

2 = pδij − 2µSij + µSkk δij , 3 µ ¶ 1 = pδij − 2µ Sij − Skk δij , 3 = pδij + τij ,

(5)

where τij is known as the deviatoric stress. The Navier-Stokes equations can then be written as, ρ

Dui ∂σij =− Dt ∂xj

(6)

For variable density flows it is not permissible to argue that the trace of the strain rate is small, and thence ignore it when computing the viscous stress tensor. This term is balanced by an accumulation of density in the continuity equation, which as we well know in combustion simulations is not negligible. Hence, 2

unless one is simulating incompressible flow in the sense that the divergence of the velocity field is identically zero, it is necessary to retain the full form of the deviatoric stress. One further comment: the pressure which is obtained from a projection method is the hydrodynamic pressure. This pressure field has the same gradient as the mechanical pressure field, but its absolute value must be obtained by defining a reference pressure within the domain.

3

The deviatoric subgrid stress

For variable density flows, the filtered Navier-Stokes equations can be written as, ∂ρui ∂ρui uj ∂σ ij + =− . ∂t ∂xj ∂xj

(7)

The Favre filtered field, u e, is defined such that, ρe u = ρu .

(8)

Inserting this definition into Equation 7 gives, ∂ρe ui ∂ρug i uj + ∂t ∂xj

=

−

∂σ ij , ∂xj

∂ρe ∂ρe ui ui u ej + ∂t ∂xj

=

−

sgs ∂τij ∂σ ij − , ∂xj ∂xj

(9)

where the subgrid stress is defined as, sgs τij ≡ ρ (ug ei u ej ) . i uj − u

(10)

Now, we wish to model the subgrid stress by direct analogy to Newton’s law of viscosity. This is called the gradient diffusion hypothesis. However, at this point, based on the definition (10) the subgrid stress is not strictly deviatoric. Therefore, we will break this stress up into two tensors, one deviatoric (denoted superscript D), and one isotropic (denoted superscript I): sgs D τij = τij + τ I δij ,

(11)

1 sgs sgs D τij = τij − τkk δij . 3

(12)

where,

We can model the deviatoric subgrid stress by, D τij

¶ µ 1e e ≈ −2ρνt Sij − Skk δij , 3

(13)

where the filtered strain rate is given by, 1 Seij ≡ 2

µ

∂e ui ∂e uj + ∂xj ∂xi 3

¶ .

(14)

Note that this is not yet the “Smagorinsky model” for the subgrid stress. The Smagorinsky model is a model for the eddy viscosity [8], e , νt = (Cs ∆) |S| 2

(15)

e ≡ (2Seij Seij )1/2 . where the magnitude of the strain rate is, |S| Putting this all together, in anticipation of developing the dynamic procedure, let us write the model for the deviatoric stress as, 2

D τij ≈ −2 (Cs ∆) βij ,

where,

4

¶ µ 1e e e βij = ρ|S| Sij − Skk δij . 3

(16)

(17)

The dynamic procedure

Dynamic determination of the model constant, Cs , is possible by employing an additional filter with a filter width larger than the mesh spacing, and then presuming that the constant should be the same for the b > ∆, to the filtered solution at both scales. Let us apply a filter (call this the “test” filter) of width, ∆ Navier-Stokes equations,

ci d bij ∂ ρu ∂ ρu ∂σ i uj + =− . ∂t ∂xj ∂xj

(18)

In analogy to Favre filtering, let us define the following decomposition (The number of over-symbols required to describe all the different filters gets ridiculous. The ˘ symbol is adopted here from Martin et al. Unfortunately, this symbol has no “wide” version in LATEX), ˘ ≡ ρc b ρu e u.

(19)

We can now write the test filtered Navier-Stokes equations as, ˘ ˘ei ∂b ρu ρug ∂b i uj + ∂t ∂xj ˘ ˘ei u ˘ej b b ∂ ρu ei ∂ ρu + ∂t ∂xj

= =

bij ∂σ , ∂xj bij ∂σ ∂Tij − − , ∂xj ∂xj −

(20)

where the “subtest” stress is defined as, ³ ´ ˘ ˘ei u ˘ej . Tij ≡ b ρ ug i uj − u The deviatoric part of the subtest stress should be modeled as, ¶ µ ³ ´2 1˘ ˘e ˘e b b TijD ≈ −2 Cs ∆ ρ|S| S ij − Sekk δij . 3

4

(21)

(22)

Now then, the Germano identity [2] gives us something that looks like a Leonard stress, sgs Lij = Tij − τd ij

³ ´ ˘ d ei u ˘ei u ˘ej − ρ (ug ρ ug ej ) , = b i uj − u i uj − u ³ ´ ³ ´ ˘ ˘ ˘ei u ˘ej − b = b ρ ug ρ ug ei˘u ej , i uj − u i uj − u ³ ´ ˘ei u ˘ej . = b ρ u ei˘u ej − u

(23)

We will explain why we do not need to worry about just taking the deviatoric part of Lij in a moment. For now, let us get (23) into a more usable form. Utilizing the Favre definitions, (8) and (19), we get the following for the Leonard term, Lij

= = =

ci ρd ρud ρu uj i ρuj −b ρ , b ρ ρ ρ b ρ ci ρd ρud ρu uj i ρuj − , b ρ ρ c c ρe ui ρe uj d ρe ui u ej − , b ρ ρ

(24)

which is the form typically seen in the literature. LATEX has a hard time covering the entire term with the “wide” version of the hat, but please note that the entire first term of Equation 24 is test filtered. The Leonard term is a defined quantity that is computable from resolved LES values. If we now look at the model for the Germano identity (the deviatoric part) we have, µ ¶ ³ ´2 1 ˘e ˘e ˘e 2 b D c D ≈ −2 C ∆ b b − τ = T LD ρ | S| S − S δ s ij kk ij + 2 (Cs ∆) βij . ij ij ij 3

(25)

Note that in actuality the entire last term should be test filtered, since Cs is not constant in space (which is a main feature of the dynamic model). However, without pulling the length scale out of the filter operation, it is difficult, if not impossible, to compute a value for the model parameter. Let us now rearrange (25) into the following form to facilitate coding, 2

D LD ij = (Cs ∆) Mij ,

where, D Mij

µ µ ¶¶ 1 ˘ ˘ ˘ b e e e = 2 βij − αb ρ|S| S ij − S kk δij , 3

(26)

(27)

2 b and α = (∆/∆) . So, for a test filter width that is two times the grid width we should have, α = 4. However,

please note that the method of discrete quadrature used in computing the test filtered values can have a significant effect on the results [4]. In this author’s work with incompressible decaying isotropic turbulence, in which the second order kinetic energy conserving scheme of Harlow and Welch [3] is used for the convective operator, and the trapezoid rule is used for quadrature, a value of α = 4 was found to work quite well.

5

The strain rates which should be used in Equation 27 are tricky, and it is not obvious that these have been computed correctly in the previous literature [6]. The Favre filtered strain rate at the test scale should be computed as follows, ˘ Seij

=

=

! ˘ei ˘ej ∂u ∂u + , ∂xj ∂xi µ ¶  µ ¶ c c ρ˜ uj ρ˜ ui ∂ ∂ bρ bρ  1  . + 2  ∂xj ∂xi  1 2

Ã

(28)

Then we define the magnitude of this strain rate as usual, ³ ´ ˘e ˘ ˘ 1/2 |S| ≡ 2Seij Seij .

(29)

D We now have a way to compute all the necessary values for Lij and Mij from known LES quantities. If D we right multiply Equation 26 by Mij we have a way to solve for the length scale, 2

(Cs ∆) =

D LD ij Mij . DM D Mij ij

(30)

We mentioned earlier that we would discuss the fact that it is unnecessary to explicitly compute the deviatoric part of the Leonard term. This is because, as luck would have it, there exists a tensor identity D D = LD that says the following, Lij Mij ij Mij . Here is the proof (due to Stas): D Lij Mij

= = =

D LD ij Mij

µ ¶ 1 Lij Mij − Mkk δij , 3 1 Lij Mij − Lij δij Mkk , 3 1 Lij Mij − Lqq Mkk . 3

µ ¶µ ¶ 1 1 Lij − Lqq δij Mij − Mkk δij , 3 3 1 1 1 = Lij Mij − Lij δij Mkk − Mij δij Lqq + δij δij Lqq Mkk , 3 3 9 1 1 1 = Lij Mij − Mkk Lqq − Lqq Mkk + Mkk Lqq , 3 3 3 1 = Lij Mij − Lqq Mkk . 3

(31)

=

(32)

Equations 31 and 32 are equal and so there is no need to go to the trouble of subtracting the isotropic part out of Lij . Further notes on implementation, which apply for both the incompressible and variable density formulation: It has been found from experience that, to maintain stability, the length scale should be averaged

6

over some homogeneous domain. This author suggests using an average over the test filter width. Hence, the model parameter is ultimately computed in the following way, 2

(Cs ∆) =

D hLij Mij i . DM Di hMij ij

(33)

Here the brackets h · i denote an arbitrary spatial average. This is the standard notation of the literature. In practice, it seems that these brackets are license to do anything within one’s power to maintain a stable code. In our codes we are explicit about using the test filter for the averaging. In channel flow it is standard practice to average over a plane normal to the wall, and hence Cs = Cs (y, t). Other authors performing LES of decaying isotropic turbulence have averaged over the entire domain [1]. In this case Cs = Cs (t). It is also standard practice to “clip” the eddy viscosity. In theory a negative value of the eddy viscosity would indicate a backscatter of energy from the unresolved scales to the resolved scales. If this sounds D dangerous from a stability perspective. . . it is. The simple way around this is to set Cs = 0; if hLij Mij i < 0.

5

The isotropic subgrid stress

Choosing not to model the isotropic part of the subgrid stress amounts to the following LES formulation, ui ui u ej ∂ρe ∂ρe + ∂t ∂xj

= −

D ∂τij ∂(p + τ I ) ∂τ ij − − , ∂xi ∂xj ∂xj

ui ui u ej ∂ρe ∂ρe + ∂t ∂xj

= −

D ∂τij ∂ph ∂τ ij − − . ∂xi ∂xj ∂xj

(34)

This is perfectly legitimate, and in fact is recommended for low Mach flows, especially if a “pressure equation” type of projection is used (i.e. ph = 0 for the predictor step). One may find, however, that the “pressure correction” method helps the speed of the linear pressure solve, and so the better the guess for the pressure I in the predictor step, the better the correction algorithm will work. In this case it makes sense to model τkk .

It is a rather simple procedure, and so we include the model here for completeness, although we do not employ this in our codes. Yoshizawa’s model is the following, sgs e2. τkk ≈ 2CI ∆2 ρ|S|

(35)

sgs Note that as we have written things above, τ I = (1/3)τkk . Also note that for some reason the original

authors have decided not to square the constant, as is done in the Smagorinsky model. No matter, we are just to be aware of it. At the test scale the model for the trace of the stress tensor is, ˘e 2 b 2b Tkk ≈ 2CI ∆ ρ|S| .

(36)

c c ρe uk ρe uk sgs d Lkk = Tkk − τd = ρe u u e − . k k kk b ρ

(37)

The Germano identity gives,

7

The model gives, ˘e 2 d b 2b e 2 = CI Mkk , Lkk ≈ 2CI ∆ ρ|S| − 2CI ∆2 ρ|S|

(38)

˘e 2 d b 2b e 2 . Then, where, Mkk = 2∆ ρ|S| − 2∆2 ρ|S| CI =

hLkk i . hMkk i

(39)

Acknowledgements Funding for this work was provided through the Department of Energy Computational Science Graduate Fellowship (DE-FG02-97ER25308). Thanks to Stas Borodai for educational discussions.

References [1] Stephen M. de Bruyn Kops. Numerical simulation of non-premixed turbulent combustion. PhD thesis, The University of Washington, 1999. [2] M. Germano, U. Piomelli, P. Moin, and W. Cabot. A dynamic subgrid-scale eddy viscosity model. Phys. Fluids A, 3(7):1760–1765, 1991. [3] F.H. Harlow and J.E. Welch. Numerical calculation of time-dependent viscous incompressible flow of fluid with free surface. Phys. Fluids, 8:2182, 1965. [4] T.S. Lund. On the use of discrete filters for large-eddy simulation. Technical report, Center for Turbulence Research – Annual Research Briefs, 1997. [5] M. Pino Mart´in, U. Piomelli, and G. Candler. Subgrid-scale models for compressible large-eddy simulations. Theoret. Comput. Fluid Dynamics, 13:361–376, 2000. [6] P. Moin, K. Squires, W. Cabot, and S. Lee. A dynamic subgrid-scale model for compressible turbulence and scalar transport. Phys. Fluids A, 3(11):2746–2757, 1991. [7] Ronald L. Panton. Incompressible Flows. John Wiley and Sons, 2nd edition, 1996. [8] J. Smagorinsky. General circulation experiments with the primitive equations. Monthly Weather Review, 91(3):99–106, 1963. [9] A. Yoshizawa. Statistical theory of compressible turbulent shear flows, with the application of subgrid modeling. Phys. Fluids A, 29:2152–2164, 1986.

8