Aerodynamic effects of compressibility for wind turbines at high tip speeds

. In the present work two dimensional airfoil computations are used to investigate the eﬀects of compressibility in the tip region of large scale wind turbines of 20 MW+ size. In the past application of incompressible CFD solvers have been wide spread for wind turbine aerodynamics, due to their eﬃciency and robustness at the near incompressible conditions experienced near the rotor center. With the increasing size of modern wind turbines and the desire to approach high tip speeds, the incompressible assumption might be violated in the tip region of the turbine. To investigate the eﬀects of compressibility and the possibility of correcting incompressible ﬂow solutions using explicit compressibility corrections, a CFD study of 2D airfoil aerodynamics at conditions of a large scale wind turbine is performed. The present study show that classical compressibility corrections can be successfully applied as a post-processing step to incompressible solutions, reducing the error in the predicted lift and drag to within a few percent for attached ﬂow conditions where viscous eﬀects are limited at Mach numbers upto 0.3.


Introduction
For large wind turbines with rotor radii of the order of 100 meters, the blade tip might reach high relative velocities approaching to thirty percent of the speed of sound.Most engineering codes are relying on airfoil data taken at Mach numbers below 0.2 intended for incompressible conditions, and several of the advanced CFD codes applied to wind turbines are based on the assumption of incompressibility.It is therefore important to evaluate the effect of the Mach numbers that surpasses the value of 0.2 which is normally taken as the limit below which the incompressible assumption is valid.
The reason why incompressible CFD solvers are of interest when solving the wind turbine rotor aerodynamics, is tied to the large variation of the local Mach number as function of radius, see Figure 1.For Mach number substantially below 0.1 the numerical approach in many compressible CFD solvers might require small time-steps, a high number of sub-iterations or advanced pre-conditioning.The incompressible flow solvers are tailored for zero Mach number, and do not face these numerical issues, but are due to the incompressible assumption incapable of resolving the weak compressible effects observed at the tip connected to Mach numbers close to 0.3.Some studies of compressibility in connection with wind turbines have been performed in the past.The active use of shock induced stall as a passive means of over-speed protection for smaller turbines were investigated early on by Wood [24], using a Blade Element Momentum and empirical formulas for the Mach number effect on the airfoil data.In Hossain et al. [5], the S809 airfoil were studied at transonic conditions at a Mach number of 0.8, which is not especially relevant for the envisioned operation of large scale wind turbines.In the European INNWIND and AVATAR projects spanning the years 2013-2017, focus were put on the effects of compressibility for large scale turbines, see [19].Recently, a study by Yan and Archer [25] from 2018 also investigated the effects of compressibility using an Actuator Line (AL) formulation coupled to both an incompressible and a compressible solver for the NREL 5MW turbine [7].The airfoil data for the AL model were obtained from the incompressible data using the Prandtl-Glauert correction, [4].In their work they reported the necessity of correction for wind speeds above 15 m/s and tip speeds ratios above 12, which corresponds to Mach numbers above 0.5 which are on the high side of the present interest.
The main objective of the present study is to investigate the effects of compressibility on the airfoil aerodynamics of relevance to modern type large scale wind turbine rotors, with diameters above 200 meters.

Methodology
In the AVATAR project, see [19], it was illustrated that an incompressible CFD solution or airfoil polar can be partially corrected for the compressible effects at the tip, by applying classical corrections for compressibility as discussed by Anderson [1], Prandtl-Glauert [4], Laitone [10], or Kármán-Tsien [8] and [22].Unfortunately, the correction is done as a post-processing of the incompressible flow solution, and therefore there are no direct coupling to the viscous development of the boundary layer.In general, compressibility effects increase the adverse pressure gradient which causes stronger boundary layer thickening or premature separation compared to the incompressible case.The associated decambering effect opposes the liftincreasing impact of the classical post-processing pressure corrections.When evaluating the effect of compressibility in a full 3D rotor CFD simulation, the problem is obscured by the fact that the effect of compressibility will influence not only rotor load but also the rotor induction, making it difficult to directly compare the compressible and incompressible flow solutions.
To circumvent this induction issue a 2D investigation is initiated for an airfoil at conditions identical to the situation at the blade tip of the AVATAR rotor, with respect to Reynolds number and Mach number.Comparisons are performed between the compressible FLOWer solver, the incompressible EllipSys2D solver, and a special compressible version of the EllipSys2D solver.Using identical grids, and computational setup we aim at isolating the effects of compressibility, which are then compared with the classical compressibility corrections.

Numerical Methods
Two different CFD codes are used in the present version, the FLOWer code and the EllipSys2D code.
The EllipSys2D is an in-house incompressible finite volume Reynolds-Averaged Navier-Stokes (RANS) flow solver based on the pressure correction approach, see [15,13,14], [17].The convective terms are discretized using a second order upwind scheme, while the diffusive terms are discretized by second order central differencing.The turbulent simulations are carried out using Menter's k − ω SST model described in [11].While the original EllipSys2D is purely incompressible, the present work include an extended capacity following the approach of [23], [3] to account for compressible effects.
For incompressible simulations, the pressure is extrapolated at all external boundaries, inlets, outlets and walls.For the pressure correction, zero gradient is assumed at all external boundaries.For the remaining variables, velocities and turbulent quantities, Dirichlet conditions are used at the inlet parts of the domain, while a zero gradient von Neumann conditions is applied at the outlet.At the wall no-slip adiabatic conditions are applied.
For the compressible simulations, the inlet conditions are changed to prescribed total conditions (pressure and temperature), while the static pressure is prescribed at the outlet.No changes are made for the other variables, see [3].
The CFD code FLOWer is a compressible RANS solver originally developed at the German Aerospace Center(DLR) [9] in the late 1990s.Over the last years, it was continuously enhanced for wind turbine application at the Institute of Aerodynamics and Gas Dynamics, University of Stuttgart [16].FLOWer's finite-volume formulation is based on block-structured grids with an in most cases cell-centered scheme.Convective fluxes are computed by a second order central discretization with artificial damping, also called the Jameson-Schmidt-Turkel method [6].Dummy layers are used around each block in order to maintain the spatial order over the block boundaries.Time integration is accomplished by an explicit multi-stage Runge-Kutta scheme.In case of steady computations convergence can be accelerated by implicit residual smoothing, local time stepping and the multi-grid algorithm.Transient simulations make use of the semi-implicit dual-time-stepping scheme.Several turbulence models are available to close the equation system.In the present study, the model by [12] is used.The simulations were performed steady mode and if needed for converged loads, continued in unsteady mode.A vortex correction for the far-field boundary is applied in order to eliminate potential influences.

Grid generation
The computational grid is generated with the DTU in-house HypGrid2D program, using an o-mesh topology, see [18].The grid consist of 320 cells in the chordwise direction and 320 in the normal direction, see Figure 2 and 3.The cell size at the trailing edge is equal to 1/4000 times the chord and 1/500 times the chord at the leading edge.The grid is cluster towards the wall in the normal direction, having a cell height of 5 × 10 −7 of the chord at the wall, assuring

Studied cases
The present study looks at three different scenarios.First we compute the actual lift polar corresponding to the AVATAR rotor setup at a realistic Mach and Reynolds number, along with incompressible simulations at an identical Reynolds number, see [2].Secondly, we look at the effect of varying the Mach number from 0 to 0.5 at selected angles of attack.Finally, a situation with varying Mach numbers at selected Reynolds numbers at a fixed AoA of 6 degrees is investigated.
The basic case is based on the DU91-W2-250 airfoil [21] scaled to have a thickness of 24 percent, according to the setup for the r/R=0.97section of the AVATAR rotor.The operational conditions for the airfoil study is derived from the AVATAR operational conditions, having an airfoil of 2 The second investigation, of the Mach number effects at the original Reynolds number of 14 millions, is performed by varying Mach number through the velocity and adjusting the viscosity to keep the Reynolds number constant.
The final investigation is obtained by simultaneously varying the viscosity and velocity to obtain the desired Reynolds and Mach number for a fixed AoA of 6 degrees.

Polar computations for the AVATAR rotor setup
Comparing the three solvers, FLOWer, EllipSys2D and EllipSys2D(comp) with respect to lift and drag for the basic AVATAR setup, the results are shown in Figures 4 and 5. From previous investigation we would expect a variation due to the different solvers to be less than 0.5% in lift and 2-3% in drag using identical grids, see [20].The origin of these differences are caused e.g. by differencing schemes and boundary condition implementation.From the Prandtl-Glauert correction we should expect a variation of the order of ∼ 5%.Based on previous experiences we should expect to be able to see this type of differences in the lift where the compressibility correction is 10 times larger than the expected solver differences, while the difference in drag is of similar order to the expected agreement between the solutions.
Examining the lift in Figure 4, we observe that in the linear region between [-8:8] degrees there is an excellent agreement between the two compressible solvers while the incompressible solution predicts a lower slope in accordance with expectations based on the Prandtl-Glauert correction.When reaching into the onset of separation both at high positive and negative values, the two compressible solver start to deviate.
For the drag polars the picture is not consistent, from Figure 5 it is observed that the drag predicted by the compressible FLOWer solver, is in agreement with the drag predicted by the incompressible EllipSys2D simulations.As the difference between the incompressible EllipSys2D and compressible EllipSys2D is around the expected 4.8 percent, the missing agreement between the two compressible solvers can be explained by the expected error band from previous investigations.Looking at the viscous part of the drag also shown in Figure 5, it is obvious that the main difference in the linear region of the drag polar is originating in the inviscid part of the drag.To avoid effects of code differences, the two versions of the EllipSys2D running with the exact same discritization are further investigated.Applying the simple Prandtl-Glauert correction to the lift and drag: , and we observe excellent agreement in the linear region both for the lift and the drag, see Figure 6 and Figure 7. Outside of the linear region, the viscous effects have a larger effect and the Prandtl-Glauert correction is no longer sufficient as should be expected based on its origin in potential thin airfoil theory.As a consequence the corrected incompressible solution does not predict the correct stalling behavior seen in the fully compressible solution.We must expect that applying the correction as a post-processing of incompressible solutions will only be valid in the attached region.

Effects of varying the Mach number at fixed Reynolds number
To investigate in more detail the accuracy of the compressibility correction under influence of viscous effects, a study is performed at a fixed Reynolds number at selected angles of attack   (0, 2, 4, 6) degrees, for Mach numbers varying between 0 and 0.5 with increments of 0.05.In Figure 8 and 9 the difference between the incompressible lift and the actual fully compressible lift is shown.From these figures the smallest error is observed at the lowest angles of attack in both lift and drag.Additionally, we observe that the error increase with increasing Mach number, and reaches a max error of -14% in lift and -18% percent in drag at 6 degrees AoA and a Mach number of 0.5.Applying the Prandtl-Glauert correction as a post-processing to the incompressible solutions and comparing with the compressible solutions, the results are shown in Figure 10 and 11.From these figures it is obvious that at moderate angles of attack a considerable improvement in the agreement is obtained by applying the compressibility correction, lowering the maximum error in lift from -14% to 6% and for drag from -18% to -6%.Generally, the correction is over correcting at low angles of attack and gives to low correction at higher angles of attack in both lift and drag.
Applying the Prandtl-Glauert compressibility correction directly to the incompressible pressure distribution an improved agreement with the actual compressible solution can be obtained, see Figure 12     include local Mach effects, the Kármán-Tsien correction provides superior agreement compared to the other two correction approaches for the present case.The downside of both the Kármán-Tsien and the Laitone corrections is that they can not be applied to the integral quantity but must be applied directly to the pressure distribution.
Looking at the viscous effects for varying Mach numbers, it is observed that the viscous effects of compressibility is very limited at Mach numbers up-to 0.3, see Figure 16 and 17.For a Mach number of 0.5 the effect of compressibility on the skin friction is cleary visible.The effect of compressibility at high Mach numbers are also clearly seen in the velocity profiles in the boundary layer at the trailing edge part of the airfoil, see Figure 3.2.For the lower Mach numbers the change in drag must thus be attributed to the change in the pressure drag.

Effects of varying the Mach number at different Reynolds numbers
Finally, the effect of the viscosity are illustrated at fixed Reynolds numbers [1,5,10,15] millions and an angle of attack of 6 degrees, by varying the Mach number from 0 to 0.5, see Figure 19 and 20.At the highest Mach number the lowest Reynolds number clearly gives the highest deviation for both lift and drag.At Mach numbers below 0.15 where the difference between the Prandtl-Glauert corrected incompressible lift and the actual compressible lift are below 0.25%, the picture is not consistent and here the lowest Reynolds number causes the smallest error.For the drag, generally the picture is not consistent, but the error in drag is quite low, and comparable to the prediction accuracy in most studies.

Discussion and Conclusions
In the present work a 2D study of the effects of compressibility is performed corresponding to the operational conditions experienced at the tip of a 20MW wind turbine having a tip speed of 100 [m/s].The initial 2D polar simulations showed good agreement between the two applied compressible solvers with respect to lift.Additionally, it is shown that the lift predicted by a 2D incompressible simulations can be brought to agree with the compressible simulations within the linear region using a simple compressibility correction.In agreement with expectation, the deviation in drag due to compressibility is comparable with the difference observed in the results from changing between solvers.showed that the standard compressibility correction could correct the lift and drag to within 2.5 percent for Mach number upto 0.3.Additionally, the final Reynolds number investigation shows that the high Reynolds numbers typically experienced in the tip region are beneficial with respect to limiting the viscous effects that deteriorates the accuracy of the compressibility correction.
The investigations of the pressure and skin friction distributions clearly indicated that the main effect with respect to compressibility originates from the change to the pressure distribution, while the viscous effects are very limited at moderate AoA's.The decreasing accuracy of the compressibility correction outside the linear region of the lift curve is clearly connected to viscous effects, as observed by the large changes in skin friction at high Mach numbers for the 6 degrees AOA.Additionally, the velocity profile at 65% chord also indicate that the change to the velocity profile is very limited at Mach 0.3 while at Mach 0.5 there is a substantial change.
The present investigation supports that incompressible flow solvers together with explicit 1234567890 ''"" The Science of Making from Wind (TORQUE 2018) IOP Publishing IOP Conf.Series: Journal of Physics: Conf.Series 1037 (2018) 022003 doi :10.1088/1742-6596/1037/2/022003 compressibility correction can be a viable option for rotor simulations.As the angles of attack is low in the tip region during normal operation, the standard compressibility correction could with advantages be used in connection with load and power curve computations when the tip velocity is high.For high AoA's, which are typically not present at the tip during normal operation, the explicit compressibility corrections are insufficient and full compressible formulations must be applied.

Figure 1 .
Figure 1.Mach number variation with radius for the AVATAR turbine at selected wind speeds excluding the effect of induction.
1234567890 ''"" The Science of Making Torque from Wind (TORQUE 2018) IOP Publishing IOP Conf.Series: Journal of Physics: Conf.Series 1037 (2018) 022003 doi :10.1088/1742-6596/1037/2/022003 1234567890 ''""The Science of Making Torque from Wind (TORQUE 2018) IOP Publishing IOP Conf.Series: Journal of Physics: Conf.Series 1037 (2018) 022003 doi :10.1088/1742-6596/1037/2/022003 y + well below one at the high Reynolds numbers investigated.The far-field boundary is placed 50 chords away from the airfoil surface to minimize effects of the far-field boundary conditions.The far-field boundary is specified as an inlet except for the +/-45 degrees downstream of the airfoil, which is specified as an outlet.

Figure 2 .
Figure 2. Total view of the computational grid.

Figure 3 .
Figure 3. Detail of the grid close to the airfoil surface.
[m] chord, a relative velocity of 102 [m/s], a temperature of 288.15 [K], a static pressure of 101300 [Pa], a viscosity equal to 1.7879 × 10 −5 [kgm −1 s −1 ]and a density of 1.22467 [kg/m −3 ].The polar is simply obtained by varying the angle of attack.The corresponding Mach number is 0.3 and the Reynolds number is 14 millions.

Figure 4 .
Figure 4. Computed lift polar with the three different formulations.

Figure 5 .
Figure 5. Computed drag polar with the three different formulations.

Figure 8 .
Figure 8. Difference between incompressible lift and actual lift in percent of actual lift.

Figure 9 .
Figure 9. Difference between incompressible drag and actual drag in percent of actual drag.

Figure 10 .
Figure 10.Difference between Prandtl-Glauert corrected incompressible lift and actual lift in percent of actual lift.

Figure 11 .
Figure 11.Difference between Prandtl-Glauert corrected incompressible drag and actual drag in percent of actual drag.

Figure 17 .
Figure 17.Skin friction distribution at three Mach numbers for 6 degrees angle of attack.

Figure 19 .
Figure 19.Difference between Prandtl-Glauert corrected incompressible lift and actual lift in percent of actual lift at 6 degrees AOA for varying Reynolds numbers.

Figure 20 .
Figure 20.Difference between Prandtl-Glauert corrected incompressible drag and actual drag in percent of actual drag at 6 degrees AOA for varying Reynolds numbers.
Figure 18.Wall normal profiles of normalized velocity profiles at x/chord=0.65 for three different Mach numbers at 6 degrees AoA.
The study at varying Mach numbers at fixed AoA's and a fixed Reynolds number clearly