- the Application Album item concerned with a two-dimensional simulation of the transient flow behind a square-sectioned obstacle.
- and a paper and presentation concerning three-dimensional simulation of the same flow, by Ochoa and Fueyo, presented at the Melbourne PHOENICS User Conference of 2004.

Subgrid-scale (SGS) turbulence models are used in Large-Eddy Simulations (LES) of turbulent flows. The LES approach ( see for example Schumann [1991] ) is a combination of Direct Numerical Simulation (DNS) of turbulence and conventional turbulence modelling.

The large scales (eddies) are simulated directly, and the small (subgrid) scales are modelled with a SGS turbulence model. The flow variables are therefore decomposed in terms of:

- resolved (grid) scales and
- unresolved (subgrid) scales, rather than time-mean and fluctuating quantities.

The SGS stresses are similar to the Reynolds stresses which result from time- or ensemble- averaging; but they differ in that they are consequences of spatial averaging and go to zero if the mesh interval goes to zero.

PHOENICS provides the simplest and most commonly-used SGS model, namely the Smagorinsky [1963] eddy-viscosity model.

Since turbulence is inherently 3-dimensional, strictly speaking LES simulations must also be 3-dimensional, even though the time-mean flow may be of lower dimensionality.

Nevertheless, some workers have used SGS turbulence models in 2- dimensional problems with reasonable results (see for example, the results of Petersen [1992] for transient 2-phase flow in bubble columns, and those of Murakami et al [1993] and Ochoa and Fueyo [2004] for flow past square-sectioned obstacles).

In addition, time-dependence is an essential feature of LES, even when the time-average flow is independent of time.

In summary, true LES is a 3-dimensional time-dependent calculation in which the large-scale motions are calculated directly, and the SGS turbulence is modelled.

The size of such a calculation is inevitably large, and therefore computationally very expensive for practical engineering calculations. Consequently, on practical grounds, its use should only be considered in flows that are inherently 3D and transient, and in challenging applications for statistical turbulence models, such as for example, bluff-body aerodynamics, vortex shedding and aeroacoustics.

Another drawback is that problems arise from the need to introduce appropriate randomness into the initial and boundary conditions ( see for example Richter et al [1987] ). The advantage of LES is that it is largely free from arbitrary closure assumptions.

**Spatially-Averaged Navier Stokes Equations**

The LES model in PHOENICS resolves each flow variable into a large-scale (resolved) component, which is still time-dependent even in flows that are stationary in the average, and into a subgrid (unresolved) component.

Following Schumann [1975] and Ciofalo [1994],
explicit filtering of the local and instantaneous quantities is replaced by volume
averaging on each grid cell, so that the equations describing the resolved flow and scalar
(*e.g.*thermal) fields are:

RHO,t + (RHO*Ui),i = 0 .... (5.1.2.1)

(RHO*Ui),t + (RHO*Uj*Ui),j = - p,i + (2.*RHO*ENUL*Sij),j - (TAUij),j .... (5.1.2.2)

(RHO*Cp*T),t + (RHO*Cp*Uj*T),j = (k*T,j),j - (Qj),j ....(5.1.2.3)

in which TAUij are residual stresses and Qj are residual heat fluxes which have to be expressed by means of a closure model (for details, see Ciofalo [1994]).

**The Subgrid-Scale Stresses**

The SGS stresses TAUij are calculated from:

TAUij = - 2.*RHO*ENUT*Sij .... (5.1.2.4)

where RHO is the density and ENUT is the SGS kinematic eddy viscosity, and Sij is the mean rate-of-strain tensor for the resolved scales:

Sij = 0.5*(dUi/dxj+dUj/dxi) .... (5.1.2.5)

**The Subgrid-Scale Model**

PHOENICS is equipped with three variants of the Smagorinsky model for the SGS eddy viscosity ENUT: namely:

- the basic Smagorinsky model;
- the Smagorinsky Wall-Damping (WD) Model; and
- the Smagorinsky Van-Driest Wall-Damping (VDWD) Model.

ENUT = (LM)**2 * (2*Sij*dUi/dxj)**0.5 .... (5.1.2.6)

where LM is the characteristic length scale of the resolved eddies. The three model variants differ in their prescription for LM, as described below.

The __basic Smagorinsky model__ computes LM from:

LM = Cs*H .... (5.1.2.7)

where Cs is the Smagorinsky constant (=0.17), and H is a representative mesh interval.

The Smagorinsky constant Cs is known to vary with application, with values reported in the literature ranging from 0.1 to 0.25. The default value in PHOENICS is 0.17.

For the __Smagorinsky WD Model__,

LM = MIN(Cs*H,k*WDIS) .... (5.1.2.8)

where k is von Karman's constant (=0.41); and WDIS is the normal distance to the nearest wall,
calculated *via* the solution of an elliptic differential equation for the variable LTLS (see the
PHOENICS Encyclopaedia entry DISWAL).

According to Schumann [1991], the MIN(...) operator makes the model valid for coarse grids, in particular near walls, where the mixing length (k*WDIS) may become smaller than the mesh interval H.

The __Smagorinsky VDWD model__ calculates LM from:

LM = Cs*H*FMU .... (5.1.2.9)

where FMU is the Van Driest [1956] damping function, given by:

FMU = (1.-EXP[-(Y+/26.0)**m])**n .... (5.1.2.10)

The dimensional wall distance y+ is computed using the normal distance to the nearest wall, WDIS. The empirical exponents m and n are defaulted to unity in PHOENICS.

Equation (5.1.2.9) is an alternative to equation (5.1.2.8) for partially taking into account wall effects by appropriately reducing the length scale Cs*H in the proximity of walls.

**The Grid Filter Width, H**

The representative mesh interval H is computed from:

H = V**1/3 .... (5.1.2.11)

where V is the local cell volume. An alternative commonly encountered in the literature is also provided as an option in PHOENICS:

H = [ S (DXi**2)/N ]**(1./N) .... (5.1.2.12)

where the summation is over the three coordinate directions and DXi is the cell width in the direction i.

**The Subgrid-Scale Heat Fluxes**

The SGS heat fluxes are modelled by:

Qj = - Kturb*T,j .... (5.1.2.13)

where Kturb=RHO*ENUT*Cp/PRT(TEM1) and PRT(TEM1) is the subgrid Prandtl number which is the LES analog of the turbulent Prandtl number.

In the LES literature, a broad range of values has been proposed for this parameter, from 0.25 to 0.9 ( see for example Ciofalo [1994]).

The default in PHOENICS is the **Smagorinsky WD model**, which may be activated from the
Q1 input file *via* the PIL command TURMOD(SGSMOD), which is equivalent to:

ENUT=GRND2;EL1=GRND10;GENK=T;EL1A=0.17;EL1B=0.41;DISWAL

Here, EL1A equals Smagorinsky's constant and EL1B von Karman's constant.

The **basic Smagorinsky model** is activated by setting

TURMOD(SGSMOD);EL1B=0.0

and if no walls are present, by deactivating the PIL command DISWAL *via*
SOLUTN(LTLS,N,N,N,N,N,N) and SOLUTN(WDIS,N,N,N,N,N,N).

The **Smagorinsky VDWD model** is activated by setting

TURMOD(SGSMOD);EL1B=0.0;EL1C=1.0; EL1D=1.0

where EL1C=n and EL1D=m in equation (5.1.2.10). If the user puts STORE(FMU) in the Q1 input file, the field values of the damping function FMU are written to the PHOENICS RESULT and PHI files.

The default grid filter width H is given by (5.1.2.11); but the alternative (5.1.2.12)
can be activated by setting the Smagorinsky coefficient to a negative value, *e.g.*
EL1A=-0.17 .

Compared with conventional
turbulence modelling, LES requires finer spatial resolution and smaller time steps, and consequently significantly more computer time.

Typically, the time step should be in the range 1/200 to 1/50 of the large-eddy turnover
time.

PHOENICS uses fully-implicit backward time-stepping (1st order Euler) and hybrid differencing for the combined effects of advection and diffusion.

This scheme reduces to the central-differencing scheme recommended by most workers for LES when the cell Peclet number is less than 2. Otherwise, it reverts to the upwind scheme. Therefore, smaller time steps and finer meshes are desirable.

Higher-order convection schemes are available in PHOENICS as built-in options, such as for example MINMOD and the Linear Upwind Scheme.

There is no built-in option for a higher-order temporal differencing scheme,
although PHOENICS is flexible enough for the user to introduce his own temporal scheme
such as Crank- Nicholson.

For example,
Ochoa and Fueyo [2004] implemented the 3rd-order
implicit Adam-Moulton scheme in PHOENICS via GROUND coding.

The only purpose of using LES is to simulate turbulent flows more realistically than can conventional models; but they can do so only if:

- the Smagorinsky constant is not too high (for example, in duct flows values in excess of 0.12 have been known to suppress turbulence completely);
- the spatial grid is not too coarse,
- the time step is not too large.

What 'not too coarse' and 'not too large' imply, in any newly-considered flow, can determined only by trial-and-error, *i.e.* by performing the calculations with a range of mesh sizes, time steps and Smagorinsky constants, and comparing the results.

Such systematic grid-refinement is rarely performed, because of its expense; indeed the need for it is often disregarded by LES practitioners. There is not even any consensus as to how the results with different numerical inputs should be compared.

However, the following advice appears to be reasonable:

- If, as the calculation proceeds, fluctuations disappear, the mesh is
**certainly**too coarse, or the Smagorinsky constant too large. - Since it is the fluctuating nature of turbulence which LES seeks to capture, there is little point in judging its success in terms of
**time-average**quantities. Yet this is the common practice. - The most informative outcomes of LES predictions would be the
**'probability-density functions'**, at various locations in the flow, of such quantities as temperature, concentration, and velocity components. - When these functions cease to vary significantly with increasing refinement, it can be presumed that the 'not too ..' requirements have been met.
- What is the correct value of the Smagorinsky constant can then be deduced by comparing the computed probability-density functions with such reliable experimental or DNS data as exist for the case in question.

M.Ciofalo, 'Large-Eddy Simulation: A Critical Survey of Models and Applications', In Advances in Heat Transfer, Vol.25, p321, Academic Press, (1994).

J.S.Ochoa and N.Fueyo, 'Large-Eddy Simulation of the flow past a square cylinder', PHOENICS 10th International User Conference, Melbourne, Australia, CHAM Limited, London (2004).

K.O.Petersen, 'Etude experimentale et numerique des ecoulements diphasiques dans les reacteurs chimiques', PhD Thesis, L'Universite Claude Bernard - Lyon, (1992).

K.Richter, R.Friedrich and R.Schmitt, 'LES of turbulent wall boundary layers with pressure gradient', p22-3-1, 6th Symp. on Turbulent Shear Flows, Sept. 7-9, Toulouse, France, (1987).

S.Murakami, A.Mochida, W.Rodi and S.Sakamoto, 'Large-Eddy Simulation of turbulent vorttex shedding past 2D square cylinders', p113-120, FDE-Vol.162, Engineering applications of LES, ASME, USA, (1993).

U.Schumann, 'Subgrid scale model for finite-difference simulations of turbulent flows in plane channels and annuli', J.Comp.Phys., Vol.18, p376, (1975).

U.Schumann, 'Direct and Large-Eddy simulation of turbulence', Introduction to the modelling of turbulence, VKI Lecture Series 1991-02, March 18-21, VKI, Belgium, (1991).

J.Smagorinsky, 'General circulation experiments with the primitive equations', Monthly Weather Review, Vol.93, No.3, p99, (1963).

E.R.Van Driest, 'On turbulent flow near a wall', J.Aerospace.Sci., Vol.23, No.11, p1007-1011, (1956).

See also numerous papers in:

Hanjalic K, Nagano Y and Jakirlic S (Eds.)

"Turbulence, Heat and Mass Transfer",

Begell House Inc., New York and Walingford (UK),

2009 .