Application Of Two-equation Turbulence Models In Aircraft Design.

.

Application of Two-Equation Turbulence Models in Aircraft Design. G. Kalitzin, A.R.B. Gould British Aerospace (Operations) Ltd, Sowerby Research Centre, Filton, Bristol, UK.

J.J. Benton British Aerospace (Regional Aircraft) Ltd, AVRO International Aerospace Division, Woodford, Cheshire, UK.

Abstract

been created. This is because the Multiblock topology description embodies a generic speci cation of line-grid stretching so that small changes in geometry require no change to the topology. The primary ow solver within the Multiblock system is a 3D, Reynolds-averaged, compressible, cell-centred, central dierencing method known as RANSMB which has been developed by the third author. This makes use of Runge-Kutta multigrid time marching and arti cial viscosity to produce steady-state solutions. Both Euler and full Navier-Stokes calculations can be made and several turbulence closures are available, from the algebraic models of Baldwin-Lomax and Johnson-King to one and two-equation transport models. Whilst the underlying solution method was rst described 10-15 years ago1; 2 it provides an unparalleled combination of eciency and maturity. Development of the solver continues in many areas. However, one of the pacing items in its development and wider application is turbulence modelling. The traditional algebraic models, including BaldwinLomax3, have been applied with a surprising degree of success to two and three dimensional ows. Although such models are cheap and robust their direct dependence on empirical mixing length assumptions severely limit the complexity of the ows which can be adequately predicted. Transport models of turbulence would appear to be more appropriate for modelling ows around multi-component geometries in three-dimensions. Second Moment Closure is an advanced manifestation of this class of models. By providing a transport equation for each Reynolds-stress component, SMC models have the potential to provide predictions of complex ows where anisotropic `non-linear' eects are important. The penalties that accompany such a capability include high cost, due to the large number of additional equations, and potentially poor numerical properties which may lead to implementation diculties. Whilst SMC remains the long term goal for a generalised, ow modelling system,

CFD methods have been used successfully for many years in aircraft design. The primary Navier-Stokes method in use at British Aerospace is an explicit, structured, Multiblock system. One of the key targets in the development programme for this system is the implementation of a Second Moment turbulence closure. Before proceeding to implement such a model it is essential to have an ecient and robust two-equation transport capability. This paper describes how this capability has been achieved. The selection and assessment of twoequation turbulence models for use within this system is discussed. Two models are considered in detail, one from the k-" family of models, and one based on the k-! class. To provide stable and accurate solutions for complex three-dimensional ow elds various modi cations are required to both types of model. Several 2D and 3D civil aircraft applications are used to validate the system. These also serve to demonstrate the use of transport turbulence modelling in the design environment.

Introduction For many years British Aerospace has been developing CFD methods for use in the design of aerospace vehicles. One suite of programs which has emerged from this development is based around a 3D, structured, Multiblock mesh generation system. This system is targeted at the modelling of external ows over complex geometries and is now being used for design applications related to military aircraft, civil aircraft and high speed weapons. This forms part of a wider initiative to integrate CFD with CAD packages which exploits the fact that mesh re-generation for small geometry changes can be automated once the initial topology and grid have  Permanent address: Otto-von-Guericke-Universit at, Institut fur Stromungstechnik und Thermodynamik, Universitatsplatz 2, 39106 Magdeburg, Germany. c 1996 British Aerospace plc. Published by the Copyright American Institute of Aeronautics and Astronautics, Inc. with permission.

1

1. No use of normal-to-wall distance. Most three dimensional calculations contain areas where two or more boundary layers are interacting, such as wing-body junctions. It is dicult, therefore, to locate a suitable wall distance for cells in such a region. For a Multiblock scheme it is also possible that the upper part of a boundary layer will be separated from the solid surface by a block boundary. Locating the relevant surface in a generalised Multiblock mesh is a complicated and expensive process. 2. Simple source terms. The turbulence model source terms should be functions of readily available quantities. For example, models which make use of second derivatives are undesirable. Source terms which become extremely large very close to walls cause stiness problems. 3. Straightforward boundary conditions. The turbulence model variables should have boundary conditions that can be readily applied by setting `halo'-cell values. The use of more complex conditions, such as forcing the dependent variables to match a function over the near-wall region is undesirable since this imposes severe constraints on the grid generation process.

the technical challenges posed by such a model mean that an interim capability is required. Two-equation models can be seen as an advance on algebraic models, since they are transport in nature and therefore can take account of history eects. On the other hand, by making use of the Boussinesq approximation there is no longer a need for a six-component Reynolds-stress closure and the number of equations is reduced. The present study discusses how a robust and ecient two-equation solution method has been achieved and demonstrates the capabilities of the system using recent BAe AVRO research on a new regional jet.

Solution Method As its name suggests, a two-equation model requires the solution of two transport equations in addition to the Navier-Stokes equations. The turbulence equations are dierent from the mean ow equations in that they include important source terms which can cause the system to become numerically very sti. The approach adopted by BAe is to solve the turbulence equations explicitly in the same manner as the mean ow. The system is kept stable through the use of relaxation on the multigrid, limiters applied to certain quantities, and careful application of the arti cial viscosity. In RANSMB the dissipation takes the form of blended second and fourth dierences. For the mean ow the blending uses a shock sensor to switch on the second dierence at shocks. For the turbulence equations the second dierencing is invoked at discontinuities via a non-linear second dierence sensor on the turbulence quantities, and has a coecient large enough to give a rst order TVD property. An alternative method of solving such equations in such an explicit scheme is given by Davidson4, who proposes a semi-implicit method for the turbulence equations, which is coupled in some way to the explicit mean

ow solution. This method has proved successful, but due to its implicit nature becomes complex to implement, especially in a Multiblock environment.

The degree to which a transport model can represent key physical phenomena is also very important. For attached ows, most traditional models perform very well. However, these models often fail to predict more complex ows containing features such as boundary layer separation, either shock induced or pressure gradient induced5; 6 . Work is underway to analyse and implement some of the numerous corrections and modi cations which have been proposed to overcome such limitations. The present study is concerned with the choice of the basic model on which to proceed with such work.

All two-equation models use the Boussinesq approximation to model the turbulent stresses:   @ui @uj 2 @uk 0 0 ? ui uj = t ( @x + @x ) ? 3 @x ij ? 23 kij ; (1) Choice of Model j i k In order to install a transport turbulence model into an where the eddy viscosity is given by: industrial CFD code the developer faces a number of issues which may determine whether or not a certain k2 t = C  : (2) model is viable. In addition to the numerical stiness, " there are further features of the turbulence equations which are unique to each new model. When assessing a The widely accepted transport equation for turbulent particular turbulence model for use in RANSMB there kinetic energy, k, is modelled as:   are certain properties which are desirable in the context @ (ui k) of a structured, Multiblock method. These include: = @ ( + t ) @k + P ? "; (3) @xi

2

@xi

k @xi

k

where Pk is the production of k, and " is the speci c a one-equation model in the near-wall region, which is dissipation. matched to the two-equation model at a certain distance from the wall. The most commonly used near-wall model To achieve closure of the system an expression for the is that of Wolfshtein10, which solves Eqn.(3) in conjuncdissipation, ", is required. A large number of closures tion with an algebraic expression for ". have been proposed over the last 20-30 years. Most of This approach has been shown to be very eective6 . these give a transport equation for " itself, based on the A major drawback in the context of the present study `standard' form: is the need to nd a suitable location for the transition   @ (ui ") " @ t @" "2 = @x ( +  ) @x + C"1 k Pk ? C"2 k : (4) from the near-wall model to the two-equation model. In @xi i " i a structured grid system it is usual to select a particular parallel-to-wall grid line which represents a y+ value of An alternative closure method is provided by around 200. For three-dimensional ows with several Wilcox7. This is written in terms of ! which represents intersecting surfaces it becomes dicult to locate a grid surface which will be satisfactory. the ratio between " and k: " most appropriate method would appear to be the (5) use The !=  : k of a single model which is modi ed through the addition of wall damping terms to allow its use in both high The transport equation for ! is given as: and low Reynolds number regions. A large number of   ! such models exist in the literature. Many of these have @ t @! @ (ui !) = ( + ) +  Pk ? !2 : (6) been tested including Jones-Launder11 , Chien12, Lam@xi @xi ! @xi k Bremhorst13, Lien-Leschziner14 and Speziale9 . In practice it has been found that all these models give very Near-wall asymptotic analysis8; 9 shows that these similar solutions for most ow problems. They do, howbasic forms of k-" and k-! models give incorrect pre- ever, dier substantially with regard to the ease with dictions of key turbulent quantities approaching a solid which they can be made to produce stable solutions. surface. In the case of the k-" equations, it is necessary to add damping terms or use wall functions to reproduce Implementation of k-" experimentally observed behaviour. In the case of the k- One key issue encountered when applying low Reynolds ! model it has been found that despite the near-wall number k-" models is the boundary condition for ". inconsistencies, good predictions of boundary layer prok tends to zero as the wall is approached, " does les and skin friction can be obtained. The suitability of Whilst not vanish. The low Reynolds number models fall into a large number of both the k-" and k-! type of models two categories, those which apply a speci c boundary have been assessed for use within the BAe solver. condition for ", and those which replace " with new variables, such as "~ or  , which are de ned such that they To test the basic implementations, and to illustrate go to zero at the surface. the following discussions, a at plate boundary layer test case is used. This is calculated for a ow at M =0.5, One of the most popular of this latter category, and Re=7  106, using a grid with cell distributions similar one which has been in use at BAe for many years, is the to those typically used for aerofoil or wing calculations. Chien k-~" model. It is worthwhile repeating the de niThe cell nearest the wall is within y+  1 and there are tion of Chien's model here, to facilitate the description around 30 cells below y+ =1000. Transition is applied at of the implementation. x=l = 0:1 and the pro les are extracted at x=l = 0:9. The new dissipation is related to " by: 2k ; Assessment of k-" "~ = " ? (7) y2

As already mentioned, wall-eects are not suciently where y is the normal-to-wall distance. well modelled by the standard k-" equations, and these equations are only valid away from a surface, in high The new dissipation equation is: turbulent Reynolds number areas of the ow. The low   @  @ "~ @ (ui "~) Reynolds number region can be bridged with wall func= ( + t) + C"1 k"~ Pk @x @x  @x i i " i tions, which allow relatively coarse meshes to be used. However, no signi cant increase in capability over the "~2 2 "~ ?C4 y+ : (8) ? C"2 k f2 ? y2 e algebraic models is obtained. An alternative is to use 3

cell. This procedure is used for all blocks adjacent to walls. This requires that the grid be generated in such a way that blocks adjacent to walls are suciently large to contain the whole low Reynolds number region. For situations such as wing-body junctions where more than one surface is present, y is set to the smaller of the two + ? C y 3 f = 1 ? e : (10) wall distances. Whilst the boundary conditions for "~ are now straightforward the damping functions contain the wall-distance, y. In addition, the wall coordinate, y+ , is used. The de nition:y+ = u y= , contains the friction velocity, u , which is itself a function of velocity gradient at the wall. At separation points the resultant y+ can be zero. This / U+ is undesirable in the context of the damping functions, and an alternative de nition is used. This is derived using the assumption that Reynolds shear stress is proportional to kinetic energy and is written as15 : Here, f2 is a wall damping function, 2 f2 = 1 ? 0:22e?(Re =6) : (9) One additional damping function is used to factor the eddy viscosity in Eqn.(2), t

35

160

140

30

120

25

100

20

t

80

15

60

10

40

y

+

p

1 = C4 

ky : 

5

(11)

20

0

0 0

1

2

3

log(y+)

Comparing this model with the three criteria listed above it is clear that, although the boundary conditions appear to be straightforward (k=0, "~=0), the normal-towall distance measure, y, is used. For mesh blocks which do not have a wall boundary the damping functions are not applicable, and the model reverts to the `standard' high Reynolds number " model, Eqn.(4). This is a valid assumption, except immediately downstream of a trailing edge, where although there is no solid surface, the boundary layer(s) from the wing surface create a ow, which for a short distance, approximates a wall bounded

ow. Without a special treatment for wake regions the sudden transition from "~ to " here can lead to convergence diculties.

4

5

0

1

2

3

4

5

log(y+)

Figure 1: Log law comparison (left) and eddy viscosity pro le for the at plate boundary layer using the Chien(+Wolfshtein) model.

Towards the outer edge of the boundary layer, a new problem is introduced by the ratio of k2 =" which is present in the eddy viscosity, Eqn.(2). As the boundary layer edge is approached, both k and " fall to zero. Any slight error in the evaluation of these quantities may cause the ratio to become signi cant. This eect can cause large values of eddy viscosity to be predicted, even though the turbulent quantities themselves are small. To control the production of `spurious' eddy viscosity outside the boundary layer a minimum limiter for " is Further, it can be seen that Eqn.(8) contains source 16 introduced , terms which become very large in the near-wall region k 32 as y ! 0. Small imbalances in these terms can lead to : (12) "min = lmax large production or destruction of "~, with undesirable consequences for the convergence of the solution. This is derived from dimensional considerations, and To overcome the near-wall and wake problems a small the maximum length scale, lmax , is typically set to 0:1 layer of cells next to the surface is solved using the Wolf- aerofoil chord. This limiter improves stability and also shtein model. This approach resembles the method men- helps to reduce the production of spurious eddy viscosity. tioned above to bridge the gap between high Reynolds However, this eddy viscosity is not eliminated altogether. number models and the wall. However, since the Chien The at plate solution for the Chien-plus-Wolfshtein model is low Reynolds number compatible, the choice of grid line for separating the two regions is not critical. In model, together with the corresponding pro le of eddy RANSMB four to eight cells next to a wall are used for viscosity, t , is shown in Fig. 1. A satisfactory comparison with the log-law velocity pro le is obtained. Howthe one-equation model. ever, it can be seen that the eddy viscosity remains at a The distance, y, is found by calculating the distance level which is still much higher than the laminar viscosity along a normal grid line from the solid surface to a given outside the boundary layer. 4

Assessment of k-!

35

The k-! model appears to meet two of the three criteria immediately. The ! equation, Eqn.(6), can be integrated all the way to the wall without the requirement for damping functions or normal-to-wall distance, and the source terms appear to be straightforward.

1.2 Eqn. (14) use of R limiter

30

1.0

25 0.8 20

U+

The boundary conditions, however, pose a signi cant challenge. Since ! / "=k it tends to in nity as a wall is approached. This is entirely unsatisfactory. Various approximations and alternative boundary conditions have been tested. Even when using boundary conditions which appear entirely inappropriate, such as (@!=@n)wall = 0, good solutions can occasionally be obtained. However, such approaches are very dependent on the near-wall mesh spacing and are not appropriate for general use. The method proposed by Wilcox8 is to t ! to a function which is derived from near-wall analysis of Eqn.(6). Wilcox suggests that this function should be applied over seven to ten cells next to a surface, and that these cells should be in the region y+  2:5. For all practical applications considered in the present study this is too restrictive since there are typically only one or two cells inside this region.

0.6 15 0.4 10 0.2

5

0

0.0 0

1

2

3

log(y+)

4

5

0

1

2

3

4

5

log(y+)

Figure 2: Log-layer comparison (left) and  pro les for the k- model with and without the R limiter applied to the  -equation source terms.

local turbulence levels are negligible. The value of  in the freestream has been found to have a strong in uence on the turbulent quantities in the boundary layer. In Menter's study of the freestream dependency of the ! equation, it is suggested that in uence on the boundary is minimised provided that the freestream value of Analysis of Eqn.(6) highlights a further shortcoming layer ! is high of this model. In areas of the ow eld where the turbu- is required.enough. This implies that a small value of 1 lent energy, k, is negligible, the equation for ! de-couples from the k equation, and a non-zero solution for ! can The at-plate solution which results from using exist everywhere. This issue manifests itself as a sensi- Eqn.(14) in conjunction with a very small value of 1 tivity to freestream, or in nity, conditions for !. It has at the in ow boundary is shown in Fig. 2, together with been demonstrated by Menter17 that the level of eddy the corresponding  pro le. It can be seen that  conviscosity in a shear layer can be profoundly aected by tinues to grow outside the boundary layer, reaching a the speci ed value of !1 . maximum well away from the surface, where k is known to be negligible. High levels of  are now present in the freestream. Due to the convective and diusive terms Implementation of k-! these high levels of  are maintained further downstream. To overcome the boundary condition problem, the ! This is undesirable in situations where a new shear layer equation was inverted15 and expressed in terms of: is developing. One relevant example of this would be the case of a ap deployed behind a wing. For the boundary 1 (13) =  : layer on the ap, the large levels of  produced by the ! shear layers on the wing give an eective 1 which is The variable  now has units of time, and goes to zero much too high, and the growth of turbulence on the ap at a solid surface. The new transport equation for  is: is now incorrect.    @ (ui  ) @  @ To suppress these high levels of  , the source terms = ( + t) ?  k Pk in Eqn.(14) are modi ed by introducing a limiter based @xi @xi  @xi   on local eddy viscosity. For the k- model, the equation t 2 @ @ +   ?  +   @x @x : (14) for t is:  i i t =  k; (15) This model would now appear to ful ll all three of the re- and the limiter is de ned as: quirements listed above, although the freestream dependency is still present. As with the ! equation, a non-zero R = max(0:01;  k ): (16) solution of the  equation can be obtained even when the 5

The source terms in Eqn.(14) are multiplied by the ratio t =R and the transport equation for  is now written as:     2 @ (ui  ) @  @ = ( + t) ?  R Pk @xi @xi  @xi    @ @ : (17) +  Rt ?  + t 2 Rk @x @x 

i

approached,  approaches zero as a function of y2 . Since the method used to calculated local gradients uses rst order averages (as in most CFD methods), the quadratic nature of  is not properly accounted for. The last term in Eqn.(17) contains a product of the gradients of  , and it is the incorrect evaluation of this source term which leads to the poor modelling of the ow very close to the wall.

i

The limiter, R, is chosen such that in regions of the ow where the eddy viscosity is larger than 1% of the laminar viscosity the source terms are unaected. A useful sideeect of such a limiter is that the primary turbulence variables are now only used in the denominator when the local turbulence level is above the 1% cut-o. This precludes the possibility of division by zero.

To eliminate this discretisation error a higher order gradient calculation method could be employed near the surface. A simpler alternative is to modify the model to eliminate the quadratic behaviour. A further transformation is carried out by the rst two authors to re-cast the  equation in terms of a new variable, g, which represents p . Using the same source term limiter, R, the The solution for this modi ed model on the at plate new model can be written as: is also given in Fig. 2. Once the edge of the boundary   layer is reached, the eddy viscosity decreases rapidly, and @ (ui k) = @ ( + t ) @k + Pk ?  2 k2 ; (18) @xi @xi k @xi R the limiter is activated, causing  to be reduced smoothly to a very low level in the freestream. Meanwhile, the    g3 velocity pro le is seen to remain largely unaected. @ (ui g) @ t @g = @x ( +  ) @x ?  2R Pk @xi Care is required in selecting the appropriate value  i   3g  kgi @g @g 2  kg for the limiter. Very close the wall the eddy viscosity + 2R ?  +  t : (19) R @xi @xi is again small. All the calculations presented here use g meshes where the minimum wall cell height has a y+ value between 0.5 and 2. For these meshes the limit of The eddy viscosity is: 0:01 is found to work well. For ner meshes, a lower t =  kg2 ; (20) limit would be required. and the constants are: 35   C = 0:09; g    ! = 2:0; 30 k = 5=9; = 0:075 k-g 25 The new solution is shown in Fig. 3 and it can be seen 20 that the correct velocity pro le is now obtained. U+

A similar transformation has been carried out by Gibson and Dafa'Allap18 to convert a k-" model to q10  , where q represents k. However, this transformation adds new source terms to the k equation. Since the k 5 diusion term is the only term involving gradients of k in the k-g model the discretisation errors at the wall have 0 a much reduced eect. Hence, a transformation of the k -1 0 1 2 3 4 5 equation does not appear to be necessary in the context log(y+) Figure 3: Flat plate velocity pro les for the k- of the present study. model, Eqn.(18) and the k-g model, Eqn.(19). 15

2D Validation

Concerning the velocity pro le predicted by the k model, it is apparent in Fig. 2 that the model gives an incorrect prediction of U + implying that there is an over-prediction of skin friction. This can been traced to an error in the calculation of the gradient of  very close to the solid surface. It is known that as the wall is

Both the k-~" and k-g implementations described above have been validated on a large number of external ows for both 2D and 3D geometries. To demonstrate the capabilities of both models a smaller number of test cases will be described. First, a transonic aerofoil, and a low 6

4.5 0.008 -1.0 0.006 -0.5

4.0

0.004

Cp

Lift Coefficient

3.5

Cf 0.002

0.0

0.0 0.5

Experiment kk-g

1.0

0.0

0.2

0.4

0.6

X/c

0.8

2.5

-0.002

-0.004 1.0

3.0

Experiment kk-g

2.0 0.0

0.2

0.4

0.6

0.8

1.0

X/c

1.5

Figure 4: Pressure (left) and friction coecients for RAE2822 Case 9.

0

5

10

15

Incidence

20

25 -0.05

0.0

0.05

0.1

0.15

Drag coefficient

Figure 6: Lift/Incidence and Lift/Drag polars for the L1T2 at M = 0:197, Re = 3:52  106 .

Mach number, high-lift, multi-element aerofoil are presented and discussed. Following this, some examples of To test the models for ows with high pressure grathree-dimensional applications related to current BAe dients, and greater geometrical complexity, a detailled AVRO research work are given. study has been carried out for a three-element high-lift A well known, and geometrically simple test case is con guration. The geometry was de ned as part of the the RAE2822 aerofoil. Test ow `Case 9' is compared National High-Lift Programme (NHLP) and is known as with experimental data19 for surface pressure and fric- `L1T2'. Experiments carried out on this geometry are tion in Fig. 4. The corrected conditions are M = 0:734, described by Moir20. Calculations have been performed  = 2:54o and Re = 6:5  106 . It can be seen that both with both k-~" and k-g over the complete incidence range models give very similar predictions, which correspond at M = 0:197, Re = 3:52  106. The geometry and well with the experimental data. In particular, the fric- streamlines for a representative ow are shown in Fig. 5. tion coecient plot demonstrates that the k-g model is Of particular interest are the interactions between the performing well in the near-wall region. wake from the slat and the boundary layer on the upper surface of the wing, and the wing wake and ap boundary layer. There are also complex re-circulation regions in the `cove' areas which exist on the lower surface of the slat, and in the ap cut-out on the lower surface of the wing. Polar plots of lift against drag and lift against incidence are given in Fig. 6. The predictions of lift from the k-g model correspond well with experiment, while those of the k-~" implementation are slightly lower. The lift/drag curve shows similar dierences. The k-g solution is closer to experiment while the Chien model overpredicts drag at every point. A closer examination of the solutions at 20:18o shows that there is little dierence in the predictions of surface pressure coecient, shown in Fig. 7. The lower lift predicted by k-~" results Figure 5: Streamlines over the L1T2 three-element from a slightly lower suction over the upper surface of aerofoil at 21.69o predicted by the k-g model. the wing and slat. It is well documented8 that k-! tends

7

shortcoming of the Chien implementation. A certain amount of `spurious' eddy viscosity outside the boundary layers can be seen, which aect a large area of the ow above the geometry. As observed in the at plate test case, this value of eddy viscosity is controlled, but not eliminated, by the use of limiters. The presence of such eddy viscosity can trigger undesirable diusion or even production of k which in turn can aect the levels of turbulence at the boundary layer edge further downstream. Meanwhile, the k-g implementation behaves well at the boundary layer edge, and the turbulence is seen to return to negligible levels smoothly.

-14

Experiment kk-g

-12 -10 -8 Cp -6 -4 -2 0 0.0

0.2

0.4

0.6 X/c

0.8

1.0

1.2 8

Figure 7: Surface pressure coecient for L1T2 at M = 0:197,  = 20:18o , Re = 3:52  106 .

8

Wing,35%

Flap, l.e.

7

7

Exp.

6

to give better predictions for adverse pressure gradient

ows. It would appear that this behaviour is preserved in the present k-g formulation. Pro les of total pressure at the same ow conditions are also available and are given in Fig. 8. These are presented for the ow above the wing at 35% chord, and at the leading edge, mid chord and trailing edge of the ap. Here, larger dierences between the models can be seen. The k-~" solution appears to give incorrect spreading of the wakes, with a rather more diuse total pressure distribution. The k-g model is consistently closer to experiment.

6

Inches above surface

k-ε k-g

5

5

4

4

3

3

2

2

1

1

0

0 -0.5

0.0

0.5

1.0

-0.5

0.0

0.5

1.0

For an explanation of these dierences, it is useful to examine the way in which the turbulent ow is hanFlap, 50% Flap, t.e. dled by these two models. Fig. 9 shows contours of noro malised eddy viscosity for the 21:18 case. Immediately evident is the large number of block boundaries which are present, this being a consequence of the mesh generation procedure which demands a single boundary condition type on a block face. Some of the block boundaries are however utilised for grid control by means of user speci ed line-grid stretching functions. The block boundaries, when they are positioned close to shear layers, can cause problems for models which require wall distance measures. For many test cases it is possible to ensure that all low Reynolds number regions are contained within the block nearest the wall. However, for moderately complex geometries, such as the L1T2, this may not be feasible. As previously stated, the k-~" model implementation reverts to the high Reynolds number form once the opposite face of a wall block is reached. It is likely for that, at such high incidences, the model is switching too Figure 8: Pro les of total o pressure coecient 6. L1T2 at M = 0 : 197 ,  = 20 : 18 , Re = 3 : 52  10 early. This aects the region above the ap in particular, where the wake and boundary layer mixing regions are still at low turbulent Reynolds number conditions. The contours of normalised eddy viscosity highlight another Inches above surface

Cpo

Cpo

8

8

7

7

6

6

5

5

4

4

3

3

2

2

1

1

0

0

-0.5

0.0

0.5

Cpo

8

1.0

-0.5

0.0

0.5

Cpo

1.0

number, or the addition of extra components such as winglets and ap track fairings. Also required is the ability to reliably detect the buet boundary. Primary areas of interest in Navier-Stokes results will be associated with shock locations and degradation of the boundary layer towards the trailing edge including 3D eects, separations, and trailing edge pressure rise. Related design areas where CFD has a particular contribution to make are in the details of ow in wing-fuselage and pylon-wing junctions, these being very sensitive to Reynolds number dependent separations, often in the presence of quite strong shocks in the wing-nacelle gap. Thus validation of the code for representative cases is vital. Figs. 10 and 11 show results for a BAe AVRO regional jet wing-body research con guration for which tunnel results are available. Flow conditions are transonic cruise at M = 0:79,  = 1:0o, transition at 15% chord, CL  0:42, but at a tunnel Reynolds number based on c of 2.7 million. The mesh has 1.2 million cells in 230 blocks with a rst cell depth of y+  2 and around 25 cells in the boundary layer. Pressure contours in Fig. 10 are obtained by use of the k-g turbulence model. Fig. 11 shows surface streamlines using both k-g and k-~" models. These streamlines can be compared to wind tunnel oil- ow visualisations. The k-~" results show a separation at the wing root trailing edge which is not present in the wind tunnel. This separation is not shown by the k-g results, for which the streamlines match the behaviour seen in the wind-tunnel. Note that the wing is attached to a belly fairing which presents a vertical wall at the wing root thus avoiding acute angles at the wing-fuselage junction, and that the leading edge root fairing further suppresses any tendency to separation. Elsewhere the two sets of streamlines show very similar properties, notably the curvature and incipient separation near the trailing edge crank which is also suggested in the tunnel oil ow. Although not clear in these gures, the surface ow pattern around the leading edge fairing is well predicted by both models. It is signi cant that trailing edge root separations similar to that seen with k-~" are also produced for this case when Baldwin-Lomax and Johnson-King turbulence models are used in RANSMB, the former model being noted for its resistance to separation. The one common feature in these three models that is not shared with k-g is the need for wall distance, which may be one source of error. Fig. 12 shows pressure contours for a dierent wing/body plus pylon/nacelle research con guration using the k-g model. Conditions are M = 0:8,  = 1:0o ,

Figure 9: Contours of eddy viscosity (t =) for the Chien model (top) and the k-g model. Values range from 0 (Red) to 800 (Violet).

3D Applications Two applications are described from current BAe AVRO research studies on regional jet con gurations. At present, the primary requirement for this ow solver is transonic cruise conditions with very accurate drag predictions for a wide spread of lift coecients around the design point extending well into the wave drag rise. The wing designer would routinely compute CL/CD polars over a CL range of 0.1 to 0.7. Initial accuracy of 1% of total aircraft drag is sought together with drag increments to substantially less than one count. These increments may be due to variation of shape or Reynolds 9

Figure 10: Contours of pressure coecient on a typical regional jet transport wing/body predicted by the model. M = 0:79,  = 1:0o , Rec = 2:7  106 .

k-g

Figure 11: Skin friction lines on the same wing/body for k-"~ (left), and k-g.

10

Figure 12: Pressure coecient contours for a wing/body/pylon/nacelle research con guration calculated using the k-g model. M = 0:8,  = 1:0o , Re = 6:0  106 .

11

transition again at 15% chord, Rec = 6  106, and CL  0:41. The mesh has 2 million cells in 900 blocks with a rst cell depth of y+ = 2. This is part of an ongoing study into nacelle-pylon installation and is shown here with a very closely coupled nacelle and overwing pylon. Flow disturbance can be seen at the wingfuselage junction. A small horseshoe vortex appears at the leading edge root. No belly or wing root fairings are present in this simulation and the wing-fuselage junction presents an acute angle at the trailing edge. This could be expected to lead to incipient separation in the wing-fuselage junction. Tunnel tests are planned.

while both models give robust convergence on moderately complex geometries, the solutions given by the k-g model are at least as good, and often better than those of the k-~" model. When applied to three-dimensional design tests, stable solutions can be obtained with both transport models. The k-g model, with its lack of dependency on wall distance, appears particularly useful for complex aircraft con gurations. Initial results indicate that accurate and reliable solutions can be obtained. The increased con dence in the solutions which are now being obtained is helping to increase the use of CFD by the design teams.

Computational Cost

Acknowledgement

The wing-body solutions required 23 Mwords and 10 The rst author's stay at BAe Sowerby was sponsored by hours on a Cray YMP machine at 60 M op/s as part the European Commission under speci c training agreeof a CL/CD polar calculation. The wing-body-pylon- ment No. AER2-CT93-5006H. nacelle solution required 41 Mwords and 75 hours on a Cray J90 machine at 32 M op/s from a free-stream start. References Although wing-body CL /CD polars are routinely computed on the YMP, work is underway to implement the

ow solver on the Cray T3D massively parallel proces- [1] Jameson A., Schmidt W. & Turkel E. Numerical solutions of the Euler equations by nite volume sor to facilitate this and extend the capability to nacelle methods using Runge-Kutta time-stepping schemes. case polars. AIAA Paper 81-1259, 1981. [2] Jameson A. & Baker T. Multi-grid solutions of the Conclusions Euler equations for aircraft con gurations. AIAA Paper 84-0093, 1984. The way in which particular two-equation transport models are chosen for implementation in an industrial [3] Baldwin B.S. & Lomax H. Thin Layer Approximadesign system has been described. This process is govtion and Algebraic Model for Separated Turbulent erned almost entirely by numerical considerations. One Flows. AIAA Paper 78-257, 1978. k-" model and one k-! model have been selected and the properties of these two models are measured against [4] Davidson L. Implementation of a Semi-Implicit three criteria: no normal-to-wall distance, simple source k- Turbulence Model Into an Explicit Rungeterms and straightforward boundary conditions. Both Kutta Navier-Stokes Code. CERFACS Report, models require manipulation before a stable implemenTR/RF/90/25, 1990. tation is achieved. The Chien k-~" model has simple boundary conditions, but makes use of wall distance and [5] Gould A.R.B. Validation of Turbulence Models for Two- and Three-dimensional Flows. In: ECARP also requires damping functions which become quite sti. (Validation Area) Report, to be published in Notes These issues are overcome using a small one-equation on Numerical Fluid Mechanics Series, by Vieweg, sub-layer next to wall surfaces, and by making use of 1996. a limiter for "~. Despite these precautions the solutions obtained are often adversely aected by the numerical [6] Lien F.S. & Leschziner M.A. Modelling 2D Sepdiculties encountered. The Wilcox k-! model has no aration from High-Lift Aerofoil with Non-Linear damping functions and simpler source terms, but the Eddy-Viscosity Models and Second-Moment Cloboundary condition for ! makes it dicult to use. This sure. UMIST Report, TFD/94/05, 1994. model is rewritten to use a new variable, g, and all three criteria can now be satis ed. The dependency on [7] Wilcox D.C. Reassessment of the scale determinfreestream conditions is controlled through the careful ing equation for advanced turbulence models. AIAA use of a source term limiter which has no eect in boundJournal, 26 pp1311-1320, 1988. ary layers and wakes. Results of two-dimensional transonic and low Mach number test cases demonstrate that, [8] Wilcox D.C. Turbulence modelling for CFD. DCW Industries Inc., 1993. 12

[9] Speziale C.G., Abid, R. & Anderson E.C. Criti-

cal Evaluation of two-equation models for near-wall turbulence. AIAA Journal, 30 pp324-331, 1992. [10] Wolfshtein M.W. The Velocity and Temperature Distribution in One-Dimensional Flow with Turbulence Augmentation and Pressure Gradient. Int.

Journal of Heat and Mass Transfer, 12 pp301-312, 1969. [11] Jones W.P. & Launder B.E. The prediction of laminarisation with a two-equation model of turbulence.

Int. Journal of Heat and Mass Transfer, 15 pp301314, 1972. [12] Chien K-Y. Predictions of Channel and Boundary-

Layer Flows with a Low-Reynolds-Number Turbulence Model. AIAA Journal, 20 pp33-38, 1982. [13] Lam C.K.G. & Bremhost K.A. Modi ed Form of k- Model for Predicting Wall Turbulence. ASME

Journal of Fluids Engineering, 103 456-460, 1981. [14] Lien F.S. & Leschziner M.A. Computational Mod-

elling of 3D Turbulent Flow in S-Diuser and Transition Ducts. Engineering Turbulence Modelling

[15]

[16]

[17] [18] [19]

and Measurements, 2 Elsevier, 1993. Kalitzin G. Validation and Development of TwoEquation Turbulence Models. In: ECARP (Validation Area) Report, to be published in Notes on Numerical Fluid Mechanics Series, by Vieweg, 1996. Benton J.J. Description of Methods Used by British Aerospace. In: `EUROVAL, A European Initiative on Validation of CFD Codes.' Haase, Brandsma, Elsholz, Leschziner, Schwamborn (eds.) Notes on Numerical Fluid Mechanics Vol 42, Vieweg, 1993. Menter F.R. In uence of Freestream Values on k-! Turbulence Model Predictions. AIAA Journal, 30 pp1657-1659, 1992. Gibson M.M. and Dafa'Alla A.A. A Two-Equation Model for Turbulent Wall Flow. AIAA Journal, 33 pp1514-1518, 1995. Cook P.H., MacDonald M.A. & Firmin M.C.P. Aerofoil 2822 - Pressure Distributions, Boundary Layer and Wake Measurements. Appendix A6,

AGARD AR-138 , 1979. [20] Moir I.R.M. Measurements on a Two-Dimensional Aerofoil With High Lift Devices. Appendix A2, AGARD AR-303 , 1994.

13