**2. Development of RST modelling**

Early work leading to the development of Reynold-stress transport (RST) modelling was mainly theoretical, due to the relative complexity of this level of modelling compared to the available computational capabilities of the time. Chou (1945) constructed a formal solution to the fluctuating pressure Poisson equation that is the basis for current models of the pressure–strain-rate correlation. Later Rotta (1951), laid the foundation for Reynolds stress transport modelling by being the first to develop a closed model of all the terms in the exact equation (Speziale, 1991). Because of limited computational capability at the time, successful computations were not carried out until several decades later (Speziale, 1991). Another important development came when the continuum mechanics community speculated on the potential similarity between turbulent flow and the flow of non-Newtonian fluids (Gatski,

The modelling approach used for the various terms appearing in the exact Reynolds stress

Reynolds Stress Transport Modelling 5

Turbulent flows are characterised by highly fluctuating velocity, pressure, and other field variables. One approach for dealing with this fluctuating nature of the flow, the one most widely used by engineers, is to work with an averaged form of the basic equations. In *Reynolds averaging* the instantaneous flow variables are decomposed into an average quantity and a

*U<sup>i</sup>* = *Ui* + *ui*

where capital letters denote averaged quantities, and small letters denote purely fluctuating quantities. The averaging can be either over time or over a repeated realisation of an experiment with the same nominal conditions. The latter, *ensemble* averaging, will be implied in the following, to allow for temporal variations of mean quantities. When this decomposition is substituted into the Navier Stokes equations for incompressible flow, and the result is ensemble averaged, one obtains the Reynolds-averaged Navier-Stokes (RANS)

*<sup>P</sup>* <sup>=</sup> *<sup>P</sup>* <sup>+</sup> *<sup>p</sup>* , (1)

<sup>−</sup> *<sup>∂</sup>uiuj ∂xj*

= 0 . (3)

= 0 , (4)

. (2)

transport equation are briefly reviewed.

*∂Ui <sup>∂</sup><sup>t</sup>* <sup>+</sup> *Uj*

and averaging is applied, one obtains for the mean flow

the fluctuating velocity is obtained

fluctuating field, the result can be written as

*Dui Dt* <sup>=</sup> <sup>−</sup><sup>1</sup> *ρ ∂p ∂xi* − *uj ∂Ui ∂xj*

*∂Ui ∂xj*

<sup>=</sup> <sup>−</sup><sup>1</sup> *ρ ∂P ∂xi* + *ν ∂*<sup>2</sup>*Ui ∂x*<sup>2</sup> *j*

When the decomposition is substituted into the continuity equation for incompressible flow,

If this is subtracted from the instantaneous continuity equation, the continuity condition for

meaning that both the mean and fluctuating velocity fields are individually divergence free. The last term in the RANS equation (2) contains the Reynolds stress tensor *uiuj*. Thus the averaging process introduced a new unknown tensor term, and the set of equations is no longer closed. This is called the *closure* problem of averaging approaches. The task of turbulence modelling is to construct appropriate models for these stresses that relate them to the mean flow quantities, and thus to construct a closed set of equations allowing numerical solutions to be obtained. An additional implied objective in the engineering context is for the models to be as computationally inexpensive as possible while being able to reproduce the behaviour and phenomena of relevance to the problem in question, at the required level of

A transport equation for the fluctuating velocity can be obtained by subtracting the RANS equation (2) from the Navier-Stokes equation. Using the divergence-free property of the

> <sup>−</sup> *<sup>∂</sup> ∂xj*

*uiuj* − *uiuj* − *ν*

*∂ui ∂xj*

(5)

*∂Ui ∂xi*

*∂ui ∂xi*

**3.1 Basic equations of turbulent flow**

fluctuation. Thus,

equations

accuracy.

2004). This meant that tensor representation results from the continuum mechanics literature could be used to formulate expressions for the Reynolds-stress tensor, as first proposed by Rivlin (1957). These ideas were then expanded by Crow (1968; 1967), and Lumley (1967; 1970). Computational work accelerated in the 1970's with the works of Daly & Harlow (1970), Reynolds (1970), Donaldson (1971), Naot et al. (1972), Hanjali´c & Launder (1972), and Lumley & Khajeh-Nouri (1974). In a landmark paper, Launder, Reece & Rodi (1975), developed a hierarchy of Reynolds-stress transport models by consolidating the work of various separate groups into a unified framework. They were able to successfully apply the models to a variety of free-shear and wall-bounded flows of practical interest (Launder et al., 1975). Their model, particularly the simple version (the 'Basic' model), has since been one of the most widely used RST models in engineering applications because of the combined advantage of being simple in form, yet retaining the ability to overcome many of the weaknesses of eddy-viscosity formulations (Hanjali´c & Jakirli´c, 2002).

Later Schumann (1977) introduced the concept of *realisability* as a constraint to guide model formulation. By this it is meant that models should be designed to prevent certain unphysical solutions, such as negative normal stress components, or a stress tensor that violates the Cauchy-Schwartz inequality. Lumley (1978) extensively discussed the significance and implementation of realisability requirements. He devised and used anisotropy invariant maps, or 'Lumley triangles', to illustrate the limiting states of turbulence with respect to values of the second and third invariants of the Reynolds-stress anisotropy tensor. Lumley pointed out that to prevent a negative normal stress component from arising during computations, the time derivative of the component must be made to vanish at the instant when the component itself vanishes, thereby preventing a negative value to arise as time progresses. Such a situation can arise near a wall or a free-surface, where the interface-normal component decays much faster than the other components as the interface is approached, thus approaching a *two-component* limit. Shih & Lumley (1985) later used these arguments to devise a realisable model for the pressure–strain-rate correlation. Their model, however, did not perform well in simple shear flows, and higher order corrections were later added to achieve better agreement with these flows (Craft & Launder, 2002). Speziale (1985; 1987) used arguments of material-frame indifference in the limit of two-dimensional turbulence to develop a model for the rapid pressure–strain-rate correlation. Speziale et al. (1991) later considered the simplest topologically equivalent form (returning the same equilibrium states) to that of the Speziale (1987) model, to arrive at a more simplified, similarly performing version (Speziale et al., 1991). This latter model is also in relatively common use in engineering RST computations. The UMIST group, starting with the work of Fu et al. (1987), Fu (1988), Craft et al. (1989), and Craft (1991) developed a model also based on ensuring realizability in the two-component limit, but using an approach slightly different from that used by Shih & Lumley. This model (the 'TCL' model, in what follows) uses a cubic expansion of the rapid pressure–strain-rate correlation in *k* and *aij*. It was shown to achieve significant improvements over previous

## **3. Modelling practises**

models in a wide range of flows.

In this section the basic equations for the mean flow of incompressible fluids are presented, along with the equations for the relevant turbulence statistics. At the level of Reynolds stress transport modelling, the Reynolds averaged Navier–Stokes (RANS) equations are solved, along with separate equations for each independent component of the Reynolds stress tensor, as well as a transport equation for the scalar rate of dissipation of turbulent kinetic energy. 2 Numerical Simulations

2004). This meant that tensor representation results from the continuum mechanics literature could be used to formulate expressions for the Reynolds-stress tensor, as first proposed by Rivlin (1957). These ideas were then expanded by Crow (1968; 1967), and Lumley (1967; 1970). Computational work accelerated in the 1970's with the works of Daly & Harlow (1970), Reynolds (1970), Donaldson (1971), Naot et al. (1972), Hanjali´c & Launder (1972), and Lumley & Khajeh-Nouri (1974). In a landmark paper, Launder, Reece & Rodi (1975), developed a hierarchy of Reynolds-stress transport models by consolidating the work of various separate groups into a unified framework. They were able to successfully apply the models to a variety of free-shear and wall-bounded flows of practical interest (Launder et al., 1975). Their model, particularly the simple version (the 'Basic' model), has since been one of the most widely used RST models in engineering applications because of the combined advantage of being simple in form, yet retaining the ability to overcome many of the weaknesses of eddy-viscosity

Later Schumann (1977) introduced the concept of *realisability* as a constraint to guide model formulation. By this it is meant that models should be designed to prevent certain unphysical solutions, such as negative normal stress components, or a stress tensor that violates the Cauchy-Schwartz inequality. Lumley (1978) extensively discussed the significance and implementation of realisability requirements. He devised and used anisotropy invariant maps, or 'Lumley triangles', to illustrate the limiting states of turbulence with respect to values of the second and third invariants of the Reynolds-stress anisotropy tensor. Lumley pointed out that to prevent a negative normal stress component from arising during computations, the time derivative of the component must be made to vanish at the instant when the component itself vanishes, thereby preventing a negative value to arise as time progresses. Such a situation can arise near a wall or a free-surface, where the interface-normal component decays much faster than the other components as the interface is approached, thus approaching a *two-component* limit. Shih & Lumley (1985) later used these arguments to devise a realisable model for the pressure–strain-rate correlation. Their model, however, did not perform well in simple shear flows, and higher order corrections were later added to achieve better agreement with these flows (Craft & Launder, 2002). Speziale (1985; 1987) used arguments of material-frame indifference in the limit of two-dimensional turbulence to develop a model for the rapid pressure–strain-rate correlation. Speziale et al. (1991) later considered the simplest topologically equivalent form (returning the same equilibrium states) to that of the Speziale (1987) model, to arrive at a more simplified, similarly performing version (Speziale et al., 1991). This latter model is also in relatively common use in engineering RST computations. The UMIST group, starting with the work of Fu et al. (1987), Fu (1988), Craft et al. (1989), and Craft (1991) developed a model also based on ensuring realizability in the two-component limit, but using an approach slightly different from that used by Shih & Lumley. This model (the 'TCL' model, in what follows) uses a cubic expansion of the rapid pressure–strain-rate correlation in *k* and *aij*. It was shown to achieve significant improvements over previous

In this section the basic equations for the mean flow of incompressible fluids are presented, along with the equations for the relevant turbulence statistics. At the level of Reynolds stress transport modelling, the Reynolds averaged Navier–Stokes (RANS) equations are solved, along with separate equations for each independent component of the Reynolds stress tensor, as well as a transport equation for the scalar rate of dissipation of turbulent kinetic energy.

formulations (Hanjali´c & Jakirli´c, 2002).

models in a wide range of flows.

**3. Modelling practises**

The modelling approach used for the various terms appearing in the exact Reynolds stress transport equation are briefly reviewed.

#### **3.1 Basic equations of turbulent flow**

Turbulent flows are characterised by highly fluctuating velocity, pressure, and other field variables. One approach for dealing with this fluctuating nature of the flow, the one most widely used by engineers, is to work with an averaged form of the basic equations. In *Reynolds averaging* the instantaneous flow variables are decomposed into an average quantity and a fluctuation. Thus,

$$\begin{aligned} \overline{\mathcal{U}}\_{i} &= \mathcal{U}\_{i} + \boldsymbol{u}\_{i} \\ \widetilde{\boldsymbol{P}} &= \boldsymbol{P} + \boldsymbol{p}\_{\prime} \end{aligned} \tag{1}$$

where capital letters denote averaged quantities, and small letters denote purely fluctuating quantities. The averaging can be either over time or over a repeated realisation of an experiment with the same nominal conditions. The latter, *ensemble* averaging, will be implied in the following, to allow for temporal variations of mean quantities. When this decomposition is substituted into the Navier Stokes equations for incompressible flow, and the result is ensemble averaged, one obtains the Reynolds-averaged Navier-Stokes (RANS) equations

$$\frac{\partial \mathcal{U}\_{l}}{\partial t} + \mathcal{U}\_{\dot{j}} \frac{\partial \mathcal{U}\_{l}}{\partial \mathbf{x}\_{\dot{j}}} = -\frac{1}{\rho} \frac{\partial P}{\partial \mathbf{x}\_{i}} + \nu \frac{\partial^{2} \mathcal{U}\_{l}}{\partial \mathbf{x}\_{\dot{j}}^{2}} - \frac{\partial \overline{\mathcal{U}\_{l} \mathcal{U}\_{\dot{j}}}}{\partial \mathbf{x}\_{\dot{j}}}.\tag{2}$$

When the decomposition is substituted into the continuity equation for incompressible flow, and averaging is applied, one obtains for the mean flow

$$\frac{\partial \mathcal{U}\_l}{\partial \mathbf{x}\_l} = \mathbf{0}.\tag{3}$$

If this is subtracted from the instantaneous continuity equation, the continuity condition for the fluctuating velocity is obtained

$$\frac{\partial \mu\_{\bar{l}}}{\partial \mathbf{x}\_{\bar{l}}} = \mathbf{0},\tag{4}$$

meaning that both the mean and fluctuating velocity fields are individually divergence free. The last term in the RANS equation (2) contains the Reynolds stress tensor *uiuj*. Thus the averaging process introduced a new unknown tensor term, and the set of equations is no longer closed. This is called the *closure* problem of averaging approaches. The task of turbulence modelling is to construct appropriate models for these stresses that relate them to the mean flow quantities, and thus to construct a closed set of equations allowing numerical solutions to be obtained. An additional implied objective in the engineering context is for the models to be as computationally inexpensive as possible while being able to reproduce the behaviour and phenomena of relevance to the problem in question, at the required level of accuracy.

A transport equation for the fluctuating velocity can be obtained by subtracting the RANS equation (2) from the Navier-Stokes equation. Using the divergence-free property of the fluctuating field, the result can be written as

$$\frac{\overline{D}u\_i}{\overline{D}t} = -\frac{1}{\rho} \frac{\partial p}{\partial \mathbf{x}\_i} - u\_j \frac{\partial \mathcal{U}\_i}{\partial \mathbf{x}\_j} - \frac{\partial}{\partial \mathbf{x}\_j} \left[ \mu\_i \mu\_j - \overline{\mu\_i \mu\_j} - \nu \frac{\partial u\_i}{\partial \mathbf{x}\_j} \right] \tag{5}$$

An equation for the kinetic energy associated with the turbulent fluctuations, *k* = *uiui*

Reynolds Stress Transport Modelling 7

<sup>−</sup> *uiuiuk* <sup>−</sup> <sup>1</sup>

*<sup>ρ</sup> puiδik*

− *ν ∂ui ∂xk*

<sup>D</sup>*<sup>t</sup>* <sup>=</sup> *<sup>P</sup><sup>κ</sup>* <sup>+</sup> D − *<sup>ε</sup>* (10)

*∂ui ∂xk*

<sup>3</sup> *<sup>δ</sup>ij* . (11)

(*uiuj* − *uiuj*) (12)

<sup>|</sup>**<sup>x</sup>** <sup>−</sup> **<sup>x</sup>**�| <sup>+</sup> Surface integral. (13)

obtained by taking half the contraction of (7). The resulting equation is

 *ν ∂k ∂xk*

D*k*

*aij* <sup>=</sup> *uiuj*

The first term on the right hand side of (9) is the production of turbulent kinetic energy by mean velocity gradients. The next term is the diffusion of turbulent kinetic energy by various mechanisms. Finally, the last term is the scalar dissipation rate of turbulent kinetic energy. The short form (10) defines the notation that will be used in the following for the respective terms in (9). It is often convenient to work with the deviatoric Reynolds stress anisotropy tensor *aij*

*<sup>k</sup>* <sup>−</sup> <sup>2</sup>

Modelling of the pressure–strain rate correlation is to a large extent guided by consideration of the exact equation for it. An equation for the fluctuating pressure can be obtained by taking

> *∂uj ∂xi*

A formal solution to this Poisson equation can be constructed using the method of Green's functions, as first demonstrated by Chou (1945). The Green's function of the Laplacian

) = <sup>−</sup><sup>1</sup>

(*uiuj* − *uiuj*)

It can be seen from this equation that the fluctuating pressure can be decomposed into three components (Pope, 2000), corresponding to the three terms appearing on the right-hand side of (13). The first term is linear in the turbulent fluctuations, and responds directly to changes in mean velocity gradient. It is thus called the *rapid* pressure, *pr*. The second is a turbulence-turbulence interaction term, that does not respond directly to changes in the mean flow, but through the turbulent cascade process, and is thus called the *slow* pressure, *ps*. The last term is the solution to the homogeneous (Laplace) equation and satisfies appropriate boundary conditions that ensure the superposition of the three parts, *p*, satisfies its own boundary conditions (Pope, 2000). This final term is only significant close to a wall or a free surface, and, since the emphasis here is on modelling regions away from walls, it will be

<sup>−</sup> *<sup>∂</sup>*<sup>2</sup> *∂xi∂xj*

4*π*|**x** − **x**�|

**x**� ,*t*

.

d**x**�

*∂Ui ∂xk* + *∂ ∂xk*

the divergence of (5), and invoking continuity (4). This gives

= −2

*g*(**x**|**x**�

*∂Ui ∂xj*

*∂*<sup>2</sup> *p ∂xk∂xk*

> *∂uj ∂x*� *i*

neglected. Wall effects on *φij* are considered in Section 3.3.

<sup>−</sup> *<sup>∂</sup>*<sup>2</sup> *∂x*� *i ∂x*� *j*

1 *ρ*

The fluctuating pressure is thus given by

 −2 *∂Ui ∂x*� *j*

 *V*

D*k*

In short form, this can be written

**3.2 Pressure–strain rate correlation**

defined as,

operator is

*p <sup>ρ</sup>* <sup>=</sup> <sup>1</sup> 4*π* <sup>D</sup>*<sup>t</sup>* <sup>=</sup> <sup>−</sup>*uiuk*

<sup>2</sup> , can be

. (9)

The operator D/*Dt* is used to denote the material derivative following the mean flow

$$\frac{\overline{D}}{\overline{\mathbf{D}}t} = \frac{\partial}{\partial t} + \mathcal{U}\_{\not{p}} \frac{\partial}{\partial \mathbf{x}\_{j}} \,. \tag{6}$$

Since this interpretation will be used exclusively here, the over-bar on this mean-flow material derivative will subsequently be dropped. An exact equation for the Reynolds stresses can be obtained by using (5) to construct

$$\frac{\mathbf{D} \overline{u\_i u\_j}}{\mathbf{D}t} = \overline{u\_j \frac{\mathbf{D} u\_i}{\mathbf{D}t} + u\_i \frac{\mathbf{D} u\_j}{\mathbf{D}t}} \,\prime$$

where it has been assumed that averaging and taking the material derivative (6) commute. The result is

$$\begin{split} \frac{\overline{\partial \overline{u}\_{i} \boldsymbol{u}\_{j}}}{\mathbf{D}t} &= -\left(\overline{\boldsymbol{u}\_{i} \boldsymbol{u}\_{k}} \frac{\partial \boldsymbol{L}\_{j}}{\partial \mathbf{x}\_{k}} + \overline{\boldsymbol{u}\_{j} \boldsymbol{u}\_{k}} \frac{\partial \boldsymbol{L}\_{i}}{\partial \mathbf{x}\_{k}}\right) \\ &+ \overline{\frac{p}{\rho} \left(\frac{\partial \boldsymbol{u}\_{i}}{\partial \mathbf{x}\_{j}} + \frac{\partial \boldsymbol{u}\_{j}}{\partial \mathbf{x}\_{i}}\right)} \\ &+ \frac{\partial}{\partial \mathbf{x}\_{k}} \left[\nu \frac{\overline{\partial \boldsymbol{u}\_{i} \boldsymbol{u}\_{j}}}{\partial \mathbf{x}\_{k}} - \overline{\boldsymbol{u}\_{i} \boldsymbol{u}\_{j} \boldsymbol{u}\_{k}} - \overline{\frac{p}{\rho} \left(\boldsymbol{u}\_{i} \delta\_{jk} + \boldsymbol{u}\_{j} \delta\_{ik}\right)}\right] \\ &- 2\nu \frac{\overline{\partial \boldsymbol{u}\_{i}}}{\overline{\partial \mathbf{x}\_{k}}} \frac{\partial \boldsymbol{u}\_{j}}{\partial \mathbf{x}\_{k}}. \end{split} \tag{7}$$

The first term on the right hand side above is the production rate of Reynolds stresses by mean velocity gradients. This term is closed at the RST level since it is given in terms of quantities that are being solved for at this level. All the remaining terms in the equation, except for viscous diffusion, require modelling. The second term is a correlation between the fluctuating pressure and the fluctuating strain rate. From continuity this term is traceless, so it does not contribute directly to the kinetic energy of the turbulence. Its effect is to redistribute the energy between the stress components, so it plays a very important role in determining the degree of anisotropy of the stresses. Accordingly, it has received much attention from researchers, and continues to do so. The third term in (7) is a combination of several diffusion terms, all having the effect of spatial redistribution of the Reynolds stresses. Finally, the last term is the dissipation rate of Reynolds stresses by viscous action at the smallest scales of turbulence. Since the smallest scales of motion are assumed to be isotropic, the dissipation rate tensor is frequently modelled as *εij* = <sup>2</sup> <sup>3</sup> *ε δij*, where *ε* is the scalar dissipation rate of turbulent kinetic energy. This approximation is not applicable near walls or free surfaces, where the dissipation tensor becomes markedly anisotropic. Equation (7) can be written in short form as

$$\frac{\mathbf{D}\overline{u\_{i}u\_{j}}}{\mathbf{D}t} = P\_{ij} + \phi\_{ij} + \mathcal{D}\_{ij} - \varepsilon\_{ij} \tag{8}$$

where it is understood that each term above defines the notation for the corresponding term in (7).

4 Numerical Simulations

Since this interpretation will be used exclusively here, the over-bar on this mean-flow material derivative will subsequently be dropped. An exact equation for the Reynolds stresses can be

> D*ui* <sup>D</sup>*<sup>t</sup>* <sup>+</sup> *ui*

where it has been assumed that averaging and taking the material derivative (6) commute.

+ *ujuk*

<sup>−</sup> *uiujuk* <sup>−</sup> *<sup>p</sup>*

The first term on the right hand side above is the production rate of Reynolds stresses by mean velocity gradients. This term is closed at the RST level since it is given in terms of quantities that are being solved for at this level. All the remaining terms in the equation, except for viscous diffusion, require modelling. The second term is a correlation between the fluctuating pressure and the fluctuating strain rate. From continuity this term is traceless, so it does not contribute directly to the kinetic energy of the turbulence. Its effect is to redistribute the energy between the stress components, so it plays a very important role in determining the degree of anisotropy of the stresses. Accordingly, it has received much attention from researchers, and continues to do so. The third term in (7) is a combination of several diffusion terms, all having the effect of spatial redistribution of the Reynolds stresses. Finally, the last term is the dissipation rate of Reynolds stresses by viscous action at the smallest scales of turbulence. Since the smallest scales of motion are assumed to be isotropic, the dissipation rate tensor is

energy. This approximation is not applicable near walls or free surfaces, where the dissipation

where it is understood that each term above defines the notation for the corresponding term

tensor becomes markedly anisotropic. Equation (7) can be written in short form as

D*uiuj*

*∂Ui ∂xk*

*ρ*

(*uiδjk* + *ujδik*)

<sup>3</sup> *ε δij*, where *ε* is the scalar dissipation rate of turbulent kinetic

<sup>D</sup>*<sup>t</sup>* <sup>=</sup> *Pij* <sup>+</sup> *<sup>φ</sup>ij* <sup>+</sup> <sup>D</sup>*ij* <sup>−</sup> *<sup>ε</sup>ij* , (8)

(7)

D*uj* <sup>D</sup>*<sup>t</sup>* ,

. (6)

The operator D/*Dt* is used to denote the material derivative following the mean flow

<sup>=</sup> *<sup>∂</sup> ∂t* + *Uj ∂ ∂xj*

*D* D*t*

D*uiuj* <sup>D</sup>*<sup>t</sup>* <sup>=</sup> *uj*

 *uiuk ∂Uj ∂xk*

 *∂ui ∂xj* + *∂uj ∂xi*

> *ν ∂uiuj ∂xk*

> > *∂uj ∂xk* .

+ *p ρ*

+ *∂ ∂xk*

−2*ν ∂ui ∂xk*

obtained by using (5) to construct

frequently modelled as *εij* = <sup>2</sup>

in (7).

D*uiuj* <sup>D</sup>*<sup>t</sup>* <sup>=</sup> <sup>−</sup>

The result is

An equation for the kinetic energy associated with the turbulent fluctuations, *k* = *uiui* <sup>2</sup> , can be obtained by taking half the contraction of (7). The resulting equation is

$$\frac{\mathrm{D}k}{\mathrm{D}t} = -\overline{u\_{i}u\_{k}}\frac{\partial U\_{i}}{\partial \mathbf{x}\_{k}} + \frac{\partial}{\partial \mathbf{x}\_{k}} \left[ \nu \frac{\partial \mathbf{k}}{\partial \mathbf{x}\_{k}} - \overline{u\_{i}u\_{i}u\_{k}} - \frac{1}{\rho} \overline{p u\_{i}} \delta\_{ik} \right] - \nu \frac{\overline{\partial u\_{i}}}{\partial \mathbf{x}\_{k}} \frac{\partial \overline{u\_{i}}}{\partial \mathbf{x}\_{k}} \,. \tag{9}$$

In short form, this can be written

$$\frac{\text{D}k}{\text{D}t} = P\_{\text{k}} + \mathcal{D} - \varepsilon \tag{10}$$

The first term on the right hand side of (9) is the production of turbulent kinetic energy by mean velocity gradients. The next term is the diffusion of turbulent kinetic energy by various mechanisms. Finally, the last term is the scalar dissipation rate of turbulent kinetic energy. The short form (10) defines the notation that will be used in the following for the respective terms in (9). It is often convenient to work with the deviatoric Reynolds stress anisotropy tensor *aij* defined as,

$$a\_{ij} = \frac{\overline{\mu\_i \mu\_j}}{k} - \frac{2}{3} \delta\_{ij} \,. \tag{11}$$

#### **3.2 Pressure–strain rate correlation**

Modelling of the pressure–strain rate correlation is to a large extent guided by consideration of the exact equation for it. An equation for the fluctuating pressure can be obtained by taking the divergence of (5), and invoking continuity (4). This gives

$$\frac{1}{\rho} \frac{\partial^2 p}{\partial \mathbf{x}\_k \partial \mathbf{x}\_k} = -2 \frac{\partial \mathcal{U}\_i}{\partial \mathbf{x}\_j} \frac{\partial \boldsymbol{u}\_j}{\partial \mathbf{x}\_i} - \frac{\partial^2}{\partial \mathbf{x}\_i \partial \mathbf{x}\_j} (\boldsymbol{u}\_i \boldsymbol{u}\_j - \overline{\boldsymbol{u}\_i} \overline{\boldsymbol{u}\_j}) \tag{12}$$

A formal solution to this Poisson equation can be constructed using the method of Green's functions, as first demonstrated by Chou (1945). The Green's function of the Laplacian operator is

$$g(\mathbf{x}|\mathbf{x}') = \frac{-1}{4\pi|\mathbf{x} - \mathbf{x}'|} \dots$$

The fluctuating pressure is thus given by

$$\frac{p}{\rho} = \frac{1}{4\pi} \iiint\_{V} \left[ -2 \frac{\partial L\_{l}}{\partial \mathbf{x}\_{j}^{\prime}} \frac{\partial u\_{j}}{\partial \mathbf{x}\_{i}^{\prime}} - \frac{\partial^{2}}{\partial \mathbf{x}\_{i}^{\prime} \partial \mathbf{x}\_{j}^{\prime}} (u\_{i} u\_{j} - \overline{u\_{i} u\_{j}}) \right]\_{\mathbf{x}^{\prime}, t} \frac{d\mathbf{x}^{\prime}}{|\mathbf{x} - \mathbf{x}^{\prime}|} + \text{Surface integral.} \tag{13}$$

It can be seen from this equation that the fluctuating pressure can be decomposed into three components (Pope, 2000), corresponding to the three terms appearing on the right-hand side of (13). The first term is linear in the turbulent fluctuations, and responds directly to changes in mean velocity gradient. It is thus called the *rapid* pressure, *pr*. The second is a turbulence-turbulence interaction term, that does not respond directly to changes in the mean flow, but through the turbulent cascade process, and is thus called the *slow* pressure, *ps*. The last term is the solution to the homogeneous (Laplace) equation and satisfies appropriate boundary conditions that ensure the superposition of the three parts, *p*, satisfies its own boundary conditions (Pope, 2000). This final term is only significant close to a wall or a free surface, and, since the emphasis here is on modelling regions away from walls, it will be neglected. Wall effects on *φij* are considered in Section 3.3.

where *Dij* is given by

expressed in terms of *k* and *aij*

As for the slow pressure–strain-rate term, *φ<sup>s</sup>*

strain term have been suggested in the literature.

**3.3 Wall effects on** *φij*

*Dij* = −*uiuk*

et al. (1972), is sometimes termed the *isotropization of production* model (LRR-IP),

*ij* <sup>=</sup> <sup>−</sup>*c*2(*Pij* <sup>−</sup> <sup>2</sup>

Various other models have been proposed following similar lines of reasoning, in which M is modelled as a tensor-polynomial function of the Reynolds stress tensor or, equivalently,

It is worth pointing out at this stage that there is an intrinsic weakness in all such models of the form (25). The tensor M, as defined by (17), contains two kinds of directional information – the direction of the energetic velocity components, and the direction of variation or dependence of the two-point correlation (Pope, 2000). Only the former type of information is contained in the Reynolds stress tensor, so two fields having the same Reynolds stresses can have different M tensors. More explicitly put, the evolution of the Reynolds stresses is not uniquely determined by the Reynolds stresses (Pope, 2000). This is an intrinsic limitation in RST modelling, that is difficult to overcome without significantly complicating the modelling approach and/or computational cost (Johansson & Hallbäck, 1994; Kassinos & Reynolds, 1994). This limitation is known to cause poor results in flows where the velocity gradient has a strong rotational component, such as in pure (or dominant) rotation, and in high shear rate flows (Johansson & Hallbäck, 1994). However, in many other flows, including ones with significant rotational effects, RST models have been shown to produce very good results.

expression, pertaining to the non-linear turbulence–turbulence interaction part of (13). Most early models followed Rotta's (1951) linear return to isotropy model for the slow term

This model is motivated by the decay of homogeneous anisotropic turbulence in the absence of mean velocity gradients. It is generally observed that in such cases turbulence progressively

Experimental evidence shows that the return-to-isotropy process is in fact non-linear in *aij* (Chung & Kim, 1995). When plotted on anisotropy invariant maps, the paths taken during return-to-isotropy experiments are not straight lines, and have different behaviour depending on the sign of the third invariant (Pope, 2000). It is also found that the rate of return is highly dependent on the Reynolds number. A number of nonlinear models for the slow pressure

The presence of a wall alters pressure fluctuations by viscous effect through the no-slip condition, and by inviscid effect through the impermeability condition. DNS results show that

*φs*

tends towards an isotropic state, hence the negative sign in (26).

*φr*

*∂Uk ∂xj*

Reynolds Stress Transport Modelling 9

and *D* = *Dii*/2. This is the first of the two Launder-Reece-Rodi (LRR) models in Launder et al. (1975), called the *Quasi-Isotropic* model (LRR-QI). A simplified version of (22) was also suggested in Launder et al. (1975) by observing that the dominant term in this equation is the first one appearing on the right hand side. The model thus obtained, first proposed by Naot

− *ujuk*

*∂Uk ∂xi*

, (23)

<sup>3</sup> *P δij*). (24)

M = M(*k*, **a**) . (25)

*ij*, it is difficult to extract anything from the exact

*ij* = −*C*1*εaij* . (26)

Based on the above decomposition, the pressure–strain rate correlation can similarly be decomposed into rapid, slow, and *wall influence* terms. The rapid part can be constructed as follows

$$\phi\_{ij}^r = \frac{p^r}{\rho} \left( \frac{\partial u\_i}{\partial \mathbf{x}\_j} + \frac{\partial u\_j}{\partial \mathbf{x}\_i} \right) \tag{14}$$

$$\begin{split} \frac{p^{r}}{\rho} \left( \frac{\partial u\_{i}}{\partial \mathbf{x}\_{j}} \right) &= \frac{1}{2\pi} \iiint\_{-\infty}^{\infty} \left( \frac{\partial \mathcal{U}\_{k}}{\partial \mathbf{x}\_{l}^{\prime}} \frac{\partial u\_{l}}{\partial \mathbf{x}\_{k}^{\prime}} \right)\_{\mathbf{x}^{\prime},t} \left( \frac{\partial u\_{i}}{\partial \mathbf{x}\_{j}} \right)\_{\mathbf{x},t} \frac{\mathbf{dx}^{\prime}}{|\mathbf{x} - \mathbf{x}^{\prime}|} \\ &= \frac{1}{2\pi} \frac{\partial \mathcal{U}\_{k}}{\partial \mathbf{x}\_{l}} \iiint\_{-\infty}^{\infty} \frac{\partial^{2}}{\partial \mathbf{x}\_{j} \partial \mathbf{x}\_{k}^{\prime}} (u\_{i}(\mathbf{x})u\_{l}(\mathbf{x}^{\prime})) \frac{\mathbf{dx}^{\prime}}{|\mathbf{x} - \mathbf{x}^{\prime}|}. \end{split} \tag{15}$$

In taking *<sup>∂</sup>Uk ∂x*� *l* outside the integral it is assumed that this term is reasonably constant over the volume integral. In homogeneous flows, that is of course exact, but is an approximation in inhomogeneous ones. One can thus write:

$$\boldsymbol{\phi}^{r} = \frac{\partial \mathcal{U}\_{k}}{\partial \boldsymbol{\chi}\_{l}} \left( \mathcal{M}\_{iljk} + \mathcal{M}\_{jllk} \right) \,. \tag{16}$$

where the fourth rank tensor M*iljk* is given by

$$\mathcal{M}\_{\text{iljk}} = \frac{-1}{2\pi} \iiint\_{\infty}^{\infty} \frac{\partial^2 \overline{u\_i(\mathbf{x}) u\_l(\mathbf{x} + \mathbf{r})}}{\partial r\_j \partial r\_k} \frac{\mathbf{dr}}{|\mathbf{r}|},\tag{17}$$

using **r** = **x**� − **x** for the separation distance. The M*iljk* tensor is symmetric in the first two indices, and in the last two

$$\mathcal{M}\_{\rm iljk} = \mathcal{M}\_{\rm lijk} = \mathcal{M}\_{\rm illkj}.\tag{18}$$

The divergence-free velocity condition means that contraction over the middle indices results in the quantity vanishing:

$$\mathcal{M}\_{ijjk} = 0 \,,\tag{19}$$

and contraction over the last two indices can be shown to yield (twice) the Reynolds stress tensor

$$\mathcal{M}\_{\text{ilk}} = 2 \,\overbrace{\mathcal{U}\_i \mathcal{U}\_l} \,. \tag{20}$$

The last of these kinematic conditions (20) suggested to workers that the M tensor could be modelled as a function of the Reynolds stresses (Launder et al., 1975). The approach taken was to model M as a polynomial function in the stresses. The most general fourth-rank tensor linear in the Reynolds stresses satisfying the symmetry conditions (18) is

$$\begin{split} \mathcal{M}\_{\text{ijkl}} &= \mathbf{a}\delta\_{\text{kl}}\overline{\mathbf{u}\_{\text{i}}}\overline{\mathbf{u}\_{\text{j}}} + \beta(\delta\_{i\mathbf{k}}\overline{\mathbf{u}\_{\text{j}}}\overline{\mathbf{u}\_{\text{l}}} + \delta\_{i\mathbf{l}}\overline{\mathbf{u}\_{\text{j}}}\overline{\mathbf{u}\_{\text{k}}} + \delta\_{j\mathbf{k}}\overline{\mathbf{u}\_{\text{i}}}\overline{\mathbf{u}\_{\text{l}}} + \delta\_{j\mathbf{l}}\overline{\mathbf{u}\_{\text{i}}}\overline{\mathbf{u}\_{\text{k}}}) \\ &+ \gamma\delta\_{i\mathbf{j}}\overline{\mathbf{u}\_{\text{k}}}\overline{\mathbf{u}\_{\text{l}}} + [\eta\delta\_{i\mathbf{j}}\delta\_{\text{kl}} + \upsilon(\delta\_{i\mathbf{k}}\delta\_{\mathbf{j}\mathbf{l}} + \delta\_{i\mathbf{l}}\delta\_{\mathbf{j}\mathbf{k}})] \mathbf{k} . \end{split} \tag{21}$$

where the coefficients *α*, *β*, *γ*, *η*, *v* are constants (or functions of the invariants of *aij*). The continuity condition (19), and the normalisation condition (20) can be used to reduce the number of undetermined constants to one. When this is done, and the resulting modelled M*ijkl* is substituted into (16) the resulting linear rapid pressure–strain rate model is

$$\phi\_{ij}^r = -\frac{\gamma+8}{11}(P\_{ij} - 2/3P\delta\_{ij}) - \frac{30\gamma-2}{55}k\left(\frac{\partial \mathcal{U}\_i}{\partial \mathbf{x}\_j} + \frac{\partial \mathcal{U}\_j}{\partial \mathbf{x}\_i}\right) - \frac{(8\gamma-2)}{11}(D\_{ij} - 2/3D\delta\_{ij}), \quad \text{(22)}$$

where *Dij* is given by

6 Numerical Simulations

Based on the above decomposition, the pressure–strain rate correlation can similarly be decomposed into rapid, slow, and *wall influence* terms. The rapid part can be constructed

> *∂ul ∂x*� *k*

*∂*2 *∂xj∂x*� *k*

volume integral. In homogeneous flows, that is of course exact, but is an approximation in

using **r** = **x**� − **x** for the separation distance. The M*iljk* tensor is symmetric in the first two

The divergence-free velocity condition means that contraction over the middle indices results

and contraction over the last two indices can be shown to yield (twice) the Reynolds stress

The last of these kinematic conditions (20) suggested to workers that the M tensor could be modelled as a function of the Reynolds stresses (Launder et al., 1975). The approach taken was to model M as a polynomial function in the stresses. The most general fourth-rank tensor

M*ijkl* =*αδkluiuj* + *β*(*δikujul* + *δilujuk* + *δjkuiul* + *δjluiuk*)

where the coefficients *α*, *β*, *γ*, *η*, *v* are constants (or functions of the invariants of *aij*). The continuity condition (19), and the normalisation condition (20) can be used to reduce the number of undetermined constants to one. When this is done, and the resulting modelled

M*ijkl* is substituted into (16) the resulting linear rapid pressure–strain rate model is

 *∂Ui ∂xj* + *∂Uj ∂xi*

<sup>55</sup> *<sup>k</sup>*

**x**� ,*t*

outside the integral it is assumed that this term is reasonably constant over the

M*iljk* + M*jlik*

*∂*<sup>2</sup>*ui*(**x**)*ul*(**x** + **r**) *∂rj∂rk*

 *∂ui ∂xj*

(*ui*(**x**)*ul*(**x**�

**x**,*t*

d**r** |**r**|

M*iljk* = M*lijk* = M*ilkj*. (18)

M*ijjk* = 0 , (19)

M*ilkk* = 2 *uiul* . (20)

<sup>+</sup>*γδijukul* + [*ηδijδkl* <sup>+</sup> *<sup>v</sup>*(*δikδjl* <sup>+</sup> *<sup>δ</sup>ilδjk*)]*k*, (21)

<sup>−</sup> (8*<sup>γ</sup>* <sup>−</sup> <sup>2</sup>)

<sup>11</sup> (*Dij* <sup>−</sup> 2/3*Dδij*), (22)

d**x**� |**x** − **x**�|

)) <sup>d</sup>**x**� |**x** − **x**�| .

, (16)

, (17)

(14)

(15)

 *∂ui ∂xj* + *∂uj ∂xi*

 *∂Uk ∂x*� *l*

 ∞ −∞

*φr ij* <sup>=</sup> *<sup>p</sup><sup>r</sup> ρ*

 ∞ −∞

*∂Uk ∂xl*

*<sup>φ</sup><sup>r</sup>* <sup>=</sup> *<sup>∂</sup>Uk ∂xl* 

2*π*

linear in the Reynolds stresses satisfying the symmetry conditions (18) is

<sup>11</sup> (*Pij* <sup>−</sup> 2/3*Pδij*) <sup>−</sup> <sup>30</sup>*<sup>γ</sup>* <sup>−</sup> <sup>2</sup>

 ∞ ∞

<sup>M</sup>*iljk* <sup>=</sup> <sup>−</sup><sup>1</sup>

as follows

In taking *<sup>∂</sup>Uk*

*∂x*� *l*

indices, and in the last two

in the quantity vanishing:

tensor

*φr*

*ij* <sup>=</sup> <sup>−</sup>*<sup>γ</sup>* <sup>+</sup> <sup>8</sup>

*pr ρ*

 *∂ui ∂xj*

inhomogeneous ones. One can thus write:

where the fourth rank tensor M*iljk* is given by

 <sup>=</sup> <sup>1</sup> 2*π*

> <sup>=</sup> <sup>1</sup> 2*π*

$$D\_{ij} = -\overline{\mu\_i \mu\_k} \frac{\partial \mathcal{U}\_k}{\partial \mathbf{x}\_j} - \overline{\mu\_j \mu\_k} \frac{\partial \mathcal{U}\_k}{\partial \mathbf{x}\_i} \,. \tag{23}$$

and *D* = *Dii*/2. This is the first of the two Launder-Reece-Rodi (LRR) models in Launder et al. (1975), called the *Quasi-Isotropic* model (LRR-QI). A simplified version of (22) was also suggested in Launder et al. (1975) by observing that the dominant term in this equation is the first one appearing on the right hand side. The model thus obtained, first proposed by Naot et al. (1972), is sometimes termed the *isotropization of production* model (LRR-IP),

$$
\phi\_{ij}^r = -\mathbf{c}\_2(P\_{i\bar{j}} - \frac{2}{3}P\_i\delta\_{i\bar{j}}).\tag{24}
$$

Various other models have been proposed following similar lines of reasoning, in which M is modelled as a tensor-polynomial function of the Reynolds stress tensor or, equivalently, expressed in terms of *k* and *aij*

$$\mathcal{M} = \mathcal{M}(k, \mathbf{a}) \,. \tag{25}$$

It is worth pointing out at this stage that there is an intrinsic weakness in all such models of the form (25). The tensor M, as defined by (17), contains two kinds of directional information – the direction of the energetic velocity components, and the direction of variation or dependence of the two-point correlation (Pope, 2000). Only the former type of information is contained in the Reynolds stress tensor, so two fields having the same Reynolds stresses can have different M tensors. More explicitly put, the evolution of the Reynolds stresses is not uniquely determined by the Reynolds stresses (Pope, 2000). This is an intrinsic limitation in RST modelling, that is difficult to overcome without significantly complicating the modelling approach and/or computational cost (Johansson & Hallbäck, 1994; Kassinos & Reynolds, 1994). This limitation is known to cause poor results in flows where the velocity gradient has a strong rotational component, such as in pure (or dominant) rotation, and in high shear rate flows (Johansson & Hallbäck, 1994). However, in many other flows, including ones with significant rotational effects, RST models have been shown to produce very good results.

As for the slow pressure–strain-rate term, *φ<sup>s</sup> ij*, it is difficult to extract anything from the exact expression, pertaining to the non-linear turbulence–turbulence interaction part of (13). Most early models followed Rotta's (1951) linear return to isotropy model for the slow term

$$
\phi\_{ij}^s = -\mathsf{C}\_1 \varepsilon a\_{ij} \,. \tag{26}
$$

This model is motivated by the decay of homogeneous anisotropic turbulence in the absence of mean velocity gradients. It is generally observed that in such cases turbulence progressively tends towards an isotropic state, hence the negative sign in (26).

Experimental evidence shows that the return-to-isotropy process is in fact non-linear in *aij* (Chung & Kim, 1995). When plotted on anisotropy invariant maps, the paths taken during return-to-isotropy experiments are not straight lines, and have different behaviour depending on the sign of the third invariant (Pope, 2000). It is also found that the rate of return is highly dependent on the Reynolds number. A number of nonlinear models for the slow pressure strain term have been suggested in the literature.

### **3.3 Wall effects on** *φij*

The presence of a wall alters pressure fluctuations by viscous effect through the no-slip condition, and by inviscid effect through the impermeability condition. DNS results show that

of energy at the large scales. This assumption is an obvious weakness in the model when the turbulence is not in equilibrium, as when unsteady solutions are sought, or where the time-scale of the mean flow is of the same order or smaller than the characteristic time-scale of turbulence. In such cases the small-scale turbulence may not have enough time to adjust to the large-scale scale variations, and the instantaneous link implied by the production term in

Reynolds Stress Transport Modelling 11

The destruction term in (29) is motivated by consideration of the decay of homogeneous isotropic turbulence in the absence of production (Pope, 2000). In such a flow one expects that the turbulence will decay in a self-similar form in which the rates of decay of *k* and *ε* are

> <sup>=</sup> <sup>−</sup>*k*/*<sup>ε</sup> ε*/ *<sup>d</sup><sup>ε</sup> dt*

If this proportionality constant is labelled *Cε*2, the following destruction term is implied

*dt* <sup>=</sup> <sup>−</sup>*Cε*<sup>2</sup>

There are three diffusive transport terms on the right hand side of (7). The first is the viscous

which is closed and does not require modelling. The following two terms are the pressure diffusion and turbulent convection, respectively. Most commonly these are modelled together as a combined turbulent diffusion term, *Tij*, using the generalised gradient diffusion

*∂*<sup>2</sup>*uiuj ∂xk∂xk*

> *∂uiuj ∂xk*

*∂ukui ∂xl*

+ *ukul*

*∂uiuj ∂xl*

. (33)

*Vij* = *ν*

 *Cs ε k uluk*

*∂ujuk ∂xl*

A deficiency of this model is that it does not preserve the symmetry under cyclic permutation of indices that is exhibited by the triple velocity moments *uiujuk*. This is only significant when the triple moments and pressure diffusion are modelled separately. In such case an improved

+ *ujul*

More elaborate models exist in the literature, as in (Craft, 1998) for example, but the models

*Tij* <sup>=</sup> *<sup>∂</sup> ∂xl*

model that has been suggested by Hanjali´c & Launder (1972) is often used,

*k ε uiul* *dε*

= *C*

*ε*2

*<sup>k</sup>* (30)

(31)

(32)

*k*/ *dk dt ε*/ *<sup>d</sup><sup>ε</sup> dt*

(29) is questionable.

**3.5 Diffusion modelling**

where *Cs* is typically 0.22.

hypothesis (GGDH) of Daly & Harlow (1970),

*uiujuk* = −*Cs*

mentioned above are the ones more commonly used.

diffusion term

proportional

the viscous effect is confined to a region within *<sup>y</sup>*<sup>+</sup> <sup>≈</sup> 15 from the wall (Mansour et al., 1988). The inviscid wall-blocking effect on the other hand is significant where the distance from the wall is of the same order as the turbulent length scale. Wall blocking causes two opposing effects; *wall reflection* of the fluctuating pressure field increases the energy-redistributing pressure fluctuations, which pushes turbulence towards isotropy, while it also causes selective damping of the wall-normal fluctuating velocity component in turbulent eddies, thereby increasing anisotropy. The latter effect dominates, and turbulence anisotropy near a wall is higher than that in a free shear flow at a similar rate of shear. To account for this, Gibson & Launder (1978) proposed two additive corrections to *φij* using the unit normal vector to the wall, *ni*. The first, based on the proposal of Shir (1973), is an additive correction to the slow part

$$\boldsymbol{\phi}\_{ij}^{\mathbf{s},\mathbf{w}} = \mathbf{C}\_1^{\mathbf{w}} \frac{\varepsilon}{k} \left( \overline{\boldsymbol{\mu}\_k \boldsymbol{\upmu}\_m} \boldsymbol{n}\_k \boldsymbol{n}\_m \boldsymbol{\updelta}\_{ij} - \frac{3}{2} \overline{\boldsymbol{\upmu}\_i \boldsymbol{\upmu}\_k} \boldsymbol{n}\_k \boldsymbol{n}\_j - \frac{3}{2} \overline{\boldsymbol{\upmu}\_j \boldsymbol{\upmu}\_k} \boldsymbol{n}\_k \boldsymbol{n}\_i \right) \boldsymbol{f}\_{\mathbf{w}} \tag{27}$$

and the second, is a correction to the rapid part

$$\boldsymbol{\Phi}^{r,w}\_{ij} = \mathbb{C}^w\_2 \left( \boldsymbol{\phi}^r\_{km} \boldsymbol{n}\_k \boldsymbol{n}\_m \boldsymbol{\delta}\_{ij} - \frac{3}{2} \boldsymbol{\phi}^r\_{ik} \boldsymbol{n}\_k \boldsymbol{n}\_j - \frac{3}{2} \boldsymbol{\phi}^r\_{jk} \boldsymbol{n}\_k \boldsymbol{n}\_i \right) f\_{iw} \tag{28}$$

where *C<sup>w</sup>* <sup>1</sup> <sup>=</sup> 0.5, *<sup>C</sup><sup>w</sup>* <sup>2</sup> <sup>=</sup> 0.3, and *fw* <sup>=</sup> 0.4*k*3/2/(*εxn*) is a damping function based on the ratio of the turbulence length scale to the normal distance to the wall, *xn*.

#### **3.4 Modelling dissipation**

While modelling of the turbulent kinetic energy, and of the pressure–strain rate correlation, has been to at least some degree guided by consideration of their exact equations, the same is not true for the standard dissipation rate model (Pope, 2000). Dissipation of turbulent kinetic energy is associated with the smallest scales of the fluctuating field, while the kinetic energy itself is mostly contained in the largest scales of fluctuations. The exact dissipation rate equation is comprised of a large number of terms that are all related to dissipative-scale processes, and all but one of the source-terms require modelling. It is thus not a useful starting point for modelling the dissipation rate. Instead the more empirical approach taken is motivated by the spectral energy transfer view of dissipation. The kinetic energy of the larger energy containing eddies is transferred by vortex-stretching in the presence of mean velocity gradients to smaller eddies, and the same process occurs at the 'next' smaller scales, and so on to the smallest dissipative scales, where kinetic energy is finally converted to heat by viscous (molecular) action. If the molecular viscosity is somehow changed, all that happens is that the size of the dissipative scales change to accommodate the rate of energy they receive, but the rate itself is not affected. Thus even though the *mechanism* of dissipation is governed by processes that occur at the smallest scales, dissipation can also be viewed as an energy-transfer rate that readjusts itself with the amount of energy it receives. In this sense, the amount (as opposed to the mechanism) of dissipation is in fact determined by the energy in larger scales. Under the assumption of spectral equilibrium, the transfer rate of energy across the spectrum of turbulence scales is constant and determined by the rate of energy input. Based on this assumption, and the preceding arguments, the conventional equation for dissipation is assumed to be of the form

$$\frac{\text{D}\varepsilon}{\text{D}t} = \text{C}\_{\varepsilon1}\frac{\varepsilon}{k}P + \text{D}\_{\varepsilon} - \text{C}\_{\varepsilon2}\frac{\varepsilon^2}{k} \,' \,. \tag{29}$$

where D*<sup>ε</sup>* is the diffusion of *ε*. The modelled production term above reflects the assumed direct link between a single rate of transfer of energy across the spectrum and production 8 Numerical Simulations

the viscous effect is confined to a region within *<sup>y</sup>*<sup>+</sup> <sup>≈</sup> 15 from the wall (Mansour et al., 1988). The inviscid wall-blocking effect on the other hand is significant where the distance from the wall is of the same order as the turbulent length scale. Wall blocking causes two opposing effects; *wall reflection* of the fluctuating pressure field increases the energy-redistributing pressure fluctuations, which pushes turbulence towards isotropy, while it also causes selective damping of the wall-normal fluctuating velocity component in turbulent eddies, thereby increasing anisotropy. The latter effect dominates, and turbulence anisotropy near a wall is higher than that in a free shear flow at a similar rate of shear. To account for this, Gibson & Launder (1978) proposed two additive corrections to *φij* using the unit normal vector to the wall, *ni*. The first, based on the proposal of Shir (1973), is an additive correction to the slow

*ukumnknmδij* <sup>−</sup> <sup>3</sup>

*kmnknmδij* <sup>−</sup> <sup>3</sup>

of the turbulence length scale to the normal distance to the wall, *xn*.

D*ε* <sup>D</sup>*<sup>t</sup>* <sup>=</sup> *<sup>C</sup>ε*<sup>1</sup> *ε k*

where D*<sup>ε</sup>* is the diffusion of *ε*. The modelled production term above reflects the assumed direct link between a single rate of transfer of energy across the spectrum and production

*P* + D*<sup>ε</sup>* − *Cε*<sup>2</sup>

*ε*2

*<sup>k</sup>* , (29)

<sup>2</sup> *uiuknknj* <sup>−</sup> <sup>3</sup>

*iknknj* <sup>−</sup> <sup>3</sup>

<sup>2</sup> <sup>=</sup> 0.3, and *fw* <sup>=</sup> 0.4*k*3/2/(*εxn*) is a damping function based on the ratio

2*φr jknkni* 

2*φr*

While modelling of the turbulent kinetic energy, and of the pressure–strain rate correlation, has been to at least some degree guided by consideration of their exact equations, the same is not true for the standard dissipation rate model (Pope, 2000). Dissipation of turbulent kinetic energy is associated with the smallest scales of the fluctuating field, while the kinetic energy itself is mostly contained in the largest scales of fluctuations. The exact dissipation rate equation is comprised of a large number of terms that are all related to dissipative-scale processes, and all but one of the source-terms require modelling. It is thus not a useful starting point for modelling the dissipation rate. Instead the more empirical approach taken is motivated by the spectral energy transfer view of dissipation. The kinetic energy of the larger energy containing eddies is transferred by vortex-stretching in the presence of mean velocity gradients to smaller eddies, and the same process occurs at the 'next' smaller scales, and so on to the smallest dissipative scales, where kinetic energy is finally converted to heat by viscous (molecular) action. If the molecular viscosity is somehow changed, all that happens is that the size of the dissipative scales change to accommodate the rate of energy they receive, but the rate itself is not affected. Thus even though the *mechanism* of dissipation is governed by processes that occur at the smallest scales, dissipation can also be viewed as an energy-transfer rate that readjusts itself with the amount of energy it receives. In this sense, the amount (as opposed to the mechanism) of dissipation is in fact determined by the energy in larger scales. Under the assumption of spectral equilibrium, the transfer rate of energy across the spectrum of turbulence scales is constant and determined by the rate of energy input. Based on this assumption, and the preceding arguments, the conventional equation for dissipation

<sup>2</sup> *ujuknkni*

*fw* (27)

*fw* (28)

part

where *C<sup>w</sup>*

*φs*,*<sup>w</sup> ij* <sup>=</sup> *<sup>C</sup><sup>w</sup>* 1 *ε k* 

<sup>1</sup> <sup>=</sup> 0.5, *<sup>C</sup><sup>w</sup>*

**3.4 Modelling dissipation**

is assumed to be of the form

and the second, is a correction to the rapid part

*φr*,*<sup>w</sup> ij* <sup>=</sup> *<sup>C</sup><sup>w</sup>* 2 *φr* of energy at the large scales. This assumption is an obvious weakness in the model when the turbulence is not in equilibrium, as when unsteady solutions are sought, or where the time-scale of the mean flow is of the same order or smaller than the characteristic time-scale of turbulence. In such cases the small-scale turbulence may not have enough time to adjust to the large-scale scale variations, and the instantaneous link implied by the production term in (29) is questionable.

The destruction term in (29) is motivated by consideration of the decay of homogeneous isotropic turbulence in the absence of production (Pope, 2000). In such a flow one expects that the turbulence will decay in a self-similar form in which the rates of decay of *k* and *ε* are proportional

$$\frac{k/\frac{dk}{dt}}{\varepsilon/\frac{d\varepsilon}{dt}} = \frac{-k/\varepsilon}{\varepsilon/\frac{d\varepsilon}{dt}} = \mathcal{C}$$

If this proportionality constant is labelled *Cε*2, the following destruction term is implied

$$\frac{d\varepsilon}{dt} = -\mathbb{C}\_{\varepsilon 2} \frac{\varepsilon^2}{k} \tag{30}$$

#### **3.5 Diffusion modelling**

There are three diffusive transport terms on the right hand side of (7). The first is the viscous diffusion term

$$\mathcal{V}\_{ij} = \nu \frac{\partial^2 \overline{u\_i u\_j}}{\partial \mathbf{x}\_k \partial \mathbf{x}\_k} \tag{31}$$

which is closed and does not require modelling. The following two terms are the pressure diffusion and turbulent convection, respectively. Most commonly these are modelled together as a combined turbulent diffusion term, *Tij*, using the generalised gradient diffusion hypothesis (GGDH) of Daly & Harlow (1970),

$$T\_{ij} = \frac{\partial}{\partial \mathbf{x}\_l} \left( \mathbf{C}\_s \frac{\varepsilon}{k} \overline{\overline{u\_l} \overline{u\_k}} \frac{\partial \overline{u\_l} \overline{u\_j}}{\partial \mathbf{x}\_k} \right) \tag{32}$$

where *Cs* is typically 0.22.

A deficiency of this model is that it does not preserve the symmetry under cyclic permutation of indices that is exhibited by the triple velocity moments *uiujuk*. This is only significant when the triple moments and pressure diffusion are modelled separately. In such case an improved model that has been suggested by Hanjali´c & Launder (1972) is often used,

$$\overline{u\_{l}u\_{j}u\_{k}} = -\mathsf{C}\_{s}\frac{k}{\varepsilon} \left( \overline{u\_{l}u\_{l}}\frac{\partial \overline{u\_{j}}\overline{u\_{k}}}{\partial \mathbf{x}\_{l}} + \overline{u\_{j}u\_{l}}\frac{\partial \overline{u\_{k}}\overline{u\_{l}}}{\partial \mathbf{x}\_{l}} + \overline{u\_{k}u\_{l}}\frac{\partial \overline{u\_{l}}\overline{u\_{j}}}{\partial \mathbf{x}\_{l}} \right). \tag{33}$$

More elaborate models exist in the literature, as in (Craft, 1998) for example, but the models mentioned above are the ones more commonly used.

geometries. One way to avoid it is based on the observation that the quantity <sup>∇</sup>*k*1/2, evaluated

Reynolds Stress Transport Modelling 13

*<sup>n</sup>* <sup>=</sup> <sup>∇</sup>*k*1/2 |∇*k*1/2|

> *∂k*1/2 *∂x*<sup>2</sup>

*∂k*1/2 *∂xi*

Following Hanjali´c & Launder (1976), when considering the implications of Low-Re effects on dissipation rate modelling, it is instructive to consider the exact transport equation for the

> 2 <sup>−</sup> *<sup>∂</sup> ∂xk*

All the terms on the right hand side above are unclosed, with the exception of viscous diffusion. The first two terms on the right hand side of (42) are the dominant ones in high Re flows. Respectively they represent generation and destruction of *ε*. The third term, which represents a combination of diffusive processes, can be of the same order as the difference of the first two, and must therefore be retained. These three terms are modelled by the three terms that typically appear in high-Re *ε* transport models, as in section 3.4. The fourth and

Launder, 1976), and are thus neglected in high-Re model versions. In low-Re models these terms need to be reconsidered and accounted for if necessary. The last term is often modelled

> *kujuk ε*

This term is present in several Low-Re models developed by the Manchester group. As for the fourth term, initial proposals meant to account for it by allowing the coefficient of the production and destruction terms, *Cε*1and *Cε*2, to be functions of Ret. Similarly, possible viscous effects on the diffusion terms were to be accounted for by allowing the term *C<sup>ε</sup>* to depend on Ret (Hanjali´c & Launder, 1976). However, computations revealed that adding the term in (43) alone was sufficient in producing good agreement between computed energy profiles and available data to within experimental accuracy. Thus dependence of the coefficients *Cε*1, *Cε*2, *C<sup>ε</sup>* on the turbulence Reynolds number is often (not always) abandoned. Finally the viscous diffusion term, neglected in high-Re models, is retained in its exact form.

 *∂*<sup>2</sup>*Ui ∂xj∂xl*

= *Cε*3*ν*

− 2*νuk*

1/2

*∂k*1/2 *∂xj*

*x*2=0

 *ukε* +

*∂ui ∂xl* 2*ν ρ ∂uk ∂xl*

1/2 and Ret smaller than the other terms (Hanjali´c &

 *<sup>∂</sup>*<sup>2</sup>*Ui ∂xk∂xl*

*∂*<sup>2</sup>*Ui ∂xk∂xl* .

= *ε* 2*ν*

, (39)

. (41)

*∂p <sup>∂</sup><sup>x</sup>* <sup>−</sup> *<sup>ν</sup>*

. (40)

*∂ε ∂xk* . (43)

(42)

near a wall, is a vector that points in the wall-normal direction. Thus


energy dissipation rate. This is given by (Daly & Harlow, 1970)

*∂uk ∂xl* − 2 *ν ∂*<sup>2</sup>*ui ∂xk∂xl*

*∂ul ∂xi* + *∂ul ∂xi*

*∂ui ∂xk*

− 2*ν ∂ui ∂xl*

fifth terms are respectively of order Ret

− 2*νuk*

*∂ui ∂xl*

*∂*<sup>2</sup>*Ui ∂xk∂xl*

*∂ui ∂xl*

D*ε* <sup>D</sup>*<sup>t</sup>* <sup>=</sup> <sup>−</sup> <sup>2</sup>*<sup>ν</sup>*

as

and using the value of the dissipation at the wall for a wall with *n* = (0, 1, 0),

Quantities of the form *ninj* appearing in (38) can therefore be replaced by

 *∂k*1/2 *∂x*<sup>2</sup>

*ninj* <sup>=</sup> <sup>2</sup>*<sup>ν</sup> ε*

> *∂ul ∂xk*

 *∂Ui ∂xk*

#### **3.6 Accounting for low-Re effects**

Viscous effects on turbulence properties and their implications on modelling are considered next. The absence of viscous terms in the equation for fluctuating pressure (12) suggests that viscous effects on the fluctuating pressure will be of secondary importance compared to the inviscid effects due to impermeability, considered in section 3.3. The focus of the discussion is thus directed to the dissipation rate tensor, and the transport equation for the scalar dissipation rate. When discussing low-Re effects, reference is frequently made to the turbulent Reynolds number, Ret, defined as

$$\text{Re}\_{\text{fl}} = \frac{k^2}{\nu \varepsilon}.\tag{34}$$

As previously mentioned, at high Reynolds numbers the dissipation rate tensor is assumed to be isotropic, *εij* = <sup>2</sup> <sup>3</sup> *ε δij*. This, however, will cease to be true near a wall where the high anisotropy of the turbulence is expected to be increasingly felt at the smaller scales as the wall is approached. The simplest model accounting for this effect is that of Rotta (1951), which is based on the idea that the anisotropy of the dissipation rate tensor is similar to the stress anisotropy, thus

$$
\varepsilon\_{ij} = \frac{\overline{u\_i u\_j}}{k} \varepsilon.\tag{35}
$$

This model was used by Hanjali´c & Launder (1976) to give the following blending approximation for the dissipation rate tensor

$$
\varepsilon\_{\rm ij} = \frac{2}{3} \varepsilon \left[ (1 - f\_s) \,\delta\_{\rm ij} + f\_s \frac{3}{2} \frac{\overline{u\_i u\_j}}{k} \right], \tag{36}
$$

where *fs* is a function of Ret whose value ranges from 1 to 0 as Ret ranges from 0 to ∞, ensuring the desired behaviour of *εij* in these limits. The near-wall model (35) is the simplest form accounting for near-wall anisotropy of the dissipation tensor. Launder & Reynolds (1983) have shown that this form does not give the correct near-wall asymptotic behaviour of the individual tensor elements, which are rather given by

$$\begin{aligned} \frac{\varepsilon\_{ij}}{\varepsilon} &= \frac{\overline{u\_i u\_j}}{k}, \quad i \neq 2, j \neq 2\\ \frac{\varepsilon\_{12}}{\varepsilon} &= 2\frac{\overline{u\_i u\_2}}{k}, \quad i \neq 2\\ \frac{\varepsilon\_{22}}{\varepsilon} &= 4\frac{\overline{u\_2 u\_2}}{k}. \end{aligned} \tag{37}$$

What is needed then is a term to replace the Rotta model in (36) which yields the correct asymptotic behaviour described by (37), and which contracts to 2*ε*. One possible form that satisfies these requirements is

$$\varepsilon\_{ij}^{\*} = \frac{\varepsilon/k\left(\overline{u\_i u\_j} + \overline{u\_i u\_k}\ n\_j \ n\_k + \overline{u\_j u\_k}\ n\_i \ n\_k + \overline{u\_k u\_l}\ n\_k \ n\_l \ \delta\_{ij}\right)}{\left(1 + \frac{5}{2} n\_p n\_q \overline{u\_p u\_q}/k\right)}\tag{38}$$

where *ni* represents a component of the wall-normal unit vector (Pope, 2000). The use of the wall vector in a model is undesirable because of the ambiguity it introduces in complex 10 Numerical Simulations

Viscous effects on turbulence properties and their implications on modelling are considered next. The absence of viscous terms in the equation for fluctuating pressure (12) suggests that viscous effects on the fluctuating pressure will be of secondary importance compared to the inviscid effects due to impermeability, considered in section 3.3. The focus of the discussion is thus directed to the dissipation rate tensor, and the transport equation for the scalar dissipation rate. When discussing low-Re effects, reference is frequently made to the

Ret <sup>=</sup> *<sup>k</sup>*<sup>2</sup>

As previously mentioned, at high Reynolds numbers the dissipation rate tensor is assumed

anisotropy of the turbulence is expected to be increasingly felt at the smaller scales as the wall is approached. The simplest model accounting for this effect is that of Rotta (1951), which is based on the idea that the anisotropy of the dissipation rate tensor is similar to the stress

*<sup>ε</sup>ij* <sup>=</sup> *uiuj*

This model was used by Hanjali´c & Launder (1976) to give the following blending

(1 − *fs*) *δij* + *fs*

where *fs* is a function of Ret whose value ranges from 1 to 0 as Ret ranges from 0 to ∞, ensuring the desired behaviour of *εij* in these limits. The near-wall model (35) is the simplest form accounting for near-wall anisotropy of the dissipation tensor. Launder & Reynolds (1983) have shown that this form does not give the correct near-wall asymptotic behaviour of the

*uiu*<sup>2</sup>

*u*2*u*<sup>2</sup> *k* .

What is needed then is a term to replace the Rotta model in (36) which yields the correct asymptotic behaviour described by (37), and which contracts to 2*ε*. One possible form that

where *ni* represents a component of the wall-normal unit vector (Pope, 2000). The use of the wall vector in a model is undesirable because of the ambiguity it introduces in complex

*<sup>k</sup>* , *<sup>i</sup>* �<sup>=</sup> 2, *<sup>j</sup>* �<sup>=</sup> <sup>2</sup>

*uiuj* + *uiuk nj nk* + *ujuk ni nk* + *ukul nk nl δij*

<sup>2</sup> *npnqupuq*/*k*

*<sup>k</sup>* , *<sup>i</sup>* �<sup>=</sup> <sup>2</sup>

<sup>3</sup> *ε δij*. This, however, will cease to be true near a wall where the high

3 2 *uiuj k* 

*νε* . (34)

*<sup>k</sup> <sup>ε</sup>*. (35)

, (36)

, (38)

(37)

**3.6 Accounting for low-Re effects**

to be isotropic, *εij* = <sup>2</sup>

anisotropy, thus

turbulent Reynolds number, Ret, defined as

approximation for the dissipation rate tensor

individual tensor elements, which are rather given by

satisfies these requirements is

*ε* ∗ *ij* <sup>=</sup> *<sup>ε</sup>*/*<sup>k</sup>* *<sup>ε</sup>ij* <sup>=</sup> <sup>2</sup> 3 *ε* 

> *εij <sup>ε</sup>* <sup>=</sup> *uiuj*

*ε*12 *<sup>ε</sup>* <sup>=</sup> <sup>2</sup>

*ε*22 *<sup>ε</sup>* <sup>=</sup> <sup>4</sup>

> 1 + <sup>5</sup>

geometries. One way to avoid it is based on the observation that the quantity <sup>∇</sup>*k*1/2, evaluated near a wall, is a vector that points in the wall-normal direction. Thus

$$
\vec{n} = \frac{\nabla k^{1/2}}{|\nabla k^{1/2}|},\tag{39}
$$

and using the value of the dissipation at the wall for a wall with *n* = (0, 1, 0),

$$|\nabla k^{1/2}|\_{\mathbf{x}\_2=0} = \left(\frac{\partial k^{1/2}}{\partial \mathbf{x}\_2} \frac{\partial k^{1/2}}{\partial \mathbf{x}\_2}\right)^{1/2}\_{\mathbf{x}\_2=0} = \sqrt{\frac{\varepsilon}{2\nu}}.\tag{40}$$

Quantities of the form *ninj* appearing in (38) can therefore be replaced by

$$m\_i n\_j = \frac{2\nu}{\varepsilon} \frac{\partial k^{1/2}}{\partial x\_i} \frac{\partial k^{1/2}}{\partial x\_j}. \tag{41}$$

Following Hanjali´c & Launder (1976), when considering the implications of Low-Re effects on dissipation rate modelling, it is instructive to consider the exact transport equation for the energy dissipation rate. This is given by (Daly & Harlow, 1970)

$$\begin{split} \frac{\mathrm{D}\varepsilon}{\mathrm{D}t} &= -2\nu \overline{\frac{\partial u\_{i}}{\partial \mathrm{x}\_{k}} \frac{\partial \boldsymbol{u}\_{i}}{\partial \mathrm{x}\_{l}} \frac{\partial \boldsymbol{u}\_{k}}{\partial \mathrm{x}\_{l}}} - 2\left(\nu \frac{\partial^{2} u\_{i}}{\partial \mathrm{x}\_{k} \partial \mathrm{x}\_{l}}\right)^{2} - \frac{\partial}{\partial \mathrm{x}\_{k}} \left[ \overline{\boldsymbol{u}\_{k}} \varepsilon + \frac{2\nu}{\rho} \overline{\frac{\partial \boldsymbol{u}\_{k}}{\partial \mathrm{x}\_{l}} \frac{\partial \boldsymbol{p}}{\partial \mathrm{x}}} - \nu \frac{\partial \varepsilon}{\partial \mathrm{x}\_{k}} \right] \\ &- 2\nu \left( \overline{\frac{\partial \boldsymbol{u}\_{i}}{\partial \mathrm{x}\_{l}} \frac{\partial \boldsymbol{u}\_{l}}{\partial \mathrm{x}\_{i}}} + \overline{\frac{\partial \boldsymbol{u}\_{l}}{\partial \mathrm{x}\_{i}} \frac{\partial \boldsymbol{u}\_{l}}{\partial \mathrm{x}\_{k}}} \right) \frac{\partial \boldsymbol{U}\_{i}}{\partial \mathrm{x}\_{k}} - 2\nu \overline{\boldsymbol{u}\_{k} \overline{\frac{\partial \boldsymbol{u}\_{i}}{\partial \mathrm{x}\_{l}}}} \frac{\partial^{2} \boldsymbol{U}\_{i}}{\partial \mathrm{x}\_{k} \partial \mathrm{x}\_{l}}. \end{split} \tag{42}$$

All the terms on the right hand side above are unclosed, with the exception of viscous diffusion. The first two terms on the right hand side of (42) are the dominant ones in high Re flows. Respectively they represent generation and destruction of *ε*. The third term, which represents a combination of diffusive processes, can be of the same order as the difference of the first two, and must therefore be retained. These three terms are modelled by the three terms that typically appear in high-Re *ε* transport models, as in section 3.4. The fourth and fifth terms are respectively of order Ret 1/2 and Ret smaller than the other terms (Hanjali´c & Launder, 1976), and are thus neglected in high-Re model versions. In low-Re models these terms need to be reconsidered and accounted for if necessary. The last term is often modelled as

$$-2\nu \overline{u\_k \frac{\partial \boldsymbol{u}\_l}{\partial \mathbf{x}\_l}} \frac{\partial^2 \mathcal{U}\_l}{\partial \mathbf{x}\_k \partial \mathbf{x}\_l} = \mathsf{C}\_{\varepsilon 3} \nu \frac{k \overline{u\_j \boldsymbol{u}\_k}}{\varepsilon} \left(\frac{\partial^2 \mathcal{U}\_l}{\partial \mathbf{x}\_j \partial \mathbf{x}\_l}\right) \left(\frac{\partial^2 \mathcal{U}\_l}{\partial \mathbf{x}\_k \partial \mathbf{x}\_l}\right). \tag{43}$$

This term is present in several Low-Re models developed by the Manchester group. As for the fourth term, initial proposals meant to account for it by allowing the coefficient of the production and destruction terms, *Cε*1and *Cε*2, to be functions of Ret. Similarly, possible viscous effects on the diffusion terms were to be accounted for by allowing the term *C<sup>ε</sup>* to depend on Ret (Hanjali´c & Launder, 1976). However, computations revealed that adding the term in (43) alone was sufficient in producing good agreement between computed energy profiles and available data to within experimental accuracy. Thus dependence of the coefficients *Cε*1, *Cε*2, *C<sup>ε</sup>* on the turbulence Reynolds number is often (not always) abandoned. Finally the viscous diffusion term, neglected in high-Re models, is retained in its exact form.

**The Shima low-Re model**

where *Cs* = 0.22.

as:

The dissipation equation is given by

D*ε* <sup>D</sup>*<sup>t</sup>* <sup>=</sup> *<sup>C</sup>ε*<sup>1</sup>

*λ*∗ = *∂ ∂xi*

In its original form, the Launder & Shima (1989) model is a low-Re version of the Basic model that uses wall reflection terms and includes Ret-based damping coefficients to return the correct near-wall behaviour. Shima (1998) later proposed a low-Re model, based on the QI pressure-strain rate model, that does away with the wall reflection terms in the interest of more general applicability to complex geometries. The model admittedly gives stress anisotropy results in steady channel flow that are inferior to his previous low-Re formulation, but this is a compromise made in order to discard the wall reflection terms with their associated difficulties related to complex geometries. The pressure-strain rate coefficients are

Reynolds Stress Transport Modelling 15

<sup>2</sup> *<sup>A</sup>*0.75[<sup>1</sup> <sup>−</sup> exp(−49*A*2)] × {<sup>1</sup> <sup>−</sup> exp[−(*Ret*/60)2]} (52a)

*A*<sup>2</sup> = *aijaji A*<sup>3</sup> = *aijajk aki* . (53)

<sup>2</sup> *C*<sup>3</sup> (52d)

<sup>8</sup> (*A*<sup>2</sup> − *A*3). (54)

. (57)

. (58e)

(55)

. (56)

*C*<sup>2</sup> = 0.7*A* (52b) *C*<sup>3</sup> = 0.3*A*0.5 (52c)

> *∂uiuj ∂xl*

<sup>2</sup>

*Cε*<sup>1</sup> = 1.44 + *β*<sup>1</sup> + *β*<sup>2</sup> , (58a) *β*<sup>1</sup> = 0.25*A* min(*λ*/2.5 − 1, 0) − 1.4*A* min(*P*/*ε* − 1, 0), (58b) *<sup>β</sup>*<sup>2</sup> <sup>=</sup> 1.0*Aλ*<sup>2</sup> max(*λ*/2.5 <sup>−</sup> 1, 0), (58c) *λ* = min(*λ*∗, 4), (58d)

no longer constant, and are given by the following expressions:

*<sup>C</sup>*<sup>4</sup> <sup>=</sup> 0.65*A*(0.23*C*<sup>1</sup> <sup>+</sup> *<sup>C</sup>*<sup>2</sup> <sup>−</sup> <sup>1</sup>) + 1.3*A*0.25

and *A* is the 'flatness' parameter first defined by Lumley (1978),

*ε k P* − *Cε*<sup>2</sup>

 *k*1.5 *ε*

 *∂ ∂xi*  *k*1.5 *ε*

where *ε*˜ is the homogeneous dissipation rate, defined as:

where *A*2, *A*<sup>3</sup> are the second and third invariants of the stress anisotropy tensor:

*<sup>A</sup>* <sup>=</sup> <sup>1</sup> <sup>−</sup> <sup>9</sup>

correlation, is modelled using the simple gradient diffusion of Daly & Harlow (1970)

 *Cs k ε ukul*

*Tij* <sup>=</sup> *<sup>∂</sup> ∂xk*

> *εε*˜ *<sup>k</sup>* <sup>+</sup>

*ε*˜ = *ε* − 2*ν*

Turbulent diffusion, comprising the triple velocity correlation and the pressure velocity

*∂ ∂xk*

> *∂k*1/2 *∂xi*

The coefficients *Cε*2,*C<sup>ε</sup>* retain their typical values 1.92, 0.15 respectively, but *Cε*<sup>1</sup> is prescribed

 *Cε k ε ukul ∂ε ∂xl* + *ν ∂ε ∂xk*

*C*<sup>1</sup> = 1 + 2.45*A*0.25

#### **3.7 The Launder–Reece–Rodi models**

In their seminal 1975 paper, Launder, Reece & Rodi laid out a hierarchy of RST models based on arguments presented in section 3.2. Two rapid pressure-strain rate models were proposed. The first is the quasi-isotropic model (LRR-QI), which has the most general linear tensorial form satisfying the required symmetry conditions, and is given by

$$\Phi\_{\rm ij}^{r} = -\mathsf{C}\_{2}(P\_{\rm ij} - \frac{2}{3}\delta\_{\rm ij}P\_{\rm k}) - \mathsf{C}\_{3}(D\_{\rm ij} - \frac{2}{3}\delta\_{\rm ij}P\_{\rm k}) - 2\mathsf{C}\_{4}k \,\mathsf{S}\_{\rm ij} \tag{44}$$

where *Sij* is the mean strain rate tensor, defined as:

$$S\_{i\bar{j}} = \frac{1}{2} \left( \frac{\partial \mathcal{U}\_i}{\partial \mathbf{x}\_{\bar{j}}} + \frac{\partial \mathcal{U}\_{\bar{j}}}{\partial \mathbf{x}\_{\bar{i}}} \right),\tag{45}$$

and the coefficients have the following values

$$\mathbf{C}\_2 = 0.764, \quad \mathbf{C}\_3 = 0.182, \quad \mathbf{C}\_4 = 0.109 \, . \tag{46}$$

The second rapid pressure-strain rate model is the isotropization of production model (LRR-IP), which is also referred to as the 'Basic' model, and simply retains the first term of the QI model and neglects the other two. Thus,

$$\Phi\_{\rm ij}^r = -\mathsf{C}\_2(P\_{\rm ij} - \frac{2}{3}\delta\_{\rm ij}P\_{\rm k}),\tag{47}$$

where the coefficient *C*<sup>2</sup> is now set at 0.6. Both models use the Rotta return-to-isotropy model for the slow pressure-strain rate term,

$$
\phi\_{\rm ij}^{\rm s} = -\mathbb{C}\_1 \varepsilon a\_{\rm ij} \,. \tag{48}
$$

but the coefficient *C*<sup>1</sup> is 1.5 for the QI model and 1.8 for the IP model.

In the original proposal turbulent diffusion *Tij* is modelled using (33) for the triple velocity moments (pressure diffusion is usually neglected). In many later implementations this is replaced by the simpler GGDH. Thus the models can be written as

$$\frac{\mathbf{D}\overline{u\_{l}u\_{j}}}{\mathbf{D}t} = \mathbf{P}\_{lj} - \mathbf{C}\_{l}\varepsilon a\_{lj} + \boldsymbol{\phi}\_{lj}^{r} + \boldsymbol{\phi}\_{lj}^{s,w} + \boldsymbol{\phi}\_{lj}^{r,w} + \frac{\partial}{\partial \mathbf{x}\_{l}} \left( \mathbf{C}\_{s}\frac{\varepsilon}{\overline{k}} \overline{u\_{l}u\_{k}} \frac{\partial \overline{u\_{l}u\_{j}}}{\partial \mathbf{x}\_{k}} \right) - \frac{2}{3}\delta\_{lj}\varepsilon \tag{49}$$

where *φ<sup>r</sup> ij* is replaced by either the QI or IP models, and the wall-reflection terms *<sup>φ</sup>s*,*<sup>w</sup> ij* , *<sup>φ</sup>r*,*<sup>w</sup> ij* are given by (27) and (28), respectively. Since these models are intended as high Re models, the viscous diffusion term is neglected and an isotropic dissipation rate tensor is assumed. Finally, closure is completed with the standard high-Re dissipation rate equation, given by

$$\frac{\partial \mathbf{D}\varepsilon}{\partial \mathbf{D}t} = \mathbf{C}\_{\varepsilon 1} \frac{\varepsilon}{k} P\_{\mathbf{k}} - \mathbf{C}\_{\varepsilon 2} \frac{\varepsilon^2}{k} + \frac{\partial}{\partial \mathbf{x}\_k} \left( \mathbf{C}\_{\varepsilon} \frac{k}{\varepsilon} \overline{\mu\_k \overline{\mu\_l}} \frac{\partial \varepsilon}{\partial \mathbf{x}\_l} \right) \,, \tag{50}$$

where

$$\mathbb{C}\_{\varepsilon1} = 1.44, \quad \mathbb{C}\_{\varepsilon2} = 1.92, \quad \mathbb{C}\_{\varepsilon} = 0.15. \tag{51}$$

#### **The Shima low-Re model**

12 Numerical Simulations

In their seminal 1975 paper, Launder, Reece & Rodi laid out a hierarchy of RST models based on arguments presented in section 3.2. Two rapid pressure-strain rate models were proposed. The first is the quasi-isotropic model (LRR-QI), which has the most general linear tensorial

<sup>3</sup> *<sup>δ</sup>ij <sup>P</sup><sup>κ</sup>* ) <sup>−</sup> *<sup>C</sup>*3(*Dij* <sup>−</sup> <sup>2</sup>

 *∂Ui ∂xj* + *∂Uj ∂xi*

The second rapid pressure-strain rate model is the isotropization of production model (LRR-IP), which is also referred to as the 'Basic' model, and simply retains the first term of

where the coefficient *C*<sup>2</sup> is now set at 0.6. Both models use the Rotta return-to-isotropy model

In the original proposal turbulent diffusion *Tij* is modelled using (33) for the triple velocity moments (pressure diffusion is usually neglected). In many later implementations this is

> *ij* <sup>+</sup> *<sup>φ</sup>r*,*<sup>w</sup> ij* +

*ij* is replaced by either the QI or IP models, and the wall-reflection terms *<sup>φ</sup>s*,*<sup>w</sup>*

given by (27) and (28), respectively. Since these models are intended as high Re models, the viscous diffusion term is neglected and an isotropic dissipation rate tensor is assumed. Finally, closure is completed with the standard high-Re dissipation rate equation, given by

> *ε*2 *<sup>k</sup>* <sup>+</sup>

*∂ ∂xk*  *Cε k ε ukul ∂ε ∂xl* 

*Cε*<sup>1</sup> = 1.44, *Cε*<sup>2</sup> = 1.92, *C<sup>ε</sup>* = 0.15. (51)

*∂ ∂xl*  *Cs ε k uluk*

*ij* <sup>=</sup> <sup>−</sup>*C*2(*Pij* <sup>−</sup> <sup>2</sup>

*φs*

*ij* <sup>+</sup> *<sup>φ</sup>s*,*<sup>w</sup>*

*P<sup>κ</sup>* − *Cε*<sup>2</sup>

*C*<sup>2</sup> = 0.764, *C*<sup>3</sup> = 0.182, *C*<sup>4</sup> = 0.109 . (46)

<sup>3</sup> *δij P<sup>κ</sup>* ) − 2*C*4*k Sij*, (44)

, (45)

<sup>3</sup> *δij P<sup>κ</sup>* ), (47)

*ij* = −*C*1*εaij* , (48)

*∂uiuj ∂xk*

 − 2

<sup>3</sup> *δijε* (49)

*ij* , *<sup>φ</sup>r*,*<sup>w</sup> ij* are

, (50)

form satisfying the required symmetry conditions, and is given by

*Sij* <sup>=</sup> <sup>1</sup> 2

*φr*

but the coefficient *C*<sup>1</sup> is 1.5 for the QI model and 1.8 for the IP model.

replaced by the simpler GGDH. Thus the models can be written as

*ε k*

*ij* <sup>=</sup> <sup>−</sup>*C*2(*Pij* <sup>−</sup> <sup>2</sup>

where *Sij* is the mean strain rate tensor, defined as:

and the coefficients have the following values

the QI model and neglects the other two. Thus,

<sup>D</sup>*<sup>t</sup>* <sup>=</sup> *Pij* <sup>−</sup> *<sup>C</sup>*1*εaij* <sup>+</sup> *<sup>φ</sup><sup>r</sup>*

D*ε* <sup>D</sup>*<sup>t</sup>* <sup>=</sup> *<sup>C</sup>ε*<sup>1</sup>

for the slow pressure-strain rate term,

D*uiuj*

where *φ<sup>r</sup>*

where

**3.7 The Launder–Reece–Rodi models**

*φr*

In its original form, the Launder & Shima (1989) model is a low-Re version of the Basic model that uses wall reflection terms and includes Ret-based damping coefficients to return the correct near-wall behaviour. Shima (1998) later proposed a low-Re model, based on the QI pressure-strain rate model, that does away with the wall reflection terms in the interest of more general applicability to complex geometries. The model admittedly gives stress anisotropy results in steady channel flow that are inferior to his previous low-Re formulation, but this is a compromise made in order to discard the wall reflection terms with their associated difficulties related to complex geometries. The pressure-strain rate coefficients are no longer constant, and are given by the following expressions:

$$C\_1 = 1 + 2.45A\_2^{0.25}A^{0.75}[1 - \exp(-49A^2)] \times \left\{1 - \exp[-(\text{Re}\_l/\text{60})^2] \right\} \tag{52a}$$

$$C\_2 = 0.7A \tag{52b}$$

$$C\_3 = 0.3A^{0.5} \tag{52c}$$

$$\mathbf{C\_4} = 0.65A(0.23\mathbf{C\_1} + \mathbf{C\_2} - 1) + 1.3A\_2^{0.25}\mathbf{C\_3} \tag{52d}$$

where *A*2, *A*<sup>3</sup> are the second and third invariants of the stress anisotropy tensor:

$$A\_2 = a\_{i\bar{j}} a\_{j\bar{l}} \quad A\_3 = a\_{i\bar{j}} a\_{j\bar{k}} a\_{k\bar{l}} \,. \tag{53}$$

and *A* is the 'flatness' parameter first defined by Lumley (1978),

$$A = 1 - \frac{9}{8}(A\_2 - A\_3). \tag{54}$$

Turbulent diffusion, comprising the triple velocity correlation and the pressure velocity correlation, is modelled using the simple gradient diffusion of Daly & Harlow (1970)

$$T\_{lj} = \frac{\partial}{\partial \mathbf{x}\_k} \left( \mathbf{C}\_s \frac{k}{\varepsilon} \overline{\mu\_k \mu\_l} \frac{\partial \overline{u\_l u\_j}}{\partial \mathbf{x}\_l} \right) \tag{55}$$

where *Cs* = 0.22.

The dissipation equation is given by

$$\frac{\mathrm{D}\varepsilon}{\mathrm{D}t} = \mathsf{C}\_{\varepsilon1}\frac{\varepsilon}{k}\mathrm{P} - \mathsf{C}\_{\varepsilon2}\frac{\varepsilon\widetilde{\varepsilon}}{k} + \frac{\partial}{\partial \mathsf{x}\_{k}} \left( \mathsf{C}\_{\varepsilon}\frac{k}{\varepsilon} \overline{\mu\_{k}\mu\_{l}} \frac{\partial \varepsilon}{\partial \mathsf{x}\_{l}} + \nu \frac{\partial \varepsilon}{\partial \mathsf{x}\_{k}} \right) \,. \tag{56}$$

where *ε*˜ is the homogeneous dissipation rate, defined as:

$$
\tilde{\varepsilon} = \varepsilon - 2\nu \left( \frac{\partial k^{1/2}}{\partial x\_i} \right)^2 \,. \tag{57}
$$

The coefficients *Cε*2,*C<sup>ε</sup>* retain their typical values 1.92, 0.15 respectively, but *Cε*<sup>1</sup> is prescribed as:

$$\mathbb{C}\_{\varepsilon 1} = 1.44 + \beta\_1 + \beta\_2. \tag{58a}$$

$$\beta\_1 = 0.25A \min(\lambda/2.5 - 1, 0) - 1.4A \min(P/\varepsilon - 1, 0) \,\,\,\tag{58b}$$

$$
\beta\_2 = 1.0A\lambda^2 \max(\lambda/2.5 - 1, 0) \,\text{.}\tag{58c}
$$

$$
\lambda = \min(\lambda^\*, 4),
\tag{58d}
$$

$$
\lambda^\* = \left[ \frac{\partial}{\partial \mathbf{x}\_i} \left( \frac{k^{1.5}}{\varepsilon} \right) \frac{\partial}{\partial \mathbf{x}\_i} \left( \frac{k^{1.5}}{\varepsilon} \right) \right]. \tag{58e}
$$

The coefficients are specified by:

*C<sup>w</sup>*

D*ε* <sup>D</sup>*<sup>t</sup>* <sup>=</sup> *<sup>C</sup>ε*<sup>1</sup>

and

where *ε*∗

The modelled dissipation rate transport equation is given by:

*P<sup>κ</sup>* − *Cε*<sup>2</sup> *f<sup>ε</sup>*

*<sup>f</sup><sup>ε</sup>* <sup>=</sup> <sup>1</sup> <sup>−</sup> *<sup>C</sup>ε*<sup>2</sup> <sup>−</sup> 1.4

*Cl*

*εij* = *fsε* ∗

*Cε*<sup>2</sup>

*∂l ∂xn* <sup>2</sup> − 1

*ε k*

+ *Cε*<sup>3</sup> *ν k ε uiuj*

The coefficients have the following specified values:

The length-scale growth correction, *Sl*, is given by:

where *l* = *k*3/2/*ε*, and *Cl* = 2.5.

*ij* is given by:

*ε* ∗ *ij* <sup>=</sup> *<sup>ε</sup> k*

*Sl* <sup>=</sup> max <sup>1</sup>

The anisotropic stress dissipation rate tensor is modelled as:

*<sup>C</sup>*<sup>1</sup> <sup>=</sup> *<sup>C</sup>* <sup>+</sup> <sup>√</sup>*AE*2, *<sup>C</sup>* <sup>=</sup> 2.5*AF*1/4 *<sup>f</sup>* , *<sup>F</sup>* <sup>=</sup> min(0.6, *<sup>A</sup>*2), (66)

*<sup>f</sup>* <sup>=</sup> min Ret

The damping coefficient appearing in the wall correction terms (27) and (28) is given by:

*fw* <sup>=</sup> min *<sup>k</sup>*3/2

*εε*˜ *<sup>k</sup>* <sup>+</sup>

*∂*<sup>2</sup>*Uk ∂xi∂xl*

<sup>2</sup> <sup>=</sup> max(<sup>1</sup> <sup>−</sup> 0.7*C*, 0.3), *<sup>C</sup><sup>w</sup>*

150

Reynolds Stress Transport Modelling 17

2.5*εxn*

*∂ ∂xk*

*∂*<sup>2</sup>*Uk ∂xj∂xl*

> exp − Ret 6

*ij* + (1 − *fs*)

*uiuj* + (*uiuk nj nk* + *ujuk ni nk* + *ukul nk nl ni nj*)*fd*

1 + <sup>3</sup> <sup>2</sup> *npnq*

 <sup>1</sup> *Cl*

*∂l ∂xn*

2 3

*upuq <sup>k</sup> fd*

*fs* <sup>=</sup> <sup>1</sup> <sup>−</sup> <sup>√</sup>*AE*2, *fd* = (<sup>1</sup> <sup>+</sup> 0.1Ret)<sup>−</sup>1. (77)

 *Cε k ε ukul ∂ε ∂xl* + *ν ∂ε*˜ *∂xk*

3/2 , 1

, 1.4

*Cε*<sup>1</sup> = 2.6 *Cε*<sup>2</sup> = 1.92 *Cε*<sup>3</sup> = 0.25 *Cε*<sup>4</sup> = 0.1 *C<sup>ε</sup>* = 0.18 , (72)

− *Cε*<sup>4</sup> *f*4*k* Ω*<sup>k</sup>* Ω*<sup>k</sup>* + *Sl*.

<sup>2</sup> 

<sup>2</sup> , 0 *εε*˜

*<sup>C</sup>*<sup>2</sup> <sup>=</sup> 0.8√*<sup>A</sup>* , (67)

, (68)

. (70)

. (73)

*<sup>k</sup> <sup>A</sup>* , (74)

, (76)

*δijε*, (75)

(71)

<sup>2</sup> = min(*A*, 0.3). (69)

#### **3.8 The Speziale–Sarkar–Gatski model**

Speziale et al. (1991) developed a pressure-strain rate model that is quadratic in *aij* by first considering the most general form for *φij* (slow and rapid) that is linear in the mean strain and rotation tensors and quadratic in *aij*. Then they obtained their model by considering the simplest subset of that general form that has an equivalent structural equilibrium in plane homogeneous flows. The resulting model has a rapid part that is linear in *aij*, and a quadratic slow part, given by

$$\begin{split} \phi\_{ij} &= -\left(2d\_1\varepsilon + d\_1^\* P\_\mathbf{k}\right) \frac{a\_{ij}}{2} + \frac{d\_2}{4} \varepsilon (a\_{ik} a\_{kj} - \frac{1}{3} a\_{kl} a\_{kl} \delta\_{ij}) \\ &+ \left(d\_3 - d\_3^\* \frac{\sqrt{A\_2}}{2}\right) k S\_{ij} + \frac{d\_4}{2} k (a\_{ik} S\_{jk} + a\_{jk} S\_{ik} - \frac{2}{3} a\_{kl} S\_{kl} \delta\_{ij}) \\ &+ \frac{d\_5}{2} k (a\_{ik} \Omega\_{jk} + a\_{jk} \Omega\_{ik}) \end{split} \tag{59}$$

where Ω*ij* is the mean vorticity tensor defined as:

$$
\Omega\_{ij} = \frac{1}{2} \left( \frac{\partial \mathcal{U}\_i}{\partial \mathbf{x}\_j} - \frac{\partial \mathcal{U}\_j}{\partial \mathbf{x}\_i} \right),
\tag{60}
$$

and the coefficients have the following values

$$d\_1 = 1.7, \quad d\_1^\* = 1.8, \quad d\_2 = 4.2, \quad d\_3 = \frac{4}{5}, \quad d\_3^\* = 1.3, \quad d\_4 = 1.25, \quad d\_5 = 0.4. \tag{61}$$

The rapid part of the SSG model, aside from the nonlinear dependence on *A*<sup>2</sup> in third term of (59), is tensorially equivalent to the QI model.

Diffusion is modelled using the GGDH, and the standard high-Re version of the *ε* equation (50) is used, but the coefficient *Cε*<sup>2</sup> is assigned the slightly lower value of 1.83.

#### **3.9 The Hanjalic–Jakirli ´ c low-Re model ´**

Jakirli´c & Hanjali´c (1995) developed a low-Re RSTM that is based on the LRR-IP model, and the Gibson & Launder (1978) wall corrections (27) and (28), making modifications to handle Low-Re and near-wall effects. The modifications are expressed in terms of Ret, the stress anisotropy invariants, *A*2, *A*3, in addition to invariants of the stress dissipation rate anisotropy tensor, *E*2, *E*3, defined as:

$$E\_2 = e\_{\rm ij} e\_{\rm ji} \quad E\_3 = e\_{\rm ij} e\_{\rm jk} e\_{\rm ki} \,\,\,\,\,\tag{62}$$

$$
\varepsilon\_{ij} = \frac{\varepsilon\_{ij}}{\varepsilon} - \frac{2}{3} \delta\_{ij} \,. \tag{63}
$$

A 'flatness' parameter based on the stress dissipation rate anisotropy invariants is also used:

$$E = 1 - \frac{9}{8}(E\_2 - E\_3) \,. \tag{64}$$

The modelled RST equation is given by:

$$\begin{split} \frac{\partial \overline{u\_{i}u\_{j}}}{\mathbf{D}t} &= \mathbf{P}\_{ij} - \mathbf{C}\_{1}\varepsilon a\_{ij} - \mathbf{C}\_{2} \left( P\_{ij} - \frac{2}{3} \delta\_{ij} P\_{\mathbf{k}} \right) + \boldsymbol{\phi}\_{ij}^{s,\overline{\mathbf{w}}} + \boldsymbol{\phi}\_{ij}^{r,\overline{\mathbf{w}}} \\ &+ \frac{\partial}{\partial \mathbf{x}\_{l}} \left( \mathbf{C}\_{s} \frac{\varepsilon}{\overline{\mathbf{u}\_{l}} \overline{\mathbf{u}\_{k}}} \frac{\partial \overline{u\_{l}} \overline{\mathbf{u}\_{j}}}{\partial \mathbf{x}\_{k}} \right) - \varepsilon\_{ij} \,. \end{split} \tag{65}$$

14 Numerical Simulations

Speziale et al. (1991) developed a pressure-strain rate model that is quadratic in *aij* by first considering the most general form for *φij* (slow and rapid) that is linear in the mean strain and rotation tensors and quadratic in *aij*. Then they obtained their model by considering the simplest subset of that general form that has an equivalent structural equilibrium in plane homogeneous flows. The resulting model has a rapid part that is linear in *aij*, and a quadratic

<sup>4</sup> *<sup>ε</sup>*(*aik akj* <sup>−</sup> <sup>1</sup>

*d*4

 *∂Ui ∂xj*

The rapid part of the SSG model, aside from the nonlinear dependence on *A*<sup>2</sup> in third term of

Diffusion is modelled using the GGDH, and the standard high-Re version of the *ε* equation

Jakirli´c & Hanjali´c (1995) developed a low-Re RSTM that is based on the LRR-IP model, and the Gibson & Launder (1978) wall corrections (27) and (28), making modifications to handle Low-Re and near-wall effects. The modifications are expressed in terms of Ret, the stress anisotropy invariants, *A*2, *A*3, in addition to invariants of the stress dissipation rate anisotropy

> *<sup>ε</sup>* <sup>−</sup> <sup>2</sup> 3

A 'flatness' parameter based on the stress dissipation rate anisotropy invariants is also used:

8

 *Pij* <sup>−</sup> <sup>2</sup> 3 *δijP<sup>κ</sup>* + *φs*,*<sup>w</sup>*

*∂uiuj ∂xk*

 − *εij* .

*eij* <sup>=</sup> *<sup>ε</sup>ij*

*<sup>E</sup>* <sup>=</sup> <sup>1</sup> <sup>−</sup> <sup>9</sup>

<sup>−</sup> *<sup>∂</sup>Uj ∂xi*

<sup>5</sup> , *d*<sup>∗</sup>

<sup>3</sup> *akl aklδij*)

<sup>3</sup> *aklSklδij*)

, (60)

<sup>3</sup> = 1.3, *d*<sup>4</sup> = 1.25, *d*<sup>5</sup> = 0.4 . (61)

*δij* . (63)

(*E*<sup>2</sup> − *E*3). (64)

*ij* <sup>+</sup> *<sup>φ</sup>r*,*<sup>w</sup> ij*

(59)

(65)

<sup>2</sup> *<sup>k</sup>*(*aikSjk* <sup>+</sup> *ajkSik* <sup>−</sup> <sup>2</sup>

*E*<sup>2</sup> = *eijeji E*<sup>3</sup> = *eijejkeki* , (62)

**3.8 The Speziale–Sarkar–Gatski model**

*φij* = − (2*d*1*ε* + *d*<sup>∗</sup>

where Ω*ij* is the mean vorticity tensor defined as:

and the coefficients have the following values

(59), is tensorially equivalent to the QI model.

**3.9 The Hanjalic–Jakirli ´ c low-Re model ´**

The modelled RST equation is given by:

D*uiuj*

+ *∂ ∂xl*

<sup>D</sup>*<sup>t</sup>* <sup>=</sup>*Pij* <sup>−</sup> *<sup>C</sup>*1*εaij* <sup>−</sup> *<sup>C</sup>*<sup>2</sup>

 *Cs ε k uluk*

*d*<sup>1</sup> = 1.7, *d*<sup>∗</sup>

tensor, *E*2, *E*3, defined as:

*d*<sup>3</sup> − *d*<sup>∗</sup> 3 <sup>√</sup>*A*<sup>2</sup> 2

+ 

+ *d*5 <sup>1</sup>*P<sup>κ</sup>* ) *aij* <sup>2</sup> <sup>+</sup> *d*2

<sup>2</sup> *<sup>k</sup>*(*aik*Ω*jk* <sup>+</sup> *ajk*Ω*ik*),

<sup>1</sup> = 1.8, *<sup>d</sup>*<sup>2</sup> = 4.2, *<sup>d</sup>*<sup>3</sup> = <sup>4</sup>

 *kSij* +

<sup>Ω</sup>*ij* <sup>=</sup> <sup>1</sup> 2

(50) is used, but the coefficient *Cε*<sup>2</sup> is assigned the slightly lower value of 1.83.

slow part, given by

The coefficients are specified by:

$$\mathbb{C}\_1 = \mathbb{C} + \sqrt{A}\mathbb{E}^2, \quad \mathbb{C} = 2.5 A \mathbb{F}^{1/4} f, \quad F = \min(0.6, A\_2) \,\,\,\tag{66}$$

$$C\_2 = 0.8\sqrt{A} \,\, \tag{67}$$

$$f = \min\left[\left(\frac{\mathbf{Re}\_t}{150}\right)^{3/2}, 1\right],\tag{68}$$

$$\mathbb{C}\_2^w = \max(1 - 0.7\mathbb{C}, 0.3), \quad \mathbb{C}\_2^w = \min(A, 0.3) \,. \tag{69}$$

The damping coefficient appearing in the wall correction terms (27) and (28) is given by:

$$f\_w = \min\left[\frac{k^{3/2}}{2.5\varepsilon\alpha\_n}, 1.4\right].\tag{70}$$

The modelled dissipation rate transport equation is given by:

$$\begin{split} \frac{\mathrm{D}\varepsilon}{\mathrm{D}t} &= \mathsf{C}\_{\varepsilon 1} \frac{\varepsilon}{k} P\_{\mathrm{X}} - \mathsf{C}\_{\varepsilon 2} f\_{\varepsilon} \frac{\varepsilon \overline{\varepsilon}}{k} + \frac{\partial}{\partial \mathrm{x}\_{k}} \left( \mathsf{C}\_{\varepsilon} \frac{k}{\varepsilon} \overline{\mu\_{k} \mu\_{l}} \frac{\partial \varepsilon}{\partial \mathrm{x}\_{l}} + \nu \frac{\partial \overline{\varepsilon}}{\partial \mathrm{x}\_{k}} \right) \\ &+ \mathsf{C}\_{\varepsilon 3} \nu \frac{k}{\varepsilon} \overline{\mu\_{i} \mu\_{j}} \frac{\partial^{2} \mathrm{U}\_{k}}{\partial \mathrm{x}\_{i} \partial \mathrm{x}\_{l}} \frac{\partial^{2} \mathrm{U}\_{k}}{\partial \mathrm{x}\_{j} \partial \mathrm{x}\_{l}} - \mathsf{C}\_{\varepsilon 4} f\_{4} k \, \Omega\_{k} \, \Omega\_{k} + \mathrm{S}\_{l}. \end{split} \tag{71}$$

The coefficients have the following specified values:

$$\mathbb{C}\_{\varepsilon1} = 2.6 \quad \mathbb{C}\_{\varepsilon2} = 1.92 \quad \mathbb{C}\_{\varepsilon3} = 0.25 \quad \mathbb{C}\_{\varepsilon4} = 0.1 \quad \mathbb{C}\_{\varepsilon} = 0.18 \, \text{A} \tag{72}$$

and

$$f\_{\varepsilon} = 1 - \frac{\mathbf{C}\_{\varepsilon 2} - 1.4}{\mathbf{C}\_{\varepsilon 2}} \exp\left[-\left(\frac{\mathbf{R} \mathbf{e}\_{\mathrm{t}}}{6}\right)^{2}\right].\tag{73}$$

The length-scale growth correction, *Sl*, is given by:

$$S\_{l} = \max\left\{ \left[ \left( \frac{1}{\mathcal{C}\_{l}} \frac{\partial l}{\partial \mathbf{x}\_{n}} \right)^{2} - 1 \right] \left( \frac{1}{\mathcal{C}\_{l}} \frac{\partial l}{\partial \mathbf{x}\_{n}} \right)^{2}, 0 \right\} \frac{\overleftarrow{\varepsilon}\varepsilon}{k} A\_{l} \tag{74}$$

where *l* = *k*3/2/*ε*, and *Cl* = 2.5.

The anisotropic stress dissipation rate tensor is modelled as:

$$
\varepsilon\_{\rm ij} = f\_{\rm s} \varepsilon\_{\rm ij}^\* + (1 - f\_{\rm s}) \frac{2}{3} \delta\_{\rm ij} \varepsilon\_{\rm \prime} \tag{75}
$$

where *ε*∗ *ij* is given by:

$$\varepsilon\_{ij}^{\*} = \frac{\varepsilon}{k} \frac{\overline{u\_i u\_j} + \left(\overline{u\_i u\_k} \, n\_j \, n\_k + \overline{u\_j u\_k} \, n\_l \, n\_k + \overline{u\_k u\_l} \, n\_k \, n\_l \, n\_i \, n\_j\right) f\_d}{1 + \frac{3}{2} n\_p n\_q \frac{\overline{u\_p u\_q}}{k} f\_d} \, \tag{76}$$

$$f\_{\rm s} = 1 - \sqrt{A}E^2, \quad f\_d = (1 + 0.1 \text{Re}\_{\rm t})^{-1}. \tag{77}$$

**Low-Re TCL model**

Constructing *φ*∗

The quantities *φ*∗,*<sup>s</sup>*

*ij* , *<sup>φ</sup>*∗,*<sup>r</sup>*

respectively, but the coefficients *C*1, *C*<sup>2</sup> and *C*�

*C*�

as

where

and

A low-Re version of the TCL model was presented by Craft (1998). This version adopts a slightly different decomposition of the velocity-pressure gradient correlation Π*ij* (which appears in the exact RST equation *before* it is decomposed, as in (7)). The alternate decomposition was found to be more appropriate when modelling inhomogeneous flows. Where this correlation is typically decomposed into the pressure–strain-rate correlation and

Reynolds Stress Transport Modelling 19

pressure diffusion, an alternative decomposition is obtained by defining:

*φ*∗ *ij* <sup>=</sup> *<sup>φ</sup>*∗,*<sup>s</sup>*

*<sup>C</sup>*<sup>1</sup> <sup>=</sup> 3.1 *fA <sup>f</sup>*Ret *<sup>A</sup>*1/2

0.55 �

*f* �

*fA* ==

*Sij* <sup>=</sup> <sup>1</sup> 2 � *∂Ui ∂xj* + *∂Uj ∂xi*

<sup>2</sup> <sup>=</sup> min(0.6, *<sup>A</sup>*) + 3.5(*S*<sup>∗</sup> <sup>−</sup> <sup>Ω</sup>∗)

⎧ ⎪⎨

⎪⎩

*SI* <sup>=</sup> <sup>2</sup>

�

*<sup>C</sup>*<sup>2</sup> <sup>=</sup> min �

*φ*∗

*ij* <sup>=</sup> <sup>Π</sup>*ij* <sup>−</sup> <sup>1</sup>

*ij* <sup>+</sup> *<sup>φ</sup>*∗,*<sup>r</sup>*

<sup>1</sup> <sup>−</sup> exp �

*<sup>A</sup>* <sup>=</sup> <sup>√</sup>*A f*Ret <sup>+</sup> *<sup>A</sup>*(<sup>1</sup> <sup>−</sup> *<sup>f</sup>*Ret

(*A*/14)1/2 *A* < 0.05 *A*/0.8367 0.05 < *A* < 0.7 *A*1/2 *A* > 0.7 ,

<sup>√</sup>2*SijSjkSki*

, <sup>Ω</sup>*ij* <sup>=</sup> <sup>1</sup> 2 � *∂Ui ∂xj*

thus cannot contribute to the level of kinetic energy. This redistributive quantity is modelled

*ij* <sup>+</sup> *<sup>φ</sup>*inh,*<sup>s</sup>*

*ij* in this way ensures that it is redistributive in nature, since it is traceless and

*ij* <sup>+</sup> *<sup>φ</sup>*inh,*<sup>r</sup>*

<sup>2</sup> are prescribed by

<sup>−</sup>*A*1.5Ret

*ij* have the same form as their homogeneous counterparts (79) and (78),

<sup>100</sup> �� , 3.2*<sup>A</sup>*

<sup>2</sup> , (82a)

1 + *S*∗

<sup>3</sup> <sup>+</sup> *<sup>S</sup>*<sup>∗</sup> <sup>+</sup> <sup>Ω</sup><sup>∗</sup> <sup>−</sup> <sup>2</sup>*SI* , (82c)

*<sup>f</sup>*Ret = min[(Ret/160)2, 1] , (84)

*S*∗ = *Sk*/*ε*, Ω∗ = Ω*k*/*ε*, (86) *S* = (2*SijSji*)1/2, Ω = (2Ω*ij*Ω*ji*)1/2, (87)

(*SlmSml*)3/2 , (88)

�

. (89)

<sup>−</sup> *<sup>∂</sup>Uj ∂xi*

�

<sup>3</sup> *δij*Π*kk* . (80)

*ij* . (81)

), (83)

, (82b)

(85)

#### **3.10 The Two-Component-Limit model**

Researchers at UMIST, starting with the work of Fu et al. (1987), and Craft et al. (1989), developed a stress transport model that satisfies the constraint of realizability in the limit of two component turbulence. An outline of the derivation of the model is presented in Craft & Launder (2002). Using similar arguments as in (21), but retaining up to cubic terms in *aij*, and using the additional constraint of realizability, the following model for *φ<sup>r</sup> ij* was obtained

*φr ij* <sup>=</sup> <sup>−</sup> 0.6 *Pij* − 2/3*δijP* + 0.6*aijP* <sup>−</sup> 0.2 *ukuj ului k ∂Uk ∂xl* + *∂Ul ∂xk* <sup>−</sup> *uluk k uiuk ∂Uj ∂xl* + *ujuk ∂Ui ∂xl* − *c*<sup>2</sup> *A*2(*Pij* − *Dij*) + 3*amianj*(*Pmn* − *Dmn*) + *c*� 2 <sup>7</sup> <sup>15</sup> <sup>−</sup> *<sup>A</sup>*<sup>2</sup> 4 (*Pij* − 2/3*δijP*) + 0.2[*aij* − 1/2(*aik akj* − 1/3*δijA*2]*P* − 0.05*aijalkPkl* <sup>+</sup> 0.1 *uium <sup>k</sup> Pmj* <sup>+</sup> *ujum <sup>k</sup> Pmi* <sup>−</sup> 2/3 *<sup>δ</sup>ij ulum <sup>k</sup> Pml* <sup>+</sup> 0.1 *ului ukuj <sup>k</sup>*<sup>2</sup> <sup>−</sup> 1/3 *<sup>δ</sup>ij ulum ukum k*2 · 6*Dlk* + 13*k ∂Ul ∂xk* + *∂Uk ∂xl* + 0.2 *ului ukuj <sup>k</sup>*<sup>2</sup> (*Dlk* <sup>−</sup> *Plk*) (78)

where *A*<sup>2</sup> is the second invariant of the stress anisotropy tensor defined in (53). In the earliest, high-Re, version of the model the recommended values of the coefficients, *C*2, *C*� <sup>2</sup>, are

$$\mathbf{C}\_2 = \mathbf{0.55}, \quad \mathbf{C}\_2' = \mathbf{0.6} \dots$$

As for the slow pressure–strain-rate term, a second-order expression in *aij* is used, where the coefficients are allowed to depend on the stress anisotropy invariants in such a way as to satisfy realizability (Craft & Launder, 2002). Dependency on the third invariant, *A*3, is introduced through the flatness parameter, *A*, defined in (54). The flatness parameter becomes zero when one stress component vanishes; thus using the form

$$\phi\_{i\dot{j}}^s = -\mathsf{C}\_1 \varepsilon [a\_{i\dot{j}} + c\_1'(a\_{i\dot{j}}a\_{j\dot{k}} - \frac{1}{3}A\_2\delta\_{i\dot{j}})] - f\_A' \varepsilon a\_{i\dot{j}},\tag{79}$$

where the coefficients are given by

$$\mathbf{C}\_1 = 3.1(A\_2 A)^{1/2} \quad \mathbf{C}'\_1 = 1.1 \quad f'\_A = A^{1/2} \text{ .} $$

ensures that *φ<sup>s</sup> ij* drops to zero when the turbulence is two-component.

#### **Low-Re TCL model**

16 Numerical Simulations

Researchers at UMIST, starting with the work of Fu et al. (1987), and Craft et al. (1989), developed a stress transport model that satisfies the constraint of realizability in the limit of two component turbulence. An outline of the derivation of the model is presented in Craft & Launder (2002). Using similar arguments as in (21), but retaining up to cubic terms in *aij*,

*ij* was obtained

+ *ujuk*

*ulum <sup>k</sup> Pml*

 · 

<sup>2</sup> = 0.6 .

<sup>1</sup> = 1.1 *f* �

<sup>3</sup> *A*2*δij*)] − *f* �

*<sup>A</sup>* <sup>=</sup> *<sup>A</sup>*1/2 ,

6*Dlk* + 13*k*

 *∂Ul ∂xk* + *∂Uk ∂xl*

<sup>2</sup>, are

*<sup>A</sup>εaij* , (79)

(78)

*∂Ui ∂xl*

and using the additional constraint of realizability, the following model for *φ<sup>r</sup>*

+ 0.6*aijP*

*A*2(*Pij* − *Dij*) + 3*amianj*(*Pmn* − *Dmn*)

+ 0.2[*aij* − 1/2(*aik akj* − 1/3*δijA*2]*P* − 0.05*aijalkPkl*

*ujum*

*<sup>k</sup>*<sup>2</sup> <sup>−</sup> 1/3 *<sup>δ</sup>ij*

*<sup>k</sup>*<sup>2</sup> (*Dlk* <sup>−</sup> *Plk*)

zero when one stress component vanishes; thus using the form

*ij* = −*C*1*ε*[*aij* + *c*�

*C*<sup>1</sup> = 3.1(*A*2*A*)1/2 *C*�

*ij* drops to zero when the turbulence is two-component.

*φs*

where the coefficients are given by

ensures that *φ<sup>s</sup>*

(*Pij* − 2/3*δijP*)

*<sup>k</sup> Pmi* <sup>−</sup> 2/3 *<sup>δ</sup>ij*

high-Re, version of the model the recommended values of the coefficients, *C*2, *C*�

*C*<sup>2</sup> = 0.55, *C*�

*ulum ukum k*2

where *A*<sup>2</sup> is the second invariant of the stress anisotropy tensor defined in (53). In the earliest,

As for the slow pressure–strain-rate term, a second-order expression in *aij* is used, where the coefficients are allowed to depend on the stress anisotropy invariants in such a way as to satisfy realizability (Craft & Launder, 2002). Dependency on the third invariant, *A*3, is introduced through the flatness parameter, *A*, defined in (54). The flatness parameter becomes

<sup>1</sup>(*aijajk* <sup>−</sup> <sup>1</sup>

**3.10 The Two-Component-Limit model**

*φr*

*ij* = − 0.6

− 0.2

− *c*<sup>2</sup> 

+ *c*� 2 <sup>7</sup>

+ 0.1 *uium <sup>k</sup> Pmj* <sup>+</sup>

+ 0.1

+ 0.2

*Pij* − 2/3*δijP*

<sup>15</sup> <sup>−</sup> *<sup>A</sup>*<sup>2</sup> 4 

*ului ukuj*

*ului ukuj*

 *∂Uk ∂xl* + *∂Ul ∂xk* <sup>−</sup> *uluk k uiuk ∂Uj ∂xl*

 *ukuj ului k*

A low-Re version of the TCL model was presented by Craft (1998). This version adopts a slightly different decomposition of the velocity-pressure gradient correlation Π*ij* (which appears in the exact RST equation *before* it is decomposed, as in (7)). The alternate decomposition was found to be more appropriate when modelling inhomogeneous flows. Where this correlation is typically decomposed into the pressure–strain-rate correlation and pressure diffusion, an alternative decomposition is obtained by defining:

$$
\phi\_{\rm ij}^\* = \Pi\_{\rm ij} - \frac{1}{3} \delta\_{\rm ij} \Pi\_{kk} \,. \tag{80}
$$

Constructing *φ*∗ *ij* in this way ensures that it is redistributive in nature, since it is traceless and thus cannot contribute to the level of kinetic energy. This redistributive quantity is modelled as

$$
\phi\_{\rm ij}^{\*} = \phi\_{\rm ij}^{\*,s} + \phi\_{\rm ij}^{\*,r} + \phi\_{\rm ij}^{\rm inh,s} + \phi\_{\rm ij}^{\rm inh,r}.\tag{81}
$$

The quantities *φ*∗,*<sup>s</sup> ij* , *<sup>φ</sup>*∗,*<sup>r</sup> ij* have the same form as their homogeneous counterparts (79) and (78), respectively, but the coefficients *C*1, *C*<sup>2</sup> and *C*� <sup>2</sup> are prescribed by

$$\mathbf{C}\_1 = \mathbf{3}.1 f\_A f\_{\text{Re}} A\_2^{1/2} \,\text{,}\tag{82a}$$

$$\mathbf{C}\_{2} = \min\left\{ 0.55 \left[ 1 - \exp\left( \frac{-A^{1.5} \mathbf{Re}\_{\mathrm{t}}}{100} \right) \right], \frac{3.2A}{1 + \mathbf{S}^{\*}} \right\},\tag{82b}$$

$$\mathbf{C}'\_2 = \min(0.6, A) + \frac{3.5(\mathbf{S}^\* - \Omega^\*)}{3 + \mathbf{S}^\* + \Omega^\*} - 2\mathbf{S}\_{I, \prime} \tag{82c}$$

where

$$f\_A' = \sqrt{A}f\_{\text{Re}\_\text{t}} + A(1 - f\_{\text{Re}\_\text{t}}) \, \tag{83}$$

*<sup>f</sup>*Ret = min[(Ret/160)2, 1] , (84)

$$f\_A = = \begin{cases} (A/14)^{1/2} & A < 0.05\\ A/0.8367 & 0.05 < A < 0.7\\ A^{1/2} & A > 0.7 \end{cases} \tag{85}$$

$$\mathbb{S}^\* = \mathbb{S}k/\varepsilon, \quad \Omega^\* = \Omega k/\varepsilon,\tag{86}$$

$$\mathcal{S} = (2\mathcal{S}\_{\vec{i}\vec{j}}\mathcal{S}\_{\vec{j}\vec{i}})^{1/2}, \quad \Omega = (2\Omega\_{\vec{i}\vec{j}}\Omega\_{\vec{j}\vec{i}})^{1/2},\tag{87}$$

$$S\_I = \frac{2\sqrt{2}S\_{\bar{i}\bar{j}}S\_{\bar{j}k}S\_{\bar{k}i}}{(S\_{lm}S\_{ml})^{3/2}},\tag{88}$$

and

$$S\_{ij} = \frac{1}{2} \left( \frac{\partial \mathcal{U}\_i}{\partial \mathbf{x}\_j} + \frac{\partial \mathcal{U}\_j}{\partial \mathbf{x}\_i} \right), \Omega\_{ij} = \frac{1}{2} \left( \frac{\partial \mathcal{U}\_i}{\partial \mathbf{x}\_j} - \frac{\partial \mathcal{U}\_j}{\partial \mathbf{x}\_i} \right). \tag{89}$$

where

*ε*���

*ε* � *ij* =*ε*

*ε* �� *ij* =*ε* 2 *uluk <sup>k</sup> <sup>d</sup><sup>A</sup> <sup>l</sup> <sup>d</sup><sup>A</sup>*

*ε* ��� *ij* =*Cεsνk*

as discussed in Section 3.6. The term *ε*��

**Dissipation rate equation**

is solved, which takes the form

D*ε*˜ <sup>D</sup>*<sup>t</sup>* <sup>=</sup> *<sup>C</sup>ε*<sup>1</sup>

and is given by

and *F* in turn is given by

*F* = *∂l ∂xj*

*ε*˜ *k*

> *∂l ∂xj*

+ *Cε*<sup>3</sup> *ν k ε uiuj*

*P<sup>κ</sup>* − *Cε*<sup>2</sup>

*ε*˜ 2 *<sup>k</sup>* <sup>−</sup> *<sup>C</sup>*� *ε*2

*∂*<sup>2</sup>*Uk ∂xi∂xl*

*YE* = *Cε<sup>l</sup>*

*<sup>D</sup>* <sup>=</sup>*ε*�

*uiuj <sup>k</sup>* <sup>+</sup> <sup>2</sup>*<sup>ν</sup>*

> *∂* <sup>√</sup>*<sup>A</sup> ∂xk*

*kk* + *ε*��

without significant viscous effects (Craft & Launder, 1996).

+2*ν ului k ∂* √*k ∂xj ∂* √*k ∂xl*

*ulun k ∂* √*k ∂xl ∂* √*k ∂xn δij*

> *∂* <sup>√</sup>*<sup>A</sup> ∂xk*

*kk* + *ε*��� *kk* <sup>2</sup>*<sup>ε</sup>* ,

and the coefficients are taken as *f<sup>ε</sup>* = *A*3/2, *Cε<sup>s</sup>* = 0.2. The term *ε*�

et al., 1999; Craft, 1998), an equation for the homogeneous dissipation rate,

*ε*˜ = *ε* − 2*ν*

*∂*<sup>2</sup>*Uk ∂xj∂xl*

> *ε*˜ 2

+ 2*ν*

Reynolds Stress Transport Modelling 21

*<sup>k</sup> <sup>δ</sup>ij* <sup>−</sup> *ului*

*uluj k ∂* √*k ∂xi ∂* √*k ∂xl* ,

*<sup>k</sup> <sup>d</sup><sup>A</sup> <sup>l</sup> <sup>d</sup><sup>A</sup>*

*δij* + 2 *∂* <sup>√</sup>*<sup>A</sup> ∂xi*

to the model in (38), and its purpose is to ensure the correct wall-limiting behaviour of *εij*,

in *ε*<sup>12</sup> near *y*/*δ* = 0.1 observed in DNS studies of plane channel flow, and finally the term

*ij* improves the behaviour of *εij* at a free surface where there is strong inhomogeneity even

Early high-Re implementations of the TCL model used the same transport equation for the scalar dissipation rate (50) as in the LRR models. In later versions of the TCL model (Batten

> *∂k*1/2 *∂xi*

(*ε* − *ε*˜)*ε*˜ *<sup>k</sup>* <sup>+</sup>

+ *YE*.

The term *YE* is a length-scale correction based on the proposal of Iacovides & Raisee (1997),

<sup>2</sup>

*∂ ∂xk*  *Cε k ε ukul ∂ε*˜ *∂xl* + *ν ∂ε*˜ *∂xk*

*<sup>j</sup>* <sup>−</sup> *uluj <sup>k</sup> <sup>d</sup><sup>A</sup> <sup>l</sup> <sup>d</sup><sup>A</sup> i* ,

> ,

*ij* serves the specific purpose of producing the dip

*<sup>k</sup>* max[*F*(*<sup>F</sup>* <sup>+</sup> <sup>1</sup>)2, 0], (102)

− *Cl*{[1 − exp(−*Bε*Ret)] + *BεCl*Ret exp(−*Bε*Ret)} , (103)

*<sup>l</sup>* = *<sup>k</sup>*3/2/*ε*, *<sup>B</sup><sup>ε</sup>* = 0.1069, *Cl* = 2.55 . (104)

, (100)

(101)

(99)

*ij* is similar in nature

*∂* <sup>√</sup>*<sup>A</sup> ∂xj*

The inhomogeneous corrections are independent of the wall-normal vector, and are given by

*φ*inh,*<sup>s</sup> ij* =*fw*<sup>1</sup> *ε k* (*ulukd<sup>A</sup> <sup>l</sup> <sup>δ</sup>ij* <sup>−</sup> <sup>3</sup> <sup>2</sup> *uiukd<sup>A</sup> <sup>j</sup>* <sup>−</sup> <sup>3</sup> <sup>2</sup> *ujukd<sup>A</sup> <sup>i</sup>* )*d<sup>A</sup> k* +*fw*<sup>2</sup> *ε <sup>k</sup>*<sup>2</sup> *ulun*(*unukd<sup>A</sup> <sup>k</sup> <sup>δ</sup>ij* <sup>−</sup> <sup>3</sup> <sup>2</sup> *uiund<sup>A</sup> <sup>j</sup>* <sup>−</sup> <sup>3</sup> <sup>2</sup> *ujund<sup>A</sup> <sup>i</sup>* )*d<sup>A</sup> l* +*fw*3*ν ail ∂* √*k ∂xl ∂* √*k ∂xj* + *ajl ∂* √*k ∂xl ∂* √*k ∂xi* − 2 3 *anl ∂* √*k ∂xl ∂* √*k ∂xn <sup>δ</sup>ij* <sup>−</sup> <sup>4</sup> 3 *aij ∂* √*k ∂xl ∂* √*k ∂xl* +*f* � *w*1 *k*2 *ε ukul ∂* <sup>√</sup>*<sup>A</sup> ∂xk ∂* <sup>√</sup>*<sup>A</sup> ∂xl <sup>δ</sup>ij* <sup>−</sup> <sup>3</sup> 2 *uiuk ∂* <sup>√</sup>*<sup>A</sup> ∂xk ∂* <sup>√</sup>*<sup>A</sup> ∂xj* − 3 2 *ujuk ∂* <sup>√</sup>*<sup>A</sup> ∂xk ∂* <sup>√</sup>*<sup>A</sup> ∂xi* , (90) *φ*inh,r *ij* = *fIk ∂Ul ∂xn dldn*(*didj* <sup>−</sup> <sup>1</sup> <sup>3</sup> *dkdkδij*), (91)

where the 'normalised length-scale gradients', *di*, *d<sup>A</sup> <sup>i</sup>* , introduced by Craft & Launder (1996), are used indicate the direction of strong inhomogeneity, when present, without the use of a wall-normal vector. These are defined by

$$d\_{\bar{l}} = \frac{N\_{\bar{l}}}{0.5 + (N\_k N\_k)^{0.5}}, \quad \text{where} \quad N\_{\bar{l}} = \frac{\partial (k^{1.5}/\varepsilon)}{\partial x\_{\bar{l}}},\tag{92a}$$

$$d\_i^A = \frac{N\_i^A}{0.5 + (N\_k^A N\_k^A)^{0.5}}, \quad \text{where} \quad N\_i^A = \frac{\partial (k^{1.5} A^{0.5} / \varepsilon)}{\partial \mathbf{x}\_i}. \tag{92b}$$

The coefficients appearing in the inhomogeneous corrections are given by:

$$f\_{w1} = 0.4 + 1.6 \min\left\{1, \max\left[0, 1 - \frac{\text{Re}\_l - 55}{20}\right]\right\},\tag{93}$$

$$f\_{w2} = 0.1 + 0.8A\_2 \min\left\{1, \max\left[0, 1 - \frac{\text{Re}\_l - 50}{85}\right]\right\},\tag{94}$$

$$f\_{w3} = 2.5\sqrt{A} \,\text{\AA} \tag{95}$$

$$f\_{w1}' = 0.22\,\text{V} \tag{96}$$

$$f\_{\mathbf{I}} = 2.5 \, f\_{\mathbf{A}} \, . \tag{97}$$

As discussed in Section 3.6, the dissipation tensor near a wall or free surface is anisotropic, and the low-Re TCL accordingly prescribes the following anisotropic model for the dissipation rate tensor,

$$
\varepsilon\_{\rm ij} = (1 - f\_{\varepsilon}) \frac{\varepsilon\_{\rm ij}^{\prime} + \varepsilon\_{\rm ij}^{\prime\prime} + \varepsilon\_{\rm ij}^{\prime\prime\prime}}{D} + \frac{2}{3} f\_{\varepsilon} \varepsilon \delta\_{\rm ij} \,\,,\tag{98}
$$

where

18 Numerical Simulations

The inhomogeneous corrections are independent of the wall-normal vector, and are given by

<sup>2</sup> *ujukd<sup>A</sup>*

*<sup>j</sup>* <sup>−</sup> <sup>3</sup>

− 2 3 *anl ∂* √*k ∂xl ∂* √*k ∂xn*

*dldn*(*didj* <sup>−</sup> <sup>1</sup>

are used indicate the direction of strong inhomogeneity, when present, without the use of a

0.5 + (*NkNk*)0.5 , where *Ni* <sup>=</sup> *<sup>∂</sup>*(*k*1.5/*ε*)

*<sup>k</sup>* )0.5 , where *<sup>N</sup><sup>A</sup>*

 1, max 

> 1, max

*f* �

*ε*� *ij* + *ε*��

As discussed in Section 3.6, the dissipation tensor near a wall or free surface is anisotropic, and the low-Re TCL accordingly prescribes the following anisotropic model for the dissipation rate

> *ij* + *ε*��� *ij <sup>D</sup>* <sup>+</sup> <sup>2</sup>

*<sup>i</sup>* )*d<sup>A</sup> k*

<sup>2</sup> *ujund<sup>A</sup>*

*∂* <sup>√</sup>*<sup>A</sup> ∂xj*

*<sup>i</sup>* )*d<sup>A</sup> l*

> − 3 2 *ujuk ∂* <sup>√</sup>*<sup>A</sup> ∂xk*

> > *∂xi*

0, 1 <sup>−</sup> Ret <sup>−</sup> <sup>55</sup> 20

0, 1 <sup>−</sup> Ret <sup>−</sup> <sup>50</sup> 85

*<sup>i</sup>* <sup>=</sup> *<sup>∂</sup>*(*k*1.5*A*0.5/*ε*) *∂xi*

*fw*<sup>3</sup> <sup>=</sup> 2.5√*<sup>A</sup>* , (95)

*<sup>w</sup>*<sup>1</sup> = 0.22 , (96)

*fI* = 2.5 *fA* . (97)

<sup>3</sup> *fεεδij* , (98)

*<sup>δ</sup>ij* <sup>−</sup> <sup>4</sup> 3 *aij ∂* √*k ∂xl ∂* √*k ∂xl*

> *∂* <sup>√</sup>*<sup>A</sup> ∂xi*

<sup>3</sup> *dkdkδij*), (91)

*<sup>i</sup>* , introduced by Craft & Launder (1996),

(90)

 ,

, (92a)

. (92b)

, (93)

, (94)

*φ*inh,*<sup>s</sup> ij* =*fw*<sup>1</sup>

*ε k*

> *ail ∂* √*k ∂xl ∂* √*k ∂xj*

> > *ukul ∂* <sup>√</sup>*<sup>A</sup> ∂xk*

wall-normal vector. These are defined by

*dA*

tensor,

+*fw*<sup>2</sup> *ε*

+*fw*3*ν*

+*f* � *w*1 *k*2 *ε*

(*ulukd<sup>A</sup>*

*<sup>k</sup>*<sup>2</sup> *ulun*(*unukd<sup>A</sup>*

*<sup>l</sup> <sup>δ</sup>ij* <sup>−</sup> <sup>3</sup>

<sup>2</sup> *uiukd<sup>A</sup>*

*<sup>k</sup> <sup>δ</sup>ij* <sup>−</sup> <sup>3</sup>

+ *ajl ∂* √*k ∂xl ∂* √*k ∂xi*

*∂* <sup>√</sup>*<sup>A</sup> ∂xl*

*φ*inh,r *ij* = *fIk*

where the 'normalised length-scale gradients', *di*, *d<sup>A</sup>*

*di* <sup>=</sup> *Ni*

*<sup>i</sup>* <sup>=</sup> *<sup>N</sup><sup>A</sup>*

*i* 0.5 + (*N<sup>A</sup>*

*fw*<sup>1</sup> = 0.4 + 1.6 min

*fw*<sup>2</sup> = 0.1 + 0.8*A*<sup>2</sup> min

*εij* = (1 − *fε*)

*<sup>k</sup> <sup>N</sup><sup>A</sup>*

The coefficients appearing in the inhomogeneous corrections are given by:

*<sup>j</sup>* <sup>−</sup> <sup>3</sup>

<sup>2</sup> *uiund<sup>A</sup>*

*<sup>δ</sup>ij* <sup>−</sup> <sup>3</sup> 2 *uiuk ∂* <sup>√</sup>*<sup>A</sup> ∂xk*

> *∂Ul ∂xn*

$$\begin{split} \varepsilon'\_{ij} &= \varepsilon \frac{\overline{u\_i} \overline{u\_j}}{k} + 2\nu \frac{\overline{u\_l} \overline{u\_l}}{k} \frac{\partial \sqrt{k}}{\partial x\_l} \frac{\partial \sqrt{k}}{\partial x\_n} \delta\_{ij} \\ &+ 2\nu \frac{\overline{u\_l} \overline{u\_i}}{k} \frac{\partial \sqrt{k}}{\partial x\_j} \frac{\partial \sqrt{k}}{\partial x\_l} + 2\nu \frac{\overline{u\_l} \overline{u\_j}}{k} \frac{\partial \sqrt{k}}{\partial x\_i} \frac{\partial \sqrt{k}}{\partial x\_l} \end{split} \tag{99}$$
 
$$\begin{split} \varepsilon''\_{ij} &= \varepsilon \left( 2 \frac{\overline{u\_l} \overline{u\_k}}{k} d^A\_l d^A\_k \delta\_{ij} - \frac{\overline{u\_l} \overline{u\_i}}{k} d^A\_l d^A\_j - \frac{\overline{u\_l} \overline{u\_i}}{k} d^A\_l d^A\_i \right), \\ \varepsilon''\_{ij} &= \mathbb{C}\_{\varepsilon \varepsilon} \nu k \left( \frac{\partial \sqrt{A}}{\partial x\_k} \frac{\partial \sqrt{A}}{\partial x\_k} \delta\_{ij} + 2 \frac{\partial \sqrt{A}}{\partial x\_i} \frac{\partial \sqrt{A}}{\partial x\_j} \right), \\ D &= \frac{\varepsilon'\_{kk} + \varepsilon''\_{kk} + \varepsilon''\_{kk}}{2\varepsilon}. \end{split} \tag{99}$$

and the coefficients are taken as *f<sup>ε</sup>* = *A*3/2, *Cε<sup>s</sup>* = 0.2. The term *ε*� *ij* is similar in nature to the model in (38), and its purpose is to ensure the correct wall-limiting behaviour of *εij*, as discussed in Section 3.6. The term *ε*�� *ij* serves the specific purpose of producing the dip in *ε*<sup>12</sup> near *y*/*δ* = 0.1 observed in DNS studies of plane channel flow, and finally the term *ε*��� *ij* improves the behaviour of *εij* at a free surface where there is strong inhomogeneity even without significant viscous effects (Craft & Launder, 1996).

#### **Dissipation rate equation**

Early high-Re implementations of the TCL model used the same transport equation for the scalar dissipation rate (50) as in the LRR models. In later versions of the TCL model (Batten et al., 1999; Craft, 1998), an equation for the homogeneous dissipation rate,

$$\vec{\varepsilon} = \varepsilon - 2\nu \left( \frac{\partial k^{1/2}}{\partial x\_i} \right)^2 \,, \tag{100}$$

is solved, which takes the form

$$\begin{split} \frac{\mathbf{D}\tilde{\varepsilon}}{\mathbf{D}t} &= \mathsf{C}\_{\varepsilon 1} \frac{\tilde{\varepsilon}}{k} P\_{\mathbf{k}} - \mathsf{C}\_{\varepsilon 2} \frac{\tilde{\varepsilon}^2}{k} - \mathsf{C}\_{\varepsilon 2}' \frac{(\varepsilon - \tilde{\varepsilon})\tilde{\varepsilon}}{k} + \frac{\partial}{\partial \mathbf{x}\_{\hat{k}}} \left( \mathsf{C}\_{\varepsilon} \frac{k}{\varepsilon} \overline{\underline{u}\_{\hat{k}} \underline{u}\_{\hat{l}}} \frac{\partial \tilde{\varepsilon}}{\partial \mathbf{x}\_{\hat{l}}} + \nu \frac{\partial \tilde{\varepsilon}}{\partial \mathbf{x}\_{\hat{k}}} \right) \\ &+ \mathsf{C}\_{\varepsilon 3} \nu \frac{k}{\varepsilon} \overline{\underline{u}\_{\hat{l}} \underline{u}\_{\hat{l}}} \frac{\partial^2 \mathsf{U}\_{\hat{k}}}{\partial \mathbf{x}\_{\hat{l}} \partial \mathbf{x}\_{\hat{l}}} \frac{\partial^2 \mathsf{U}\_{\hat{k}}}{\partial \mathbf{x}\_{\hat{l}} \partial \mathbf{x}\_{\hat{l}}} + \mathsf{Y}\_{\mathcal{E}}. \end{split} \tag{101}$$

The term *YE* is a length-scale correction based on the proposal of Iacovides & Raisee (1997), and is given by

$$Y\_E = \mathcal{C}\_{\mathfrak{sl}} \frac{\mathfrak{z}^2}{k} \max[F(F+1)^2, 0],\tag{102}$$

and *F* in turn is given by

$$F = \left(\frac{\partial l}{\partial \mathbf{x}\_{\dot{j}}} \frac{\partial l}{\partial \mathbf{x}\_{\dot{j}}}\right) - \mathbb{C}\_{l} \{ [1 - \exp(-B\_{\ell} \mathbf{Re}\_{l})] + B\_{\ell} \mathbb{C}\_{l} \text{Re}\_{l} \exp(-B\_{\ell} \mathbf{Re}\_{l}) \},\tag{103}$$

$$l = k^{3/2} / \varepsilon, \quad B\_{\varepsilon} = 0.1069, \quad \mathcal{C}\_{l} = 2.55 \,\text{.} \tag{104}$$

Thus a suitable choice for *ν*<sup>11</sup> is

and relating the shear stress *uv* to *<sup>∂</sup><sup>U</sup>*

*u*2 *<sup>E</sup>* <sup>=</sup> *HE aE*

> <sup>2</sup> (*u*<sup>2</sup> *<sup>P</sup>* <sup>+</sup> *<sup>u</sup>*<sup>2</sup> *E*)

 linear interpolation

Using linear interpolation for *ν<sup>e</sup>*

+ 1 2Δ*x νP* <sup>11</sup> <sup>+</sup> *<sup>ν</sup><sup>E</sup>* 11

**5. Concluding remarks**

*u*2 *<sup>e</sup>* = <sup>1</sup>

Similarly,

terms.

discussed.

*u*2 *<sup>P</sup>* <sup>=</sup> <sup>1</sup> *aP* ∑ *i aiu*<sup>2</sup> *<sup>i</sup>* + *Su* + *Su aP*

> + *ν<sup>E</sup>* 11

*<sup>ν</sup>*<sup>11</sup> <sup>=</sup> <sup>2</sup> <sup>−</sup> <sup>4</sup>

*<sup>ν</sup>*<sup>22</sup> <sup>=</sup> <sup>2</sup> <sup>−</sup> <sup>4</sup>

*<sup>ν</sup>*<sup>12</sup> <sup>=</sup> <sup>1</sup> <sup>−</sup> *<sup>C</sup>*<sup>2</sup> *C*1

be accomplished through a Rhie-Chow-type interpolation (Leschziner & Lien, 2002):

 *HP*/*ap*

(*UP* <sup>−</sup> *UE*) <sup>−</sup> *<sup>ν</sup><sup>P</sup>*

 velocity smoothing

Similar expressions can be constructed for the remaining faces, and for the remaining stress

This chapter has provided an introduction to the subject of Reynolds stress transport modelling. A brief historical account of the development of this class of RANS models was presented. This was followed by an account of the theoretical background, assumptions, approximations, as well as the rationale behind the most commonly adopted RST modelling practises. Finally, some numerical implementation issues specific to RST models were briefly

The account served to illustrate areas of strength of this class of RANS models, such as the exact form of the stress production terms, and the abandoning of the incorrectly assumed direct link between stress and strain that characterises eddy–viscosity formulations. The presentation also serves to illustrate some inherent weaknesses of present RST models,

<sup>Δ</sup>*<sup>x</sup>* , *<sup>u</sup>*<sup>2</sup>

(*Uw* − *Ue*)*<sup>E</sup>*

Maintaining the required coupling between the velocity and Reynolds stress components can

Similar consideration of the *v*<sup>2</sup> transport equation leads to the specification

<sup>3</sup>*C*<sup>2</sup> *C*1

Reynolds Stress Transport Modelling 23

<sup>3</sup>*C*<sup>2</sup> *C*1

*k ε*

> *k ε*

*k ε*

*<sup>∂</sup><sup>y</sup>* leads to the following specification for *ν*<sup>12</sup>

+*ν<sup>P</sup>* 11

*<sup>e</sup>* <sup>=</sup> *He ae* + *ν<sup>e</sup>* 11

<sup>11</sup> and *He*/*ae*, one obtains for the value at face *e*:

<sup>11</sup>(*Uw* <sup>−</sup> *Ue*)*<sup>P</sup>* <sup>−</sup> *<sup>ν</sup><sup>E</sup>*

(*Uw* − *Ue*)*<sup>P</sup>*

*u*<sup>2</sup> . (110)

*v*2, (111)

*v*<sup>2</sup> . (112)

(*UP* − *UE*) <sup>Δ</sup>*<sup>x</sup>* .

<sup>11</sup>(*Uw* − *Ue*)*<sup>E</sup>*

. (114)

<sup>Δ</sup>*<sup>x</sup>* . (113)

The remaining coefficients are given by

$$\begin{aligned} \mathsf{C}\_{\varepsilon1} &= 1.0, \quad \mathsf{C}\_{\varepsilon2} = \frac{1.92}{1 + 0.7A\_d\sqrt{A\_2}}, \quad A\_d = \max(A, 0.25), \\ \mathsf{C}\_{\varepsilon2}' &= 1.0, \quad \mathsf{C}\_{\varepsilon3} = 0.875, \\ \mathsf{C}\_{\varepsilonl} &= 0.5, \quad \mathsf{C}\_{\varepsilon} = 0.15. \end{aligned} \tag{105}$$

#### **4. Numerical issues specific to RST modelling**

There are a number of numerical difficulties associated with the use of RST models that are not present when using eddy viscosity formulations. In particular, the use of RST models results in relatively large source terms that increase the stiffness of the algebraic equation system, in addition to the fact that the equation set becomes highly non-linear and strongly coupled (Leschziner & Lien, 2002; Lien & Leschziner, 1994). When using a collocated grid, there is also the issue of odd-even decoupling of the velocities and the Reynolds stresses.

The use of an eddy-viscosity approach adds to the momentum equations a momentum diffusion term that can be treated implicitly, thus enhancing stability. Since no such term is present in RST model equations, one approach to improve stability when applying RST models is to add and subtract a gradient-diffusion term based on an effective viscosity, *ν*eff. Considering the stress term *u*2, for example, one may write

$$
\overline{u^2} = \left(\overline{u^2} + \nu\_{\rm eff} \frac{\partial U}{\partial \mathbf{x}}\right) - \nu\_{\rm eff} \frac{\partial U}{\partial \mathbf{x}}\,\boldsymbol{\chi}
$$

allowing the unbracketed term to be treated implicitly in the *U*-momentum equation.

Since the effective viscosity does not affect the final converged solution, it is not uniquely specified. One would, in general, simply be trying to significantly reduce the residual stress term that must be treated explicitly in the source term. One way to specify the effective viscosity is by reference to a simplified form of the Basic Reynolds stress model equations. What is needed is to construct a relation between *u*<sup>2</sup> and *<sup>∂</sup><sup>U</sup> <sup>∂</sup><sup>x</sup>* , between *<sup>v</sup>*<sup>2</sup> and *<sup>∂</sup><sup>V</sup> <sup>∂</sup><sup>x</sup>* , and so on. Take *u*<sup>2</sup> for example, and start by assuming its transport equation is source dominated:

$$P\_{11} + \phi\_{11} - \frac{2}{3}\varepsilon\delta\_{ij} = 0 \,. \tag{106}$$

Substituting for *φ*<sup>11</sup> from the Basic model,

$$P\_{11} - C\_1 \varepsilon \left(\frac{\mu^2}{k} - \frac{2}{3}\right) - C\_2 \left[P\_{11} - \frac{1}{3}(P\_{11} + P\_{22} + P\_{33})\right] - \frac{2}{3} \varepsilon \delta\_{ij} = 0\,\tag{107}$$

This leads to

$$\begin{aligned} -2\overline{u^2}\frac{\partial U}{\partial \mathbf{x}}\left(1 - \frac{2}{3}\mathbf{C}\_2\right) - \mathbf{C}\_1\frac{\varepsilon}{k}\overline{u^2} + \\ \left(\text{other terms not containing } \overline{u^2} \text{ or } \frac{\partial U}{\partial \mathbf{x}}\right) &= 0, \quad \text{(108)} \end{aligned}$$

or

$$
\overline{u^2} = \frac{(2 - \frac{4}{3}C\_2)\overline{u^2}}{C\_1}\frac{k}{\varepsilon}\frac{\partial U}{\partial x} + \text{O.T.}\tag{109}
$$

20 Numerical Simulations

There are a number of numerical difficulties associated with the use of RST models that are not present when using eddy viscosity formulations. In particular, the use of RST models results in relatively large source terms that increase the stiffness of the algebraic equation system, in addition to the fact that the equation set becomes highly non-linear and strongly coupled (Leschziner & Lien, 2002; Lien & Leschziner, 1994). When using a collocated grid, there is also

The use of an eddy-viscosity approach adds to the momentum equations a momentum diffusion term that can be treated implicitly, thus enhancing stability. Since no such term is present in RST model equations, one approach to improve stability when applying RST models is to add and subtract a gradient-diffusion term based on an effective viscosity, *ν*eff.

> *∂U ∂x* − *ν*eff

Since the effective viscosity does not affect the final converged solution, it is not uniquely specified. One would, in general, simply be trying to significantly reduce the residual stress term that must be treated explicitly in the source term. One way to specify the effective viscosity is by reference to a simplified form of the Basic Reynolds stress model equations.

3

(*P*<sup>11</sup> + *P*<sup>22</sup> + *P*33)

other terms not containing *u*<sup>2</sup> or

*∂U ∂x* ,

*<sup>∂</sup><sup>x</sup>* , between *<sup>v</sup>*<sup>2</sup> and *<sup>∂</sup><sup>V</sup>*

*εδij* = 0 . (106)

 − 2 3 *<sup>∂</sup><sup>x</sup>* , and so on.

*εδij* = 0 . (107)

= 0 , (108)

*∂U ∂x* 

+ O.T. (109)

*u*<sup>2</sup> + *ν*eff

allowing the unbracketed term to be treated implicitly in the *U*-momentum equation.

Take *u*<sup>2</sup> for example, and start by assuming its transport equation is source dominated:

*<sup>P</sup>*<sup>11</sup> <sup>+</sup> *<sup>φ</sup>*<sup>11</sup> <sup>−</sup> <sup>2</sup>

3*C*2)*u*<sup>2</sup> *C*1

*k ε ∂U ∂x*

*<sup>u</sup>*<sup>2</sup> <sup>=</sup> (<sup>2</sup> <sup>−</sup> <sup>4</sup>

<sup>√</sup>*A*<sup>2</sup>

, *Ad* = max(*A*, 0.25),

(105)

1 + 0.7*Ad*

the issue of odd-even decoupling of the velocities and the Reynolds stresses.

The remaining coefficients are given by

*C*�

*<sup>C</sup>ε*<sup>1</sup> <sup>=</sup> 1.0, *<sup>C</sup>ε*<sup>2</sup> <sup>=</sup> 1.92

*<sup>ε</sup>*<sup>2</sup> = 1.0, *Cε*<sup>3</sup> = 0.875, *Cε<sup>l</sup>* = 0.5, *C<sup>ε</sup>* = 0.15.

**4. Numerical issues specific to RST modelling**

Considering the stress term *u*2, for example, one may write

What is needed is to construct a relation between *u*<sup>2</sup> and *<sup>∂</sup><sup>U</sup>*

Substituting for *φ*<sup>11</sup> from the Basic model,

 *u*<sup>2</sup> *<sup>k</sup>* <sup>−</sup> <sup>2</sup> 3 − *C*<sup>2</sup> *<sup>P</sup>*<sup>11</sup> <sup>−</sup> <sup>1</sup> 3

*P*<sup>11</sup> − *C*1*ε*

This leads to

<sup>−</sup> <sup>2</sup>*u*<sup>2</sup> *<sup>∂</sup><sup>U</sup> ∂x* <sup>1</sup> <sup>−</sup> <sup>2</sup> 3 *C*2 − *C*<sup>1</sup> *ε k u*2+

or

*u*<sup>2</sup> =  Thus a suitable choice for *ν*<sup>11</sup> is

$$\nu\_{11} = \frac{2 - \frac{4}{3}\mathcal{C}\_2}{\mathcal{C}\_1} \frac{k}{\varepsilon} \overline{u^2}. \tag{110}$$

Similar consideration of the *v*<sup>2</sup> transport equation leads to the specification

$$\nu\_{22} = \frac{2 - \frac{4}{3} \mathcal{C}\_2}{\mathcal{C}\_1} \frac{k}{\varepsilon} \overline{v^2} \,, \tag{111}$$

and relating the shear stress *uv* to *<sup>∂</sup><sup>U</sup> <sup>∂</sup><sup>y</sup>* leads to the following specification for *ν*<sup>12</sup>

$$\nu\_{12} = \frac{1 - \mathcal{C}\_2}{\mathcal{C}\_1} \frac{k}{\varepsilon} \overline{v^2}. \tag{112}$$

Maintaining the required coupling between the velocity and Reynolds stress components can be accomplished through a Rhie-Chow-type interpolation (Leschziner & Lien, 2002):

$$\overline{u\_P^2} = \underbrace{\frac{1}{a\_P} \left(\sum\_i a\_i \overline{u\_i^2} + S\_{\mathcal{U}}\right) + \frac{S\_{\mathcal{U}}}{a\_P}}\_{\text{H}p/a\_P} + \nu\_{11}^P \frac{(\mathcal{U}\_{\mathcal{U}} - \mathcal{U}\_{\mathcal{E}})\_P}{\Delta \mathbf{x}}.\tag{113}$$

Similarly,

$$
\overline{u\_E^2} = \frac{H\_E}{a\_E} + \nu\_{11}^E \frac{(\mathcal{U}\_w - \mathcal{U}\_\varepsilon)\_E}{\Delta \mathfrak{x}}, \qquad \overline{u\_\varepsilon^2} = \frac{H\_\varepsilon}{a\_\varepsilon} + \nu\_{11}^\varepsilon \frac{(\mathcal{U}\_P - \mathcal{U}\_E)}{\Delta \mathfrak{x}}.
$$

Using linear interpolation for *ν<sup>e</sup>* <sup>11</sup> and *He*/*ae*, one obtains for the value at face *e*:

$$\begin{split} \overline{u\_{\varepsilon}^{2}} &= \underbrace{\frac{1}{2} (\overline{u\_{P}^{2}} + \overline{u\_{E}^{2}})}\_{\text{linear interpolation}} \\ &+ \underbrace{\frac{1}{2\Delta x} \left\{ \left[ \nu\_{11}^{P} + \nu\_{11}^{E} \right] (\mathcal{U}\_{P} - \mathcal{U}\_{E}) - \nu\_{11}^{P} (\mathcal{U}\_{W} - \mathcal{U}\_{\varepsilon})\_{P} - \nu\_{11}^{E} (\mathcal{U}\_{W} - \mathcal{U}\_{\varepsilon})\_{E} \right\}}\_{\text{velocity smoothing}}. \end{split} \tag{114}$$

Similar expressions can be constructed for the remaining faces, and for the remaining stress terms.

#### **5. Concluding remarks**

This chapter has provided an introduction to the subject of Reynolds stress transport modelling. A brief historical account of the development of this class of RANS models was presented. This was followed by an account of the theoretical background, assumptions, approximations, as well as the rationale behind the most commonly adopted RST modelling practises. Finally, some numerical implementation issues specific to RST models were briefly discussed.

The account served to illustrate areas of strength of this class of RANS models, such as the exact form of the stress production terms, and the abandoning of the incorrectly assumed direct link between stress and strain that characterises eddy–viscosity formulations. The presentation also serves to illustrate some inherent weaknesses of present RST models,

Donaldson, C. d. (1971). A progress report on an attempt to construct an invariant model of

Reynolds Stress Transport Modelling 25

Fu, S. (1988). *Computational modelling of turbulent swirling flows with second-moment closures*,

Fu, S., Launder, B. E. & Tselepidkis, D. P. (1987). Accomodating the effect of high strain rates in

Gatski, T. B. (2004). Constitutive equations for turbuelnt flows, *Theoret. Comput. Fluid Dynamics*

Gibson, M. M. & Launder, B. E. (1978). Ground effects on pressure fluctuations in the atmospheric boundary layer, *Journal of Fluid Mechanics* 86(03): 491–511. Hanjali´c, K. & Jakirli´c, S. (2002). Second-moment turbulence closure modelling, *in* B. Launder

Hanjali´c, K. & Launder, B. E. (1972). A reynolds stress model of turbulence and it's application

Hanjali´c, K. & Launder, B. E. (1976). Contribution towards a reynolds-stress closure for low-reynolds number turbulence, *Journal of Fluid Mechanics* 74: 593–610. Iacovides, H. & Raisee, M. (1997). Computation of flow and heat transfer in 2d rib roughened

Jakirli´c, S. & Hanjali´c, K. (1995). A second-moment closure for non-equilibrium and separating

Kassinos, S. C. & Reynolds, W. C. (1994). A structure-based model for the rapid distortion

Launder, B. E. & Reynolds, W. C. (1983). Asymptotic near-wall stress dissipation rates in a

Launder, B., Reece, G. & Rodi, W. (1975). Progress in the development of a reynolds-stress

Launder, B. & Shima, N. (1989). Second-moment closure for the near-wall sublayer:

Leschziner, M. A. & Lien, F. S. (2002). Numerical aspects of applying second-moment closure

Lien, F. S. & Leschziner, M. a. (1994). Upstream monotonic interpolation for scalar transport

Lumley, J. (1967). Rational approach to relations between motions of differing scales in

Lumley, J. & Khajeh-Nouri, B. (1974). Computational modelling of turbulent transport, *Adv.*

to complex flows, *in* B. Launder & N. Sandham (eds), *Closure Strategies for Turbulent*

with application to complex turbulent flows, *International Journal for Numerical*

PhD thesis, Faculty of Technology, University of Manchester.

to thin shear flows, *Journal of Fluid Mechanics* 52(4): 609–638.

*Symposium on Turbulence, Heat and Mass Transfer*, Delft.

closures, *Journal of Fluid Mechanics* 269: 143–168.

turbulent flow, *Physics of Fluids* 26(5): 1157–1158.

*and Transitional Flows*, Cambridge University Press.

turbulent flows, *Physics of Fluids* 10: 1405–1408.

*Methods in Fluids* 19: 527–548.

*Geophys.* 18Aa: 169–192.

turbulence closure, *Journal of Fluid Mechanics* 68(3): 537–566.

development and appilication, *AIAA Journal* 27: 1319–1325.

London.

Eng., UMIST.

18: 345–369.

University Press.

University.

turbulent shear flows, *Proc. AGARD Conf. on Turbulent Shear Flows*, number 1 in *B*,

modelling the pressure-strain correlation, *Technical Report TFD/87/5*, Dept. of Mech.

& N. Sandham (eds), *Closure Strategies for Turbulent and Transitional Flows*, Cambridge

passges, *in* K. Hanjali´c & T. Peeters (eds), *Proceedings of the Second International*

high- and low-re number flows, *Proceedings of the* 10th *Symposium on Turbulent Shear Flows*, Vol. 3, The Pennsylvania State University, Pennsylvania, pp. 23–25–23–30. Johansson, A. V. & Hallbäck, M. (1994). Modelling of rapid pressure–strain in reynolds-stress

of homogeneous turbulence, *Technical Report TF-61*, Mech. Eng. Dept., Stanford

which might also be thought of as areas for potential improvement. These weaknesses are a natural result of the complexity of turbulent phenomena, and of the persistent *closure* problem—transport equations for any level of statistical moments will always contain unclosed higher moment terms.

The realizability constraint is ultimately a kinematic constraint that serves to prevent certain un-physical results. Aside from that, it does not prescribe any particular dynamic *stimulus–response* type of link between the strain field and inter-component redistribution processes. Therefore there is no reason to expect that redistributive models, in the form of tensor polynomial expansions in stress anisotropy and velocity gradient, satisfying such constraints should return the correct response to all possible strain fields and histories, particularly ones far removed from those for which the models were calibrated.

This does not diminish the value of RST models, but rather serves to emphasise the importance of testing and validation in order to understand the limits of validity and accuracy for intended applications. As discussed earlier, there is always a trade-off between accuracy and computational cost, and the need for reliable RANS models for many types engineering simulations is not likely to be replaced by LES or DNS in the near future. More importantly, these arguments emphasise the need to strive for a deeper and more general understanding of the complex turbulent phenomena described by the unclosed terms in the transport equations, with the aim of building better models.

#### **6. References**


22 Numerical Simulations

which might also be thought of as areas for potential improvement. These weaknesses are a natural result of the complexity of turbulent phenomena, and of the persistent *closure* problem—transport equations for any level of statistical moments will always contain

The realizability constraint is ultimately a kinematic constraint that serves to prevent certain un-physical results. Aside from that, it does not prescribe any particular dynamic *stimulus–response* type of link between the strain field and inter-component redistribution processes. Therefore there is no reason to expect that redistributive models, in the form of tensor polynomial expansions in stress anisotropy and velocity gradient, satisfying such constraints should return the correct response to all possible strain fields and histories,

This does not diminish the value of RST models, but rather serves to emphasise the importance of testing and validation in order to understand the limits of validity and accuracy for intended applications. As discussed earlier, there is always a trade-off between accuracy and computational cost, and the need for reliable RANS models for many types engineering simulations is not likely to be replaced by LES or DNS in the near future. More importantly, these arguments emphasise the need to strive for a deeper and more general understanding of the complex turbulent phenomena described by the unclosed terms in the transport equations,

Batten, P., Craft, T. J., Leschziner, M. A. & Loyau, H. (1999). Reynolds-stress-transport

Chou, P. Y. (1945). On velocity correlations and the solutions of the equations of turbulence

Chung, M. K. & Kim, S. K. (1995). A nonlinear return-to-isotropy model with reynolds number

Craft, T. (1991). *Second-moment modelling of turbulent scalar transport*, PhD thesis, Faculty of

Craft, T. J. (1998). Developments in a low-reynolds-number second moment closure and it's

Craft, T. J., Fu, S., Launder, B. E. & Tselepidakis, D. P. (1989). Developments in modelling the

Craft, T. J. & Launder, B. E. (1996). A reynolds stress closure designed for complex geometries,

Craft, T. J. & Launder, B. E. (2002). Closure modelling near the two-component limit, *in*

Crow, S. (1968). Visco-elastic properties of fine-grained incompressible turbulence, *Journal of*

Crow, S. C. (1967). Visco-elastic character of fine-grained isotropic turbulence, *Physics of Fluids*

Daly, B. J. & Harlow, F. H. (1970). Transport equations in turbulence, *Physics of Fluids*

application to separating and reattaching flows, *International Journal of Heat and Fluid*

turbulent second-moment pressure correlations, *Technical Report TFD/89/1*, Dept. of

B. Launder & N. Sandham (eds), *Closure Strategies for Turbulent and Transitional Flows*,

and anisotropy dependency, *Physics of Fluids* 7: 1425–1436.

*International Journal of Heat and Fluid Flow* 17(3): 245–254.

modelling for compressible aerodynamic applications, *AIAA Journal* 37(7): 785–797.

particularly ones far removed from those for which the models were calibrated.

unclosed higher moment terms.

with the aim of building better models.

*Flow* 19: 541–548.

Mech. Eng., UMIST.

Cambridge University Press.

*Fluid Mechanics* 33: 1–20.

10: 1587–1589.

13: 2634–2649.

fluctuation, *Quart. Appl. Math.* 3: 38–54.

Technolgy, University of Manchester.

**6. References**


**2** 

Xiaolong Yang *Hunan University* 

*China* 

**Study of Some Key Issues for Applying LES** 

Most of nature and industry flows are turbulence. There are three kinds of numerical simulation methods for turbulent flows (Lesieur 1990; Pope 2000; Sagaut 2000, 2006): direct numerical simulation (DNS), Reynolds-averaged Navier-Stokes equations (RANS) and large eddy simulation (LES). DNS is a straightforward way to simulate turbulent flows. Full Navier-Stokes equations are discretized and solved numerically without any model, empirical parameter or approximation. Theoretically speaking, results of DNS exactly reflect the real flow and the whole range of turbulence scales are computed. With DNS, people can compute and visualize any quantity of interest, including some that are too difficult or impossible to be measured by experiments. But as we all know the computation cost is very high. For high Reynolds number flow, even modern computer technology can not satisfy the

In RANS, the flow quantities are decomposed into two parts: the average or mean term and the fluctuating term by applying Reynolds averaging. The effect of the fluctuating quantities on the mean flow quantities is described by the so called Reynolds stress tensor, which is must be modelled in terms of the mean velocities. Typical models can be grouped loosely into three categories: algebraic models, one-equation models and two-equation models. RANS is simple and robust. It is widely used in engineering problem. The general limitation of RANS is the fact that the model must represent a very wide range of scales. While the small scales tend to be universal, and depend on viscosity, the larger scales depend largely on flow condition and boundaries. So there is no one universal model for all flows. For different flows, the model must be modified to obtain good results. Another issue is that usually a time averaging is adopted in RANS. So RANS has difficult to handle unsteady

In LES, a filter is applied to separate the large scales from small scales. Then only the large, energy carrying scales (or called resolved scales) of turbulence are computed exactly by solving the governing equations. While the small, fluctuating scales are modelled, which is also called subgrid scales (SGS). Compared to RANS, LES has several advantages: 1) LES can capture the large scales directly which are the main energy container of turbulence and response for the momentum and energy transfer. 2) The dissipation of turbulence energy is believed to be done by small scales. Since small scales are thought to be homogenous,

**1. Introduction** 

computation requirement.

flows.

**to Real Engineering Problems** 

