12th International Symposium on Turbulence and Shear Flow Phenomena (TSFP12)
Osaka, Japan (Online), July 19-22, 2022

TOWARDS A GENERALIZABLE DATA-DRIVEN APPROACH TO PREDICTSEPARATION-INDUCED TRANSITION

Vishal Srivastava
Department of Aerospace Engineering
University of Michigan
Ann Arbor, MI 48105
vsriv@umich.edu

Karthik Duraisamy
Department of Aerospace Engineering
University of Michigan
Ann Arbor, MI 48105
kdur@umich.edu

ABSTRACTof Expressions (EVE) by Waschkowski et al. (2022) among
others. It should be noted, however, that model augmentations
extracted using several of these frameworks struggle to improve predictive accuracy across a wide range of geometries
and physical conditions.
Srivastava & Duraisamy (2021) presented Learning
and Inference assisted by Feature-space Engineering (LIFE)
framework which is based on the IIML framework by Holland et al. (2019) and provides expert modelers with guiding
principles, tools and techniques to extract generalizable, robust and model-consistent augmentations from limited data.
This work presents a hierarchical model augmentation strategy
via the LIFE framework that involves consecutive introduction
and inference of two augmentation functions, β1 and β2 , into abaseline turbulence model to enable prediction of bypass transition caused by freestream disturbances and flow separation,
respectively. The robustness offered by the LIFE framework(by reverting to baseline model behavior within feature-space
regions where no training data is available) makes it easier
to infer hierarchical augmentations. Designing an appropriate
blending function is key to extend the applicability of hierarchically augmented models to a wide range of geometries andoperating conditions. A blending function, σ , is designed inthis work to isolate the effects of the hierarchical augmentation to cases involving flow separation.

This work introduces a hierarchical data-driven model
augmentation technique to alleviate structural inadequacies ina transition model. This technique is demonstrated by introducing two different model augmentations, β1 and β2 within
Wilcox’s 1988 k-ω turbulence model. These functions are consecutively inferred in a model-consistent fashion as functionsof local flow quantities to improve predictions for bypass transition caused by freestream disturbances and flow separation,
respectively. The key step described in the process is designing a physics-informed blending function to isolate the effectsof the β2 only to cases involving flow separation. β1 is inferred from only two flat plate cases from the ERCOFTAC
dataset, whereas LES data for only one compressor cascadecase is used to infer β2 . While both, β1 and β2 , consistently
improve predictive accuracy across all unseen test cases, it isobserved that transition is predicted significantly upstream forcases involving transition of attached flow under adverse pressure gradients, and some cases involving flow separation.

INTRODUCTION
While Reynolds-Averaged Navier-Stokes (RANS) simulations have been the industry workhorse for applications involving design and optimization of flow paths, turbulence closures used in RANS models consist of model-form inadequacies that can result in significant prediction errors for cases
involving complex geometries or unseen physical conditions.
Over the past decade, the emergence of machine learning asa scientific tool has led to the development of several datadriven frameworks that can extract corrections to reducedfidelity models from high-fidelity data. These include theField Inversion and Machine Learning (FIML) framework by
Parish & Duraisamy (2016), Integrated Inference and Machine
Learning (IIML) by Holland et al. (2019), Sparse Regressionof Turbulence Stress Anisotropy (SpaRTA) by Schmelzer et al.
(2020), Deep learning PDE Model (DPM) by Sirignano et al.
(2020), Data Assimilation and Machine Learning by Volpiani
et al. (2021), End-to-end Differentiable Learning by Ströfer &Xiao (2021), and EVolutionary algorithm for the development

METHODOLOGY
The LIFE framework is explained as follows. Assume
that a baseline model is given as R(u; ξ ) = 0 where u and ξ
represent the model states and the inputs to the model (mesh,boundary conditions, etc.), respectively. An augmentation
function β (η(u; ξ ); w) is then carefully introduced within themodel such that it can modify model behavior for physical conditions affected by the inadequacy in consideration and has
a minimal influence on other physical conditions. The augmented model can then be written as R(u, β (η(u; ξ )); w)) =0. Here η refers to features which are local geometricallyinvariant quantities derived from model states. The additional dependence of η on ξ is to account for features be-1

12th International Symposium on Turbulence and Shear Flow Phenomena (TSFP12)
Osaka, Japan (Online), July 19-22, 2022
as an automatic switch to activate/deactivate the β2 augmentation depending on whether the flow exhibits flow separation.
The baseline value for β1 , β2 and σ is unity which makes
γ = 1 everywhere in the physical domain, hence ensuring baseline k-ω model predictions. It should be noted here that two
augmentation functions are needed because a single augmentation function was unable to establish a functional relationship consistent with bypass transition originating due to both,freestream disturbances as well as flow separation. Given how
the augmentation term is introduced in the model, the spatial
field for intermittency (γ) closely resembles the β1 β2σ field. A
transport equation is used, instead of directly multiplying thesource term in the k-equation with β1 β2σ , because it provides
a smoother field for the solver to work with, hence preserving
stability and also because it brings in the effects of flow history
within the intermittency.

ing functions of quantities like wall distance or gradients ofmodel states. The augmentation function β is characterized
by function parameters w. The LIFE framework offers several
guiding principles for feature design, the major ones among
which include choice of features relevant to the inadequacy inconsideration, physics-based (instead of statistics-based) nondimensionalization, a bounded functional form for features,and ensuring that physical conditions that require different
augmentation values are manifested in different regions within
the feature-space. In addition to these guidelines, LIFE also
introduces localized learning techniques to ensure that the augmentation behavior remains unchanged in the feature-space
regions that do not contain data. This offers robustness andensures that the augmented model falls back to its baseline behavior for unseen physical conditions. To infer the optimal values of these parameters from available data, the inverse problem shown in Eqn. 1 is solved.
N∗

w = arg min ∑ [Cw

Feature Designk(ykdata , ykpred (uk ; ξ k ))k=1

The features used in this work are identical to those proposed by Srivastava & Duraisamy (2021). Both, β1 and β2 arefunctions of the same three features. These features are briefly
described as follows.

(1)+ λ k T k (β (η(uk ; ξ k ); w))]
Here, k is the index assigned to individual training datasets,
C refers to the cost function defining the discrepancy between sparsely available data and corresponding predictions
for quantities y, and T refers to a regularization function.
Once the minimization problem is solved, the optimal weightsw∗ can be used within the augmented solver for predictions as
R(u, β (η(u; ξ )); w∗ )) = 0. Note that this work does not use
any explicit regularization as there is ample implicit regularization as a consequence of inference from multiple cases, theaugmentation function being a function of specified featuresand even the underlying functional form of the augmentation
being constrained in the feature-space.

dw2 Ω
2.188U∞η1 = min

r

9ω∞
,3
7ν

!
(3)

The first feature is the ratio of the vorticity Reynolds numberReΩ = dw2 Ω/(2.188ν) (see Durbin (2012)) and Reθ ,tr , whereθtr is the correlation for transition momentum
thickness givenpby Praisner & Clark (2004) as θtr = 7ν/9ω∞ . All freestream
quantities are extracted from a user-specified distance awayfrom the wall. A conservative upper-bound of 3 is also explicitly applied.

Augmenting Wilcox’s 1988 k-ω modelη2 =A bare-bones intermittency transport equation, inspired
from the model by Durbin (2012), was added to Wilcox’s 1988k-ω model (Wilcox (1988)) and augmented as shown in Eqn.
2.

dw ω
√
dw ω + k(4)νν + νt

(5)

η3 =

D(ρk)
= ∇ · ((µ + σk µt )∇k) + γPk −Cµ ρkω
Dt
ω
D(ρω)
= ∇ · ((µ + σω µt )∇ω) +Cω1 Pk −Cω2 ρω 2 (2)
Dt
k

√
D(ργ)
= ∇ · (σγ,ℓ µ + σγ,t µt )∇k + ρΩ(β1 β2σ − γ) γ
Dt

The second and third features distinguish the laminar boundary layer upstream of the transition location from near-wall regions in a turbulent boundary layer and turbulent flow regions,
respectively. It should be noted here that it is in these laminarboundary layer regions where the augmentation values need tobe significantly lower than their baseline value of unity.

Here, the production of turbulent kinetic energy Pk is given as
follows.

Localized Learning
An interpolation-based C0 -continuous functional form
was chosen to perform localized learning for this application.
The feature-spaces corresponding to both, β1 and β2 were discretized into uniform grids. The augmentation values on thenodes of these grids act as function parameters, while the augmentation at any point in the feature-space is evaluated via
multi-linear interpolation based on these nodal values (see Fig.
1). Note here that since the functional form is constrained to
be C0 -continuous, sudden jumps in the augmentation values
are not possible within the feature-space. On the other hand,
the C0 -continuity provides for excellent convergence characteristics of the flow solver.

2 ∂ uiPk = µt Ω − ρk
3 ∂ xi2

The main differences compared to Durbin’s model are that theγmax term is replaced by β1 β2σ , and that the function Fγ isabsent from the source term. β1 is inferred using data fromonly two flat plate cases to enable predictions for bypass transition due to freestream disturbances. β2 is inferred from onlya single compressor cascade case (after β1 is inferred) to enable predictions of separation-induced bypass transition. The
blending function σ can vary between zero and unity and acts2

12th International Symposium on Turbulence and Shear Flow Phenomena (TSFP12)
Osaka, Japan (Online), July 19-22, 2022

Figure 2.

Figure 3.

Figure 1. Schematic detailing C0 -continuous interpolationbased augmentation functions

Designing the Blending Function
To create a blending function σ that varies between zeroand unity, one or more sigmoid functions can be applied to an
appropriate quantity√ fσ . Ge et al. (2014), used the function
ωdw (nw · ∇)|S|/ 2|S|2 to distinguish regions that exhibit flow
separation. Here nw is the wall-normal direction corresponding to the nearest point on the wall. This can be calculated byevaluating the gradient of the wall distance (dw ) and normalizing the result (to ensure that the magnitude is unity). Notethat for body-fitted grids with low cell skewness, this approximation is fairly accurate close to the walls and that is exactly
where the blending function is needed. |S| denotes the magnitude of the strain rate tensor. Three changes were made tothis function to formulate an appropriate fσ that was found toeffectively detect separated flow. The changes include the useof vorticity
p magnitude Ω instead of |S|, use of laminar lengthscale ν/Ω instead of dw , and using a bounded functional
form. The final functional form of fσ is shown in Eqn. 6.

fσ =

√


ν(nw · ∇)Ωω√√ν|(nw · ∇)Ω| + Ω1.5
2Ω + ω

1
1 + exp(−( fσ + 0.05)/0.003)

Mesh used for T3C1, T3C2, T3C3 and T3C5 cases

Table 1.

Inflow conditions for the T3 test cases

Cases

Tuin

νtin

L (in m)

ReL,in

T3A

0.035

14.0

1.5

520000

T3C1

0.1

50.0

1.65

660000

T3B

0.065

100.0

1.5

940000

T3C2

0.037

12.0

1.65

550000

T3C3

0.034

8.0

1.65

418000

T3C5

0.043

17.0

1.65

946000Inferring β1
For this step, β2 was constrained to unity for all regionsin the feature space, which means all function parameters forβ2 were set to unity. Flat plate cases T3A and T3C1 fromthe ERCOFTAC dataset are used to train β1 whereas other flat
plate cases viz., T3B, T3C2, T3C3 and T3C5, were used to test
the inferred augmentation. Meshes used for these cases can beseen in Figs. 2 and 3. The inflow conditions for all these cases
are listed in Table 1.

(6)

This function mainly targets the region within a separated flowwhere Ω increases with wall distance i.e., where fσ is positive. However, it was seen that the σ needs to remain activefor very small negative values of fσ as well, to predict transition accurately for cases exhibiting separated boundary layers.
This is probably because fσ is positive in a very thin region
and intermittency needs to be lowered within a slightly thicker
region for k-production to sufficiently drop in order to predict
transition at the correct location. The final functional form of
σ is given in Eqn. 7.

σ=Mesh used for T3A and T3B cases

Training Only cases T3A and T3C1 were used to infer function parameters β1 . Note here that transition occurs
in a zero pressure gradient region for T3A and a favorablepressure gradient for T3C1. The training procedure is similar to that used by Srivastava & Duraisamy (2021). The
cost functions for both these cases were chosen to be the sum
squared discrepancy in the skin friction coefficients i.e., C =∥C f ,pred − C f ,data ∥22 . As stated previously, no regularization
was used as ample implicit regularization is already present.
The functional form corresponding to the interpolation-based
C0 -continuous augmentation function was characterized bydiscretizing the feature-space into a grid with 45 × 15 × 15
cells. Steepest gradient descent was used to solve the inference problem and the method of adjoints (along with automatic differentiation) was used to evaluate sensitivities. A step
size of 0.1 was used. The intermittency was restricted within
its physical bounds of 0 and 1 after every steepest gradient descent step. The corresponding predictions for both these cases(7)

RESULTSThis section describes the results obtained after inferring
β1 (β2 assumed to be unity everywhere in the feature-space),
inferring β2 (β1 frozen to its inferred form and σ = 1) andintroducing the functional form for σ in a step-wise manner.

3

12th International Symposium on Turbulence and Shear Flow Phenomena (TSFP12)
Osaka, Japan (Online), July 19-22, 2022

Figure 7.
(a) T3A

Blade geometries for compressor cascade(b) T3C1

Figure 4. Skin friction predictions for T3C1. Blue dots Data; Red curves - Predictions using baseline model; Green
curves - Predictions with parameters inferred for β1 only

Figure 5.

Residuals vs solver iterations

are shown in Fig. 4. The transition location is predicted wellfor both these cases while achieving better laminarization upstream to the location compared to that observed in Srivastava& Duraisamy (2021). The residual convergence plot when the
T3A case is solved with the augmented model is shown in Fig.
5. The convergence characteristics are far superior compared
to what was seen in Srivastava & Duraisamy (2021). This isbecause of the continuous nature of the functional form.

Figure 8.
casecurate, transition is predicted significantly upstream of the true
locations for T3C2 and T3C3. Even though the model was
not trained on cases where transition occurs in adverse pressure gradients (which is the case for T3C2 and T3C3), thesepredictions are worse compared to those in Srivastava & Duraisamy (2021). This is because a continuous augmentation
does not allow for sudden jumps in augmentation values within
the feature space and hence intermittency has a more smeared
growth. As soon as this smeared growth results in intermittency values high enough to facilitate k-production, transition
is predicted. Hence, while C0 -continuous augmentation functions provide more laminarization and significantly better convergence characteristics, they also lead to a lower predictive
accuracy for transition in adverse pressure gradient regions.
This can be corrected in the future by using a more adaptivefunctional form for the augmentation function.

Testing The prediction results of this inferred augmentation on other flat plate cases are shown in Fig. 6. While(a) T3B

Mesh used for NACA 65-010 compressor cascade(b) T3C2

Inferring β2(c) T3C3

LES data for six compressor cascade cases was provided
by RTRC. The geometries along with their blade stagger angles have been shown in Fig. 7. The mesh used for theNACA65-010 compressor cascade case is shown in Fig. 8.
Similar meshes were used for other geometries. The inflow
and outflow conditions for the compressor cascade cases havebeen shown in Tables 2 and 3. Fig. 9 shows the transition
predictions made using the augmented model on a NACA 65010 compressor cascade geometry. In this case, transition waspredicted at the separation location which is significantly upstream of the actual transition location. Holding function parameters of β1 constant at their inferred values, and assuming
σ a constant baseline value of 1 for this step, function parameters for β2 are inferred as follows.

(d) T3C5

Figure 6. Skin friction predictions for T3 testing cases. Blue
dots - Data; Red curves - Predictions using baseline model;
Green curves - Predictions with parameters inferred for β1
only

the transition predictions for T3B and T3C5 are reasonably ac-412th International Symposium on Turbulence and Shear Flow Phenomena (TSFP12)
Osaka, Japan (Online), July 19-22, 2022
Table 2. Inflow conditions common for all compressor cascade cases

p0,in

T0,in

Tuin

νt,in /ν∞14.7705724 PSI

288.5672892 K

1%

10Figure 10. Wall shear stress predictions for NACA 65-010.
Black dots - Data; Red curves - Predictions with parametersinferred for β1 only; Green curves - Predictions with parameters inferred for both, β1 and β2

Table 3. Outlet back-pressure and flow angles for compressor cascade cases

Cases

pb (in PSI)

Flow angles

NACA65-010

14.695946

45◦

NACA65-410

14.715946

45◦

NACA65-1210

14.725946

45◦

NACA65-1810

14.725946

45◦

NACA65-1510

14.735946

60◦

NACA65-2110

14.735946

60◦(a) NACA 65-410(b) NACA 65-1210(c) NACA 65-1510(d) NACA 65-1810

Figure 9. Wall shear stress predictions for NACA 65-010.
Black dots - Data; Red curves - Predictions using baselinemodel; Green curves - Predictions with parameters inferredfor β1 only(e) NACA 65-2110

Figure 11. Wall shear stress predictions for unseen compressor cascade cases. Black dots - Data; Red curves - Predictions
with parameters inferred for β1 only; Green curves - Predictions with parameters inferred for both, β1 and β2

Training The training process for β2 was identical tothat used for β1 with three differences.
• Only the NACA 65-010 compressor cascade case wasused to infer β2 .
• The feature-space for β2 was divided into 30 × 10 × 10
cells.
• The cost function C was set to the sum squared discrepancy between the predictions and LES data for wall shear
stress.

using only β1 . The corresponding plots for NACA 65-410,
NACA 65-1210, NACA 65-1510, NACA 65-1810 and NACA
65-2110 are shown in Fig. 11. In some cases transition is predicted considerably upstream of the true transition location.
However, the model shows consistent and significant improvements in predictive accuracy even for these cases when compared to using just β1 .

Fig. 10 shows the wall shear stress profile predicted using thehierarchical augmentation. The predicted transition locationshows excellent agreement with data. However, regions downstream of the transition location show a lower wall shear stress.
This discrepancy is due to other inadequacies in the turbulencemodel and cannot be alleviated within the bounds set for thecurrent augmentation.

Introducing σ
In the absence of any blending function to isolate effectsof β2 , the predictive accuracy on the flat plate cases is abysmal
when both β1 and β2 are used as shown by red curves in Fig.
12. However, when the blending function is used, the transition location predictions are restored (green curves) close towhere they are predicted when using only β1 (black dashed

Testing When tested on other compressor cascadeconfigurations, the hierarchically augmented model showeda significant improvement compared to the results obtained
5

12th International Symposium on Turbulence and Shear Flow Phenomena (TSFP12)
Osaka, Japan (Online), July 19-22, 2022
cies. In addition to hierarchical augmentations, this work also
provides a brief commentary on continuous functional forms
for augmentations. As is demonstrated in the results, while
such functional forms provide excellent solver convergence,
they could slightly hamper accuracy for applications like bypass transition where sharp jumps in augmentation values in
the feature-space might provide better results.
(a) T3A(b) T3C1ACKNOWLEDGEMENTS
The authors acknowledge funding for this work from the
ARPA-E Differentiate program (Award no. DE-AR0001201).
The authors also thank RTRC for providing LES data for compressor cascade geometries.

(c) T3B(d) T3C2REFERENCES(e) T3C3

Durbin, Paul 2012 An intermittency model for bypass transition. International Journal of Heat and Fluid Flow 36, 1–6.
Ge, Xuan, Arolla, Sunil & Durbin, Paul 2014 A bypass transition model based on the intermittency function. Flow, turbulence and combustion 93 (1), 37–61.
Holland, Jonathan R, Baeder, James D & Duraisamy, Karthik
2019 Towards integrated field inversion and machine learning with embedded neural networks for rans modeling. In
AIAA Scitech 2019 Forum, p. 1884.
Parish, Eric J & Duraisamy, Karthik 2016 A paradigm fordata-driven predictive modeling using field inversion andmachine learning. Journal of Computational Physics 305,
758–774.
Praisner, T. J. & Clark, J. P. 2004 Predicting Transition in Turbomachinery—Part I: A Review and New Model Development. Journal of Turbomachinery 129 (1), 1–13.
Schmelzer, Martin, Dwight, Richard P & Cinnella, Paola
2020 Discovery of algebraic reynolds-stress models using
sparse symbolic regression. Flow, Turbulence and Combustion 104 (2), 579–603.
Sirignano, Justin, MacArt, Jonathan F & Freund, Jonathan B
2020 Dpm: A deep learning pde augmentation method withapplication to large-eddy simulation. Journal of Computational Physics 423, 109811.
Srivastava, Vishal & Duraisamy, Karthik 2021 Generalizable
physics-constrained modeling using learning and inference
assisted by feature-space engineering. Phys. Rev. Fluids 6,124602.
Ströfer, Carlos A Michelén & Xiao, Heng 2021 End-to-end
differentiable learning of turbulence models from indirect
observations. Theoretical and Applied Mechanics Letters
11 (4), 100280.
Volpiani, Pedro Stefanin, Meyer, Morten, Franceschini, Lucas,
Dandois, Julien, Renac, Florent, Martin, Emeric, Marquet,
Olivier & Sipp, Denis 2021 Machine learning-augmented
turbulence modeling for rans simulations of massively separated flows. Physical Review Fluids 6 (6), 064607.
Waschkowski, Fabian, Zhao, Yaomin, Sandberg, Richard &Klewicki, Joseph 2022 Multi-objective cfd-driven development of coupled turbulence closure models. Journal ofComputational Physics 452, 110922.
Wilcox, David C 1988 Reassessment of the scale-determining
equation for advanced turbulence models. AIAA journal
26 (11), 1299–1310.

(f) T3C5

Figure 12. Skin friction predictions for T3 cases. Blue dots
- Data; Black dashed curves - Predictions with parameters inferred for β1 only; Red curves - Predictions with parametersinferred for both, β1 and β2 , with σ = 1; Green curves - Predictions with parameters inferred for both, β1 and β2 , withσ = 1/(1 + exp(−( fσ + 0.05)/0.003))

curves). The predictions for transition under adverse pressure gradients (cases T3C2 and T3C3) actually improves whenβ1 β2σ is used instead of just β1 . It should be noted here thatthe blending function has negligible effect on the compressor
cascade cases.

CONCLUSIONSThis work describes a hierarchical augmentation technique to infer multiple augmentation functions that can address
various model-form inadequacies manifested in different physical regimes. The process is demonstrated by consecutively inferring two augmentation functions within a bare-bones intermittency function that was added to Wilcox’s 1988 k-ω model
to predict transition caused by freestream disturbances andflow separation, respectively. An appropriate blending function is designed to isolate the effects of the second augmentation to cases involving separation-induced transition. The augmentation consistently improves predictions across a range of
geometries and flow conditions. However, significant inaccuracies are observed for cases involving transition of attachedflows in adverse pressure gradients and some cases involving
separation-induced transition as well. In such cases, the augmentation predicts an early transition. Future work involves
development of better features and a formal framework to optimize grid resolution in feature-space to reduce these inaccura-6