**2. Governing equations and numerical methods**

In Large Eddy Simulation (LES) a filtering operation is applied to separate the large scales from the small scales (Leonard 1974). In general, a filtered variable can be written as

$$\overline{f}(\mathbf{x}) = \int\_{D} f(\mathbf{x'}) G(\mathbf{x}, \mathbf{x'}; \overline{\mathbf{A}}) d\mathbf{x'} \tag{2.1}$$

where *G* is the filter kernel and *D* is the filtering domain. The filter is characterized by a filter width . The corresponding wave number *<sup>c</sup> k* is called as the filter cut-off wave number.

For our study, the fluid is assumed to be incompressible; the viscosity is constant; there are no body forces; and the flow is initially homogenous, isotropic, i.e. there are no mean velocity gradients. So the incompressible Navier-Stokes equations after applying a low-pass filter can be written as

$$
\frac{\partial \overline{u\_i}}{\partial t} + \frac{\partial \overline{u\_i u\_j}}{\partial \mathbf{x}\_j} = -\frac{1}{\rho} \frac{\partial \overline{P}}{\partial \mathbf{x}\_i} + \frac{\partial}{\partial \mathbf{x}\_j} \left| \nu \left( \frac{\partial \overline{u\_i}}{\partial \mathbf{x}\_j} + \frac{\partial \overline{u\_j}}{\partial \mathbf{x}\_i} \right) \right| - \frac{\partial \mathbf{r}\_{ij}}{\partial \mathbf{x}\_j} \frac{\partial \overline{u\_i}}{\partial \mathbf{x}\_i}.\tag{2.2}
$$

Above equations are also called the incompressible LES equations. The *u*, *P*, , are the velocity, pressure, density and kinematic viscosity, respectively. *ij* is the subgrid stress (SGS) tensor

$$
\pi\_{i\dot{j}} = \overline{u\_i u\_j} - \overline{u\_i}\,\overline{u\_j}\,. \tag{2.3}
$$

It represents the effect of the unresolved (small) scales. It is the only unclosed term in the above LES equations (2.2) and should be parameterized in terms of the resolved (large) scales.

In order to isolate other effects, the simplest homogenous, isotropic turbulence is chose as our simulation case. The advantage is that we can obtain the statistical quantities of this turbulence easily in spectral space, such s energy spectrum, total kinetic energy etc. So in

Study of Some Key Issues for

**3. SGS model** 

calculate the *ν*t.

*al.* (Piomelli *et al.*1988)).

resolved scales. The small ones, *ui*

the eddy viscosity model briefly. The eddy viscosity models assume:

the equation and needed to be modelled, i.e. *ij*

Applying LES to Real Engineering Problems 31

*<sup>n</sup> u u <sup>n</sup> <sup>N</sup> t*

2 *n n u u <sup>n</sup> N N*

where *n* and *n*+1 represent the different time steps, and \* represent the middle step variable.

In LES, only the large, energy carrying scales of turbulence are computed exactly. Specify in LES equation (2.2), the large scales are the filtered velocities, *ui* , which are also called the

The SGS model is the key issue in LES. Since only large scales are resolved in LES, the energy transfer from large scales to small scales is cut off. The energy will accumulate at the cut-off wave number and lead to the unphysical solution. So the main role of SGS model is to provide necessary small scales dissipation and thus remove the accumulated energy. There are many different approaches for the modelling of the SGS stress tensor. Traditionally they are divided into three main categories: eddy viscosity models, similarity models and mixed models. Discussion of standard LES models can be found in some review paper, such as Piomelli (Piomelli 1999), Mathew (Mathew 2010) etc. Below we only discuss

2 13 *ij uu uu S <sup>i</sup> <sup>j</sup> <sup>i</sup> <sup>j</sup> t ij ij kk*

and *ν*t is the eddy viscosity. Like RANS, equation (3.1) was developed by analogizing to the molecular viscosity. So different eddy viscosity models are actually different methods to

The Smagorinsky model (Smagorinsky 1963) is perhaps the most successful SGS eddy viscosity models, which takes eddy viscosity proportional to the product of *Δ*2 and *s* ,

<sup>2</sup>

magnitude of the strain-rate tensor. By choosing different *C*s for various flows, Smagorinsky model has been used with considerable success. For isotropic decaying turbulence, the value of the Smagorinsky constant is taken to be around 0.18∼0.23 (Lilly (Lilly 1996)), but in shear flow or near boundaries, *C*s must be decreased and values 0.06∼0.1 are preferred (Piomelli *et* 

*j i*

*j i u u*

*x x* 

1 2

which relate the SGS stresses to the large scale strain-rate tensor *S*ij , where *S*ij is

*ij*

where *C*s is called the Smagorinsky constant, *Δ* is the grid size and

*s*

<sup>1</sup> 1

*t*

Eq. (2.8) is also called the predictor step and Eq. (2.9) is called the corrector step.

(2.8)

(2.9)

, (unresolved, or subgrid scales) have been removed from

in equation (2.2).

 

*t s C s* (3.3)

(3.1)

(3.2)

1/2 2 *ij ij s ss* is the

such case it is convenience to write the governing equations in spectral space. Note the continuity equation can be combined with the pressure term through the projecting operation (Lesieur 1990). So the governing equation in spectral space can be simplified as

$$
\hat{\mathbf{i}}\left(\frac{\partial}{\partial t} + \nu \boldsymbol{k}^2\right) \hat{\boldsymbol{\mu}}\_i(\mathbf{k}, t) = -P\_{im}(\mathbf{k}) i k\_j \widehat{\boldsymbol{\mu}\_j \hat{\boldsymbol{\mu}}\_m}(\mathbf{k}) \,. \tag{2.4}
$$

where ' ' means the Fourier transform, the tensor <sup>2</sup> ( ) *P k im im i m* **k** *k k* is called the projection operator, which ensures the continuity equation automatically satisfied. And the *k* **k** .

For spatial discretization, a computation method similar to Rogallo's (Rogallo 1981) is adapted here. For the viscous term in the left hand side of equation (2.4), Rogallo proposed an integrated factor method which can solve it analytically.

$$\frac{\partial}{\partial t}\bigg[e^{\nu k^2 t}\hat{u}\_i(\mathbf{k}, t)\bigg] = -e^{\nu k^2 t}P\_{im}(\mathbf{k})N\_m(\mathbf{k})\,. \tag{2.5}$$

So the only term needed to be discretized is the nonlinear term in the right hand side, () () *m jj <sup>m</sup> N ik u u* **k k** , which usually is solved by high order spectral scheme. But for engineering problem, spectral method is not available at most cases. Finite difference scheme or finite volume scheme is used instead. Among them, Padé compact scheme is widely adapted due to its flexibility in handling complex geometry and to obtaining high order. For one dimensional derivative, the Padé scheme proposed by Lele (Lele 1992) can be expressed as

$$a x\_{i-1}^{'} + f\_i^{'} + a f\_{i+1}^{'} = a \frac{f\_{i+1} - f\_{i-1}}{2 \Lambda} + b \frac{f\_{i+2} - f\_{i-2}}{4 \Lambda} \,. \tag{2.6}$$

Different coefficient defines different order of compact scheme. The highest order is 6th for 3 point stencil. The parameters of Lele (Lele 1992) are shown in Table 2.1.


Table 2.1. Parameters for Padé compact scheme.

For the temporal advancement of the nonlinear term, an explicit second-order Runge-Kutta scheme, also known as predictor-corrector scheme, is used. It's simple and efficient. Briefly, equation (2.5) with only nonlinear term on the right hand side can be seen as the following form

$$\frac{\partial}{\partial t}\mu = \mathbf{N} \tag{2.7}$$

where *N* represents the nonlinear term. By applying second-order R-K scheme to above equations, we get

$$\frac{\mu^\*-\mu^n}{\Delta t} = N^n \tag{2.8}$$

$$\frac{\mu^{n+1} - \mu^n}{\Delta t} = \frac{1}{2} \left( N^n + N^\star \right) \tag{2.9}$$

where *n* and *n*+1 represent the different time steps, and \* represent the middle step variable. Eq. (2.8) is also called the predictor step and Eq. (2.9) is called the corrector step.
