Direct numerical simulation of the flow around a wing section at moderate Reynolds number

ARTICLE IN PRESS

JID: HFF

[m5G;March 26, 2016;19:24]

International Journal of Heat and Fluid Flow 0 0 0 (2016) 1–12

Contents lists available at ScienceDirect

International Journal of Heat and Fluid Flow journal homepage: www.elsevier.com/locate/ijheatfluidflow

Direct numerical simulation of the flow around a wing section at moderate Reynolds number S.M. Hosseini a, R. Vinuesa a, P. Schlatter a,∗, A. Hanifi a,b, D.S. Henningson a a b

KTH Royal Institute of Technology, Department of Mechanics, Linné FLOW Centre, Swedish e-Science Research Centre (SeRC), SE-100 44 Stockholm, Sweden Swedish Defence Research Agency, FOI, SE-164 90 Stockholm, Sweden

a r t i c l e

i n f o

Article history: Available online xxx Keywords: Turbulent boundary layer Vortex shedding Wake Incipient separation Pressure gradient NACA4412

a b s t r a c t A three-dimensional direct numerical simulation has been performed to study the turbulent flow around the asymmetric NACA4412 wing section at a moderate chord Reynolds number of Rec = 40 0, 0 0 0, with an angle of attack of AoA = 5◦ . The mesh was optimized to properly resolve all relevant scales in the flow, and comprises around 3.2 billion grid points. The incompressible spectral-element Navier–Stokes solver Nek50 0 0 was used to carry out the simulation. An unsteady volume force is used to trip the flow to turbulence on both sides of the wing at 10% of the chord. Full turbulence statistics are computed in addition to collection of time history data in selected regions. The Reynolds numbers on the suction side reach Reτ ࣃ 373 and Reθ = 2, 800 with the pressure-gradient parameter ranging from β ≈ 0.0 to β ≈ 85. Similarly, on the pressure side, the Reynolds numbers reach Reτ ≈ 346 and Reθ = 818 while β changes from β ≈ 0.0 to β ≈ −0.25. The effect of adverse pressure gradients on the mean flow is consistent with previous observations, namely a steeper incipient log law, a more prominent wake region and a lower friction. The turbulence kinetic energy profiles show a progressively larger inner peak for increasing pressure gradient, as well as the emergence and development of an outer peak with stronger APGs. The present simulation shows the potential of high-order (spectral) methods in simulating complex external flows at moderately high Reynolds numbers. © 2016 Elsevier Inc. All rights reserved.

1. Introduction A clear evolutionary path can be observed looking back at the aircraft from early days and comparing them to the state of the art modern civil aircraft. This is owed to our deeper knowledge of different flow phenomena around such bodies, i.e., laminar-turbulent transition, wall-bounded turbulence under pressure gradients, flow separation, turbulent wake flow, etc. Despite the conspicuous advances, there still remain major challenges in terms of understanding the complex flow phenomena and a design procedure that can efficiently exploit the interacting features on an airplane. Traditionally such procedures relied heavily on experimental findings. Recent advances in supercomputers and massive parallelization have enabled designers to determine more and more parameters in aircraft design via computational approaches in the early stages of the design, effectively minimizing costs. Recently NASA issued a report by Slotnick et al. (2014) laying out a number of findings and recommendations regarding the role of computational fluid dynamics (CFD) in aircraft design. The main revelations point out ∗

Corresponding author. Tel.: +46 87907176. E-mail address: [email protected] (P. Schlatter).

the necessity of accurate prediction of turbulent flows with significantly separated flow regions, more robust, fast and reliable mesh generation tools and development of more multidisciplinary simulations for both analysis and design optimization procedures. Johnson et al. (2005) also discuss the changing role of CFD from a mere curiosity to a significant role in efficient designs. For instance, in a Boeing–777 significant part of high-speed wing design and propulsion/airframe has been done via CFD simulations. Despite numerous efforts in keeping the flow laminar on aerodynamic surfaces, turbulence flow dominates over a large portion of modern aircraft. In particular pressure gradients (PG) on different surfaces play a major role in the development of turbulent boundary layers (TBL). Some of the early numerical studies of PG TBLs on flat plates using direct numerical simulations (DNS) were the work by Spalart and Watmuff (1993), and Skote et al. (1998), which focused on the role of adverse pressure gradients (APG). This problem has attracted growing attention over the last decade in the recent works by Lee and Sung (2008), Gungor et al. (2014) and Bobke et al. (2016) . Favorable pressure gradients (FPG) have also been recently studied numerically by Piomelli and Yuan (2013). One of the first structure-resolving numerical studies of the flow around a wing was the large-eddy simulation (LES) performed

http://dx.doi.org/10.1016/j.ijheatfluidflow.2016.02.001 S0142-727X(16)30016-9/© 2016 Elsevier Inc. All rights reserved.

Please cite this article as: S.M. Hosseini et al., Direct numerical simulation of the flow around a wing section at the moderate Reynolds number, International Journal of Heat and Fluid Flow (2016), http://dx.doi.org/10.1016/j.ijheatfluidflow.2016.02.001

JID: HFF 2

ARTICLE IN PRESS

[m5G;March 26, 2016;19:24]

S.M. Hosseini et al. / International Journal of Heat and Fluid Flow 000 (2016) 1–12

by Jansen (1996), who considered a NACA4412 profile at a Reynolds number based on freestream velocity U∞ and chord length c of Rec = 1.64 × 106 . The idea behind that study was to compare his results with three different available experimental campaigns of the same flow configuration; those experiments were carried out by Coles and Wadcock (1979), Hastings and Williams (1987) and Wadcock (1987). The flow around the symmetric NACA0012 wing profile was studied by means of DNS by Shan et al. (2005), who considered a Reynolds number of Rec = 10 0, 0 0 0, and an angle of attack of AoA = 4◦ . They found that the backward effect of the disturbed flow on the separated region may be connected to the self-sustained turbulent flow and the self-excited vortex shedding on the suction side of the wing. In addition to this, they found that the vortex shedding from the separated free shear layer was due to a Kelvin–Helmholtz instability. At lower Reynolds number, Rec = 50, 0 0 0, and larger angles of attack of 9.25° and 12°, the same case was analyzed through DNS by Rodríguez et al. (2013). They used a second-order conservative scheme, and found that a combination of leading-edge and trailing-edge stall caused the massive separation on the suction side of the wing. Besides, they also reported that the vortices formed after the initial shear layer undergoes transition are shed forming a von Kármán vortex street in the wake. It should be put in perspective that this range of Reynolds number is significantly smaller than a conventional passenger plane cruising at Re ≈ 108 . Note that the resolution requirements increase with Re37/14 (cf. Choi and Moin, 2012). For instance, increasing the Reynolds number with a factor of 103 dictates a resolution of approximately 80 million times higher. Therefore despite the invaluable information gained from performing well-resolved DNSs at higher Reynolds numbers, they are accompanied by numerous challenges, such as mesh generation, computational resources, and data handling. Here we investigate the turbulent flow around the asymmetric NACA4412 airfoil, with a chord Reynolds number of Rec = 40 0, 0 0 0 and AoA = 5◦ . The airfoil geometry includes a sharp trailing edge, obtained by modifying the last coefficient in the NACA airfoil equation to −0.1036 (see for instance Abbot and von Doenhoff, 1959). This is a big leap in Reynolds number compared to similar previous studies. The Reynolds number corresponds to that of a small cruising glider. Moreover, the use of high-order methods (as opposed to other numerical studies) ensures high-fidelity simulations, especially at the moderate Re under consideration where a significant scale separation starts to emerge. This paper describes the different steps required in performing such large-scale simulations. In addition to full turbulence statistics, time history data is also collected to monitor the different transient behavior of turbulent structures in different regions of the flow.

is mapped isoparametrically from the reference element to the actual grid. The nonlinear terms are treated explicitly by third-order extrapolation (EXT3), whereas the viscous terms are treated implicitly by a third-order backward differentiation scheme (BDF3). In order to avoid quadrature errors, the nonlinear terms are oversampled by a factor of 3/2 in each direction. Furthermore, a low amplitude filter of the highest expansion coefficient (Amp f ilter = 1%) is used in order to avoid errors from potential geometrical aliasing (cf. Fischer and Mullen, 20 01). Nek50 0 0 has previously been used to study turbulent flow in a pipe by El Khoury et al. (2013) and flow around a wall-mounted square cylinder by Vinuesa et al. (2015) at moderately high Reynolds-numbers. The following results have been obtained using comparably high order, N = 11. The present SEM code is optimized for modern computers with thousands of processors Tufo and Fischer (2001). Here, we have performed parallel computations on 16,384 processors. Scaling tests performed on the Cray-XC40 computer Beskow at PDC (KTH) are shown in Fig. 1, for the case under consideration in the present study (left), which has a total of 3.2 billion grid points (and around 1.85 million spectral elements), and for a similar configuration with a total of 120 million grid points (right). In order to establish a good measure of scaling, we report the time required to perform one GMRES (generalized minimal residual method) iteration for the pressure solve. Doing so, we can characterize code performance by isolating it from other factors contributing to the total time per time-step, such as I/O operations, different tolerances etc. Note that minor superlinear scaling is observed in Fig. 1 (left) up to 32,768 cores. This can be explained by the fact that Nek50 0 0 exhibits best parallel performance when the number of cores is a power of 2, which was not the case for the simulation with the lowest number of processes (12,0 0 0 cores) due to the required memory. Interestingly, the scaling tests shown in Fig. 1 (right), which were performed on numbers of cores which were powers of 2, exhibit linear scaling up to 8192 cores. It is also interesting to note that in the case under consideration in this study the best parallel efficiency is achieved when running on 32,768 cores (with around 10 0,0 0 0 grid points per core), and in the smaller test case the best performance is observed on 4,096 cores, with a total of 30,0 0 0 grid points per core. A total of 35 million core-hours were spent to collect full turbulence statistics, time history data, and flow snapshots for visualization purposes. The flow was initially run for three flow-over times with N = 5, and after this point the polynomial order was progressively increased up to the final value of N = 11.

2. Direct numerical simulation

Initially a RANS simulation was performed by means of EDGE code developed at FOI, the Swedish Defence Research Agency (Eliasson, 2002). The code uses the explicit algebraic Reynolds stress model (EARSM) by Wallin and Johansson (20 0 0). The RANS domain extends up to 200c in every direction. The solution from the RANS simulation is used to extract an accurate velocity distribution in the near field corresponding to the averaged flow at a given angle of attack. Consequently these values are imposed as Dirichlet boundary conditions on the DNS domain. Note that the natural stress-free condition was used at the outlet, and periodicity was imposed in the spanwise direction. We considered a Cmesh topology of radius c centered at the leading edge of the airfoil, with total domain lengths of 6.2c in the horizontal (x), 2c in the vertical (y) and 0.1c in the spanwise (z) directions, see Fig. 2. The resulting spectral-element mesh (without the GLL points) and the computational domain are shown in Fig. 3. In this figure it can be observed that a local smoothing was applied to a very small region close to the sharp trailing edge.

2.1. Nek5000 Direct numerical simulations were performed using the incompressible Navier–Stokes solver ‘Nek50 0 0’ by Fischer et al. (2008), which uses the spectral-element method proposed (SEM) by Patera (1984). Enabling the geometrical flexibility characteristic of finiteelement methods combined with the accuracy provided by spectral methods are the main advantages of using SEM. The equations are cast into weak form. The spatial discretization is then done by means of the Galerkin approximation, following the PN − PN−2 formulation. The solution to the Navier–Stokes equations is approximated element-wise as a weighted sum of Lagrange interpolants defined by an orthogonal basis of Legendre polynomials up to degree N (of degree N − 2 in the case of the pressure). Inside each element, Gauss–Lobatto–Legendre (GLL) quadrature points are used, and extended to 3D in a tensor-product formulation. The geometry

2.2. Boundary conditions, mesh design and simulation procedure

Please cite this article as: S.M. Hosseini et al., Direct numerical simulation of the flow around a wing section at the moderate Reynolds number, International Journal of Heat and Fluid Flow (2016), http://dx.doi.org/10.1016/j.ijheatfluidflow.2016.02.001

JID: HFF

ARTICLE IN PRESS

[m5G;March 26, 2016;19:24]

S.M. Hosseini et al. / International Journal of Heat and Fluid Flow 000 (2016) 1–12

3

Fig. 1. Strong scaling for problem sizes of (left) 3.2 billion and (right) 120 million grid points vs. the number of cores on a Cray XC40 system. The dashed line shows the linear scaling.

As will be discussed below, the methodology proposed here using the RANS solution as a Dirichlet boundary condition shows excellent results in this configuration, with low angle of attack and far from stall. However, in cases closer to stall with large-scale unsteady separation it is important to carefully assess the impact of the RANS solution on the flow dynamics developing on the suction side of the wing. Future work will be devoted to the development of a general methodology applicable to such cases. Mesh generation is one of the key factors of the present simulation, where the computational grid must meet the requirements for a fully-resolved DNS by allowing sufficient geometrical flexibility. The mesh is structured, and is designed to meet the following criteria around the wing section: x+ < 10 (defined in the direction tangential to the wing surface), y+ w < 0.5 (at the wing, in the direction normal to its surface) and z+ < 5. Note that here ‘+’ denotes scaling with the viscous length ∗ = ν /uτ (where ν is the √ fluid kinematic viscosity and uτ = τw /ρ is the friction velocity, ρ being the fluid density). Since the wall shear stress τ w is not constant over the wing surface, and not know a-priori, we considered the RANS results to design the mesh according to the following guidelines: Fig. 2. Schematic three-dimensional layout of the set up for direct numerical simulation. The chord length is denoted as c. The domain extends to 5c downstream, and 1c upstream, top and bottom from the leading edge. The red lines mark the boundary where the Dirichlet boundary condition has been used. The blue line indicate the boundaries with the stress-free boundary condition. The incoming flow has an angle of attack of AoA = 5◦ . The Cartesian coordinate system was aligned with the chord. (For interpretation of the references to color in this figure legend, the reader is referred to the web version of this article.)

I Between x/c = 0.1 and 0.6, the uτ from the RANS was used to compute x, both on suction (ss) and pressure sides (ps). II For x/c > 0.6, the uτ from the RANS on the pressure side was used to design the mesh in x, both on the suction and pressure sides. The strong APG on the suction side in this region leads to a progressively smaller value of uτ , which in principle would yield a coarser mesh. However, we considered that at least the

Fig. 3. Spectral element mesh used for the direct numerical simulation. Note that the mesh was designed based on the variation of the Kolmogorov scale in the turbulent boundary layers and the wake.

Please cite this article as: S.M. Hosseini et al., Direct numerical simulation of the flow around a wing section at the moderate Reynolds number, International Journal of Heat and Fluid Flow (2016), http://dx.doi.org/10.1016/j.ijheatfluidflow.2016.02.001

ARTICLE IN PRESS

JID: HFF 4

[m5G;March 26, 2016;19:24]

S.M. Hosseini et al. / International Journal of Heat and Fluid Flow 000 (2016) 1–12

resolution of the pressure side should be used in this region of the upper surface. III For x/c < 0.1, the respective values of uτ at x/c = 0.1 on suction and pressure sides were used to design the mesh in x. IV The spectral element mesh is uniform in z, and the value of uτ at x/c = 0.2 on the suction side was used to define the constant spanwise width of the spectral elements. We considered the value of uτ at this location because it is the maximum value in the turbulent region. V The mesh around the wing is orthogonal, and the yw value was computed based on the value of uτ close to the trailing edge on the pressure side. Note that as observed in Fig. 3 a local smoothing was applied to a small region close to the sharp trailing edge, and therefore at this particular location the mesh is not orthogonal to the wing surface. Fig. 4 shows the final mesh resolution around the wing in the directions tangential and normal to the wing surface, as well as in the spanwise direction, obtained following the previous guidelines. Note that the friction velocity distribution used in this figure was obtained from the time and spanwise averaged DNS data and the wall-shear stress was expressed in the tangential direction. In Fig. 4 (top) the maximum and minimum x+ from each spectral element is represented, for both suction and pressure sides. The flow is laminar from x/c = 0 (the leading edge of the wing) to x/c = 0.1, the location where the flow is tripped. After this point the flow transitions to turbulence, and for x/c ࣡ 0.2 it can be considered to be fully turbulent. Note the peak in the x+ curves at the tripping location on both sides of the wing, and how the criterion x+ < 10 is fulfilled throughout the whole turbulent region. It is also interesting to note that on the suction side the friction velocity decreases significantly due to the effect of the APG, whereas the very mild FPG on the pressure side leads to a slightly increasing uτ trend, and both effects are reflected in all the resolution plots. Fig. 4 (middle) shows how the condition for the first point away from the wall is also satisfied in the turbulent regions of both surfaces, as well as the condition in the spanwise direction which is shown in Fig. 4 (bottom). Here we show the minimum and maximum spanwise spacings (zmin and zmax ), which are constant throughout the whole domain, scaled with the local friction velocity values. The location of the tripping and the uτ trends can also be observed in Fig. 4 (middle) and (bottom). So far we have discussed typical mesh requirements for DNS, which have been developed for canonical flow cases where it is easy to define the wall shear stress and the parallel/orthogonal directions. However, they do not necessarily apply to the mesh around a wing section, especially farther from the wall, and therefore we considered additional criteria to ensure that the mesh requirements were met everywhere in the domain. To this end, we used distributions of the Kolmogorov scale η obtained from the RANS to optimize the mesh. Note that the Kolmogorov scale is



1/4

computed as η = ν 3 /ε , where ε is the local isotropic dissipation. State-of-the-art high-Reynolds number DNSs (Schlatter and Örlü, 2010; Sillero et al., 2013) typically follow y(y ) ∼ 4–8η (y ) of boundary layers, which together with the previous restrictions in x and z provide sufficient mesh resolution for a fully resolved DNS. In our case, our mesh was designed in order to satisfy the condition h ≡ (x · y · z)1/3 < 5η everywhere in the domain, which ensures that the mesh is fine enough to capture the smallest turbulent scales. A similar criterion was also adopted in the wing DNSs on unstructured meshes by Rodríguez et al. (2013). In Fig. 5 we show the wall-normal evolution of the ratio hmax /η (obtained by considering the maximum spanwise spacing zmax ) at the wing locations x/c = 0.3 and 0.95, for both suction and pressure sides. Note that the Kolmogorov scale profiles were also obtained from the time- and spanwise- averaged DNS data, and that

Fig. 4. Streamwise evolution of mesh resolution in (top) direction tangential to the wing surface, (middle) direction normal to the wing surface at the wall and (bottom) spanwise direction. Values are shown for the suction and the pressure sides of the wing.

the coordinate yn is defined in the direction normal to the wing surface. The first important observation is that this ratio reaches maximum values of around 3, which in fact implies that the mesh resolution is very fine around the wing on both sides. The growth of the Kolmogorov scale with yn is similar on both sides of the wing at x/c = 0.3, whereas at x/c = 0.95 the hmax /η ratio starts to decrease much farther away from the wall than on the pressure side. This is of course due to the fact that the very strong APG on the suction side significantly thickens the boundary layer, which is also observed in the value of the 99% boundary layer thickness δ 99 . At x/c = 0.95 the value of δ 99 is around two times larger on the top surface than on the bottom. The oscillations observed in these

Please cite this article as: S.M. Hosseini et al., Direct numerical simulation of the flow around a wing section at the moderate Reynolds number, International Journal of Heat and Fluid Flow (2016), http://dx.doi.org/10.1016/j.ijheatfluidflow.2016.02.001

JID: HFF

ARTICLE IN PRESS

[m5G;March 26, 2016;19:24]

S.M. Hosseini et al. / International Journal of Heat and Fluid Flow 000 (2016) 1–12

5

Fig. 5. Wall-normal evolution of the hmax /η ratio at (left) x/c = 0.3 and (right) x/c = 0.95, for both suction and pressure sides of the wing. The corresponding 99% boundary layer thickness values δ 99 are also shown for the four profiles.

Fig. 6. Vertical evolution of the hmax /η ratio at different locations in the turbulent wake downstream of the wing section.

profiles are produced by the particular GLL distribution of grid points within spectral elements. The quality of the mesh is also assessed in the wake, where in Fig. 6 (left) we show the hmax /η ratios at x/c ࣃ 1.1 and 2, and in Fig. 6 (right) at x/c ࣃ 2.9 and 4.7. The maximum value of the ratio is found at x/c ࣃ 1.1 close to y/c ࣃ 0, and it is around 3. Beyond this point the maximum values of this ratio further decrease, which again indicates that the mesh is very fine in the wake and therefore able to properly represent the turbulent wake. Another interesting observation inferred from this figure is the fact that all the hmax /η profiles show a value close to zero at approximately y/c ࣃ 0, then this ratio increases both above and below y/c ࣃ 0 reaching two relative maxima, and eventually goes towards zero

for large enough values of y. The point where hmax /η ࣃ 0 between the two relative maxima shows the location of maximum velocity defect in the wake, and therefore its streamwise evolution can be inferred from this figure. Interestingly, the hmax /η ࣃ 0 locations form a straight line in the xy plane, with an angle of ࣃ 3.2° with respect to the horizontal plane. Although the incidence angle was 5°, the fact that the boundary layer on the suction side is much thicker than the one on the pressure side tilts down the angle of the wake after they merge. This is indeed a manifestation of the downwash effect at the trailing edge of the wing. Also note the asymmetry in hmax /η: the relative maximum below the zero value is larger than the one above, which implies that the Kolmogorov scale is smaller and therefore the turbulent structures on the wake

Please cite this article as: S.M. Hosseini et al., Direct numerical simulation of the flow around a wing section at the moderate Reynolds number, International Journal of Heat and Fluid Flow (2016), http://dx.doi.org/10.1016/j.ijheatfluidflow.2016.02.001

ARTICLE IN PRESS

JID: HFF 6

[m5G;March 26, 2016;19:24]

S.M. Hosseini et al. / International Journal of Heat and Fluid Flow 000 (2016) 1–12

Fig. 7. Visualization of the instantaneous vertical forcing field located at x/c = 0.1. The blue and red iso-surfaces show the positive and negative vertical force. The Gaussian shape of the forcing field is depicted in a two-dimensional slice via the iso-contours. The pseudo-colors on each slice show the streamwise velocity. (For interpretation of the references to color in this figure legend, the reader is referred to the web version of this article.)

originated in the pressure side of the wing are smaller than the corresponding ones in the thicker boundary layer from the suction side. Note that the vertical grid spacing y is slightly larger on the lower part of the curve, which also contributes to the larger hmax /η ratio, but both x and z are the same on both sides, and these constitute 2/3 of the value of hmax . 2.3. Tripping In order to initiate the turbulent boundary layer two tripping strips are placed at x/c = 0.1 on both sides of the wing. This is achieved by applying a weak random volume forcing as reported in Schlatter and Örlü (2012). The direction of the forcing is only in the vertical direction. The volume forcing reproduces the same effect that tripping strips have in wind-tunnel experiments. The volume forcing F2 enters into the momentum equation, whereby amplitude (At ), spanwise length scale (zs ) and temporal frequency (ts ) are the relevant parameters. Hence F2 is defined as



 2

F2 = exp [(x − x0 )/lx ]2 − [y/ly ] g(z, t ),

(1)

where the forcing g(z, t) varies in time and spanwise direction, and it is mapped into a Gaussian shape, with lx and ly being the Gaussian attenuation of the forcing region. The distribution in time and space of g(z, t) reads

g(z, t ) = At



 {1 − b(t )}hi (z ) + b(t )hi+1 (z ) ,

(2)

with b(t ) = 3 p2 − 2 p3 , p = t/ts − i and i = int(t/ts ), where int(·) represents the integer part of the argument. The temporal scale is represented by ts , which is in effect the cut-off frequency of the perturbations. Third-order Lagrangian interpolation is used in order to interpolate g(z, t) in time. The values of hi (z) are random harmonic functions having a unit amplitude for wavenumbers below 2π /zs and zero amplitudes for wavenumbers above that range. A pseudo-random number generator, which enables accurate restart of the simulations, was used to generate the spanwise distribution. The Gaussian attenuation of the forcing function was optimally chosen to be lx = 4δx∗t and ly = δx∗t , with xt denoting the location of the trip forcing at xt /c = 0.1. This location corresponds to Reθ ࣃ 180 on the suction side and Reθ ࣃ 100 on the pressure side. Similarly a temporal cut-off scale of ts = 4δx∗t /U∞ and the spanwise cut-off scale of zs = 1.7δx∗t were chosen. The shape of the forcing is depicted in Fig. 7. On both sides of the wing the main

Fig. 8. Instantaneous visualization showing coherent vortices in the region close to the location of the tripping strip on the suction side of the wing. Note the emergence of hairpin vortices visualized by the λ2 -criterion (Jeong and Hussain, 1995), = −0.0 0 01. The structures are colored with chordwith the isosurfaces set at λ+ 2 wise velocity, ranging from -0.1 (dark blue) to 1.5 (red). (For interpretation of the references to color in this figure legend, the reader is referred to the web version of this article.)

axis of the forcing shape is aligned with the surface tangent. The amplitude was adjusted in order to have immediate transition at the location of the trip. This prevents the transition location moving back and forward along the wing, giving a consistent transition point. Fig. 8 shows vortex visualization of coherent vortices obtained from an instantaneous field by means of the λ2 criterion by Jeong and Hussain (1995). It clearly shows the emergence of hairpin vortices immediately after the tripping strip, which lead to the development of a turbulent boundary layer. 3. Results Fig. 9 shows an instantaneous visualization of vortical structures using the λ2 criterion. The smoothness of the hairpin vortices, characteristic of transitional flows, and of the smallest structures in the turbulent region shows that the mesh is appropriate to capture all the relevant flow features. It is interesting to observe how the turbulent vortical structures occupy a significantly

Please cite this article as: S.M. Hosseini et al., Direct numerical simulation of the flow around a wing section at the moderate Reynolds number, International Journal of Heat and Fluid Flow (2016), http://dx.doi.org/10.1016/j.ijheatfluidflow.2016.02.001

JID: HFF

ARTICLE IN PRESS

[m5G;March 26, 2016;19:24]

S.M. Hosseini et al. / International Journal of Heat and Fluid Flow 000 (2016) 1–12

7

Fig. 9. Instantaneous vortical structures visualized by the λ2 criterion, colored with chordwise velocity, ranging from −0.1 (dark blue) to 1.5 (red). The flow is tripped at 10% chord on both sides. The angle of attack is AoA = 5◦ and the chord Reynolds number is Rec = 40 0, 0 0 0. The spectral element mesh is also shown, but not the individual grid points within elements. (For interpretation of the references to color in this figure legend, the reader is referred to the web version of this article.)

Fig. 10. Instantaneous vortical structures visualized by λ2 criterion, colored with chordwise velocity, ranging from −0.1 (dark blue) to 1.1 (red). (For interpretation of the references to color in this figure legend, the reader is referred to the web version of this article.)

Please cite this article as: S.M. Hosseini et al., Direct numerical simulation of the flow around a wing section at the moderate Reynolds number, International Journal of Heat and Fluid Flow (2016), http://dx.doi.org/10.1016/j.ijheatfluidflow.2016.02.001

JID: HFF 8

ARTICLE IN PRESS

[m5G;March 26, 2016;19:24]

S.M. Hosseini et al. / International Journal of Heat and Fluid Flow 000 (2016) 1–12

Mach numbers (Ma ≈ 0.04). Since no pressure input from RANS is delivered to DNS, a comparison of pressure coefficient obtained from RANS and DNS gives further verification on the validity of the set up. This comparison is achieved by adjusting the pressure level such that the pressure coefficient at the stagnation (C pst ) is unity point for both cases using

C pst =

Fig. 11. Comparison of the time- and spanwise-averaged chordwise velocity from ) and the RANS (− − −) solution. the DNS (

larger volume on the suction side than on the pressure side, due to the fact that the APG on the top surface significantly thickens the boundary layer, and also leads to more energetic large-scale structures in the outer region of the boundary layer. A close up view of the wake region up to two chord lengths downstream of the airfoil is shown in Fig. 10. Near the trailing edge the structures on the upper side travel through a low speed region, colored in blue. The near-wall structures from both sides start to interact as they leave the trailing edge, and a region of reverse flow right after the trailing edge of the wing is observed. The reverse flow also appears instantaneously near the very end of the wing, while in the mean the flow remains attached. This relates to the similar observation by Shan et al. (2005) where they found a relation between the backward effect of the separated region on the suction side with the self-excited vortex shedding. Further downstream the emergence of a von Kármán-type of vortex street can be observed. The accuracy of the set up is confirmed by comparing the streamwise mean velocity from the RANS and the DNS (in this case obtained from the time- and spanwise-averaged field), as shown in Fig. 11. The excellent agreement obtained around the wing and in the wake is remarkable, and also highlights the quality of the setup used in the present study. In Nek50 0 0 the pressure and velocity are decoupled in the incompressible Navier–Stokes formulation. At each time step a non-divergence free velocity field is computed using the pressure field from the previous time step, followed by a pressure correction. On the other hand, the RANS is performed using the compressible Navier–Stokes equations at extremely low

1 Pst − P∞ 2 = 1 ⇒ P∞ = Pst − ρU∞ , 2 2 1/2ρU∞

(3)

where P∞ and Pst are the ambient and stagnation  point pressure. 2 is then computed The pressure coefficient C p = (P − P∞ )/ 1/2ρU∞ using the new pressure levels. The comparison is shown in Fig. 12a. A remarkably good agreement can be observed. Moreover Fig. 12b shows the friction coefficient C f = 2 ) from the DNS results. A sudden rise of friction τw /(1/2ρU∞ coefficient near the tripping strip can be clearly observed and this value later settles to a turbulent state at x/c = 0.2 on the suction side. Subsequently it continues to decline as expected from an APG turbulent boundary layer, whereas on the pressure side the value stays approximately constant with a very slight increase due to the FPG as discussed in Section 2.2. The friction coefficient value drops two orders of magnitude on the suction side while on the pressure side it exhibits a small variation. This behavior appears to be consistent with the reported values of friction coefficient of a symmetric NACA0012 airfoil by Wolf et al. (2012). A comparison with different RANS solutions with different prescribed transition locations xtr is also presented in Fig. 12b. The first interesting observation is the excellent agreement between DNS and RANS in the laminar region from both sides before the tripping location in the DNS, x/c = 0.1. At this location, the skin friction exhibits a sharp peak on both sides associated with the tripping, after which both boundary layers undergo different transitional stages. Note that the various tripping locations prescribed in the RANS allow to observe the evolution of the laminar Cf , and interestingly the skin friction curve from the pressure side returns to the laminar trend after the tripping location up to around x/c ࣃ 0.4. After this point, the Cf on the pressure side shows excellent agreement with the RANS, and therefore it can be argued that in this region the boundary layer in the DNS exhibits fully-developed turbulence features. With respect to the suction side, the boundary layer does not exactly follow the laminar curve after the tripping location, but it can be considered to be a fully-developed TBL after x/c ࣃ 0.2, i.e. earlier than the one on the pressure side. Although there is a small deviation between DNS and RANS for x/c > 0.2 on the suction side, the trends from both simulation techniques exhibit a very similar behavior. Note that small deviations are observed

Fig. 12. Pressure coefficient Cp (left) and skin friction coefficient Cf (right). The pressure side and suction side are denoted by ps and ss on both panels. Skin friction coefficients obtained from various RANS solutions at different prescribed transition locations xtr = 10%, 15% and 20% chord length are also presented.

Please cite this article as: S.M. Hosseini et al., Direct numerical simulation of the flow around a wing section at the moderate Reynolds number, International Journal of Heat and Fluid Flow (2016), http://dx.doi.org/10.1016/j.ijheatfluidflow.2016.02.001

ARTICLE IN PRESS

JID: HFF

[m5G;March 26, 2016;19:24]

S.M. Hosseini et al. / International Journal of Heat and Fluid Flow 000 (2016) 1–12

9

Fig. 13. Spanwise and time-averaged streamwise velocity. The black rectangle shows the control volume bound ( ) over which the momentum equation is integrated.

in the RANS solutions with the various prescribed transition locations, which points out the dependence of RANS on the initial parameters and highlights the importance of finely resolved simulations such as the one presented here for accurate predictions of relevant aerodynamic quantities. The resulting lift and drag coefficients from the DNS obtained by integration  of 2the pressure and friction forces on the surface are Cl = L/ 1/2ρU∞ c = 0.905,





2 c = 0.01681 (where L and D denote lift and and Cd = D/ 1/2ρU∞ drag forces). With respect to the RANS simulations, the Cl values are 0.893, 0.898 and 0.902 for the prescribed transition locations of 10%, 15% and 20% respectively, whereas the corresponding Cd values are 0.01720, 0.01681 and 0.01642. An alternative approach can also be considered by computing the forces acting on a control volume as indicated in Fig. 13. The time-averaged velocity components are first interpolated onto an equidistant Cartesian grid. The pressure and momentum flux forces are integrated around the bound and the corresponding lift and drag coefficients are calculated using

D∗ = − − L∗ = − −



ρ (u¯ · dl )u¯ −  







ρ u u dl



−





x

ρ u v dl



x

p¯ dl , x

ρ (u¯ · dl )v¯ −    



p¯ dl ,

(4)

 

ρv v dl

 y

−

 

ρ u v dl

 y

(5)

y

where u¯ is the mean velocity vector (u¯ , v¯ ) in Cartesian coordinates with u u and u v the corresponding Reynolds stresses. The mean pressure is denoted by p¯ and the integral bound  is marked in Fig. 13. The vector normal and outward on each segment on  is represented by dl. The subscripts, [·]x and [·]y denote the projection in horizontal and vertical Cartesian coordinates respectively. The resulting D and L forces are then obtained by rotating D∗ and L∗ such that the forces are aligned with the incoming flow. The resulting lift and drag coefficients from this integral method are Cl = 0.904, and Cd = 0.01682, i.e., the deviations with respect to the values obtained by integrating the forces on the wing surface are Cderr = 0.1% and Cl err = 0.06% for the drag and lift coefficients respectively. As is clear from Fig. 12a, the boundary layers on both sides of the wing are subjected to streamwise pressure gradients, adverse in the case of the suction side, and slightly favorable in the pressure side. In order to characterize the pressure gradient at any particular location of the wing, we consider the Rotta– Clauser pressure gradient parameter β = δ ∗ /τw dPe /dxt , where δ ∗ is the displacement thickness and dPe /dxt is the gradient of the pressure at the boundary layer edge in the direction tangential to the wing surface. The boundary layer edge was determined using the diagnostic-plot concept (Alfredsson et al., 2011), which as shown

Fig. 14. Reynolds number based on momentum thickness Reθ for the suction and pressure sides of the wing.

by Vinuesa et al. (2016) is a robust method to determine the 99% boundary-layer thickness δ 99 in pressure gradient TBLs, and will be considered in the present study. The suction side of the wing experiences a significant increase in pressure gradient magnitude, with β values ranging from around 0 up to 85. On the pressure side the pressure gradient range is quite narrow, and β changes from around 0 to −0.25. The Reynolds number based on momentum thickness (denoted by θ ) is shown in Fig. 14 as a function of the chordwise coordinate x. The different evolutions on both sides of the wing are clear from this figure, where the progressively more intense adverse pressure gradient on the suction side leads to a more rapid growth rate in Reynolds number. In fact, the Reθ curve on the suction side follows a roughly exponential growth beyond x/c ࣃ 0.8, where the value of β is around 4.1. This figure also shows that the maximum Reynolds number based on momentum thickness is Reθ = 2, 800 on the suction side, and Reθ = 818 on the pressure side. With respect to the friction Reynolds numbers Reτ = δ99 uτ /ν, the maxima are ࣃ373 and ࣃ346 for the suction and pressure sides, respectively. Fig. 15 (left) shows displacement and momentum thicknesses as a function of the streamwise coordinate for both suction and pressure sides of the wing. The velocity profiles used to evaluate these quantities were projected on the directions tangential (t) and normal (n) to the wing surface, at a total of 80 locations throughout both sides. Note that the method used to determine δ 99 is based on definitions valid for turbulent boundary layers, and therefore we do not show results for x/c < 0.2 due to the fact that in this region the flow is laminar or transitional. It can be observed how for x/c ࣠ 0.4 the rate of growth of δ ∗ and θ is very similar on both sides of the wing, due to the fact that the pressure gradient is very mild on both surfaces. After this point, the boundary layer

Please cite this article as: S.M. Hosseini et al., Direct numerical simulation of the flow around a wing section at the moderate Reynolds number, International Journal of Heat and Fluid Flow (2016), http://dx.doi.org/10.1016/j.ijheatfluidflow.2016.02.001

JID: HFF 10

ARTICLE IN PRESS

[m5G;March 26, 2016;19:24]

S.M. Hosseini et al. / International Journal of Heat and Fluid Flow 000 (2016) 1–12

Fig. 15. (Left) momentum (red) and displacement thickness (black) for suction ( ) and pressure sides (− − −) . (Right) Shape factor for suction ( (− − −). (For interpretation of the references to color in this figure legend, the reader is referred to the web version of this article.)

) and pressure sides

Fig. 16. Inner-scaled mean velocity profiles, where results from the present simulation are represented by the solid (suction side) and dashed (pressure side) lines. Circles denote the DNS of ZPG turbulent boundary layer by Schlatter and Örlü (2010). (Left) x/css = 0.4 (suction side) and x/c ps = 0.65 (pressure side), Reτ ࣃ 250. (Right) x/css = 0.7 and x/c ps = 0.98, Reτ ࣃ 360.

Fig. 17. Turbulence kinetic energy profiles, where results from the present simulation are represented by the solid (suction side) and dashed (pressure side) lines. Circles denote the DNS of ZPG turbulent boundary layer by Schlatter and Örlü (2010). (Left) x/css = 0.4 (suction side) and x/c ps = 0.65 (pressure side), Reτ ࣃ 250. (Right) x/css = 0.7 and x/c ps = 0.98, Reτ ࣃ 360.

developing on the pressure side does not exhibit a significant change in the rate of growth, whereas on the suction side both thicknesses increase exponentially due to the progressively increasing magnitude of the APG. For x/c ࣡ 0.9 the APG is so strong that the boundary layer is very close to detachment, with the corresponding significant decrease in wall-shear stress as discussed in Fig. 12. In this region the flow detaches instantaneously, but the positive Cf value indicates that the averaged profiles are attached up to the trailing edge. Also note that the APG strongly reduces wall shear by lifting up the profile and therefore an APG boundary layer would carry more momentum than a ZPG one

overcoming the same friction. This can also be observed in Fig. 15 (right), where the shape factor H = δ ∗ /θ is shown for both sides of the wing. Whereas H exhibits a decreasing trend on the pressure side, similar to the one characteristic of ZPG boundary layers, the rate of growth of H on the suction side shows a significant increase for x/c ࣠ 0.4, which is also noticeable in Fig. 15 (left) due to the much faster growth of δ ∗ than of θ . This is another manifestation of the effect of the pressure gradient, which leads to a much thicker boundary layer for basically the same momentum. At x/c ࣃ 0.9, the boundary layer on the suction side exhibits a shape factor of H ࣃ 2.03 and a pressure gradient parameter of β ࣃ 14.1; the

Please cite this article as: S.M. Hosseini et al., Direct numerical simulation of the flow around a wing section at the moderate Reynolds number, International Journal of Heat and Fluid Flow (2016), http://dx.doi.org/10.1016/j.ijheatfluidflow.2016.02.001

ARTICLE IN PRESS

JID: HFF

[m5G;March 26, 2016;19:24]

S.M. Hosseini et al. / International Journal of Heat and Fluid Flow 000 (2016) 1–12 Table 1 Summary of H and β cases compared with the literature. Note that ss and ps stand for suction and pressure sides respectively. Reference

Reθ

Reτ

H

β

Present x/css = 0.9 Present x/c ps = 0.9 Present x/css = 0.7 Schlatter and Örlü (2010) Skåre and Krogstad (1994) Monty et al. (2011) Bobke et al. (2016)

2255 785 1415 1007 39,120 11,860 4300

328 317 356 359 4,300 2,860 600

2.03 1.48 1.66 1.45 2.00 1.50 1.70

14.1 −0.16 2.4 ࣃ0 19.9 1.9 2.0

corresponding values on the pressure side are H ࣃ 1.48 and β −0.16 (note that the negative sign indicates FPG). The Reynolds numbers at this location are Reτ ࣃ 328 and Reθ ࣃ 2, 255 on the suction side, and Reτ ࣃ 317 and Reθ ࣃ 785 on the pressure side. For similar Reynolds numbers, the ZPG boundary layer simulated by Schlatter and Örlü (2010) has a shape factor of H ࣃ 1.45, which is very similar to the FPG value on the pressure side due to the very mild PG. The strong APG on the suction side leads to a value of H approximately 38% larger than the ZPG value, which further highlights the significant thickening of the boundary layer. The previous comparison of shape factor values H in terms of β , together with the moderate APG case shown in Fig. 16 (right) below, is summarized in Table 1. The results from this study are also compared with data available in the literature, namely the ZPG DNS from Schlatter and Örlü (2010), the strong APG experiment by Skåre and Krogstad (1994), the moderate APG experiments by Monty et al. (2011) and the LES of moderate APG by Bobke et al. (2016). The first conclusion from Table 1 is that, as stated above, the mild FPG boundary layer on the pressure side of the wing exhibits a very similar shape factor as the ZPG by Schlatter and Örlü (2010). The strong APG found at x/c ࣃ 0.9 on the suction side of the wing exhibits almost the same shape factor H ࣃ 2 as the one from the experiment by Skåre and Krogstad (1994) (despite the difference in Reynolds number of an order of magnitude), with a comparable β value of 19.9. Regarding the moderate APG at x/c ࣃ 0.7 on the suction side, the agreement in shape factor with the simulation by Bobke et al. (2016) is quite remarkable (H ࣃ 1.7), for a Reynolds number around 3 times larger. The profile at x/c ࣃ 0.7 also agrees well with the experimental data from Monty et al. (2011), with H ࣃ 1.5, at a Reynolds number around 8 times larger. The inner-scaled mean flow is shown in Fig. 16 at two locations on the suction and pressure sides, and compared with the ZPG results by Schlatter and Örlü (2010). Following the approach by Monty et al. (2011), we selected profiles that roughly exhibited the same friction Reynolds numbers Reτ : Fig. 16 (left) shows profiles with Reτ ࣃ 250 (at x/c ࣃ 0.4 on the suction side and 0.65 on the pressure side), whereas in Fig. 16 (right) the Reτ is approximately 360 (obtained at x/c ࣃ 0.7 on the top surface and 0.98 on the bottom one). In Fig. 16 (left) the value of β on the suction side is 0.6, which is associated with a moderate adverse pressure gradient. Localized pressure gradients of similar magnitude were studied experimentally and computationally by Vinuesa et al. (2014) over a higher Reynolds number range up to Reθ ࣃ 35, 0 0 0. The effect of this moderate pressure gradient on the flow is noticeable especially in the outer part of the boundary layer, where a more prominent wake is observed, together with a higher value of inner-scaled edge velocity. Interestingly, the mild FPG on the pressure side leads to an inner-scaled profile in very good agreement with the corresponding ZPG boundary layer. The profile from the suction side shown in Fig. 16 (right) has a higher value of β ࣃ 2.4, representative of some of the intermediate pressure gradients analyzed by Monty et al. (2011) on boundary layers with Reθ up to 18,500. The stronger pressure gradient shows a significant effect in

11

the outer flow, where the Ue+ value is around 20% larger, but is also noticeable in the incipient logarithmic region and all the way down to the buffer layer. In this case the shape factor from the profile on the suction side of the wing is 1.66, and the one from the ZPG boundary layer is 1.45. The steeper logarithmic region exhibited by this profile is also characteristic of APG boundary layers, which is associated with lower values of the von Kármán coefficient κ as discussed by Nagib and Chahuan (2008). In this case the boundary layer on the pressure side differs slightly from the ZPG case due to the effect of the mild FPG it is subjected to, defined by a value of β −0.24. Additional insight on the effect of APGs on the boundary layers developing around the wing can be gained by analyzing the turbulence kinetic energy (TKE), which is shown in Fig. 17 for the same cases as in Fig. 16. Note that here we use the components projected on the t and n directions to

compute the turbulence kinetic energy, i.e., k = 1/2 ut2 + v2n + w2 , and therefore the spanwise tur-

bulence intensity w2 remains unchanged. It can be observed how at x/c ࣃ 0.4 the APG on the suction side has a large impact on the TKE : the inner peak is increased, and the effect on the outer region is quite noticeable, as also observed by Monty et al. (2011). As in the case of the mean flow, the mild FPG on the pressure side leads to a very similar TKE profile compared with the corresponding ZPG boundary layer. At x/c ࣃ 0.7, the APG greatly affects the TKE profile on the suction side, where the moderately strong pressure gradient leads to a larger inner peak, and more interestingly to a prominent outer peak which was also observed by Monty et al. (2011) and is due to the fact that the APG energizes the largest and most energetic turbulent structures in the flow. To a lower extent, the opposite effect is observed in the FPG acting on the boundary layer in the pressure side, which exhibits slightly lower values of TKE, including a less prominent inner peak. The decrease in turbulent fluctuations was also observed in the FPG boundary layers simulated by Piomelli and Yuan (2013). 4. Conclusions A direct numerical simulation of the flow around a wing section at a moderate Reynolds number of Rec = 40 0, 0 0 0 at an angle of attack of 5° has been performed. Computations based on the Reynolds-averaged Navier–Stokes (RANS) equations were initially carried out to determine the velocity distribution around the wing. The simulations are then done using the spectral-element code Nek50 0 0. This is the first time well-resolved direct numerical simulations of the flow around a wing are performed at such high Reynolds number. Turbulence statistics are computed over the suction and pressure sides of the wing, and expressed in local tangential and normal directions. The Reynolds numbers at x/c = 0.9 are Reτ ࣃ 328 and Reθ ࣃ 2, 255 on the suction side, and Reτ ࣃ 317 and Reθ ࣃ 785 on the pressure side. Statistical analysis of the turbulent boundary layer at two different stations on the upper and lower side of the wing were compared with the DNS of ZPG turbulent boundary layer from Schlatter and Örlü (2010) at matching friction Reynolds numbers. The effect of the APG on the mean flow is consistent with previous observations (Monty et al., 2011; Vinuesa et al., 2014), namely a steeper incipient log law, a more prominent wake region and a lower Cf . This is connected with the fact that the APG decelerates the flow and thickens the boundary layer, therefore for similar Reτ values the APG boundary layer would carry more momentum than the corresponding ZPG one. This is also associated with larger shape factors, which relate displacement thickness of the boundary layer with the corresponding momentum thickness. It is also interesting to note that the effect of the APG at a

Please cite this article as: S.M. Hosseini et al., Direct numerical simulation of the flow around a wing section at the moderate Reynolds number, International Journal of Heat and Fluid Flow (2016), http://dx.doi.org/10.1016/j.ijheatfluidflow.2016.02.001

JID: HFF 12

ARTICLE IN PRESS

[m5G;March 26, 2016;19:24]

S.M. Hosseini et al. / International Journal of Heat and Fluid Flow 000 (2016) 1–12

moderate value of β ࣃ 2.4 is reflected all the way down to the buffer layer which suggests that there is a different momentum transport mechanism across the boundary layer in the presence of APGs. We also assess the effect of the APG on the TKE profiles, which show a progressively larger inner peak value for increasing β , as well as the emergence and development of an outer peak when stronger APGs are present. These findings are consistent with the observations of Monty et al. (2011) in the streamwise turbulence intensity, and show the interaction of the APG with the large-scale structures of the flow, which become stronger as β increases. As a final remark, we would like to highlight the fact that the boundary layers developing around the wing section are influenced by the particular angle of attack of 5°, which basically determines the characteristics of the flow around the wing and the particular pressure-gradient distributions. Different angles of attack will lead to substantially different β distributions and TBL characteristics. Moreover, the boundary layers under consideration are also strongly influenced by history effects, i.e., the particular evolution of β (x), and the curvature of the wing surface. We are currently investigating the influence of such effects on the boundary layers by performing numerical simulations on flat-plate boundary layers subjected to various pressure-gradient configurations (Bobke et al., 2016). The simulations were performed on resources provided by the Swedish National Infrastructure for Computing (SNIC) at the Center for Parallel Computers (PDC), in KTH, Stockholm. RV and PS acknowledge the funding provided by the Knut and Alice Wallenberg Foundation; PS acknowledges financial support from the Swedish Research Council (VR). References Abbot, I., von Doenhoff, A., 1959. Theory of Wing Sections, 2nd ed. Dover Publications. Alfredsson, P.H., Segalini, A., Örlü, R., 2011. A new scaling for the streamwise turbulence intensity in wall-bounded turbulent flows and what it tells us about the “outer” peak. Phys. Fluids 23, 041702. Bobke, A., Vinuesa, R., Örlü, R., Schlatter, P., 2016. Large-eddy simulations of adverse pressure gradient turbulent boundary layers. J. Phys. : Conf. Ser. (in press). Choi, H., Moin, P., 2012. Grid-point requirements for large eddy simulation: Chapman’s estimates revisited. Phys. Fluids 24 (1), 011702. Coles, D., Wadcock, A.J., 1979. Flying-hot-wire study of flow past an NACA4412 airfoil at maximum lift. AIAA J. 17, 321–329. El Khoury, G., Schlatter, P., Noorani, A., Fischer, P.F., Breethouwer, G., Johansson, A., 2013. Direct numerical simulation of turbulent pipe flow at moderately high Reynolds numbers. Flow Turb. Combust. 91 (3), 475–495. Eliasson, P., 2002. EDGE, a Navier-Stokes solver for unstructured grids. In: Kroner, D., Herbin, R. (Eds.), Proceedings to Finit Volumes for Complex Applications III. Hemre Penton Science London, pp. 527–534. Fischer, P., Mullen, J., 2001. Filter-based stabilization of spectral element methods. C. R. Acad. Sci. Paris 332 (3), 265–270.

Fischer, P. F., Lottes, J. W., Kerkemeier, S. G., 20 08. nek50 0 0 Web page. http://nek50 0 0.mcs.anl.gov. Gungor, A.G., Maciel, Y., Simens, M.P., Soria, J., 2014. Analysis of a turbulent boundary layer subjected to a strong adverse pressure gradientoundary layer subjected to a strong adverse pressure gradient. J. Phys: Conf. Ser. 506, 012007. Hastings, R.C., Williams, B.R., 1987. Studies of the flow field near a NACA 4412 aerofoil at nearly maximum lift. Aero. J. 91, 29. Jansen, K., 1996. Large-eddy simulation of flow around a NACA 4412 airfoil using unstructured grids. CTR. Annu. Res. Br. 225–232. Jeong, J., Hussain, F., 1995. On the identification of a vortex. J. Fluid Mech. 285, 69–94. Johnson, F.T., Tinoco, E.N., Yu, N.J., 2005. Thirty years of development and application of CFD at Boeing commercial airplanes. Comput. Fluids 34 (10), 1115–1151. Lee, J.H., Sung, J., 2008. Effects of an adverse pressure gradient on a turbulent boundary layer. Int. J. Heat Fluid Flow 29 (3), 568–578. Monty, J., Harun, Z., Marusic, I., 2011. A parametric study of adverse pressure gradient turbulent boundary layers. Int. J. Heat Fluid Flow 32 (3), 575–585. Nagib, H.M., Chahuan, K.A., 2008. Variations of von Kármán coefficient in canonical flows. Phys. Fluids 20, 101518. Patera, A.T., 1984. A spectral element method for fluid dynamics: laminar flow in a channel. J. Comput. Phys. 54, 468–488. Piomelli, U., Yuan, J., 2013. Numerical simulations of spatially developing, accelerating boundary layers. Phys. Fluids 25, 101304. Rodríguez, I., Lehmkuhl, O., Borrell, R., Oliva, A., 2013. Direct numerical simulation of a NACA0012 in full stall. Int. J. Heat Fluid Flow 43, 194–203. Schlatter, P., Örlü, R., 2010. Assessment of direct numerical simulation data of turbulent boundary layers. J. Fluid Mech. 659, 116–126. Schlatter, P., Örlü, R., 2012. Turbulent boundary layers at moderate reynolds numbers: inflow length and tripping effects. J. Fluid Mech. 710, 5–34. Shan, H., Jiang, L., Liu, C., 2005. Direct numerical simulation of flow separation around a NACA0012 airfoil. Comput. Fluids 34, 1096–1114. Sillero, J.A., Jiménez, J., Moser, R.D., 2013. One-point statistics for turbulent wall-bounded flows at Reynolds numbers up to δ + 20 0 0. Phys. Fluids 25, 105102. Skote, M., Henningson, D.S., Henkes, R.A.W.M., 1998. Direct numerical simulation of self-similar turbulent boundary layers in adverse pressure gradients. Flow Turb. Comb. 60, 47–85. Skåre, P.E., Krogstad, P.-r., 1994. A turbulent equilibrium boundary layer near separation. J. Fluid Mech. 272, 319–348. Slotnick, J., Khodadoust, A., Alonso, J., Darmofal, D., Gropp, W., Laurie, E., Marviplis, D., 2014. CFD Vision 2030 Study: A Path to Revolutionaty Computational Aerosciences. Technical Report 2014-218178. NASA. Spalart, P., Watmuff, H., 1993. Experimental and numerical study of a turbulent boundary layer. J. Fluid Mech. 249, 337–371. Tufo, H.M., Fischer, P.F., 2001. Fast Parallel Direct Solvers for Coarse Grid Problems. J. Parallel Distrib. Comput. 61 (2), 151–177. Vinuesa, R., Örlü, R., Schlatter, P., 2016. On determining characteristic length scales in pressure gradient turbulent boundary layers. J. Phys.: Conf. Ser. (in press). Vinuesa, R., Rozier, P.H., Schlatter, P., Nagib, H.M., 2014. Experiments and computations of localized pressure gradients with different history effects. AIAA J. 52, 368–384. Vinuesa, R., Schlatter, P., Malm, J., Mavriplis, C., Henningson, D.S., 2015. Direct numerical simulation of the flow around a wall-mounted square cylinder under various inflow conditions. J. Turbul. 16 (6), 555–587. Wadcock, A. J., 1987. Investigation of low-speed turbulent separated flow around airfoils, NACA CR 177450. Wallin, S., Johansson, A.V., 20 0 0. An explicit algebraic Reynolds stress model for incompressible and compressible turbulent flows. J. Fluid Mech. 403, 89–132. Wolf, W.R., Azevedo, J.L.F., Lele, S.K., 2012. Convective effects and the role of quadrupole sources for aerofoil aeroacoustics. J. Fluid Mech. 708, 502–538.

Please cite this article as: S.M. Hosseini et al., Direct numerical simulation of the flow around a wing section at the moderate Reynolds number, International Journal of Heat and Fluid Flow (2016), http://dx.doi.org/10.1016/j.ijheatfluidflow.2016.02.001