Study of laminar separation bubble on low Reynolds number operating airfoils: RANS modelling by means of an high-accuracy solver and experimental verification

This work is devoted to the Computational Fluid-Dynamics (CFD) simulation of laminar separation bubble (LSB) on low Reynolds number operating airfoils. This phenomenon is of large interest in several fields, such as wind energy, and it is characterised by slow recirculating flow at an almost constant pressure. Presently Reynolds Averaged Navier-Stokes (RANS) methods, due to their limited computational requests, are the more efficient and feasible CFD simulation tool for complex engineering applications involving LSBs. However adopting RANS methods for LSB prediction is very challenging since widely used models assume a fully turbulent regime. For this reason several transitional models for RANS equations based on further Partial Differential Equations (PDE) have been recently introduced in literature. Nevertheless in some cases they show questionable results. In this work RANS equations and the standard Spalart-Allmaras (SA) turbulence model are used to deal with LSB problems obtaining promising results. This innovative result is related to: (i) a particular behaviour of the SA equation; (ii) a particular implementation of SA equation; (iii) the use of a high-order discontinuous Galerkin (DG) solver. The effectiveness of the proposed approach is tested on different airfoils at several angles of attack and Reynolds numbers. Numerical results were verified with both experimental measurements performed at the open circuit subsonic wind tunnel of Università Politecnica delle Marche (UNIVPM) and literature data.


Introduction
In this paper an innovative approach for the study of low Reynolds numbers operating airfoils involved in the blade sections of wind turbines is presented. In particular, great efforts are addressed to Laminar Separation Bubble (LSB) prediction, which is particularly difficult to be accurately computed without using large amount of computational resources, those required by the Direct Numerical Simulation (DNS) or by the Large Eddy Simulations (LES) approaches, as well known in literature.
The LSB, illustrated in figure 1, is a local boundary layer separation phenomenon that may occur on aerodynamic bodies operating at Re ≤ 10 6 . It can occur in presence of some Figure 1. Laminar Separation Bubble features, [1]. conditions that are briefly described: a laminar separation of the laminar boundary layer due to an adverse pressure gradient; a turbulent transition within the separated shear layer; a turbulent reattachment. Under these conditions a region characterised by a slow recirculating flow and by an almost constant pressure is formed. The presence of a LSB may introduce two classes of problems: (i) a decrease in the efficiency of the airfoil, mainly due to the airfoil drag rise; (ii) the appearance of large pressure fluctuations in the case of bubble bursting. This complex phenomenon is a challenging task of aerodynamics and it has already been widely studied by several authors with both experimental [2,3,4,5,6,7,8,9,10,11] and numerical techniques [12,13,14,15,16,17]. However, a computational efficient prediction of the LSBs still remains a difficult task. From the engineering point of view, thanks to their limited computational requests, the main simulation technique feasible to deal with complex configuration, such as wind turbines blades, appears to be the Reynolds Averaged Navier-Stokes (RANS) class of methods. A critical point in applying the RANS approach to low-Reynolds number flows is the presence of laminar regions. In fact in these cases the flow is laminar at least up to the separation point and the transition is difficult to be accurately detected. Widely used turbulence models are designed only for fully turbulent flows, thus they cannot be easily adopted in this context. However in this work we have used an incompressible high-order discontinuous Galerkin (DG) RANS solver, described in [18,19] and based on the standard Spalart-Allmaras (SA) model, to simulate the occurrence of LSBs on airfoils. As recognised by Rumsey [20] the SA model exhibits an apparent transition behaviour if the boundary free-stream condition forν is sufficiently small. In our numerical results, reported in [18], the artificial transition appears even if the trip term, considered by Rumsey the main responsible for this undesirable behaviour, is not implemented. Furthermore we have also verified that accurate numerical solutions, either on refined grids or with higher polynomial orders, produce larger laminar regions. Our approach for the LSBs prediction uses the previous observations, in other words we adopt the standard SA model with the ideal free-stream conditionν ∞ = 0, beingν ∞ the freestream turbulent viscosity, without any modification or specific equations to model the turbulent transition. If a fine enough resolution is guaranteed the artificial transition is delayed. Hence the computed solution is induced to remain laminar up to the separation point. Downstream the flow separation the production of vorticity activates the SA production term and in this way the turbulent transition is modelled. This behaviour allows, as it will be clear below, to capture the LSB phenomena almost correctly. It is particularly important to remark that the use of the ideal free stream condition,ν ∞ = 0, is not trivial, due to the numerical stability problems induced by the condition itself. Our computations can benefit from this condition thanks to our robust SA model implementation characterised by a non-standard r closure function definition, see §2.
The numerical results obtained with our DG solver were verified with the experimental data available in literature and with the measurements performed at the open circuit subsonic wind tunnel of Università Politecnica delle Marche (UNIVPM). More in depth, two airfoils were investigated: (i) the Eppler 387 (E387) airfoil, to validate both the numerical code and the experimental facilities with literature data; (ii) the WT2 airfoil. The later airfoil was designed at UNIVPM as a medium blade airfoil for a small wind turbine, having a nominal power of 3.5 kW, for possible installations in green docks. The airfoil maximum thickness is 0.139c, placed at 0.318c, while the maximum camber is equal to 0.03c and is located at 0.448c, being c the chord length.

Governing equations
The complete set of the incompressible RANS-SA Partial Differential Equations (PDEs) system can be written as where u is the velocity vector, p = P/ρ is pressure divided by the density, ν is the kinematic viscosity and ν t is the turbulent viscosity, which is computed according to theν variable as The SA equation diffusion coefficient ξ and the source term s will be defined below. Even if the SA model was developed only for positiveν values, its numerical solution can admit un-physical negative values which often produce numerical instabilities and the blow-up of the computation. For this reason the equation forν is manipulated with the same approach used for the turbulent kinetic energy in the k − ω turbulence model implementation proposed by Bassi et al. [21]. Whenever negativeν values are reached the source terms are set to zero, while the diffusion coefficient is fixed to the kinematic viscosity value, this means where χ =ν/ν is the non dimensional turbulent variable.
Note that with positive turbulent viscosity the standard SA equation is recovered. To completely define the PDEs system, equation (1), the following closure functions are introduced.
where d is the minimum distance from wall, r max is constant positive value andS is a function of both the vorticity magnitude, S, and of theν variablẽ Finally, in order to completely define the SA model, the standard closure constants are introduced The resulting turbulence model is the standard SA model without the trip term, that in the standard nomenclature is referred to as f t1 ∆u 2 . Moreover also the ad-hoc numerical fix term is missing, f t2 c b1 Sν − (ν/kd) 2 . This term was originally introduced in order to avoid the turbulence transition upstream from where the numerical trip was set. However, it has the drawback of promoting an undesired apparent transition behavior if it is run in a f t simulation, [20]. This model is usually classified as fully turbulent but, as it will be clear below, it can effectively perform even in laminar regions. Note that in our high-order SA model implementation the r closure function, equation (5), is modified by introducing a kind of realizable condition proposed in [18]. The new definition consists in limiting r to a positive value, here r max = 10, and handles negative and large values in the same way because the original r, for a positive χ, is negative only if it goes through a singular point with r → ±∞.

Numerical Discretization
RANS equations and Spalart-Allmaras turbulence model were solved, as introduced in §1, by using the DG method, which is a finite-element method of Galerkin characterized by discontinuous solutions across the cell boundaries. During the last years, DG methods have received increasing interest from the Computational Fluid-Dynamics (CFD) community because they combine flexibility, stability, robustness and accuracy. As in continuous Galerkin methods, DG methods can increase their order of accuracy by rising the degree k of the piecewise polynomial approximation P k describing the solution. Moreover, the resulting high-order accuracy is retained even using unstructured and anisotropic hybrid meshes usually adopted in real-life applications. The discontinuous approximation between neighbouring elements is treated, as in finite volume methods, using numerical fluxes. In these methods the removal of the continuity constraint enhances the compactness of the scheme, as elements share information only with their closest neighbours through surface flux terms, establishing the implicit time integration to block matrices of lower density. This compactness of the computational stencil also implies a straightforward parallelization on clusters of CPUs.

Implementation issues
The DG code here employed is able to solve the incompressible RANS equation up to the seventh order in space and up to the fourth order in time. A peculiarity of our implementation is related to the shape finite-element functions, which are defined in the physical space rather than in the reference one. Another distinctive trait of our solver is connected to the convective flux evaluation, which is defined using the artificial compressibility flux introduced in [22]. In this scheme the incompressibility constraint is relaxed by adding an artificial compressibility, but only at the interface flux level, in order to recover the hyperbolicity of the equations. This algorithm allows, as in the compressible case, the definition of a Riemann problem for the inviscid numerical flux evaluation. One advantage of this method, differently by the standard artificial compressibilty algorithm, is that the resulting scheme is time consistent, see [23] for two-dimensional and [24] for three-dimensional unsteady applications. Lastly the viscous fluxes are computed using the well established BR2 scheme, [25]. The space discretized non-linear RANS equations are here solved with a pseudo time-stepping strategy based on a fully-coupled backward Euler implicit scheme. The linear system derived from the DG space discretization was solved using the PETSc library [26] adopting its restarted GMRES preconditioned algorithm in both the matrix-based and matrix-free versions [27]. All the numerical solutions here reported were obtained in parallel on a small Linux Cluster, with 8 AMD Opteron based nodes for a total of 64 CPU cores operating at 2.3 GHz.

Experimental Apparatus
Aerodynamic tests were carried out in the open circuit subsonic wind tunnel available at the DIISM Department of UNIVPM and having the following features: Flow velocity may be varied from 0 to 35 m/s by controlling the fan motor RPMs with an AC inverter. Several measurements, periodically performed by means of a calibrated hot wire anemometer, showed a mean turbulence level lesser than 0.3% for a section area greater than 90%. Both wing sections (the one equipped with the E387 airfoil and the one with the WT2 airfoil) were realised by joining the two moulds CNC built: in this way it was possible to reduce the uncertainties related to the manufacturing process, especially for the airfoil nose. Before assembling the two moulds, a complete set of pressure taps were realised in the middle section of the wing section, both on the airfoil extrados and intrados. Each tap was equipped with a plenum chamber, flush mounted below the wing surface, so to reduce the dynamic pressure component and obtain a more reliable static pressure reading. Velocity and turbulence measurements were carried out for each experimental test by a Dantec CTA anemometric system equipped with a pneumatic calibration unit; this is used before every measurement session in order to reduce the velocity overall error, estimated as 0.08 m/s. Pressure measurements were carried out with a Druck LPM 9481 differential pressure transducer having an operating range between 10 Pa and 1000 Pa and a typical error less than 0.1% FS BSL. In order to obtain the pressure coefficient distribution on the airfoil, this device was sequentially connected to all the pressure taps drilled in the middle section of the wing.

Eppler 387 airfoil
The Eppler 387 airfoil was extensively investigated at Re = 10 5 and Re = 3 · 10 5 at several angles of attack in order to verify the accuracy of both the high-order CFD solver and the experimental facilities available at UNIVPM. Reference values for this test case are: pressure and force coefficient measurements, performed at NASA Langley Research Center, [28], and force coefficients measurements obtained at University of Illinos, [29]. The numerical solutions were computed, up to P 5 polynomial approximation, on two grids. The first one has a C-type topology and it is named G1 in the following. It consists of 2880 computational elements with a piecewise cubic representation of the edge elements. More in depth this grid was obtained by agglomerating 8 elements of multi-blocks structured finer grid. The second grid, named G2 and characterised by 9239 elements and by a circular domain shape, is unstructured, triangular elements, with a boundary layer refinement consisting of quadrilaterals. It was generated using a code developed by Ghidoni et al. [30]. In this case the edges of the elements have a piecewise quadratic representation. Both G1 and G2 grids are able to guarantee O(y + c ) = 1, where y + c is the viscous sub-layer scaled non-dimensional first centroid distance of the elements next to the wall. Both the grids here used are represented in figure 2. At Re = 10 5 our numerical simulations, performed using both the G1 and the G2 grids, and the experimental measurements are in good agreement with the reference data, figure 3(a). It is worth mentioning that the numerical computations do not exhibit completely satisfactory results in the stall region due to RANS modelling weakness for large separated flows. At α = 0 • , figure 3(b), the laminar separation point, predicted by numerical computations, is located upstream with respect to the experimental data. At α = 4 • , figure 3(c), the UNIVPM numerical and experimental data agree pretty well, however in the NASA reference data the LSB is shifted downstream. In figure 3(d) the results at α = 8 • are reported. Numerical and reference data are in very good agreement, while our experimental results do not show evidence of LSB. Nevertheless, as already showed in [31], our wind tunnel can be considered reliable only up to the angles prior to the stall.
The flow field at Re = 3 · 10 5 has been up to now computed using only the G2 grid. Figure 4 shows the comparison between the numerical results and the reference one and they clearly agree very well. These data highlight what we have observed many times during our computational experience about this specific topic: the LSB prediction capabilities of our approach increase with the Reynolds number and seem to be "optimal" at Re = 2 · 10 5 ÷ 3 · 10 5 . However, it is important to mention that at α = 8 • there is no experimental evidence of LSB and that the related simulation was performed adopting a small free-streamν ∞ value, 10 −3 ν, otherwise the solver did not converge to a steady state solution. This behaviour seems to establish in some way a further link between the laminar separation and the transitional capabilities, unintentionally inserted in the PDE of the SA model. If the laminar separation does not take place in the flow field and if the apparent transition of the SA model does not move up the natural transition location, in this case experimentally detected at 0.4c, serious convergence problems of the solver are observed. This circumstance apparently prevents to compute un-physical solutions signalling the need of a larger free-streamν ∞ boundary condition to encourage the numerical transition in the attached boundary layer. However, this subject still requires a more accurate investigation.

WT2 airfoil
A second airfoil was here considered in order to fully prove the almost surprising, since the SA model was not designed for this purpose, effectiveness in handling LSB of the proposed approach. This profile, called WT2, was designed for wind turbine applications and tested at UNIVPM. Numerical computations were performed, at Re = 2 · 10 5 , up to P 5 polynomial approximation, see figure 5 for a representation of the flow-field, while in figure 6 the solution behaviour is depicted close and inside the LSB. A C-type structured grid with 2688 elements and a piecewise cubic representation of the edges was used for these computations. Obviously the computational elements are clustered near the airfoil surface in order to accurately solve the boundary layer and, even in this case, the first cell height ensures a y + value of about 1. Another time the high-order grid was obtained by agglomerating 32 elements of multi-blocks structured linear grid. Comparisons between numerical and experimental pressure distribution at α = −2 • , 1 • , 3 • and 5 • are here reported in figure 7. Also in this case the results appear very convincing and seem to confirm the reliability of the proposed approach for the simulation of the LSB phenomena at moderate Reynolds numbers. Note that at α = −2 • laminar separation bubbles are correctly detected at both the pressure and the suction sides.

Conclusions
This work assesses the reliability the results of a high-order CFD solver for LSBs simulations. The code is based on a DG method and implements, as described in [18], the SA equation for turbulence modelling. The RANS-SA system is employed in the fully-turbulent mode using the ideal free-stream condition,ν ∞ = 0. In this way the flow behaves as a laminar at least up to the separation point. However, the separated flow exhibits high level of vorticity activating the source terms of the SA equation and producing significative values of turbulent viscosity. This behaviour allows to promote the turbulent transition, without introducing a specific model for this purpose, inducing in the here investigated cases the turbulent boundary layer reattachment. The numerical approach was tested on two airfoils: the E387 and the WT2. In the first case the CFD solver and the wind tunnel facility, available at UNIVPM, were tested with respect to literature results. Lastly the airfoil for small wind turbines WT2, developed within UNIVPM,  was tested both numerically and experimentally. Very promising results were obtained for both the profiles, especially for the Re = 2 · 10 5 and Re = 3 · 10 5 flow regimes. Moreover note that the main object of this work was to develop a computational efficient simulation tool suitable for the study and the design of airfoils and to validate its effectiveness. Indeed the aim was not to obtain a "perfect" numerical solution for LSB problems, in this occurrence the use of more computationally demanding approaches like DNS or LES are almost unavoidable. Our goal was, in our opinion, completely reached since each of the computations here reported required no more than 8 CPU cores, nowadays available in almost all workstations, and only few hours of wall clock time to achieve a converged, up to machine precision, P 5 solution which is clearly even uselessly accurate. In fact, the most of the cases could be computed, with no evident changes in the obtained results and with a clear advantage in the CPU time, only  up to the P 2 or P 3 orders of approximation. A paper with a more in-depth numerical and theoretical analysis of the SA behaviour related to the laminar separation inducing the turbulent transition is under preparation, while we are actually considering the application of our model to three-dimensional problems.