Precise drag prediction of airfoil flows by a new algebraic model

We report the results of accurate prediction of lift ($$C_L$$CL) and drag ($$C_D$$CD) coefficients of two typical airfoil flows (NACA0012 and RAE2822) by a new algebraic turbulence model, in which the eddy viscosity is specified by a stress length (SL) function predicted by structural ensemble dynamics (SED) theory. Unprecedented accuracy of the prediction of $$C_D$$CD with error of a few counts (one count is $$10^{-4}$$10-4) and of $$C_L$$CL with error under 1%-2% are uniformly obtained for varying angles of attack (AoA), indicating an order of magnitude improvement of drag prediction accuracy compared to currently used models (typically around 20 to 30 counts). More interestingly, the SED-SL model is distinguished with fewer parameters of clear physical meaning, which quantify underlying turbulent boundary layer (TBL) with a universal multi-layer structure, and is thus promising to be more easily generalizable to complex TBL. The use of the new model for the calibration of flow condition in experiment and the extraction of flow physics from numerical simulation data of aeronautic flows are discussed.


Introduction
Drag prediction continues to be the main challenge for computational fluid dynamics (CFD) simulation of aerodynamic applications, which has inspired such major international effort as drag prediction workshop (DPW) over more than a decade. Among all factors affecting the drag prediction accuracy, turbulence model in the Reynolds-averaged Navier-Stokes (RANS) framework plays a central role [1], due to the lack of a complete turbulent boundary layer (TBL) theory based on physical understanding of wall turbulence. Current models typically involve the modeling of correlation terms by dimensional argument, which is too simplistic and then inevitably introduces numerous coefficients/functions having no physical interpretation. Consequently, the determination of the parameters rely heavily on the experience of engineers, so that high-quality CFD becomes rather "artistic" and the generalization of a working program to new flows (with new features) becomes increasingly a formidable task (since model formulation becomes increasingly more complex). Theoretically speaking, the bottleneck lies on the lack B Zhen-Su She she@pku.edu.cn 1 State Key Laboratory for Turbulence and Complex Systems, College of Engineering, Peking University, Beijing 100871, China of a understanding of wall turbulence, in terms of its universal structure. The mixing length theory of Prandtl [2] and the log law of von Karman [3], and later the similarity hypothesis of Townsend, represent the milestones in this pursuit, but unfortunately they have not yet been developed to form a complete TBL theory, able to guide the development of turbulence model. The present work explores such a path, namely, developing accurate algebraic model based on well-founded physics theory of wall turbulence.
Recently proposed structural ensemble dynamics (SED) theory [4][5][6] is a symmetry-based approach, which expresses the (universal) wall constraint on turbulence eddies in terms of a (Lie-group) dilation invariance of a stress length (SL), 12 , quantifying eddies' scales. The theory predicts a multilayer analytic form (see below) for 12 [4][5][6] which includes a linear region corresponding to the log law in the mean velocity, but extends to the entire flow domain with four layers including viscous sublayer, buffer layer and a bulk flow region. The four-layer structure characterizing the entire TBL is physically sound and then capture the true similarity form of wall turbulence, for increasing Reynolds numbers (Re). The SED theory is able to predict its universal analytic form, based on a dilation invariance assumption, which is presumably universal whenever a wall is present. More interestingly, the variation of the (multi-layer) parameters may represent the variation among different TBLs of a variety of physical conditions. Indeed, the theory has yielded unified description of mean velocity as well as mean kinetic energy profiles in a variety of wall flows [4,5], from channel/pipe (which is internal flow) to TBL (which is external flow), from incompressible to supersonic and hypersonic TBL, and from smooth to rough wall. On the other hand, the theory, at the very beginning, aimed at developing a closure for the RANS equations: an order-function-based model was reported in 2009 [6], and a SED-k-ω model in 2017 [7]. The latter yields a uniformly accurate description of mean velocity and streamwise kinetic energy profiles in turbulent channel and pipe, covering a wide range of Re.
Recently, this theory is further extended to derive an algebraic model [8] for compressible turbulent boundary layers (CTBL), with accurate prediction of mean profiles of velocity, temperature and density, valid from incompressible to hypersonic flow regimes (thus achieving a Mach number (Ma) invariant description). More recently, a transition model [9] for transitional flat plate flow is shown to yield accurate prediction, for the first time, of the entire streamwise profile of the friction coefficient, C f , of all seven sets of experimental (e.g., of T3 series) and simulations data with varying incoming turbulence intensities (T u) (or scales), superior to existing closure models. This work explores its applicability to industrial turbulent boundary layers such as airfoil flows, to assess its potential merit for aeronautic applications.
Here, we focus on two transonic airfoil flows, NACA0012 and RAE2822, which stands as common benchmark flow for testing turbulence model for aeronautic applications. In these two airfoil flows, shock may form and interact with TBL, with notable pressure gradient. Two quantities of interest are the coefficients of lift and drag, C L and C D ; the latter is particularly challenging for currently used turbulence models, which have limited accuracy, due to inadequate description of turbulent boundary layer in the presence of varying pressure gradient of airfoil flow, compared to flat plate TBL for which the model parameters are calibrated. Specifically, the drag coefficient predicted by popular models for the two standard airfoil configurations always have (absolute) error of 20 to 30 counts (each count is 10 −4 ), compared to experimental measurements, despite extensive effort by the international CFD community. The remarkable message of the present work is that, with a refined physical (multi-layer) description of wall turbulence, this bottleneck has been overcome.

Governing equations
Below, we present our new algebraic model of eddy viscosity specified by SL of the SED theory, thus hereafter referred to as SED-SL model. The SED-SL model starts from compressible RANS equations [10], which express the mean conservation of mass, momentum and energy. Under the Boussinesq approximation and assuming constant turbulent Prandtl number, the RANS equations is written as follows where ρ is fluid density, U j is the jth component of velocity and x j represents Cartesian coordinates. p = ρ RT is the static pressure (R is the gas constant) and entropy H = E + p ρ , where total energy E = C ν T +k with C v the specific heat at constant volume, T the static temperature and k the kinetic energy of the flow. The laminar Prandtl number Pr is 0.7, and the turbulent Prandtl number Pr τ is set to be 0.9. Meanwhile, the the Boussinesq approximation is expressed by where δ i j is the Kronecker delta function defined as δ i j = 1 for i = j and δ i j = 0 for i = j. S i j denotes the strain rate tensor, expressed as follows

SED-SL model
While the molecular viscosity μ is computed, using Sutherlands law: μ where T S = 110.4 K and subscript ∞ indicates the free-stream property, the socalled turbulent viscosity μ T is described according to the SED theory, by the stress length function l 12 : The SED theory predicts that the stress length l 12 has a multi-layer form in fully developed turbulent boundary layer [4]. It reads where y + is the distance from the wall in wall unit (i.e. which is well-defined for attached flows and mild separated flows, where du/dy| w > 0 ), r = 1 − y/δ = 1 − y + /Re τ (Re τ = δu τ /ν, friction Reynolds number) is the distance from the edge of TBL whose location is at δ. Equation (7) describes a canonical four-layers structure of wall turbulence, including viscous sublayer (with thickness 9.7), buffer layer (with thickness y + bu f ), log layer and bulk flow (with 1−r 4 structure); it is easy to verify that for y + y + bu f , we obtain the celebrated linear law +inner 12 ≈ κ y + for the matching function between the inner and outer region, with 0 = 9.7 2 κ/y + bu f , where κ is the Karman constant. The outer region is described using the variable r , and δ is defined as the position of 0.99U e for flat plate TBL (U e is free stream velocity), but generalized here for airfoil flows, following Baldwin-Lomax model [11] by where y max = location max +inner 12 |ω| .
Note that there are two important physical parameters in the multi-layer Eq. (7), 0 and y bu f , which have clear physical interpretation: the former is the size of very near-wall eddy (e.g. at y + ≈ 9.7) and the latter the buffer layer thickness, where near-wall coherent vortices are formed. They together characterize TBL and thus must evolve in the streamwise direction. In order to describe the streamwise variation and, in particular, the laminar-turbulent transition, we extend our multi-layer formalism [4] by postulating a simple two-state ("laminar" versus "turbulent") structure [9]: 0 and y bu f take a similar multi-layer form, provided that a dilation invariance is valid in the streamwise direction x, with respect to the leading edge location (i.e. x = 0): where x * denotes the middle transition location from laminar to fully developed turbulent regime, which is a single parameter, set to be the location where the transition is forced in experiment for airfoil flows studied below. The exponent γ = 0.5 is the Blasius scaling describing the growth rate of the very near-wall eddy size following the growth of the laminar boundary layer thickness (δ ∼ x 0.5 for x x * ) for transition due to incoming turbulence, but is different for triggering by blow and suction disturbance at a specific location. In the latter case, γ = 20 is used to describe a sudden transition. On the other hand, the exponent 5.5 in Eq. (10) is obtained from fitting direct numerical simulation (DNS) data of flat plate TBL, which yields an empirical κ ∼ x 6 (for x x * ), which seems to be universal for all cases studied below. Note that the SED-SL transition model is significantly simpler than the popular approach adding the transition to RANS model by introducing an intermittency transport equation [12,15,16]; for details, see Ref. [9].
Finally, the SED-SL model has two parameters, ∞ 0 and y ∞ bu f , characterizing the fully developed TBL, which may be slightly different for different physical conditions such as wall curvature, pressure gradient, or Ma. For flat plate TBL (presented in Sect. 3.1), they are set to be 1.1 and 41, respectively, (yielding a Karman constant κ ≈ 0.45), according to the original SED theory [4,5]. For the two airfoil flows (presented in Sects. 3.2 and 3.3), we empirically determine that ∞ 0 ≈ 0.3, while y ∞ bu f ≈ 85 (for yielding a smaller κ of 0.27), which are fixed for all cases and all angle of attack (AoA). In the future, more simulation studies need to be conducted to establish a systematic variation of ∞ 0 and y ∞ bu f with Ma and wall curvature, in order to develop a wide range of engineering predictions.

Numerical method and model implementation
The above model is implemented on a so-called OpenCFD-EC platform, which is a parallel, multi-blocks, structuredgrid, finite-volume compressible CFD software developed for engineering computation by Li et al. [17] at the Institute of Mechanics of Academy of Science of China. In this code, a third-order MUSCL-scheme and a van-Leer method are adopted to discretize the inviscid terms and the viscous terms, respectively; an implicit LU-SGS method is used for time advance. Our model is also successfully implemented in CFL3D platform [18], which is a structured-grid, cellcentered, upwind-biased RANS code, developed in the early 1980's in the Computational Fluids Laboratory at NASA Langley Research Center used to support numerous NASA programs since then and continues to be used today. The results (e.g. accuracies) are similar using both codes, with slightly adjustment of the model parameter ∞ 0 (e.g. changing from 0.3 to 0.5 owing to different numerical viscosity associated with different numerical schema).
In order to assess the quality of the predictions, we compare our model with the popular γ − Re θt transition model [12] for the flat plate transitional flow case, and with two wellknown turbulence models, namely the Spalart-Allmaras (SA) one-equation model [13] and the shear stress transport (SST) two-equation model [14], for the airfoil flows; the results clearly show that the SED-SL model perform the best. But we do not claim the general superiority of the SED-SL model over the popular models, since more flows need to be tested.

Results and discussion
Below, we first present the results of the prediction by the SED-SL model for transitional flat plate flow for a basic validation [9], following the validation of this model for the description of fully developed flat plate flows from incompressible to hypersonic regimes [8]. Then, we report the results of supercritical airfoil RAE2822 and symmetric airfoil NACA0012. For the former, we calculate the flow for a wide range of Re, and for the latter, the computation is carried out for different AoA, with significant improvement of the prediction accuracy of the lift (C L ) and drag (C D ) coefficients compared to the two popular models (e.g. SA and SST). Note that for latter, experimentally reported incoming flow conditions (Ma and Ao A) are slightly modified in a systematic way, which is explained in Appendix.

Transitional flow over flat plate
The first test case is chosen to be an incompressible flow over ERCOFTAC T3B flat plate [19] with free stream turbulence intensity of 6.5%. As the transition is induced by freestream turbulence, the present model parameter γ is set to 0.5, and the transition Re by a simple correlation formula, following [9]. The flow is computed on an 500 × 150 grid over the calculation domain 1.0 m × 0.1 m, with a free stream boundary condition applied at the upper, inlet and outlet boundaries, and no-slip condition at the (bottom) flat plate wall. The grids are homogeneous in the streamwise direction, while stretched exponentially in the wall-normal direction with a small value for the closest mesh point to the wall: y + 1 ≈ 0.3. Figure 1(a) reports the variation of the skin friction (C f ) along Re x = U in f x/ν (U in f is the free stream velocity, x denotes the streamwise coordinate and ν represents kinetic viscosity) for T3B flat-plate test case [19] with zero pressure gradients, showing much better agreement of the present model, compared with the γ − Re θt transition model [12]. As shown in Fig. 1(b), the errors of the predictions by the SED-SL model are bounded within 5%, which is clearly below that by Menter's model, especially in the transition region. wall normal direction y + , at different streamwise locations for T3B transitional flat plate flow with our DNS data. Note that the incoming turbulence intensity is 6.5%, which is rather large, so that a bypass transition takes place, which is difficult to describe for most models. Nevertheless, the predicted mean velocity profiles by the SED-SL model are in almost perfect agreement with DNS results; the Reynolds stress profiles have slight deviations, only near the transition region.

Transonic flow over RAE2822 airfoil
Then, consider flow passing RAE2822 airfoil for case 6, whose experimental condition is well documented in Ref. [20]. Note that previous studies showed that the experimentally reported flow conditions, mainly Ma and Ao A, need to be rectified, in order for the computation to reproduce accurately the measured lift and drag coefficients. For case 6, the used flow conditions was suggested in NASA turbulence a b Fig. 2 Predicted mean velocity and Reynolds stress profiles compared with DNS data of T3B transitional flow at three streamwise locations with different Re x . Black symbols: DNS (squares for Re x = 0.94×10 5 ; circles for Re x = 1.27 × 10 5 ; triangles for Re x = 4.40 × 10 5 ); Red lines: SED-SL model validation and verification website: Re = 6.5 × 10 6 , Ma = 0.729, Ao A = 2.31 • , while for case 9, we use following flow conditions: Re = 6.5 × 10 6 , Ma = 0.732, Ao A = 2.65 • , as discussed in Appendix. Using these values, we compute the flow passed RAE2822 airfoil, with a standard mesh of 395 × 65 with 306 points lying on the airfoil surface, provided by NASA turbulence modeling validation page. Grid convergence is fully achieved as shown in Appendix (in Fig.  6(a)). Figure 3 plots the pressure coefficient distribution C p around the whole surface (x/x c , x is the streamwise locations of surface points, x c denotes the chord length of airfoil), and the friction coefficient C f distribution along the upper surface. It is seen that the shock position is well captured, significantly closer to experimentally measured values than the prediction by (standard) SA and SST model; the latter show a too early shock formation with a under-estimation of pre-shock pressure level. Same results can be seen for case 9, as shown in Fig. 4, where the results predicted by CKDO model [21] with Ma = 0.73, Ao A = 3.19 • (same as experimental flow condition without correction) are also included for comparison (black solid lines). The prediction of preshock skin friction on upper surface by the CKDO model are yet much higher than experimental measurements, while pressure coefficient distribution is captured better than SA and SST models.
Then, remarkably, the predicted lift and drag coefficient by the SED-SL almost exactly matches the experimental measurements, sharply distinct from the two other popular models, as shown in Table 1. The deviation of the SED-SL model for case 6 and case 9 are bounded within 1% for C L , and within 2 counts for C D (below 1 count in case 6). The drag accuracy is clearly superior to the SA, SST and CKDO a b Fig. 4 Pressure coefficient and friction coefficient along the upper surface of case 9 for RAE2822 airfoil flows: comparison between experimental data (black symbols) from AGARD report and computation results of SED-SL (red solid lines), SA (blue dash lines), SST (brown dot lines) and CKDO (black solid lines) models (whose maximum deviation up to nearly 30, 10 and 36 counts, respectively); the errors are reduced by nearly one order of magnitude. Note that our reported results by the SA and SST model are obtained using the same OpenCFD software with the same grid resolution, so the comparison is meaningful.
One would naturally wonder if this astonishingly small error of a few counts is meaningful, under experimental uncertainty. On one hand, experimental flow condition may need to be corrected, due to wall interference. The present accuracy is obtained under a systematic flow condition rectification; see Appendix B. We will show below that such high accuracy can be achieved for other cases of flow passing RAE2822, in a future communication. On the other hand, the evidence provided by the results of another airfoil flow passing NACA0012 for varying Ao A further enhance the conclusion that the SED-SL model delivers significantly better predictions than the existing models.

Flow over NACA0012 airfoil
Flows passing NACA0012 airfoil model were performed in the ONERA S3MA transonic wind tunnel by Thibert et al. [22]. Test cases selected for a range of Ao A under Ma ≈ 0.7, are computed. The computation is performed on 449 × 129 grids with 257 points lying on the airfoil surface, provided also by NASA Langley Research Center Turbulence Modeling Resource [23]. For this set of grids, the far field is located approximately 500 chords from the airfoil surface. Grid convergence is shown in Fig. 6 Under a systematic adjustment of Ao A explained in the Appendix, the predictions of C L and C D are compared to experimental data, as shown in Fig. 5. For all cases, the relative errors of C L for the SED-SL model is bounded within 2%, while the absolute errors of C D is within 6 counts, which is close to the experimental uncertainty [22]. On the other hand, the SA and SST model have errors for C L up to 20%, with, in addition, a systematic variation with increasing Ao A, showing that both models do not capture the right similarity TBL structure for varying Ao A. The errors of SA and SST for C D also widely spread in the range of 10 to hundreds counts, with again a wild variation of Ao A, with mostly significant over-predictions. In sharp contrast, the SED-SL model's errors remain unchanged for increasing Ao A, indicating that the multi-layer description by the SED theory, and hence the SED-SL model, capture the right similarity property of the TBL in airfoil flows, versus Ao A. Note that the SA and SST model used here are the standard version without any refinement, e.g. with no modeling of streamwise transition [17], and so they may not be the best version used in the industry. But it is suffice here for the purpose to demonstrate the effectiveness of the SED-SL model.
One may wonder how much the transition model contributes to the accurate drag prediction. We have investigated this question by removing the transition from the SED-SL model; only very small changes in C L and C D are observed. For instance, for RAE2822 airfoil, a full turbulence SED-SL model yields changes of 0.2% and 0.3 counts in C L and C D , respectively. Therefore, adding the transition model is not a major factor for the accurate drag prediction in the present set of computations. Instead, the simultaneously accurate lift and drag predictions are due to refined multi-layer description of TBL structure over the airfoil surface, much better than existing models.
In summary, the SED-SL model has achieved an improvement of drag prediction by nearly one order of magnitude compared to the popular models. Thus, we assert that the SED-SL model has achieved a breakthrough in the drag prediction.

Conclusion
In summary, the SED-SL yields accurate prediction of C L and C D for airfoil flows over RAE2822 and NACA0012 in a range of Re and Ao A. The improvement of drag prediction accuracy by an order of magnitude over current popular models (approaching the limit of experimental uncertainty) is rather striking, in two aspects. First, it is a sign that the essence of wall turbulence is captured, namely, the universal multi-layer description of the stress length defines the invariant (in x) wall-normal distribution of the eddy viscosity, which is behind the success of the computation. It is remarkable that the extension from flat plate to airfoil flow is so straightforward that only two parameters slightly change: y +∞ bu f from 41 to 85, and ∞ 0 from 1.1 to 0.3. Secondly, it offers a more optimistic view for the future of RANS model: accurate prediction is possible when fundamental understanding of wall turbulence is taken into consideration.
Furthermore, three specific merits of the present algebraic model are noteworthy. First, it may produce an alternative and better (more smooth and accurate) wall constraint for unsteady simulations such as detached-eddy or large-eddy simulations (LES) (e.g. to be used in the recently developed constrained LES [24,25]. Second, it may provide a tool to calibrate the experimental flow conditions. Our study shows that the suggested corrections of Ao A for RAE2822 and NACA0012 are reasonable (see Appendix). Thirdly, for new applications such as flows at higher Ma (e.g. hypersonic flow) or flows over new geometrical configurations, for which simulation data (from DNS or LES) is available, the present model may be used to determine the corresponding multi-layer parameters (i.e. ∞ 0 and y +∞ bu f ) which specify flow physics of TBL. Note that the Ma-dependence of ∞ 0 and y +∞ bu f is generally smooth, due to the fact that the multilayer structure is universal. This consideration highlights the interest of carrying out a combined DNS and LES study with the SED-SL model: it would yield a refined description of the corresponding TBL structures, and, at the same time, an accurate turbulence model. Thus, it is highly recommended to carry out such study for a larger class of wall flows. Note that some very different wall flows may display new multilayer structure. An example is flow with separation under strong adverse pressure gradient, for instance, the recently concerned dynamic stall problem [26]; our study shows that additional layers need to be added, so that stall at higher Ao A and unsteady separation can be correctly predicted, the results of which will be reported elsewhere. the past fifteen years. MJX likes to thank Xin-Liang Li for the help in implementing the model to OpenCFD. We also thank Xi Chen for helpful discussions, and to Shiyao Li for assistance.
Open Access This article is distributed under the terms of the Creative Commons Attribution 4.0 International License (http://creativecomm ons.org/licenses/by/4.0/), which permits unrestricted use, distribution, and reproduction in any medium, provided you give appropriate credit to the original author(s) and the source, provide a link to the Creative Commons license, and indicate if changes were made.   Fig. 6, the results of the medium grids (red solid lines) of both two airfoils are almost identical to that of the fine grids (blue dotted lines). Thus, the medium grid is chosen as the standard grid for simulation with balanced efficiency and accuracy.

A.2 Correction on flow conditions
Previous studies show that experimentally reported flow conditions for flows over airfoils, mainly Ma and Ao A, need to be rectified, in order for computation to match more accurately experimentally measured lift and drag coefficients.
Here, we discuss a procedure using the SED-SL model for this rectification, which appears to show a systematic feature, attributable to flow condition corrections due to wall interference and model deformation during experiment, which has been widely discussed by numerous researchers (see Refs. [28,29]). Our procedure consists in examining the polar curve (C D versus C L ), for the measured and computed results, with a systematic variation of Ma and Ao A. It turns out that for the SED-SL model, there always exists an optimal set of flow conditions (Ma and Ao A) such that the computed C D and C L are very close to the measured ones. Then, we take the corresponding flow condition to be the correct one. Note that such optimal flow condition (with the same C D and C L ) can not be found for any other model like SA or SST (because simultaneously accurate predictions of C D and C L are impossible for SA or SST).
The above procedure is validated by the case 6 flow over airfoil RAE2822, for which our suggested correction (i.e. Re = 6.5 × 10 6 , Ma = 0.729, Ao A = 2.31 • ) exactly coincides with the ones suggested by NASA turbulence validation and verification website [27]. Furthermore, we carry out the procedure for a series of flow conditions such as cases 7 to 9 of RAE2822 airfoil flows, as shown in Table 2. In both cases, the corrections seem to be sound, with a nearly constant Ma correction ( Ma ∼ 0.004) and Ao A correction ( Ao A ∼ −0.6 • ). Such corrections suggest a possible systematic mis-calculation of the flow conditions in the experiment; note that we hold Re unchanged, to be 6.5 × 10 6 , since it is generally more accurately determined.
Similarly, for flows over NACA0012 under Ma = 0.7, our suggested corrections on Ao A (with Ma keeping unchanged) are shown in Table 3. Ao As of experiment before and after previously considered wall interference corrections [30] are also shown, indicating that our suggested corrections further extend the previously suggested correction magnitude (by wall interference), excepted at the largest Ao A. In most cases, we believe that it may be attributed to model deformation. But for Ao A approaching to the largest angle, some dynamic feature occurs, indicating a deformation of the TBL multilayer description. Note that the present suggested correction needs to be considered together with other modeling, which will be discussed elsewhere.