Transition at Low-Re Numbers for some Airfoils at High Subsonic Mach Numbers

‡High altitude long endurance UAV flight regime imposes certain difficult testing conditions for the ground based test facilities, such as low density, low freestream turbulence, high alpha, low-Reynolds, and highsubsonic-Mach numbers. High altitude flight testing would be required to collect actual experimental data at an expense of much higher costs. Computational approach is a viable alternative to support generation of design data base for the laminar and transitional boundary layers at high-subsonic-Mach-numbers. In the present study, a two-equation γ-Reθt correlation-based transition model is further developed to predict some airfoils that are frequently used in the design of UAVs. Firstly, the empirical correlations already validated for low to high Mach number flat plate cases is validated for the thin NACA64A006 airfoil at low subsonic speed and high-Reynolds-number. Secondly, the present methodology is successfully demonstrated for the E387 and SD7037 moderately thick UAV type airfoils in the low-Reynolds and low-subsonic-Mach-number conditions. Finally, the relatively thicker APEX-16 airfoils at high-altitude, low-Reynolds-number conditions for high-subsonic Mach numbers are simulated. Results are compared with the available numerical data in the literature obtained through Reynolds Averaged Navier-Stokes and viscous-inviscid interaction methods using the e N method.


Introduction
High altitude long endurance UAVs with the role of being low cost satellites due to their long on station times bring renewed interest in Low Reynolds Number Airfoils. These airfoils' behaviors are quite different from their high Reynolds number counterparts [1]. Transition to turbulence and separation bubbles play important roles for these airfoils from low to high Mach numbers. Besides, compressibility makes the stability and transition problems more complex and realistic that are encountered in the transonic regime. As for the density effects, air gets thinner and the Reynolds number starts decreasing as the altitude increases [2]. For instance, at 10,000 m, despite the reduction of gravity of 0.3%, the reduction of air density from 1.225 kg/m 3 to 0.413 kg/m 3 is quite disadvantageous. At 21,000 m, the air density drops significantly to 0.0757 kg/m 3 . Meantime, the high altitude UAVs must either fly faster or increase the coefficient of lift to carry the weight. The net result is either increasing the airspeed within the power consumption limits and/or increasing the angle of attack as Re number decreases. Therefore, high altitude flight puts a lot of pressure on designers as to balance the power consumption, high angles of attack nearing stall angles, growing separation bubbles as Re number gets smaller and high subsonic Mach numbers adding the possibility of lambda shocks accelerating the instability in the bubbles. Furthermore, increasing the coefficient of lift moves the cruise point out of the so-called drag bucket limits where the drag increases very quickly. A good selection of a UAV airfoil must account for all these factors for the reason that a good airfoil at sea level may turn out to be a worse selection at high altitude and high alpha conditions. Predicting low-Reynolds number airfoil performance is a difficult task that requires correctly modeling several flow phenomena such as inviscid flow field with the presence of shock waves, laminar separation regions with presence of separation bubbles, transition to turbulence in the free shear layer and turbulent boundary layer. Especially the presence of the separation bubble may affect the results significantly. Constant pressure assumption across the bounday layer may not be valid across the bubble. Thus, correct modeling of the flow around the airfoils operating at low Reynolds numbers becomes a challenging research problem. Experimental study for low-Reynolds number flows also have some certain difficulties. In a low Reynolds number experimental work regarding the drag coefficient measurements on a Wortmann FX63-137 airfoil at the same test conditions at three different research facilities, it is reported that the results show differences of more than 50% [3]. High altitude long endurance vehicles' flight conditions impose difficult testing conditions for the ground based test facilities to emulate the low density, low freestream turbulence and high subsonic Mach numbers at the same time. Flight testing would be required to collect the actual experimental data. For instance, a high altitude sailplane project called APEX was started by NASA that aimed for collecting boundary layer data after being released from a very high altitude of 108 K ft. after a balloon launch at that altitude [4]. However, it was clear for the designers of the aircraft that computational methodology was very much required to improve our know-how on laminar and transitional boundary layers at high Mach numbers. In reference [4], a viscous-inviscid interaction methodology [5,6] and time-accurate RANS [7,8] were used in the numerical predictions. Today, state of the art Reynolds Averaged Navier-Stokes (RANS) solvers are widely available for numerically predicting fully turbulent part of flow fields, but none of these models are adequate to handle flows with significant transition effects because of lack of practical transition modeling. Nevertheless, transition predictions have shown certain progress and utility by means of the well-known e N method [5], some two-equation low Re-number turbulence models [9], and some methods based on experimental correlations [10].
Yet, great strides have been taken by the recent introduction of what is called as the "engineering transition modeling" by Menter et al. [11] that rely on local data to circumvent some complicated procedures in other methods. Further work demonstrated the viability and practicality of the transition correlation-based model that warranted further investigation. Since the present model [11] does not attempt to model the physics of the transition process, but rather to form a framework for implementation of transition correlations into general purpose CFD methods, different correlations for each user either remain propriatery [11] or user-dependent [12,13] based on the specific experimental data set. Therefore, those correlations that are used in the model are not universal, but reflect each user's own data base. A number of applications for the Menter et al [11] model were done by Genç [14] for a thin airfoil at high Reynolds numbers and Genç, Kaynak and Yapici [15] for flow control around an airfoil by jet blowing and suction at low Reynolds numbers. Following the Menter et al. [11] model, a number of successful two-or three-equation models also appeared in the literature such as k-k L -ω model of Walters and Leylek [16], near/freestream intermittency model of Lodefier et al. [17], and k-ω-γ model of Fu and Wang [18]. High speed applications of the Menter et al. [11] model was done by Kaynak [19] for flat plates up to supersonic Mach number of 2.7. Other high speed calculations were done by Fu and Wang [18] for supersonic flow past a straight cone and hypersonic flow over a flared cone at zero angle of attack.

Numerical method
In this study, 2-D computational results were obtained using the FLUENT software [20] for NACA64A006 and APEX 16 [4] airfoils. Although 3D effects are present, the aim of this study is to investigate the prediction of the boundary layer and stall characteristics of the airfoils. The k-ω SST turbulence model, k-ω SST transition model and k-k L -ω transition model are used in conjunction with the built-in RANS solver. The results were compared with the results gathered from the MSES code explained in the study of Drela et al. [4]. The MSES code is a viscous-inviscid interaction code that adapts the Euler equations coupled with the boundary layer code including the e N method for transition prediction. For both airfoils, free stream boundary conditions were used in the upstream, downstream and outer boundaries. Pressure-farfield boundary condition was used for both cases, with different Mach and Reynolds numbers.

Solution grid
C-Type structured solution grids were generated by the ICEM CFD software [20] for both NACA64A006 and APEX 16 airfoils as shown in Figures 2 and 3 respectively. The grid extends from -12 chords from upstream to 16 chords downstream and the upper and lower boundary extends 12 chords from the airfoil. Both grids include 61588 cells and 62058 nodes.

Turbulence and transition model
The k-ω SST turbulence model as implemented within the RANS equations blends the formulation of the Wilcox k-ω model [9] in the near-wall region with the formulation of the k-ε model [21] in the far-field developed by Menter [22]. The k-ω SST transition model [11] solves for four transport equations, such as the turbulent kinetic energy (k), specific turbulence dissipation rate (ω), intermittency (γ), and the transition onset momentum thickness Reynolds number (Re θt ) equations in addition to the basic RANS equations. The γ transport equation and Re θt transport equation are used to initiate the transition process and for establishing link with experimental correlations, respectively. The correlations are based on the free stream turbulence intensity (Tu) the Re θt at transition onset. The k-k L -ω model [16] solves three transport equations for the turbulent kinetic energy (k), laminar kinetic energy (k L ), and specific turbulence dissipation rate (ω) in addition to the RANS equations. The k L is based on the non-turbulent fluctuations in the laminar boundary layer, as defined in the work of Mayle and Schulz [23].  Figure 4 shows the numerical data obtained by using the FLUENT [20] showing experimental C L and C D values for different angle of attack values for the NACA64A006 airfoil. In this figure, all models give reasonably good results against the experiment in the linear region, but the calculation starts to differ from the experiment [24] after 5 o angle of attack. It is observed that after 5 o , the k-ω SST and k-ω-SST transition models underpredict the lift coefficient, whereas the k-k L -ω transition model predicts the lift coefficient better than those two models. On the other hand, Figure 4 shows that there is a quite good agreement between the experiment and computational results for the drag coefficient.   Figure 5 shows that there are significant differences between the experiment and numerical predictions in the pressure coefficient distribution. Comparing Figs. 4 and 5, although the k-k L -ω transition model appears to yield the best result for the lift coefficient, detailed pressure distributions do not wholly support this finding as there are large discrepancies for the local pressures.

APEX 16
Numerical results for APEX 16 airfoil are divided into two categories such as one high subsonic Mach number at a range of low Reynolds numbers and a range of low Reynolds numbers at range of high subsonic Mach numbers. All the numerical data for this airfoil is gathered using k-k L -ω transition model and they are compared with the numerical results obtained using the MSES code mentioned in the work of Drela et al [4].
i. High Subsonic Mach Number Case: In this section, a constant Mach number of 0.6 is selected and the simulations are made at Reynolds numbers of 200,000, 300,000 and 500,000. The aim of this investigation is to understand the effect of Reynolds number for the aerodynamic characteristics of the APEX www.intechopen.com Transition at Low-Re Numbers for some Airfoils at High Subsonic Mach Numbers 85 16 airfoil. Also in this section, prediction performances of the MSES code and k-k L -ω transition model are compared. Figure 6 shows the numerical data of C L values against C D values using the k-k L -ω transition model and MSES code for different Reynolds numbers. This figure shows that the MSES code and k-k L -ω transition model reasonably agree for the lift coefficient; but, there are some differences in the drag coefficient. Especially for low Reynolds numbers, the k-k L -ω transition model predicts smaller drag coefficients than the MSES code. As the Reynolds number increases, the agreement between two models start improving. For low Reynolds number flow conditions, predicting the drag coefficient becomes even more difficult, since low Reynolds number airfoils typically exhibit laminar separation bubbles, which are known to significantly affect the performance of an airfoil. As for increasing the Reynolds number, the ability of prediction of the drag and lift coefficient becomes easier for the flow solvers because the flow is no longer laminar, and turbulent boundary layer is effective on the surface of the airfoil beginning from the leading edge.
From Figure 6, it can be conluded that as the Reynolds number gets higher, lift coefficient gets higher whereas the drag coefficient gets lower. Although the predictions of MSES code and k-k L -ω transition model do not fully agree, the relative agreement is still reasonable as both models at least agree on the trend of the lift and drag coefficients as the Reynolds number changes. Taking Ma = 0.6 and Re = 300,000, further information on the characteristics of the APEX 16 airfoil is obtained using FLUENT's k-k L -ω transition model and the results are compared with the MSES code. For this condition, pitching moment predictions, lift coefficient predictions and transition location on the upper surface of the airfoil are compared. Figure 7 shows the lift coefficient predictions versus angle of attack comparison for the Ma = 0.6 and Re = 300,000 case. There is a good agreement for the lift coefficient between the kk L -ω transition model and MSES code predictions until the angle of attack reaches around 5 degrees; but after 5 degrees, the MSES code predicts lift coefficient lower than k-k L -ω transition model.  Figure 8 shows the pitching moment predictions for the same case. The pitching moment results obtained using FLUENT are all based on the pitching moment center at 25% of the chord length. The results show that, there is a good agreement between the two models in the angle of attack range from 0 to 4 degrees. However, the k-k L -ω transition model calculates a larger pitching moment for angles of attack larger than 4 degrees than the MSES code. Finally, Figure 9 shows the predictions of the location of the transition location on the upper surface of the airfoil against the lift coefficient. The transition locations obtained from the FLUENT results are all assumed that the transition occurs when turbulent to laminar viscosity ratio reaches about 80-100. Looking at this figure, the agreement between the two models are very good at low angles of attack. The reason why both curves are not matching perfectly is that, the lift coefficient predictions are harder to obtain at high angles of attack.

ii. High Subsonic Mach Numbers at a Range of Low Reynolds Numbers
In this section, low Reynolds numbers in the range (200,000-500,000) are kept constant, and Mach numbers are changed in the high subsonic range of 0.60, 0.65 and 0.70 at each Re number in order to understand the effect of changing Mach number.  Figure 10 shows the high subsonic Mach number data for the C L against C D values gathered from k-k L -ω transition model and the MSES code at Re=200,000. As seen in this figure, the drag coefficients calculated by the k-k L -ω transition model are much higher than the MSES code. Also, for Mach number 0.7, lift coefficient calculations of k-k L -ω transition model and the MSES code differs a lot. However, the k-k L -ω transition model which is based on Reynolds-averaged Navier-Stokes equations seems to produce smoother drag polars which more look like normal trend. On the contrary, the MSES code predicts smaller C D values as C L gets larger which should not be a normal pattern. Figure 11 shows the high subsonic Mach number data for the C L against C D values gathered from k-k L -ω transition model and MSES code at Re=300,000. As seen in the figure, the drag coefficients calculated by k-k L -ω transition model were again much higher than the MSES code. Also, for Mach number 0.7, lift coefficient calculations of k-k L -ω transition model and the MSES code differs a lot. Yet, as Reynolds number gets larger in this case, the shape of the drag polars slightly start normalizing as they look more like standard high Re number drag polars. Some kind of kink is apperant in the MSES simulations around C L =0.4.  Figure 12 shows the high subsonic Mach number data for the C L against C D values gathered from k-k L -ω transition model and the MSES code at Re=500,000. The drag coefficients calculated by the k-k L -ω transition model were again much higher than the MSES code. The shape of the drag polars quickly start normalizing as they look more like standard high Re number drag polars. The kink which is apperant in the MSES simulations around C L =0.4 is www.intechopen.com more pronounced for this Re number. From all these numerical analysis for the APEX 16 airfoil, it can be concluded that the drag coefficient decreases as Reynolds number increases for constant Mach numbers, and also drag coefficient increases as Mach number increases for the same Reynolds number. It appears that as the Mach number increases, the adverse pressure gradient causes the drag coefficient to increase while keeping the Reynolds number constant.

Re = 500,000 case
The separation bubble stability has an important effect on the airfoil predicted performance [4]. The MSES code, based on the stable bubble calculations, predicts a lift coefficient of 0.96 at the flight condition of Ma = 0.65, Re = 200,000 and an angle of attack of 4 degrees, whereas k-k L -ω transition model predicts an average section lift coefficient of 0.82 for the same flight condition. In the following, the shape, location and extent of the separation bubbles for different flow cases are introduced. Figure 13 shows the velocity contours and the velocity vectors of the APEX 16 airfoil at 0 o angle of attack, Ma = 0.6 and Re = 300,000 based on k-k L -ω transition model. From this figure, it is seen that the flow separation occurs at around 0.60c away from the leading edge. At this angle of attack, flow reattachment is not observed until the trailing edge. Figure 14   shows the velocity contours and the velocity vectors of the APEX 16 airfoil at 6 o angle of attack, Ma = 0.6 and Re = 300,000 based on k-k L -ω transition model. From this figure, it is seen that the flow separation occurs at around 0.10c away from the leading edge. Comparing this figure with Figure 13, it is concluded that the flow separation occurs closer to the leading edge as the angle of attack increases, as expected. Figure 15 shows the turbulent kinetic energy distributions for 0 o and 6 o angle of attack for the APEX 16 airfoil at the same flow conditions. This figure supports the observations and comments made for Figures 13 and 14. As seen in Figure 15, the flow transition into turbulence occurs closer to the leading edge as the angle of attack increases. Also, there is a laminar boundary layer separation first, followed by a shear layer transition into turbulence, and finally there is turbulent reattachment.

Conclusions
In this study, firstly, aerodynamic characteristics of the NACA64A006 airfoil is investigated using the k-ω SST turbulence model, k-ω SST transition model and k-k L -ω transition model using FLUENT. The results obtained from FLUENT are compared with the experiment [24], and it is observed that all numerical approaches give reasonably good results in the linear region, although the results began to differ as the angle of attack gets larger. Especially after 5 o of angle of attack, whereas the k- SST turbulence and k-ω SST transition models greatly underpredict the lift coefficient, the k-k L -ω transition model yields the best result. With regards to the drag coefficient, it is reasonable to say that all numerical methods agree quite well with the experimental data. For the pressure coefficient, it is observed that the k-k L -ω transition model also fares better than the other models.
In the second part of this study, aerodynamic characteristics of the APEX 16 airfoil is investigated using the k-k L -ω transition model and results are compared with the e N based results of the MSES code by Drela et al [4]. The first apparent characteristic for the APEX 16 is that increasing the Mach number results in a decrease in the maximum lift coefficient. Although the lift coefficient predictions of the MSES code and k-k L -ω transition model slightly differ, the general trends in both results are similar. As the Reynolds number decreases, the separation bubbles become larger, which is the reason for the increase in the drag coefficient. Another observation is that the slope of the lift curve is relatively unaffected by the Mach and Reynolds numbers except near the stall. The predicted transition location versus the lift coefficient is also presented for the Re = 300,000 and Ma = 0.6 case where both models agree on the transition location but the lift coefficient predictions differ at high angles of attack. It is shown that the transition location on the upper surface moves forward with the increasing angle of attack, as expected.