# A high-order local correlation based transition model for transonic airfoil flows

Matthias Plath, Florent Renac, Olivier Marquet, Christian Tenaud

► To cite this version:

Matthias Plath, Florent Renac, Olivier Marquet, Christian Tenaud. A high-order local correlation based transition model for transonic airfoil flows. Computers and Fluids, 2024, 285, pp.106461. 10.1016/j.compfluid.2024.106461. hal-04857496

HAL Id: hal-04857496

https://hal.science/hal-04857496v1

Submitted on 16 Jan 2025

HAL is a multi-disciplinary open access
archive for the deposit and dissemination of sci-
entific research documents, whether they are pub-
lished or not. The documents may come from
teaching and research institutions in France or
abroad, or from public or private research centers.

L'archive ouverte pluridisciplinaire HAL, est
destinée au dépôt et à la diffusion de documents
scientifiques de niveau recherche, publiés ou non,
émanant des établissements d'enseignement et de
recherche français ou étrangers, des laboratoires
publics ou privés.

<figure><img src="image_1.png" /></figure>

Distributed under a Creative Commons Attribution 4.0 International License
<figure><img src="image_2.png" /></figure>

<figure><img src="image_3.png" /></figure>

# A high-order local correlation based transition model for transonic airfoil flows

Matthias Plath<sup>a,*</sup>, Florent Renac<sup>a</sup>, Olivier Marquet<sup>b,c</sup>, Christian Tenaud<sup>d</sup>

## ARTICLE INFO

MSC:

0000

1111

Keywords:

Transition model

RANS equations

Discontinuous Galerkin

Transonic

## ABSTRACT

Accurate transition modeling is a key ingredient towards the all-time goal of improved aerodynamic performance in the aircraft industry. While several different approaches exist to model transition, local correlation based transition models (LCTM) are one of the most versatile strategies. Many versions of LCTM were proposed in the context of lower order finite volume discretizations, but so far only a few transition models have been brought to higher-order numerical schemes. The present work represents a first attempt of including a one equation LCTM into the high-order discontinuous Galkerin (DG) framework. Several modifications are accessed to improve the model's robustness and tailor it for transonic airfoil flows. Validation is conducted for a series of test cases, starting with zero pressure gradient flat plates, progressing through high Reynolds number airfoil flows, and concluding with a complex unsteady transonic buffet scenario. All results are compared against experiments and previous numerical references.

## 1. Introduction

In recent years, the development of different approaches to predict laminar-to-turbulent transition is receiving more and more attention. One main driving force behind this development, are plans by the aircraft industry to include laminar flow wings into their next generation aircraft designs. Laminar flow techniques are encouraging candidates to significantly decrease the aerodynamic drag and hence improve aerodynamic efficiency and fuel consumption for civil aircraft. While the benefits are very promising, including accurate transition prediction capabilities into an industrial computational fluid dynamics (CFD) framework are not straight forward. Since the transition process is inherently nonlinear and can be triggered through several different mechanisms, there exists no unified physical approach to model it. The only exception are computationally highly demanding direct numerical simulations [1], which are up to now not close to be applicable for most industrial configurations. Considering that chord Reynolds numbers for commercial airliners can reach values up to several tens of million, even LES simulations are out of reach for aircraft in cruise configuration. During cruise flight conditions, the main transition mechanisms are natural transition due to Tollmien-Schlichting (T-S) waves and cross-flow transition. Both mechanisms can be sufficiently well predicted by linear stability based approaches such as the $e^N$ method, first

introduced by van Ingen [2,3]. The coupled usage of the $e^N$ method and an industrial RANS solver is the current standard in the aircraft industry, as described in Krumbein et al. [4]. While such an approach can give reliable results even for full aircraft configurations [5,6], its implementation for unstructured CFD solvers is inherently complex.

Since stability based methods are not well suited for unstructured CFD, different classes of transport equation based transition models were developed. The transport equation based models can be separated into the groups of correlation based [7,8] and physical/phenomenological models [9,10]. Phenomenological models were first introduced by Walters et al. [9,10] and usually apply the strategy of a laminar kinetic energy (LKE) by Mayle & Schulz [11]. In case of LKE models, a modeled transport equation for the laminar kinetic energy is solved, where trigger functions inside the source terms control the transition process. The laminar kinetic energy is then coupled to a turbulence model by injecting it into the equation for the turbulent kinetic energy. While being a very promising approach, LKE models momentarily mainly suffer from the complexity to integrate several different transition mechanisms at once. The group of correlation based transport equation models tries to alleviate this problem and was introduced by Menter et al. [7]. Commonly, those models are formulated within the framework of a modeled intermittency transport

<footnote><sup>a</sup> ONERA, Université Paris Saclay, DAAA, 29 Avenue de la Division Leclerc, Châtillon, F-92322, Île-de-France, France</footnote>

<footnote><sup>b</sup> ONERA, Université Paris Saclay, DAAA, 8 Rue des Vertugadins, Meudon, F-92190, Île-de-France, France</footnote>

<footnote><sup>c</sup> ONERA, Institut Polytechnique de Paris, DAAA, 8 Rue des Vertugadins, Meudon, F-92190, Île-de-France, France</footnote>

<footnote><sup>d</sup> Université Paris-Saclay, CNRS, CentraleSupélec, Laboratoire EM2C, 8-10 Rue Joliot Curie, Gif-sur-Yvette, F-91190, Île-de-France, France</footnote>

<footnote>* Corresponding author.</footnote>

<footnote>E-mail address: matthias.plath@onera.fr (M. Plath).</footnote>
equation. Inside the source terms, correlation functions are used to predict transition based on the critical momentum thickness Reynolds number $Re_{\theta c}$. In the beginning, these type of models were still requiring integrated quantities, making them non-local. However, since Langtry & Menter [8,12] proposed the first fully local version, many different local versions were developed [13–16]. Because a correlation function is used to predict transition, those models can be explicitly tailored for various transition mechanisms such as natural and bypass transition [13,17], crossflow transition [18,19] as well as compressibility effects [20].

For modern industrial CFD solvers, novel higher-order discretization schemes suitable to complex unstructured meshes are envisioned to be employed. One such promising technique is the discontinuous Galerkin (DG) method [21–25]. The DG method offers a flexible higher-order scheme, characterized by an arbitrary selectable high-order and hence low numerical dissipation. Furthermore, the DG method allows the use of unstructured meshes, which make it well suited for complex geometries. Additionally, due to its discontinuous formulation and if used with suitable numerical fluxes, it is ideal for compressible flows showing shocks [26]. So far, only limited efforts have been made to bring RANS transition models into higher-order discretizations. A first attempt was made by Crivellini et al. [27] for laminar separation bubbles, by employing the Spalart–Allmaras model subject to very low inflow turbulence boundary conditions. Lorini et al. [28] introduced the first non-local version of a transition model discretized by a DG method. Following Lorini et al. [29] proposed a fully local version of a LKE based transition model. Recently, Crivellini et al. [30] integrated modified versions of the local correlation based algebraic SABCM [31] and $\gamma$-k-$\omega$ models into a high-order DG framework. This approach is particularly promising as it eliminates the need to solve additional transport equations for transition variables. The aim of the present work is to introduce a local correlation based transition model (LCTM) into the high-order DG framework Aghora [32] for external aerodynamic applications. Therefore, Menter et al.’s. [14] one equation model concept is coupled to a regularized Spalart–Allmaras turbulence model. A logarithmic transformation [33] is applied to the intermittency equation to achieve positivity preservation of the model and increase robustness. A critical momentum Reynolds number correlation according to Fehrs [34] is employed, tailored for high Reynolds number, transonic flows. A main focus is put into making the model robust, capable of simulating transonic flows, while achieving deep and fast convergence to steady state. Therefore, the model is made fully differentiable and subsequently analytically linearized for an implicit quasi-Newton method. Model validation is performed for: several zero pressure gradient flat plate cases offering natural T–S (Tollmien–Schlichting) and bypass transition, subsonic laminar NLF0414F airfoil at high Reynolds number $Re = 9.5 \cdot 10^6$, transonic laminar CAST10-2 airfoil and transonic laminar OALT25 airfoil. Finally, the model is even validated for an unsteady transonic buffet case around the OALT25 airfoil, while comparing against LES results [35].

The paper is organized as follows. Section 2 introduces the transition model and its discretization with the DG method is described in Section 3. Numerical results about steady solutions are provided in Section 4 to validate the present approach, while steady and unsteady simulations of the transonic OALT25 airfoil are given in Section 5. Concluding remarks about this work are summarized in Section 6.

## 2. SA-neg-log($\gamma$) transition model

In this section, we introduce a transition model, which is based on the one equation $\gamma$ transition model by Menter et al. [14] coupled to the well-known Spalart–Allmaras (SA-neg) turbulence model [36,37]. The starting point for the transition model is a transport equation for the intermittency factor $\gamma$ [14] to which we apply a logarithmic transformation $\tilde{\gamma} = \ln \gamma$. This technique was originally introduced by Ilinca & Pelletier [33] in the context of turbulence modeling in

order to ensure positivity preservation and hence increase robustness. Furthermore, the concept of logarithmic variables was first applied in a high-order DG framework by Bassi et al. [38] for solving the RANS k-$\omega$ equations. The transport equation for $\tilde{\gamma}$ is therefore derived by taking the natural logarithm of the classical $\gamma$ equation and is given as

$$ \partial_t(\rho\tilde{\gamma}) + \nabla \cdot (\rho\tilde{\gamma}\mathbf{v}) = P_{\tilde{\gamma}} - \epsilon_{\tilde{\gamma}} + \nabla \cdot \left( \left(\mu + \frac{\mu_t}{\sigma_{\tilde{\gamma}}}\right) \nabla\tilde{\gamma} \right) + \left(\mu + \frac{\mu_t}{\sigma_{\tilde{\gamma}}}\right) \nabla\tilde{\gamma} \cdot \nabla\tilde{\gamma}, \quad (1) $$

where $\rho$ defines the density, $\mathbf{v}$ the velocity vector, $\mu_t$ the turbulent viscosity and $\sigma_{\tilde{\gamma}} = 1$ is a model constant. The last term on the right hand side of Eq. (1) is a result of the logarithmic transformation. The production $P_{\tilde{\gamma}}$ and destruction term $\epsilon_{\tilde{\gamma}}$ are given by

$$ P_{\tilde{\gamma}} = F_{\text{length}} \rho S (\gamma - 1) F_{\text{onset}} \quad (2a) $$

$$ \epsilon_{\tilde{\gamma}} = \rho c_{a2} \Omega F_{\text{turb}} (c_{e2}\gamma - 1), \quad (2b) $$

where $S$ defines the strain rate magnitude and $\Omega$ the vorticity magnitude. Following the simplification process applied by several authors [14,34,39] the transition length $F_{\text{length}} = 14$ is taken as a model constant. The model constants $c_{a2} = 0.06$ and $c_{e2} = 50$ are taken similar to Menter et al. [14]. The trigger functions controlling the on and offset of transition are defined according to Menter et al. [14] as

$$ F_{\text{onset1}} = \frac{Re_V}{2.2 Re_{\theta c}}, \quad (3a) $$

$$ F_{\text{onset2}} = \min(F_{\text{onset1}}, 2), \quad (3b) $$

$$ F_{\text{onset3}} = \max\left(1 - \left(\frac{R_T}{3.5}\right)^3, 0\right), \quad (3c) $$

$$ F_{\text{onset}} = \max(F_{\text{onset2}} - F_{\text{onset3}}, 0), \quad (3d) $$

$$ F_{\text{turb}} = e^{-\left(\frac{R_T}{2}\right)^4}, \quad (3e) $$

$$ R_T = \frac{\mu_t}{\mu}, \quad (3f) $$

where $Re_V$ depicts the vorticity Reynolds number and $Re_{\theta c}$ the critical momentum thickness Reynolds number. In equation Eq. (3f) $R_T$ differs from the original formulation by Menter et al. [14] in order to couple the transition model to the SA turbulence model as introduced by several other authors [39,40]. The vorticity Reynolds number $Re_V$ in Eq. (3a) is computed as

$$ Re_V = \frac{\rho_e d_w^2 S}{\mu_v}, \quad (4) $$

where $d_w$ describes the wall distance and the subscript $(\cdot)_e$ indicates boundary layer edge values. The use of boundary layer edge values instead of simple local values was introduced by Fehrs [34], who pointed out that relation Eq. (3a) follows from similarity solutions of compressible Falkner-Skan profiles where $Re_\theta$ is only defined in terms of boundary layer edge quantities and hence $Re_V$ needs to be computed accordingly. By using isentropic relations and assuming a constant pressure normal across the boundary layer $\rho_e = \rho$ (local pressure value) one can compute the density $\rho$ and temperature $T$ at the boundary layer edge as

$$ \rho_e = \rho_\infty \left( \frac{p_e}{p_\infty} \right)^{\frac{1}{\kappa_{PG}}} \quad (5a) $$

$$ T_e = T_\infty \left( \frac{p_e}{p_\infty} \right)^{\frac{\kappa_{PG}-1}{\kappa_{PG}}}, \quad (5b) $$

where $\kappa_{PG}$ defines the specific heat ratio. The use of Sutherland’s law allows then to compute $\mu_e$ as

$$ \mu_e = \mu_\infty \left( \frac{T_e}{T_\infty} \right)^{\frac{3}{\kappa_{PG}}} \frac{T_\infty + 110.4}{T_e + 110.4}. \quad (6) $$

The correlation for $Re_{\theta c}$ in Eq. (3a) is one of the key parameters for accurate transition prediction in local correlation based models. Airfoil
<figure><img src="image_4.png" /><figcaption>Fig. 1. Comparison of modified Fehrs correlation [34] to correlations by Menter at al. [14], Langtry & Menter [12] for critical momentum thickness Reynolds number $Re_{\theta_c}$ as a function of the turbulence intensity $Tu$.</figcaption></figure>

flows are most often characterized by high Reynolds numbers and low turbulence intensities <1%. Hence, here we select a correlation from Fehrs [34], which was specifically developed for transonic, low turbulence environment, high Reynolds number airfoil flows. Accordingly, $Re_{\theta_c}$ is given by,

$$Re_{\theta_c} = C_{Tu1} + C_{Tu2} k_t Tu^{-C_{Tu3}}, \qquad (7)$$

where the parameter

$$k_t = e^{(3.2-3.2 \min(Tu,1))k_t} \qquad (8)$$

includes pressure gradient effects through the acceleration parameter $K_t$. Here $Tu$ defines the turbulence intensity, which is a constant parameter in this model. The acceleration parameter $K_t$ is defined by

$$K_t = - \frac{\nu_e}{\rho_e U_e^3} \frac{dp}{ds} Re, \qquad (9)$$

where the pressure gradient along the streamline is given by

$$\frac{dp}{ds} = \frac{v}{|v|} \cdot \nabla p. \qquad (10)$$

The boundary layer edge velocity $U_e$ is found by using isentropic relations and the compressible Bernoulli equation [20,41] as

$$U_e = \sqrt{U_{\infty}^2 + \frac{2\kappa_{PG}}{\kappa_{PG} - 1} \frac{p_{\infty}}{\rho_{\infty}} \left[ 1 - \left( \frac{p_e}{p_{\infty}} \right)^{1-\frac{1}{\kappa_{PG}}} \right]} \qquad (11)$$

The correlation constants in Eq. (7) were calibrated for low turbulence environment, high Reynolds number airfoil flows as

$$C_{Tu1} = 120, C_{Tu2} = 380, C_{Tu3} = 2/5. \qquad (12)$$

For bypass transition cases with $Tu > 1\%$ those constants give too high values for the critical momentum thickness Reynolds number as can be seen from Fig. 1. This leads to delayed transition onset for bypass transition cases and therefore the constants are modified as follows,

$$C_{Tu1} = 120, C_{Tu2} = 380, C_{Tu3} = 2/5 \quad Tu \le 1\% \qquad (13a)$$

$$C_{Tu1} = 37.8, C_{Tu2} = 462.2, C_{Tu3} = 1.073 \quad Tu > 1\% \qquad (13b)$$

leading to the modified correlation presented in Fig. 1. From Fig. 1 it can be seen that the modified constants for $Tu > 1\%$ were selected such that the correlation closely resembles the ones by Menter et al. [14] and Langtry & Menter [12] as they are proven to work well for bypass transition.

The coupling of the transition and turbulence model is achieved in the same way as previously done by Nichols [40] or Piotrowski &

Zingg [16] by multiplying production and destruction terms in the SA-neg model [36] with the intermittency function. Hence, the transport equation for the modified turbulent viscosity $\tilde{\nu}$ reads

$$\partial_t (\rho \tilde{\nu}) + \nabla \cdot (\rho \tilde{\nu} \mathbf{v}) = \gamma_s P_{\tilde{\nu}} - \max(\gamma_s, 0.1) \epsilon_{\tilde{\nu}} + \frac{1}{\sigma_{\tilde{\nu}}} [\nabla \cdot (\nu + \tilde{\nu}) \nabla \tilde{\nu} + c_{b2} \nabla \tilde{\nu} \nabla \tilde{\nu}] - \gamma_s C_5 \frac{\tilde{\nu}^2}{c_s^2} \nabla \tilde{\nu} \nabla \tilde{\nu}. \qquad (14)$$

It is important to note that Nichols [40] found that the production term of the SA model $P_{\tilde{\nu}}$ needs to be exactly zero in laminar regions. By default, this is not the case for Eq. (1), since it has a theoretical limit of $\tilde{\gamma} = \ln(1/c_{e2})$ at the wall as most of the intermittency transport models [14]. To solve this issue Nichols [40] proposed a modified function $\gamma_s$ defined by

$$\gamma_s = \frac{\min(\gamma, 1) - 1/c_{e2}}{1 - 1/c_{e2}} \qquad (15a)$$

$$\gamma_s = \max(\min(\gamma_s, 1), 0) \qquad (15b)$$

which is utilized here as well. Keeping in mind, that in the present case $\gamma = e^{\tilde{\gamma}}$. The limiting of $\gamma_s$ in front of the destruction term $\epsilon_{\tilde{\nu}}$ ensures remaining destruction of turbulent viscosity even if $\gamma$ is very close to zero. Here the limiting value 0.1 is used as applied by most of the authors [12,14,40]. It should be noted that Coder and Maughmer [17] suggested a value of 0.5, but without Nichols [40] correction Eq. (15). Throughout this work, it was found that a limiting value of 0.1 gives a reliable behavior in combination with Nichols correction term Eq. (15). Additionally, since in this work transonic, separated flows are considered the last term on the right-hand side of Eq. (14) includes the compressible mixing layer correction [42]. Since this term should only be active in turbulent separation regions, it is multiplied by $\gamma_s$ as well. The constant $C_5 = 3.5$ is kept unchanged with respect to the default value [42]. The eddy viscosity $\mu_t$ follows from the modified eddy viscosity $\tilde{\nu}$ as

$$\mu_t = f_{v1} \mu \psi \qquad (16)$$

where $f_{v1} = \psi^3/(\psi^3 + c_{v1}^3)$. The limiter $\psi = \psi(\chi)$ is a function of, $\chi = \rho\tilde{\nu}/\mu$ similar to the limiter defined in [37]. All instances of the max/min functions in the source terms are replaced with the logarithmic P-mean function, as proposed by Piotrowski & Zingg [16], to eliminate discontinuities and ensure differentiability. The logarithmic P-mean function [16] for a max($x_1, x_2$) operation with two arguments is defined as:

$$\max(x_1, x_2) \approx \varphi_p(x_1, x_2) = x_1 + \frac{\log(1 + \exp(P(x_2 - x_1)))}{P} \qquad (17)$$

where the parameter $P$ controls the degree of smoothing around the discontinuity. In this work, the default value $P = 300$ is used, as suggested by Piotrowski & Zingg [16]. The proposed transition model is calibrated, such that free stream values for the modified turbulent viscosity of $\tilde{\nu}/\nu = 0.05 - 0.1$ should be applied.

### 3. Discretization

In this work, a higher-order discontinuous Galerkin (DG) discretization is adopted to discretize the system of partial differential equations. Therefore, we gather all equations in a vectorial form as

$$\partial_t q + \nabla \cdot f_c(q) - \nabla \cdot f_v(q, \nabla q) - s(q, \nabla q) = 0, \qquad (18)$$

where

$$f_c(q) = \begin{pmatrix} \rho \tilde{v}^T \\ \rho \tilde{v} \tilde{v}^T + \rho I \\ (\rho E + \rho) \tilde{v}^T \\ \rho \tilde{\nu} \tilde{v}^T \\ \rho \tilde{\gamma} \tilde{v}^T \end{pmatrix} \qquad (19)$$
contains the convective fluxes and

$$f_v(q, \nabla q) = \begin{pmatrix} \mathbf{0} \\ \tau + \tau_t \\ \nu^T (\tau + \tau_t) - \Phi - \Phi_t \\ \frac{1}{\sigma_v} (\mu + \rho \bar{v}) \nabla \bar{v}^T \\ \left( \mu + \frac{\mu_t}{\sigma_y} \right) \nabla \bar{\gamma}^T \end{pmatrix} \quad (20)$$

contains the viscous fluxes. Here $\tau = \mu(-\frac{2}{3}(\nabla \cdot \nu) + 2S)$ is the viscous stress tensor and $\tau_t = \mu_t(-\frac{2}{3}(\nabla \cdot \nu) + 2S)$ the turbulent stress tensor, where $S = (\nabla \nu + \nabla \nu^T)/2$ describes the deformation tensor. The lami-nar $\Phi = -\mu C_p/\text{Pr}\nabla T$ and turbulent $\Phi_t = -\mu_t C_p/\text{Pr}\nabla T$ heat flux vectors are given with a Prandtl number $\text{Pr} = 0.72$ and a turbulent Prandtl number $\text{Pr}_t = 0.9$. The vector of conserved variables $q$ is given by

$$q = (\rho, \rho v^T, \rho E, \rho \bar{v}, \rho \bar{\gamma})^T, \quad (21)$$

where $E = e + v^2/2$ represents the total energy and $e$ denotes the internal energy. The pressure is defined by a polytropic ideal gas law: $p = (\kappa_{PG}-1)\rho e$. For a DG method we define a discrete weak formulation of Eq. (18). We divide the physical domain $\Omega$ into $N_e$ non-overlapping elements $\kappa$ forming the computational mesh $\Omega_h$. Inside element $\kappa$ we search an approximate solution in the function space of discontinuous polynomials $V_h^p = \{\phi \in L^2(\Omega_h) : \phi|_{\kappa} \in P^p \text{ for } \kappa \in \Omega_h\}$ as

$$q_h(x, t) = \sum_{l=1}^{N_p} \phi_{\kappa}^{l}(x) Q_{\kappa}^{l}(t) \quad \forall x \in \kappa, \kappa \in \Omega_h, \quad (22)$$

where $Q_{\kappa}^l$ are the $N_p = (\rho+1)(\rho+2)/2$ degrees of freedom inside element $\kappa$ in two space dimensions, while the $\phi_{\kappa}^l$ constitute an orthonormal basis [43] of $V_h^p$ inside $\kappa$. Multiplication of Eq. (18) by test functions $v_h$ of the same space $V_h^p$ and integration by part leaves us with

$$\int_{\Omega_h} v_h \partial_t q_h dV + \mathcal{F}_c(q_h, v_h) + \mathcal{F}_v(q_h, v_h) + \mathcal{F}_s(q_h, v_h) = 0, \quad (23)$$

where $\mathcal{F}_c$, $\mathcal{F}_v$ and $\mathcal{F}_s$ are operators gathering the discretization of convective, viscous and source terms respectively. The discretization of convective fluxes is given by

$$\mathcal{F}_c(q_h, v_h) = - \int_{\Omega_h} f_c(q_h) \cdot \nabla v_h dV + \int_{\partial \kappa} \left( \left[ \left[ v_h \right] \right] h_c(q_h^+, \bar{q}_h^-, n) \right) dS + \int_{\partial \Omega_h} v_h^+ f_c(q_h(q_h^+, n)) \cdot n dS, \quad (24)$$

where $\left[ \left[ \psi \right] \right] = \psi^+ - \psi^-$ denotes a jump operator of any quantity $\psi$ along an interior interface $\partial \kappa$ between master element $\kappa^+$ and neighboring element $\kappa^-$. The normal vector pointing outward with respect to $\kappa^+$ is expressed as $n$. For the convective flux function $h_c$ we use a Roe flux [44] incorporating Harten's entropy fix [45]. The operator $q_h$ allows setting boundary conditions such as the non-reflective Riemann farfield boundary condition and the no-slip and adiabatic condition at walls. The viscous term is discretized by utilizing Bassi & Rebay's BR2 method [38]

$$\begin{align} \mathcal{F}_v(q_h, v_h) = & \int_{\Omega_h} f_v(q_h, \nabla q_h + \mathbf{R}_h) \cdot \nabla v_h dV \nonumber \\ & - \int_{\partial \kappa} \left( \left[ v_h \right] h_v(q_h^+, \bar{q}_h^-, \nabla q_h^+, \nabla \bar{q}_h^-) \cdot n \right) dS \tag{25} \\ & - \int_{\partial \Omega_h} v_h^+ f_v(q_h(q_h^+, n), \nabla q_h(q_h^+, n) + \eta_{BR2} \mathbf{r}_h^e) \cdot n dS \nonumber \end{align}$$

where the operator $\nabla q_h$ allows applying Neumann type boundary conditions and $\eta_{BR2}$ is a stabilization parameter. Throughout this work, a value of $\eta_{BR2} = 2$ has been applied. The numerical viscous face flux function $h_v$ is taken as the average between two adjacent elements,

$$h_v(q_h^+, \bar{q}_h^-, \nabla q_h^+, \nabla \bar{q}_h^-, n) = \{f_v(q_h, \nabla q_h + \eta_{BR2} \mathbf{r}_h^e)\} \cdot n, \quad (26)$$

<figure><img src="image_5.png" /><figcaption>Fig. 2. Specification of geometry and applied boundary conditions for flat plates cases. Displayed mesh corresponds to coarse level for T3A case.</figcaption></figure>

where $\{\psi\} = 1/2 (\psi^+ + \psi^-)$ denotes the average of an arbitrary quantity across an internal face $\partial \kappa$. The global lifting operator $\mathbf{R}_h$ is defined as the sum of the local lifting operator $\mathbf{r}_h^e$

$$\mathbf{R}_h = \sum_{e \in \partial \kappa} \mathbf{r}_h^e \quad \forall \kappa \in \Omega_h, \quad (27)$$

where the local lifting operator is the solution of

$$\int_{\kappa \cup \kappa^{-}} v_h \mathbf{r}_h^e dV = - \int_{\partial \kappa} \{v_h\} \left[ \left[ q_h \right] \right] \mathbf{n}^T dS, \quad (28)$$

defined on internal faces. The source term is discretized in a similar way by utilizing the BR2 method [38]

$$\mathcal{F}_s(q_h, v_h) = \int_{\Omega_h} v_h s(q_h, \nabla q_h + \mathbf{R}_h) dV. \quad (29)$$

Specifically, the gradients within the source term are computed by incorporating the global lifting operator as suggested by Bassi et al. [38]. Finally, we can cast (23) into matrix form as

$$\mathbf{M} \partial_t \mathbf{Q} = \mathbf{R}(\mathbf{Q}), \quad (30)$$

where $\mathbf{M}$ is a diagonal mass matrix by applying a Gram-Schmidt orthonormalization [43] to the DG basis in (22), while $\mathbf{R}(\mathbf{Q})$ is a vector of discrete residuals gathering all convective, viscous and source terms. Regarding unsteady computations Eq. (30) is advanced in time using a second-order implicit backward differencing scheme (BE-BDF2), while steady state solutions are obtained by employing a quasi-Newton approach with local time stepping in combination with a transient continuation technique [46]. Resulting linear systems are solved by a restarted GMRES method utilizing a block ILU(0) preconditioner. The linear operator or Jacobian matrix is obtained by carefully analytically linearizing the discrete residuals including all model terms.

## 4. Results: steady solutions

This section is dedicated to the transition model validation for steady state solutions. Therefore, several test cases from the literature are considered. To start with, the canonical zero pressure gradient flat plate cases Schubauer & Klebanoff (SK) [47] and ERCOFTAC T3A, T3B, T3A- [48] are used to show the model capabilities for T-S and bypass transition at different turbulence intensities and Reynolds numbers. In the following, the subsonic NLF0414F airfoil [49] case at high chord Reynolds number, is used to investigate natural transition at low turbulence environment under the influence of favorable and adverse pressure gradients. Finally, to access the models performance for transonic flows, subject to weak shocks and varying pressure gradients, the CAST10-2 airfoil is investigated.
<table><caption>**Table 1**
Flow conditions for the zero pressure gradient flat plate cases SK [47], T3A, T3B and T3A- [48].</caption><thead><tr><th></th><th>SK</th><th>T3A</th><th>T3B</th><th>T3A-</th></tr></thead><tbody><tr><td><i>Re</i></td><td>3.36 · 10<sup>6</sup></td><td>3.7 · 10<sup>5</sup></td><td>6.3 · 10<sup>5</sup></td><td>1.32 · 10<sup>6</sup></td></tr><tr><td>Tu[%]</td><td>0.18</td><td>3.0</td><td>6.0</td><td>0.6</td></tr></tbody></table>

### 4.1. Zero pressure gradient flat plate

To validate the transition model for zero pressure gradient flows, the Schubauer & Klebanoff (SK) [47] and ERCOFTAC T3A, T3B, T3A- [48] cases are simulated. The cases offer a wide range of Reynolds numbers based on the length of the plate and turbulence intensities, as can be seen from Table 1. The higher turbulence intensity for the ERCOFTAC cases is leading to bypass transition, whereas the low turbulence intensity SK case shows natural (T-S) transition. A full spacial convergence analysis is performed for the SK and the T3A case. For both cases, solutions for polynomial degrees between $p1 - 4$ are computed on three consecutively refined meshes. For the T3A case, quadrilateral meshes of 1551 (coarse), 4524 (medium) and 13593 (fine) elements are used. Since the SK case has a much higher Reynolds number, slightly finer meshes of 4356 (coarse), 9018 (medium) and 21420 (fine) elements are utilized. Each grid consist of linear straight sided quadrilateral elements. For both the SK and T3A cases, the mesh levels are designed to achieve specific $\Delta y^+$ values: coarse $\Delta y^+ = 5$, medium $\Delta y^+ = 3$, and fine $\Delta y^+ = 1$. In the context of the current DG method, these values are scaled by the polynomial degree $p$ as $\Delta y_p^+ = \Delta y^+/(p+1)$. Similarly, the wall tangential spacing is scaled as $\Delta x_p^+ = \Delta x^+/(p+1)$. The scaled wall-normal and tangential mesh spacings near the transition locations are provided in Table 2. The general domain geometry as well as the applied boundary conditions are the same for all cases, as displayed in Fig. 2. Along the plate, an adiabatic no-slip wall boundary condition is applied. Subsonic inlet condition in combination with a static outlet pressure specify a low Mach number of $M = 0.2$. On the top side of the 2D domain, a non-reflective Riemann farfield condition is employed. Upstream of the plate, a short symmetry patch is employed. In Fig. 3 skin friction distributions $C_f$ are displayed for the SK case for each mesh up to polynomial degree $p4$. On the fine mesh starting from $p2$ the results can be considered fully converged. In comparison to the experiment, the SA-neg-log($\gamma$) shows good agreement. When comparing to the $\gamma$-SST model by Menter et al. [14], the SA-neg-log($\gamma$) model predicts transition earlier, matching better with the transition onset from the experimental data. However, with both models the predicted transition is quicker than in the experiment, where Menter et al.'s $\gamma$-SST model [14] transitions more towards the end of the transition process. Hysteresis effects were tested for the SK case by using both laminar and fully turbulent (pure SA model) initial conditions. The results showed that no differing solutions emerged, as both initial conditions converged to the same final solution.

The T3A case displayed in Fig. 4 is much less demanding in terms of spacial resolution compared to the SK case, because of the much lower Reynolds number of the T3A case. The results are already fully spatially converged for $p2$ on the coarse mesh by only 9306 DoFs per equation. Lorini et al. [29] found similar convergence properties for the T3A case by their high-order DG $k_L - k - \log(\omega)$ model. The results using the SA-neg-log($\gamma$) model compare well against the experiment, giving similar behavior as the $\gamma$-SST model by Menter et al. [14]. In Fig. 5 the skin friction distribution for the T3A case is displayed without using the modified coefficients (Eq. (12)). It can be seen that transition is erroneously predicted too far downstream, which highlights the necessity for the modified coefficients in case of bypass transition. In Fig. 6 the skin friction distributions for the T3B case are displayed computed on the medium mesh for $p1 - 4$. It can be seen that the transition location is well predicted compared to the experiment. Furthermore, the SA-neg-log($\gamma$) model predicts the jump between minimum and maximum

<table><caption>**Table 2**
Mesh parameters for SK and T3A case for the different mesh resolutions and polynomial degrees. The scaled parameters follow from: $\Delta y_p^+ = \Delta y^+/(p+1)$ and $\Delta x_p^+ = \Delta x^+/(p+1)$.</caption><thead><tr><th></th><th></th><th>Δy<sub>p</sub><sup>+</sup></th><th>Δx<sub>p</sub><sup>+</sup></th></tr></thead><tbody><tr><td rowspan="4">Coarse<br>Δy<sup>+</sup> = 5</td><td>p1</td><td>2.50</td><td>400.00</td></tr><tr><td>p2</td><td>1.67</td><td>266.67</td></tr><tr><td>p3</td><td>1.25</td><td>200.00</td></tr><tr><td>p4</td><td>1.00</td><td>160.00</td></tr><tr><td rowspan="4">Medium<br>Δy<sup>+</sup> = 3</td><td>p1</td><td>1.50</td><td>200.00</td></tr><tr><td>p2</td><td>1.00</td><td>133.33</td></tr><tr><td>p3</td><td>0.75</td><td>100.00</td></tr><tr><td>p4</td><td>0.60</td><td>80.00</td></tr><tr><td rowspan="4">Fine<br>Δy<sup>+</sup> = 1</td><td>p1</td><td>0.50</td><td>75.00</td></tr><tr><td>p2</td><td>0.33</td><td>50.00</td></tr><tr><td>p3</td><td>0.25</td><td>37.50</td></tr><tr><td>p4</td><td>0.20</td><td>30.00</td></tr></tbody></table>

$C_f$ better than the $\gamma$-SST model by Menter et al. [14]. In Fig. 7 the skin friction distribution is displayed for the T3A- case. For this case, the SA-neg-log($\gamma$) model predicts transition too early as well as too steep. The $\gamma$-SST model by Menter et al. shows a similar behavior and is attributed to the fact that the transition length is a simple model constant in both models [14]. Therefore, such a slow transition process is not well captured by both models, whereas the $\gamma - Re_0$-SST model [8,12] captures this flow better since it uses a correlation for the transition length.

### 4.2. NLF0414F airfoil

The NLF0414F airfoil case is selected to test the ability of the SA-neg-log($\gamma$) model to predict natural transition at high Reynolds number under the influence of different adverse and favorable pressure gradient effects. McGhee et al. [49] performed multiple experiments at various Reynolds and Mach numbers. Here, the experimental run 15 is selected at a Mach number of $M = 0.4$ and a chord Reynolds number of $9.5 \cdot 10^6$. McGhee et al. [49] gave a low turbulence intensity of $Tu \approx 0.08$, which has a slight dependency on the Reynolds and Mach number. In order to compare the present results also with the ones by Fehrs [34] we set $Tu = 0.08$. This test case is characterized by a so-called laminar drag bucket, which is a sudden decrease of drag for a certain range of angles of attack. Typically, a laminar drag bucket is formed, due to a large change in transition location. More specifically, in the present case, the transition locations change by up to 70% with respect to the chord length. Since the transition in case of laminar drag buckets is usually quite sensitive, those cases are very well suited to test transition models. It is important to note that the original experiments were carried out on an NLF0414F airfoil model, which trailing edge is below the center line. An updated version of this airfoil is publicly available in the UIUC database, where a correction has been applied such that the trailing edge lies on the center line, as can be seen in Fig. 8. The original coordinates given by the experiment were found to be insufficient to create properly defined high-order curved meshes. Hence, the UIUC coordinates were used, by applying a rotation correction of $\Delta\alpha_{corr} = -1.512^\circ$. From Fig. 8 can be seen, that when applying this correction both geometries are the same. In the following, we always refer to angles of attack with respect to the UIUC airfoil. However, for comparison with the experiment [49] the rotation correction is taken into account. Three consecutively refined CH-topology meshes with 54 821 (coarse), 64 884 (medium) and 76 518 (fine) elements are used to prove spatial convergence. Each mesh consists of 2nd order curved quadrilaterals to accurately represent the airfoil shape. The medium-resolution mesh is shown in Fig. 9, with the CH topology specifically designed to accurately capture the wake behind the blunt trailing edge. The non-dimensional, scaled wall normal and tangential mesh spacings for the three mesh levels are given in Table 3. The given tangential spacing represents values in the vicinity of the transition location. The
<figure><img src="image_6.png" /><figcaption>Fig. 3. Schubauer & Klebanoff [47] flat plate test case. Mesh and polynomial degree sensitivity analysis for skin friction coefficient $C_f$. Comparison with experimental data [47] and $\gamma$-SST results by Menter et al. [14].</figcaption></figure>

<figure><img src="image_7.png" /><figcaption>Fig. 4. ERCOFTAC T3A flat plate test case. Mesh and polynomial degree sensitivity analysis for skin friction coefficient $C_f$. Comparison with experimental data and $\gamma$-SST results by Menter et al. [14].</figcaption></figure>

<figure><img src="image_8.png" /><figcaption>Fig. 5. ERCOFTAC T3A flat plate test case. Skin friction coefficient $C_f$ using Fehrs [34] correlation with (Eq. (13)) and without (Eq. (12)) correction for higher turbulence intensities $Tu > 1\%$.</figcaption></figure>

computational domain extends 100 chord lengths in every direction, and a non-reflective farfield condition is applied at the outer boundary. In Fig. 10 the skin friction distributions are displayed on each mesh for polynomial degrees $p1-3$. At this angle of attack, the flow on the suction side stays laminar until $x/c = 0.71$, whereas on the pressure side a tiny suction peak at the leading edge is just strong enough to cause transition. The flow is spatially converged from $p2$ on the medium mesh corresponding to 389 304 DoFs per equation, which shows how demanding this high Reynolds number flow is. On the coarse mesh, no spatial convergence is achieved for $p1$ and $p2$. For a polynomial degree of $p3$ the flow is converged on all three meshes. In Fig. 11 the convergence history of the steady state quasi-Newton method is presented for the $\alpha = -1.525^{\circ}$ case. It can be seen that even for these stiff conditions, posed by the high Reynolds number and low turbulence intensity, the solver converges quite fast and deep. During the last stage of convergence, almost quadratic convergence rate is achieved, a key
<figure><img src="image_9.png" /><figcaption>**Fig. 6.** ERCOFTAC T3B flat plate test case. Skin friction coefficient $C_f$ on medium mesh for $p1-4$. Comparison with experimental data and $\gamma$-SST results by Menter et al. [14].</figcaption></figure>

<figure><img src="image_10.png" /><figcaption>**Fig. 7.** ERCOFTAC T3AM flat plate test case. Skin friction coefficient $C_f$ on medium mesh for $p1-4$. Comparison with experimental data and $\gamma$-SST results by Menter et al. [14].</figcaption></figure>

feature of the quasi-Newton method. In Fig. 12(a) the aforementioned laminar drag bucket is visible in the $C_D-C_L$ curve. It can be seen that the SA-neg-log($\gamma$) model captures the laminar drag bucket well compared to the experiment. Overall for this case the SA-neg-log($\gamma$) model gives a similar performance as the $\gamma$-SST model by Fehrs [34], which mainly attributes to the same $Re_{\theta_L}$ correlation used in both models. The minimum drag inside the laminar drag bucket is caused by a vanishing suction peak at the leading edge for certain angles of attack, leading to delayed transition on both sides of the airfoil. From Fig. 13 it can be seen that for an angle of attack in the center of the laminar drag bucket the flow on the suction side stays laminar until $x/c = 0.58$ and on the pressure side until $x/c = 0.65$. The lift curves in Fig. 12(b) are generally matching well between experiment and the present computations. However, for larger angles of attack $\alpha > 4^\circ$ discrepancies are visible. It is not entirely clear where those differences are stemming from, since the drag values are comparing well, indicating that an erroneous predicted transition is not the origin. One aspect that was found to be very important to accurately predict the transition for this flow is the order of the mesh. The effect of high-order mesh curving on the transition prediction inside the laminar drag bucket can be seen from Fig. 14. Therefore, simulations were performed for minimum drag configuration at $\alpha = 0.542^\circ$ by using coarse and medium mesh resolutions, each one with either quadratic curved or linear non-curved quadrilateral elements. Apart from the mesh curving, the meshes are exactly the same. Both the coarse and medium mesh in quadratic curved form at a polynomial degree of $p3$ converge to the same late transition on both sides of the airfoil. Contrary, the linear meshes do not converge to a single transition location. This transition behavior is considered to be erroneous, since the predicted drag values on the linear meshes are much larger than those of the experiment and the curved meshes. On the pressure side the differences are largest, where the transition location using the linear meshes occurs upstream at $x/c = 0.05$ instead of $x/c = 0.58$ using the curved meshes. On the suction side, the differences are smaller, but still amount to about 10% on the medium

mesh. The cause of this behavior was traced down to local curvature oscillations in case of the linear meshes. Since the laminar boundary layer at this high Reynolds number is very sensitive, already small geometrical imperfections can be enough to trigger transition much earlier. Moreover, this can not only be attributed to a single linear mesh generation software, as this was tested with linear meshes generated using GMSH and Pointwise software. The high sensitivity of this case was also reported by Langel et al. [50]. A more detailed discussion can be found in Appendix. Additionally, it was recognized that the usage of linear meshes drastically slowed down the convergence rate of the steady state quasi-Newton solver. In order to avoid this and create high quality curved airfoil meshes the airfoil geometry coordinates were re-splined using a high-order polynomial and based on this improved geometry the mesh was curved using GMSH. In Fig. 15, local Mach number fields are presented for $\alpha = -1.525^\circ$ using the medium mesh resolution and polynomial degrees $p1$ to $p3$. No oscillations are observed in the subsonic flow. Similarly, Fig. 16 shows the turbulent variable field, normalized by the free-stream viscosity, with a close-up view near the downstream boundary layer edge for each polynomial degree. Weak spurious oscillations are visible at $p1$ (Fig. 16(a)), but these artifacts diminish with higher polynomial degrees, as seen in $p3$ (Fig. 16(c)), where the oscillations have practically disappeared. Furthermore, since the eddy viscosity is limited according to Eq. (16), these artifacts do not affect the rest of the solution. In Fig. 17, a close-up view of the transition variable $\tilde{\gamma}$ is shown near the suction side transition location for the same configuration, with increasing polynomial degrees. At $p1$ (Fig. 17(a)), a slight overproduction of $\tilde{\gamma}$ occurs, reaching values slightly above zero just after the transition. This is not unexpected, as the transition model formulation used here does not explicitly prevent such overshoots. However, for higher polynomial degrees (Figs. 17(b) and 17(c)), the overshoots diminish significantly and nearly vanish. Additionally, since the limiter in Eq. (15) is applied, these overshoots do not affect the rest of the solution. No impact on the steady-state convergence has been observed as a result of these weak overshoots.

<figure><img src="image_11.png" /><figcaption>**Fig. 8.** Illustration of the angle of attack correction applied to the NLF0414F airfoil from the UIUC database to match the original geometry by McGhee et al. [49].</figcaption></figure>

<figure><img src="image_12.png" /><figcaption>**Fig. 9.** Computational CH-type mesh (medium) for the NLF0414F airfoil. The outer farfield boundary condition is of non-reflective Riemann type and the airfoil surface uses an adiabatic no-slip condition.</figcaption></figure>
<figure><img src="image_13.png" /><figcaption>**Fig. 10.** NLF0414F airfoil case at $\alpha = -1.525^{\circ}$. Mesh and polynomial degree sensitivity analysis for skin friction coefficient $C_f$. Solid lines represent suction side distribution and dashed lines pressure side distributions.</figcaption></figure>

**Table 3**

Mesh parameters for the NLF0414F case for the different mesh resolutions and polynomial degrees. The scaled parameters follow from: $\Delta y_p^+ = \Delta y^+/(p+1)$ and $\Delta x_p^+ = \Delta x^+/(p+1)$.

<table><thead><tr><th></th><th></th><th>Δy<sub>p</sub><sup>+</sup></th><th>Δx<sub>p</sub><sup>+</sup></th></tr></thead><tbody><tr><td rowspan="4">Coarse</td><td>p1</td><td>1.50</td><td>500.00</td></tr><tr><td>p2</td><td>1.00</td><td>333.33</td></tr><tr><td>p3</td><td>0.75</td><td>250.00</td></tr><tr><td>p4</td><td>0.60</td><td>200.00</td></tr><tr><td rowspan="3">Medium</td><td>p1</td><td>1.00</td><td>375.00</td></tr><tr><td>p2</td><td>0.67</td><td>250.00</td></tr><tr><td>p3</td><td>0.50</td><td>187.50</td></tr><tr><td rowspan="3">Fine</td><td>p1</td><td>0.50</td><td>250.00</td></tr><tr><td>p2</td><td>0.33</td><td>83.30</td></tr><tr><td>p3</td><td>0.25</td><td>62.50</td></tr></tbody></table>

<figure><img src="image_14.png" /><figcaption>**Fig. 11.** Convergence history of the DG SA-neg-log($\gamma$) model for the NLF0414F airfoil using the medium mesh at $Re = 9.5 \cdot 10^6$ and $\alpha = -1.525$ for polynomial degrees $p1-3$, where $p2$ was initialized with the prior $p1$ solution and $p3$ with the prior $p2$ solution.</figcaption></figure>

### 4.3. CAST10-2 airfoil

The next case is the transonic flow around the CAST10-2 airfoil at a Mach number of $M = 0.7$ and a Reynolds number of $2 \cdot 10^6$. This case is used to demonstrate the transition models capabilities for transonic flows involving weak shocks. Furthermore, again a laminar drag bucket is present for this flow, which is a suitable test for the transition models sensitivity in the transonic regime. The results will be compared against an experiment by Hebler [51] and computational results by Ströer et al. [52] and Fehrs [34]. The turbulence intensity was not exactly specified for the experiment, however values between $Tu = 0.2-0.4\%$ are suggested in the literature [20,52,53]. Therefore, in this work $Tu = 0.35$ is selected. Hebler [51] suggested a wind tunnel correction of $\Delta\alpha = -0.3^{\circ}$ and $\Delta M = -0.01$, which was applied by Ströer et al. [52]. At first, a mesh and polynomial sensitivity analysis is performed for an angle of attack of $\alpha = 1.077^{\circ}$, while applying the aforementioned correction. The selected domain extends 100 chords in all directions, and again a non-reflective condition is applied at the farfield. In Fig. 18 skin friction distributions can be seen for three CH meshes consisting of 27808 (coarse), 38752 (medium) and 53917 (fine) quadratic curved quadrilateral elements. The scaled wall mesh spacings for the three mesh levels are given in Table 4. The given tangential spacing represents values in the vicinity of the transition location. It can be seen that already on the coarse mesh utilizing a polynomial degree of $p2$ gives practically converged results, amounting to just 166848 DOFs. In Fig. 19 the pressure distribution at $\alpha = 1.077^{\circ}$ is compared against results by Ströer et al. [52] obtained by a compressible AHD (Arnal-Habiballah-Delcourt criterion [54]) transition model, which shows good agreement. At this angle of attack, the transition on the suction side is caused by an adverse pressure gradient induced by a weak shock around $x/c \approx 0.18$. On the pressure side, the flow is accelerated up to $x/c \approx 0.5$, followed by an adverse pressure gradient, leading to transition at $x/c = 0.6$. In Fig. 20 the predicted transition locations are displayed for several angles of attack, by comparing against the experiment [51] and the predictions by Ströer et al. [52]. Presented are predictions with and without the suggested wind tunnel correction. Overall, a good agreement with respect to the
<figure><img src="image_15.png" /><figcaption>**Fig. 12.** Drag and lift polar for the NLF0414F airfoil. Comparison against experimental data by McGee et al. [49] and γ-SST results by Fehrs [34].</figcaption></figure>

<figure><img src="image_16.png" /><figcaption>**Fig. 13.** Skin friction distribution for NLF0414F airfoil at $\alpha = 0.542^{\circ}$ showing its minimum drag configuration inside the laminar drag bucket. Solid lines represent suction side distribution and dashed lines pressure side distributions.</figcaption></figure>

<figure><img src="image_17.png" /><figcaption>**Fig. 14.** Skin friction distribution for NLF0414F airfoil at $\alpha = 0.542^{\circ}$. Effect of high-order mesh curving on transition onset. Solid lines represent suction side distribution and dashed lines pressure side distributions.</figcaption></figure>

experiment is achieved. Only the exact angle of attack at which the transition location suddenly jumps is slightly shifted compared to the experiment. The computational results by Ströer et al. [52] show a similar slight shift with respect to the transition location jump. The transition location on the pressure side is practically unaffected for this range of angles of attack and stays at $x/c \approx 0.6$. The sudden transition change around $\alpha \approx 0^{\circ}$ leads to a formation of a laminar drag bucket and a non-linear lift plateau, which can be seen in Fig. 21. The drag curve Fig. 21(a) shows good overall agreement against the experiment. Especially the laminar drag bucket is well predicted. Furthermore, the lift plateau around $\alpha \approx 0^{\circ}$ compares generally well between experiment and present computations. At higher angles of attack, further away from the laminar drag bucket, slightly larger deviations for the lift values are observed. Results by Piotrowski & Zingg [20] with the 3-equation transition model SA-sLM2015 show lift deviations of similar

magnitude compared to the experiments, although at a slightly larger Mach number. For this case, the transition locations were found to be much less sensitive to high-order mesh curving than for the NLF0414F case from Section 4.2. In fact, even linear meshes were able to reproduce the here presented results. This is most likely attributed to the fact, that the Reynolds number for the CAST10-2 case is more than four times lower than for the NLF0414F from the previous section. In Fig. 22, the local Mach number fields for polynomial degrees $p1-p3$ on the medium mesh at $\alpha = 1.077^{\circ}$ ($\alpha_{corr} = 0.777$) are shown. A small “shocklet” forms on the suction side above the transition location at $x/c \approx 0.18$, which is resolved across all polynomial degrees. Additionally, no solution oscillations are visible near the weak shock wave in the Mach number field for any polynomial degree. In Fig. 23, the turbulent variable $\bar{v}$ is displayed for the same configuration and polynomial degrees on the medium mesh. In general, no prominent oscillations are present in the turbulent solution. However, a slight oscillation around the boundary layer edge’s downstream part can be observed for the lowest polynomial degree $p1$ in Fig. 22(a). As the polynomial degree increases, these oscillations disappear, as seen for $p3$ in Fig. 23(c), where they are almost entirely removed. Due to the eddy viscosity limiter Eq. (16), these oscillations do not negatively affect the overall solution. Furthermore, no significant impact on steady-state convergence was observed. In Fig. 24, the solution field for the logarithmic transition variable $\bar{\gamma}$ is shown for polynomial degrees $p1-p3$ on the medium mesh, with a close-up of the region near the transition location below the weak shock. For polynomial degree $p1$, as seen in Fig. 23(a), a slight oscillation along the boundary layer edge is visible. However, similar to the turbulent variable, these oscillations diminish with increasing polynomial degree, as shown for $p2$ and $p3$ in Figs. 24(b) and 24(c). This behavior may be attributed to the model formulation, which preserves positivity but is not explicitly designed to prevent overshoots as already mentioned for the NLF0414F case in the previous section. These overshoots can occur in regions of rapid transition, such as beneath the shock wave in the current configuration, especially in spatially under-resolved simulations. However, due to the limiter in Eq. (15), these minor oscillations in under-resolved cases do not affect the overall solution. Additionally, no significant impact on steady-state convergence was observed in the presence of such weak overshoots.

## 5. Results: steady and unsteady transonic OALT25 airfoil

The last case is the transonic flow around the supercritical laminar airfoil ONERA OALT25. This case is introduced to demonstrate the SA-neg-log($\gamma$) model’s capabilities to capture the flow transition in a laminar separation bubble occurring below a shock. It is important to
<figure><img src="image_18.png" /><figcaption>Fig. 15. Local Mach number fields at $\alpha = -1.525$ on the medium mesh resolution for the polynomial degrees $p1$ (a), $p2$ (b) and $p3$ (c).</figcaption></figure>

<caption>**Table 4**</caption>

Mesh parameters for the CAST10-2 case for the different mesh resolutions and polynomial degrees. The scaled parameters follow from: $\Delta y_p^+ = \Delta y^+/(p+1)$ and $\Delta x_p^+ = \Delta x^+/(p+1)$.

<table><thead><tr><th></th><th></th><th>Δy<sub>p</sub><sup>+</sup></th><th>Δx<sub>p</sub><sup>+</sup></th></tr></thead><tbody><tr><td rowspan="3">Coarse</td><td>p1</td><td>1.50</td><td>450.00</td></tr><tr><td>p2</td><td>1.00</td><td>300.00</td></tr><tr><td>p3</td><td>0.75</td><td>225.00</td></tr><tr><td rowspan="3">Medium</td><td>p1</td><td>1.00</td><td>300.00</td></tr><tr><td>p2</td><td>0.67</td><td>200.00</td></tr><tr><td>p3</td><td>0.50</td><td>150.00</td></tr><tr><td rowspan="3">Fine</td><td>p1</td><td>0.50</td><td>200.00</td></tr><tr><td>p2</td><td>0.33</td><td>133.33</td></tr><tr><td>p3</td><td>0.25</td><td>100.00</td></tr></tbody></table>

note that the SA-neg-log(γ) model does not incorporate a specific treatment for separation-induced transition. Instead, transition in separated flows is indirectly addressed through the correlation for $Re_{\theta c}$, which is calibrated against adverse pressure gradients. This approach is similar

to the methods used in the γ model by Menter et al. [14] and the γ-Reθ model [8,12]. The chord based Reynolds number is $Re = 500000$ at a fixed angle of attack of $\alpha = 4^{\circ}$. The results are compared against high resolved LES simulations by Zauner et al. [35]. A turbulence intensity of $Tu = 0.1\%$ is estimated, based on free stream fluctuations and comparison against several steady state results. However, a perfectly accurate value for the turbulence intensity in LES simulations is typically difficult to obtain, and therefore the value should be taken carefully. At first, a grid convergence study is performed for the steady flow at an angle of attack of $\alpha = 4^{\circ}$ and a free stream Mach number of $M = 0.67$. Secondly, unsteady URANS simulations are performed at a higher Mach number of $M = 0.735$, where the airfoil is experiencing transonic buffet [35]. Transonic buffet is an unsteady low frequency aerodynamic phenomenon occurring at certain flow conditions around airfoils and wings [55], where shock waves are periodically oscillating inline with flow separation. While the transonic buffet phenomenon is studied in depth for fully turbulent conditions using RANS models [55–57], this is not the case for free-transitional scenarios. To the author's best knowledge, this is the first time applying a RANS based transition model to free-transitional transonic buffet. All simulations are carried
<figure><img src="image_19.png" /><figcaption>**Fig. 16.** Turbulent variable $\bar{v}$ field normalized by the free-stream viscosity $\nu_{\infty}$ at $\alpha = -1.525^{\circ}$ on the medium mesh resolution for the polynomial degrees $p1$ (a), $p2$ (b) and $p3$ (c).</figcaption></figure>

out on a mesh consisting of 48 736 curved elements for polynomial degrees $p1$ – 3. The scaled wall-normal mesh spacings, based on the steady-state case at $M = 0.67$, are: $\Delta y_{p1}^{+} = 0.5$, $\Delta y_{p2}^{+} = 0.33$, and $\Delta y_{p3}^{+} = 0.25$. The corresponding scaled tangential spacings near transition locations, separation bubbles, and shock waves are approximately: $\Delta x_{p1}^{+} = 225$, $\Delta x_{p2}^{+} = 150$, and $\Delta x_{p3}^{+} = 112.5$. It should be noted here, that the LES domain is selected too small to achieve completely domain independent results. However, for reasons of comparability, the present domain size is selected exactly as for the LES [35].

### 5.1. Steady flow at $M = 0.67$

The convergence history of the steady state quasi-Newton solver is displayed in Fig. 25, proving that even for complicated transonic flows with separation, a robust, fast and deep convergence is achieved. In Fig. 26 the obtained pressure distributions alongside the skin friction distributions are presented, in comparison to mean LES results by

Zauner et al. [35]. The pressure distribution matches very well with the LES, except for small deviations behind the shock. When comparing the skin friction distribution, slightly larger differences are visible. In the laminar region, the RANS matches the LES results, by also showing the same separation locations. Furthermore, the transition locations are almost identical. Downstream of the transition location, larger discrepancies are visible. First of all the skin friction inside the main vortex region is predicted too low by the SA-neg-log($\gamma$) model. Additionally, the transition is slower compared to the LES, leading to a delayed reattachment. However, the error for the reattachment location with respect to the chord length is only about 3%. This is very much acceptable considering that most of the common transition models perform worse or at most similar for separation induced transition [58]. The main reason why most transition models encounter problems for laminar separation relates to an insufficient over production of turbulence inside the turbulent shear layer. The commonly observed result is a too low skin friction shortly after reattachment, as can be seen
<figure><img src="image_20.png" /><figcaption>Fig. 17. Logarithmic transition variable $\tilde{\gamma}$ solution at $\alpha = -1.525^{\circ}$ on the medium mesh resolution for the polynomial degrees $p1$ (a), $p2$ (b) and $p3$ (c).</figcaption></figure>

<figure><img src="image_21.png" /><figcaption>Fig. 18. CAST10-2 airfoil case at $\alpha = 1.077^{\circ}$ ($\alpha_{corr} = 0.777$). Mesh and polynomial degree sensitivity analysis for skin friction coefficient $C_f$. Solid lines represent suction side distribution and dashed lines pressure side distributions.</figcaption></figure>
<figure><img src="image_22.png" /><figcaption>**Fig. 19.** Pressure distribution $C_p$ for CAST10-2 airfoil case at $\alpha = 1.077^\circ$ ($\alpha_{corr} = 0.777^\circ$). Comparison against results by Ströer et al. [52]. Solid lines represent suction side distribution and dashed lines pressure side distributions.</figcaption></figure>

<figure><img src="image_23.png" /><figcaption>**Fig. 20.** Comparison of transition locations on the suction side of the CAST10-2 airfoil. Results are compared against experiments by Hebler [51] and computations by Ströer et al. [52] obtained by a compressible AHD transition model.</figcaption></figure>

in Fig. 26. In Fig. 27 a close up view of the resulting intermittency distribution $\gamma$ is displayed. It can be seen, that the transition model correctly starts to produce intermittency along the shear layer, above the laminar separation. At the apex of the separation bubble, rapid transition is occurring, when $\gamma$ quickly reaches values of $\gamma = 1$. Hence, Fig. 27 shows, that the transition process itself is predicted well by the model, however the subsequent production of turbulence by the turbulence model is too slow. The reason why the SA-neg-log($\gamma$) model still achieves quite reasonable results in that regard, is attributed to its ability to transition very rapidly between $\gamma = 0$ and $\gamma = 1$. Hence, even without any additional tweaks to the underlying turbulence model, [14, 58, 59] solid results can be obtained. The local Mach number fields, displayed in Fig. 28, are matching almost perfectly between RANS and LES. Nonetheless, a closer view of the $u$ velocity field in Fig. 29 shows, that the main differences are very localized around the laminar separation bubble. The maximum reverse flow velocity inside the main vortex is higher in case of the RANS prediction, as can already be anticipated through the skin friction distribution Fig. 26.

## 5.2. Unsteady transonic buffet at $M = 0.735$

Now the Mach number is increased to $M = 0.735$, while keeping all the remaining parameters fixed and unsteady URANS simulations are performed. For the unsteady $p3$ simulation presented in the following, additional shock-capturing was applied using the entropy-based method introduced by Guermond et al. [60] to stabilize the computation. At first in Fig. 30 it can be seen from the lift history, how a low frequency limit cycle oscillation is predicted at this Mach number for all polynomial degrees. Starting from $p2$ the results can be considered converged in space and time, when comparing the lift history in Fig. 30 as well as the mean skin friction and pressure distributions in Fig. 31. When comparing the lift histories in Fig. 30 against the LES

<table><caption>**Table 5**
Comparison of predicted low frequency buffet between SA-neg-log($\gamma$) model and a LES simulation by Zauner et al. [35].</caption><thead><tr><th></th><th>St</th><th>C&#773;<sub>L</sub></th><th>&Delta;C<sub>L</sub></th><th>min C<sub>L</sub></th><th>max C<sub>L</sub></th></tr></thead><tbody><tr><td>p1</td><td>0.072</td><td>0.718</td><td>0.317</td><td>0.542</td><td>0.859</td></tr><tr><td>p2</td><td>0.070</td><td>0.703</td><td>0.298</td><td>0.533</td><td>0.832</td></tr><tr><td>p3</td><td>0.070</td><td>0.700</td><td>0.296</td><td>0.535</td><td>0.831</td></tr><tr><td>LES [35]</td><td>0.082</td><td>0.752</td><td>0.39</td><td>0.53</td><td>0.92</td></tr></tbody></table>

results, it can be seen, that the minimum lift is quite well captured by the transition model. Controversially, the high lift phase is generally underpredicted. As a result, the overall lift variation $\Delta C_L$ is smaller in the URANS simulation. In Table 5 a direct comparison between LES and URANS can be seen for selected quantities of the predicted buffet instability. The Strouhal number $St = fc/U_\infty$, where c represents the chord length and f the buffet frequency, is slightly underpredicted by the present model. However, a slight dependency of the buffet frequency to the turbulence intensity was recognized. Since the selected turbulence intensity is only an estimation, a more accurate value could potentially further improve the frequency prediction. A comparison of the mean local Mach number fields is displayed in Fig. 32. While being of course only a qualitative measure, both fields show good overall agreement in terms of extent of the supersonic region and their shape. A more quantitative comparison is shown in Fig. 31 by means of the mean skin friction and pressure distributions. Although some discrepancies are visible, both are showing quite good agreement compared to the reference LES. In case of the LES higher frequency content is present in the lift signal. More precisely, an intermediate frequency separation bubble instability ($St = 0.393$) and a vortex shedding ($St = 2$) were identified by Zauner et al. [35]. In case of the present URANS simulations, no higher frequencies are present in the lift signal. Nonetheless, it should be noted here that a “shedding like” transient phenomenon is visible in the instantaneous flow fields. Even in case of the LES those higher frequencies are much less pronounced than the low frequency buffet. Therefore, it is possible that the higher frequencies in the URANS are not influencing the lift strongly enough to be visible. A more detailed analysis of those phenomena in case of the URANS is required, which is not the scope of the present work. Given the high complexity of this unsteady buffet case, it can be concluded that the SA-neg-log($\gamma$) model captures the main flow features to an acceptable degree. Moreover, considering the computational cost and time consumption of the highly resolved LES compared to the URANS, the obtained results are a promising indication that sufficiently accurate free-transitional buffet prediction is possible using cheaper models.

## 6. Conclusion

The present work has presented an implementation of a one equation local correlation based transition model into a high-order DG solver. This marks one of the first attempts to introduce the class of LCTM into a high-order numerical framework. The one equation intermittency equation in combination with a correlation function tailored for tranonic, high Reynolds number airfoil flows allows accurate prediction of natural transition subject to sensitive pressure gradients. A simple modification of the correlation constants is proposed to improve prediction capabilities of bypass transition. A special attention was paid to increase the model robustness. Therefore, several modifications are included to make the model positivity preserving, smooth and achieve full differentiability. Validation has been performed starting from zero pressure gradient flat plate cases for natural and bypass transition up to complex transonic airfoil test cases. Through mesh and polynomial degree convergence studies, performed for each case, spatial convergence of the model has been proven. Overall, for all cases, accurate transition prediction of the model is shown. Furthermore, it has been presented how sensitive high-order transition models are against proper
<figure><img src="image_24.png" /><figcaption>Fig. 21. Comparison of drag and lift polar for the CAST10-2 airfoil against experiments by Hebler [51] and results by Fehrs [34] obtained using a γ-SST model.</figcaption></figure>

<figure><img src="image_25.png" /><figcaption>Fig. 22. Local Mach number fields at $\alpha = 1.077^{\circ}$ ($\alpha_{corr} = 0.777$) on the medium mesh resolution for the polynomial degrees $p1$ (a), $p2$ (b) and $p3$ (c).</figcaption></figure>
<figure><img src="image_26.png" /><figcaption>**Fig. 23.** Turbulent variable $\bar{v}$ field normalized by the free-stream viscosity $\nu_\infty$ at $\alpha = 1.077^\circ$ ($\alpha_{corr} = 0.777$) on the medium mesh resolution for the polynomial degrees $p1$ (a), $p2$ (b) and $p3$ (c).</figcaption></figure>
<figure><img src="image_27.png" /><figcaption>**Fig. 24.** Logarithmic transition variable $\tilde{\gamma}$ solution at $\alpha = 1.077^{\circ}$ ($\alpha_{corr} = 0.777$) on the medium mesh resolution for the polynomial degrees $p1$ (a), $p2$ (b) and $p3$ (c). The black line indicates the sonic line $M = 1$.</figcaption></figure>
<figure><img src="image_28.png" /><figcaption>**Fig. 25.** Convergence history of the DG SA-neg-log(γ) model for the OALT25 airfoil at *M* = 0.67 for polynomial degrees *p*1 − 3, where *p*2 was initialized with the prior *p*1 solution and *p*3 with the *p*2 solution.</figcaption></figure>

<figure><img src="image_29.png" /><figcaption>**Fig. 26.** Spatial convergence of skin friction *C*<sub>f</sub> (dashed) and pressure distribution *C*<sub>p</sub> (solid) around the suction side of the OALT25 at *M* = 0.67. Comparison against mean LES results by Zauner et al. [35].</figcaption></figure>

<figure><img src="image_30.png" /><figcaption>**Fig. 27.** Close up view of the intermittency distribution γ = e^β in the vicinity of the laminar separation bubble for the OALT25 airfoil at M = 0.67.</figcaption></figure>

high-order mesh generation. In case of the NLF0414F airfoil at high Reynolds number, it has been revealed that linear meshes can give large errors in transition location compared to quadratic meshes. The superior convergence to steady state of the model in combination with a quasi-Newton solver is presented. Even for complex transonic and highly stiff problems, a rapid and deep convergence is achieved. Finally, the SA-neg-log(γ) model's ability to predict free-transitional transonic buffet is accessed. This marks one of the first attempts to predict transonic buffet in a free-transitional context by a RANS based transition model. Despite the high complexity of this unsteady flow,

the model accurately predicts the main flow features compared to a high resolved LES simulation. In a next step, the model's behavior for laminar separation prediction will be improved. Therefore, possible ways to increase turbulence production inside the turbulent shear layer will be addressed. Due to the high model accuracy and the ability to converge deeply even for complex cases, future work will be devoted to predict transonic buffet. Thus, the model will be used in combination with global linear stability analysis for buffet prediction.

### CRediT authorship contribution statement

**Matthias Plath:** Writing – original draft, Visualization, Validation, Software, Methodology, Conceptualization. **Florent Renac:** Writing – review & editing, Software, Methodology, Conceptualization. **Olivier Marquet:** Writing – review & editing, Methodology, Conceptualization. **Christian Tenaud:** Writing – review & editing, Methodology, Conceptualization.

### Declaration of competing interest

The authors declare that they have no known competing financial interests or personal relationships that could have appeared to influence the work reported in this paper.

### Acknowledgments

This work received funding from the European Union's Horizon 2020 research and innovation program under the Marie Skłodowska-Curie grant agreement No 955923. The authors extend special thanks to Markus Zauner and Pradeep Moise for generously sharing their data and expertise for the unsteady transonic buffet case around the OALT25 profile.

### Appendix. Effect of mesh curving on local curvature

We here illustrate the importance of using a high-order mesh representation for obtaining accurate transition modeling. In Fig. A.33 the signed local surface curvature at the surface Gauss integration points is displayed for the NLF0414F airfoil case (discussed in Section 4.2). Both cases displayed in Fig. A.33 have the same polynomial degree *p*3 and the same number of Gauss integration points. However, one case uses a second order curved mesh and the other one a simple linear mesh. It is obvious that the linear mesh case shows strong oscillations of local curvature at the leading edge Fig. A.33(a), while the second order mesh avoids most of the large amplitude oscillations. Further downstream, the oscillations decrease in amplitude for the linear case, as the overall curvature of the airfoil also changes much more gradually. Nonetheless, the oscillations still persist over the complete airfoil surface. Contrarily, the second order mesh almost completely avoids curvature oscillations apart from the leading edge. For the NLF0414F case from Section 4.2 the Reynolds number is quite high. Hence, the boundary layer is very sensitive, because the local momentum thickness Reynolds number *Re*<sub>θ</sub> is very close to criticality for the whole chord [50]. As a result, the curvature oscillations visible in Fig. A.33 can be strong enough to create small pressure gradients, which lead to early transition. In general, it has been found for the NLF0414F case at high Reynolds number, that the higher the polynomial degree of the solution approximation, the higher the sensitivity to geometrical imperfections.

### Data availability

Data will be made available on request.
<figure><img src="image_31.png" /><figcaption>**Fig. 28.** Comparison of local Mach number fields for p2 solution of SA-neg-log(γ) model and LES mean results by Zauner et al. [35] at a free stream Mach number of M = 0.67. The black solid line indicates the so-called sonic line where M = 1.</figcaption></figure>

<figure><img src="image_32.png" /><figcaption>**Fig. 29.** Close view of $u$ velocity fields in the vicinity of the laminar separation bubble for p2 solution of SA-neg-log(γ) model and LES mean results by Zauner et al. [35] at a free stream Mach number of M = 0.67. The black solid line indicates the so-called sonic line where M = 1 and the dashed line indicates the separation line.</figcaption></figure>

<figure><img src="image_33.png" /><figcaption>**Fig. 30.** Comparison of low frequency lift oscillation for an OALT25 airfoil experiencing transonic buffet at a Mach number of M = 0.735. Comparison against high-resolved LES results by Zauner et al. [35].</figcaption></figure>

<figure><img src="image_34.png" /><figcaption>**Fig. 31.** Comparison of mean skin friction $C_f$ (dashed) and mean pressure distribution $C_p$ (solid) around the suction side of the OALT25 at M = 0.735 experiencing strong transonic buffet. Comparison against mean LES results by Zauner et al. [35].</figcaption></figure>
<figure><img src="image_35.png" /><figcaption>Fig. 32. Comparison of mean local Mach number fields for URANS p2 solution of SA-neg-log(γ) model and LES results by Zauner et al. [35] at a free stream Mach number of M = 0.735 experiencing strong transonic buffet. The black solid line indicates the so-called sonic line where M = 1.</figcaption></figure>

<figure><img src="image_36.png" /><figcaption>Fig. A.33. Effect of mesh representation on local curvature for NLF0414F airfoil. Displayed is the signed local surface curvature $k_s$ at the Gauss integration points along the airfoil surface for a polynomial degree p3 using linear and curved meshes. Solid lines are indicating suction side values and dashed lines pressure side values.</figcaption></figure>

References

[1] Zauner M, De Tullio N, Sandham ND. Direct numerical simulations of transonic flow around an airfoil at moderate Reynolds numbers. AIAA J 2019;57(2):597–607. http://dx.doi.org/10.2514/1.J057335.

[2] van Ingen JL. A suggested semi-empirical method for the calculation of the boundary Layer Transition Region. Technische hogeschool deflt, vliegtuigbouwkunde, rapport VTH-74, 1956.

[3] Van Ingen JL. The eN method for transition prediction: Historical review of work at TU Delft. In: Proceedings of the 38th AIAA fluid dynamics conference and exhibit, Seattle, Washington, June 23-26, 2008, pp. 1–49; AIAA paper 2008-3830. AIAA; 2008.

[4] Krumbein A, Krimmelbein N, Schrauf G. Automatic transition prediction in hybrid flow solver, Part 2: Practical Application. J Aircr 2009;46(4):1191–9. http://dx.doi.org/10.2514/1.39738.

[5] Krumbein A, Krimmelbein N, Schrauf G. Automatic transition prediction in hybrid flow solver, Part 1: Methodology and sensitivities. J Aircr 2009;46(4):1176–90. http://dx.doi.org/10.2514/1.39736.

[6] Schrauf G. Large-scale laminar flow tests evaluated with linear stability theory. J Aircr 2004;41(2):224–30. http://dx.doi.org/10.2514/1.9280.

[7] Menter FR, Esch T, Kubacki S. Transition modelling based on local variables. In: Rodi W, Fueyo N, editors. Engineering turbulence modelling and experiments 5. Oxford: Elsevier Science Ltd; 2002, p. 555–64. http://dx.doi.org/10.1016/B978-008044114-6/50053-3.

[8] Menter FR, Langtry R, Völker S. Transition modelling for general purpose CFD codes. Flow Turbul Combust 2006;77(1):277–303. http://dx.doi.org/10.1007/s10494-006-9047-1.

[9] Walters DK, Leylek JH. A new model for boundary layer transition using a single-point RANS approach. J Turbomach 2004;126(1):193–202. http://dx.doi.org/10.1115/1.1622709.

[10] Walters DK, Cokljat D. A three-equation eddy-viscosity model for Reynolds-averaged Navier-Stokes Simulations of transitional flow. J Fluids Eng 2008;130(121401). http://dx.doi.org/10.1115/1.2979230.

[11] Mayle RE, Schulz A. Heat transfer committee and turbomachinery committee best paper of 1996 Award: The Path to predicting bypass transition. J Turbomach 1997;119(3):405–11. http://dx.doi.org/10.1115/1.2841138.

[12] Langtry RB, Menter FR. Correlation-based transition modeling for unstructured parallelized computational fluid dynamics codes. AIAA J 2009;47(12):2894–906. http://dx.doi.org/10.2514/1.42362.

[13] Durbin P. An intermittency model for bypass transition. Int J Heat Fluid Flow 2012;36:1–6. http://dx.doi.org/10.1016/j.ijheatfluidflow.2012.03.001.

[14] Menter FR, Smirnov PE, Liu T, Avancha R. A one-equation local correlation-based transition model. Flow Turbul Combust 2015;95(4):583–619. http://dx.doi.org/10.1007/s10494-015-9622-4.

[15] Medida S, Baeder J. Application of the correlation-based Gamma-re theta transition model to the spalart-allmaras turbulence model. In: 20th AIAA computational fluid dynamics conference. American Institute of Aeronautics and Astronautics; 2011, p. 3979. http://dx.doi.org/10.2514/6.2011-3979.

[16] Piotrowski MGH, Zingg DW. Smooth local correlation-based transition model for the Spalart-Allmaras Turbulence Model. AIAA J 2021;59(2):474–92. http://dx.doi.org/10.2514/1.J059784.

[17] Coder J, Maughmer M. One-equation transition closure for eddy-viscosity turbulence models in CFD. In: 50th AIAA aerospace sciences meeting including the new horizons forum and aerospace exposition. Aerospace sciences meetings, American Institute of Aeronautics and Astronautics; 2012, p. 672. http://dx.doi.org/10.2514/6.2012-672.

[18] Medida S, Baeder J. A new crossflow transition onset criterion for RANS turbulence models. In: 21st AIAA computational fluid dynamics conference. Fluid dynamics and co-located conferences, American Institute of Aeronautics and Astronautics; 2013, p. 3081. http://dx.doi.org/10.2514/6.2013-3081.

[19] Grabe C, Krumbein A. Extension of they-Re$_{6}$t model for prediction of crossflow transition. In: 52nd aerospace sciences meeting. American Institute of Aeronautics and Astronautics; 2014, p. 1269. http://dx.doi.org/10.2514/6.2014-1269.

[20] Piotrowski MGH, Zingg DW. Compressibility corrections to extend a smooth local correlation-based transition model to transonic flows. Aeronaut J 2023;127(1313):1141–70. http://dx.doi.org/10.1017/aer.2022.105.

[21] Reed WH, Hill TR. Triangular mesh methods for the neutron transport equation. Tech. rep. LA-UR-73-479; CONF-730414-2, Los Alamos Scientific Lab., N.Mex. (USA); 1973.

[22] Cockburn B, Shu C-W. The local discontinuous Galerkin method for time-dependent convection-diffusion systems. SIAM J Numer Anal 1998;35(6):2440–63. http://dx.doi.org/10.1137/S0036142997316712.

[23] Arnold DN, Brezzi F, Cockburn B, Marini LD. Unified analysis of discontinuous Galerkin methods for elliptic problems. SIAM J Numer Anal 2002;39(5):1749–79. http://dx.doi.org/10.1137/S0036142901384162.

[24] Bassi F, Rebay S. A high-order accurate discontinuous finite element method for the numerical solution of the compressible Navier-Stokes Equations. J Comput Phys 1997;131(2):267–79. http://dx.doi.org/10.1006/jcph.1996.5572.

[25] Bassi F, Rebay S. A high order discontinuous Galerkin method for compressible turbulent flows. In: Cockburn B, Karniadakis GE, Shu C-W, editors. Discontinuous Galerkin methods. Lecture notes in computational science and engineering,
Berlin, Heidelberg: Springer; 2000, p. 77–88. http://dx.doi.org/10.1007/978-3-642-59721-3_4.

[26] Renac F. Stationary discrete shock profiles for scalar conservation laws with a discontinuous Galerkin method. SIAM J Numer Anal 2015;53(4):1690–715. http://dx.doi.org/10.1137/14097906X.

[27] Crivellini A, D'Alessandro V. Spalart–Allmaras model apparent transition and RANS simulations of laminar separation bubbles on airfoils. Int J Heat Fluid Flow 2014;47:70–83. http://dx.doi.org/10.1016/j.ijheatfluidflow.2014.03.002.

[28] Lorini M, Bassi F, Colombo A, Ghidoni A. High-order implementation of a non-local transition model in a DG solver for turbomachinery applications. Comput & Fluids 2016;127:115–30. http://dx.doi.org/10.1016/j.compfluid.2015.12.009.

[29] Lorini M, Bassi F, Colombo A, Ghidoni A, Noventa G. Discontinuous Galerkin solution of the RANS and kl-k-log(ω) equations for natural and bypass transition. Comput & Fluids 2021;214:104767. http://dx.doi.org/10.1016/j.compfluid.2020.104767.

[30] Crivellini A, Ghidoni A, Noventa G. Algebraic modifications of the K-Ω and Spalart–Allmaras turbulence models to predict bypass and separation-induced transition. Comput & Fluids 2023;253:105791. http://dx.doi.org/10.1016/j.compfluid.2023.105791.

[31] Cakmakcioglu SC, Bas O, Mura R, Kaynak U. A revised one-equation transitional model for external aerodynamics. In: AIAA aviation 2020 forum. AIAA aviation forum, American Institute of Aeronautics and Astronautics; 2020, http://dx.doi.org/10.2514/6.2020-2706.

[32] Renac F, de la Llave Plata M, Martin E, Chapelier JB, Couaillier V. Aghora: a high-order DG solver for turbulent flow simulations. In: Kroll N, Hirsch C, Bassi F, Johnston C, Hillewaert K, editors. IDIHOM: industrialization of high-order methods - a top-down approach: results of a collaborative research project funded by the European union, 2010 - 2014. Notes on numerical fluid mechanics and multidisciplinary design, Cham: Springer International Publishing; 2015, p. 315–35. http://dx.doi.org/10.1007/978-3-319-12886-3_15.

[33] Ilinca F, Pelletier D. Positivity preservation and adaptive solution for the K-ε model of turbulence. AIAA J 1998;36(1):44–50. http://dx.doi.org/10.2514/2.350.

[34] Fehrs M. Boundary layer transition in external aerodynamics and dynamic aeroelastic stability (Ph.D. thesis), Technische Universität Braunschweig; 2018.

[35] Zauner M, Moise P, Sandham ND. On the co-existence of transonic buffet and separation-bubble modes for the OALT25 laminar-flow wing section. Flow Turbul Combust 2023;110(4):1023–57. http://dx.doi.org/10.1007/s10494-023-00415-4.

[36] Allmaras SR, Johnson FT. Modifications and clarifications for the implementation of the Spalart-Allmaras turbulence model. In: Seventh international conference on computational fluid dynamics (ICCFD7). Vol. 1902, Big Island, HI; 2012.

[37] Burgess N, Mavriplis D. Robust computation of turbulent flows using a discontinuous Galerkin method. In: 50th AIAA aerospace sciences meeting including the new horizons forum and aerospace exposition. American Institute of Aeronautics and Astronautics; 2012, p. 457. http://dx.doi.org/10.2514/6.2012-457.

[38] Bassi F, Crivellini A, Rebay S, Savini M. Discontinuous Galerkin solution of the Reynolds-averaged Navier–Stokes and k-ω turbulence model equations. Comput & Fluids 2005;34(4):507–60. http://dx.doi.org/10.1016/j.compfluid.2003.08.004.

[39] Lee B, Baeder JD. Prediction and validation of laminar-turbulent transition using SA-γ transition model. In: AIAA scitech 2021 forum. American Institute of Aeronautics and Astronautics; 2021, p. 1532. http://dx.doi.org/10.2514/6.2021-1532.

[40] Nichols RH. Addition of a local correlation-based boundary layer transition model to the CREATE\textsuperscript{scripttm-}AV <italic>kestrel</ italic> unstructured flow solver. In: AIAA scitech 2019 forum. American Institute of Aeronautics and Astronautics; 2019, http://dx.doi.org/10.2514/6.2019-1343.

[41] Coder JG, Maughmer MD. Computational fluid dynamics compatible transition modeling using an amplification factor transport equation. AIAA J 2014;52(11):2506–12. http://dx.doi.org/10.2514/1.J052905.

[42] Spalart P. Trends in turbulence treatments. In: Fluids 2000 conference and exhibit. Fluid dynamics and co-located conferences, American Institute of Aeronautics and Astronautics; 2000, p. 2306. http://dx.doi.org/10.2514/6.2000-2306.

[43] Bassi F, Botti L, Colombo A, Di Pietro DA, Tesini P. On the flexibility of agglomeration based physical space discontinuous Galerkin discretizations. J Comput Phys 2012;231(1):45–65. http://dx.doi.org/10.1016/j.jcp.2011.08.018.

[44] Roe PL. Approximate Riemann solvers, parameter vectors, and difference schemes. J Comput Phys 1997;135(2):250–8. http://dx.doi.org/10.1006/jcp.1997.5705.

[45] Harten A, Hyman JM. Self adjusting grid methods for one-dimensional hyperbolic conservation laws. J Comput Phys 1983;50(2):235–69. http://dx.doi.org/10.1016/0021-9991(83)90066-9.

[46] Crivellini A, Bassi F. An implicit matrix-free discontinuous Galerkin solver for viscous and turbulent aerodynamic simulations. Comput & Fluids 2011;50(1):81–93. http://dx.doi.org/10.1016/j.compfluid.2011.06.020.

[47] Schubauer GB, Klebanoff PS. Contributions on the mechanics of boundary-layer transition. 1956.

[48] Coupland J. ERCOFTAC special interest group on laminar to turbulent transition and retransition: T3A and T3B test cases. 1990, A309514.

[49] McGhee RJ, Viken JK, Pfenninger W, Beasley WD, Harvey WD. Experimental results for a flapped natural-laminar-flow airfoil with High Lift/Drag Ratio. Tech. rep. NASA-TM-85788, 1984.

[50] Langel CM, Chow R, Van Dam CCP. A comparison of transition prediction methodologies applied to high Reynolds number external flows. In: 54th AIAA aerospace sciences meeting. AIAA sciTech forum, American Institute of Aeronautics and Astronautics; 2016, p. 0551. http://dx.doi.org/10.2514/6.2016-0551.

[51] Hebler A. Aeroelastik am Laminarflügel - Versuchsbericht ALF-3. 2013.

[52] Ströer P, Krimmelbein N, Krumbein A, Grabe C. Stability-based transition transport modeling for unstructured computational fluid dynamics at transonic conditions. AIAA J 2021;59(9):3585–97. http://dx.doi.org/10.2514/1.J059892.

[53] Fehrs M, van Rooij ACL, Nitzsche J. Influence of boundary layer transition on the flutter behavior of a supercritical airfoil. CEAS Aeronaut. J. 2015;6(2):291–303. http://dx.doi.org/10.1007/s13272-014-0147-7.

[54] Arnal D, Habiballah M, Coustols E. Laminar instability theory and transition criteria in two and three-dimensional flow. Rech Aerospat (Engl Ed) 1984;2:45–63, (ISSN 0379-380X).

[55] Crouch JD, Garbaruk A, Magidov D, Travin A. Origin of transonic buffet on aerofoils. J Fluid Mech 2009;628:357–69. http://dx.doi.org/10.1017/S0022112009006673.

[56] Crouch JD, Garbaruk A, Magidov D. Predicting the onset of flow unsteadiness based on global instability. J Comput Phys 2007;224(2):924–40. http://dx.doi.org/10.1016/j.jcp.2006.10.035.

[57] Sartor F, Mettot C, Sipp D. Stability, receptivity, and sensitivity analyses of buffeting transonic flow over a profile. AIAA J 2015;53(7):1980–93. http://dx.doi.org/10.2514/1.J053588.

[58] Bernardos L, Richez F, Gleize V, Gerolymos GA. Algebraic nonlocal transition modeling of laminar separation bubbles using K-ω turbulence models. AIAA J 2019;57(2):553–65. http://dx.doi.org/10.2514/1.J057734.

[59] Langtry RB. A correlation-based transition model using local variables for unstructured parallelized CFD codes (Ph.D. thesis), University of Stuttgart; 2006, http://dx.doi.org/10.18419/opus-1705.

[60] Guermond J-L, Pasquetti R, Popov B. Entropy viscosity method for nonlinear conservation laws. J Comput Phys 2011;230(11):4248–67. http://dx.doi.org/10.1016/j.jcp.2010.11.043.