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

INVESTIGATION OF REYNOLDS STRESSES PRIOR TO VORTEXBREAKDOWN ON A TRIPLE-DELTA WING AT TRANSONIC CONDITION

Tony Di Fabbio, Eike Tangermann and Markus Klein
Department of Aerospace Engineering
University of the Bundeswehr Munich
85577 Neubiberg, Germany
tony.difabbio@unibw.de

ABSTRACT

been analysed in the current work. The simulations have been
performed employing the DLR TAU-Code (Gerhold (2005))
at transonic conditions with Ma = 0.85, Re = 12.53 · 106 for
α = 20° and β = 5°. The results from three different simulations featuring different turbulence model approaches are presented. Results from URANS based on the linear Boussinesq
assumption (LBA), the QCR.2020 nonlinear model version of
Spalart-Allmaras (QCR) and hybrid RANS/LES (IDDES) are
discussed. The SAneg-based IDDES model is applied in thescale-resolving computations, whereas the Spalart-Allmaras
One-Equation Model with corrections for negative turbulent
viscosity and rotation/curvature correction (SA-negRC) is employed to close the URANS equations within the QCR and
LBA approaches.
Since the flow separation, which forms the initial stageof vortex formation, is fixed by the sharp leading edge, themain challenge within the simulation is to correctly produce
formation and further development of the vortical flow system along the wing surface, which is primarily affected by
the treatment of turbulence in terms of modeling or resolving turbulent eddies. The IDDES serves as a reference, providing the highest affordable fidelity of flow field prediction.
The LBA URANS model, even without the necessity to predict
the flow separation, fails to produce several features (including
vortex breakdown) of the flow field correctly, as it has been already shown in previous work by the authors (Di Fabbio et al.,2022a,b). Finally, the QCR.2020 nonlinear model version of
Spalart-Allmaras is able to capture the vortex breakdown over
the windward wing resulting in a clear improvement in the predicted aerodynamic coefficients.
The present work provides insight to the understandingof the relationship between vortex breakdown and Reynolds
stresses. When used in conjunction with the SA turbulencemodel with rotation/curvature correction, the QCR equation
improves the Reynolds stress predictions as will be shown bycomparison with IDDES results. This improvement is significant in the observed case, though it is noteworthy, that the momentum equation incorporates the divergence of the Reynolds
stresses. Therefore, the gradients of the Reynolds stresses impact the stress-induced vortex breakdown and capturing theirabsolute levels typically is less significant as it has been discussed by Rumsey et al. (2020).

The behavior of leading-edge vortices on a triple-delta
wing configuration has been analysed numerically at a sideslip condition. The flow is transonic, therefore several shockwaves appear over the fuselage with multiple vortex-shock interactions occurring. The paper focuses on the investigationof the Reynolds stress tensor with a detailed analysis of individual components and their relevance for the prediction of
the flow features. Three different approaches of turbulencetreatment and of constitutive relation have been compared inorder to understand the correlation between Reynolds stressfield further upstream and the consecutive appearance of vortex breakdown on the windward wing. The Quadratic Constitutive Relation appears in the end a promising approach inorder to mitigate the deficiencies of the linear Boussinesq assumption and to predict the flow physics over such configuration with a better accuracy.

INTRODUCTION
Solving the Reynolds averaged Navier-Stokes equations
employing one-equation eddy viscosity turbulence models remains the most widely used approach for the prediction of
turbulent flows. The eddy-viscosity turbulence models, commonly used to close the RANS equations, are linear models.
These models derived from the Boussinesq hypothesis appear
to lack capability for prediction of vortical flows. A more realistic effective-viscosity formulation has been proposed byPope (1975). Its advantage over the isotropic-viscosity hypotheses is that the whole velocity-gradient tensor affects thepredicted Reynolds stresses. Two notable consequences of thisare that the Reynolds-stress tensor is realistically modelled andthe influence of streamline curvature on the Reynolds stresses
is incorporated. The failure of isotropic viscosity hypotheses
to provide correct predictions in many flow situations either isdue to inapplicability of an effective-viscosity approach or toinadequacy of the isotropic-viscosity hypotheses. To remedy
some of the shortcomings of the linear eddy-viscosity models,the Quadratic Constitutive Relation (QCR) for eddy viscosity
has been proposed by Spalart (2000). QCR uses a nonlinear
turbulent stress/strain equation, in contrast to the most commonly used linear Boussinesq relation. Therefore, a key driveris the QCR’s ability to modify the incorrect behavior of the
Reynolds normal stresses trying to take its anisotropy property
into account.
The vortex dominated flow around a triple-delta wing
ADS-NA2-W1 as described by Hövelmann et al. (2020) has

TURBULENCE MODELThe SAneg-IDDES model according to Spalart et al.
(1997) and Shur et al. (2008) is based on the following trans-

1

12th International Symposium on Turbulence and Shear Flow Phenomena (TSFP12)
Osaka, Japan (Online), July 19-22, 2022
port equation for the modified eddy viscosity∂(ρ ν̃) +⃗u · ⃗∇(ρ ν̃) = ⃗∇ ·∂tµ + ρ ν̃ ⃗
∇ν̃ +
σ

c
+ ρ b2 (⃗∇ν̃)2 + cb1 ρ S̃ν̃ − cw1 fw ρ
σ

 2
ν̃
.
d˜

(1)

In the original SA-model the length scale d = LRANS = dw (distance to the nearest wall), whereas in the IDDES model d is
replaced with d˜ = LIDDES , defined as follows(a) Slice plane η = 0

d˜ = f˜d (1 + fe )d + (1 + f˜d )CDES Ψ∆
featuring several fitted functions and coefficients of great complexity which shall not be presented further in this context.
Besides, Shur & al. (2000) proposed a streamline curvaturecorrection (SA-RC), which alters the production term with arotation function, applied in the current work just within the
URANS computations as in IDDES the vortices are located in
the LES zone.
Finally, the QCR can be applied as an extension to every
linear eddy-viscosity model τiLj in which the Boussinesq assumption approximates the Reynolds stress tensor as follows

1 ∂ uk2
τiLj = 2µt Si j −δi j − ρkδi j .
3 ∂ xk3(b) Slice plane ξ = 0.3
Figure 1: ADS-NA2-W1 geometry and mesh.

(2)
sections. The computational mesh is of unstructured topology
with 35 prism layers close to the aircraft surface and tetrahedral cells around. The domain size is 50 times the characteristic length L, which is the fuselage length as illustrated inFig. 1a. It consists of about 40 million cells. Figure 1b shows
the inhomogeneous mesh refinement, in which the cells size
varies through the computational domain. The mesh is refined
most close to the leading edge, where the generation of inboard and the outboard vortices starts. The cell refinement
roughly follows the vortices to capture the strong gradients aswell as for the IDDES case the resolved turbulent fluctuations.
In previous work (Di Fabbio et al., 2022b) the mesh has beenvalidated by a convergence study.

The objective is to correct the linear Boussinesq relation byadding quadratic, cubic, or even higher-order terms to obtain anonlinear relation. In the current work, the QCR.2020 method
proposed by Rumsey et al. (2020), which is given in Eq. 3, isused in conjunction with the SA-negRC model.
τi j,QCR2020 = τiLj +
h
i
p
“
“
−Ccr1
Oi j τ Ljk + O jk τikL −Ccr2µt 2WmnWmn δi j

(3)

NUMERICAL APPROACH
Unsteady simulations have been performed using the
DLR TAU-Code (Gerhold (2005)). An implicit dual-time
stepping approach has been selected employing a BackwardEuler/LUSGS implicit smoother. The computation of thefluxes has been performed with a central scheme and the matrix dissipation model has been selected. However, in hybrid
RANS/LES the artificial dissipation should be reduced in orderto prevent excessive damping of the resolved turbulent structures. A (hybrid) low-dissipation low-dispersion discretization
scheme (HLD2) and a vorticity-sensitive sub-grid filter scale,
which enhances the development of turbulent structures on
anisotropic meshes, have been used (Probst & Reuß (2016)).
More details regarding the numerical approach and the time
length series of the simulation to collect the statistics variables
for the different approaches can be found in previous work
(Di Fabbio et al., 2022b).

RESULTSAs a primary indicator of result accuracy, the curves ofrolling and pitching moment coefficients are shown in Fig. 2.
The experimental data according to Hövelmann et al. (2020)are plotted in comparison with the CFD results: URANS
based on linear Boussinesq assumption (LBA); URANS based
on QCR.2020 nonlinear model version of Spalart-Allmaras(QCR); hybrid RANS/LES (IDDES). For a wider overview ofthe predictions, Fig. 2 also includes aerodynamic coefficientsfor other angles of attack.
The QCR.2020 method considerably improves the prediction of the aerodynamic coefficients in comparison to the linear
SA results especially at α = 20°, on which the present work is
focused, therefore. The main reason behind the improvement
is the prediction of a vortex breakdown above the windward
wing which affects the suction footprint over the aircraft as illustrated in Fig. 3. It generates a sudden increase of pressurein the rear part and consequently a drop of local lift. The integral moments of pitching and rolling are very sensitive to suchvariations of the flow pattern as indicated by the gradients in
Fig. 2 which appear at changes of the vortex topology.

GEOMETRY AND MESH
The ADS-NA2-W1 geometry illustrated in Fig. 1 is a
1:30-scaled version of a generic combat aircraft. It is characterized by a triple-delta wing with three different leading-edge

2

12th International Symposium on Turbulence and Shear Flow Phenomena (TSFP12)
Osaka, Japan (Online), July 19-22, 2022(a) Rolling moment(a) IDDES(b) Pitching moment
Figure 2: Aerodynamic coefficients over angle of attack

with Ma = 0.85, Re = 12.53 · 106 and β = 5°.

Fig. 3a shows the pressure root-mean-square (rms) on a
Q-criterion iso-surface and on the aircraft surface from IDDES. prms has been used to identify the magnitude of thefluctuations. They clearly highlight regions of vortex burst
and separated shear-layer. Two well-distinguished vortices arepresent on the leeward wing and two vortices, closer to eachother, are captured on the windward wing breaking down in
the rear part. The vortex breakdown appears as an abruptchange in the flow topology where the flow decelerates and diverges. As it has been discussed by Escudier (1988), locationand mode of breakdown depend on various parameters such as
adverse pressure gradients, type of delta wing planform, angleof attack, sweep angle and interaction with shock waves. As
it can be seen in Fig. 3b, the QCR extension also is able topredict the vortex breakdown. The onset location is not fixed,
it shows a buffeting behavior. This buffeting is visualized byincreased values of prms . Fig. 3 further shows the iso-surfacesof Ma = 1 over the fuselage and highlights the presence of fiveshock waves (denoted by the numbers in Fig. 3b) between thewing apex and the trailing edge.
The exact shock locations over the fuselage can be seenin Fig. 4, where the ratio p/p∞ in the symmetry plane η = 0
is plotted along the surface of the aircraft from the different
simulations. The LBA results do not predict the fourth shock
wave and consequently the vortex breakdown does not appear.
Therefore, the aerodynamic coefficients are predicted wrong,
as it has been already discussed before. The prms plot of the
LBA simulation has not been included in the current work as
no visible fluctuations have been captured over the wing. The
shock waves highlighted in Fig. 4 are caused by the aircraft
shape. The fuselage diameter is not constant, as can be seenin Fig. 1 or more precisely in Hövelmann et al. (2020). The
transonic flow coming from the windward side aligns with thefuselage, follows its shape while being accelerated up to sonic
speed and consequently forms a shock wave. The first threeshocks in Fig. 3b are caused by this phenomenon. The fifth
and final shock is caused by the need for the flow to return to(b) QCR
Figure 3: Rms pressure normalized by free-stream pres-sure prms /p∞ and Q-criterion iso-surface, IDDES andQCR results for α = 20° and β = 5°. Black: isosurfaceof Ma = 1 and vortex center lines.

its initial conditions and adjust its pressure field. These shocks
due to geometric and boundary conditions are well predictedby all the numerical results. However, it is not possible to givea complete explanation for the presence of the fourth shock. It
may be caused by the increase in the sweep angle at the third
section of the leading edge. In the following, shock-vortex
interaction and Reynolds stresses are analyzed in order to understand the reasons behind the occurrence and prediction of
the vortex breakdown on the windward wing.
Fig. 3 also shows indications of the vortex center lines(in black) that have been extracted using an automated criterion. The method considers the maximum of swirl (S =(⃗ω ·⃗u)/(ρ|⃗u|2 )) in a delimited vortex region in order to determine the vortex core. Fig. 5 shows the extracted vortex center
lines from the three different simulation approaches. The simulations predict a similar path. Besides, the vortex breakdown
in the IDDES and QCR results can be identified. Fig. 3 and 5
show that in this case of a triple-delta wing configuration withtwo primary vortices in each side, the inboard vortex is the first
one to break down. Since the present study primarily focuseson investigating the Reynolds stresses prior to vortex breakdown, the following analysis considers only the flow physics
over the windward wing.
Fig. 6 shows the mean Mach number in the vortex core3

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

Figure 4: Shock wave locations on the aircraft fuselage.

Black line indicates aircraft geometry.

(b) Outboard windward vortex
Figure 6: Numerical evolution of the longitudinal mean

Mach number Ma in the leading-edge vortex core.
Dashed lines indicate instantaneous, solid lines averaged
data.

less visible because the flow field is dominated by the inboard
vortex burst, as it is evinced by the turbulent fluctuations of theinstantaneous Mach number plot. However, a drop of the Ma
can be seen at around ξ = 0.66 from the QCR results and theinstantaneous IDDES plot. It indicates the interaction between
the inboard vortex breakdown onset as well as shock number
4 and the outboard vortex.
Since the QCR approach significantly improves the results, the anisotrpy of the Reynolds stress is suspected be a
driving cause behind the LBA mispredictions. Therefore, theReynolds stresses are investigated along the vortex core lines
in proximity to the vortex breakdown in order to better understand their correlation with this crucial phenomenon. Moreover, the comparison between the different model approaches
will help to comprehend the negative effects of the linear
Boussinesq assumption. The specific Reynolds stresses havebeen computed taking into account the resolved and the modeled parts as follows.

Figure 5: Vortex core line.

lines where the different locations of shock-vortex interaction
can be identified. Regarding the IDDES results, the second
and third shock in Fig. 4 interact with the inboard windward
vortex at around ξ = 0.4 and ξ = 0.5, respectively, causing a
reduction in the streamwise velocity and increasing the vulnerability of the vortex (Delery, 1994). As it can be seen in Fig 6a,
the inboard windward vortex breaks at around ξ = 0.62 where
the interaction with the fourth shock occurs. Without the prediction of shock and vortex breakdown, the LBA produces arather different Mach number profile. The QCR results arecloser to the IDDES high fidelity data and improve the overall
prediction of the flow physics. However, the first two shockvortex interactions are captured slightly upstream at ξ = 0.4
and ξ = 0.45 and the breakdown at ξ = 0.65, which is further
downstream over the wing with respect to the IDDES as wellas to the experimental results.
Besides, Fig. 6 illustrates the instantaneous Mach number
with the dashed lines. The fourth shock wave seems to oscillate in a buffeting way. This oscillation causes a displacementof the interaction plane between the vortex and the shock. Thisis the primary cause for the unsteadiness of the vortex breakdown which appears in Fig. 3. A clear onset point for the vortex breakdown cannot be determined. The buffeting will haveto be discussed in further work dedicated to the unsteady behaviour of such flow physics phenomena.
Fig. 6b shows the mean Mach number in the outboard
windward vortex core line. The vortex breakdown location is

mod
res
τitot
j = τi j + τi j

(4)where the two contributions have been calculated using the expressions
′

′

τires
j = −ui u jτimod= τiLj /ρ
j

or

τimod
= τi j,QCR2020 /ρ
j

(5)

(6)

according to the respective approach. Considering the mean
flow field, the modeled Reynolds stresses τimod
are computedj4

12th International Symposium on Turbulence and Shear Flow Phenomena (TSFP12)
Osaka, Japan (Online), July 19-22, 2022
velocity and local mean density. τ11 shows the turbulent behavior of the transported turbulent shear-layer. The turbulent
motion becomes more intense once the vortex breaks down
and the turbulence kinetic energy is then transported downstream. The covariance τ12 can be mainly considered to visualize the vortex cores and the region of separated shear-layer
close to the leading edge, where its value is negative. τ22 indicates that the fluctuations are generated as the flow departs
from the leading-edge. It is also the primary origin of thehigh turbulence kinetic energy in the vortex core. Finally, as it
can be seen from the legend scale, τ23 is the strongest covariance component and it mainly acts in the shear-layer where thecomplex process of separation and rolling-up appears.
It is important to highlight how the QCR extension improve the URANS results, even approaching in magnitude theresolved Reynolds stresses from IDDES. The anisotropy of the
Reynolds stresses has not been predicted by the LBA which
underlines the shortcomings of linear eddy viscosity models
applied for vortical flows. In particular, as it can be seen in
the τ12 plot, QCR allows to model the Reynolds stresses betterclose to the leading edge which are then involved in the rollup process feeding the vortex with turbulence kinetic energy.
This could explain the higher Mach number in the vortex core
seen above and the occurrence of the forth shock which is notpredicted by the LBA. However, both the URANS approaches
overestimate the prediction of τ23 . Consequently, the modeled
vortex is stronger, less vulnerable and the probability of vortex
breakdown occurrence decreases.

(a) Inboard Windward Vortex(b) Outboard Windward Vortex
Figure 7: Numerical evolution of the normalized specificReynolds stress τ11 /U∞ in the leading-edge vortex core.

with the Boussinesq assumption in Eq. 2 or with the QCR.2020
in Eq. 3. On the other hand, the resolved fluctuations have beencollected to compute the resolved Reynolds stresses τiresj .
Fig. 7 shows the normalized specific Reynolds stressτ11 /U∞ along the vortex core lines on the windward side. Velocity fluctuations appear in the leading-edge vortex core even
upstream of the breakdown location and it has been demonstrated by Menke & Gursul (1997) that these fluctuations arecaused by a random displacement of the vortex core. The
averaged Navier Stokes equations only incorporate the divergence of the Reynolds stress tensor. The turbulence effects inthe averaged momentum equation can be represented by the
body force f = ∇ · τi j . It means that a negative value of the
Reynolds stress tensor gradient produces a negative force that
decelerates the flow field over the aircraft. As it can be seen in
Fig. 7, the breakdown is characterized by a strong decrease ofτ11 (the derivative is negative) representing a rapid deceleration of the turbulent flow caused by the aforementioned shock
wave. Further, taking a look at the slope of τ11 before theinboard vortex breakdown (0.5 < ξ < 0.6), the vortex is becoming increasingly vulnerable. The inboard windward vortex
is already losing energy and consequently decreasing its axial
velocity before the breakdown. This process is not capturedby the LBA simulation, where τ11 remains almost constant.
This phenomenon generates an unsustainable flow condition
that leads to the fourth shock over the fuselage and consequently the breakdown. τ11 in the outboard vortex core line
experiences the effect of the inboard vortex breakdown and themagnitude of the captured fluctuations prior to its breakdown
is higher. For this reason, the slope is less pronounced and thebreakdown less visible.
Fig. 8 shows components of the normalized Reynolds
stress tensor τi j = u′i u′j /U∞2 on the windward wing. τi j represents the intensity of the turbulent fluctuations along thethree directions and their covariance. The components of the
Reynolds stress tensor have been normalized by free-stream

CONCLUSIONS AND OUTLOOKS
Numerical simulations have been performed to analyze
the leading-edge vortex behavior on a triple-delta wing configuration in a side slip condition. Results from URANS basedon linear Boussinesq assumption, QCR.2020 nonlinear modelversion of Spalart-Allmaras and hybrid RANS/LES have been
compared.
The simulation results show that the application of the
Quadratic Constitutive Relation improves the accuracy compared to standard linear eddy-viscosity models not only interms of integral coefficients of forces and moments but also
regarding specific flow features like shocks and vortex breakdown. This observation has been supplemented with a discussion of the Reynolds stress components modeled and resolved by the different methods, respectively. The results underline the shortcomings of linear eddy viscosity models dueto the lack of considering the aspects of anisotropy, namely
the differences between streamwise, wall-normal, and lateral
Reynolds stresses. This affects not only the accuracy of theresults but also the resolved flow physics over the wing. The
present study serves as a basis and motivation for the development of a modified constitutive relation for specific applications in future work.

REFERENCES
Delery, Jean M. 1994 Aspects of vortex breakdown. Progress
in Aerospace Sciences 30 (1), 1–59.
Di Fabbio, T., Tangermann, E. & Klein, M. 2022a Flow pattern analysis on a delta wing at transonic speed. AIAASCITECH 2022 Forum, San Diego, CA.
Di Fabbio, T., Tangermann, E. & Klein, M. 2022b Investigation of transonic aerodynamics on a triple-delta wing in sideslip conditions. CEAS Aeronautical Journal .
Escudier, Marcel 1988 Vortex breakdown: Observations

5

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

Figure 8: Normalized specific Reynolds stress components, IDDES (top), QCR (middle) and LBA (bottom). Vortex corelines illustrated in black over the aircraft.

and explanations. Progress in Aerospace Sciences 25 (2),
189–229.
Gerhold, Thomas 2005 Overview of the hybrid rans code tau.
In MEGAFLOW - Numerical Flow Simulation for Aircraft
Design (ed. Norbert Kroll & Jens K. Fassbender), p. 81–92.
Berlin, Heidelberg: Springer Berlin Heidelberg.
Hövelmann, A, Winkler, A, Hitzel, S M, Richter, K & Werner,
M 2020 Analysis of vortex flow phenomena on generic delta
wing planforms at transonic speeds. In New Results in Numerical and Experimental Fluid Mechanics XII, pp. 307–316. Cham: Springer International Publishing.
Menke, Matthew & Gursul, Ismet 1997 Unsteady nature ofleading edge vortices. Physics of Fluids 9 (10), 2960–2966.
Pope, S. B. 1975 A more general effective-viscosity hypothesis. Journal of Fluid Mechanics 72 (2), 331–340.
Probst, A & Reuß, S 2016 Progress in scale-resolving simulations with the dlr-tau code. Deutscher Luft- und Raumfahrtkongress, Braunschweig, Deutschland.

Rumsey, C. L., Carlson, J.-R., Pulliam, T. H. & Spalart, P. R.
2020 Improvements to the quadratic constitutive relation
based on nasa juncture flow data. AIAAJ 58.
Shur, M L & al. 2000 Turbulence modeling in rotating andcurved channels:assessing the spalart-shur correction. AIAA
Journal 38 (5).
Shur, Mikhail L., Spalart, Philippe R., Strelets, Mikhail Kh.
& Travin, Andrey K. 2008 A hybrid RANS-LES approachwith delayed-DES and wall-modelled LES capabilities. International Journal of Heat and Fluid Flow 29 (6), 1638–
1649.
Spalart, P. R. 2000 Strategies for turbulence modelling andsimulations. International Journal of Heat and Fluid Flow
21 (3), 252–263.
Spalart, P R, Jou, W H, Strelets, M & Allmaras, S R 1997
Comments on the feasibility of les for wings, and on a hybrid rans/les approach. First AFOSR International Conference on DNS/LES, Ruston, LA.

6