# 计算流体力学(CFD)文档——10. Advanced turbulence modelling

10. ADVANCED TURBULENCE MODELLING 1. Turbulence models for general-purpose CFD 2. Linear eddy-viscosity models 3. Non-linear eddy-viscosity models 4. Differential stress models 5. Implementation of turbulence models in CFD References Examples

SPRING 2011

1. Turbulence Models For General-Purpose CFD Turbulence models for general-purpose CFD must be frame-invariant – i.e. independent of any particular coordinate system – and hence must be expressed in tensor form. This rules out simpler models of boundary-layer type. Turbulent flows are computed either by solving the Reynolds-averaged Navier-Stokes equations with suitable models for turbulent fluxes or by computing the fluctuating quantities directly. The main types are summarised below.

Reynolds-Averaged Navier-Stokes (RANS) Models Linear eddy-viscosity models (EVM) – (deviatoric) turbulent stress proportional to mean strain; – eddy viscosity based on turbulence scalars (usually k + one other), determined by solving transport equations. Non-linear eddy-viscosity models (NLEVM) – turbulent stress is a non-linear function of mean strain and vorticity; – coefficients depend on turbulence scalars (usually k + one other), determined by solving transport equations; – mimic response of turbulence to certain important types of strain. Differential stress models (DSM) – aka Reynolds-stress transport models (RSTM) or second-order closure (SOC); – solve (modelled) transport equations for all turbulent fluxes.

Models That Compute Fluctuating Quantities Large-eddy simulation (LES) – compute time-varying fluctuations, but model sub-grid-scale motion. Direct numerical simulation (DNS) – no modelling; resolve the smallest scales of the flow.

1

David Apsley

2. Linear Eddy-Viscosity Models 2.1 General Form Stress-strain constitutive relation: U U j 2 k ui u j = t i + x xi 3 j

ij

,

t

=

t

(1)

The eddy viscosity t is derived from turbulent quantities such as the turbulent kinetic energy k and its dissipation rate . These quantities are themselves determined by solving scalartransport equations (see below). A typical shear stress and normal stress are given by U V uv = t y + x U 2 u2 = 2 t k x 3 From these the other stress components are easily deduced by inspection/cyclic permutation. Most eddy-viscosity models in general-purpose CFD codes are of the 2-equation type; (i.e. scalar-transport equations are solved for 2 turbulent scales). The commonest types are k- and k- models, for which specifications are given below.

2.2 k- Models Eddy viscosity:
t

Scalar-transport equations (non-conservative form): k Dk + (P (k ) ) = ( (k ) ) xi xi Dt D ( () ) + (C 1 P ( k ) C 2 ) = k Dt xi xi rate of change diffusion production dissipation + advection The diffusivities of k and are related to the eddy-viscosity:
(k )

The rate of production of turbulent kinetic energy (per unit mass) is U i P ( k ) ≡ u i u j x j In the standard k- model (Launder and Spalding, 1974) the coefficients take the values (k) ( ) C 2 = 1.44, C = 0.09, C 1 = 1.92, = 1.0, = 1.3
¤ ¤ ¤

2

t (k )

,

( )

=

+

=

+

t ( )

=C

k2

(2)

(3)

(4)

(5)

David Apsley

Other important variants include RNG k- (Yakhot et al., 1992) and low-Re models such as Launder and Sharma (1974), Lam and Bremhorst (1981), and Lien and Leschziner (1993). Modifications are employed in low-Re models (see later) to incorporate effects of molecular viscosity. Specifically, C , C 1 and C 2 are multiplied by viscosity-dependent factors f , f1 and f2 respectively, and an additional source term S( ) may be required in the equation. The damping factor f is necessary because t ∝ y3 as y → 0, but k ∝ y2 and ~ constant, so that k2/ yields the wrong power of y. Some models (notably Launder and Sharma, 1974) solve for the homogeneous dissipation rate ~ which vanishes at solid boundaries and is related to by =~+D, D = 2 (k 1 / 2 ) 2 (6) 2 This reflects the theoretical near-wall behaviour of (i.e. 2 k / y ) in a form which avoids using a geometric distance y explicitly.
¤ ¤ ¤

2.3 k-

Models

) is sometimes known as the specific dissipation rate and has C k dimensions of frequency or (time)–1. (nominally equal to

Eddy viscosity: k t = Scalar-transport equations: k Dk = ( (k ) ) + ( P (k ) xi Dt xi
( )

(7)

*

k)

Again, the diffusivities of k and
(k )

The original k- model was that of Wilcox (1988a) where the coefficients take the values 9 5 3 (k) ( ) * = 2.0, = 2.0 (9) = , = , = , 100 9 40 (see Wilcox, but in later versions of the model the coefficients become functions of k 2 / 1998). Menter (1994) devised a shear-stress-transport (SST) model. The model blends the kmodel (which is – allegedly – superior in the near-wall region), with the k- model (which is less sensitive to the level of turbulence in the free stream) using wall-distance-dependent blending functions. Transport equations are solved for k and , but this is an odd choice Advanced turbulence modelling 3 David Apsley

t (k )

,

( )

=

+

D ( = Dt xi

)+ ( xi

(8)

2

P
t

(k )

)

are related to the eddy-viscosity:
= +
t ( )

because in a free flow with no wall boundaries (e.g. a jet) the model is simply a transformed k- model. All models of k- type suffer from a problematic wall boundary condition ( which is routinely fudged! as y → 0)

2.4 Behaviour of Linear Eddy-Viscosity Models in Simple Shear In simple shear flow the shear stress is U uv = t y The three normal stresses are predicted to be equal: 2 u 2 = v 2 = w2 = k 3 whereas, in practice, there is considerable anisotropy; e.g. in the log-law region: u 2 : v 2 : w 2 ≈ 1.0 : 0.4 : 0.6
y

v

u

U(y)

Actually, in simple shear flows, this is not a problem, since only the gradient of the shear stress uv plays a dynamically-significant role in the mean momentum balance. However, it is a warning of more serious problems in complex flows.

4

David Apsley

3. Non-Linear Eddy-Viscosity Models 3.1 General Form The stress-strain relationship for linear eddy-viscosity models gives for the deviatoric Reynolds stress (i.e. subtracting the trace): U U j 2 u i u j k ij = t i + x 3 xi j 2 Dividing by k and writing t = C k / gives
(10) k We define the LHS of (10) as the anisotropy tensor aij, the dimensionless and traceless form of the Reynolds stress: ui u j 2 aij ≡ (11) ij k 3 For the RHS of (10), the symmetric and antisymmetric parts of the mean-velocity gradient are called the mean strain and mean vorticity tensors, respectively: U j 1 U 1 U i U j S ij ≡ ( i + (12) ), ( ) ij ≡ xi xi 2 x j 2 x j Their dimensionless forms, scaled by the turbulent timescale k/ , are written in lower case: k k sij ≡ S ij , (13) ij ≡ ij
ui u j
ij

Then equation (10) can be written a ij = 2C s ij or, (14) a = 2C s Hence, for a linear eddy-viscosity model the anisotropy tensor is proportional to the dimensionless mean strain. The main idea of non-linear eddy-viscosity models is to generalise this to a non-linear relationship between the anisotropy tensor and the mean strain and vorticity: a = 2C s + NL(s, ) (15) Additional non-linear components cannot be completely arbitrary, but must be symmetric and traceless. For example, a quadratic NLEVM must be of the form a = 2C s (16) + 1 (s 2 1 {s 2 }I) + 2 ( s s ) + 3 ( 2 1 { 2 }I) 3 3 where {.} denotes a trace and I is the identity matrix; i.e. {M} ≡ trace(M) ≡ M ii , (I) ij ≡ ij (17) We shall see below that an appropriate choice of the coefficients 1, 2 and 3 allows the model to reproduce the correct anisotropy in simple shear.

Theory (based on the Cayley-Hamilton Theorem – a matrix satisfies its own characteristic equation) predicts that the most general possible relationship involves ten independent tensor

2 3

= C

k U i U j + x xi j

5

David Apsley

bases and includes terms up to the 5th power in s and
a = ∑ C T (s, )
=1 10

: (18)

where all T are linearly-independent, symmetric, traceless, second-rank tensor products of s and and the C are, in general, functions of their invariants. One possible choice of bases (but by no means the only one) is

T1 = s Linear: Quadratic: T2 = s 2 1 {s 2 }I 3 T3 = s s
T4 =
2

Cubic: Quartic:

Quintic:

Exercise. (i) Prove that all these bases are symmetric and traceless. (ii) Show that the bases T5 – T10 vanish in 2-d incompressible flow. The first base corresponds to a linear eddy-viscosity model and the next three to the quadratic extension in equation (16). T5, T7, T8, T9 clearly contain multiples of earlier bases and hence could be replaced by simpler forms; however, the bases chosen here ensure that they vanish in 2-d incompressible flow. A number of routes have been taken in devising such NLEVMs, including: assuming the form of the series expansion to quadratic or cubic order and then calibrating against important flows (e.g. Speziale, 1987; Craft, Launder and Suga, 1996); simplifying a differential stress model by an explicit solution (e.g. Speziale and Gatski, 1993) or by successive approximation (e.g. Apsley and Leschziner, 1998); renormalisation group methods (e.g. Rubinstein and Barton, 1990); direct interaction approximation (e.g. Yoshizawa, 1987). In devising such NLEVMs, model developers have sought to incorporate such physicallysignificant properties as realisability: u2 ≥ 0

u u

2

≤ u2 u2

1{ 3

2

}I

T5 = T6 = T7 = T9 = T10 =

2

s+s s2 + s2
2 2

2

{
2

2

}s 2 { s }I 3
2

s2 s2
2

2 {s 2 3 1{ 2 s
2

}I {

2

}( s 2 1 {s 2 }I) 3

T8 = s 2 s s s s
2

s 2 1 {s 2 }( s s ) 2
2

s
2

2

}( s s )

(positive normal stresses) (Cauchy Schwartz inequality)

(19)

6

David Apsley

3.2 Cubic Eddy-Viscosity Models The preferred level of modelling in MACE is a cubic eddy viscosity model, which can be written in the form a = 2C f s + 1 (s 2 1 {s 2 }I) + 2 ( s s ) + 3 ( 2 1 { 2 }I) 3 3 1{s 2 }s 2 { 2 }s 3 ( 2 s + s 2 { 2 }s 2 { s }I) 4 ( s 2 s 2 ) 3 (20) Note the following properties (some of which will be developed further below and in the Examples). (i) A cubic stress-strain relationship is the minimum order with at least the same number of independent coefficients as the anisotropy tensor (i.e. 5). In this case it will be precisely 5 if we assume 3 = 0 (see (vi) below) and note that the 1 and 2 terms are tensorially similar to the linear term (see (iv) below). (ii) The first term on the RHS corresponds to a linear eddy-viscosity model. (iii) The various non-linear terms evoke sensitivities to specific types of strain: – the quadratic ( 1, 2, 3) terms evoke sensitivity to anisotropy; – the cubic 1 and 2 terms evoke sensitivity to curvature; – the cubic 4 term evokes sensitivity to swirl. (iv) The 1 and 2 terms are tensorially proportional to the linear term; however they (or rather their difference) provide a sensitivity to curvature, so have been kept distinct. (v) The 3 and 4 terms vanish in 2-d incompressible flow. (vi) Theory and experiment indicate that pure rotation generates no turbulence. This implies that 3 should either be zero or tend to 0 in the limit S → 0 . As an example of such a model we cite the Craft et al. (1996) model in which coefficients are functions of the mean-strain invariants and turbulent Reynolds number: 0.3[1 exp(0.36e 0.75 )] C = 1 + 0.35 3 / 2 (21) Rt 1 / 2 Rt 2 k2 f = 1 exp[( ) ( Rt = ~ ) ], 90 400 where k S = 2 S ij S ij , = 2 ij ij , = ~ max(S , ) (22) The coefficients of non-linear terms are (in the present notation): ( 1 , 2 , 3 ) = (0.4, 0.4, 1.04)C f ( 1 , 2 , 3 , 4 ) = (40, 40, 0, 80)C 3 f (23) Non-linearity is built into both tensor products and strain-dependent coefficients – notably C . The model is completed by transport equations for k and ~ where ~ = 2 (k 1 / 2 ) 2 ~ is called the homogeneous dissipation rate, which vanishes at solid boundaries because of the near-wall behaviour of k and (see later): k ~ k0 y 2 , ~ 2 k0 ~ 2 k / y 2 In this model, mean strain and vorticity are also non-dimensionalised using ~ rather than .

7

David Apsley

3.3 General Properties of Non-Linear Eddy-Viscosity Models (i) 2-d Incompressible Flow The non-linear combinations of s and incompressible flow. In such a flow: 0 s11 s12 0 = 21 s = s 21 s 22 0 , 0 0 0 0 have particularly simple forms in 2-d
12

0 0

0 0 0
ij

Incompressibility ( s11 = s 22 ) and the symmetry and antisymmetry properties of sij and ( s 21 = s12 , 21 = 12 ) reduce these to 0 0 s11 s12 0 12 = 12 0 0 s = s12 s11 0 , 0 0 0 0 0 0 From these we find
1 0 0 s = ( s + s ) 0 1 0 , 0 0 0 1 0 0 s s = 2 12 s12 0 1 0 2 0 0 0
2 2 11 2 12

1 0 0 = 0 1 0 0 0 0 0 1 0 0 0 12 s11 1 0 0 0
2 2 12

(24)

PROPERTY 1 In 2-d incompressible flow: 2 2 s 2 = ( s11 + s12 )I 2 = 1 {s 2 }I 2 2 (25) 2 2 1 = 12I 2 = 2 { 2 }I 2 where, I2 = diag(1,1,0). In particular, taking tensor products of s2 or 2 with matrices whose third row and third column are all zero has the same effect as multiplication by the scalars 2 2 1 1 } respectively. 2 {s } or 2 {

PROPERTY 2 P (k ) = aij s ij = {as}

(26)

Moreover, in 2-d incompressible flow the quadratic terms do not contribute to the production of turbulent kinetic energy. Proof. U i P ( k ) = u i u j = k (aij + 2 ij )( S ij + ij ) 3 x j

8

David Apsley

Now (aij +

2 3

ij

)

ij

= 0 since the first factor is symmetric whilst the
ij

ij

is antisymmetric.

Also, incompressibility implies P ( k ) = kaij S ij or P (k ) = aij s ij = {as}

S ij = S ii = 0 . Hence,

This is true for any incompressible flow, but, in the 2-d case, multiplying (20) by s, taking the trace and using the results (25) it is found that the contribution of the quadratic terms to {as} is 0.

PROPERTY 3 In 2-d incompressible flow the 3- and 4-related terms of the non-linear expansion (20) vanish. Proof. Substitute the results (25) for s2 and 2 into (20).

(ii) Particular Types of Strain The non-linear constitutive relationship (20) allows the model to mimic the response of turbulence to particular important types of strain.
PROPERTY 4 The quadratic terms yield turbulence anisotropy in simple shear:

u2 2 = +( 3 k v2 2 = +( k 3 w2 2 = ( k 3

1 +6

2

2

2 3)

12
2

1 6

2

3)

12

where

=

k U y

(27)

1

3)

6
450

This may be deduced by substituting the results (24) into (20), noting that s11 = 0, whilst 1 k U 1 = s12 = 12 = 2 y 2
y+

400

350

300

250

uu+ vv+ ww+ -uv+

As an example the figure right shows application of the Apsley and Leschziner (1998) model to computing the Reynolds stresses in channel flow.

200

150

100

50

0 0 1 2 3 4 5 6 7 8

9

David Apsley

PROPERTY 5 The 1 and 2-related cubic terms yield the correct sensitivity to curvature. U U U s V = = s , where Rc is radius of curvature. From (24), , In curved shear flow, Rc y R x

{s 2 } + { where s12 = Hence,

2

2 } ≡ 2( s12

2 12

) = 1 U s U s + 2 R Rc

1 U s U s , Rc 2 R
2

12

k U s U s } ≡ 2( ) 2 R Rc Inspection of the production terms in the stress-transport equations (Section 4) shows that curvature is stabilising (reducing turbulence) if Us increases in the direction away from the centre of curvature (Us/R > 0) and destabilising 'stable' curvature (increasing turbulence) if Us (reducing turbulence) decreases in the direction away from the centre of curvature (Us/R < 0). In the constitutive relation (20) the response is correct if positive. {s 2 } + {

'unstable' curvature (increasing turbulence)

1

and

2

are both

PROPERTY 6 In 3-d flows, the swirl.

4-related

term evokes the correct sensitivity to

W

U

10

David Apsley

4. Differential Stress Modelling
Differential stress models (aka Reynolds-stress transport models or second-order closure) solve separate scalar-transport equations for each stress component. By writing a momentum equation for the fluctuating ui component, multiplying by uj, adding a corresponding term with i and j reversed and then adding we obtain: 1 (u i u j ) + U k (u i u j ) = [ (u i u j ) p (u i jk + u j ik ) u i u j u k ] t x k x k x k
(u i u k U j x k + u j uk u i u j U i p u u j ) + (u i f j + u j f i ) + ( i + )2 x j xi x k x k x k

or, multiplying by : D(u i u j ) Dt = d ijk x k + ( Pij + Fij +
ij

ij

)

(28)

For a full derivation see the Appendix. Term Name and role Exact or Model EXACT ( ui u j ) + ( U k ui u j ) t x k

D(u i u j ) RATE OF CHANGE (time derivative + advection) Dt Transport with the mean flow.
Pij

PRODUCTION (mean strain) Generation by mean-velocity gradients. PRODUCTION (body forces) Generation by body forces. DIFFUSION Spatial redistribution.

EXACT
Pij ≡ u i u k

U j x k

u j uk

U i x k

Fij

EXACT (in principle) Fij ≡ u i f j + u j f i MODEL
k u k ul (u i u j ) xl

d ijk

d ijk = (

kl

+ Cs
( 2) ij

)

MODEL (1) ij = ij +
(1) ij

+

( w) ij

Φ ij

PRESSURE-STRAIN Redistribution between components.

= C1

( 2) ij
( w) ij

ij ) k = C 2 ( Pij 1 Pkk ij ) 3 ~ ~ ~ = ( kl nk nl ij 3 ik n j nk 3 jk ni nk ) f 2 2

(u i u j 2 k 3

ij

=C

( w) 1

k

ui u j + C

( w) 2

( 2) ij

,

f =

ε ij

DISSIPATION Removal by viscosity

MODEL 2 ij = 3 ij

11

David Apsley

~

C 3/ 4k 3/ 2/ yn

Typical values of the constants are: C1 = 1.8 , C 2 = 0.6 , C1( w ) = 0.5 ,

( C 2w ) = 0.3

(29)

The stress-transport equations must be supplemented by a means of specifying – typically by its own transport equation, or a transport equation for a related quantity such as . The most significant term requiring modelling is the pressure-strain correlation. This term is traceless (the sum of the diagonal terms 11 + 22 + 33 = 0 ) and its accepted role is to promote isotropy – hence models for
(1) ij

and

( 2) ij

. Near walls this isotropising tendency
(w ) ij

must be over-ridden, necessitating a "wall-correction" term

.

General Assessment of DSMs For:

Include more turbulence physics than eddy-viscosity models. Advection and production terms (the "energy-in" terms) are exact and do not need modelling.

Against: Models are very complex and many important terms (particularly the redistribution and dissipation terms) require modelling. Models are very expensive computationally (6 stress-transport equations in 3 dimensions) and tend to be numerically unstable (only the small molecular viscosity contributes to any sort of gradient diffusion term).

DSMs of Interest The nearest thing to a standard model is a high-Re closure based on that of Launder et al. (1975), with wall-reflection terms from Gibson and Launder (1978). Other models of interest include: Speziale et al. (1991) – non-linear ij formulation, eliminating wall-correction terms; Craft (1998) – low-Re DSM, attempting to eliminate wall-dependent parameters; Jakirli and Hanjali (1995) – low-Re DSM admitting anisotropic dissipation; Wilcox (1988b) – low-Re DSM, with rather than as additional turbulent scalar.

12

David Apsley

5. Implementation of Turbulence Models in CFD 5.1 Transport Equations
The implementation of a turbulence model in CFD requires: (1) a means of specifying the turbulent stresses by either: – a constitutive relation (eddy-viscosity models), or – individual transport equations (differential stress models); (2) the solution of additional scalar-transport equations.

Special Considerations for the Mean Flow Equations Only part of the stress is diffusive. u i u j represents a turbulent flux of Ui-momentum in the xj direction. For linear or non-linear eddy-viscosity models only a part of this can be treated implicitly as a diffusion-like term; e.g. for the U equation through a face normal to the y direction: U V uv = t ( + ) + (non linear terms ) y 14444244443 x 123
diffusive part transferred to source

The remainder of the flux is treated as part of the source term for that control volume. Nevertheless, it is still treated in a conservative fashion; i.e. the mean momentum lost by one cell is equal to that gained by the adjacent cell. The lack of a turbulent viscosity in differential stress models leads to numerical instability. This can be addressed by the use of "effective viscosities" – see below.

Special Considerations for the Turbulence Equations They are usually source-dominated; i.e. the most significant terms are production, redistribution and dissipation (this is sometimes invoked as an excuse to use a loworder advection scheme). Variables such as k and must be non-negative. This demands: – care in discretising the source term (see below); – use of an unconditionally-bounded advection scheme.

Source-Term Linearisation For Non-Negative Quantities The general discretised scalar-transport equation for a control volume centred on node P is a P φ P ∑ a F φ F = bP + s P φ P
F

For stability one requires sP ≤ 0 To ensure non-negative φ one requires, in addition,

13

David Apsley

bP ≥ 0 You should, by inspection of the k and source term is linearised in this way.

transport equations (3), be able to identify how the

If bP < 0 for a quantity such as k or which is always non-negative (e.g. due to transfer of non-linear parts of the advection term or non-diffusive fluxes to the source term) then, to ensure that the variable doesn't become negative, the source term should be rearranged as bP sP → sP + * (30) φP
bP → 0 where * denotes the current value of a variable.

5.2 Wall Boundary Conditions

At walls the no-slip boundary condition applies, so that both mean and fluctuating velocities vanish. At high Reynolds numbers this presents three problems: there are very large flow gradients; wall-normal fluctuations are suppressed (i.e., selectively damped); viscous and turbulent stresses are of comparable magnitude. There are two main ways of handling this in turbulent flow: low-Reynolds-number turbulence models – resolve the flow right up to the wall with a very fine grid and viscous modifications to the turbulence equations; wall functions – use a coarser grid and assume theoretical profiles in the unresolved near-wall region.

5.2.1 Low-Reynolds-Number Turbulence Models

Aim to resolve the flow right up to the boundary. Have to include effects of molecular viscosity in the coefficients of the eddy-viscosity formula and (or ) transport equations. Try to ensure the theoretical near-wall behaviour (see the Examples): 2 k 3 ~ 2 ~ constant , k ∝ y2, ( y → 0) t ∝ y y Full resolution of the flow requires the near-wall node to satisfy y+ ≤ 1, where u y y+ ≡ , u = w/ This can be very computationally demanding, particularly for high-speed flows.

(31)

(32)

14

David Apsley

5.2.2 High-Reynolds-Number Turbulence Models

Bridge the near-wall region with wall functions; i.e., assume profiles (based on boundary-layer theory) between near-wall node and boundary. OK if the equilibrium assumption is reasonable (e.g. slowly-developing boundary layers), but dodgy in highly non-equilibrium regions (particularly near impingement, separation or reattachment points).

control volume

near-wall node

Up
assumed velocity profile

w (wall shear stress)

τ

The near-wall node should ideally (if seldom in practice) be placed in the log-law region (30 < y+ < 150). This means that numerical meshes cannot be arbitrarily refined close to solid boundaries.

In the finite-volume method, various quantities are required from the wall-function approach. Values may be fixed on the wall (w) itself or by forcing a value at the near-wall node (P).
Variable Mean velocity (U,V,W) Wall boundary condition (relative) velocity =0 at the wall Required from wall function Wall shear stress

k, u i u j

u i u j = k = 0 at the wall; zero flux
P

Cell-averaged production and dissipation Value at the near-wall node

fixed at near-wall node

The means of deriving these quantities are set out below.

Mean-Velocity Equation: Wall Shear Stress
Method 1

The friction velocity u is defined in terms of the wall shear stress: 2 w = u If the near-wall node lies in the logarithmic region then y u UP 1 + + = ln( Ey P ) , (33) yP = P u where subscript P denotes the near-wall node. Given the value of UP this could be solved (iteratively) for u and hence the wall stress w. This works well only for near-equilibrium boundary layers.

Method 2

A better approach when the turbulence is far from equilibrium (e.g. near separation or reattachment points) is to estimate an "equivalent" friction velocity from the turbulent kinetic Advanced turbulence modelling 15 David Apsley

energy at the near-wall node: 1 u 0 = C 1 / 4 k P/ 2 and integrate the mean-velocity profile assuming an eddy viscosity t. If we adopt the log-law version: t = u0 y and solve for U from U w = t y we get u 0U P (34) w/ = u0 y P ln( E )

(If the turbulence were genuinely in equilibrium, then u0 would equal u and (33) and (34) would be equivalent).

Method 3

An even more advanced approach which depends less on the assumption that the near-wall node is well within the log layer is to assume a total viscosity (molecular + eddy) which matches both the viscous ( eff = ) and log-layer ( eff u y) limits: + max{0, u 0 ( y y )} (35) eff = where y is a matching height. Similar integration to that above (see the examples) leads to a complete mean-velocity profile satisfying both linear viscous sublayer and log-law limits: y+, y+ ≤ y+ yu U w = 2 × + 1 y+ ≡ 0 (36) + + + + , y ≥y u0 u 0 y + ln{1 + ( y y )}, where we note that y+ is based on u0 rather than the unknown u . A similar approach can be applied for rough-wall boundary layers (Apsley, 2007). A typical value of the nondimensional matching height (for smooth walls) is y + = 7.37 . As far as the computational implementation is concerned the required output for a finitevolume calculation is the wall shear stress in terms of the mean velocity at the near-wall node, yp, not vice versa. To this end, (36) is conveniently rearranged in terms of an effective wall viscosity eff,wall such that Up (37) w = eff , wall yp where + 1, yP ≤ y + + yP + = × (38) yP ≥ y + , eff , wall 1 y + + ln{1 + ( y + y + )} P (The reason for using the form (37) is that the code will "see" the mean velocity gradient as Up/yp.)

16

David Apsley

k Equation: Cell-Averaged Production and Dissipation

Because, in a finite-volume calculation, it represents the total production and dissipation for the near wall cell the source term of the k transport equation requires cell-averaged values of production P(k) and dissipation rate . These are derived by integrating assumed profiles for these quantities: y≤y 0 U (k ) U 2 where t = eff – P ≡ uv = (39) y>y ( ) y t y
(y ≤ y ) w 3 u0 (40) = (y > y ) (y y ) d For smooth walls, the matching height y and offset yd are given in wall units by (see Apsley, 2007): + y + = 27.4 , y d = 4 .9

Integration over a cell (see the examples) leads to cell averages ( / )2 ( + y+ ) 1 ( ln[1 + ( + y + )] Pavk ) ≡ P ( k ) dy = w u0 1+ ( + y+ ) 0
av
¤

Equation: Boundary Condition on is fixed from its assumed profile (equation (40)) at the near-wall node. A particular value at a cell centre can be forced in a finite-volume calculation by modifying the source coefficients: sP → , bP → P where is a large number (e.g. 1030). The matrix equations for that cell then become ( + a P )φ P ∑ a F φ F = P or ∑ aF φ F + φP = P + aP + aP Since is a large number this effectively forces φP to take the value P.
P

Reynolds-Stress Equations For the Reynolds stresses, one method is to fix the values at the near-wall node from the nearwall value of k and the structure functions u i u j /k , the latter being derived from the differential stress-transport equations on the assumption of local equilibrium. For the standard model this gives (see the examples):

17

¤

¤

1 ≡ 0

3 yd u0 y )+ dy = ln( y yd y yd

¤

(41) (42)

David Apsley

( v 2 2 1 + C1 + C 2 2C 2 w ) C 2 = k 3 C1 + 2C1( w) ( u 2 2 2 + C1 2C 2 + C 2 w ) C 2 C1( w) v 2 + = C k k C1 3 1 ( w 2 2 1 + C1 + C 2 + C 2 w ) C 2 C1( w ) v 2 + = C k k C1 3 1

(43)

( 1 C 2 + 3 C 2w) C 2 v 2 uv 2 = C + 3 C ( w) k k 1 2 1

With the values for C1, C2, etc. from the standard model this gives u2 = 1.098 , k v2 = 0.248 , k w2 = 0.654 , k uv = 0.255 k (44)

When the near-wall flow and wall-normal direction are not conveniently aligned in the x and y directions respectively, the actual structure functions can be obtained by rotation. However, for 3-dimensional and separating/reattaching flow the flow-oriented coordinate system is not fixed a priori and can swing round significantly between iterations if the mean velocity is small, making convergence difficult to obtain. A second – and now my preferred – approach (Apsley, 2007) is to use cell-averaged production and dissipation in the Reynolds-stress equation in the same manner as the k-equation.

5.3 Effective Viscosity for Differential Stress Models
DSMs contain no turbulent viscosity and have a reputation for numerical instability. An artificial means of promoting stability is to add and subtract a gradient-diffusion term to the turbulent flux: U U u u = (u u + ) (45) x x with the first part averaged between nodal values and the last part discretised across a cell face and treated implicitly. (This is analogous to the Rhie-Chow algorithm for pressurevelocity coupling in the momentum equations).

The simplest choice for the effective viscosity k2 = t =C
¤

A better choice is to make use of a natural linkage between individual stresses and the corresponding mean-velocity gradient which arise from the actual stress-transport equations. Assuming that the stress-transport equations (with no body forces) are source-dominated then Pij + ij ij ≈ 0 or, with the basic DSM (without wall-reflection terms),

is just (46)

18

David Apsley

Pij C1 (

ui u j k

2 3

ij

) C 2 ( Pij 1 Pkk 3

ij

) 2 3

ij

≈0
¤

Expand this, identifying the terms which contain only u u or
¤

For the normal stresses u 2 :
k

P C1

(u 2 K) C 2 2 P + K = 0 3

Hence, u2 =

(1 2 C 2 ) k (1 2 C 2 ) k U 3 3 P +K = (2u 2 + K) C1 C1 x

Similarly for the shear stresses u u :
k

P C1

u u C2 P + K = 0

whence

u u =

U (1 C 2 ) k (1 C 2 ) k P +K = (u 2 + K) C1 C1 x

Hence, from the stress-transport equations, U +K = u2 x U +K u u = x where the effective viscosities (both for the Uα component of momentum) are:

1 2 C2 k u 2 , =2 3 C 1

1 C2 k u 2 = C 1

Note that the effective viscosities are anisotropic, being linked to particular normal stresses. A more detailed analysis can accommodate wall-reflection terms in the pressure-strain model, but the extra complexity is not justified.

19

U x

as follows.

¤

¤

¤

¤¤

¤

¤

¤

(47)

(48)

David Apsley

References
Apsley, D.D., 2007, CFD calculation of turbulent flow with arbitrary wall roughness, Flow, Turbulence and Combustion, 78, 153-175. Apsley, D.D. and Leschziner, M.A., 1998, A new low-Reynolds-number nonlinear twoequation turbulence model for complex flows, Int. J. Heat Fluid Flow, 19, 209-222 Craft, T.J., 1998, Developments in a low-Reynolds-number second-moment closure and its application to separating and reattaching flows, Int. J. Heat Fluid Flow, 19, 541-548. Craft, T.J., Launder, B.E. and Suga, K., 1996, Development and application of a cubic eddyviscosity model of turbulence, Int. J. Heat Fluid Flow, 17, 108-115. Gatski, T.B. and Speziale, C.G., 1993, On explicit algebraic stress models for complex turbulent flows, J. Fluid Mech., 254, 59-78. Gibson, M.M. and Launder, B.E., 1978, Ground effects on pressure fluctuations in the atmospheric boundary layer, J. Fluid Mech., 86, 491-511. Hanjali , K., 1994, Advanced turbulence closure models: a view of current status and future prospects, Int. J. Heat Fluid Flow, 15, 178-203. Jakirli , S. and Hanjali , K., 1995, A second-moment closure for non-equilibrium and separating high- and low-Re-number flows, Proc. 10th Symp. Turbulent Shear Flows, Pennsylvania State University. Lam, C.K.G. and Bremhorst, K.A., 1981, Modified form of the k-e model for predicting wall turbulence, Journal of Fluids Engineering, 103, 456-460. Launder, B.E., 1989, Second-Moment Closure and its use in modelling turbulent industrial flows, Int. J. Numer. Meth. Fluids, 9, 963-985. Launder, B.E., Reece, G.J. and Rodi, W., 1975, Progress in the development of a Reynoldsstress turbulence closure, J. Fluid Mech., 68, 537-566. Launder, B.E. and Sharma, B.I., 1974, Application of the energy-dissipation model of turbulence to the calculation of flow near a spinning disc, Letters in Heat and Mass Transfer, 1, 131-138. Launder, B.E. and Spalding, D.B., 1974, The numerical computation of turbulent flows, Computer Meth. Appl. Mech. Eng., 3, 269-289. Lien, F-S. and Leschziner, M.A., 1993, Second-moment modelling of recirculating flow with a non-orthogonal collocated finite-volume algorithm, in Turbulent Shear Flows 8 (Munich, 1991), Springer-Verlag. Menter, F.R., 1994, Two-equation eddy-viscosity turbulence models for engineering applications, AIAA J., 32, 1598-1605. Rubinstein, R. and Barton, J.M., 1990, Nonlinear Reynolds stress models and the renormalisation group, Phys. Fluids A, 2, 1472-1476. Speziale, C.G., 1987, On nonlinear K-l and K- models of turbulence, J. Fluid Mech., 178, 459-475. Speziale, C.G., Sarkar, S. and Gatski, T.B., 1991, Modelling the pressure-strain correlation of turbulence: an invariant dynamical systems approach, J. Fluid Mech., 227, 245-272. Wilcox, D.C., 1988, Reassessment of the scale-determining equation for advanced turbulence models, AIAA J., 26, 1299-1310. Wilcox, D.C., 1988, Multi-scale model for turbulent flows, AIAA Journal, 26, 1311-1320. Wilcox, D.C., 1998, Turbulence modelling for CFD, 2nd Edition, DCW Industries. Yakhot, V., Orszag, S.A., Thangam, S., Gatski, T.B. and Speziale, C.G., 1992, Development of turbulence models for shear flows by a double expansion technique, Phys. Fluids A, 7, 1510. Yoshizawa, A., 1987, Statistical analysis of the derivation of the Reynolds stress from its eddy-viscosity representation, Phys. Fluids, 27, 1377-1387.

20

David Apsley

Appendix: Derivation of the Reynolds-Stress Transport Equations (For Reference and Hangover-Cure Only!)
Restrict to constant-density fluids for simplicity. fi are the components of problem-dependent body forces (buoyancy, Coriolis forces, ...). Use prime/overbar notation (e.g. u = u + u ′ ) and summation convention throughout. Continuity Instantaneous: Average: Subtract: u k =0 x k u k =0 x k ′ u k =0 x k (A1) (A2) (A3)

Result (I): Both mean and fluctuating quantities satisfy the incompressibility equation.

Momentum Instantaneous: Average: u i 2ui u 1 p + + uk i = + fi xi t x k x k x k (A4)

u i 2ui u u ′ 1 p ′ + + uk i + uk i = + fi (A5) xi t x k x k x k x k Du i u i 1 p (A6) Rearrange: ( u i′u ′ ) + f i + = k Dt xi x k x k ′ (Note: u k can be "taken through" the /xk derivative whenever required, due to the incompressibility condition)

Result (II): Mean flow equation is the same as the instantaneous equation except for additional apparent stresses u i′u ′j . Note: this is actually the net transport of momentum by turbulent fluctuations, not a force.
2 u i′ u u i′ u i′ u ′ u ′ 1 p ′ ′ i + uk i uk i = ′ ′ + f i′ + + uk + uk x k x k xi x k t x k x k x k

Subtract: Similarly j:

(A7)i (A7)j

u ′j

2 u ′j 1 p ′ ′ + + u′ + uk uk + f j′ = + u′ k k x j x k t x k x k x k x k x k u ′j u j u ′j

u ′j

Form u i′ (A7) j + u ′j (A7) i :

21

David Apsley

u j u ′ ′ (u i′u ′j ) + u k (u i′u ′j ) + u i′u ′ (u i′u ′j u k ) + u ′j u k i + k x k x k x k x k t 2 u ′j 2 u i′ p ′ p ′ ) + (u i′ ) + (u i′ f j′ + u ′j f i′) + u ′j + u ′j = (u i′ x k x k x k x k xi x j Rewrite the pressure terms and rearrange: 1 ′ (u i′u ′j ) + u k (u i′u ′j ) = [ (u i′u ′j ) p ′(u i′ jk + u ′j ik ) u i′u ′j u k ] t x k x k x k 1
(u i′u ′ k u j x k + u ′j u ′ k u i u i′ u ′j p ′ u ′ u ′j ) + (u i′ f j′ + u ′j f i′) + ( i + )2 x k x j xi x k x k

Result (III): Reynolds-stress transport equation: d ijk D (u i′u ′j ) = + Pij + Fij + ij Dt x k where: D (u i′u ′j ) = (u i′u ′j ) + u k (u i′u ′j ) Dt t x k u j u + u ′j u ′ i ) Pij = (u i′u ′ k k x k x k Fij = u i′ f j′ + u ′j f i′
ij

ij

(A8)

advection (by mean flow) production (by mean strain) production (by body forces) pressure-strain correlation dissipation

= =2

p ′ u i′ u ′j + ( ) x j xi u i′ u ′j x k x k (u i′u ′j ) x k 1 p ′(u i′
jk

ij

d ijk =

+ u ′j

ik

′ ) u i′u ′j u k

diffusion

Contract (A8), then divide by 2. Change summation subscript k to i to minimise confusion. Result (IV): turbulent kinetic energy equation: Dk d i( k ) = + P (k ) + F (k ) Dt xi where: u P ( k ) = u i′u ′j i production (by mean strain) x j F ( k ) = u i′ f i′ = (
d i( k )

(A9)

production (by body forces) dissipation diffusion

u i′ 2 ) x j k 1 p ′u i′ 1 u ′j u ′j u i′ = 2 xi

22

David Apsley

Examples
RANS Models - General Q1. (a) (b)

What is meant by a "Reynolds-averaged Navier-Stokes" solver? Describe the main principles of: (i) eddy-viscosity models; (ii) non-linear eddy-viscosity models; (iii) differential stress models; (iv) large-eddy simulation.

Q2. in terms of the fundamental physical Write down the dimensions of , t, k, and dimensions M (mass), L (length) and T (time). Hence, show that any expression for t in terms of , k and must be of the form k2 = constant × t whilst any expression for
t

t

in terms of , k and

must be of the form

= constant ×

k

Summation Convention Q3. U i U 1 U 2 U 3 U V W + + . Expand the following. or + is shorthand for + x y x1 xi x 2 x3 z U i U i DU i (a) (= +Uk ) when i = 1 and when i = 2. Dt t x k U i P ( k ) = u i u j (b) x j (c) (d) U i ) when i = 1, j = 1 and when i = 1, j = 2. x k x k For a general matrix M, what quantities are represented by Mii, MijMji and MijMjkMki? Pij = (u i u k + u j uk U j

Q4. What do the following quantities reduce to in a simple shear flow (in which U/y is the only non-zero velocity component)? (a) P(k) (use the definition in Question 3(b)) (b) Pij (for each combination of i and j; use the definition in Question 3(c)) (c) S = (2S ij S ij )1 / 2 (d) = (2
ij ij

)1 / 2

23

David Apsley

Production of Turbulent Kinetic Energy Q5. The rate of production of turbulent kinetic energy (per unit mass) is given by U i P ( k ) = u i u j x j Show that, in an incompressible flow: (a) only the symmetric part Sij of the mean velocity gradient affects production; (b) only the anisotropic part of the turbulent stress tensor kaij = u i u j 2 k 3 production; (c) with the eddy-viscosity hypothesis, P(k) is greater than or equal to zero.

ij

affects

Linear Eddy-Viscosity Models Q6. Two-equation turbulence models require the solution of two scalar-transport equations, typically for the turbulent kinetic energy k and a second dimensionally-independent variable φ = km n, where is the rate of dissipation of turbulent kinetic energy and m and n are constants. These transport equations take the form: Dk = ( ( k ) k ) + P ( k ) Dt Dφ φ = ( ( φ ) φ) + (C φ1 P ( k ) C φ 2 ) + S ( φ ) Dt k where P(k) is the rate of production of turbulent kinetic energy. (a) If the kinematic eddy viscosity k2 =C t

t

is given by

find an equivalent expression for (b)

t

in terms of k and φ.

S(

)

Show that the transport equation for φ can be transformed into a transport equation for in the form D = ( ( ) ) + (C 1 P ( k ) C 2 ) + S ( ) Dt k where 2 2 S ( φ) m k k (φ) (k ) (φ) )k + = + ( + 2mn + n(n 1) m(m 1) k n φ k k and find expressions for ( ), C 1 and C 2.

[

]

(A personal view: if the transport equation for any second variable can always be transformed into a transport equation for – or vice versa – then it seems that there is no fundamental reason for preferring one particular pair of turbulence variables over any other. Rather, there are pragmatic reasons for the choice, such as the ease of setting boundary conditions or the relationship to measurable physical quantities.)

24

¤

¤

¤

David Apsley

Non-Linear Eddy-Viscosity Models Q7. A quadratic stress-strain relationship can be written in the form a = 2C s + 1 (s 2 1 {s 2 }I) + 2 ( s s ) + 3 ( 2 1 { 2 }I) 3 3 where the anisotropy a, non-dimensional mean-strain s and vorticity are defined by ui u j 2 aij ≡ ij , k 3 U j U j k k 1 U 1 U S ij = ( i + sij = S ij , = ( i ), ), ij ij = ij xi 2 x j 2 x j xi Prove that with this model the non-zero stresses in a simple shear flow are given by uv = C k
u2 2 = +( 3 k 2 v = +( 3 k w 2 = ( k 3
2 2

Q8. Find similar results to those of Q7 for the cases of: V U (a) plane strain ( = , W = 0) y x 1 U V W ) = = (b) axisymmetric strain ( 2 x y z k U In each case, use = . x

1

+6 6

2

2

2 3

) )

12
2

where

=

k U . y

1

2

3

12

1

3

)

6

25

David Apsley

Differerential Stress Models Q9. (CFD 2 Exam, April 2006 – modified) A basic differential stress model (DSM) for turbulence closure in wall-bounded incompressible fluid flow is given (in standard suffix notation) by d ijk D (u i u j ) = + Pij + ij ij Dt x k where U j U i u j uk Pij ≡ u i u k x k x k
ij

=

(1) ij

+

(2) ij

+

( w) ij

( Φ ij1) = C1ε(
( w) ij

ui u j
ij

k = ( kl nk nl ~ ~
ij

2 δ ij ) , 3 3 2 k ~
ik

( Φ ij2) = C 2 ( Pij 1 Pkk δ ij ) 3

n j nk 3 2
( w) 2 ( 2) ij

~
jk

ni n k ) f , f = C 3/ 4k 3/ 2 /

=C

( w) 1

ui u j

+C

,

yn

ε ij = 2 εδ ij 3

dijk is a diffusive flux, ni are the components of a unit wall-normal vector and yn is the distance to the nearest wall. k is the turbulent kinetic energy (per unit mass) and ε is its dissipation rate. C1, C2, C1(w) and C2(w) are model constants. (a) (b) What are the accepted physical roles of the terms denoted Pij, Φij and εij? What is the purpose of the wall correction equilibrium wall boundary layer?
( w) ij

and what is the value of f in an

(c)

Write down expressions for all components of Pij, Φij and εij in simple shear flow (where U1/x2 ≡ U/y is the only non-zero mean-velocity derivative). Assume that n = (0,1,0). Show that, in an equilibrium turbulent flow (where Pij + Φ ij ε ij = 0 ),
( v 2 2 1 + C1 + C 2 2C 2 w ) C 2 = k 3 C1 + 2C1( w) ( u 2 2 2 + C1 2C 2 + C 2 w ) C 2 C1( w) v 2 + = C k k C1 3 1 ( w 2 2 1 + C1 + C 2 + C 2 w ) C 2 C1( w ) v 2 + = C k k C1 3 1 ( 1 C 2 + 3 C 2w) C 2 v 2 uv 2 = k k C1 + 3 C1( w ) 2

(d)

26

David Apsley

Q10. (a)

In the transport equation for u i u j the production term is given in suffix notation by U i x k x k Write this out in full for P11, P22 and P12. Pij ≡ u i u k u j uk U j

(b)

In a simple shear flow, U/y is the only non-zero mean-velocity gradient and the production term simplifies dramatically. Write expressions for all Pij in this case. Give two reasons why the streamwise velocity variance u 2 tends to be larger than the wall-normal component v 2 in flow along a plane wall y = 0.

(c)

Near-Wall Behaviour Q11. By expanding the fluctuating velocities in the form u = a1 + b1 y + c1 y 2 + L v = a 2 + b2 y + c 2 y 2 + L w = a 3 + b3 y + c3 y 2 + L show that u 2 = b12 y 2 + L and derive similar expressions for v 2 , w 2 , uv , k and t.

Q12. Use the turbulent kinetic energy equation and the near-wall behaviour of k from Q11 to show that the near-wall behaviour of is 2 k ( y → 0) ~ 2 ~ constant y (Note that this implies that some modification is required in the equation as y → 0).

27

David Apsley

Wall Functions Q13. (Cell-averaged production and dissipation) An effective viscosity profile is given by + max{0, u 0 ( y y )} eff = where y is the zero-eddy-viscosity height. (a) Show that the mean velocity profile and cell-averaged rate of production of turbulent kinetic energy are given by y+, y+ ≤ y+ U w = 2 × + 1 + + y+ ≥ y+ u0 u 0 y + ln{1 + ( y y )}, ( / )2 ( + y+ ) ( Pavk ) = w ln[1 + ( + y + )] u0 1+ ( + y+ )

respectively, where y + ≡ (b)

yu 0

and

(assumed greater than y ) is the depth of cell.

If the dissipation rate is given by (y ≤ y ) w 3 u0 = (y > y ) (y y ) d where w is such as to make the profile continuous at y = y , show that the cellaveraged dissipation rate is given by 3 yd u0 y )+ ln( av = y yd y yd
¤

28

David Apsley

Q14. (Rough and smooth walls) A generalised form of the log-law mean-velocity profile which satisfies both smooth- and rough-wall limits can be written in wall units as y+ 1 )+B U + = ln( (*) 1 + ck s+ where , B and c are constants and u y u ks U , k s+ = U + = , y+ = u y is the distance from the wall, u is the friction velocity and ks is the Nikuradse roughness. (a) Assuming that (*) holds from the wall to the centreline of a pipe of diameter D, integrate to find an implicit relationship between the skin-friction coefficient cf, pipe Reynolds number Re and relative roughness ks/D. By comparing your results with the Colebrook-White formula: k 1 1.26 = 4.0 log 10 s + 3.7 D Re c cf f deduce values for , B and c. Show that (*) can also be written in the form 1 ~ U + = ln y + + B (k s+ ) (**) ~ and deduce the functional form of B (k s+ ) . (d) In the fully-rough limit ( k s+ >> 1 ) (**) can be written as y 1 U + = ln + Bk ks Use your answers to parts (b) and (c) to deduce a value for Bk.

(b)

(c)

29

David Apsley

10页 免费 计算流体力学(CFD)文档... 18页 免费 计算流体力学(CFD)文档...means trivial see the lectures on turbulence modelling (Sections 7 and ...

produce a second mass-consistent flow field but with time-advanced pressure....10页 免费 计算流体力学(CFD)文档... 20页 免费 计算流体力学(CFD)文档...

2018年计算流体力学(CFD)市场现状与发展趋势预测 (目录).doc
2018-2025 年全球及中国计算流体力学 (CFD)市场调查研究及发展前景趋势分析 报告 报告编号:2271659 中国产业调研网 www.cir.cn 计算流体力学(CFD) 2018-2025 年...

【系统仿真学报】_计算流体力学(cfd)_期刊发文热词逐年....xls
10 科研热词 计算流体力学 自航发射auv 网格收敛性研究 湍流边界层 极限流域 turbulence boundary layer narrow domain grid convergence analysis cfd auv swim-out...

CFD商业软件综合介绍.doc
CFD商业软件综合介绍_能源/化工_工程科技_专业资料 暂无评价|0人阅读|0次下载|举报文档CFD商业软件综合介绍_能源/化工_工程科技_专业资料。计算流体力学(CFD)是...
CFD介绍.pdf
CFD的概念 { 计算流体力学(Computational fluid ...3 (末尾空一行) 湍流模型简介 What is Turbulence?...文档贡献者 新月1921 贡献于2015-05-10 专题推荐...
Fluent计算流体力学软件.pdf