ENCM 7910
Advanced Turbulence Modeling
Reynolds-Averaged Navier-Stokes Flow Solver: Implementation and Potential Coupling with Acoustics Weiyang Lin, Lafayette K. Taylor, and Kidambi Sreenivas University of Tennessee at Chattanooga March 30, 2014
1
Introduction
One of the essences in CFD is its capability of accurately predicting complex flows. Since many flows of engineering interest are turbulent, a suitable treatment of turbulent flows is therefore critical. While the flow field is sufficiently described by the Navier-Stokes equations, the accurate solution to them is usually not practical. Reynolds-averaged Navier-Stokes (RANS) models have been proven successful in the last few decades. The goal of this report is to describe the details of the implementation of two well-known turbulence models in couple with the mean flow in three dimensional grids. One model is the Spalart-Allmaras standard (SA) model, where the turbulent molecular viscosity is modeled in one equation. The other is the Menter Shear Stress Transport (SST) model, in which the kinetic energy (per unit mass) of the turbulent fluctuations and the dissipation per unit turbulence kinetic energy are modeled in two equations. Boussinesq eddy-viscosity approximation is employed to couple the mean flow with the solutions to the turbulence models. This report consists of two major parts. The first part covers through chapters 2-6, including the descriptions of the governing equations, their nondimensionalization, turbulence models, numerical discretization, results on a flat plate, and related discussions. The second part, chapter 7, discusses about some aspects of acoustic research that might be related to the current topic. One is the possibility of coupling of turbulent flow solvers with acoustics. Finite element method and other high order schemes are considered to solve the linearized Euler equations (LEE) in the context of hybrid methods with acoustic analogies.
2
Governing Equations
2.1
Navier-Stokes Equations Since air is mostly common compressible flow in engineering problems, the flow solver is set to solve the air in normal situations. Several constants and characteristics of air at standard atmosphere are listed in Table 1 below. Table 1. Constant and characteristics of air (standard atmosphere). The main dimension used are mass M , length L , time T , and absolute temperature . Quantity Dimension SI 2 2 1 R LT 287.053J .kg 1.K 1 1.2250kg.m3 ML3
p T
ML1T 2
ML1T 1
1.0132 105 Pa 288.15K 1.78935 105 kg.m1.s 1
The governing equations are the time-dependent compressible Navier-Stokes equations. The equations of conservation laws can be written in the tensor form as [1] Conservation of mass: (1) u 0 , t
Page 1 of 24
ENCM 7910
Conservation of momentum: Conservation of total kinetic energy: In Eq. (2),
Advanced Turbulence Modeling
u uu f p t et b f u t
(2) (3)
u I 2 D ,
and 1 T u u , 2 2 , 3 the bulk coefficient of viscosity is usually ignorable. In Eq. (3), D
b et p u u k T , and in both equations the body force f can be ignored in the current context. By expanding Eqs. (1)-(3), we obtain the vector form of the equations as Q fi gi hi f v gv hv 0 t x y z x y z where Q , f i , g i , hi , f v , g v and hv are vectors given by
(4)
u v w u 2 p uv uw u Q v , fi uv , gi v 2 p , hi vw , 2 uw vw w p w et p u et p v et p w et 0 0 0 xx xy xz , gv , hv , xy yy yz fv xz yz zz u xx v xy w xz qx u xy v yy w yz q y u xz v yz w zz qz where the components of the viscous stress tensor ij are given by 2 u v w 2 u v w 2 u v w , yy 2 , zz 2 , 3 x y z 3 x y z 3 x y z u v v w u w xy yx , xz zx , yz zy , z x y x z y and Fourier's law for heat transfer by conduction will be assumed so that the heat transfer q can be expressed as q k T , where k is the coefficient of thermal conductivity and T is the temperature. The equations are closed by introducing the Sutherland's law:
xx 2
2
T 3 T0 S , 0 T0 T S Page 2 of 24
ENCM 7910
Advanced Turbulence Modeling
where the constants are T0 273K , 0 1.716 105 N
s , S 111K . m2
2.2
Non-Dimensionalization The equations are nondimensionalized by free stream density , speed of sound c , temperature T , viscosity , thermal conductivity k , and a reference length L . The nondimensionalization yields the two groups of nondimensionalized free-stream variables, the first has * 1 , c* 1 , T* 1 , * 1 , k* 1 , and the second has 1 c2 V V* p 1 V * L * , , , M V* . p* 2 V t c c* c c2 c V The nondimensionalization leads to several flow constants as p* 1 1 1 R* * * , C *p , Cv* . 1 1 T The substitution of the nondimensionalized variables to the governing equations leads to Q* fi * gi* hi* f v* gv* hv* 0 (5) t x y z x y z where the vectors in the equation have exactly the same form as those in Eq. (4), the only difference is the construction of the viscous fluxes, which can be expressed as 2 M u * v* w* xx* * 2 * * * , 3 Re x y z
2 M u * v* w* *yy * * 2 * * , 3 Re x y z 2 M u * v* w* zz* * * * 2 * , 3 Re x y z
xy* *
* * M u * v* * w* * * w* * * * M u * M v * * yx , xz * * zx , yz * * zy , Re y x Re z x Re z y * * * * * * M M M T T T q*x , q*y , q*z , * * Re 1 Pr x Re 1 Pr y Re 1 Pr z *
where the Reynolds number
Re
c L ,
the Mach number
M
V*
RT
V* .
The Sutherland's formula is also nondimensionalized as
* T Tref 0 T0 or it will take the following form if Tref
2
3 T0 S , * T T S ref T0 and ref 0 : *
Page 3 of 24
ENCM 7910
Advanced Turbulence Modeling
* T * 3 2
1 C* , T * C*
where S 111K . T0 273K For future convenience, the asterisks will be dropped from the nondimensional equations. C*
2.3
Equations of State Several equations of state are listed here (with duplicates) for reference: C h C pT , p 1.4 , R C p Cv , e Cv T , p RT , Cv
c RT
p
,
M
V , c
h
c2 , 1
p 1 e , Pr
et e
1 2 u , 2
ht h
1 2 u , 2
Cp
0.72 . k The nondimensionalization leads to an additional equation for temperature: M2p p 2 T c .
3
Turbulence Models
The concept of Reynolds averaging varies from time averaging, space averaging to ensemble averaging [2]. Time averaging is appropriate when considering a stationary turbulence. The main idea of Reynolds averaging is to decompose the flow to averaged and fluctuating components: ui ui ui , p p p and t , where the variables with bars are averaged values and the ones with primes represent turbulent fluctuations. Using the Boussinesq approximation and the assumption of a constant turbulent Prandtl number, the Reynolds-averaged Navier-Stokes equations can be written as u j 0 , t x j
ui p ij ui u j , t x j xi x j et u j ht ui ij t x j x j x j
t Pr Prt
T , x j
uj 1 u 1 1 1 where et e ui ui ui' ui' , ij 2 t Sij Skk ij and Sij i . The 2 x j xi 2 2 3 nondimensionalization of the equations above looks very similar to Eq. (5), except for the 1 additional turbulent kinetic energy ui' ui' (which is often ignored in the definition of total energy), 2 eddy viscosity t and the assumption of a constant turbulent Prandtl number Prt 0.90 . RANS models based on the Boussinesq eddy-viscosity approximation have succeeded to predict many flow properties, however, they can be inaccurate in some circumstances. Such situations include [3] 1. Flows with sudden changes in mean strain rate Page 4 of 24
ENCM 7910
Advanced Turbulence Modeling
2. Flow over curved surfaces 3. Flow in ducts with secondary motions 4. Flow in rotating fluids 5. Three-dimensional flows. In addition, corner region is where turbulence models would potentially fail [4]. Second order models, such as algebraic stress model (ASM) and the Reynolds stress model (RSM), could be implemented to overcome the plausible Boussinesq hypothesis and the limitations of the first order models in dealing with the isotropy of turbulence and the extra strains. A general form for a PDE that governs transport turbulence models is given by Z Z 1 Z Ui T P Z D Z , t xi xi xi where Z represents the transported turbulence quantity. Typically closure is obtained by a relationship between the eddy viscosity T and the turbulence quantity Z , in addition of course to the mean flow variables. 3.1
Spalart-Allmaras Model The one-equation model is given by the following equations (NASA Langley) 2 cb1 1 uj cb1 1 ft 2 S cw1 f w 2 f t 2 cb 2 x j x j xi xi , t d x j production diffusion convection destruction t f v1 where
3 , , 3 3 cv1 and is the molecular kinematic viscosity, and is the molecular dynamic viscosity. In addition, S 2 2 f v 2 , where 2WijWij is the magnitude of vorticity, d is the distance d f v1
from the field point to the nearest wall, and 1
1 c6 6 , f w g 6 w36 , g r cw2 r 6 r , fv 2 1 1 f v1 g cw3 1 u u j r min 2 2 ,10 , ft 2 ct 3 exp ct 4 2 , Wij i 2 x j xi S d
. The term f t 2 is sometimes neglected since it is originally derived from the trip term. The wall boundary condition is set to be wall 0 . Spalart and Allmaras [5] suggested the farfield boundary condition as 10 , but later Spalart [6] suggested 3 ,5 , if the trip term is being left out. The constants are cb1 0.1355 , 2 3 , cb 2 0.622 , 0.41 , c 1 cb 2 cw 2 0.3 , cw3 2.0 , cv1 7.1 , ct 3 1.2 , ct 4 0.5 , cw1 b21 .
The nondimensionalized equation can be expressed by M * * u *j * cb1 1 ft 2 S * * t x j Re Page 5 of 24
cb1 * c f w1 w 2 ft 2 d *
2
ENCM 7910
Advanced Turbulence Modeling
* * * * * * * cb 2 * * , x j xi xi x j the equation above can also be rearranged to 2 * * cb1 * M * u j * * cb1 1 f t 2 S cw1 f w 2 f t 2 * t x*j Re d
1 M Re
1 M Re
* * 2 * * 1 cb 2 * * cb 2 * *2 x j xi x j
(6)
(7)
note that
V * * 2 2 *2 f v 2 , L L d but after it is plugged into the equation, this term will end up with the following form M * S * * 2 *2 f v 2 , Re d also, M r min ,10 . 2 2 Re S d The Spalart-Allmaras turbulence model is a very stable and generally reasonably accurate model for a wide range of turbulent flows. One drawback is the model requires the calculation of the distance to the nearest wall for all field points. This can be an expensive computation, especially for unstructured grids. A few hints of implementing this model are given by [4] 1. The eddy viscosity should be limited so that it will not run away in some complex applications. Generally a limit of t 200,000 is acceptable. 2. This model tends to smear out three-dimensional vortical flows. Rotation and curvature corrections (SARC and ASARC) can significantly improve the results. 3. This model can over-damp some unsteady flows. 4. This model contains no corrections for compressibility and will over-predict the growth rate of high speed shear layers. S*
3.2
Menter SST Standard Model The standard Menter SST (Shear Stress Transport) equations below [7] (written in conservation form) k u j k P * k t x j x j
u j
two-equation model is given by the
k k t , x j
2 k , t 2 1 F1 t x j x j x j x j after nondimensionalization (drop the asterisks), k u j k M M k P * k k t , t x j Re Re x j x j t
u j x j
P 2 t x j
(8)
2 k M P 2 , (9) t 2 1 F1 t Re x j x j x j x j
Page 6 of 24
ENCM 7910
where P ij
Advanced Turbulence Modeling
ui , and x j
ij t 2Sij
2 uk Re 2 ij k ij . 3 xk M 3
The strain Sij is defined like before. Note that it is generally recommended to use a production limiter
Re P min P, 20 * k . M The turbulent eddy viscosity is computed from a1k Re t , M max a1 , F2
(10)
The magnitude of vorticity is defined same way like before, and it is replaced by the strain magnitude in some applications. Each of the constants ( , k , , ) is calculated via
F11 1 F1 2 ,
where 1 and 2 represents constant 1 and 2, respectively. Additional functions are given by
k M 500 4 2 k , F1 tanh arg14 , arg1 min max * , 2 , 2 d Re d CDk d
CDk
k CDk max 2 2 ,1020 , F2 tanh arg 22 , x j x j k M 500 arg 2 max 2 * , 2 , t t . d Re d can also be plugged into Eq. (9) during the calculation. The other constants are given by
1
1 1 2 2 2 2 , , k1 0.85 , k 2 1.0 , 1 0.5 , 2 0.856 , 2 * * * * 1 0.075 , 2 0.0828 , * 0.09 , 0.41 , a1 0.31 .
The boundary conditions recommended in the original reference are 105 M 2 101 M 2 k farfiled M farfiled 10M , , Re L Re L M 6 wall 10 , kwall 0 . Re 1 d1 2
4
Numerical Discretization
4.1
Finite Volume Method for the Navier-Stokes Equations The finite volume discretization takes the integral form of Eq. (5) and transform it with the divergence theorem. The final form of the discretized equations can be expressed as Qi Vi Fij nˆij Sij 0 , ti j where F here represents both inviscid and viscous fluxes. A typical 2D control volume constructed by surrounding medium duals is shown in Fig. 1 below. The formulation of 3D control volumes follows the same method, though it involves more efforts on the construction of Page 7 of 24
ENCM 7910
Advanced Turbulence Modeling
medium dual surfaces. It should be noted that the potential division by a small number should be avoided when one calculates the normal of a control surface, especially when dealing with grids with tight spacing. Roe fluxes are used here to calculate the inviscid flux; the flux through an edge i, j is defined as 1 ij Fij nˆij Fi Qi ; nˆ Fi Q j ; nˆ A Qi , Q j ; nˆ Q j Qi , 2 where i and j are the nodes that make up the edge. A is the Jacobian matrix evaluated using Roe averaged values, which can be practically either treated as constant or calculated using automatic differentiation wherever flux Jacobians are needed. The constant treatment may be written down as ij Fij nˆij 1 Fi Qi ; nˆ A Qi , Q j ; nˆ , Qi Qi 2 Qi ij Fij nˆij 1 Fj Q j ; nˆ A Qi , Q j ; nˆ . Q j Q j 2 Q j
Figure 1.
A typical 2D node-centered control volume with the control surfaces illustrated by slashed lines
A second-order spatial accuracy can be achieved by the reconstruction of conservative variables on the control surfaces. Consider the truncated Taylor series expansions in space from the points xi and x j to the midpoint of the edge xe , a more accurate solution could be achieved by using available gradient information at a node to approximate the left and right states: U L U i U i xe xi , U R U j U j xe x j .
The gradients can be approximately solved by the least squares method, which can be summarized by 1. Precompute N
N
N
e 1 N
e 1 N
s11 xe , s12 xe ye , s13 xe ze , e 1 N
2
s22 ye , s23 ye ze , s33 ze , e 1
2
e 1
2
e 1
where xe j xe . For unweighted least squares, the constant j is unity; for inverse distance
Page 8 of 24
ENCM 7910
Advanced Turbulence Modeling
weighted least squares, j j s j (where s j s ). The unweighted least squares method is often used for the inviscid fluxes. 2. Let r11 s11 , r12 s12 , r13 s13 , r22o r22 , r33o r33 ,
r22 s22
r2 r2 r122 r , r23 s23 12 r13 , r33 s33 13 23 , r11 r11 r22 r11
and
r 1 1 xe r12Wey r13Wez , Wey o ye 12 xe r23Wez , r22 r12 r11 r r r 1 Wez o ze 13 xe 23 ye 12 xe . r33 r11 r22 r11 3. Construct the gradients Wex
N
N
N
u xi Wex u j e ui , u yi Wey u j e ui , u zi Wez u j e ui , e 1
e 1
e 1
and reconstruct the left and right states. Note that points 3 and 4 can be combined into one step during implementation. The viscous fluxes are computed with a finite volume technique with a direct approximation for the gradients at the quadrature points, referred to as "directional derivative" technique. The formula for the edge gradients is sij , U ij U ij U j U i U ij sij s ij
which follows the evaluation of nodal gradients with the weighted least squares method. The average of nodal gradients U ij , which significantly increases the stencil in the calculation of Jacobian matrix, will be treated as constant. 4.2
Numerical Steps of Coupling with Turbulence Models In this implementation, both models are solved in a manner that is loosely-coupled with the mean flow; that is, at each time step, the solver would attempt to obtain the mean flow solution first, and then to solve the turbulence models separately. 4.2.1 Spalart-Allmaras Standard Model Equation (6) is a convective-diffusive equation which can be discretized by two steps of calculation: one is for the convective part, and the other is for the source terms. The equation can be recast into the following form (drop the asterisks) 1 M U t Re A
B
M cb1 1 f t 2 S Re
cb1 1 M 2 cw1 f w 2 f t 2 d Re cb 2 , 2
C
this leads to the discretization of the integral form of the equation with corresponding terms i 1 M Vi H ij nˆij Sij G nˆ S V ti Re j ij ij ij i C i j A
B
where convective term A is formulated using the single upwind method
Page 9 of 24
ENCM 7910
Advanced Turbulence Modeling
H ij nˆij i j , where
,
and U ij nˆij . 2 2 This gives the Jacobian contribution from term A as H ij nˆij , H ij nˆij . j i Diffusive term B can be calculated in the same way as the viscous fluxes; that is, can be computed using the directional derivative sij , ij ij j i ij sij 2 sij and the multiplier ij 1 cb 2 ij can be evaluated using average. This gives the Jacobian contribution from term B as
Gij nˆij i i j
G
ij
nˆij
i j i j 2 2
sij nˆij ij i , ij nˆij 2 s ij
sij nˆij i j i j ij j . ij nˆij 2 j 2 2 s ij
The source term C is integrated using the local value at the node i . Its calculation is relatively straightforward. To calculate the Jacobian, note that M Si 2 i 2 f v 2 , Re d substitute to the term, 2 c M M 1 M 2 i cb1 1 ft 2 2 i 2 f v 2 i cw1 f w b21 ft 2 i cb 2 , Re d Re d Re the Jacobian will only contribute to the main diagonal if is treated as constant in this case. This gives i c M f M cb1 1 ft 2 2 v22 2i 2 cw1 f w b21 f t 2 i2 . i Re d Re d For stability concern, only non-negative contributions are placed on the main diagonal of the global flux Jacobian. An alternative discretization can be done based on Eq. (7). The difficulty comes from the 2 discretization of the term cb 2 2 . Since it is not in divergence (or conservative) form, it cannot xi be integrated by the traditional finite volume way. A plausible trick can be applied here; that is, to approximate this term as , cb 2 * xi xi where * is taken as the bias nodal value rather than the edge-associated average. 4.2.2 Menter SST Standard Model The Mentor SST model of Eqs. (8) and (9) is discretized via the very similar way as the SpalartAllmaras model. The convective terms are calculated by simple upwind scheme, the diffusive terms Page 10 of 24
ENCM 7910
Advanced Turbulence Modeling
are approximated with directional derivatives, and the productions and destructions are treated as source terms. Since the conservative variables of k and are not in direct usage during implementation, the primitive variables k and are stored instead. Thus the update of the solution yields
qT
QT
,
where QT k , , and qT and are the primitive variables and local density, respectively. Besides, the flux Jacobians are again forced to be diagonal dominant by adding non-negative contributions for stability concerns. T
4.3
Time Step Integration In this report, only the 1st-order time step (finite difference) is implemented. For transient problems, an implicit Runge-Kutta scheme may be considered.
5
Results
5.1
Laminar Flow Over a Flat Plate Laminar flow over a flat plate is typically the first benchmark problem to solve in order to validate the code. With the dimensionless similarity variable set to be U Re , y y* 2 x 2x the normalized velocity profile u U can be compared with the analytic Blasius flat-plate solution [8]. The case is ran on a flat plate with length of 1.0. It has infinitely long width and infinitely small thickness (due to the symmetric boundary condition). The grid being used has 50 points on the x-direction (parallel to the flow) of the plate and 2 points on the z-direction; the upstream and downstream directions are extended towards the x-coordinates of -6.0 and 6.0 with 11 grid points, and all of the points are extruded towards the y-coordinates of 5.0 for 31 points in total with the initial spacing of 1E-3 to resolve the laminar boundary layers. A side look of the grid is shown in Fig. 2. In this case, the free-stream flow has the Reynolds number of 1E4 and a low Mach number of 0.2, this results in the general contour of u-velocity as is plotted in Fig. 3. The normalized velocity profile is compared with the Blasius solution below in Fig. 4.
Figure 2.
A side look of the grid being used to simulate the laminar flow over a flat plate. Page 11 of 24
ENCM 7910
Advanced Turbulence Modeling
Spacing off the solid wall is 1E-3
Figure 3.
The u-velocity contour from the solution to the laminar flow over a flat plate. Flow parameters are: Re = 1E4, Mach = 0.2
Figure 4.
Normalized u-velocity against plot in comparison with the Blasius solution
Besides, the skin friction, which is a more sensitive parameter to verify the solution, can also be plotted. The wall shear stress is given by a u u * w * * , y w L y w and the skin friction can be calculated by * 2 w 2 * u Cf U 2 a L * U *2 y *
w
* M 2 * u . Re * U *2 y * w
It is plotted against Re x in comparison with the Blasius solution 0.664
Page 12 of 24
Re x in Fig. 5 below.
ENCM 7910
Figure 5.
Advanced Turbulence Modeling
Skin friction C f against Re x plot in comparison with the Blasius solution (to be modified)
5.2
Turbulent Flow Over a Flat Plate The Spalding velocity profile is derived from a dimensionless analysis based on the theory of Ludwig Prandtl and Theodore von Karman. The theory consists of three laws. The inner law states that the profile would depend upon wall shear stress, fluid properties, and distance y from the wall. The outer law states that the wall acts merely as a source of retardation, reducing local velocity U below the stream velocity U e in a manner independent of viscosity but dependent upon wall shear stress, layer thickness, and freestream pressure gradient. The overlap law simply specifies the inner and outer functions to merge together smoothly over some finite intermediate region. Define friction velocity by Cf v* U , 2 and it gives the normalized u and y by (nondimensionalized)
Re yv* U y and . M v* Then the inner and outer laws may be written in U 1 yv* B * ln v , U e U 1 y v* k ln A or in a compact form (Spalding 1961) 2 3 u u B u , y u e e 1 u 2 6 where the constants can be either 0.40 and B 5.5 , or 0.41 and B 5.0 . u
Page 13 of 24
ENCM 7910
Advanced Turbulence Modeling
The cases are ran on the same flat plate as the laminar case. The initial spacing of 1E-5, 5E-6, 1E-6 and 5E-7 are chosen to resolve the turbulent boundary layers. The free-stream flow has the Reynolds number of 5E6 and a low Mach number of 0.2. For the SA model, the farfield boundary condition is chosen to be 1.341946 , and the trip function f t 2 is neglected; for the SST
2 103 M 2 model, the farfield boundary condition is chosen to be k , 4 M , . This results Re L in the general contour of u-velocity as is plotted in Fig. 6, and a typical contour of eddy viscosity is plotted in Fig. 7. Both pictures are from the SST model.
Figure 6.
The u-velocity contour from the solution to the turbulent flow over a flat plate, solution from SST model. Flow parameters are: Re = 5E6, Mach = 0.2
Figure 7. The turbulent eddy viscosity t contour from the solution to the turbulent flow over a flat plate, solution from SST model. Flow parameters are: Re = 5E6, Mach = 0.2 The normalized velocity profile is compared with the Spalding solution below in Figs. 8-10. Both models are ran in 4 grids with different grid spacing off the wall, as stated before. Although not sure if it should be the way, it is found that the SA model has relatively more consistent than the SST model subject to the change of grids (Figs. 8 and 9). The SST model tends to "converge" to a solution as the wall spacing becomes tighter. The velocity profiles of the SA solution at different locations of the plate, range from 1st quarter, 2nd quarter to 3rd quarter, are plotted in Fig. 10. The varying plots are due to the change of the skin friction as a denominator in the equation. Skin friction is a more sensitive parameter to validate the solution. The White-Christoph formula 0.455 Cf 2 ln 0.06 Re x is used for comparison to the skin friction coefficient results. Results on 4 different grids are plotted in Figs. 11 (SA model) and 12 (SST model). Again, the SA solution is more stable in terms of grids, and the SST solution "converges" towards the White-Christoph solution. Figure 13 shows the comparison of C f from SA and SST solutions on a same grid with off-wall spacing of 5E-7. The SST solution is closer to the analytic solution than SA.
Page 14 of 24
ENCM 7910
Figure 8.
Advanced Turbulence Modeling
Velocity profile u against log10 y plot (solutions from SA model) in comparison with the Spalding solution
Figure 9.
Velocity profile u against log10 y plot (solutions from SST model) in comparison with the Spalding solution
Page 15 of 24
ENCM 7910
Figure 10.
Advanced Turbulence Modeling
Velocity profile u against log10 y plots (solutions from SA model) at three locations on the flat plate, in comparison with the Spalding solution
Figure 11.
Skin friction coefficient plots on different grids (solutions from SA model) in comparison with the White-Christoph solution
Page 16 of 24
ENCM 7910
Figure 12.
Figure 13.
Advanced Turbulence Modeling
Skin friction coefficient plot on different grids (solutions from SST model) in comparison with the White-Christoph solution
Skin friction coefficient plot on the same grid (with 5e-7, solutions from SA and SST models) in comparison with the White-Christoph solution
5.3
Turbulent Flow Over an Airfoil The RANS solver is ran on NACA0012 airfoil out of interest. Results from both SA and SST models are presented in Figs. 14 and 15, and the SST results are compared with the corresponding ones obtained from Tenasi (Fig. 16). It is shown that the two models give similar mean flow results (u-velocity), though the eddy viscosities t differ from each other. Page 17 of 24
ENCM 7910
Advanced Turbulence Modeling
Figure 14. Results of turbulent flow over NACA0012, solution obtained from SA model. Left is the u-velocity contour, and the right is the eddy viscosity one. Flow parameters are: Re = 5E6, Mach = 0.2
Figure 15. Results of turbulent flow over NACA0012, solution obtained from SST model. Left is the u-velocity contour, and the right is the eddy viscosity one. Flow parameters are: Re = 5E6, Mach = 0.2
Figure 16. Results of turbulent flow over NACA0012, solution obtained from SST model by Tenasi. Left is the u-velocity contour, and the right is the eddy viscosity one. Flow parameters are: Re = 5E6, Mach = 0.2
6
Discussions
In the first part of the report, a basic detailed implementation with two RANS turbulence models, Spalart-Allmaras one equation model and Menter SST two equation model, is shown with corresponding results on a flat plate. Future work includes the investigation of the behavior turbulence models in predicting turbulent flows in more complicated situations, such as for airfoils. Turbulence models are often being adjusted to particular flow cases. People come up with various ideas of modeling certain behaviors of turbulent flows and then try to fit the results by tweaking certain constants and/or functions, usually of those appear in the diffusive and destruction terms, since the dissipation of the turbulent behavior can sometimes be arbitrary. For instance, the magnitude of vorticity is exchangeable with the magnitude of strain in some engineering applications. It is also reported that small changes (5-10%) in modeling constants can lead to significant modifications on the simulation predictions [7]. Besides, even if the turbulence model is not asymptotically consistent, the mean flow profile and the wall skin friction can still be predicted correctly. There was probably no major conceptual break done since the introduction of RANS. An alternative way of solving turbulent flows is direct numerical simulation (DNS). DNS is of great importance in that the results obtained could be considered equivalent to experiments. On Page 18 of 24
ENCM 7910
Advanced Turbulence Modeling
the other hand, DNS is usually not affordable. For example, the number of uniformly distributed t 49 grid points would be Nuniform 110ReT , and number of time steps is given by Ntime total , t 0.003 L where t [3]. * ReT v Large-Eddy simulation (LES) is currently an active research area. LES is a computation where large vortexes (eddies) are computed directly, while small scale eddies are modeled, hence the LES is less expensive. Another important concept regarding LES is filtering. Sub-grid fluctuations should be modeled by averaging in order to be removed. LES is, however, approximately 20 times slower than RANS models. A combination of LES and RANS gives yields two types of hybrid methods. The first one models the near-wall structures using a RANS-type model, and the outer effects with LES. The second one is the so-called detached eddy simulation (DES), which uses RANS for the entire attached boundary layer, while LES is used to treat the separated region of the flow field. It is commented that DES can be easily implemented in an existing RANS code with a small modification of the turbulence model, but the strong numerical dissipation inherent in typical RANS algorithms is detrimental to the LES region, and the artificial transition from RANS to LES requires special attention [9].
Page 19 of 24
ENCM 7910
7
Advanced Turbulence Modeling
Aeroacoustics
Aeroacoustics is a branch of acoustics that studies noise generation via either turbulent fluid motion or aerodynamic forces interacting with surfaces (Wikipedia). Prompted by the need for quieter jet engines, noise reduction is also in demand in marine propellers, hydrofoils, and transitional and turbulent boundary layers on sonar domes [9]. The hierarchy of noise prediction methods is summarized in Fig. 17 below .
Figure 17.
A hierarchy of noise prediction methods [9]
While this first literature review provides some preliminary investigation of potential numerical solution with linearized Euler equations (LEE) to the aeroacoustic problems, several cutting-edge aspects of research shall also be noted for future reference. Regarding the challenges in CAA, Kaltenbacher, et al. [10] summarized some main difficulties which have to be considered for the simulation of flow noise problems including: Energy disparity and acoustic inefficiency. Length scale disparity. Preservation of multipole character. Dispersion and dissipation. Flows with high Mach and Reynolds number. Simulation of unbounded domains. This remains to be a very active filed of research. 7.1
Introduction to the Hybrid Approach The basis of numerical solutions with hybrid methods can be traced back to Lighthill's well-
Page 20 of 24
ENCM 7910
Advanced Turbulence Modeling
known analogy for the L-S decomposition1 [11] 2Tij 2 2 2 , 2 c0 xi xi xi x j t L
S
where the Lighthill's stress tensor is given by Tij ui u j ij p c02 ij .
One of the most important part of this analogy is the transformation of the acoustic sources from the computed flow data to the acoustic field; that is, the interpretation of the source term where mean flow effects on the wave propagation are included. The most general form of Lighthill's analogy is the extension developed by Ffowcs Williams & Hawkings (FW-H formulation), which incorporates the effect of surfaces in arbitrary motion [12], and some of the difficulties involved in Lighthill's analogy motivated a generalized acoustic theory by Goldstein [13]. A illustrative discussion about acoustic analogy can be found in [9]. Figure 18 is helpful to illustrate the strategy of coupling aeroacoustic analogies with the nonlinear solution obtained by CFD solvers. On one hand, integral methods are widely used for open domain problems. Once the unsteady flow solution is obtained, one can evaluate acoustic sources using the fluid field, by assuming that there is no significant physical influence occurs from the acoustic propagation to the flow field. On the other hand, for confined aeroacoustic problems, or if structural/acoustic effects are present, a volume discretization method to account for the interactions between the solid surfaces and the flow-induced noise shall be employed. This includes finite differences (FD), discontinuous Galerkin (DG) or finite volume methods (FVM), for solving the systems of equations like LEE or acoustic perturbation equations, and so on. Traditional finite element methods (FEM) are also used to solve the variational formulation of Lighthill's acoustic analogy [10]. Regard the FVM, a 2nd-order time conservative FVM can solve relevant aeroacoustic problems with high–fidelity [14]. Two issues remain a concern regarding the numerical methods. The first is the accurate boundary conditions, and the other is the spatial resolution of numerical schemes, often quantified as the number of points per wavelength needed to represent advection or wave propagation to within some error threshold. With the exception of Fourier spectral methods, discretizations include artificial dispersion and, in some cases, dissipation due to spatial discretization. For example, it is reported in [15] that FW-H calculations produce results that are equivalent to the more direct but costlier computational approaches; besides, CFD failed to capture higher frequencies. DNS with computation of sound generation was first studied by [16]. Recently a promising heterogeneous domain decomposition technique has been proposed [17], as a way to overcome the scales disparities and still allow a direct simulation. Later LES and DES become a tool to facilitate the simulation. Wang and Moin [18] applied incompressible LES in conjunction with Lighthill's theory to simulate the trailing-edge aeroacoustic experiments. Hedges [19] computed noise radiated from an aircraft landing gear model using DES and FW-H equation and compare the results with those from URANS computations.
1
An acoustic analogy is formulated by rearranging the Navier-Stokes solution R Q 0 into LQ S Q ,
where L is some (usually linear) wave propagation operator and S Q is its corresponding nominal (nonlinear) sound source.
Page 21 of 24
ENCM 7910
Figure 18.
Advanced Turbulence Modeling
Schematic depicting some of the possible strategies when using an aeroacoustic hybrid approach [10]
7.2
Linearized Euler Equations LEE has the advantage that it makes possible the computation of noise propagation effects in non-uniform flows outside the turbulent region. In addition, as stated in section 7.1, the effects due to the presence of complex geometries in the flow are implicitly taken into account in the acoustic solution in LEE, since refraction effects are considered in it. The disadvantage is the higher computational costs. LEE can be derived from the Navier-Stokes equations by linearized the mean flow, turbulent and acoustic variables, by assuming that the acoustic components are much smaller than mean and turbulent counterparts. The equation may be written as [20] uai u U i 1 pa a P U j ai uaj t x j x j xi 2 xi U j
xi U i u 1 p ui u j u j i u j ui , x j x j x j xi t x j
where a refers to acoustic components and the prime superscript refers to turbulent ones. Besides, the turbulent velocity field needed to compute the LEE source terms is obtained using the method of stochastic noise generation and radiation (SNGR). Successful applications of LEE have been found in literatures, which implement specific source-terms [21] or projected source-terms [22]. The numerical difficulty of this approach stem from the fact that the full LEE set (say about a transversely sheared flow) admits non-trivial instability wave solution of the homogeneous equations, which would contaminate the near acoustic field and the errors inherent in the boundary schemes or mesh non-uniformities spuriously reflect them into radiating acoustic waves. Several approaches have been proposed to overcome this problem [22, 23]. Page 22 of 24
ENCM 7910
Advanced Turbulence Modeling
Shen and Tam [24, 25] computed jet screech by URANS with Euler equations. URANS models usually are, however, unable to provide broadband source information because of some modeling ambiguities such as the concept of eddy viscosity (2nd-order turbulence models may relieve the problem), whereas RANS calculations are insufficient by themselves but they may be used to provide some statistical descriptions of acoustic source terms. Applications with broadband range of frequencies of sound energy is distributed in turbulent flows may be solved with the assistance of broadband noise source models including source terms in LEE, others include Proudman's formula, jet noise source model, boundary layer noise source model, source terms in Lilley's equation. However, these source models can only be used to determine which portion of the flow is primarily responsible for the noise generation (and even the low frequency components), but not to predict the sound at the receivers [20]. To validate the numerical solutions, one can either compare them with a DNS solution, or a closed form solution to the Lighthill equation for an unsteady flow in an unbounded domain. 7.3
Applications Wei [26] employed DNS and an adjoint method to explore means of active flow control for noise control in a free shear flow. In comparison, passive control by modifying the shapes of solid boundaries in the flow has been employed in many practical applications; for example, the use of serrated nozzles is effective in reducing jet noise [27]. See also [28-31] for shape optimizations for noise reduction.
References 1. 2. 3. 4. 5. 6. 7. 8. 9. 10.
11.
12.
13. 14.
Warsi, Z. U. Fluid dynamics: theoretical and computational approaches: CRC press, 2005. Sodja, J. "Turbulence models in CFD," University of Ljubljana, 2007. Wilcox, D. C. Turbulence modeling for CFD: DCW industries La Canada, CA, 1998. Briley, W. R. "ENCM 7340 Viscous Flow Computation," 2012. Spalart, P. R., and Allmaras, S. R. "A one equation turbulence model for aerodinamic flows," AIAA journal Vol. 94, 1992. Spalart, P. "Strategies for turbulence modelling and simulations," International Journal of Heat and Fluid Flow Vol. 21, No. 3, 2000, pp. 252-263. Menter, F. R. "Two-equation eddy-viscosity turbulence models for engineering applications," AIAA journal Vol. 32, No. 8, 1994, pp. 1598-1605. White, F. M., and Corfield, I. Viscous fluid flow: McGraw-Hill New York, 1991. Wang, M., Freund, J. B., and Lele, S. K. "Computational prediction of flow-generated sound," Annu. Rev. Fluid Mech. Vol. 38, 2006, pp. 483-512. Kaltenbacher, M., Escobar, M., Becker, S., and Ali, I. "Computational aeroacoustics based on Lighthill’s acoustic analogy," Computational Acoustics of Noise Propagation in Fluids-Finite and Boundary Element Methods. Springer, 2008, pp. 115-142. Lighthill, M. J. "On sound generated aerodynamically. I. General theory," Proceedings of the Royal Society of London. Series A. Mathematical and Physical Sciences Vol. 211, No. 1107, 1952, pp. 564-587. Williams, J. F., and Hawkings, D. L. "Sound generation by turbulence and surfaces in arbitrary motion," Philosophical Transactions of the Royal Society of London. Series A, Mathematical and Physical Sciences Vol. 264, No. 1151, 1969, pp. 321-342. Goldstein, M. "A generalized acoustic analogy," Journal of Fluid Mechanics Vol. 488, 2003, pp. 315-333. Aybay, O. "Time-conservative finite-volume method with large-eddy simulation for computational aeroacoustics." Durham University, 2010. Page 23 of 24
ENCM 7910
15.
16.
17.
18. 19. 20. 21. 22. 23.
24. 25. 26. 27. 28. 29.
30.
31.
Advanced Turbulence Modeling
Khorrami, M. R., and Lockard, D. P. "Effects of geometric details on slat noise generation and propagation," International Journal of aeroacoustics Vol. 9, No. 4, 2010, pp. 655-678. Mitchell, B. E., Lele, S. K., and Moin, P. "Direct computation of the sound from a compressible co-rotating vortex pair," Journal of Fluid Mechanics Vol. 285, No. 1, 1995, pp. 181-202. Utzmann, J., Schwartzkopff, T., Dumbser, M., and Munz, C.-D. "Heterogeneous domain decomposition for computational aeroacoustics," AIAA journal Vol. 44, No. 10, 2006, pp. 2231-2250. Wang, M., and Moin, P. "Computation of trailing-edge flow and noise using large-eddy simulation," AIAA journal Vol. 38, No. 12, 2000, pp. 2201-2209. Hedges, L., Travin, A., and Spalart, P. "Detached-eddy simulations over a simplified landing gear," Journal of Fluids Engineering Vol. 124, No. 2, 2002, pp. 413-423. Fluent, I. "FLUENT 6.3 user’s guide," Fluent documentation, 2006. Bogey, C., Bailly, C., and Juvé, D. "Computation of flow noise using source terms in linearized Euler's equations," AIAA journal Vol. 40, No. 2, 2002, pp. 235-243. Ewert, R., Meinke, M., and Schroeder, W. "Comparison of source term formulations for a hybrid CFD/CAA method," AIAA paper Vol. 2200, No. 7, 2001. Ewert, R., and Schröder, W. "Acoustic perturbation equations based on flow decomposition via source filtering," Journal of Computational Physics Vol. 188, No. 2, 2003, pp. 365-398. Shen, H., and W. Tam, C. K. "Effects of jet temperature and nozzle-lip thickness on screech tones," AIAA journal Vol. 38, No. 5, 2000, pp. 762-767. Shen, H., and W. Tam, C. K. "Three-dimensional numerical simulation of the jet screech phenomenon," AIAA journal Vol. 40, No. 1, 2002, pp. 33-41. Wei, M., and Freund, J. B. "A noise-controlled free shear flow," Journal of Fluid Mechanics Vol. 546, 2006, pp. 123-152. Bridges, J., and Wernet, M. P. "Turbulence measurements of separate flow nozzles with mixing enhancement features," AIAA paper Vol. 2484, 2002. Marsden, A. L. "Aerodynamic noise control by optimal shape design." Stanford University, 2004. Marsden, A. L., Wang, M., Dennis Jr, J. E., and Moin, P. "Optimal aeroacoustic shape design using the surrogate management framework," Optimization and Engineering Vol. 5, No. 2, 2004, pp. 235-262. Marsden, A. L., Wang, M., Dennis Jr, J. E., and Moin, P. "Suppression of vortex-shedding noise via derivative-free shape optimization," Physics of Fluids (1994-present) Vol. 16, No. 10, 2004, pp. L83-L86. Economon, T. D., Palacios, F., and Alonso, J. J. "A coupled-adjoint method for aerodynamic and aeroacoustic optimization," AIAA Paper Vol. 5598, 2012, p. 2012.
Page 24 of 24
