- Email: [email protected]
A numerical method for solving the compressible Navier-Stokes equations
IMPACT
OF
(‘OMPU
I INCi
IN
S(‘IEU(‘E
hNI1
l:NuCilNttKING
1,
64-02
( 1Y8Y )
A Numerical Method for Solving the Compressible Navier-Stokes Equations S. BOIVIN
D6partemeni de Mathimatiques, Universitk Laval, Laval, Quebec, Canada GlK 7P4. and INRIA. Menusin. 78153 Le Chesnay, France S. Boivin, A Numerical Method for Solving the Compressible Navier-Stokes Equations, IhlPACT ofComputing in Science and Engineering 1, 64-92 (1989).
A new computer code solving the compressible Navier-Stokes equations is described. The code solves the two-dimensional (2D), time-dependent equations using the finite-element method. Filtering methods are used on gradients of convection terms to ensure stability. Two turbulence models are integrated into the code to compare, evaluate, and ultimately improve model performance. The results of numerical experiments for the flow around an NACAO012 airfoil and a simplified version ofthe Hermes space shuttle are presented to show the possibilities of the methods used.
@I1989
Academtc
Press. Inc.
INTRODUCTION
Numerical simulations of flows around airplanes and space vehicles at high Mach and Reynolds numbers are still a challenge to the numerical analyst. Most of the difficulties arise from the simultaneous presence in the flow of complex physical phenomena of very different scales. The most frequently encountered are shocks, boundary layers (both laminar and turbulent), wakes, and various interactions between these phenomena. Throughout this paper, we hypothesize that all these phenomena are adequately modeled by the compressible Navier-Stokes equations. In particular, we assume that these equations are valid for highly turbulent flows. Under this hypothesis, our main problem is how to solve the Navier-Stokes equations without overdamping some of these phenomena. The strategy of resolution presented is a step in that direction and is based on the following key ideas: 64 0899-8248/89 $3.00 Copyngbt
0
1989 by Academac
All rights of reproduction
Press. Inc.
m any form resewed
COMPRESSIBLE NAVIER-STOKESEQUATIONS
65
first, being unable to represent the smaller physical scales, we take account of them through a viscosity computed with appropriate turbulence models, ??
second, unresolved frequencies following a shock must be damped by a numerical viscosity taking into account the physical viscosity, ??
third, with no reliable a priori information about the boundary layer, we seek to resolve it directly, or eventually, by the use of special elements. ??
Although very. appealing, no wall law is used because even if the region where viscous effects are important is small, its influence is felt everywhere in the flow. The formulation, schemes, and method of solution were also chosen to achieve some particular goals. In fact, it is our opinion that a good twodimensional (2D) Navier-Stokes code should: (i) be directly generalizable to 3D case, (ii) include a turbulent closure model with transport equations, (iii) be independent
from the mesh,
(iv) enable control of the numerical diffusion, (v) have solid theoretical support. In this report, we present a resolution algorithm which includes most of these features: (i), (iii), and (v) through a finite-element approach, (ii) and (iv) through adequate use of classical methods. The results of numerical experiments are presented to show the possibilities of the methods discussed in this report. 1. FORMULATIONOF THE NAVIER-STOKES EQUATIONS FOR HIGH REYNOLDSFLOWS It is well known [l] that with few initial data, the Navier-Stokes equations possess a unique solution. But in the usual case of many stiff initial data, instabilities appear and degenerate in turbulence. We make the usual assumption that these chaotic flows are within the range of validity of the NavierStokes equations (albeit not within the capabilities of our numerical approximation of these equations). Let D E RN (N = 2 or 3 in practice) be the flow domain and F be its boundary. The conservative form of the equations is given by
$+v.(pu)=o ~+v.(puc3u)=v.(7),
(1.1) 7ij
=
7ji
(1.2)
66
S. BOIVIN
(1.3) 711=
-Ph, + adub4,
(1.4)
+ &S,,(u)
P = (Y - l)w,
(1.5)
where S,j(u) = $( u,,, + Uj,i), E = e + $u2, e = C,T, y = C,/C,, q = -kVT, x = -&L. Despite all the recent advances in computer technology, turbulent flows cannot be computed at present by the direct resolution of the Navier-Stokes equations. The reason is that the turbulent motion contains a continuous range of scales in which the ratio of the larger to the smaller scale is about 0(Re3’4). For this reason, and because one is generally not interested in complete details of the behavior of all the instantaneous variables, we work with the mass-weighted averages of those quantities. To clarify, we recall the definition of the averaging operators. If p is the instantaneous density and u the instantaneous velocity, then we define the average of p and the mass-average of u on a time interval of length to by
p”
s 10
to
0
pdt,
a=& Pto
s 10
0
pudt
and u’by u’ = u - 6. To simplify the writing, we will retain the same notation for the mean variables and use the prime for the fluctuating part of mass-weighted variables. Consequently, for the remainder of the paper p, p are mean density and pressure; u, e, Tare mass-weighted mean velocity, internal energy, and temperature; u’, T’ are the fluctuating part of the velocity and the temperature; and u”is the mass-weighted averaging operator. Now, neglecting buoyancy effects and all viscosity and thermal conductivity fluctuations, the open equations for mass, momentum, and total energy conservations are
$+v.(pu)=o
(1.6)
~+v.(,,@24)=..((7_(,uGL~)) y
+ V.(pEu)
= V.[(T
- (pu’&‘)).
u - qtol],
(1.7) (1.8)
where pE = pe + $pu2 + ~PU’~ = pe + jpu* + pk, k being the turbulent kinetic energy, and qtOlthe total (laminar + turbulent) heat transfer.
67
COMPRESSIBLE NAVIER-STOKES EQUATIONS
The system ( 1.6) - ( 1.8 ) contains only second-order moments as supplementary unknowns, higher order moments being neglected and correlation between velocity and temperature =
modeled by
qi,tot =
-k(aT/dXi)
+
G
-k&(dT/dXi).
There are two ways to express the Reynolds stress tensor R = -pu'% u': To assume that the Reynolds stresses are a local property of the mean flow and, by analogy with the viscous stresses, express them by a constitutive law. ??
To assume that the Reynolds stresses are independent variable quantities and obtain their values through the resolution of a set of transport equations. ??
In this report, we consider only the first approach, the second requiring too much empiricism [ 21. The choice of a constitutive law relating R to mean quantities is very difficult. But, under certain approximations, Chacon and Pironneau [ 31 make it clear that, in the incompressible case, R' and (VU + Vu’) generate the same subspace of second-order tensors; hence in 2D R’ = cd+
B(Vu + Vu'),
(1.9)
where (Y,/3 are only functions of the nontrivial invariants of (V u + Vu’), and possibly of the turbulent kinetic energy and mixing length. We suppose that this can be generalized to the compressible case as
R = --put&d = -$(pk
+ /.L~V.u)l+
/&VU + Vu').
(1.10)
More general hypotheses [4] seem to give more realistic approximations, but without any proof that the Reynolds stress tensor still lies in the same space of tensors. In a classical way, we completely decouple the resolution of the mean Navier-Stokes equations ( 1.6 )-( 1.8 ) and the calculation of the Reynolds stresses. The next section is devoted to the resolution of the mean Navier-Stokes equations and Section 3 to the modeling and resolution of the model giving the Reynolds stresses. 2. RESOLUTIONOFTHENAVIER-STOKESEQUATIONS FORTHEMEANFLOW
We now consider the following nonconservative the Navier-Stokes equations ( 1.6 ) - ( 1.8 ) ,
nondimensional
form of
68
S. BOIVIN
ap
dr+u.vp+pv*u=o
(2.1)
$+(u.v)u+h9-~v.(“*o)=0
(2.2)
P
~+U.vTt~v.U-Vu: P
( 1-h.(K*vT)=o. V*
-0.
P
P
(2.3)
In (2.1)-( 2.3 ), we have normalized by the subscript r (i) the density p by pr (ii) the velocity u by I ur I (iii) the internal energy e by ( ur 12 (iv) the pressurep by p,(u,l’ (v) the viscosities p and pt by pr
(vi) the kinetic turbulence energy k by I ur 12 (vii) the temperature
T by I ur 1*/C,
which imply e = T. The pressure obeys the ideal gas law; the number y and the functions u, 8, v*, K* are defined by: ??
a=vu+vut-
#V.uI,
??
0 = p + $pk (p = (y - l)pT),
??
y = C,/ C, is the ratio of the specific heats ( y E 1.4 in air),
v* = (p + pt)/Rer is the total viscosity defined from the computed laminar and turbulent viscosities divided by the reference Reynolds number Re, = P~w%IcL~, K* = y/Re,( p/Pr + pt/Prt) is the total conductivity coefficient, also defined from laminar and turbulent viscosities. ??
??
Remark. The field (v* )-’ will be referenced as the total Reynolds num-
ber field. We consider external flows around 2D geometries; the domain of computation is described in Fig. 1a. Let r, be a far-field boundary of the domain; we introduce
r, = {XIXEr,,U,vz
r; = r,.w;,
where u, denotes the free stream velocity and n the unit vector of the outward normal to r.
COMPRESSIBLE NAVIER-STOKES
a
b
c
d
EQUATIONS
69
“t U “2 “3
FIG. 1. (a) Computational domain. (b) P, -iso Pz element: @ velocity, density, temperature, q and w d.o.f.; 0 velocity, q and w d.o.f. (c) Upwind finite element T(s), to node s. (d) Upwind finite elements T(s,), to nodes s, of a given element.
We assume the flow to be uniform at infinity, and the corresponding variables to be normalized by the free stream values; then, for example, we prescribe at infinity (Yis the angle of attack, P=
1,
T= Tm.= 1lI-d~ - lMb1, where it4, denotes the free stream Mach number. The boundary conditions on the computational boundary I?, are on FL :
u=u,,T=T,,p=
1,
and on r”,:
VT.n=O,[n.(Vu+vut-fv.Uz).n]=O,
[n*(Vu
+ Vu’ -
$0. uZ)* t] = 0,
where n, I are the local normal and tangent on I’:.
S. BOIVIK
70
On the rigid boundary rB, we shall use the conditions ( no-slip condition),
u=o T = TB = T,(l
+ (y - 1)/2M:,)
(free stream total temperature).
Finally, since steady solutions are sought through time-dependent equations, initial conditions must be added; we shall take P(X, 0) = PO(X),
T(x,O)= To(x).
u(x, 0) = e(x),
Solving the compressible Navier-Stokes equations is a difficult task. Most of the existing numerical solution methods are based on finite-difference techniques, for both space and time discretizations. Following the work of Bristeau et al. [ 5, 61, we will consider a method based on finite-element techniques for space discretization while using finite differences in time. Let the time derivatives be discretized using a classical implicit Euler finitedifference formula; then at each time step, we solve the nonlinear system of variational equations ol(p - ;, N) + (U-VP, N) + (~0. U, N) = 0
(2.4)
LY(U-G,M)+((U’V)U,M)+
+(FVT,VK)=O,
(2.6)
where the solution is sought in V X W X Z, Y=
EH’(n)lp,r,
w=
{uE(H’(Q))+Q.~
Z=
{TEH’(O)lT,r,
although there is no triple (p, u, T) E V of test functions (N, of the spaces X(Q), Remark
{p
= 1) = U,, UjrB = O} = T,,
TI,., = TB),
complete theoretical justification for this. We look for a X W X Z such that (2.4)-( 2.6) is verified for all triples M, K) E Y(Q) X (X(Q))* X x(Q). For the definition Y(Q) see Appendix: Nomenclature.
1. ( , ) denotes the scalar product in L2( Q).
COMPRESSIBLE
NAVIER-STOKES
EQUATIONS
71
Remark 2. The natural boundary conditions already defined were introduced by setting the boundary integrals appearing from the integration by parts to zero. Remark 3. The 0 term of Eq. (2.5 ) is integrated in the form
1
l)vT+$vq+ (y- 1)f+2dVP.
-+(r-
3P
[
We do not integrate by parts the pressure term of the momentum equation because we cannot set p + 6,, = 0 on the boundary (doing so would perturbate the solution) and we prefer to avoid the computation of boundary integrals. To approximate (2.4)-( 2.6), by the finite-element method, we must divide 52into small elements (triangles) and replace all the functions by their interpolate oh, uh, Th. Interpolates are defined inside the elements from their values at the nodes of the elements by an interpolation formula. The choice of the approximation spaces is a critical and still partly obscure point. But, we see that when u is very small, for example, near a solid boundary, Eq. (2.4) looks like a divergence free equation. This is especially true after the change of variable: (I = log p (see [ 5 ] ). We may conjecture that to ensure stability, a necessary condition is that the element chosen should be stable for the incompressible case. This seems to be confirmed by the various experiments made by Bristeau et al. [ 61. A possible choice for such an element (and consequently spaces) is the PI -iso P2 element (Fig. 1b). If Th12is the triangulation obtained from Th by dividing each triangle of Th into four equal triangles whose vertices are the vertices { qr, q2, q3} of Th and ( h(ql + q’), f(q2 + q3), f(q’ + q’)}, we define
v, = {~/a E VIP~ E c’(Q), VTE
Th, PAITE f”(T)} C
w, = {uh E WjUh E (C0(fi))2, VTE T/,,z, Uh,TE P’(T)}
v c w
Zh= {ThEZIThEC”(n),VTETh,Th,TEP’(T)}CZ, where P’ ( T) = { set of polynomials of degree d 1 on T} . Restricted to these finite-dimensional spaces, Eqs. (2.4)-( 2.6) lead to the following nonlinear problem: Find (ph, uh, Th) in (VAx wh x z,,) a
SOhtiOIl
of Fh(P,,, u,,, Th) = 0,
Fh being the discrete version of the system (2.4)-( 2.6).
(2.7)
s. ROlVllU
72
We now consider iteration schemes for solving the nonlinear system F,?(s ) = 0, where s = (oh, ul,, T,,). Newton’s method applied to this system results in the iteration 1. Set so an initial guess 2. For IZ= 0, 1, 2, . . . , until convergence do 2.1. Solve J(s”)6” = --F,,(Y), (2.8)
2.2. Set sn+’ = sn + 6”,
where J( s”) = F’(f) is the system Jacobian. For large problems, iterative methods are frequently used to solve (2.8) only approximately, giving rise to methods which can be viewed as inexact Newton methods. The particular method we use is the Generalized Minimum Residual Method (GMRES) due to Saad and Schultz [ 71. This method has the virtue of requiring virtually no matrix storage and only the action of the Jacobian matrix J times a vector r, and not J explicitly. In the setting of nonlinear equations, this action is approximated by a difference quotient of the form J(s)r
z
F/,(s+ br) - F/,(s) h
3
where s is the current approximation of a root of (2.7) and b is a scalar. For details of the algorithm, see [ 71 and also [ 61 for a clear setting of this algorithm within the context of the conjugate gradient methods. To conclude this section we present and discuss two strategies to add some numerical diffusion to the approximation of the convective part of the equations. It is well known that when the Reynolds number is high, boundary layers and shocks become hard to resolve and centered approximations of first derivatives are dangerous. In fact, it is necessary to add some numerical diffusion where the mesh size h is such that the local Peclet number satisfies the condition Pe=@>2, C
a = P max(Iuil,
lull),
to damp unresolved frequencies. There are many different ways to do that, upwinding, addition of an artificial viscosity, etc., but the problem is to control the numerical diffusion to be sure it does not degenerate the approximation of gradients. Our approach is to filter the discrete gradients in the following way: as the discrete gradients are constant by element, we define the gradients at a node
COMPRESSIBLE
NAVIER-STOKES
EQUATIONS
73
or on an element by an upwinded combination of the gradients on the nearby elements. The first method can be viewed as a projection followed by a combination with the upwinded gradients, that is, if a node s is surrounded by elements 1 to n , with respective area ak which sum to A. Then we define the gradient of a function 4 at s by the combination of a centered gradient and an upwinded gradient
VA = (1 - /w4: + PV4F,
(2.9)
where the centered and upwinded gradients at s are respectively defined by
where T(S) is the upwind element at node S. The choice of p is discussed below. The second method is a pure filtering and can also be written as a combination of a centered and an upwinded gradient. More precisely, on an element e we define the gradient by
vdJt== (1 - P)V4c,+ PVK,
(2.10)
where the centered and upwinded gradients on e are respectively defined by
A fundamental difference between these two upwinding methods (for elements where gradients are constant by element) is that in the second there is no need for projection. We now describe the computation of p. The method, based on the analysis of the scalar steady convection diffusion equation, is to relate the coefficient to a local Peclet number by 0 = coth F - ;, ( 1
where the Peclet number is an estimation of the importance of convection compared to diffusion. Although the definition of this number is obvious in the scalar case, the generalization to multidimensional cases is not immediate. But in general we may define it by
74
s. BOIVIN
Pe = Re phi, where h is an estimation of the characteristic length and a”an estimation of the convection speed. A possible choice for A, a^is it =
I(~~~)lll~l,
a= IUI.
(2.11)
But some computational experiments demonstrated that the following choices, inspired from the work of Hughes et al. [ 81, give better results, Ji = I(~~~~)lll~~l,
6 = l(~~~~)Ill~~l.
(2.12)
In fact, with these definitions, /3 is more closely related to discontinuities. Note also that the coefficient ,f3can be multiplied by a tuning coefficient for further control of the diffision and more specifically to weight the diffusion added in each equation. The equation for p being the only equation without any physical diffusion it is fair to believe that more numerical diffusion should be added to this equation. Although less robust, the second method, the filtering method (2.9)-( 2.12), proved to be very efficient and precise, provided that the mesh was adequately refined. Remark. The value of h = (h,, h,) at a given point of the mesh is deduced from the computed value on each element using the relations
where qi are the PI basis functions on a given triangle.
3. CLOSUREMODELOFTURBULENTFLOW The basic goal in developing turbulence models is to specify closure conditions in which the unknown Reynolds stresses of turbulence are related, either algebraically or through differential relations, to known mean-flow variables. The fundamental problem is to deduce a suitable mathematical model which should be a good compromise between mathematical simplicity and the complex exact transport equations of Reynolds stresses. Our starting point will be the classical k - c model of Jones and Launder [ 91. A reason for it is that this model can be deduced rigorously in the in-
COMPRESSIBLE
NAVIER-STOKES
EQUATIONS
75
compressible case [ 31 and also because there are some efficient (but not trivial) ways to solve it. It is well known [lo] that the choice of the scheme to solve the k - c equations is very important. Thus, to ensure numerical stability we make a change of variable inspired by Coakley [ 111. The basic closure turbulence model as described under Appendix: TwoEquation Closure Model is given by
am) at
+u.v(Pq)+;pqv.u--IR:vu+~pwq 37
-v.[(,+:W]=o
ah4 +
at
u*V(pw)
+ pwv-
u - (Cl - l$R:
(3.1)
vu + (c2 - l)pw2 -v.[(p+~~~]=O,
(3.2)
where q2 = k, w = e/k,pt= c,,pq2/w, and R is the Reynolds stress tensor. In the construction of the two-equation closure model it is assumed that the flow is fully turbulent and that the Reynolds number is everywhere high. Thus, to provide predictions of the flow within the viscous layer adjacent to the wall, the above model must be extended in three ways. These are: (i) prescribe the turbulent / nonturbulent
interface,
(ii) the terms containing the cl and c, will become dependent upon the Reynolds number of turbulence, (iii) if possible, correct values using empirical formulas. The complex problem of the interface between a turbulent region and a nonturbulent region will not be discussed here. The corrected values for the constants, proposed by Coakley [ 111, are given by replacing c, by c,D and (1 - cl ) by (0.405 D + 0.045), where D is a damping factor given by
D = [l - exp-“‘I,
r=
p42 110 ’
where the value of the constant is a! = 0.0018. Initial profiles for q and w and corrected values are deduced from an estimation of the turbulent viscosity provided by an algebraic model, namely the Baldwin-Lomas [ 12 ] model, through the empirical relations
S.
76
BOlVllv
(Bradshaw hypothesis) 3 w = -s,* 10
( Hypothesis: production
= dissipation).
As in the case of the mean Navier-Stokes equations we work with the nonconservative and nondimensional form of these equations. That is, ~+~.V~-~~V.u-~c,D(a:Vu)+~wg-~V.[C,Vql=O 2
(3.3)
+ u*vw + ; (c, - 1)oV. u - (c, - I)c,D(a:
Vu)
+ (c* - l)w2 - A V*[C,Vw] P
= 0.
(3.4)
In (3.3)-( 3.4), we have normalized by the subscript r (i) p, u, T as in Section 2 (ii) the viscosities p and pt by pr (iii) the kinetic turbulence energy q* by I u, I ’ (iv) the specific dissipation rate of q*, w by 1u, I/L,. The functions D, 6, C,, C, are defined by: ??
D = 1 - exp(-apq*Re,/w),
??
a=Vu+Vu’-
??
C, = (CL+ d~y)/Rer,
??
C, = (P + dcJ/Re,,
~04,
The values of the constants c,, cl, c2, ug = uk are taken from the original paper of Launder and Spalding [ 13 J; for the value of 6, see the first appendix. The boundary conditions on smooth wall boundaries are q = 0 and VW. n = 0. At the computational boundary I’, small values are prescribed: q = qa , w = w, . More exactly, qa is set to the known free stream value of the kinetic turbulent energy and w, is defined by setting Put=
Re,c,p,&
= I
W,
We now describe the resolution strategy. As already noted, the overall resolution is split in two parts; first we solve the mean Navier-Stokes equations and then we solve the turbulence
COMPRESSIBLE
NAVIER-STOKES
EQUATIONS
77
closure equations. At each step the resolution of (3.3)-(3.4) is done with p = i, u = u^,P = ji, that is, with known values from the preceding NavierStokes step. The technique used to solve (3.3)-(3.4) is very similar to that used, in Section 2, to solve the mean Navier-Stokes equations and is based on the following weak formulation. Let the time derivative be discretized using a classical implicit Euler finitedifference formula; then at each step, we solve the nonlinear system of variational equations
+;(wq,N)+($‘q,VN)=O
a(w-i,K)+(ri*vqK)++ - (c, - l)(c,Di
(3.5)
l)(wV*ti,K)
: Vii, K) + (cz - l)(w $K)+($
VW, VK ) = 0,
(3.6)
where the solution is sought in Q X S,
Q = {QE H’W)lq,r, = q-3
s = {w E H’(O)lo,r,
qlrB =
0)
= cd,}.
Here again, we obtain an approximation of the solution using the inexact Newton-GMRES method on the nonlinear system (3.5)-( 3.6) considered on the finite-dimensional approximation spaces
Qh = {Q/, E Qlsh E cot% &=
b’7-E l-h/z, 4h,rE P’(T)}
C Q
{OhESIOhEC”(,),VT.Th,2,Wh(TEP’(T)}CS,
where P’(T) = {set of polynomial of degree <1 on T} and Th,2 is the fine grid of the P, -iso P2 grids system discussed in Section 2. To stabilize our scheme, we replace the convection terms (u^*Vq, N) and (11. VW, K) by filtered terms (see Section 2 for the filtering technique). The second model of turbulence we use is an algebraic model derived from the Baldwin-Lomax [ 121 model. We include this model as a basis for comparison and eventually to decrease the cost of computation. (Another reason was to try to combine both models but preliminary computations led us to the conclusion that the two-equation model gives better results.)
S.
78
BOIVIN
In this model, the turbulent viscosity is computed using a two-layer approach,
where Y is the normal distance from the surface and Y, is the least value of Y at which pt, = ptO. In the inner layer pti =
K’plw( Y2[1 - exp(-Y’/26)12,
where w is the vorticity, and Y + = Y ( p0 I c, 1/I;) ‘I2 with the subscript o referring to the value of the parameter at the nearest solid boundary, and mu is the local shear stress in the direction of the flow. In the outer layer
whereF,,,,,=max(Ylwl[l -exp(-Y+/26)]),andY,,,,,isthevalueofYat which F,,,,, occurs. The KIebanoff intermittency correction is given by Fkleb = [1 + 5,
~(~klebY/Ymd61-‘.
The model constants and functions are K = 0.40, k = 0.0168, 0.2933M,
C,, =
C kleb
I
ifA4, < 3; ifM,
2.08,
-0.4118A4, 0.3,
=
+ 1.2,
+ 0.65,
I
> 3,
if M, < 0.85; if M,. > 0.85.
In the wake region, the following nonclassical formulation is used,
Pt = CLtr(IWI/I%l), with the subscript r referring to a reference profile, usually the last profile before the trailing edge. This unusual wake formulation enables us to simulate the characteristic transport-diffusion behavior of the turbulence in the wake, keeping the emphasis on the evolution of the vorticity. This scheme avoids the appearance of steep gradients when moving from the wall model to the wake model of the usual algebraic models. Some details about the implementation of this model for unstructured meshes are given below.
COMPRESSIBLE
NAVIER-STOKES
EQUATIONS
79
This kind of model is one of the most widely used because of the ease of its implementation in finite-difference codes. The basic property appearing here is the orthogonality between the mesh lines and the walls. This important property will have to be replaced by an abstract structure for general finiteelement meshes. The most obvious way is to compute all the model parameters on a second mesh ensuring the orthogonality property and to use an interpolation operator between the two meshes. This procedure can be difficult to implement for complex (2D) geometries and inadequate in 3D. Another way to implement algebraic turbulence models for general meshes is to consider independently each node. That is, keeping some simple geometric characteristics of each node in a table, for example, the nearest point of each wall, the normal lines passing throughout this node can be found and so all the model parameters can be evaluated. But as usual, the simplicity of this method must be paid for; the maximization of functionals along the normal lines can be very costly. Consequently, the complexity of the process must be diminished by the use of approximate maximization techniques and eventually by restricting the process to adequately chosen nodes. These simplifications make the overall process interesting. 4. RESOLUTION
SCHEME
We summarize here the overall resolution scheme as described in the previous sections. This reads: Step 1: Read initial values of u”, p”, To or let p” = 1 and solve A.u” = 0, AT0 = 0 with appropriate boundary conditions. Step2:
Let q" = qcu,oo = w,(qpr, = 0).
Step 3: u”, p”, T”, pf, q",W" are known from the previous time step. Step 3.1: Update the laminar viscosity (P”+‘) using the Sutherland formula. Step 3.2: Compute the algebraic viscosity (&zi) using the BaldwinLomax model. Step 3.3: Make a correction to q" and J’, 4 “+*‘* = max[q”, qa],
cd"+"*= max[w”, w,]
or 4
n+1/*= max[(l(gzl$2)~'*,q.,
w n+t/*= max
+2,
[
an,
WCC
1 -
q_]
s. BOIVIN
80
Step 3.4: Solve the q - w equations F+( q" ’, w"+' ) = 0 (see Section 3). Step 3.5: Update the turbulent viscosity (this step may include a projection of q and w to eliminate negative values),
or p;+’
=
max[c,D(q”+‘)‘/o”+‘,
&,:I.
Step 3.6: Solve the N-S equations FNs( z~“+‘, p”+‘, Tnf’) = 0 (see Section 2 ) . Step 4: If Is”+ - s” ]L~cr2)< c, sn = (p”, un, T”), then stop, else go to step 3. 5. NUMERICALEXPERIMENTS Before the exposition of the results of the numerical experiments let us recall that as v* is defined from the laminar and the turbulent viscosities, it is not a fluid property. This field is of primary interest; more exactly we are interested in the field 1/v* called the total Reynolds number. Computations were performed for some flow condition around an NACAOO12 airfoil and a simplified version of the Hermes space shuttle. In each case, the steady-state solution was reached after 80 to 120 time steps depending on whether the initial flow field was uniform or a known flow field from a previous computation. Two test cases were chosen here for the computations, each case involving complex flow with shock waves. The first series of numerical experiments concerns a compressible viscous and turbulent flow around an NACAO012 airfoil. A part of the coarse mesh used is shown in Fig. 2a. We took here M, = 0.85, Re = 105, y = 1.4, Pr = 0.72, Pr, = 0.9, and a zero angle of attack. We have shown in Fig. 3 the iso-Mach contours (Fig. 3a) for a NavierStokes computation with the total Reynolds number field deduced from the Sutherland formula (Fig. 3b). The same problem was solved a second time but using the two-equation turbulence model, discussed in Section 3, for the computation of the total Reynolds number field. See Fig. 4 for the iso-Mach contours and the total Reynolds number field and Fig. 5 for the iso-k and iso-w contours. Despite the problems caused by the coarse mesh in the wake, we note the presence of vortices in both computations. In the turbulent one, the boundary layer thickened abruptly but did not damp out the time-dependent nature of the shocks and the wake vortices.
COMPRESSIBLE NAVIER-STOKES
mesh: 11~30012
- 2000
mesh: pseudo-hermes
-
nodes
1883
nodes
81
EQUATIONS
Univ. Lava1
Univ. Lava1
27 Jul 88
27 Jd 88
FIG. 2. Partial view of the principal meshes.
It is clear that this flow is dependent on the capabilities of the mesh we use. But it is interesting that most of the physical phenomena we tend to see are present. An obvious way to obtain better results is to refine the mesh in regions where steep gradients are present. This is only a matter of computer facilities, but we are now developing a strategy to obtain better results via the use of enriched finite elements in the corresponding areas. The second series of numerical experiments concerns a compressible viscous
S. BOIVIN
iso-Mach
MAX=O.l78e+Ol
MlK0.485~~19
univ. Lawtl 27 Jul88
,
FIG. 3. Mach contours and total Re number; laminar computation at A4 = 0.85, Re = IO>.
and turbulent flow around a simplified version of the Hermes space shuttle. A part of the coarse mesh used is shown in Fig. 2b. We took here M, = 2.0, Re = 104, y = 1.4, Pr = 0.72, Pr, = 0.9, and a zero angle of attack. Note that the length of the body is approximatively 0.14 of characteristic length. We have shown in Fig. 6 the iso-Mach contours (Fig. 6a) for a NavierStokes computation with the total Reynolds number field given by the Sutherland formula (Fig. 6b). The same problem was solved a second time using
COMPRESSIBLE
iso-Mach
NAVIER-STOKES
MAX=O.l67eGIl
83
EQUATIONS
MIK0.485e19
Univ. Lad 26 SW
06
\ Re_total
MAX=0.1221?+06
MIN=O.l95e*03
Univ. Lava1 26Sep66
nc.
4. Mach contours
and total Re number:
q - w computation
at M = 0.85,Re = 105.
the two-equation turbulence model (see Fig. 7 for the iso-k and iso-w contours), but all the characteristics of the flow remained similar (see, for example, a comparison of the friction coefficients for the turbulent and nonturbulent computations shown in Fig. 9b). This is encouraging because at this low Reynolds number the flow should be mainly laminar. We have conducted a sensibility experiment on the values of the constant a, of the w equation. The new iso-k and iso-w contours for the last q - w computation with a, = 13
S. BOIVIN
84 /
I
I
I
k
MAX=0.625e-01
w
MAX=O.l08eO3
MIN=O.OOOeOO
MIK0.9OOe-02
univ. 27
Lava1 Jul68
univ. Lava1 27
Jul98
FIG. 5. k and w from a q - w computation at A4 = 0.85, Re = 105.
instead of 1.3 are shown in Fig. 8. Although satisfactory for a very wide range of flow conditions [ 141 preliminary results on simple nearly incompressible cases tend to indicate that this choice for 0, leads to overdiffised values for k and q* and w. That is, the constant u, is not universal. Finally, a numerical experiment was conducted to compare the diffusion associated with the turbulence model used. The test case was the transonic flow around an NACAO012 airfoil at AI, = 0.85, Re = 104, y = 1.4, Pr
COMPRESSIBLE NAVIER-STOKES
iso-Mach
MAX=0226e+Ol
85
EQUATIONS
MIrG=O.441e-19
Univ. Laval 27 Jtd 66
////////I
I-_ .~a
MAX=O.l lOed
MltWS23e94
Univ. Lava1 27 Jul66
FIG. 6. Mach contours and Re number; q - w computation at M = 2., Re = 104.
= 0.72, Pr, = 0.9, and a zero angle of attack. Only the resulting values of the coefficient of friction (see Fig. 9a) are discussed. The Navier-Stokes computations (inner curves) show a large separation zone while the use of the algebraic turbulence model gives no separation at all (outer curves). It is clear that the diffusion caused by this turbulence model is overestimated. The estimation given by the two-equation model (see middle curve in Fig. 9a) is more realistic.
S. BOIVlN
86
! /
k
/
w
MAX=O.445~02
MIN=aQOOe03
univ. Lava1 27 Jo188
FIG. 7. k and w from a q - w computation at A4 = 2.0, Re = IO?
The choices made in steps 3.3 and 3.6 of the resolution scheme (see Section 4) were proved to be less important than expected. In fact, the solution process described to solve the turbulence transport equations with the given boundary and flow condition converge without the need of the algebraic turbulence model. The latter may be used to start the q - w computation from a given state, causing the steady state to be reached more quickly. We present here some statistics to give one an idea of the time used for these computations on a SUN3 workstation with a floating point accelerator
COMPRESSIBLE NAVIER-STOKES
87
EQUATIONS
MAX=0.621e-O2
I
w
MAX=O.566&2
MIKa900e03
W.
Lava1
30 JUl88
FIG. 8. Same as Fig. 7 but with a different constant of diffusion in the w equation.
(0.4 Mflop). For a computation on the 1883-node simplified Hem& shuttle, it takes at each time step:
space
110 s to compute the turbulent viscosity with the algebraic model, 2300 s to compute the turbulent viscosity with the two-equation model (complete resolution at each time step), 1100 to 2000 s to solve the Navier-Stokes equations.
S. BOIVIl\i
88
0
02
0
01
0 00
-0
01
-0.
02
-0
03 1 lII/,IIII
0 00 b
0. 04 0 03
0 02
0 01
0 00
-0
01
-0
02
-0
03
i
i
I,,,
0 25
0. 50
I,,,
0. 75
1 00
-.--. _=-_ _-_ ---_ _ _ .I:.... -I..._..I _..........~ _.
‘.;=
__---__-_--
_---
_--_
--
---
-
-0.04
0 00
0.03
0 07
0. 10
0. 1
FIG. 9. Comparison of friction coefficients;see text for details.
More numerical experiments and further details about the method can be found in a research report [ 151 published at INRIA, France.
6. CONCLUDING REMARKS We have presented a numerical scheme to solve the compressible NavierStokes equations written in nonconservative form with emphasis on turbulence modeling. We have obtained good results for moderate Reynolds numbers using an approximation satisfying some “inf-sup” condition and using a filtering of the gradients of the convective terms. This filtering is dependent on
COMPRESSIBLE
NAVIER-STOKES
EQUATIONS
89
the local direction of the fluid flow and on the local amount of physical viscosity. A particular form of the k - c model was solved with standard boundary conditions to obtain an estimation of the turbulent viscosity using the same techniques. Further developments are needed to obtain accurate solutions for higher Mach and Reynolds numbers where special care must be taken to not introduce too much numerical dissipation. In the near future we plan to compare computations on finer grids with experimental tests and emphasize filteringupwinding methods and also the resolution of the boundary layers. APPENDIX
Two-Equation Closure Model We present here some details about the two-equation closure model ( 3.1)(3.2) that we used. Let us recall the simplest and certainly the most popular closure model, namely, the k - Emodel, where k is the turbulence energy and e is its rate of dissipation. According to this prescription
a(Pw
-+uV(pk)=-pkV.u+R:Vu+V.
at
am at
+uV(ps)=-p67~u+c,fR:Vu+V
[(p+$Vk]-pe .[(P+37c]-c2P~?
(1) (2)
where c,, (lk, ur, cl, c2 are model constants, usually set to the recommended standard values 0.09, 1.O, 1.3, 1.44, and 1.92, respectively. It is well known that near a solid boundary (at the boundary we suppose k = t = 0) this model needs a low-Reynolds-number dissipation term to balance molecular diffusion in the sublayer. In fact, letting y tend to 0, Eq. ( 1) becomes
which is incompatible with the experimental fact that k - y2 in the sublayer. Instead of adding empirical terms to correct it, we assume k to be constant in the diffusion term and make the following change of variable k = q2, from which we deduce directly Eq. ( 3.1). Note that if we assume q2 = k - y2 we obtain the correct behavior of q when y + 0. Finally, to reduce the problems of numerical instabilities, we reduce the truncation errors due to quotients, using the change of variable w = c/k. But this must be paid by supplementary hypotheses.
S. BOIVIN
90
If we insert w = c/k in Eq. (2) and use the identity k(pw)’ = up’+ (PC)’- u(pk)’ we find the equation
ah.4 -=w~+wu.~p-~.V(pw)+(~~-l)wR:Du
at
42
+ D - (cz - l)pw2,
(3)
where the diffusion term reads D=V.[(p+;-Ve]-wV-[(p++‘k].
(4)
Thus we need to assume (pt + ~*/u,) = (p -t p,/c~k) = cte in order to write it in the simpler form
where the recommended value for Q, [ 1 l] is 1.3. But, although satisfactory for a very wide range of flow conditions [ 141 preliminary results on simple nearly incompressible cases tend to indicate that this choice for u, leads to overdiffused values for k = q* and w. The standard form for this equation follows from the use of the nonconservative mass balance equation.
Nomenclature Letters incompressible Reynolds stress tensor compressible Reynolds stress tensor viscous and pressure stress tensor viscous stress tensor strain tensor laminar viscosity turbulent viscosity reference viscosity turbulent kinetic energy
'smog
snoas!A a[q!ssaJduroau! pue aIq!ssaJdwoa JO uoyelnmp aql 01 uoyagddy
S=fOlS-Ja!AEN ‘( L861)
gLI
aql 'OJ SpOqKK4I .zpqq
p!nld T
ll23!JZNUIlN ‘XllE~?d
'f PUFJ ‘!~SU’MOI~
‘a
.suogenba
‘nl.?a~S&I~ ‘0
‘m
‘s
.aaualnqJnijo stapolu 3 - y pua I-? Jeauquou u0 ‘afe!zadS ‘S ‘3 ‘p
‘(9861) YJoA MaN/u!lJatl ‘%‘aA “XII aJt?M$jOs uoyezyqdo ‘S3!lVlUJylV~ pqddy u! SVLVA ‘(‘pa) ueuqs&~+?[s~ ‘A ‘V u[ .lapom 3 - 9 aqlJ0 uoyJpunoJ Iea!i~maqleur aql u0 ‘nsauuoJ!d ‘0 pur! uoaeq3 ‘1 ‘E
-Ja%U!JdS
3puepaqiaN aqL ‘qaJeas -aa
‘OJpilH
‘OSSV
‘WUI
'SJ!lllVlpkH U!
uo,m,qddy
.l!JlfLpUV
Slapow 't2861)
-qUoN ‘I\
'~U!.IJJU@U~
uo!lour JO suogsnba
aql
pun
SJ~UJ!~~
JOJ SLW~qOJd
pu= ‘IOPJQ ‘8 ‘nEaUUOJ!d aM ‘uo!l’+V la+J@lln==a
‘0
pqddy
uf SPO~IJM
an@AiClEpUnOq
‘!‘lSU!MO[f) [a3JEm
‘u
SUO!AV
‘U!lJOJ “OJJ
Sqndtuo3
[qpq
‘m
‘Ep!qS!N
J3UJlTlq.Vl~ ‘!pOx
'M
'2
uJePJal=W ‘PuelloH y ‘spq ‘1
pue
lEJaua;l JO
eJeuInsle~
‘v
1
3uo!ssnas!p pg!nJJ JOJ xnepad ‘1 ‘0 ‘m xueq$ LIIe!aadsa pauoddns L[&wds! ~JOM s!qL
OSlC PUP nW$S!Ja
KWJ1UO3 l? bq
*sa~gt?“~~appue saxpu! pawadaJ “03 uo~y.~a~uo~pa)y~ -0sse aqt put? uogwou losuat ~txymp ayl asn ahi ‘papaau uayM ym4.4a~
sa+w+ap
lsry
‘10 = EL”XI(0),H 3x1 = (0)x (0 = =“‘xl(D),H 3x1 = (0)x (z&7 u;[ ov a= (7~)~7 30 suogx1n3 30 axds Aaloqos pmsn (0 ) , H (8 )I 7 uoyun3 alqt?Salur an%saq~ aanbs 30 axds
asoym
16
SNOI.LVL%XI SEfXO_LS-XIIAVN 318ISS3lIdI403
92
S. BOIVIN
6. M. 0. Bristeau, R. Glowinski, B. Mantel, J. Periaux, and G. Rage, Adaptative finite element methods for three dimensional compressible viscous flow simulation in aerospace engineering. In 11th Inter. Conf: on Numer. Meth. in Fluid D_vnamics. Springer-Verlag, Berlin/New York. in press. 7. Y. Saad and M. H. Schultz, GMRES: A generalized minimal residual algorithm nonsymmetric linear systems. SIAM J. Sci. Statist. Comput. 7 (1986).
for solving
8. T. J. R. Hughes, M. Mallet, and A. Mizukami, A new finite element formulation for computational fluid dynamics. II. Beyond SUPG. Compul. Methods App/. Mech. Engrg. (1985). 9. W. P. Jones and B. E. Launder, The prediction of laminarization of turbulence. Int. J. Heat Mass Tram@ 15 (1972).
with a two equations
model
10. A. G. Hutton, R. M. Smith, and S. Hickmott, The computation of turbulent flows of industrial complexity by the finite element method-Progress and prospects. In Finite Elemenfs in Fluids, Vol. 7. Wiley, New York (1987). 11. T. J. Coakley, Turbulence modeling AIAA Paper 83-1693 (1983).
methods
for the compressible
12. B. S. Baldwin and H. Lomax, Thin layer approximation turbulent flows. AIAA Paper 78-0257 (1978). 13. B. E. Launder and D. B. Spalding, The numerical Methods Appl. Mech. Engrg. 3 ( 1974). 14. T. J. Coakley, (1987).
Numerical
simulation
and algebraic
computation
of viscous transsonic
Navier-Stokes
model for separated
of turbulent
Navier-Stokes
16. B. Escande, Modelisation de I’interaction onde de choc/couche limite turbulente. calcul-experience. These Doctorat, Ecole Centrale de Lyon (1986). Separated
flow treatment
with a new turbulence
18. M. Tabata, Finite element approximation corresponding Mem. Numer. Math. Univ. Kyoto Tokyo Y (1977 ).
flows. Comput.
airfoil flows. AIAA Paper 87-046
15. S. Boivin, A finite element method to solve the compressible turbulence modelling. INRIA Research Report (1988).
17. U. C. Goldberg, (1986).
equations.
equations
I
with
Comparaison
model. AIAA J. 24, No. 10
to the upwind
finite differencing.
OF
(‘OMPU
I INCi
IN
S(‘IEU(‘E
hNI1
l:NuCilNttKING
1,
64-02
( 1Y8Y )
A Numerical Method for Solving the Compressible Navier-Stokes Equations S. BOIVIN
D6partemeni de Mathimatiques, Universitk Laval, Laval, Quebec, Canada GlK 7P4. and INRIA. Menusin. 78153 Le Chesnay, France S. Boivin, A Numerical Method for Solving the Compressible Navier-Stokes Equations, IhlPACT ofComputing in Science and Engineering 1, 64-92 (1989).
A new computer code solving the compressible Navier-Stokes equations is described. The code solves the two-dimensional (2D), time-dependent equations using the finite-element method. Filtering methods are used on gradients of convection terms to ensure stability. Two turbulence models are integrated into the code to compare, evaluate, and ultimately improve model performance. The results of numerical experiments for the flow around an NACAO012 airfoil and a simplified version ofthe Hermes space shuttle are presented to show the possibilities of the methods used.
@I1989
Academtc
Press. Inc.
INTRODUCTION
Numerical simulations of flows around airplanes and space vehicles at high Mach and Reynolds numbers are still a challenge to the numerical analyst. Most of the difficulties arise from the simultaneous presence in the flow of complex physical phenomena of very different scales. The most frequently encountered are shocks, boundary layers (both laminar and turbulent), wakes, and various interactions between these phenomena. Throughout this paper, we hypothesize that all these phenomena are adequately modeled by the compressible Navier-Stokes equations. In particular, we assume that these equations are valid for highly turbulent flows. Under this hypothesis, our main problem is how to solve the Navier-Stokes equations without overdamping some of these phenomena. The strategy of resolution presented is a step in that direction and is based on the following key ideas: 64 0899-8248/89 $3.00 Copyngbt
0
1989 by Academac
All rights of reproduction
Press. Inc.
m any form resewed
COMPRESSIBLE NAVIER-STOKESEQUATIONS
65
first, being unable to represent the smaller physical scales, we take account of them through a viscosity computed with appropriate turbulence models, ??
second, unresolved frequencies following a shock must be damped by a numerical viscosity taking into account the physical viscosity, ??
third, with no reliable a priori information about the boundary layer, we seek to resolve it directly, or eventually, by the use of special elements. ??
Although very. appealing, no wall law is used because even if the region where viscous effects are important is small, its influence is felt everywhere in the flow. The formulation, schemes, and method of solution were also chosen to achieve some particular goals. In fact, it is our opinion that a good twodimensional (2D) Navier-Stokes code should: (i) be directly generalizable to 3D case, (ii) include a turbulent closure model with transport equations, (iii) be independent
from the mesh,
(iv) enable control of the numerical diffusion, (v) have solid theoretical support. In this report, we present a resolution algorithm which includes most of these features: (i), (iii), and (v) through a finite-element approach, (ii) and (iv) through adequate use of classical methods. The results of numerical experiments are presented to show the possibilities of the methods discussed in this report. 1. FORMULATIONOF THE NAVIER-STOKES EQUATIONS FOR HIGH REYNOLDSFLOWS It is well known [l] that with few initial data, the Navier-Stokes equations possess a unique solution. But in the usual case of many stiff initial data, instabilities appear and degenerate in turbulence. We make the usual assumption that these chaotic flows are within the range of validity of the NavierStokes equations (albeit not within the capabilities of our numerical approximation of these equations). Let D E RN (N = 2 or 3 in practice) be the flow domain and F be its boundary. The conservative form of the equations is given by
$+v.(pu)=o ~+v.(puc3u)=v.(7),
(1.1) 7ij
=
7ji
(1.2)
66
S. BOIVIN
(1.3) 711=
-Ph, + adub4,
(1.4)
+ &S,,(u)
P = (Y - l)w,
(1.5)
where S,j(u) = $( u,,, + Uj,i), E = e + $u2, e = C,T, y = C,/C,, q = -kVT, x = -&L. Despite all the recent advances in computer technology, turbulent flows cannot be computed at present by the direct resolution of the Navier-Stokes equations. The reason is that the turbulent motion contains a continuous range of scales in which the ratio of the larger to the smaller scale is about 0(Re3’4). For this reason, and because one is generally not interested in complete details of the behavior of all the instantaneous variables, we work with the mass-weighted averages of those quantities. To clarify, we recall the definition of the averaging operators. If p is the instantaneous density and u the instantaneous velocity, then we define the average of p and the mass-average of u on a time interval of length to by
p”
s 10
to
0
pdt,
a=& Pto
s 10
0
pudt
and u’by u’ = u - 6. To simplify the writing, we will retain the same notation for the mean variables and use the prime for the fluctuating part of mass-weighted variables. Consequently, for the remainder of the paper p, p are mean density and pressure; u, e, Tare mass-weighted mean velocity, internal energy, and temperature; u’, T’ are the fluctuating part of the velocity and the temperature; and u”is the mass-weighted averaging operator. Now, neglecting buoyancy effects and all viscosity and thermal conductivity fluctuations, the open equations for mass, momentum, and total energy conservations are
$+v.(pu)=o
(1.6)
~+v.(,,@24)=..((7_(,uGL~)) y
+ V.(pEu)
= V.[(T
- (pu’&‘)).
u - qtol],
(1.7) (1.8)
where pE = pe + $pu2 + ~PU’~ = pe + jpu* + pk, k being the turbulent kinetic energy, and qtOlthe total (laminar + turbulent) heat transfer.
67
COMPRESSIBLE NAVIER-STOKES EQUATIONS
The system ( 1.6) - ( 1.8 ) contains only second-order moments as supplementary unknowns, higher order moments being neglected and correlation between velocity and temperature =
modeled by
qi,tot =
-k(aT/dXi)
+
G
-k&(dT/dXi).
There are two ways to express the Reynolds stress tensor R = -pu'% u': To assume that the Reynolds stresses are a local property of the mean flow and, by analogy with the viscous stresses, express them by a constitutive law. ??
To assume that the Reynolds stresses are independent variable quantities and obtain their values through the resolution of a set of transport equations. ??
In this report, we consider only the first approach, the second requiring too much empiricism [ 21. The choice of a constitutive law relating R to mean quantities is very difficult. But, under certain approximations, Chacon and Pironneau [ 31 make it clear that, in the incompressible case, R' and (VU + Vu’) generate the same subspace of second-order tensors; hence in 2D R’ = cd+
B(Vu + Vu'),
(1.9)
where (Y,/3 are only functions of the nontrivial invariants of (V u + Vu’), and possibly of the turbulent kinetic energy and mixing length. We suppose that this can be generalized to the compressible case as
R = --put&d = -$(pk
+ /.L~V.u)l+
/&VU + Vu').
(1.10)
More general hypotheses [4] seem to give more realistic approximations, but without any proof that the Reynolds stress tensor still lies in the same space of tensors. In a classical way, we completely decouple the resolution of the mean Navier-Stokes equations ( 1.6 )-( 1.8 ) and the calculation of the Reynolds stresses. The next section is devoted to the resolution of the mean Navier-Stokes equations and Section 3 to the modeling and resolution of the model giving the Reynolds stresses. 2. RESOLUTIONOFTHENAVIER-STOKESEQUATIONS FORTHEMEANFLOW
We now consider the following nonconservative the Navier-Stokes equations ( 1.6 ) - ( 1.8 ) ,
nondimensional
form of
68
S. BOIVIN
ap
dr+u.vp+pv*u=o
(2.1)
$+(u.v)u+h9-~v.(“*o)=0
(2.2)
P
~+U.vTt~v.U-Vu: P
( 1-h.(K*vT)=o. V*
-0.
P
P
(2.3)
In (2.1)-( 2.3 ), we have normalized by the subscript r (i) the density p by pr (ii) the velocity u by I ur I (iii) the internal energy e by ( ur 12 (iv) the pressurep by p,(u,l’ (v) the viscosities p and pt by pr
(vi) the kinetic turbulence energy k by I ur 12 (vii) the temperature
T by I ur 1*/C,
which imply e = T. The pressure obeys the ideal gas law; the number y and the functions u, 8, v*, K* are defined by: ??
a=vu+vut-
#V.uI,
??
0 = p + $pk (p = (y - l)pT),
??
y = C,/ C, is the ratio of the specific heats ( y E 1.4 in air),
v* = (p + pt)/Rer is the total viscosity defined from the computed laminar and turbulent viscosities divided by the reference Reynolds number Re, = P~w%IcL~, K* = y/Re,( p/Pr + pt/Prt) is the total conductivity coefficient, also defined from laminar and turbulent viscosities. ??
??
Remark. The field (v* )-’ will be referenced as the total Reynolds num-
ber field. We consider external flows around 2D geometries; the domain of computation is described in Fig. 1a. Let r, be a far-field boundary of the domain; we introduce
r, = {XIXEr,,U,vz
r; = r,.w;,
where u, denotes the free stream velocity and n the unit vector of the outward normal to r.
COMPRESSIBLE NAVIER-STOKES
a
b
c
d
EQUATIONS
69
“t U “2 “3
FIG. 1. (a) Computational domain. (b) P, -iso Pz element: @ velocity, density, temperature, q and w d.o.f.; 0 velocity, q and w d.o.f. (c) Upwind finite element T(s), to node s. (d) Upwind finite elements T(s,), to nodes s, of a given element.
We assume the flow to be uniform at infinity, and the corresponding variables to be normalized by the free stream values; then, for example, we prescribe at infinity (Yis the angle of attack, P=
1,
T= Tm.= 1lI-d~ - lMb1, where it4, denotes the free stream Mach number. The boundary conditions on the computational boundary I?, are on FL :
u=u,,T=T,,p=
1,
and on r”,:
VT.n=O,[n.(Vu+vut-fv.Uz).n]=O,
[n*(Vu
+ Vu’ -
$0. uZ)* t] = 0,
where n, I are the local normal and tangent on I’:.
S. BOIVIK
70
On the rigid boundary rB, we shall use the conditions ( no-slip condition),
u=o T = TB = T,(l
+ (y - 1)/2M:,)
(free stream total temperature).
Finally, since steady solutions are sought through time-dependent equations, initial conditions must be added; we shall take P(X, 0) = PO(X),
T(x,O)= To(x).
u(x, 0) = e(x),
Solving the compressible Navier-Stokes equations is a difficult task. Most of the existing numerical solution methods are based on finite-difference techniques, for both space and time discretizations. Following the work of Bristeau et al. [ 5, 61, we will consider a method based on finite-element techniques for space discretization while using finite differences in time. Let the time derivatives be discretized using a classical implicit Euler finitedifference formula; then at each time step, we solve the nonlinear system of variational equations ol(p - ;, N) + (U-VP, N) + (~0. U, N) = 0
(2.4)
LY(U-G,M)+((U’V)U,M)+
+(FVT,VK)=O,
(2.6)
where the solution is sought in V X W X Z, Y=
EH’(n)lp,r,
w=
{uE(H’(Q))+Q.~
Z=
{TEH’(O)lT,r,
although there is no triple (p, u, T) E V of test functions (N, of the spaces X(Q), Remark
{p
= 1) = U,, UjrB = O} = T,,
TI,., = TB),
complete theoretical justification for this. We look for a X W X Z such that (2.4)-( 2.6) is verified for all triples M, K) E Y(Q) X (X(Q))* X x(Q). For the definition Y(Q) see Appendix: Nomenclature.
1. ( , ) denotes the scalar product in L2( Q).
COMPRESSIBLE
NAVIER-STOKES
EQUATIONS
71
Remark 2. The natural boundary conditions already defined were introduced by setting the boundary integrals appearing from the integration by parts to zero. Remark 3. The 0 term of Eq. (2.5 ) is integrated in the form
1
l)vT+$vq+ (y- 1)f+2dVP.
-+(r-
3P
[
We do not integrate by parts the pressure term of the momentum equation because we cannot set p + 6,, = 0 on the boundary (doing so would perturbate the solution) and we prefer to avoid the computation of boundary integrals. To approximate (2.4)-( 2.6), by the finite-element method, we must divide 52into small elements (triangles) and replace all the functions by their interpolate oh, uh, Th. Interpolates are defined inside the elements from their values at the nodes of the elements by an interpolation formula. The choice of the approximation spaces is a critical and still partly obscure point. But, we see that when u is very small, for example, near a solid boundary, Eq. (2.4) looks like a divergence free equation. This is especially true after the change of variable: (I = log p (see [ 5 ] ). We may conjecture that to ensure stability, a necessary condition is that the element chosen should be stable for the incompressible case. This seems to be confirmed by the various experiments made by Bristeau et al. [ 61. A possible choice for such an element (and consequently spaces) is the PI -iso P2 element (Fig. 1b). If Th12is the triangulation obtained from Th by dividing each triangle of Th into four equal triangles whose vertices are the vertices { qr, q2, q3} of Th and ( h(ql + q’), f(q2 + q3), f(q’ + q’)}, we define
v, = {~/a E VIP~ E c’(Q), VTE
Th, PAITE f”(T)} C
w, = {uh E WjUh E (C0(fi))2, VTE T/,,z, Uh,TE P’(T)}
v c w
Zh= {ThEZIThEC”(n),VTETh,Th,TEP’(T)}CZ, where P’ ( T) = { set of polynomials of degree d 1 on T} . Restricted to these finite-dimensional spaces, Eqs. (2.4)-( 2.6) lead to the following nonlinear problem: Find (ph, uh, Th) in (VAx wh x z,,) a
SOhtiOIl
of Fh(P,,, u,,, Th) = 0,
Fh being the discrete version of the system (2.4)-( 2.6).
(2.7)
s. ROlVllU
72
We now consider iteration schemes for solving the nonlinear system F,?(s ) = 0, where s = (oh, ul,, T,,). Newton’s method applied to this system results in the iteration 1. Set so an initial guess 2. For IZ= 0, 1, 2, . . . , until convergence do 2.1. Solve J(s”)6” = --F,,(Y), (2.8)
2.2. Set sn+’ = sn + 6”,
where J( s”) = F’(f) is the system Jacobian. For large problems, iterative methods are frequently used to solve (2.8) only approximately, giving rise to methods which can be viewed as inexact Newton methods. The particular method we use is the Generalized Minimum Residual Method (GMRES) due to Saad and Schultz [ 71. This method has the virtue of requiring virtually no matrix storage and only the action of the Jacobian matrix J times a vector r, and not J explicitly. In the setting of nonlinear equations, this action is approximated by a difference quotient of the form J(s)r
z
F/,(s+ br) - F/,(s) h
3
where s is the current approximation of a root of (2.7) and b is a scalar. For details of the algorithm, see [ 71 and also [ 61 for a clear setting of this algorithm within the context of the conjugate gradient methods. To conclude this section we present and discuss two strategies to add some numerical diffusion to the approximation of the convective part of the equations. It is well known that when the Reynolds number is high, boundary layers and shocks become hard to resolve and centered approximations of first derivatives are dangerous. In fact, it is necessary to add some numerical diffusion where the mesh size h is such that the local Peclet number satisfies the condition Pe=@>2, C
a = P max(Iuil,
lull),
to damp unresolved frequencies. There are many different ways to do that, upwinding, addition of an artificial viscosity, etc., but the problem is to control the numerical diffusion to be sure it does not degenerate the approximation of gradients. Our approach is to filter the discrete gradients in the following way: as the discrete gradients are constant by element, we define the gradients at a node
COMPRESSIBLE
NAVIER-STOKES
EQUATIONS
73
or on an element by an upwinded combination of the gradients on the nearby elements. The first method can be viewed as a projection followed by a combination with the upwinded gradients, that is, if a node s is surrounded by elements 1 to n , with respective area ak which sum to A. Then we define the gradient of a function 4 at s by the combination of a centered gradient and an upwinded gradient
VA = (1 - /w4: + PV4F,
(2.9)
where the centered and upwinded gradients at s are respectively defined by
where T(S) is the upwind element at node S. The choice of p is discussed below. The second method is a pure filtering and can also be written as a combination of a centered and an upwinded gradient. More precisely, on an element e we define the gradient by
vdJt== (1 - P)V4c,+ PVK,
(2.10)
where the centered and upwinded gradients on e are respectively defined by
A fundamental difference between these two upwinding methods (for elements where gradients are constant by element) is that in the second there is no need for projection. We now describe the computation of p. The method, based on the analysis of the scalar steady convection diffusion equation, is to relate the coefficient to a local Peclet number by 0 = coth F - ;, ( 1
where the Peclet number is an estimation of the importance of convection compared to diffusion. Although the definition of this number is obvious in the scalar case, the generalization to multidimensional cases is not immediate. But in general we may define it by
74
s. BOIVIN
Pe = Re phi, where h is an estimation of the characteristic length and a”an estimation of the convection speed. A possible choice for A, a^is it =
I(~~~)lll~l,
a= IUI.
(2.11)
But some computational experiments demonstrated that the following choices, inspired from the work of Hughes et al. [ 81, give better results, Ji = I(~~~~)lll~~l,
6 = l(~~~~)Ill~~l.
(2.12)
In fact, with these definitions, /3 is more closely related to discontinuities. Note also that the coefficient ,f3can be multiplied by a tuning coefficient for further control of the diffision and more specifically to weight the diffusion added in each equation. The equation for p being the only equation without any physical diffusion it is fair to believe that more numerical diffusion should be added to this equation. Although less robust, the second method, the filtering method (2.9)-( 2.12), proved to be very efficient and precise, provided that the mesh was adequately refined. Remark. The value of h = (h,, h,) at a given point of the mesh is deduced from the computed value on each element using the relations
where qi are the PI basis functions on a given triangle.
3. CLOSUREMODELOFTURBULENTFLOW The basic goal in developing turbulence models is to specify closure conditions in which the unknown Reynolds stresses of turbulence are related, either algebraically or through differential relations, to known mean-flow variables. The fundamental problem is to deduce a suitable mathematical model which should be a good compromise between mathematical simplicity and the complex exact transport equations of Reynolds stresses. Our starting point will be the classical k - c model of Jones and Launder [ 91. A reason for it is that this model can be deduced rigorously in the in-
COMPRESSIBLE
NAVIER-STOKES
EQUATIONS
75
compressible case [ 31 and also because there are some efficient (but not trivial) ways to solve it. It is well known [lo] that the choice of the scheme to solve the k - c equations is very important. Thus, to ensure numerical stability we make a change of variable inspired by Coakley [ 111. The basic closure turbulence model as described under Appendix: TwoEquation Closure Model is given by
am) at
+u.v(Pq)+;pqv.u--IR:vu+~pwq 37
-v.[(,+:W]=o
ah4 +
at
u*V(pw)
+ pwv-
u - (Cl - l$R:
(3.1)
vu + (c2 - l)pw2 -v.[(p+~~~]=O,
(3.2)
where q2 = k, w = e/k,pt= c,,pq2/w, and R is the Reynolds stress tensor. In the construction of the two-equation closure model it is assumed that the flow is fully turbulent and that the Reynolds number is everywhere high. Thus, to provide predictions of the flow within the viscous layer adjacent to the wall, the above model must be extended in three ways. These are: (i) prescribe the turbulent / nonturbulent
interface,
(ii) the terms containing the cl and c, will become dependent upon the Reynolds number of turbulence, (iii) if possible, correct values using empirical formulas. The complex problem of the interface between a turbulent region and a nonturbulent region will not be discussed here. The corrected values for the constants, proposed by Coakley [ 111, are given by replacing c, by c,D and (1 - cl ) by (0.405 D + 0.045), where D is a damping factor given by
D = [l - exp-“‘I,
r=
p42 110 ’
where the value of the constant is a! = 0.0018. Initial profiles for q and w and corrected values are deduced from an estimation of the turbulent viscosity provided by an algebraic model, namely the Baldwin-Lomas [ 12 ] model, through the empirical relations
S.
76
BOlVllv
(Bradshaw hypothesis) 3 w = -s,* 10
( Hypothesis: production
= dissipation).
As in the case of the mean Navier-Stokes equations we work with the nonconservative and nondimensional form of these equations. That is, ~+~.V~-~~V.u-~c,D(a:Vu)+~wg-~V.[C,Vql=O 2
(3.3)
+ u*vw + ; (c, - 1)oV. u - (c, - I)c,D(a:
Vu)
+ (c* - l)w2 - A V*[C,Vw] P
= 0.
(3.4)
In (3.3)-( 3.4), we have normalized by the subscript r (i) p, u, T as in Section 2 (ii) the viscosities p and pt by pr (iii) the kinetic turbulence energy q* by I u, I ’ (iv) the specific dissipation rate of q*, w by 1u, I/L,. The functions D, 6, C,, C, are defined by: ??
D = 1 - exp(-apq*Re,/w),
??
a=Vu+Vu’-
??
C, = (CL+ d~y)/Rer,
??
C, = (P + dcJ/Re,,
~04,
The values of the constants c,, cl, c2, ug = uk are taken from the original paper of Launder and Spalding [ 13 J; for the value of 6, see the first appendix. The boundary conditions on smooth wall boundaries are q = 0 and VW. n = 0. At the computational boundary I’, small values are prescribed: q = qa , w = w, . More exactly, qa is set to the known free stream value of the kinetic turbulent energy and w, is defined by setting Put=
Re,c,p,&
= I
W,
We now describe the resolution strategy. As already noted, the overall resolution is split in two parts; first we solve the mean Navier-Stokes equations and then we solve the turbulence
COMPRESSIBLE
NAVIER-STOKES
EQUATIONS
77
closure equations. At each step the resolution of (3.3)-(3.4) is done with p = i, u = u^,P = ji, that is, with known values from the preceding NavierStokes step. The technique used to solve (3.3)-(3.4) is very similar to that used, in Section 2, to solve the mean Navier-Stokes equations and is based on the following weak formulation. Let the time derivative be discretized using a classical implicit Euler finitedifference formula; then at each step, we solve the nonlinear system of variational equations
+;(wq,N)+($‘q,VN)=O
a(w-i,K)+(ri*vqK)++ - (c, - l)(c,Di
(3.5)
l)(wV*ti,K)
: Vii, K) + (cz - l)(w $K)+($
VW, VK ) = 0,
(3.6)
where the solution is sought in Q X S,
Q = {QE H’W)lq,r, = q-3
s = {w E H’(O)lo,r,
qlrB =
0)
= cd,}.
Here again, we obtain an approximation of the solution using the inexact Newton-GMRES method on the nonlinear system (3.5)-( 3.6) considered on the finite-dimensional approximation spaces
Qh = {Q/, E Qlsh E cot% &=
b’7-E l-h/z, 4h,rE P’(T)}
C Q
{OhESIOhEC”(,),VT.Th,2,Wh(TEP’(T)}CS,
where P’(T) = {set of polynomial of degree <1 on T} and Th,2 is the fine grid of the P, -iso P2 grids system discussed in Section 2. To stabilize our scheme, we replace the convection terms (u^*Vq, N) and (11. VW, K) by filtered terms (see Section 2 for the filtering technique). The second model of turbulence we use is an algebraic model derived from the Baldwin-Lomax [ 121 model. We include this model as a basis for comparison and eventually to decrease the cost of computation. (Another reason was to try to combine both models but preliminary computations led us to the conclusion that the two-equation model gives better results.)
S.
78
BOIVIN
In this model, the turbulent viscosity is computed using a two-layer approach,
where Y is the normal distance from the surface and Y, is the least value of Y at which pt, = ptO. In the inner layer pti =
K’plw( Y2[1 - exp(-Y’/26)12,
where w is the vorticity, and Y + = Y ( p0 I c, 1/I;) ‘I2 with the subscript o referring to the value of the parameter at the nearest solid boundary, and mu is the local shear stress in the direction of the flow. In the outer layer
whereF,,,,,=max(Ylwl[l -exp(-Y+/26)]),andY,,,,,isthevalueofYat which F,,,,, occurs. The KIebanoff intermittency correction is given by Fkleb = [1 + 5,
~(~klebY/Ymd61-‘.
The model constants and functions are K = 0.40, k = 0.0168, 0.2933M,
C,, =
C kleb
I
ifA4, < 3; ifM,
2.08,
-0.4118A4, 0.3,
=
+ 1.2,
+ 0.65,
I
> 3,
if M, < 0.85; if M,. > 0.85.
In the wake region, the following nonclassical formulation is used,
Pt = CLtr(IWI/I%l), with the subscript r referring to a reference profile, usually the last profile before the trailing edge. This unusual wake formulation enables us to simulate the characteristic transport-diffusion behavior of the turbulence in the wake, keeping the emphasis on the evolution of the vorticity. This scheme avoids the appearance of steep gradients when moving from the wall model to the wake model of the usual algebraic models. Some details about the implementation of this model for unstructured meshes are given below.
COMPRESSIBLE
NAVIER-STOKES
EQUATIONS
79
This kind of model is one of the most widely used because of the ease of its implementation in finite-difference codes. The basic property appearing here is the orthogonality between the mesh lines and the walls. This important property will have to be replaced by an abstract structure for general finiteelement meshes. The most obvious way is to compute all the model parameters on a second mesh ensuring the orthogonality property and to use an interpolation operator between the two meshes. This procedure can be difficult to implement for complex (2D) geometries and inadequate in 3D. Another way to implement algebraic turbulence models for general meshes is to consider independently each node. That is, keeping some simple geometric characteristics of each node in a table, for example, the nearest point of each wall, the normal lines passing throughout this node can be found and so all the model parameters can be evaluated. But as usual, the simplicity of this method must be paid for; the maximization of functionals along the normal lines can be very costly. Consequently, the complexity of the process must be diminished by the use of approximate maximization techniques and eventually by restricting the process to adequately chosen nodes. These simplifications make the overall process interesting. 4. RESOLUTION
SCHEME
We summarize here the overall resolution scheme as described in the previous sections. This reads: Step 1: Read initial values of u”, p”, To or let p” = 1 and solve A.u” = 0, AT0 = 0 with appropriate boundary conditions. Step2:
Let q" = qcu,oo = w,(qpr, = 0).
Step 3: u”, p”, T”, pf, q",W" are known from the previous time step. Step 3.1: Update the laminar viscosity (P”+‘) using the Sutherland formula. Step 3.2: Compute the algebraic viscosity (&zi) using the BaldwinLomax model. Step 3.3: Make a correction to q" and J’, 4 “+*‘* = max[q”, qa],
cd"+"*= max[w”, w,]
or 4
n+1/*= max[(l(gzl$2)~'*,q.,
w n+t/*= max
+2,
[
an,
WCC
1 -
q_]
s. BOIVIN
80
Step 3.4: Solve the q - w equations F+( q" ’, w"+' ) = 0 (see Section 3). Step 3.5: Update the turbulent viscosity (this step may include a projection of q and w to eliminate negative values),
or p;+’
=
max[c,D(q”+‘)‘/o”+‘,
&,:I.
Step 3.6: Solve the N-S equations FNs( z~“+‘, p”+‘, Tnf’) = 0 (see Section 2 ) . Step 4: If Is”+ - s” ]L~cr2)< c, sn = (p”, un, T”), then stop, else go to step 3. 5. NUMERICALEXPERIMENTS Before the exposition of the results of the numerical experiments let us recall that as v* is defined from the laminar and the turbulent viscosities, it is not a fluid property. This field is of primary interest; more exactly we are interested in the field 1/v* called the total Reynolds number. Computations were performed for some flow condition around an NACAOO12 airfoil and a simplified version of the Hermes space shuttle. In each case, the steady-state solution was reached after 80 to 120 time steps depending on whether the initial flow field was uniform or a known flow field from a previous computation. Two test cases were chosen here for the computations, each case involving complex flow with shock waves. The first series of numerical experiments concerns a compressible viscous and turbulent flow around an NACAO012 airfoil. A part of the coarse mesh used is shown in Fig. 2a. We took here M, = 0.85, Re = 105, y = 1.4, Pr = 0.72, Pr, = 0.9, and a zero angle of attack. We have shown in Fig. 3 the iso-Mach contours (Fig. 3a) for a NavierStokes computation with the total Reynolds number field deduced from the Sutherland formula (Fig. 3b). The same problem was solved a second time but using the two-equation turbulence model, discussed in Section 3, for the computation of the total Reynolds number field. See Fig. 4 for the iso-Mach contours and the total Reynolds number field and Fig. 5 for the iso-k and iso-w contours. Despite the problems caused by the coarse mesh in the wake, we note the presence of vortices in both computations. In the turbulent one, the boundary layer thickened abruptly but did not damp out the time-dependent nature of the shocks and the wake vortices.
COMPRESSIBLE NAVIER-STOKES
mesh: 11~30012
- 2000
mesh: pseudo-hermes
-
nodes
1883
nodes
81
EQUATIONS
Univ. Lava1
Univ. Lava1
27 Jul 88
27 Jd 88
FIG. 2. Partial view of the principal meshes.
It is clear that this flow is dependent on the capabilities of the mesh we use. But it is interesting that most of the physical phenomena we tend to see are present. An obvious way to obtain better results is to refine the mesh in regions where steep gradients are present. This is only a matter of computer facilities, but we are now developing a strategy to obtain better results via the use of enriched finite elements in the corresponding areas. The second series of numerical experiments concerns a compressible viscous
S. BOIVIN
iso-Mach
MAX=O.l78e+Ol
MlK0.485~~19
univ. Lawtl 27 Jul88
,
FIG. 3. Mach contours and total Re number; laminar computation at A4 = 0.85, Re = IO>.
and turbulent flow around a simplified version of the Hermes space shuttle. A part of the coarse mesh used is shown in Fig. 2b. We took here M, = 2.0, Re = 104, y = 1.4, Pr = 0.72, Pr, = 0.9, and a zero angle of attack. Note that the length of the body is approximatively 0.14 of characteristic length. We have shown in Fig. 6 the iso-Mach contours (Fig. 6a) for a NavierStokes computation with the total Reynolds number field given by the Sutherland formula (Fig. 6b). The same problem was solved a second time using
COMPRESSIBLE
iso-Mach
NAVIER-STOKES
MAX=O.l67eGIl
83
EQUATIONS
MIK0.485e19
Univ. Lad 26 SW
06
\ Re_total
MAX=0.1221?+06
MIN=O.l95e*03
Univ. Lava1 26Sep66
nc.
4. Mach contours
and total Re number:
q - w computation
at M = 0.85,Re = 105.
the two-equation turbulence model (see Fig. 7 for the iso-k and iso-w contours), but all the characteristics of the flow remained similar (see, for example, a comparison of the friction coefficients for the turbulent and nonturbulent computations shown in Fig. 9b). This is encouraging because at this low Reynolds number the flow should be mainly laminar. We have conducted a sensibility experiment on the values of the constant a, of the w equation. The new iso-k and iso-w contours for the last q - w computation with a, = 13
S. BOIVIN
84 /
I
I
I
k
MAX=0.625e-01
w
MAX=O.l08eO3
MIN=O.OOOeOO
MIK0.9OOe-02
univ. 27
Lava1 Jul68
univ. Lava1 27
Jul98
FIG. 5. k and w from a q - w computation at A4 = 0.85, Re = 105.
instead of 1.3 are shown in Fig. 8. Although satisfactory for a very wide range of flow conditions [ 141 preliminary results on simple nearly incompressible cases tend to indicate that this choice for 0, leads to overdiffised values for k and q* and w. That is, the constant u, is not universal. Finally, a numerical experiment was conducted to compare the diffusion associated with the turbulence model used. The test case was the transonic flow around an NACAO012 airfoil at AI, = 0.85, Re = 104, y = 1.4, Pr
COMPRESSIBLE NAVIER-STOKES
iso-Mach
MAX=0226e+Ol
85
EQUATIONS
MIrG=O.441e-19
Univ. Laval 27 Jtd 66
////////I
I-_ .~a
MAX=O.l lOed
MltWS23e94
Univ. Lava1 27 Jul66
FIG. 6. Mach contours and Re number; q - w computation at M = 2., Re = 104.
= 0.72, Pr, = 0.9, and a zero angle of attack. Only the resulting values of the coefficient of friction (see Fig. 9a) are discussed. The Navier-Stokes computations (inner curves) show a large separation zone while the use of the algebraic turbulence model gives no separation at all (outer curves). It is clear that the diffusion caused by this turbulence model is overestimated. The estimation given by the two-equation model (see middle curve in Fig. 9a) is more realistic.
S. BOIVlN
86
! /
k
/
w
MAX=O.445~02
MIN=aQOOe03
univ. Lava1 27 Jo188
FIG. 7. k and w from a q - w computation at A4 = 2.0, Re = IO?
The choices made in steps 3.3 and 3.6 of the resolution scheme (see Section 4) were proved to be less important than expected. In fact, the solution process described to solve the turbulence transport equations with the given boundary and flow condition converge without the need of the algebraic turbulence model. The latter may be used to start the q - w computation from a given state, causing the steady state to be reached more quickly. We present here some statistics to give one an idea of the time used for these computations on a SUN3 workstation with a floating point accelerator
COMPRESSIBLE NAVIER-STOKES
87
EQUATIONS
MAX=0.621e-O2
I
w
MAX=O.566&2
MIKa900e03
W.
Lava1
30 JUl88
FIG. 8. Same as Fig. 7 but with a different constant of diffusion in the w equation.
(0.4 Mflop). For a computation on the 1883-node simplified Hem& shuttle, it takes at each time step:
space
110 s to compute the turbulent viscosity with the algebraic model, 2300 s to compute the turbulent viscosity with the two-equation model (complete resolution at each time step), 1100 to 2000 s to solve the Navier-Stokes equations.
S. BOIVIl\i
88
0
02
0
01
0 00
-0
01
-0.
02
-0
03 1 lII/,IIII
0 00 b
0. 04 0 03
0 02
0 01
0 00
-0
01
-0
02
-0
03
i
i
I,,,
0 25
0. 50
I,,,
0. 75
1 00
-.--. _=-_ _-_ ---_ _ _ .I:.... -I..._..I _..........~ _.
‘.;=
__---__-_--
_---
_--_
--
---
-
-0.04
0 00
0.03
0 07
0. 10
0. 1
FIG. 9. Comparison of friction coefficients;see text for details.
More numerical experiments and further details about the method can be found in a research report [ 151 published at INRIA, France.
6. CONCLUDING REMARKS We have presented a numerical scheme to solve the compressible NavierStokes equations written in nonconservative form with emphasis on turbulence modeling. We have obtained good results for moderate Reynolds numbers using an approximation satisfying some “inf-sup” condition and using a filtering of the gradients of the convective terms. This filtering is dependent on
COMPRESSIBLE
NAVIER-STOKES
EQUATIONS
89
the local direction of the fluid flow and on the local amount of physical viscosity. A particular form of the k - c model was solved with standard boundary conditions to obtain an estimation of the turbulent viscosity using the same techniques. Further developments are needed to obtain accurate solutions for higher Mach and Reynolds numbers where special care must be taken to not introduce too much numerical dissipation. In the near future we plan to compare computations on finer grids with experimental tests and emphasize filteringupwinding methods and also the resolution of the boundary layers. APPENDIX
Two-Equation Closure Model We present here some details about the two-equation closure model ( 3.1)(3.2) that we used. Let us recall the simplest and certainly the most popular closure model, namely, the k - Emodel, where k is the turbulence energy and e is its rate of dissipation. According to this prescription
a(Pw
-+uV(pk)=-pkV.u+R:Vu+V.
at
am at
+uV(ps)=-p67~u+c,fR:Vu+V
[(p+$Vk]-pe .[(P+37c]-c2P~?
(1) (2)
where c,, (lk, ur, cl, c2 are model constants, usually set to the recommended standard values 0.09, 1.O, 1.3, 1.44, and 1.92, respectively. It is well known that near a solid boundary (at the boundary we suppose k = t = 0) this model needs a low-Reynolds-number dissipation term to balance molecular diffusion in the sublayer. In fact, letting y tend to 0, Eq. ( 1) becomes
which is incompatible with the experimental fact that k - y2 in the sublayer. Instead of adding empirical terms to correct it, we assume k to be constant in the diffusion term and make the following change of variable k = q2, from which we deduce directly Eq. ( 3.1). Note that if we assume q2 = k - y2 we obtain the correct behavior of q when y + 0. Finally, to reduce the problems of numerical instabilities, we reduce the truncation errors due to quotients, using the change of variable w = c/k. But this must be paid by supplementary hypotheses.
S. BOIVIN
90
If we insert w = c/k in Eq. (2) and use the identity k(pw)’ = up’+ (PC)’- u(pk)’ we find the equation
ah.4 -=w~+wu.~p-~.V(pw)+(~~-l)wR:Du
at
42
+ D - (cz - l)pw2,
(3)
where the diffusion term reads D=V.[(p+;-Ve]-wV-[(p++‘k].
(4)
Thus we need to assume (pt + ~*/u,) = (p -t p,/c~k) = cte in order to write it in the simpler form
where the recommended value for Q, [ 1 l] is 1.3. But, although satisfactory for a very wide range of flow conditions [ 141 preliminary results on simple nearly incompressible cases tend to indicate that this choice for u, leads to overdiffused values for k = q* and w. The standard form for this equation follows from the use of the nonconservative mass balance equation.
Nomenclature Letters incompressible Reynolds stress tensor compressible Reynolds stress tensor viscous and pressure stress tensor viscous stress tensor strain tensor laminar viscosity turbulent viscosity reference viscosity turbulent kinetic energy
'smog
snoas!A a[q!ssaJduroau! pue aIq!ssaJdwoa JO uoyelnmp aql 01 uoyagddy
S=fOlS-Ja!AEN ‘( L861)
gLI
aql 'OJ SpOqKK4I .zpqq
p!nld T
ll23!JZNUIlN ‘XllE~?d
'f PUFJ ‘!~SU’MOI~
‘a
.suogenba
‘nl.?a~S&I~ ‘0
‘m
‘s
.aaualnqJnijo stapolu 3 - y pua I-? Jeauquou u0 ‘afe!zadS ‘S ‘3 ‘p
‘(9861) YJoA MaN/u!lJatl ‘%‘aA “XII aJt?M$jOs uoyezyqdo ‘S3!lVlUJylV~ pqddy u! SVLVA ‘(‘pa) ueuqs&~+?[s~ ‘A ‘V u[ .lapom 3 - 9 aqlJ0 uoyJpunoJ Iea!i~maqleur aql u0 ‘nsauuoJ!d ‘0 pur! uoaeq3 ‘1 ‘E
-Ja%U!JdS
3puepaqiaN aqL ‘qaJeas -aa
‘OJpilH
‘OSSV
‘WUI
'SJ!lllVlpkH U!
uo,m,qddy
.l!JlfLpUV
Slapow 't2861)
-qUoN ‘I\
'~U!.IJJU@U~
uo!lour JO suogsnba
aql
pun
SJ~UJ!~~
JOJ SLW~qOJd
pu= ‘IOPJQ ‘8 ‘nEaUUOJ!d aM ‘uo!l’+V la+J@lln==a
‘0
pqddy
uf SPO~IJM
an@AiClEpUnOq
‘!‘lSU!MO[f) [a3JEm
‘u
SUO!AV
‘U!lJOJ “OJJ
Sqndtuo3
[qpq
‘m
‘Ep!qS!N
J3UJlTlq.Vl~ ‘!pOx
'M
'2
uJePJal=W ‘PuelloH y ‘spq ‘1
pue
lEJaua;l JO
eJeuInsle~
‘v
1
3uo!ssnas!p pg!nJJ JOJ xnepad ‘1 ‘0 ‘m xueq$ LIIe!aadsa pauoddns L[&wds! ~JOM s!qL
OSlC PUP nW$S!Ja
KWJ1UO3 l? bq
*sa~gt?“~~appue saxpu! pawadaJ “03 uo~y.~a~uo~pa)y~ -0sse aqt put? uogwou losuat ~txymp ayl asn ahi ‘papaau uayM ym4.4a~
sa+w+ap
lsry
‘10 = EL”XI(0),H 3x1 = (0)x (0 = =“‘xl(D),H 3x1 = (0)x (z&7 u;[ ov a= (7~)~7 30 suogx1n3 30 axds Aaloqos pmsn (0 ) , H (8 )I 7 uoyun3 alqt?Salur an%saq~ aanbs 30 axds
asoym
16
SNOI.LVL%XI SEfXO_LS-XIIAVN 318ISS3lIdI403
92
S. BOIVIN
6. M. 0. Bristeau, R. Glowinski, B. Mantel, J. Periaux, and G. Rage, Adaptative finite element methods for three dimensional compressible viscous flow simulation in aerospace engineering. In 11th Inter. Conf: on Numer. Meth. in Fluid D_vnamics. Springer-Verlag, Berlin/New York. in press. 7. Y. Saad and M. H. Schultz, GMRES: A generalized minimal residual algorithm nonsymmetric linear systems. SIAM J. Sci. Statist. Comput. 7 (1986).
for solving
8. T. J. R. Hughes, M. Mallet, and A. Mizukami, A new finite element formulation for computational fluid dynamics. II. Beyond SUPG. Compul. Methods App/. Mech. Engrg. (1985). 9. W. P. Jones and B. E. Launder, The prediction of laminarization of turbulence. Int. J. Heat Mass Tram@ 15 (1972).
with a two equations
model
10. A. G. Hutton, R. M. Smith, and S. Hickmott, The computation of turbulent flows of industrial complexity by the finite element method-Progress and prospects. In Finite Elemenfs in Fluids, Vol. 7. Wiley, New York (1987). 11. T. J. Coakley, Turbulence modeling AIAA Paper 83-1693 (1983).
methods
for the compressible
12. B. S. Baldwin and H. Lomax, Thin layer approximation turbulent flows. AIAA Paper 78-0257 (1978). 13. B. E. Launder and D. B. Spalding, The numerical Methods Appl. Mech. Engrg. 3 ( 1974). 14. T. J. Coakley, (1987).
Numerical
simulation
and algebraic
computation
of viscous transsonic
Navier-Stokes
model for separated
of turbulent
Navier-Stokes
16. B. Escande, Modelisation de I’interaction onde de choc/couche limite turbulente. calcul-experience. These Doctorat, Ecole Centrale de Lyon (1986). Separated
flow treatment
with a new turbulence
18. M. Tabata, Finite element approximation corresponding Mem. Numer. Math. Univ. Kyoto Tokyo Y (1977 ).
flows. Comput.
airfoil flows. AIAA Paper 87-046
15. S. Boivin, A finite element method to solve the compressible turbulence modelling. INRIA Research Report (1988).
17. U. C. Goldberg, (1986).
equations.
equations
I
with
Comparaison
model. AIAA J. 24, No. 10
to the upwind
finite differencing.