CFD code comparison for 2D airfoil flows

. The current paper presents the eﬀort, in the EU AVATAR project, to establish the necessary requirements to obtain consistent lift over drag ratios among seven CFD codes. The ﬂow around a 2D airfoil case is studied, for both transitional and fully turbulent conditions at Reynolds numbers of 3 × 10 6 and 15 × 10 6 . The necessary grid resolution, domain size, and iterative convergence criteria to have consistent results are discussed, and suggestions are given for best practice. For the fully turbulent results four out of seven codes provide consistent results. For the laminar-turbulent transitional results only three out of seven provided results, and the agreement is generally lower than for the fully turbulent case.


Introduction
In connection with the EU AVATAR project, dealing with the aerodynamics and aeroelasticity of the wind turbines of the future, a large effort is invested in the exploration of Computational Fluid Dynamics (CFD) codes.The continuous upscaling of wind turbines requires code validation for new flow regimes, where the Reynolds numbers increase above 10 × 10 6 , and the Mach numbers may be in the range [0.01:0.3].In this region, experimental data relevant to wind turbines are sparse and costly to obtain.In the EU AVATAR Project a three step procedure is used to solve this problem: • Establish the necessary requirements for several CFD solvers to provide consistent results.
• Validate the CFD solvers on the sparse data available.
• Apply the solvers to explore the target flow regime.
Finally, the results from the CFD solvers will be used to develop the engineering methods used in industry to handle the upscaled turbines.
The present work deals with a code comparison, that aims at establishing the requirement for the solvers to consistently predict both fully turbulent and transitional 2D airfoil flows.The focus is on conditions where the Mach number is 0.1 while the Reynolds number is varied from a moderate value of 3 × 10 6 to a high value of 15 × 10 6 .The conditions corresponds to a recent wind tunnel test of the DU00-W-212 airfoil, performed in the AVATAR project.Several of the codes have a long track record, and have previously been validated against data in several studies, both in connection with blind comparisons, and against data from the open literature, [29,30,28,10,23,2].
The main focus is to investigate the requirement of the codes to have consistent solutions for the target flow regimes, in terms of the necessary domain size, grid resolution, and iterative convergence criteria.Even though the issue of obtaining iterative and grid convergence of numerical methods is not new, it is still an essential problem, see [24,25,5].

Methods
Both compressible and incompressible Reynolds Averaged Navier-Stokes (RANS) methods are applied.All solvers apply the k − ω SST turbulence model by Menter [12], while different types of transition models are applied.The computations are carried out in time or iteration, until the variation over the last 10 percent of the integration time or iterations, are within 0.05 percent for the lift, and 0.1 percent for the drag.For most of the solvers, the L2 norm of the residuals are reduced by 5-6 orders of magnitude to assure this.
The CENER calculations are performed using the code WMB [6].The code is a structured multi-block cell-centred finite volume method, based on the compressible Unsteady Reynolds Averaged Navier-Stokes (URANS) equations.The calculations are done in steady state using the k − ω SST turbulence model for the fully turbulent calculations, and k − ω SST combined with the e N method for the transitional ones.
The CRES computations are performed using an incompressible implicit pressure correction solver.The algorithm introduces a matrix-free algorithm for pressure updating, which maintains the compatibility of the velocity and pressure field corrections, allowing for practically unlimited large time steps in steady state calculations.Spatial discretisation is performed on a computational domain, resulting from a body fitted coordinate transformation using finite volume schemes.The velocity is stored at grid nodes, while the pressure is computed at midcells.A linear 4 th order pressure dissipation term is added into the continuity equation to prevent velocity-pressure decoupling.Further details can be found in [3].For turbulent flow computation, the standard k − ω or the k − ω SST [13] models can be used.In the present paper all results are produced with the k − ω SST for fully turbulent flow.
The DTU computations are performed in steady state using the EllipSys2D in-house incompressible finite volume RANS flow solver, see [18,19], [26].The convective terms are discretised using the QUICK scheme, as given by [9].The turbulent simulations are carried out using Menter's k − ω SST model described in [12], while the transitional simulations are based on the e N transition model, as described in [20].
The NTUA computations are performed using the MaPFlow in-house compressible finite volume RANS flow solver.The implementation uses Roes scheme, low Mach preconditioning, and non-preconditioned far-field conditions in Riemann invariant form [22]. Menters k − ω SST turbulence closure model [12], and the e N transition model as described in [4] are used.In the present simulations limiters are suppressed.For angles of attack outside of the range [-11:11], unsteady simulations are performed.
The TUDelft calculations are performed using OpenFOAM v2.3.0, an open source CFD software distributed under the General Public Licence (GPL), see [7].OpenFoam is a segregated finite volume code able to solve compressible and incompressible flows.For this work, the steady, incompressible RANS equations are solved using the SIMPLE algorithm.The governing equations are solved using first order upwind discretisation schemes for the convective terms.Fully turbulent simulations use Menter's two-equation k-ω SST model [12].Transitional simulations are conducted using the kkl-ω model [35], which solves for the laminar and turbulent kinetic energy, as well as a third transport equation which governs low-frequency velocity fluctuations that trigger transition in the boundary layer.
The UOG computations use the Helicopter Multi-Block (HMB3) code [2], a cell-centred RANS finite volume approach.The solver uses a fully-implicit time integration.Following the pseudo-time approach, after the linearisation of the residual at the new pseudo time step, the discretised RANS result in a large system of linear equations.In HMB3, the latter is solved employing a Generalised Conjugate Gradient (GCG) method [1], pre-conditioned with an Incomplete Lower Upper (ILU) factorisation [11].To discretise the convective part of the Navier-Stokes equations, a formally third order Monotone Upstream-Centred Scheme for Conservation Laws (MUSCL) [33] is employed with the Van Albada limiter [34].In the present work, the Osher scheme [21] is used for the fluxes, and the turbulent closure is the k − ω SST model of Menter [13].When the laminar-turbulent transition is considered, the k − ω SST turbulence model is coupled with the γ−equation LCTM of Menter [17].A second order central discretisation scheme is used for the viscous fluxes.
The USTUTT simulations are performed with the CFD code FLOWer, which is developed by the German Aerospace Center (DLR), see [8].FLOWer is a compressible code that solves the RANS equations in integral form.The numerical scheme is based on a finite-volume formulation for block-structured grids.To determine the convective fluxes, a second order central discretisation with artificial damping is used, also called the Jameson-Schmidt-Turkel (JST) method.Time integration is accomplished by an explicit multi-stage scheme including local time-stepping.The time-accurate simulations make use of the implicit Dual-time-stepping method.In the present study, the simulations are performed in steady and unsteady mode if needed for converged loads.Turbulence modelling is accomplished by the Menter's k − ω SST model described in [12].

Grid Generation
A common grid is provided for all partners to serve as a reference point, generated by the hyperbolic grid generation code HypGrid2D [27] around the DU00-W-212 airfoil [32].The grid has a single block O-topology with a finite sharp trailing edge, and a total of 73728 cells.The grid has 384 cells in the chord-wise direction and 192 cells in the normal direction, placing the outer boundary 40 chords away from the airfoil.The cell spacing in normal direction has a ∆y/Chord 1.5 × 10 −6 , which should assure that y + is below 2 for both the considered Reynolds numbers.At the finite trailing edge 16 cells are distributed with a cell size of ∆s/Chord 2×10 −4 at the sharp corners.At the leading edge the cell size is set to ∆s/Chord 1 × 10 −3 .The stretching in the chord-wise and normal directions is accomplished by double sided tanh stretching functions, giving low expansion rates.Most partners are using the common grid, but two partners used their own grids, namely the UOG and CENER.Additionally, some partners tested other mesh topologies (C, OCH), resolutions, and domain sizes.
Even though CENER tested the common O-mesh, the results presented here is computed on a C-H meshes topology generated in ANSYS ICEMCFD, and coarsened by the same method as the common grid, see Table 1.The UOG computations are performed on a special C-H mesh, where the coarser grids are regenerated reducing the number of points by a factor of 2.5 (L2) and 3.5 (L3) keeping the y + constant.

Grid refinement method
As the applied common grid is a block structured grid, the grid refinement/coarsening approach is kept simple.The coarser grid is constructed from the finest grid (L1) by removing every second point in both directions, resulting in the medium grid (L2).By repeating the coarsening procedure even coarser levels can be obtained.Here the coarsening is limited to two consecutive times, resulting in three grid levels, the finest L1, the medium L2, and the coarsest L3, see Figure 1 The advantage of this approach, identical to the approach which has been used within the EllipSys solver for years to do grid-sequencing, is that it is a very tough test where all parameters are varied at the same time.The cell size, the y + value at the wall surface, and the grid expansion rate is more than doubled, while the number of cells within the boundary layer is halved.Looking at the glide ratio, focusing on the AOA region between [-5:9] degrees the discretisation error for the best performing code is around 2 percent.The discretisation error of most of the codes are 4-7 percent.For three of the codes, the discretisation error is more than 10 percent.It is important to remember that grid independence does not guarantee that the solution is approximating the physics, only that the solution is a true representation of the implemented models.

Far-field boundary condition
Initial investigations revealed that even with a domain size of 40 Chords, a point vortex correction is needed by some solvers to account for the induced velocity in the far-field and to provide a domain size independent solution on the common O-mesh, eg. the FLOWer solver.Other solvers worked well for the common O-mesh configuration, while CH-and OCH-meshes, constructed with a similar domain size, revealed the need for a vortex correction on these topologies to avoid influence on the integral quantities from the pressure disturbance in the far-field, eg. the EllipSys2D solver.The influence from a lifting airfoil on its surroundings can be approximated from simple potential flow, using the Kutta-Jukowski theorem and the Biot-Savart law, resulting in the following boundary condition at the outer boundary: where C l is the Lift coefficient per unit span, U ∞ is the free stream velocity, and r/Chord is the normalized distance from the domain centre.As an alternative to applying the point vortex correction, the expression can be used to determine the necessary domain size to reduce the neglected induced velocity at the far-field below a given threshold value.The effect of including the vortex correction in the EllipSys2D computations is illustrated in Figure 2, where it is evident that the correction limits the distortion of the iso-bars near the outer boundary.

Fully turbulent results
A direct comparison of all partner results on the finest grid level (L1) with respect to the glide ratio (Cl/Cd), is shown in Figure 3.Here all partners are shown with identical colour, with the sole intention of indicating the spread.The glide ratio is an important parameter in connection with wind turbine design, and additionally a highly sensitive parameter.At the maximum glide ratio, the difference in the glide ratio at the low and high Reynolds numbers are 8 and 18 percent, respectively.
In Figure 4, a representative example of the spread between the seven individual results is illustrated.As seen in the left side of the figure, the pressure distributions are in very good agreement, and only one of the codes can be seen to deviate slightly.For the skin-friction, one of the codes show some numerical wiggles, while the remaining codes show acceptable mutual agreement.
A representative picture can be seen when comparing the grid dependence of the glide ratio computed by the different solvers for the fully turbulent simulation at Re=15 × 10 6 , using consecutive coarsenings of the grids, see Figure 6.Several of the solvers do not provide grid independence results, even in the attached flow region of [-8:10] degrees, where RANS solvers are expected to provide good predictions.It is also evident from this figure, that grid dependence varies with the AOA.Based on this, it is concluded that the grid dependence must be investigated for the full AOA range of interest.As grid independence is a necessary prerequisite for comparing the solution between different solvers, the three solvers having the highest variability with grid refinement from level 2 to level 1 are excluded, and the four remaining solvers are compared.Comparing the four remaining solvers, see Figure 5, the deviation between the solvers is down to less than 2 percent in glide ratio.Looking at the lift and drag values at 8 degrees, close to maximum glide ratio, it can be derived from Table 2 that most codes are agreeing within 0.25 percent in lift, and 2 percent in drag.

Results including the laminar-turbulent transition process
A similar exercise is done for the transitional cases, looking again at Reynolds numbers of 3×10 6 and 15 × 10 6 .In the present study, only three codes were applied successfully to the transitional computations, all based on the e N type transition models.Results are not included for the correlation based transition model by Menter and Langtry [15], [16], available in some of the codes, as it fails to correctly predict the natural transition behaviour at high Reynolds numbers.This was first experienced by some of the authors in connection with the studies of the DTU 10 MW turbine, in the INNWIND EU project back in 2013, and illustrated in some detail in 2014 [31].
Looking at the glide ratio, as shown in Figure 7, several tendencies are similar to the observations for the fully turbulent computations.Two out of the three codes have only minor variation in glide ratio between the finest grid level L1, and the coarser L2, for angles between [-10:10].The third code exhibits large differences between the two finest grid levels, similar to its fully turbulent behaviour.Even though the third code exhibits a large variation from L2 to L1, the actual computed values on L1 is close to the values computed by the two other codes.To verify whether the result of the third code on grid L1 is grid independent, a refined version of the L1 grid was generated and tested.The results are not shown here, but the solution on the refined grid showed good agreement with the L1 solution.Generally, the transitional results are showing larger spread than the fully turbulent results, which indicate that transition prediction is not as well tested as the fully turbulent approach.

Discussion
Several lessons were re-learned by the partners during the exercise, when going into the new regime of high Re and low Mach that will be experienced on wind turbines in the future.
Grid generation has to be done very carefully.This implies using e.g.hyperbolic tangent stretching to limit the expansion rate in the high gradient regions, both chord-wise and in the normal direction.At the high Reynolds numbers dealt with here, the very small cells needed at the wall will result in slow iterative convergence for many solves.Comparing the results from the seven codes is illustrative, and it is observed that for codes formally of second order, the grid requirement to obtain grid independent results might differ.
The effect of the domain size is illustrated, and it is concluded that a large domain size of the order of 100 chords or including a vortex correction, is necessary to have consistent results between the codes.Additionally, it was observed that the individual formulations require a very different number of iteration or time-steps to converge.The fastest solver might require less than 1000 iterations, while others required more than 20000 iterations.The efficiency with respect to computing time, not discussed here, will depend on the cost of the iteration or time-step.
For the laminar-turbulent transitional computations, only three partners contributed with results.Even though the standard correlation based model is available in some of the codes, it is not applied due to the problem of predicting the natural transition on an airfoil at high Reynolds numbers.Generally, the three codes predict results in decent agreement, even though it is not as good as for the fully turbulent computations.Finally, it was illustrated that grid dependency should be investigated at all angles of attack of interest.For structured grids, a simple coarsening where every second cell is removed in the chord-wise, and the normal direction, is well suited for the purpose.

Conclusion
This paper presents the effort to make seven different RANS CFD solvers provide consistent results.Even though the results are not fully satisfactory for all solvers, four out of seven codes are capable of predicting the maximum glide-ratio with a difference below 2 percent for fully turbulent conditions.To obtain this level of agreement domain size, careful treatment of boundary conditions, grid quality, and iterative convergence must be in focus.
More work is need to fully establish the application of laminar-turbulent transition models for the flow regime of the wind turbines of the future.This conclusion is based on the fact that several codes are not capable of providing transitional results, and that the mutual agreement between the results from the codes are inferior to the fully turbulent behaviour.

Figure 1 .
Figure 1.Details of the grids near the surface for the described grid coarsening from left to right, L1, L2, and L3.

Figure 2 .
Figure 2. The effect of a point vortex correction on the pressure, illustrated by iso-bars showing the uncorrected results on an O-and OCH-mesh left and the corrected results right, based on the EllipSys2D results.The red colours represents high pressure values, while the blue colours represent low pressure values.

Figure 3 .Figure 4 .
Figure 3.The spread in the glide ratio predicted by the individual partners, fully turbulent conditions at Re=3 × 10 6 left, and Re=15 × 10 6 right.

Figure 6 .
Figure 6.Grid dependency of the fully turbulent results for the individual partners at Re=15 × 10 6 .

Table 1 .
Mesh topologies used in the simulations.

Table 2 .
Lift and drag at an angle of attack equal to 8 degrees for fully turbulent conditions at Re=15 × 10 6 .