A Study of Physical and Numerical Effects of Dissipation on Turbulent Flow Simulations

AbstrAct: The present work is primarily concerned with studying the effects of artificial dissipation and of certain diffusive terms in the turbulence model formulation on the capability of representing turbulent boundary layer flows. The flows of interest in the present case are assumed to be adequately represented by the compressible Reynolds-averaged NavierStokes equations, and the Spalart-Allmaras eddy viscosity model is used for turbulence closure. The equations are discretized in the context of a general purpose, density-based, unstructured grid finite volume method. Spatial discretization is based on the Steger-Warming flux vector splitting scheme and temporal discretization uses a backward Euler point-implicit integration. The work discusses in detail the theoretical and numerical formulations of the selected model. The computational studies consider the turbulent flow over a flat plate at 0.3 freestream Mach number. The paper demonstrates that the excessive artificial dissipation automatically generated by the original spatial discretization scheme can deteriorate boundary layer predictions. Moreover, the results also show that the inclusion of Spalart-Allmaras model cross-diffusion terms is primarily important in the viscous sublayer region of the boundary layer. Finally, the paper also demonstrates how the spatial discretization scheme can be selectively modified to correctly control the artificial dissipation such that the flow simulation tool remains robust for high-speed applications at the same time that it can accurately compute turbulent boundary layers.


INTRODUCTION
The present work is primarily interested in studying the effects of artificial dissipation and of certain diffusive terms in the turbulence model formulation on the capability of representing turbulent boundary layer flows.This interest comes from the fact that situations could arise in which one has a certain computational fluid dynamics (CFD) tool, for instance, developed elsewhere, and wants to apply this code to a particular application.It is clear that all decisions taken in the selection of the computational tool, from the choice of a specific turbulence model to numerical issues, such as the type of spatial discretization used, may have consequences on the quality of numerical results that might be obtained from the simulations.The research group, in the context of which the present effort is inserted, has recently experienced exactly this sort of situation.Therefore, an extensive study on the effects of artificial dissipation had to be performed in order to be able to correctly reproduce turbulent boundary layer flows.Similar issues with the effect of artificial dissipation terms on boundary layer flows have been previously addressed in the literature (Bigarella, 2007;Bigarella and Azevedo, 2012).However, this previous work is mostly concerned with centrally-differenced schemes and explicitly added artificial dissipation, whereas the present effort focuses on the artificial dissipation terms that arise from an upwind, flux-vector splitting-type discretization.Furthermore, the present study also addressed the decision to include, or not, some terms of the turbulence model formulation in the implemented code, since some of them, in many turbulence models, are computationally stiff and, hence, expensive.Thus, many investigators simply do not include these troublesome terms in their implementation of such particular model.
The present research group has been focusing on different aspects of CFD in the past years.For instance, the group maintains lines of work (Bigarella, 2007;Bigarella et al., 2007;Bigarella and Azevedo, 2009) aimed at creating new capabilities on numerical methods, multigrid techniques, and turbulence modeling.The study of such aspects is also a major issue in the present work.Previous effort was primarily geared towards the simulation of satellite launch vehicle (SLV) flows, which is one of the main interests of Instituto de Aeronáutica e Espaço.It resulted in a powerful Navier-Stokes solver, known as BRU3D, frequently used by the research group.
Additional effort at the group has also addressed the issue of high-order methods (Wolf, 2006;Wolf andAzevedo, 2006 and2007;Breviglieri, 2010;Breviglieri et al., 2010a and b), and successful applications of schemes such as essentially nonoscillatory (ENO) and weighted essentially nonoscillatory methods (WENO), and spectral methods have been demonstrated.Basso (1997) used preconditioning matrices to extend CFD codes for all speed applications.All these numerical technologies are applied on aeronautical and aerospace simulations, such as high lift and drag predictions, aerodynamics optimization, aeroacoustics, turbulent flows, and wind tunnel validation.The version of the BRU3D code of interest herein is a serial Navier-Stokes solver developed to simulate threedimensional (3-D) viscous turbulent flows over general aerospace configurations.The code presents different turbulent closures such as linear eddy-viscosity turbulence models, explicit algebraic Reynolds-stress models -EARSM (Wallin and Johansson, 2000;Hellsten and Laine, 2000), and Reynolds-stress models -RSM (Batten et al., 1999).A thorough study of flux computational schemes was also undertaken during the development of the code.Spatial discretization of the BRU3D code can be performed with the second-order accurate centered scheme of Jameson et al. (1981) and the Roe flux-difference splitting upwind scheme (Roe, 1981).Different artificial dissipation terms are also added for the Jameson centered spatial discretization (Jameson et al., 1981), such as the convective upwind split pressure (CUSP) scheme (Jameson, 1995a and b), scalar and matrix versions of switched seconddifference and fourth-difference models (Mavriplis, 1988;Turkel and Vatsa, 1990).The efforts of Bigarella (Bigarella, 2007;Bigarella and Azevedo, 2009) provided extensive expertise on turbulence modeling, which is a pacing item in CFD (Chapman, 1981).
On the other hand, previous work by Scalabrin and collaborators (Scalabrin, 2007;Scalabrin and Boyd, 2007;Schwartzentruber et al., 2007;Schwartzentruber et al., 2008) developed a numerical tool using upwind schemes, unstructured meshes, high-performance computing (HPC), and implicit integration for numerical simulations of weakly ionized hypersonic flows over reentry capsules.Such research has resulted in a very efficient numerical framework, called LeMANS, to simulate reentry flows over space capsules.These extreme conditions demanded the implementation of the Navier-Stokes equations coupled to nonequilibrium chemical reaction equations.The computation of these sets of equations requires very fine meshes, which makes impractical the use of serial algorithms.Therefore, message passing interface (MPI) protocols were implemented to parallelize LeMANS and, then, reduce computational costs.One should understand that high-fidelity CFD solvers have very complex algorithms and, hence, their parallelization involves advanced numerical and computational issues.Thus, a full parallel code, as LeMANS, is always welcome.Furthermore, since LeMANS already incorporates several programming issues and highspeed flow physics models, it seems to be a more suitable code for continued developments in the future.
In this context, an important motivation of the present work is to take full advantage of all scientific technology concerning turbulence modeling, boundary condition, and initial condition treatment implemented in the BRU3D solver in order to extend the LeMANS code for the research group needs, more specifically, parallel turbulent flow simulations for high dissipative spatial discretization, which are strongly recommended for high-speed configurations, such as reentry flows.However, such highly dissipative methods can strongly deteriorate boundary layer flow predictions.It is clear that the challenge is to selectively modify the discretization scheme in order to correctly control the artificial dissipation such that the flow simulation tool remains robust for high-speed applications, at the same time that it can accurately compute turbulent boundary layers.
Therefore, the present work selected the well-known problem of a turbulent flow over a flat plate in order to address the dissipation issues that are relevant in this case.Two major aspects are investigated, namely the problem of excessive artificial dissipation of the upwind spatial discretization scheme, and the inclusion of numerically stiff cross-diffusionlike terms in the formulation of the turbulence model.For the present study, the Spalart-Allmaras turbulence model was selected, primarily because it is probably the most widely used turbulence closure for realistic aerospace applications at the time the study was carried out.
This study considers the case of freestream Mach number equal to 0.3, because there are experimental and other independent computational results available.Moreover, since all computational codes here considered implement a compressible formulation, there are no issues with the incompressible limit at such Mach number.The paper demonstrates that the excessive artificial dissipation automatically generated by the original spatial discretization scheme can deteriorate boundary layer predictions.The paper also demonstrates how the spatial discretization scheme should be selectively modified to correctly control the artificial dissipation.Finally, the results show that the inclusion of Spalart-Allmaras model cross-diffusion terms is primarily important in the viscous sublayer region of the boundary layer.

THEORETICAL FORMULATION
The formulation used in the present work is based on the Reynolds-averaged Navier-Stokes set of equations, also known by the CFD community as RANS equations.They are obtained by filtering the Navier-Stokes set of equations.This process filters the fluctuation part of the fluid and maintains only the mean contribution.The filtered information needs to be recovered somehow.Turbulence models are applied to the RANS formulation to recover the effect of the fluctuating part.The levels of turbulence modeling are also discussed in this section.The most used filtering processes are based on the time, space, and ensemble averages.The filtering based on the time average is the most used for steady state applications and it is the one applied in the present work.For the sake of simplicity, the filtering process is not discussed in this work.The reader can find further details on the work of Bigarella (2007) and Junqueira-Junior (2012).
The filtered compressible Reynolds-averaged Navier-Stokes equations are written in the vector form as The conserved variables vector, Q, the inviscid flux vector, Fe, and viscous flux vector, Fv, are given by where ρ stands for the density, v = {u, v, w} is the velocity vector in Cartesian coordinates, p is the static pressure, τ is the viscous stress tensor, q H is the heat flux vector, e is the total energy per unit volume, and β i is given by The î x , î y , and î z terms are the Cartesian coordinate orthonormal vector basis.The { } ̅ terms are the averaged and weighted averaged properties.Therefore, it is very important to emphasize that field forces, such as gravity, are neglected here.
Other equations are necessary in order to close the system of equations given by Eq. ( 2), which are called constitutive relations.The first constitutive equation presented to close the Navier-Stokes set is known as the equation of state.This equation considers the perfect gas law, and it is written as in which the mean total energy per unit volume, e ̅ , is given by and i stands for the internal energy, defined as in which T stands for the mean static temperature and C v is the specific heat at constant volume.The heat flux from Eq. ( 2) is obtained from the Fourier law for heat conduction, and it is given by in which γ is the ratio of specific heats and Pr is the Prandtl number.Typically, for air, it is assumed that γ = 1.4 and Pr = 0.72.C p is the gas specific heat at constant pressure, and µ is the dynamic molecular viscosity coefficient, calculated as a function of the temperature by the Sutherland law equation (Anderson, 1991), written as In the above equation, S = 110K, and µ ∞ is the dynamic molecular viscosity coefficient of the fluid at temperature T ∞ .The components of the viscous stress tensor, for a Newtonian fluid, are given by in which δ ij stands for the Kronecker delta.
All the terms marked with the superscript { } t , in Eq. (2), appear after the time filtering processes.These terms carry important turbulent information and need to be modeled.The turbulence closures are responsible for representing them.There are two major families of turbulence models for the RANS equations, the first and second order closures.The present paper focuses on the first order closures, more specifically, on the Spalart and Allmaras (1992) model, which is an one-equation closure.This model was chosen because it is, by far, the most used turbulence model for realistic aerospace applications.Furthermore, the research group has already achieved good results using it on previous applications (Bigarella, 2007;Bigarella et al., 2007;Bigarella and Azevedo, 2009).The Spalart-Allmaras turbulence closure is a partial differential equation, which models the turbulent eddy viscosity transport.The theoretical and numerical formulations of the turbulence closure are discussed in details in the forthcoming sections.

NUMERICAL FORMULATION spAtiAl discretizAtion
The spatial discretization used is the first aspect to be discussed in this section, starting with the finite volume formulation and followed by the flux calculations.For the sake of simplicity, from here, all the averaged terms are written without the { } ̅ notation.

Finite volume formulation
The finite volume formulation is a numerical method applied to represent and evaluate partial differential equations.It is applied by the CFD community to find the solution of the RANS equations, Eq. ( 1).The method is obtained integrating the flow equations for each control volume within a given mesh, Considering a cell-centered formulation, V i is a determined cell of the given grid.After the integration, it is possible to apply the Gauss theorem over Eq. 10, resulting in in which S i is the outward-oriented area vector and it is defined as Considering the mean value of the conserved variables within the i-th control volume, one can write the first term of Eq. ( 11) as The second term of Eq. ( 11) can be written as the sum of all faces of a cell in which the k subscript is the index of the cell face, and nf indicates the number of faces of the i-th volume.Finally, the RANS equations discretized with a finite volume approximation is given by For this formulation, the fluxes are computed at the faces of the control volume, and the conserved variables are computed in the cell.

Inviscid flux calculation
The inviscid fluxes are calculated using a method based on a classical flux vector splitting formulation, the Steger-Warming scheme (Steger and Warming, 1981).The formulation implemented to compute the inviscid fluxes is explained here.This method is an upwind scheme that uses the homogeneous property of the inviscid flux vectors to write where Fe n is the normal flux at the k-th face, and A is the Jacobian matrix of the inviscid flux that can be diagonalized by the matrices of its eigenvectors from the left and from the right, L and R, respectively, as and Λ is the diagonal matrix of the eigenvalues of the Jacobian matrix.The A matrix can be split into positive and negative parts as The splitting separates the flux into two parts, the downstream and the upstream fluxes, in relation to the face orientation as where the cl and cr subscripts are the cells on the left and right sides of the face.The split eigenvalues of the Jacobian matrix are given by In order to avoid sudden sign transitions, as illustrated in Fig. 1, the split eigenvalues receive a small number, ϵ, turning Eq. ( 20) into Numerical studies performed in the present paper indicated that this flux vector splitting is too dissipative and it can deteriorate the boundary layer profiles (MacCormack and Candler, 1989;Junqueira-Junior et al., 2011).To avoid such issue a pressure switch is implemented to smoothly shift the Steger-Warming scheme into a centered one.Then, the artificial dissipation is controlled and the numerical stability is maintained as presented in in which The switch, w, is given by Therefore, for small ∇p, w = (1-w) = 0.5, the code runs with a centered scheme, and for larger values of ∇p, w = 0 and (1-w) = 1, the code runs with the Steger-Warming scheme.For Eq. ( 24), it is suggested α = 6, but some problems may require larger values (Scalabrin, 2007).The applied formulation was originally created with interest on studying flows over reentry capsules.For this particular case, with very strong shock waves, it is very common to find solutions with nonphysical numerical structures such as carbuncles (Ramalho et al., 2011).To prevent such numerical problems, artificial dissipation has necessarily to be added to the method.The dissipation was included into the split eigenvalues, Eq. ( 21), using the ϵ factor, which is given by: Junqueira-Junior, C., Azevedo, J.L.F., Scalabrin, L.C. and Basso, E.
where d k is the distance of the k-th face to the the nearest wall boundary, d 0 is set by the user and must be smaller than the boundary layer thickness and larger than the shock stand-off distance, mk → is the normal vector of the nearest wall, and nk → is the normal vector to the k-th face.Equation ( 25) applies the term (1 -|nk → •mk → |) to decrease the ϵ value at the faces parallel to the wall inside the boundary layer (Scalabrin, 2007).This artificial dissipation model has shown an important role in the prediction of boundary layer profiles, and several tests were performed in the present work with the aim of understanding its behavior (Junqueira-Junior et al., 2011).

Viscous flux calculation
The viscous terms are based on derivative of properties on the faces.To build the derivative terms, two volumes are created over the face where the derivative is being calculated.At the center of each new volume, the derivative is calculated using the Green-Gauss theorem.This computation is used to find the derivative at the desired face.
A two-dimensional (2-D) example is used in this section to better explain the derivative calculation.Considering the two cells, S1 and S2, in Fig. 2, two new cells, S3 and S4, are created using node points, P1 and P3, and cell centered points, P2 and P4, to calculate the derivative on the face 1-3.
The properties at the faces are calculated using simple averages.For the example in Fig. 2, they are given by Using the averaged properties, Q 12 , Q 23 , Q 13 , Q 34 , and Q 41 , with the normal vectors, n 12 , n 23 , n 13 , n 34, and n 41 , and the surface of the faces, S 12 , S 23 , S 13 , S 34 , and S 41 , as illustrated in Fig. 3, it is possible to calculate the derivative at the points P7 and P8 using the Green-Gauss theorem (Jawahar and Kamath, 2000), which is applied to a scalar and relates the volume integral of the gradient of its area integral over the boundary as Considering ∇Q as a constant over the cell, Eq. 27 yields in which ∇Q is the constant cell-centered gradient.Using the derivatives in the cells S3 and S4, the derivative at face 1-3 is computed using The derivative computation for other types of element, 2D or 3D, is straightforward.

time integrAtion
Simulations of turbulent flows can become very stiff.Such stiffness substantially limits the use of large time steps.One classical solution is the use of implicit time integration.This work applies an implicit integration based on the backward Euler method, which is given by One can linearize the residue at time n+1 as a function of properties at time n.
From the spatial discretization, the inviscid terms can be written as with and The formulation described assumes which is not true and can decrease the numerical stability of the method.Then, the true Jacobian matrices of the split fluxes were implemented in place of A ± to calculate the implicit operator.Issues involving the true Jacobians matrices have a great importance in the context of numerical stability for computational methods.The reader with interest in this subject must look further in references Anderson et al. (1986) and Steger and Warming (1981), and in chapter 20 from Hirsch (1990).
The viscous terms can be written in the same form as The viscous Jacobian matrices are represented by B. The splitting of these matrices is written as The true Jacobian matrices, for A ± and B ± , can be found in the work of Scalabrin (2007).One can write the system as It is, then, possible to write with ( ) As the code is an unstructured solver, this system of equations results in a sparse block matrix, where each block is a square matrix of size equal to the number of equations to be solved in each control volume.The solution of such system is typically very expensive and, depending on the size of the mesh, it is not even practical.A less expensive implicit method is applied in the present paper, the point implicit integration (Gnoffo, 2003;Venkatakrishnan, 1995;Wright, 1997).
The main idea of the point implicit integration is to move all the off-diagonal terms to the right hand side and solve the resulting system iteratively, i.e., It is assumed that ∆Q n+1,0 = 0 and four iterations are taken in the process as suggested in the literature (Wright, 1997).The sparse linear system illustrates the point implicit method: Each , in the sparse matrix, is a block matrix.The time step is computed by in which CFL (Azevedo, 1988) is a parameter set to ensure stability of the time integration method, l is the size of the cell and || -v || + a is the largest wave speed in the cell (Scalabrin, 2007).

second-order extension of inviscid fluxes
The monotone upstream-centered scheme for conservation laws, known as MUSCL approach (van Leer, 1979), is used in order to obtain second-order extension for the inviscid fluxes calculation.This section presents the classical formulation of the MUSCL approach and an extension for unstructured grids.

MUSCL approach
A first-order spatial discretization is equivalent to represent the numerical approximation of the solution as a piecewise constant.The MUSCL idea is to use a linear approximation to achieve a second-order space discretization.A linear solution is exactly resolved, which generates a truncation error of the order ∆x 2 .In order to represent the conservation laws, the discrete state variables express the average state within the cells.Then, the linear approximation has to average out these values (Hirsch, 1990).
One can consider the Taylor linearization as a onedimensional (1-D) local representation, Fig. 4, valid in a given cell "i", at a determined instant The use of backward and forward derivatives provides One can rewrite these terms at the same faces as Using limiters to treat correctly the discontinuities, the MUSCL scheme can be written as where r is the ratio of consecutive variations, given by and ψ(r) is a limiter function.Th ese functions are extremely important in the context of high-order methods.However, this discussion is not in the scope of the present work.
One can fi nd further information about limiter functions in chapter 21 from Hirsch (1990).Th ere are several limiter functions available for highorder methods of CFD applications.Two of these functions are implemented herein, the van Albada limiter (van Albada et al., 1982), given by and the minmod limiter (Hirsch, 1990), written as Th e simulations performed in the present work used mainly the van Albada limiter function (Hirsch, 1990).

Second-order extension for unstructured grids
Th e MUSCL variable extrapolation for 2-D or 3-D is straightforward for structured meshes.However, it is not very simple for unstructured solvers.Th e approach applied here is based on the work of Batina (1990) and Bibb et al. (1997).
Here, the stencils are created using only cell-centered values.Th e points "i" and "i+1" are, respectively, the center of the cell at the left and right of the face.Th e other two points, "i-1" and "i+2", are defi ned by a stencil search.
Th e search for cells "i+2" and "i-1" is limited to the ones that share at least one node with the volumes "i+1" and "i", respectively.Th e selected "i+2" point is the one that has the maximum positive value of the dot product between the face normal and the normalized vector joining the face centroid to cell centroid.Using the same principle, the chosen "i-1" point is the one with the maximum negative value of the dot product.Figure 5 is extracted from Scalabrin (2007), and illustrates the search for the "i+2" point in a given mesh.
It is possible to observe that, when the stencil search is applied to unstructured grids, the distance between the "i+1" and "i+2" points can be diff erent from the distance between the "i" and "i+1" points.Th erefore, a correction is applied on the limiter, as given by where "a" is the distance between the "i+1" and "i" points, "b" is the distance between the "i+2" and "i+1"points, and "c" is the distance between the "i" and "i-1" points.and Allmaras, 1992 and1994).The one equation model is originally written as

TURBULENCE MODELING
The kinematic eddy viscosity is defined as The classical approximation used by the CFD community to write this single equation model in a conservative form for compressible flows consists in multiplying Eq. ( 59) by ρ, which yields The new variable to be solved is μ ̃, defined μ ̃ = ρν ̃.It is very important to point here that and are not fluid properties, but flow properties.The closure coefficients and auxiliary relations are given by in which d is the distance from the closest surface and Ω stands for the terms present in the anti-symmetric part of the mean velocity gradient field.
In general terms, turbulence is modeled by transport equations in order to represent turbulent properties being carried by the mean flow.These transport equations have advection, diffusion, source, production, and destruction terms, such as the ones indicated in Eq. 63: or using the definition of total derivative: in which q is the turbulent property, Pq is the production term, Sq is the destruction term, Dq is the diffusion term, and CDq is the cross-diffusion-like term.The second term in the left-hand side is the advection term, here defined as Cpq.
In the context of the SA turbulent closure, these previously discussed terms are written as The SA turbulence model has been extensively used by the CFD community for 3D compressible flow with very good agreement to experimental data for many relevant applications (Spalart andAllmaras, 1992 and1994;Bigarella, 2007;Bigarella and Azevedo 2009).

numericAl implementAtion
In order to discretize the SA turbulence model equation in a finite volume context, it is necessary to integrate Eq. ( 64), yielding Using the mean property definition and the Green-Gauss theorem, it is possible to rewrite the finite volume equation The production, Pi, destruction, Si, and the cross-diffusion, CDi, terms are considered constant in the i-th volume.Summation over the faces forming the i-th volume needs to be performed in order to calculate the surface integral, The S k term is the outward facing normal area vector.Hence, Eq. ( 69) can be written as After the operations, the SA equation is written as where the residue is It is important out here that this formulation needs property values on the cell faces to compute the summation terms and property values on the cell centers for the computation of the source terms.

mesh requirements for turbulent simulAtions
Turbulent simulations require that a sufficiently refined mesh at the wall is provided.In particular, the parameter typically used to measure such refinement is denoted y + , which is the dimensionless distance from the point to the nearest wall (Tennekes and Lumley, 1972).Schlichting (1978) defines y + as a Reynolds number based on the friction velocity, in which y is the distance to the nearest wall, and u τ is the friction velocity, written as

√ ¬ ¯ (74)
Here, τ w is the wall shear stress.For the correct solution of wall-bounded turbulent flows, with turbulence models solved for up to the wall, it is usually required that the first point off the wall be located so as to satisfy y + ≤1 (Bigarella, 2007).However, obviously, it is not possible to know τ w value before achieving the solution of the simulations.Then, in the present work, the position of the first point off the wall is determined using an empirical approximation (Bigarella, 2007;van der Burg et al., 2000): where y + is the desired value set by the user, and Re L is the Reynolds number based on the reference length, L.

BOUNDARY CONDITIONS
The boundary conditions are implemented using ghost cells.The solver creates the ghost cells to hold properties that satisfy the correct flux calculation at the boundaries.The implementation assigns properties that satisfy the Euler boundary conditions to calculate the inviscid fluxes, and properties that satisfy the Navier-Stokes boundary conditions for calculating the viscous fluxes.Therefore, the ghost volumes store two different types of fluxes for the correct computation of the RANS equations.

Wall and symmetry boundary conditions
Ghost cells are applied for the implementation of boundary conditions.The ghost cells hold the properties in the same manner to calculate the inviscid fluxes at the wall and at the symmetry boundaries.Mass and energy fluxes should yield zero, and the momentum flux is equal to the pressure flux.This is accomplished by setting the normal velocity component to the boundary face zero.To simplify, the properties at the left side of the boundary face, which is the interior domain, are rotated to the face coordinates using As the left side of the boundary is the interior domain, the right side of the boundary is the ghost cell.
In the present section, the turbulent variable from the Spalart-Allmaras model is included in the Q vector, which now is written as T . (77) In Eq. ( 76), is the rotation matrix given by, R= 1 0 0 0 0 0 0 n x n y n z 0 0 0 t x t y t z 0 0 0 rx ry rz 0 0 0 0 0 0 1 0 0 0 0 0 0 1 ( One can write in the matrix form as in which W is the inviscid wall matrix given by W= 1 0 0 0 0 0 0 − 1 0 0 0 0 0 0 1 0 0 0 0 0 0 1 0 0 0 0 0 0 1 0 0 0 0 0 0 1 Therefore, the boundary condition can be written as in which the R -1 matrix is given by It returns the properties to the Cartesian coordinate frame.

Nonreflecting farfield boundary condition
The concept of Riemann invariants (Long et al., 1991;Bigarella, 2007) is implemented to a nonreflecting farfield boundary condition at subsonic speeds.These invariants are derived from the characteristic relations for the Euler equations.The formulation at the boundary is given by in which v n is the normal velocity component given by The subscripts ∞ and int represent the property at the freestream and in the interior domain, respectively.The normal velocity component and the speed of the sound at the boundary face can be written as in which the f subscript represents the property at the farfield computational surface.It is possible to write the velocity for a subsonic exit, 0 < v n int / a int < 1, using the tangential velocity components of the interior and the definition of normal velocity.
The other properties are given by For a subsonic entrance boundary, -1 < v n int / a int < 0, one should extrapolate the freestream properties as For supersonic fl ows, it is not necessary to use the concept of Riemann invariants, because information propagates in only one direction in inviscid supersonic fl ows.Th erefore, zeroth-order extrapolation is used for it yields a cheaper computation.Hence, for a supersonic exit boundary, v n int / a int > 1, the properties are extrapolated from the interior of the domain as On the other hand, for an entrance boundary, v n int / a int < -1, the properties are extrapolated from the freestream as The freestream properties are previously provided by the user.

Ghost cells
One obtain the properties at the boundary face using the Riemann invariants.However, in order to obtain the boundary conditions these properties need to be computed in the ghost cells.It is possible to use an average to calculate the properties in the ghost volume as in which the subscript gh stands for ghost cell.
Riemann invariants are derived for farfi eld boundary conditions.It is strongly recommended to avoid their use for other fl ow situations, such as entrance and exit boundary conditions for internal fl ow cases (Bigarella, 2007).Th e implementation has shown to be very sensitive when the farfi eld boundary is set close to solid walls.

viscous boundAry conditions
For adiabatic boundary condition, it should be assumed that the heat conduction through the boundary face yields zero, q H wall •n = 0, hence Th erefore, one can state that Tcr =T wall =T cl . (94) In order to satisfy wall pressure condition (Schilichting, 1978), and then write Th us, it is possible to extrapolate the density from the interior Th e Cartesian components of the wall velocity, u wall, v wall and w wall , are set by the user.It is possible to use these velocity components and the average procedure to calculate the velocity components, u cr , v cr and w cr , in the ghost cells.
Th e turbulent property in the ghost cell is set to -μ cl in order to force μwall = 0. Symmetry, nonrefl ecting farfi eld, and zeroth-order extrapolation boundary conditions use the same procedures applied on the computation of the corresponding inviscid boundary conditions.

implicit boundAry conditions
Implicit boundary conditions are necessary in order to obtain a truly implicit time-marching method.Th e use of explicit boundary conditions can substantially limit the stability of the numerical method in the marching procedure for the solution convergence.

bAsic implicit formulAtion
A simplified form of the implicit equation, Eq. ( 38), is written in this section in order to detail the implementation of implicit boundary conditions for flux vector splitting schemes: In such equation, the repeated k index in the second term in left hand side of the equation indicates summation over all the k faces of the control volume.Equation ( 101) is written only to present the relation between an internal cell, cl, and a boundary cell, cr, k.In the original formulation, Eq. ( 38), the cl-th cell has contributions from other faces, which may or may not be boundaries.
As presented in the beginning of this section, the ghost cells hold different values for inviscid and viscous calculations.Hence, using the splitting definition, presented in section numerical formulation, Eq. ( 101) can be written as The contributions of the boundary face can be expressed in terms of the internal cell corrections as Hence, Eq. ( 102), can be rewritten as The viscous Jacobians are calculated using primitive variables, and the corrections are set for the primitive variables and applied directly at the calculation of the viscous Jacobians.The matrix B + k+ already includes the contribution from the boundary.The A + k+ , A - k-, B + k+ and B - k-are presented in the work of Scalabrin (2007).

Inviscid boundary conditions
The matrix for an inviscid wall or for a symmetry boundary is the same presented in Eq. ( 82), The matrices are applied to ∆Q for the implicit boundary condition according to Eq. ( 104).
It is possible to use the identity matrix to represent the zeroth-order extrapolation as given by For the purpose of developing the implicit boundary condition, the farfield variables are considered constant.Hence, ∆Q at a farfield boundary is given by where 0 stands for the zero matrix.

Viscous boundary conditions
The viscous Jacobians are created using primitive variables.Therefore, the implementation of implicit viscous matrices is performed using primitive variables.They are applied directly at the calculation of the Jacobians matrices as The adiabatic wall implicit boundary condition is calculated using The same procedure used for the inviscid zeroth-order extrapolation is used here, hence A Study of Physical and Numerical Effects of Dissipation on Turbulent Flow Simulations Th e viscous farfi eld matrix is the same one used for the inviscid farfi eld boundary condition, which is given by in which 0 stands for the zero matrix.

FLOW SIMULATION RESULTS AND DISCUSSION
Th e test case used in this work is the simulation of turbulent subsonic fl ow over a fl at plate in absence of streamwise pressure gradients.Th e work studies the addition of diff erent levels of artifi cial dissipation, at diff erent distances of the wall, in an attempt to fully understand the eff ects of the high dissipative upwind spatial discretization over turbulent dimensionless boundary layer profi les and over the friction coeffi cient.All results are compared to analytical and experimental data.Th e Mach number at the freestream condition is set to M ∞ = 0.3, and the Reynolds number based on 1 m long fl at plate is 7.2 million.Th erefore, the fl ow can be considered as a turbulent compressible fl ow.Th e computational domain covers a 3 m high and 7 m long region.Th e wall fl at plate is located between 3< x <4.As illustrated in Fig. 6, the applied boundary conditions are as follows: in the lower portion of the computational domain, symmetry is used between the entry and the fl at plate, and between the trailing edge and the outlet; an adiabatic wall is assumed over the fl at plate surface; and farfi eld Riemann-type boundary conditions are used for all the other boundaries.Th e simulations are performed using two grids, one with 60,000 cells and another with 210,000 cells, as illustrated in Fig. 7.Both meshes have y + = 0.5 at the fi rst grid point off the wall, but the mesh with 210,000 cells has more points near to the wall.
Th e dimensionless boundary layer profi les obtained by the simulations are compared to the analytical formulation given by the law of wall (Schilichting, 1978), which is written by y + < viscous sublayer u + = y + 30 5 < y + <300 log-law region u + = 1 0.41 ln (y + )+5.5 .
Th e distribution of the local skin friction coeffi cient, c f , which is written as is compared to experimental data (Coles and Hirst, 1969), and to the analytical formulation (von Karman, 1934) given by  where Re x is local Reynolds number based on a distance from the current position to a reference point, and is given by in which U ∞ is the freestream velocity and x is the distance to a reference position.The upwind Steger-Warming (Steger and Warming, 1981) scheme is too dissipative and, therefore, implementing a turbulence closure to treat the turbulent flows is not enough to provide accurate results.The dissipative terms, present in the spatial discretization and in the turbulence model equations, need to be carefully treated.The effects of the spatial discretization on the results are analyzed by forcing the switch term, Eq. ( 24), to the classical Steger-Warming scheme, w = 0, and to the centered scheme, w = 0.5.The artificial dissipation term, Eq. ( 25), and the Spalart-Allmaras cross-diffusion-like term, Eq. ( 65), are added in different quantities and at different positions of the boundary layer profile in order to study their influence on the skin friction coefficient and on the dimensionless boundary layer.Table 1 presents the validation studies performed for this test case.It enumerates the simulations and indicates the chosen spatial discretization; the distance, d 0 , for the artificial dissipation; the amount of artificial dissipation at a given region; the choice of the turbulence model; and the choice of the cross-diffusion-like term.
The first three simulations are an introductory study of the numerical issues concerned in this work, and the results are illustrated in Fig. 8. Simulation No. 0 is performed without any turbulence modeling.It provides an underestimated c f distribution over the flat plate and a boundary layer profile completely different from the analytical turbulent one.For case No. 1, as an initial approximation, the cross-diffusionlike term of the SA equation, the most expensive source term, is not included in the formulation and the inviscid flux is calculated using the original scheme of the code, as presented in Eqs. ( 23) and ( 24).The results with the turbulent model have shown to be better than those without any turbulence model.The boundary layer profile matches the analytical results at the viscous sublayer, however the excessive artificial dissipation, provided from the spatial discretization, deteriorates the boundary layer profile at the log-law zone.The skin friction distribution is still underestimated, which indicates that the computation of shear-stress tensor, τ w , at least at the wall, is not correct.Simulation No. 2 uses the centered scheme without any artificial dissipation.The results present the same shortcomings as the previous ones discussed here.However, the results for test case No. 2 have shown to be closer to the analytical data.
The solution of a nonlinear partial differential equation, spatially discretized without the addition of artificial dissipation, is numerically unstable (Lomax et al., 2001), due to the frequency cascade phenomenon.Test case No. 2 achieved convergence of the solution because it provides a very simple geometry and due to the presence of the viscous fluxes, which are dissipative by nature.The numerical solver set up, which uses the centered scheme without the addition of artificial dissipation, is unstable and cannot be used for complex configurations.
Simulations including test cases No. 3 to No. 19 are performed using the centered scheme.Artificial dissipation is added to the formulation in different quantities, and at different locations within the boundary layer profile.The local error is measured as the difference between the analytical, and is given by  The c f distributions are underestimated and presented nonphysical oscillations.
The CD terms are implemented to test the full SA turbulent equation.Simulations No. 7, 8, 9 and 10 are performed to study the effects of the cross-diffusionlike terms and of the artificial dissipation on the results.These simulations use a centered scheme and the artificial dissipation in the flow direction in different quantities and positions off the wall.Figures 11 and 12 demonstrate that the numerical oscillations, discussed in the previous simulations without the CD terms, are no longer present in the skin friction coefficient distribution along the flat plate.However, the shape of the boundary layer profile in the log-law zone is very different from the analytical one, with errors superior to 100% of the analytical value.Once again, the addition of artificial dissipation in the flow direction does not significantly affect the boundary layer profile.
Simulations No. 11, 12, 13 and 14 are carried out using a centered scheme with artificial dissipation added in all directions, in different quantities and positions of the boundary layer profile.From the results, it is clear that the addition of artificial dissipation improves the boundary layer results  the analytical and experimental data, with errors of ≈ 5% for the log-law zone of the boundary layer profile and for the skin friction coefficient distribution over the flat plate.These simulations are carried out using the mesh illustrated in Fig. 7(b).Results are accurate and correspond to aerospace industry needs.On the other hand, the sublayer solutions presented maximum error of ≈10% for the simulations No. 18 and 19.These sublayer errors are not excellent but plausible, since there is good agreement of the log-law zone and the c f distribution with the literature.Therefore, results with good agreement with analytical (von Karman, 1934) and experimental (Coles and Hirst, 1969) data are achieved when one controls the artificial dissipation in the log-law zone, using 1.25 (a k +|u k |) < ϵ < 1.5 (a k +|u k |), and turn off the cross-diffusion-like-terms above the sublayer zone.experimental data.The very dissipative spatial discretization can deteriorate the boundary layer profile and the distribution of the local skin friction coefficient.Therefore, the effort is undertaken in order to adapt the numerical scheme for turbulent applications.The switch term is changed to force the upwind scheme to a centered one in boundary layer regions.Different artificial dissipation terms are added to improve the numerical stability of the method, shift the shape of the boundary layer, and change the distribution of the local skin friction coefficient.The cross-diffusion-like terms of the turbulent equation have also shown to possess an important role on the results.The study has shown that it is possible to achieve good results for turbulent flows using a discretization based on the Steger-Warming flux vector splitting scheme, provided that an appropriate control of the intrinsic artificial dissipation terms is implemented.The artificial dissipation must be correctly applied only in the log-law zone, and the cross-diffusion-like terms only in the viscous sub-layer region.
In the context of validating the code for turbulent flows application, it is possible to state that good results were achieved for the 2D flat plate flow simulation using  spatial discretization based on the Steger-Warming scheme.However, the good agreement is dependent upon the addition of artificial dissipation and CD terms at the correct position in the boundary layer profile.In the present work, these modifications are explicitly performed to provide good results for a predefined local Reynolds number.It is very important to adapt the switches in the computational implementation of the numerical scheme in order to add artificial dissipation in the correct position in the boundary layer profile.Similarly, the cross-diffusion-like terms must be switched on only in the sublayer region.Hence, the turbulence model implementation must be performed in such a way that the code could automatically recognize the correct position in the dimensionless boundary layer profile and switch on/off the CD terms and the artificial dissipation.This can certainly be accomplished using an appropriate scaling of boundary layer quantities and the normalized distance to the wall, which is already computed by the turbulence model.However, these are implementation issues, which are beyond the scope of the present work.Clearly, such implementations would be necessary in order to handle more complex configurations.More complex applications, in the sense, for instance, of treating flows with more obvious compressibility effects, are not an issue, since the turbulence model and the present CFD tool were both originally developed for compressible flows.

Figure 3
Figure 3. 2-D example of a derivative calculation.
the spAlArt-AllmArAs turbulence model Th e Spalart-Allmaras (SA) closure(Spalart andAllmaras, 1992 and1994) is a one-equation, linear eddy-viscosity turbulence model.It solves one transport equation for the modifi ed eddy viscosity coeffi cient, ν .Th e model uses eight closure coeffi cients and three closure functions derived along intuitive and empirical lines, relying heavily on calibration by reference to a wide range of experimental data(Spalart

Figure 7 .
Figure 7. Visualization of the meshes used in the zeropressure-gradient fl ow over a fl at plate.(a) Mesh with 60,000 volumes; (b) Mesh with 210,000 volumes.

Figure 6 .
Figure 6.Domain used for the zero-pressure-gradient fl ow over a fl at plate.

Figure 8 .
Figures 9 and 10 present the results and errors for simulations No. 3, 4, 5 and 6, using the SA model without the crossdiffusion-like term.These figures illustrate that the addition of artificial dissipation, in the flow direction, does not affect substantially the boundary layer profile.As with simulation No. 1, the results are in good agreement with the analytical ones at the sublayer zone and overpredicted the log-law zone by about 20% of the analytical value.

Figure 9 .
Figure 9.Comparison of the dimensionless turbulent boundary layer and skin friction coefficient to analytical and experimental data(Coles and Hirst, 1969) for simulations 3, 4, 5, and 6.
Figures 13 and 14 illustrate that this simple modification significantly affects the results.The simulation results with ϵ < 0.3 (a k +|u k |) have the same behavior as the other simulations presented so far, a good agreement of the sublayer region with analytical data, overpredicted log-law boundary layer, and an underestimated c f distribution.The simulations with ϵ > 0.3 (a k +|u k |) have shown good agreement of the skin friction coefficient with the analytical and experimental data, and an underestimated boundary layer profile.

Figure 10 .
Figure 10.Boundary layer error distribution for simulations 3, 4, 5, and 6.(a) Error in the log region; (b) Error in the viscous sublayer.

Figure 11 .
Figure 11.Comparison of the dimensionless turbulent boundary layer and skin friction coefficient to analytical and experimental data(Coles and Hirst, 1969) for simulations 7, 8, 9, and 10.

Figure 15 .
Figure 15.Comparison of the dimensionless turbulent boundary layer and skin friction coeffi cient to analytical and experimental data (Coles and Hirst, 1969) for simulations 15, 16, and 17.

Figure 17 .
Figure 17.Boundary layer and skin friction coefficient error distribution of simulations 18 and 19.(a) Error in the boundary layer log region; (b) Error in the viscous sublayer; (c) Error in the skin friction coefficient.

Figure 16 .
Figure 16.Comparison of the dimensionless turbulent boundary layer and skin friction coefficient to analytical and experimental data (Coles and Hirst, 1969) for simulations 18 and 19.

Table 1 .
Zero-pressure-gradient flow over flat plate simulations.