2D NACA 0012 airfoil validation
This is part one of a two article series on lift in 2D which uses the NACA 0012 airfoil to illustrate some concepts related to lift. The purpose of this validation is to compare our CFD results against known data to certify that we reproduce the physics correctly. Since we will rely on these data to draw conclusions in part two of this series, being able to trust the results of our model is vital.
The reference validation case is available through NASA Langley Turbulence Modeling Resource (TMR) website, while detailed results have been reported through various papers such as Diskin et al. (2015).
Contents
Setup
Here we will give only a brief outline of the case setup; for more details you can download the finest mesh case at the end of this post, here.
Grid
NASA makes available several grids used in their verification and validation studies. Among the available options, I have opted for the Family II of structured grids, from size 225 x 64 to size 1793 x 513. The grids have been generated such that the airfoil closes at chord = 1 with a sharp trailing edge. In all available grids, the farfield extends about 500c.
The minimum length —after scaling— at the wall is 7 x 10-7 and average stretching rate normal to the wall is ~1.02 for the points near the wall. As a result, the grid contains approximately 1 million cell.
The structured PLOT3D grids can be easily converted to OpenFOAM format with the plot3dToFoam tool.
Boundary conditions
I’ll be using the results from Diskin et al. (2015) and the TME website as reference for this case. There, we get that the conditions for the case are M = 0.15, Reynolds number per chord is Re = 6 million, alpha = 10 deg, reference temperature = 540 R. The Prandtl number Pr is taken to be constant at 0.72, while the turbulent Prandtl number Prt is taken to be constant at 0.9. The heat capacity ratio ($\gamma$) is 1.4.
The conditions above states are given for nondimensional CFD codes such as CFL3D or FUN3D. That’s not the case for OpenFOAM, so the conditions need to be translated. Firstly, we convert the reference temperature from Imperial to SI (540 R = 300 K). Secondly, once we know temperature, we can compute the speed of sound as
$$ a = \sqrt{\gamma R T} $$
where $R$ is the specific gas constant (287.058 J⋅kg−1⋅K−1 for dry air). As a result we get a speed of sound of $a$ = 347.2 m/s. Thirdly, now we can use the Mach number to get the flow velocity, $U_\infty$ = 52.08 m/s. Fourthly, we can get additional dry air properties such as kinematic viscosity $\nu_\infty$ = 15.68 x 10-6 m2/s. Finally, we can obtain the airfoil chord from Reynolds number, air velocity, and kinematic viscosity as
$$ L = \frac{\nu_\infty Re}{U_\infty} $$
which results in a chord of L = 1.80645 m. Therefore, the grid has to be scaled.
On the other hand, in the farfield we have also the following conditions:
$$ \tilde{\nu}_\rm{farfield} = 3\nu_\infty $$
$$ \nu_{t,\rm{farfield}} = 0.210438\nu_\infty $$
OpenFOAM
The following table summarises the boundary conditions used in OpenFOAM (front and back are of the empty type).
| airfoil | farfield | |
| U | fixedValue uniform (0 0) | freestreamVelocity* freestreamValue 52.08*(cos(10°) sin(10°)) |
| p | zeroGradient | freestreamPressure* freestreamValue 101360.2 |
| nuTilda | fixedValue uniform 0 | freestream* freestreamValue 47.04e-6 |
| nut | nutLowReWallFunction uniform 0 | freestream* freestreamValue 3.299668e-6 |
| T | zeroGradient | freestream* freestreamValue 300 |
| alphat | alphatWallFunction Prt 0.9 | calculated |
*The freestream, freestreamPressure, and freestreamVelocity boundary conditions inherit from inletOutlet and outletInlet. Those boundary conditions switch the mode of operation between fixed (free stream) value and zero gradient based on the sign of the flux.
Sutherland’s law
NASA repeats through several sections in their website —like here— that in their codes, dynamic viscosity is computed using Sutherland’s law,
$$ \mu = \mu_0 \left( \frac{T}{T_0} \right)^{3/2} \left( \frac{T_0 + S}{T + S} \right) $$
were $S$ it the Sutherland temperature (110.4 K). In OpenFOAM, Sutherland’s law is written differently, as
$$ \mu = \frac{A_s \sqrt{T}}{1 + S/T} $$
Comparing the formulas above, we can write the constant as:
$$ A_s = \frac{\mu_0}{T_0^{3/2}} (T_0 + S) $$
were $\mu_0$ = 1.176 x 10-5 kg⋅m-1⋅s-1, $T_0$ = 273.15 K, and $S$ = 110.4 K. Thus, $A_s$ = 1.458 x 10-6 kg⋅m-1⋅s-1⋅K-1/2.
We can enable Sutherland’s law in the thermophysicalProperties dictionary.
Note: I compiled a new implementation of the Sutherland transport model which uses the Prandtl number to calculate thermal conductivity. The code is accessible in GitHub: https://github.com/pfsq/mySutherland
The Spalart-Allmaras turbulence model
Ludwig Prandtl pointed out that certain fluid flows could be divided into a thin viscous layer —boundary layer— near surfaces and an effectively inviscid outer layer. In the inviscid outer layer the Euler and Bernoulli equations apply. But in boundary layers, which are predominantly turbulent, the assumptions taken by Euler and Bernoulli do not hold. Some empirical equations were developed for predicting the turbulent over a flat plate, yet the turbulent flow over arbitrarily shaped bodies still involves the solution of the continuity, momentum, and energy equations along with some model of the turbulence. Philippe Spalart and Steve Allmaras. proposed one of such models.
The one uncontroversial fact about turbulence is that it is the most complicated kind of fluid motion.
Peter Bradshaw, Imperial College
The Spalart-Allmaras model is a linear eddy viscosity that solves one additional transport equation. Because it is computationally cheaper, it is used in many codes and, for many flows, its performance is comparable to that of many more complicated models.
Many variants of the model have been proposed since its introduction in the early 1990s. Specifically, OpenFOAM’s implementation of the model ignores the $f_{t2}$ laminar suppression term; otherwise, the model is identical to the standard version. But as long as the simulation is at a relatively high Reynolds and $\tilde{\nu}_\rm{farfield} = 3\nu_\infty$, ignoring the $f_{t2}$ term should make little difference, as is the case here. In any case —for the inquiring eye— I have implemented the standard version of the Spalart-Allmaras model by reintroducing the $f_{t2}$ term: https://github.com/pfsq/standardSA; unsurprisingly, the results are almost identical.
In Diskin et al. (2015), the codes FUN3D and TAU use the SA-neg variant of the SA turbulence model that admits negative values of the Spalart turbulence variable, while for positive values of $\tilde{\nu}$ the model remains identical to the standard version. In any case, this variant should provide negligible differences in most cases.
Note: older versions of OpenFOAM (v2.x and older) need changes to the turbulence models to match the definitions found in NASA’s website.
Numerical schemes
The solutions from all three codes in Diskin et al. (2015) are second order accurate via a MUSCL scheme. Although these schemes are also available in OpenFOAM, I was getting unstable simulations, so I used linearUpwind instead, for velocity, energy, and turbulence, and linear for the rest. As with MUSCL, both schemes are second order accurate discretizations of the advective terms.
Since the solution was initiated with the coarsest mesh and then each successive solution mapped to the next grid level, it was not necessary to start with lower order schemes.
Linear solvers
I used the PCG solver with FDIC preconditioner for pressure —GAMG was very unstable— and PBiCGStab with DILU preconditioner for the asymmetric matrices.
For asymmetric matrices, we set a tight tolerance of less than 1e-10, while for pressure —as it is more expensive to calculate— we run with a relative tolerance of 0.001.
Solution
The cases were solved with a compressible solver, rhoSimpleFoam, and with the standard Spalart-Allmaras turbulence model. Each was iterated for 5000 steps —except for the finest mesh, which required 10000—, reaching convergence to at least 4 significant digits.
Convergence of residuals and variables of interest
Residuals go down nicely, although turbulence and energy seem to reach a lower limit at 1e-4 and 1e-7, respectively.