EuropeanCongressonComputationalMethodsinAppliedSciencesandEngineering ECCOMAS2000 Barcelona,11-14September2000 (cid:13)cECCOMAS LESFOIL: A EUROPEAN PROJECT ON LARGE EDDY SIMULATIONS AROUND A HIGH-LIFT AIRFOIL AT HIGH REYNOLDS NUMBER L. Davidson Dept. of Thermo and Fluid Dynamics, Chalmers University of Technology SE-412 96 G(cid:127)oteborg, Sweden, Email: [email protected] web page: www.tfd.chalmers.se/~lada Key words: LES, wall functions, airfoil, channel (cid:13)ow, LES-RANS, DES. Abstract. This paper presents work carried out in the LESFOIL project, which studies Large Eddy Simulations (LES) of the (cid:13)ow around high-lift airfoils, during its (cid:12)rst 24 6 months. The Reynolds number for the selected Aerospatiale A-airfoil is high (Re = 2:1(cid:1)10 based on the freestream velocity and the chord length). If some kind of near-wall treatment could be used, the near-wall streaks would not be resolved and a much coarser grid could be used in the streamwise and spanwise directions. Two di(cid:11)erent near-wall treatment methodologies are used in the LESFOIL project, either hybrid LES-RANS (or DES) or o wall functions. The angle of incidence is (cid:11) = 13:3 , and a small separation bubble is, according to experiments, present on the suction side in the trailing edge region. Thus the method for treating the near-wall wall region must be able to handle both attached boundary layer (cid:13)ow, including streamwise pressure gradients, and separated (cid:13)ow. Subgrid-scale (SGS) models and parallelized numerical methods are two other subjects covered in the LESFOIL project. SGS models are developed and evaluated in simple (cid:13)ows such as channel (cid:13)ow. The hill (cid:13)ow (Reynolds number 10;000 based on the hill height) is also used as a test case in which the performance of the wall treatment and SGS models can be evaluated in recirculating (cid:13)ow. The (cid:12)rst part of the paper presents LES of channel (cid:13)ow. The second part presents hill (cid:13)ow computations. In the following section, LES around the A-airfoil is shown and conclusion are drawn in the (cid:12)nal section. 1 L. Davidson 1 INTRODUCTION Computational Fluid Dynamics (CFD) is an important tool in the design of aircraft. For conventional aircraft under cruise conditions, the current generation of CFD codes allows a computer-based design to be achieved with some con(cid:12)dence. The design process of future aircraft requires extensive research and development in CFD. E(cid:14)cient and ac- curate CFD methods allow the aircraft industry to reduce extensive experimental testing, which usually is very expensive. CFD will also allow faster development of new aircraft and modi(cid:12)cation of existing aircraft, as it may be su(cid:14)cient to use measurements only in the later stages of development. To perform accurate numerical simulations of such complex (cid:13)ows, new and more ad- vanced turbulence models are needed. In the LESFOIL project, the feasibility of the application of Large Eddy Simulations (LES) to aerodynamic (cid:13)ows will be assessed, both with aspect to the reliability of the results and the computer resources needed for the calculations. The main objectives of the project are: (cid:15) to provide know-how of an advanced CFD method to the European aeronautical industry; (cid:15) to determine when LES will be feasible for the European aeronautical industry; (cid:15) to demonstate the feasibility of LES of (cid:13)ow over an airfoil; (cid:15) to assess the computational requirements to carrying out LES for the simple airfoil and for more complex con(cid:12)gurations in future; (cid:15) to develop highly e(cid:14)cient numerical methods for the LES of airfoil (cid:13)ows; (cid:15) to compare competing subgrid-scale models for these (cid:13)ows. The work is divided into (cid:12)ve main tasks. Task 1: Development of subgrid models. If reliable subgrid models are developed, it would be possible to model a larger portion of the turbulence, allowing coarser grids. Task 2: Near-walltreatment. AthighReynoldsnumbers, LESrequiresa(cid:12)negridinallthree coordinate directions. Use of wall functions or simple RANS (Reynolds Averaging Navier-Stokes) turbulence models in the near-wall region will be studied. This wouldallowustouse coarsergridspacing, especiallyinthe streamwiseand spanwise directions. Task 3: Transition. In airfoil (cid:13)ow, transition occurs on both the suction and pressure sides. For a correct development of the boundary layer on the suction side, it is crucial to trip transition in a correct way. 2 L. Davidson Task 4: Development of e(cid:14)cient numerical parallel methods. LES is computationally very expensive, and it is thus extremely important to reduce the computation time as much as possible. Task 5: LES of airfoil (cid:13)ow. In this task the experience gained in Tasks 1{4 is put together. Nine parters work in the LESFOIL project, which is coordinated by Chalmers. The partners are (with the researcher in charge in parentheses): 1. Chalmers, Sweden (L. Davidson) 2. Alenia, Italy (N. Ceresola) 3. CERFACS, France (F. Ducros) 4. Dassault, France (M. Mallet) 5. Fluent Europe Ltd, UK (D. Cokljat) 6. Karlsruhe University, Germany (W. Rodi) 7. ONERA, France (T.H. Le) 8. University of Surrey, UK (P. Voke) 9. UMIST/QMW, UK (M.A. Leschziner) The LESFOIL project started on the 1st of February 1997. This paper presents the outcome of some of the cooperative work carried out thus far in the LESFOIL project. Most of the (cid:12)gures and discussions below have been taken from progress reports written by the partners in the LESFOIL project. 2 CHANNEL FLOW Some work in the LESFOIL project has been carried out studying channel (cid:13)ow. Simu- lationsofthis(cid:13)owarefairlycheapandcanbecarriedoutonworkstationsorPCs. Because the required CPU time is reasonable, parametric studies have been done in the project investigating di(cid:11)erent grids, subgrid modelsand wallmodels. The experience gained from this simple (cid:13)ow is expected to be valuable in airfoil simulations. Although the channel (cid:13)ow is geometrically very simple, it is a wall-bounded (cid:13)ow in which the near-wall tur- bulence must be properly predicted in order to perform accurate simulations. In airfoil (cid:13)ow, the attached (cid:13)ow along the suction side prior to separation shares many features of channel (cid:13)ow, although the boundary layer in the former case experiences both favourable and adverse pressure gradients, which are absent in channel (cid:13)ow. 3 L. Davidson 2.1 Wall-resolved Computations Coarse grid simulations can be considered meaningless since large turbulent structures in the near wall layer cannot be captured. Nevertheless the resolution is representative of the resolutionused insomeregionsofthe airfoil(cid:13)ow calculation. It isthus worth checking how subgrid scale models behave on such a poor grid, which, in addition, contains high aspect ratio cells. CERFACS[1]hasusedseveralgridsinsimulationswiththelocalizeddynamicmodel[2]. Pro(cid:12)les of mean velocity and rms (cid:13)uctuations shown in Fig. 1 show the improvement in the results when the grid is re(cid:12)ned. This is rather reassuring but also proves that the subgrid scale model does not completely ful(cid:12)l its duty: one could expect results to be insensitive to the grid resolution, the subgrid scale model taking on a greater portion of the job as the grid iscoarsened. This isnot the case for the range of resolutionconsidered. Actually, the opposite occurs, as shown in Fig. 2. The subgrid scale dissipation provided by the model is actually higher in the (cid:12)ner grid simulation. Eddy-viscosity pro(cid:12)les are interesting: in the region where grids are certainly too coarse to resolve the large eddies, + (cid:23)t is stronger on the (cid:12)ner grid but, for y > 50, as the grids become more adequate, the expected behaviour is recovered, (cid:23)t being stronger on the coarser grid. 4.0 LES Piomelli D x+=412 D z+=154 3.0 D x+=412 D z+=102 20.1 D x+=412 D z+=51 s D x+=206 D z+=51 m 2.0 r + u U 10.1 u+ =2.5 ln y+ +5.2 1.0 0.1 0.0 1 10 100 1000 0.0 0.2 0.4 0.6 0.8 1.0 + y y=(cid:14) 1.5 1.5 1.0 1.0 ms ms vr wr 0.5 0.5 0.0 0.0 0.0 0.5 1.0 0.0 0.5 1.0 y=(cid:14) y=(cid:14) Figure 1: Re(cid:28) = 1050. Pro(cid:12)les of streamwise velocity and of root-mean-square velocity (cid:13)uctuations in + simulations with the localized dynamic model. Grid spacings in the wall normal direction are (cid:1)y = 2(cid:0)86. Results from Ref. [1]. 4 L. Davidson 101(cid:23)t=(cid:23) 102 100 101 10−1 n o 10−2 ati 100 p 10−3 i ss D x+=412 D z+=154 10−4 DDxx++==441122 DDzz++==115042 di10−1 DD xx++==441122 DD zz++==15012 10−5 DDDxxx+++===422105260 DDDzzz+++===555111 GS DD xx++==220560 DD zz++==5511 S 10−6 1 10 100 1000 1 10 100 1000 10000 + + y y Figure 2: Re(cid:28) = 1050. Localized dynamic model. Left: Pro(cid:12)les of time-averaged eddy-viscosity. Right: Pro(cid:12)les of time-averaged modeled subgrid scale dissipation. Results from Ref. [1]. Figure 3 shows Chalmers’ predictions. The dynamic one-equation model by David- son [3] (see also Ref. [4]) was employed, in which an equation is solved for the subgrid kinetic energy (see Eq. 2). The coe(cid:14)cients in the production (C) and dissipation (C(cid:3)) terms are computed dynamically. To ensure numerical stability, a constant value of C in space (Chom) is used in the momentum equations. The idea is to include alllocaldynamic informationthrough the source terms of the transport equation for ksgs. In this way, large (cid:13)uctuations in the dynamic coe(cid:14)cients are smoothed out in a natural way and the need to restrict or limit the dynamic coe(cid:14)cients is eliminated altogether. The sensitivitytothe gridpresented inFig.3exhibits thesametrend as thepredictions made with the Dynamic SGS model in Fig. 1. 25 4 s3.5 20 m r 3 w 15 ;2.5 s + rm 2 U 10 v 1.5 ; s m 1 5 r u0.5 0 100 102 00 0.2 0.4 0.6 0.8 1 + y y=(cid:14) Figure 3: Channel (cid:13)ow. Re(cid:28) = 1050. Velocity and resolved normal stresses using the dynamic one- equation model [3]. (cid:1)x+ =200, (cid:1)y+ =3:9(cid:0)85. Thick lines: (cid:1)z+ =50; thin lines: (cid:1)z+ =25. Results from Ref. [5]. 5 L. Davidson 2.2 Wall functions 6 The Reynolds number for the selected Aerospatiale A-airfoil is high (Re = 2:1 (cid:1) 10 based on the free stream velocity and the chord length). A wall-resolved mesh based + + + on (cid:1)x ’ 100 (streamwise direction), (cid:1)ymin ’ 1 (wall-normal direction) and (cid:1)z ’ 20 would require a mesh of approximately 2000 (cid:2) 100 (cid:2) 300 = 60 million cells, which is unrealistic. If instead some kind of near wall treatment can be used, in which the near-wall streaks would not be resolved but modelled altogether, a much coarser grid could be used in the streamwise and spanwise directions. By increasing the grid size, time step constraints are also considerably alleviated. Two di(cid:11)erent near-wall treatment methodologies are used in the LESFOIL project, either a hybrid LES-RANS (or DES) or + wall functions. For both methodologies, the required cell spacing would be 100 < (cid:1)x < + 600 and 100 < (cid:1)z < 300. The cell spacing in the wall-normal direction would be much + + the same as for RANS computations, i.e. (cid:1)y ’ 1 or (cid:1)y ’ 30 (near-wall node) for a near-wall RANS method and wall functions, respectively. In the airfoil (cid:13)ow chosen in the LESFOIL project, a small separation bubble is, ac- cording to experiments, present on the suction side in the trailing edge region. Thus the method for treating the near-wall wall region must be able to handle both attached boundary layer (cid:13)ow including streamwise pressure gradients as well as separated (cid:13)ow. While developing such a method is a formidable task that cannot be expected to be accomplished in the LESFOIL project, a (cid:12)rst attempt will be made. In the LESFOIL project, two di(cid:11)erent wall functions, the law of the wall and the Werner-Wengle wall function, are used. Both wall functions are formulated in instanta- neous quantities. The Werner-Wengle wall function uses a power (cid:12)t of the logarithmic law and allows an explicit evaluation of the friction velocity. Use of a wall function allows + placement of the (cid:12)rst point well above the wall (y ’ 30). Karlsruhe [6] evaluated the performance of the dynamic SGS model on stretched and distorted grids. The Werner-Wengle wall function [7] provides the near-wall treatment and the dynamic model the SGS stresses. DNS of Moser et al. [8] is used as reference data. A 222 (cid:2) 38 (cid:2) 66 node grid is used and is 8(cid:25) (cid:2) 2 (cid:2) (cid:25) units in size. Periodic boundary conditions are applied in the streamwise and spanwise directions. Wall normal + + and spanwise cell sizes are constant at (cid:1)y = 30 and (cid:1)z = 60, respectively. In the streamwise direction, the grid is divided into (cid:12)ve di(cid:11)erent regions. Regions AB and EF + contain cells of streamwise size (cid:1)x = 30. The cells in region CD are also of uniform size + with (cid:1)x = 220. The cells in regions BC and DE are either stretched or compressed in the streamwise direction at a rate of 5%. Regions AB and EF are (cid:25) long, while regions BC, CD and DE are each approximately 2(cid:25) in length. Figure4showssomeresultsfromthecomputationundertakenwiththe5%stretched/com- pressed grid. It can be seen that the grid stretching/compressing has an impact on the time-averaged statistical data. The mean velocity pro(cid:12)le appears to be least a(cid:11)ected. Variation at the channel centerline is in the order of 1%, increasing to approximately 6 L. Davidson + 3% at y = 63. The di(cid:11)erences between this computation and the computation with a + constant (cid:1)x = 220 are of similar magnitudes. The increase and decrease of the stream- wise stresses when increasing/decreasing (cid:1)x agree with the results in Figs. 1 and 5. A large/small (cid:1)x gives large/small streamwise stresses. It should be noted that the evolu- tion of the wall-normaland spanwise stresses is the reverse [6]. This is also seen in Figs. 1 and 5. The results presented in Fig. 4 show that the stretching in the streamwise direction of a grid in a boundary layer should be kept at a minimum. This (cid:12)nding is a con(cid:12)rmation of earlier (cid:12)ndings of Surrey [9] when doing LES of the (cid:13)ow around turbine blades. The results obtained above should be kept in mind when generating the grid for the airfoil. 21 20 y+=590 15 1290 B:stretch C:constant D:shrink E:constant +U> +U> 18 <10 5%stretching < 5%stretching LDLooNgcSa-ltaiownB 17 dyy++x+===56293200 5 LocationC LocationD 16 LocationE y+=63 15 0 100 101 102 5 10 15 20 25 y+ x 10 12 y+=63 11 9 78 5%stretching 1890 B:stretch C:constant D:shrink E:constant 6 DNS 7 +> LocationB +> uu 5 LocationC uu 6 < LocationD < 5%stretching 4 LocationE 5 dx+=220 4 y+=590 3 y+=63 3 2 2 1 1 y+=590 0 0 0 100 200 300 400 500 5 10 15 20 25 y+ x Figure 4: Channel (cid:13)ow, Re(cid:28) = 590. Velocity pro(cid:12)le and streamwise normal stresses. AB and EF: + + + (cid:1)x =30; CD: (cid:1)x =220; BC and DE: (cid:1)x stretched/compressed by 5%. Solid lines in the (cid:12)gures + to the right show predictions with constant (cid:1)x =220. Results from Ref. [6]. 7 L. Davidson Some results obtained by CERFACS [10] are shown with the Werner-Wengle wallfunc- tion in Fig. 5. The WALE subgrid-scale model[11] is used and several resolutions are con- sidered. Overall, results are satisfactory, even on the coarsest grid. Of course, the peaks of velocity (cid:13)uctuations cannot be properly captured, but the simulations give acceptable estimations of their amplitude. CERFACS [10] has compared, for one case (20(cid:2)20(cid:2)20 grid), the results obtained with a wall function based on the logarithmic law of the wall with those obtained with the Werner-Wengle wall function. Results are similar but not identical; in particular, the predicted turbulent Reynolds number shows a 6% di(cid:11)erence between the two simulations. This is however within the range of the approximation to the logarithmic law by the power law in the Werner-Wengle wall function. 4.5 4.0 LES Piomelli 20.1 3.5 50x20x50 40x20x40 3.0 20x20x20 15.1 2.5 20x10x10 +U u’ 2.0 10.1 5.1 u+=2.5 ln y++5.2 11..05 0.5 0.1 0.0 1.0 10.0 100.0 1000.0 0.0 0.5 1.0 y+ y/d 1.5 1.5 1.0 1.0 v’ w’ 0.5 0.5 0.0 0.0 0.0 0.5 1.0 0.0 0.5 1.0 y/d y/d Figure 5: Re(cid:28) = 1050. Pro(cid:12)les of streamwise velocity and of root-mean-square velocity (cid:13)uctuations in + simulations with the Werner-Wengle wall function and the WALE SGS model. Resolutions are: (cid:1)x = 164,(cid:1)y+ =105,(cid:1)z+ =32(50(cid:2)20(cid:2)50);(cid:1)x+ =206,(cid:1)y+ =105,(cid:1)z+ =41(40(cid:2)20(cid:2)40);(cid:1)x+ =412, (cid:1)y+ =105, (cid:1)z+ =82 (20(cid:2)20(cid:2)20); (cid:1)x+ =412, (cid:1)y+ =210; (cid:1)z+ =164(20(cid:2)10(cid:2)10). Results from Ref. [10]. The work carried out by QMW [12] on channel (cid:13)ow with wall functions con(cid:12)rm the results obtained by CERFACS. Figure 6 compares the velocity pro(cid:12)le obtained with the Werner-Wengle wall function with that obtained with the instantaneous logarithmic law. The predictions in Fig.6 should be compared with the predictions of CERFACS in Fig. 5. 2.3 LES-RANS As mentioned above, wall-resolved LES is not a(cid:11)ordable in airfoil (cid:13)ow. Using wall functions is one alternative and hybrid LES-RANS is another. The method presented in 8 L. Davidson 25 25 20 DNS 20 DNS WW + Wale LLI + WALE 15 15 10 10 5 5 0 0 1 10 100 1000 1 10 100 1000 + + y y Figure6: Velocitypro(cid:12)lesusingtheWerner-Wenglewallfunction(left)andtheinstantaneouslogarithmic law (right). Re(cid:28) =1050. The WALE SGS model [11] is used. Results from Ref. [12]. [13] combines the two-equation k (cid:0)! model of Peng et al. [14] @k @ @ (cid:23)t @k + (Ujk) = (cid:23) + +Pk (cid:0)ckfk!k @t @xj @xj (cid:20)(cid:18) (cid:27)k(cid:19) @xj(cid:21) @! @ (cid:23)t @! ! (cid:23)t @k @! + (Uj!) = (cid:23) + + (c!1f!Pk (cid:0)c!2k!)+c! (1) @t @xj (cid:20)(cid:18) (cid:27)!(cid:19) @xj(cid:21) k k (cid:18)@xj @xj(cid:19) k (cid:23)t = f(cid:22) ! with the one-equation ksgs model of Davidson [3] @ksgs @ @ @ksgs + (u(cid:22)jksgs) = ((cid:23) +(cid:23)sgs) @t @xj @xj (cid:20) @xj (cid:21) 3=2 +2C(cid:1)ks1g=s2S(cid:22)ijS(cid:22)ij (cid:0)C(cid:3)ksgs (2) (cid:1) 1=2 (cid:23)sgs = Chom(cid:1)ksgs; see also Section 2.1. The above combination has been investigated for channel (cid:13)ow at a Reynolds number Re(cid:28) = 1050. A mesh with 32(cid:2) 64 (cid:2) 32 (x;y;z) cells has been used. Di(cid:11)erent computational box sizes were used, as given in Table 1. The matching line near the lower wall is located at yml, see Table 1. The cell below the matching line is denoted by jmatch. At the lower matching line, the following boundary conditions are used: @k @! j = jmatch : = = 0 @y @y 2 (3) (cid:0) (cid:23)t j = jmatch +1 : (cid:23)t;jmatch+1 = (cid:23)sgs;jmatch ) ksgs = : (cid:18)Chom(cid:1)(cid:19)jmatch (cid:0) ksgs is the SGS kinetic energy, which is transported by convection-di(cid:11)usion to the LES region. It is introduced via sources and the convection-di(cid:11)usion coe(cid:14)cient connecting 9 L. Davidson + + + Case xmax zmax m.l. y=(cid:14) m.l. jmatch m.l. y (cid:1)x (cid:1)z 1 2(cid:25) 0:5(cid:25) 0.023 4 25 206 52 2 2(cid:25) 0:5(cid:25) 0.057 8 60 206 52 3 4(cid:25) (cid:25) 0.023 4 25 412 104 4 4(cid:25) (cid:25) 0.057 8 60 412 104 5 4(cid:25) (cid:25) 0 0 0 412 104 Table 1: Size of the computational domain and position of the matching line (m.l.) between the LES and RANS regions. The jmatch value represents number of cells in the RANS region. Note that in Case 5 only LES is used. the LES region to the RANS region is set to zero, i.e. aS;jmatch+1 = 0 [13]. Results for mean velocity and shear stress are given in Fig. 7. The velocity pro(cid:12)les for Cases 1 and 2 (restricted computational domain) are seen to agree well with the LES benchmark. A minorkink isvisiblenear the matchingline. When the computationaldomain isincreased (Cases 3 and 4), the agreement is not that good, although it is still much better than when only LES is used (Case 5). The total shear stress for Case 1, Fig. 7b, varies linearly with y, as it should. The sum of SGS shear stress (RANS shear stress near the wall) and the viscous shear stress increases as the wall is approaches. 30 0 25 −0.2 20 −0.4 hu(cid:22)ixzt 15 −0.6 u(cid:28) 10 5 −0.8 0 100 102 −10 0.2 0.4 0.6 0.8 1 + y y=(cid:14) (a) Solid thin line: Case 1 (b) Case 1. Solid thick line: (see Table 1); solid thick total shear stress; solid thin line: Case 2; dashed thin line: resolved shear stress; line: Case 3; dashed thick dashed line: sum of SGS or line: Case 4; dash-dotted RANS shear stress and vis- line: Case 5. cous shear stress. Figure 7: Channel (cid:13)ow. Re(cid:28) =1050, (cid:1)y+ = 2:6(cid:0)97. Markers: LES by Piomelli [15]. a) u(cid:22) pro(cid:12)les. b) shear stresses. Results from Ref. [13,16]. 10