Turbulent Round Jet Flows:

Direct Numerical Simulations

and a Symmetry Analysis
Zur Erlangung des akademischen Grades Doktor-Ingenieur (Dr.-Ing.)
Genehmigte Dissertation von Cat Tuong Nguyen aus Frankfurt am Main
Tag der Einreichung: 29. April 2024, Tag der Prüfung: 24. Juli 2024

1. Gutachten: Prof. Dr.-Ing. M. Oberlack

2. Gutachten: Prof. Dr. rer. nat. M. Schäfer
Darmstadt, Technische Universität Darmstadt

Fachbereich Maschinenbau
Fachgebiet für
Strömungsdynamik
Turbulent Round Jet Flows: Direct Numerical Simulations and a Symmetry Analysis

Genehmigte Dissertation von Cat Tuong Nguyen

Tag der Einreichung: 29. April 2024

Tag der Prüfung: 24. Juli 2024

Darmstadt, Technische Universität Darmstadt

Bitte zitieren Sie dieses Dokument als:

URN: urn:nbn:de:tuda-tuprints-278622
URL: [Link]
Jahr der Veröffentlichung auf TUprints: 2024

Dieses Dokument wird bereitgestellt von tuprints,

E-Publishing-Service der TU Darmstadt
[Link]
tuprints@[Link]

Die Veröffentlichung steht unter folgender Creative Commons Lizenz:

Namensnennung 4.0 International
[Link]
For Family & Friends
Erklärungen laut Promotionsordnung

§ 8 Abs. 1 lit. c PromO

Ich versichere hiermit, dass die elektronische Version meiner Dissertation mit der schriftlichen
Version übereinstimmt.

§ 8 Abs. 1 lit. d PromO

Ich versichere hiermit, dass zu einem vorherigen Zeitpunkt noch keine Promotion versucht
wurde. In diesem Fall sind nähere Angaben über Zeitpunkt, Hochschule, Dissertationsthema
und Ergebnis dieses Versuchs mitzuteilen.

§ 9 Abs. 1 PromO

Ich versichere hiermit, dass die vorliegende Dissertation selbstständig und nur unter Verwen-
dung der angegebenen Quellen verfasst wurde.

§ 9 Abs. 2 PromO

Die Arbeit hat bisher noch nicht zu Prüfungszwecken gedient.

Darmstadt, 29. April 2024

C.T. Nguyen

v
Abstract

A spatially evolving turbulent round jet flow is studied using Lie symmetry analysis and vali-
dated against data from direct numerical simulations (DNS) of the Navier-Stokes equations
(NSEs). The, in the literature, unprecedented simulations were performed at two Reynolds
numbers of Re = 3500 and Re = 7000 based on the orifice diameter D and the bulk velocity
at the orifice Ub with a passive scalar Θ at a Prandtl number of P r = 0.71 in a long numerical
box of z/D = 75.
To achieve self-similarity in the round jet DNS, turbulent pipe flows at the corresponding
Reynolds numbers and a length of z/D = 5 were used as upstream inflow boundary conditions.
A fast convergence of self-similarity was achieved by this approach. High quality statistics
were generated, e.g., by ensemble averaging over 200 washouts of a particle for the simulation
at Re = 3500. Analysis of mean velocities, Reynolds stresses, and turbulent kinetic energy
budgets revealed nearly perfect classical scaling based on a similarity coordinate η = r/z with
the radius r. In addition, statistical analysis of probability density functions (PDFs) for the
axial velocity Uz showed Gaussian behavior along the jet axis, with a transition to heavier tails
and skewness away from the axis.
Using the new DNS data of the high-velocity moments, the statistical behavior of the turbulent
round jet is further analyzed by applying Lie symmetry analysis for pure hydrodynamics. The
turbulent velocity scaling laws derived by Lie symmetry analysis reveal a possible variation of
the turbulent decay behavior due to the statistical symmetry. Interestingly, the statistical sym-
metry is found only in the multi-point moment equations and not in the NSEs. This variation is
possible for all moment orders n except for the second order moment. However, the present
simulations break this symmetry, leading to a classical scaling behavior characterized by η and
a scaling of the velocity moments with z −n , which has been validated up to moment order 10.
The prefactors of the scaling laws are exponential in n. Notably, Gaussian behavior is observed
in the Uz -moments, although the Gaussian exponent shows non-linear behavior in n, which
implies significant intermittency in Uzn . The statistical symmetry, which gives a measure of
intermittency, does not play a role in the base scaling laws, but is the basis for high-moment
scaling for η in turbulent jet flows.
Additionally, moments and PDF statistics of Θ are analyzed with the two simulations. Fur-
thermore, Lie symmetries applied to multi-point velocity-scalar correlation equations lead to
a generalization of the passive scalar and velocity scaling laws, where the scaling prefactors
are exponential for varying velocity and scalar moment order n and m, respectively. Further,
inflow variation is controlled solely by the velocity inlet. Gaussian distribution of instanta-
neous Θ-moments and mixed Uz -Θ-moments with η is observed. Similar to the Uz -moments,
non-linear coefficients in the Gaussian exponent are traced back to intermittency. Unlike the
velocity PDF statistics, scalar PDF statistics deviate from the Gaussian distribution on the jet
axis, but also become skewed and heavy-tailed with increasing η.

vii
Zusammenfassung

Räumlich sich entwickelnde turbulente runde Freistrahlströmungen werden mit der Lie-
Symmetrieanalyse untersucht und mit Daten aus direkten numerischen Simulationen (DNS) der
Navier-Stokes-Gleichungen (NSEs) validiert. Die in der Literatur bisher größten Simulationen
wurden bei Reynoldszahlen von Re = 3500 und Re = 7000 basierend auf dem Düsendurchmes-
ser D und der mittleren Geschwindigkeit durch die Düse Ub mit einem passiven Skalar Θ bei
einer Prandtzahl von P r = 0.71 in einer langen numerischen Box von z/D = 75 durchgeführt.
Um Selbstähnlichkeit in der DNS der Strahlströmung zu erreichen, wurden Geschwindigkeits-
profile turbulenter Rohrströmungen mit einer Länge von z/D = 5 bei den entsprechenden
Reynoldszahlen als Einströmbedingung verwendet. Statistiken von hoher Qualität wurden
generiert, indem, im Falle der Simulation bei Re = 3500, der Ensemblemittelwert über 200
Durchläufe eines Partikels aus der Box gebildet wurde. Analysen der mittleren Geschwindigkeit,
der Reynoldsspannungen und der turbulenten kinetischen Energiebudgets zeigen eine nahezu
perfekte klassische Skalierung für die Ähnlichkeitsvariable η = r/z basierend auf dem Radius
r im Bereich z/D = 25 − 65. Darüber hinaus zeigt die Analyse der Wahrscheinlichkeitsdichte-
funktionen (PDFs) der axialen Geschwindigkeit Uz auf der Mittelachse der Strahlströmung ein
Gaußsches Verhalten. Mit zunehmendem Abstand von der Mittelachse werden die PDFs für
seltene Ereignisse stark nicht-gaussisch und schief.
Weiterhin wurden die neuen DNS-Daten verwendet, um die Geschwindigkeitsmomente höherer
Ordnung mit Hilfe der Lie-Symmetrieanalyse zu untersuchen. Die aus der Symmetrieanalyse
abgeleiteten Skalengesetze der turbulenten Geschwindigkeit weisen auf eine mögliche Variation
des turbulenten Abklingverhaltens hin. Dies wird durch eine statistische Symmetrie ermöglicht,
die nicht in den NSEs, sondern in den Multipunktmomentengleichungen zu finden ist. Bis auf
das zweite Geschwindigkeitsmoment ist die Variation für alle Momentenordnungen möglich.
Die Simulationsdaten bis zur zehnten Momentenordnung brechen jedoch diese Symmetrie
und die Validierung zeigt, dass sich das klassische Skalierverhalten mit η und eine Skalierung
von den Geschwindigkeitsmomenten mit z −n ergibt. Die Vorfaktoren der Geschwindigkeits-
momente sind in n exponentiell. Auffällig ist, dass in den Momenten Uz eine Gauß-Verteilung
beobachtet werden kann und der Gauß-Exponent ein nichtlineares Verhalten in n zeigt, was
eine Intermittenz in Uzn impliziert. Die statistische Symmetrie, die ein Maß für die Intermittenz
ist, hat keinen Einfluss auf die grundlegenden Skalengesetze. Sie ist jedoch wesentlich für die
Skalierung der hohen Momente in turbulenten Strahlströmungen in η.
Lie-Symmetrien angewendet auf die Multipunktgleichung führen zu einer Verallgemeinerung
der Θ- und Geschwindigkeitsmomente und zeigt, dass die Variation des turbulenten Abkling-
verhaltens nur vom Geschwindigkeitseinlass abhängt. Gaußsches Verhalten ist sowohl für Θ
als auch für die gemischten Uz -Θ-Momente zu erkennen. Auch hier sind die Gauß-Exponenten
nichtlinear und können auf Intermittenzen zurückgeführt werden. Im Gegensatz zu Uz -PDFs
sind die Θ-PDFs auf der Mittelachse leicht nicht-gaussisch, werden aber mit zunehmendem η
auch für seltene Ereignisse stark nicht-gaussisch und schief.

ix
Acknowledgements

This dissertation represents the entirety of my efforts during the last three and a half years at
the Chair of Fluid Dynamics. I would like to take this opportunity to express my gratitude to
all those who have contributed to and supported the completion of this work.
First and foremost, I express my heartfelt appreciation to Prof. Dr.-Ing. habil. Martin Oberlack
for his guidance, supervision, and the numerous insightful discussions. His mentorship has
been indispensable, and without his support, this work would not have been possible. His
positive demeanor and encouragement served as a beacon of strength during challenging
phases of my research, and his willingness to share his expertise and connect me with the right
experts was invaluable.
Further, I would like to thank him and the Chair of Fluid Dynamics for the opportunity to
use the infrastructure of the institute, to participate in several international conferences, and
especially for the opportunity to take on the role of the exercise lecturer in “Statics” (TM1). I
fondly remember sitting in the same auditorium 10 years ago and listening to my predecessor.
This experience was particular fulfilling as it allowed me to give back to a new generation of
engineers and contributed greatly to my personal growth.
I am also grateful to Prof. Dr.-Ing. habil. Yongqi Wang for his kindness, as well as his professional
and administrative support. Further, I would like to thank Dr.-Ing. Florian Kummer for his
assistance in IT matters and advice on numerical simulations.
Special thanks are owed to Ruth Völker who has not only been irreplaceable in her administrative
assistance but has also not hesitated to give out advice wherever possible.
Additionally, I extend my appreciation to Prof. Dr. rer. nat. Michael Schäfer for agreeing to
serve as the co-referee for this dissertation.
I am indebted to Paul Hollmann for his assistance in my role as the exercise lecturer and also
for creating a video as well as a beautiful picture of my simulation of a turbulent round jet
flow, which now hangs in the institute.
Moreover, I would also like to express my deepest gratitude to Dr. Hamed Sadeghi for his
support during the initial stages of my research, particularly by sharing his knowledge on Lie
symmetry analysis of temporal evolving plane jet flows. I am also grateful to Dr. Sergio Hoyas
and Dr. Carlos B. da Silva for their expert advice on DNS and their willingness to engage
in digital meetings with me. I also thank Dr. Philipp Schlatter for his generous guidance
in extending a statistics extraction routine for the generated DNS data. It has truly been a
pleasure to meet each of you in person at international conferences.
My time at the Chair of Fluid Dynamics owes much of its pleasantness to the camaraderie and
support of my esteemed colleagues, whose professional and personal connections have greatly

xi
enriched my experience. First of all, I extend my sincere gratitude to Toni Dokoza, my office
mate and long-time friend since the beginning of our bachelor studies, where we started our
academic pursuits together. The PhD journey has been enriched by our many conversations,
not only about our research but also simply sharing anecdotes from our daily lives. I am
equally grateful to Jakob Vandergrift, a friend since elementary school, whose reconnection
during this PhD journey has been profoundly rewarding. Special thanks go to Simon Görtz for
his unwavering support and numerous professional and personal exchanges that have been
invaluable during my time at the institute. I have particularly enjoyed working with Dr.-Ing.
Alparslan Yalcin, with whom I have not only engaged in productive professional discussions,
but also shared many enjoyable conversations about chess. My gratitude extends to Felician
Putz, Jie Liu, Lara De Broeck, Lauritz Beck, Matthias Rieckmann, Schahin Akbari and Weihang
Sun for their contributions to many engaging and fruitful discussions and conversations. Each
interaction has added depth and insight to my academic journey. Finally, I would like to thank
all of my other colleagues who have contributed to the friendly and supportive atmosphere at
the institute, even if they are not explicitly mentioned. Your presence has greatly enriched my
experience and made my time at the Chair of Fluid Dynamics truly memorable.
I am deeply thankful for the extremely helpful feedback provided by Felician Putz, Jakob
Vandergrift, Johannes Conrad, Matthias Rieckmann, Simon Görtz and Toni Dokoza whose
constructive input significantly enhanced the quality of this dissertation.
I extend my deepest gratitude to my father, my mother and my two sisters for their unconditional
support, which transcends all aspects of my life, extending far beyond my academic journey.
Additionally, I am deeply thankful for the encouragement and support of my friends, without
whom this journey would not have been possible. This dissertation is dedicated to my family
and friends, whose individual contributions through their encouragement and support have
been instrumental in its completion. Their belief in me has been a constant source of inspiration
throughout this academic endeavor.
Finally, I want to express my gratitude to the German Academic Scholarship Foundation
(German: Studienstiftung des deutschen Volkes) whose generous support funded three years
of my research. Additionally, I am deeply thankful to the Gauss Centre for Supercomputing e.V.
([Link]) for funding this project under the project number pn73fu by providing
computing time on the GCS Supercomputer SUPERMUC-NG at Leibniz Supercomputing Centre
([Link]).

xii
Contents

List of Tables xiii

List of Figures xv

List of Abbreviations xix

List of Symbols xxi

1. Introduction 1
1.1. Motivation . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 1
1.2. Symmetry-based turbulence theory and round jet flows . . . . . . . . . . . . . 2
1.3. Outline of this work . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 5

2. Governing Equations 7
2.1. Conservation laws . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 7
2.1.1. Continuity equation . . . . . . . . . . . . . . . . . . . . . . . . . . . . 8
2.1.2. Momentum equation . . . . . . . . . . . . . . . . . . . . . . . . . . . . 10
2.2. Navier-Stokes equations . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 10
2.2.1. Navier-Stokes equations in cylindrical coordinates . . . . . . . . . . . . 12
2.2.2. Non-dimensional Navier-Stokes equations . . . . . . . . . . . . . . . . 13
2.3. Convection-diffusion equation . . . . . . . . . . . . . . . . . . . . . . . . . . . 14

3. Statistical Description of Turbulence 17

3.1. Description of random variables . . . . . . . . . . . . . . . . . . . . . . . . . . 17
3.1.1. Gaussian distribution . . . . . . . . . . . . . . . . . . . . . . . . . . . . 19
3.1.2. Statistics in the context of turbulence . . . . . . . . . . . . . . . . . . . 19
3.2. Reynolds-averaged Navier-Stokes equations . . . . . . . . . . . . . . . . . . . 20
3.3. Reynolds stress transport equations . . . . . . . . . . . . . . . . . . . . . . . . 21
3.4. Multi-point moment equations . . . . . . . . . . . . . . . . . . . . . . . . . . . 22
3.5. Multi-point scalar moment equation . . . . . . . . . . . . . . . . . . . . . . . . 23
3.6. Multi-point velocity-scalar correlation equations . . . . . . . . . . . . . . . . . 24

4. Lie Symmetry Theory 27

4.1. Symmetries of algebraic equations . . . . . . . . . . . . . . . . . . . . . . . . . 27
4.2. Symmetries of differential equations . . . . . . . . . . . . . . . . . . . . . . . 28
4.2.1. Lie symmetries of a one-dimensional heat equation . . . . . . . . . . . 30

5. Direct Numerical Simulations of a Turbulent Round Jet Flow at two Reynolds Numbers 33
5.1. Computational domain of the direct numerical simulations . . . . . . . . . . . 37
5.1.1. Computational domain of the pipe flow . . . . . . . . . . . . . . . . . . 38

xiii
 5.1.2. Computational domain of the jet flow . . . . . . . . . . . . . . . . . . . 39
5.2. Mean velocity statistics . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 41
5.3. Turbulent intensities . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 45
5.4. Reynolds stress . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 47
5.5. Reynolds stress transport budgets . . . . . . . . . . . . . . . . . . . . . . . . . 50
5.6. Probability density functions . . . . . . . . . . . . . . . . . . . . . . . . . . . . 55
5.7. Conclusive remarks . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 58

6. Lie Symmetry Analysis of a Turbulent Round Jet Flow 59

6.1. Constructing the turbulent velocity scaling laws . . . . . . . . . . . . . . . . . 60
6.1.1. Symmetries of the multi-point moment equations . . . . . . . . . . . . 60
6.1.2. Derivation of invariant solutions . . . . . . . . . . . . . . . . . . . . . . 61
6.2. Validation of the turbulent velocity scaling laws on the centerline . . . . . . . 64
6.2.1. Results for the axial velocity scaling law . . . . . . . . . . . . . . . . . 65
6.2.2. Results for the radial and azimuthal velocity scaling laws . . . . . . . . 65
6.3. Gaussian behavior of the mean axial velocity . . . . . . . . . . . . . . . . . . . 65
6.3.1. Symmetry reduction of the multi-point moment equations . . . . . . . 69
6.3.2. Hidden intermittency . . . . . . . . . . . . . . . . . . . . . . . . . . . . 70
6.4. Conclusive remarks . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 72

7. Passive Scalar Statistics in a Turbulent Round Jet Flow 75

7.1. Passive scalar statistics . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 77
7.1.1. Mean passive scalar statistics . . . . . . . . . . . . . . . . . . . . . . . 77
7.1.2. Second order passive scalar and mixed statistics . . . . . . . . . . . . . 79
7.1.3. Third order mixed statistics . . . . . . . . . . . . . . . . . . . . . . . . 82
7.1.4. Probability density functions . . . . . . . . . . . . . . . . . . . . . . . 83
7.2. Symmetry-based scaling laws for a passive scalar . . . . . . . . . . . . . . . . 86
7.2.1. Symmetries of the governing equations . . . . . . . . . . . . . . . . . . 86
7.2.2. Derivation of invariant solutions . . . . . . . . . . . . . . . . . . . . . . 87
7.2.3. Validation of the scaling laws on the centerline . . . . . . . . . . . . . 88
7.2.4. Gaussian behavior of the higher passive scalar moments . . . . . . . . 91
7.2.5. Hidden intermittency in passive scalars . . . . . . . . . . . . . . . . . . 96
7.2.6. Gaussian behavior of the higher mixed moments . . . . . . . . . . . . 98
7.3. Conclusive remarks . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 101

8. Conclusion 103
8.1. Closing remarks . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 103
8.2. Outlook . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 105

A. Appendix 115
A.1. Estimation of the grid sizes . . . . . . . . . . . . . . . . . . . . . . . . . . . . 115
A.1.1. Grid size of a turbulent pipe flow . . . . . . . . . . . . . . . . . . . . . 115
A.1.2. Grid size of a turbulent round jet flow . . . . . . . . . . . . . . . . . . 115
A.2. Self-preservation of additional moments . . . . . . . . . . . . . . . . . . . . . 116
A.3. Symmetry reduction of the multi-point scalar equation . . . . . . . . . . . . . 121

xiv
List of Tables

5.1. List of 75 variables that have been computed with the extension of the statistics
toolbox. Note that these variables have been computed on top of the 44 variables
already calculated by the toolbox (Rezaeiravesh et al., 2019). The position
describes where they are stored in the statistical data output. . . . . . . . . . 37
5.2. Boundary conditions (BCs) of the main computational domain . . . . . . . . . 40
5.3. Various jet parameters of DNS and experiments. . . . . . . . . . . . . . . . . . 43

7.1. BCs of the main computational domain for a passive scalar . . . . . . . . . . . 77

xv
List of Figures

5.1. Schematic view of a round jet flow. Fluid is blown through a nozzle with diameter
D. The mean axial velocity of the jet is denoted by U z while the mean axial
centerline velocity is denoted by U z,c . . . . . . . . . . . . . . . . . . . . . . . . 33
5.2. Cross-sectional view of the computational domain for the pipe at Re = 3500.
The N = 7 Gauss-Lobatto-Legendre (GLL) points has been included in the mesh. 38
5.3. A cross-sectional view of the pseudo-color visualized magnitude of the instan-
taneous velocity field at Re = 3500 (left) and Re = 7000 (right). The values
range from 0 (blue) to 1.4 (red). . . . . . . . . . . . . . . . . . . . . . . . . . 39
5.4. Cross-sectional view of the main computational box for the Re = 3500 case at
z/D = 0 (left) which is scaled linearly in z-direction to obtain the whole main
computational box (right). The GLL points are omitted for better visibility. The
DNS uses curved elements which are not apparent in this figure. . . . . . . . . 41
5.5. A cross section of q-criterion isosurfaces at q = 0.01 of the conducted jet DNS at
Re = 3500 colored with the velocity magnitude. . . . . . . . . . . . . . . . . . 41
5.6. Inverse of the mean axial centerline velocity according to Equation (5.4) plotted
over the distance from the orifice. Present DNS at Re = 3500 ( ), present
DNS at Re = 7000 ( ), Boersma et al. (1998) at Re = 2400 ( ), Babu and
Mahesh (2004) at Re = 2400 ( ), Taub et al. (2013) at Re = 2000 ( ). The
inset figure with the same axis shows a magnification at the near-field, where
the cutoff ( ) marks the end of the potential core. Both present simulations
there have an earlier onset of spreading compared to the comparative studies. 42
5.7. Mean axial velocity profiles normalized with the axial centerline velocity accord-
ing to Equation (5.7) as functions of η at different distances from the orifice for
Re = 3500 (left) and Re = 7000 (right): z/D = 25 ( ), 35 ( ), 45 ( ),
55 ( ), 65 ( ). . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 43
5.8. Mean radial velocity profiles normalized with the axial centerline velocity ac-
cording to Equation (5.7) as functions of η at different distances from the orifice
for Re = 3500 (left) and Re = 7000 (right): z/D = 25 ( ), 35 ( ), 45
( ), 55 ( ), 65 ( ). . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 44
5.9. The invariant U z (η = 0) plotted over the distance from the orifice. Present DNS
e
at Re = 3500 ( ), present DNS at Re = 7000 ( ), Boersma et al. (1998)
( ), Babu and Mahesh (2004) ( ), Taub et al. (2013) ( ). . . . . . . . 45
[Link] component of the turbulent intensity on the centerline compared to
previous studies: present DNS at Re = 3500 ( ), present DNS at Re = 7000
( ), Taub et al. (2013) ( ), Bogey and Bailly (2009) ( ), Panchapakesan
and Lumley (1993a) ( ). . . . . . . . . . . . . . . . . . . . . . . . . . . . . 46

xvii
 [Link] component of the turbulent intensity on the centerline compared to previ-
ous studies: present DNS at Re = 3500 ( ), present DNS at Re = 7000 ( ),
Taub et al. (2013) ( ), Bogey and Bailly (2009) ( ), Panchapakesan and
Lumley (1993a) ( ). . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 46
2
[Link] stresses ui uj normalized with the axial centerline velocity accord-
U z,c
ing to Equation (5.8) at different distances from the orifice for the Re = 3500
case: z/D = 25 ( ), 35 ( ), 45 ( ), 55 ( ), 65 ( ). . . . . . . . . . 48
2
[Link] stresses ui uj normalized with the axial centerline velocity accord-
U z,c
ing to Equation (5.8) at different distances from the orifice for the Re = 7000
case: z/D = 25 ( ), 35 ( ), 45 ( ), 55 ( ), 65 ( ). . . . . . . . . . 49
2
[Link] of the Reynolds stresses at Re = 3500 scaled with U z,c /z at z/D = 25:
production ( ), dissipation ( ), turbulent diffusion ( ), velocity-pressure
gradient correlation ( ), convection ( ), pressure diffusion ( ), pressure
strain ( ). . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 52
3
[Link] of the Reynolds stresses at Re = 7000 scaled with U z,c /z at z/D = 35:
production ( ), dissipation ( ), turbulent diffusion ( ), velocity-pressure
gradient correlation( ), convection ( ), pressure diffusion ( ), pressure
strain ( ). . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 53
[Link] kinetic energy budgets for the Re = 3500 case at z/D = 25 (left)
and for the Re = 7000 case at z/D = 35 (right): production ( ), dissipation
( ), turbulent diffusion ( ), convection ( ), pressure diffusion ( ),
sum of the budget terms ( ). . . . . . . . . . . . . . . . . . . . . . . . . . . 54
[Link] dissipation of the turbulent kinetic energy comparing the present DNS at
Re = 3500 and Re = 7000 to the DNS of Taub et al. (2013) at Re = 2000, the
experiments of Panchapakesan and Lumley (1993a) at Re = 11,000, Hussein
et al. (1994) at Re = 95,500 and a LES of Bogey and Bailly (2009) at Re = 11,000. 55
[Link] of Uz (η = 0, z)/U z,c (z) at Re = 3500 (left) and Re = 7000 (right) for
z/D = 15 ( ), 25 ( ), 35 ( ), 45 ( ), 55 ( ), 65 ( ) compared to
a Gaussian ( ). . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 56
[Link] of Uz (η, z)/U z,c (z) at Re = 3500 (above) and Re = 7000 (below) for
z/D = 28 ( ), 42 ( ), 56 ( ). . . . . . . . . . . . . . . . . . . . . . . . 57
[Link] K (above) and skewness S (below) of Uz (η, z)/U z,c (z) at Re = 3500
(left) and Re = 7000 (right) for z/D = 28 ( ), 42 ( ), 56 ( ). K = 3,
S = 0 (dashed) are the Gaussian values. . . . . . . . . . . . . . . . . . . . . . 57

6.1. The exponential prefactors αi,n eci n of the DNS according to the velocity scaling
law on the centerline (6.28) are shown for each moment up to order n = 10. For
Re = 3500, they are presented for each direction i = r ( ), z ( ) and determined
with the DNS data ( ). Likewise, for Re = 7000, they are displayed for each
direction i = r ( ), z ( ) and determined with the DNS data ( ). . . . . . . . 66
6.2. The radial profiles of the nth axial moment normalized with the scalings in
(6.25) at different distances from the orifice for the Re = 3500 case: z/D = 25
( ), 35 ( ), 45 ( ), 55 ( ), 65 ( ). The black solid lines indicate the
Gaussian behavior from (6.36) using γn from (6.37). . . . . . . . . . . . . . . 67

xviii
6.3. The radial profiles of the nth axial moment normalized with the scalings in
(6.25) at different distances from the orifice for the Re = 7000 case: z/D = 25
( ), 35 ( ), 45 ( ), 55 ( ), 65 ( ). The black solid lines indicate the
Gaussian behavior from (6.36) using γn from (6.37). . . . . . . . . . . . . . . 68
6.4. The radial profiles (blue, dashed) of the 1st (top) up to the 10th (bottom) axial
moment and the corresponding curves (black) from (6.36) shown in a semi-
logarithmic plot at z/D = 45 for the Re = 3500 case (left) and the Re = 7000
case (right). . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 69
6.5. Constants γn from (6.36) are shown for each moment up to order n = 10
determined by fitting to the DNS yielding γn = −1.54n2 + 61.95n + 27.10
(Re = 3500) and γn = −1.51n2 + 63.39n + 28.21 (Re = 7000). . . . . . . . . . 70
6.6. The radial profiles of for the Re = 3500 at different distances z from the
Ur Uzn−1
orifice compared to the solution in (6.40) ( ): z = 25 ( ), z = 35 ( ),
z = 45 ( ), z = 55 ( ), z = 65 ( ). . . . . . . . . . . . . . . . . . . . . 71

7.1. Mean inverse passive scalar at the centerline over the distance of the orifice:
Birch et al. (1978) at Re = 16,000 ( ), Babu and Mahesh (2005) at Re = 2400
( ), Lubbers et al. (2001) at Re = 2000 ( ), present DNS at Re = 3500
( ), present DNS at Re = 7000 ( ). In the magnified view, the potenial
core of the passive scalar is marked by a dashed line ( ). . . . . . . . . . . . 78
7.2. Mean passive scalar Θ scaled with the centerline passive scalar Θc (7.4) for
the Re = 3500 case (left) plotted as a function of the similarity coordinate η
at different distances from the orifice: z/D = 15 ( ), 25 ( ), 35 ( ), 45
( ), 55 ( ). For the Re = 7000 case (right), Θ/Θc is plotted at: z/D = 20
( ), 30 ( ), 40 ( ), 50 ( ), 60 ( ). . . . . . . . . . . . . . . . . . . 79
7.3. Variance of the passive scalar fluctuations RΘΘ at Re = 3500 (left) scaled with
the centerline passive scalar Θ2c (7.4) plotted as a function of the similarity
coordinate η at different distances from the orifice: z/D = 25 ( ), 35 ( ),
45 ( ), 55 ( ). For the Re = 7000 case (right), RΘΘ /Θ2c is plotted at:
z/D = 30 ( ), 40 ( ), 50 ( ), 60 ( ). . . . . . . . . . . . . . . . . . . 80
√
7.4. Centerline root-mean square (rms) of the passive scalar fluctuation RΘΘ scaled
with the centerline mean passive scalar Θc (7.4) over the distance z from the
orifice: Darisse et al. (2015) ( ), Birch et al. (1978) ( ), Babu and Mahesh (2005)
( ), Lubbers et al. (2001) ( ), present DNS at Re = 3500 ( ), present
DNS at Re = 7000 ( ). . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 80
7.5. Turbulent heat fluxes RrΘ and RzΘ normalized with Uz,c Θc for the Re = 3500
case at different distances from the orifice: z/D = 25 ( ), 35 ( ), 45 ( ),
55 ( ). . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 81
7.6. Turbulent heat fluxes RrΘ and RzΘ normalized with Uz,c Θc for the Re = 7000
case at different distances from the orifice: z/D = 30 ( ), 40 ( ), 50 ( ),
60 ( ). . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 81
7.7. Turbulent heat fluxes RrΘΘ and RzΘΘ at Re = 3500 normalized with Uz,c Θ2c
at different distances from the orifice: z/D = 25 ( ), 35 ( ), 45 ( ), 55
( ). . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 82
7.8. Turbulent heat fluxes RrΘΘ and RzΘΘ at Re = 7000 normalized with Uz,c Θ2c at
different distances from the orifice: z/D = 30 ( ), 40 ( ), 50 ( ), 60 ( ) 83

xix
 7.9. PDFs of Θ(η = 0, z)/Θc (z) for the Re = 3500 case (left) and Re = 7000 case
(right) at z/D = 15 ( ), 25 ( ), 35 ( ), 45 ( ), 55 ( ), 65 ( )
compared to a Gaussian ( ). . . . . . . . . . . . . . . . . . . . . . . . . . . 84
[Link] of Θ(η, z)/Θc (z) at Re = 3500 (above) and Re = 7000 (below) for z/D =
28 ( ), 42 ( ), 56 ( ). . . . . . . . . . . . . . . . . . . . . . . . . . . . 84
[Link] K (above) and skewness S (below) of Θ(η, z)/Θc (z) at Re = 3500 for
z/D = 28 ( ), 42 ( ), 56 ( ). K = 3, S = 0 (dashed) are the Gaussian
values. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 85
[Link] K (above) and skewness S (below) of Θ(η, z)/Θc (z) at Re = 7000 for
z/D = 28 ( ), 42 ( ), 56 ( ). K = 3, S = 0 (dashed) are the Gaussian
values. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 85
[Link] of Θ(η, z = 28)/Θc (z = 28) for η = 0.139 (solid) and η = 0.149 (dashed). 86
[Link] exponential prefactor αzΘ,nm ecz,nm (n+m) ( ) from Equation (7.19) deter-
mined with the DNS data at Re = 3500 for mixed moments up to n + m = 6 ( ),
for Θ moments up to m = 10 ( ) and Uz moments up to n = 10 ( ) is shown.
Additionally, Equation (7.19) is highlighted for m = 0 ( ) and n = 0 ( ). 90
[Link] exponential prefactor αrΘ,nm e c r,nm (n+m) ( ) from (7.19) determined with
the DNS data at Re = 3500 for mixed moments up to n + m = 6 ( ), for
Θ moments up to m = 10 ( ) and Ur moments up to n = 10 ( ) is shown.
Additionally, Equation (7.19) is highlighted for m = 0 ( ) and n = 0 ( ). 91
th
[Link] radial profiles of the m axial moment normalized with the scalings in
(7.17) at different distances from the orifice for the Re = 3500 case: z/D = 25
( ), 35 ( ), 45 ( ), 55 ( ). The black solid lines indicate the Gaussian
from Equation (7.26) using γm from Equation (7.27). . . . . . . . . . . . . . . 93
[Link] radial profiles of the mth axial moment normalized with the scalings in
(7.17) at different distances from the orifice for the Re = 7000 case: z/D = 30
( ), 40 ( ), 50 ( ), 60 ( ). The black solid lines indicate the Gaussian
from Equation (7.26) using γm from Equation (7.27). . . . . . . . . . . . . . . 94
[Link] radial profiles (blue, dashed) of the 1st (top) up to the 10th (bottom) axial
moment at Re = 3500 (left) and Re = 7000 (right) and the corresponding
Gaussian (black) from Equation (7.26) shown in a semi-logarithmic plot at
z/D = 45. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 95
[Link] γm from Equation (7.26) are shown for each moment up to order
n = 10 determined by fitting to the DNS yielding the following fits: γm =
−1.27m2 + 37.52m + 27.34 (Re = 3500) and γm = −1.39m2 + 39.59m + 24.04
(Re = 7000). . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 95
[Link] radial profiles (blue, dashed) of m = 1 (top) up to the n + m = 6 (bottom)
axial mixed moment Uzn Θm at Re = 3500 and the corresponding Gaussian
(black) from Equation (7.26) shown in a semi-logarithmic plot at z/D = 45. . 99
[Link] radial profiles (blue, dashed) of m = 1 (top) up to the n + m = 6 (bottom)
axial mixed moment Uzn Θm at Re = 7000 and the corresponding Gaussian
(black) from Equation (7.26) shown in a semi-logarithmic plot at z/D = 45. . 100
[Link] γnm from Equation (7.42) at Re = 3500 are shown for pure moments
up to order n, m = 10 and mixed moments up to order n + m = 6 determined
by fitting to the DNS yielding the following fit: γnm = −1.51n2 − 1.29m2 −
0.03nm + 61.2n + 37.34m + 30.75. . . . . . . . . . . . . . . . . . . . . . . . . 101

xx
A.1. The radial profiles of the nth axial moment normalized with the scalings in
(6.25) at different distances from the orifice for the Re = 3500 case: z/D = 25
( ), 35 ( ), 45 ( ), 55 ( ), 65 ( ). The black solid lines indicate the
Gaussian behavior from (6.36) using γn from (6.37). . . . . . . . . . . . . . . 117
A.2. The radial profiles of the nth axial moment normalized with the scalings in
(6.25) at different distances from the orifice for the Re = 7000 case: z/D = 25
( ), 35 ( ), 45 ( ), 55 ( ), 65 ( ). The black solid lines indicate the
Gaussian behavior from (6.36) using γn from (6.37). . . . . . . . . . . . . . . 118
A.3. The radial profiles of the mth axial moment normalized with the scalings in
(7.17) at different distances from the orifice for the Re = 3500 case: z/D = 25
( ), 35 ( ), 45 ( ), 55 ( ). The black solid lines indicate the Gaussian
from Equation (7.26) using γm from Equation (7.27). . . . . . . . . . . . . . . 119
A.4. The radial profiles of the mth axial moment normalized with the scalings in
(7.17) at different distances from the orifice for the Re = 7000 case: z/D = 30
( ), 40 ( ), 50 ( ), 60 ( ). The black solid lines indicate the Gaussian
from Equation (7.26) using γm from Equation (7.27). . . . . . . . . . . . . . . 120

xxi
List of Abbreviations

BC boundary condition
BDF2 second order backward differentiation formula

CFD computational fluid dynamics

DNS direct numerical simulation

DOF degree of freedom

GLL Gauss-Lobatto-Legendre

LES large-eddy simulation

LMN Lundgren-Novikov-Monin

MPI Message Passing Interface

MPME multi-point moment equation
MPSE multi-point scalar equation
MPVSCE multi-point velocity-scalar correlation equation

NSE Navier-Stokes equation

ODE ordinary differential equation

PDE partial differential equation

PDF probability density function

RANS Reynolds-averaged Navier-Stokes

rms root-mean square
RST Reynolds stress tensor

TKE turbulent kinetic energy

xxiii
List of Symbols

( · )(0) correlation at one point indicated by a superscript

a group parameter
(·) mean quantity
α thermal diffusivity
αi,n scaling law parameter of velocity moments
αiΘ,nm scaling law parameter of mixed velocity-scalar mo-
ments
αΘ,m scaling law parameter of passive scalar moments
( · )[·] brackets, meaning no implied summation
BΘ decay rate of the passive scalar
Bu decay rate of the axial velocity
c speed of sound
Cij convection term in the Reynolds stress transport equa-
tions
ci concentration quantity i
ci,n scaling law parameter of velocity moments
ci,nm scaling law parameter of mixed velocity-scalar mo-
ments
Ck class of functions that have a k th derivative that is
continuous in its domain
cm scaling law parameter of passive scalar moments
( · )c value at the centerline
D diameter
∆r+ radial grid spacing in wall units
δv viscous length scale
∆x Cartesian grid spacing in x-direction
∆y Cartesian grid spacing in y-direction
∆z + axial grid spacing in wall units
∆Rϕ+ grid spacing on an arc with radius R in wall units
Di diffusion coefficient of species i
dij deviatoric stress tensor
( · )!! double factorial
eij rate-of-strain tensor
ij dissipation term in the Reynolds stress transport equa-
tions
η normalized radial coordinate
η1/2 jet half-width
ηk infinitesimal k

xxv
 fD friction factor
f (V ) probability density function
γm Gaussian prefactor of the passive scalar moment pro-
files
γn Gaussian prefactor of the axial velocity moment pro-
files
γnm Gaussian prefactor of the mixed velocity-scalar mo-
ment profiles
H velocity moment based on instantaneous values
Ii momentum
κ bulk viscosity
ki volume force
Kn Knudsen number
δij Kronecker delta
K kurtosis
L length of a flat plate
Lij local temporal term in the Reynolds stress transport
equations
λ mean free path
Lc characteristic length scale
m mass
Ma Mach number
( · )max maximum value
( · )min minimum value
µ dynamic viscosity
µ00 second coefficient of viscosity
∇ nabla operator
ni surface normal vector
nc,i concentration flux of i
nS unit normal vector to a surface S
ν kinematic viscosity
O order of magnitude
P pressure
Pe Péclet number
Πij velocity-pressure gradient correlation term in the
Reynolds stress transport equations
Πdij pressure diffusion term in the Reynolds stress trans-
port equations
Πsij pressure strain term in the Reynolds stress transport
equations
Pr Prandtl number
Pij production term in the Reynolds stress transport equa-
tions
P∗ pressure rescaled by a constant density
R velocity moment based on fluctuation values. Rij is
also known as the Reynolds stress tensor
Re Reynolds number

xxvi
Reτ friction Reynolds number
ρ density
r+ wall unit
Sb boundary
σ standard deviation
S skewness
t time
Tij turbulent diffusion term in the Reynolds stress trans-
port equations
τij Cauchy stress tensor
Θ passive scalar quantity
Θi passive scalar quantity i
ti surface force
(f·) similarity variables or invariants
tr(·) trace of a square matrix
U flow velocity vector
u velocity fluctuation
U mean velocity
Ub bulk velocity
Uc characteristic velocity scale
uc,i global speed of concentration ci
Ui velocity component i
U∞ free-stream velocity
ul local velocity
uτ friction velocity
V volume
Vij viscous diffusion term in the Reynolds stress transport
equations
vn central statistical moment of order n
Vn statistical moment of order n
X infinitesimal operator
xi spatial coordinate i
ξi infinitesimal i
X (n) nth prolonged operator
z0 virtual origin

xxvii
1. Introduction

This dissertation represents the entirety of my work over the past three and a half years at
the Chair of Fluid Dynamics at the Technical University of Darmstadt. Notably, it incorporates
content from several papers, which have been slightly modified. These papers are a direct
result of my contributions as a doctoral researcher:
• Nguyen, C. T. and Oberlack, M. (2024a). “Analysis of a turbulent round jet based on
direct numerical simulation data at large box and high Reynolds number”. Physical Review
Fluids 9, 074608
• Nguyen, C. T. and Oberlack, M. (2024b). “Comparative Study of Turbulent Round Jet
Flows through Direct Numerical Simulation at Medium-High Reynolds Numbers”. Under
review with Physics of Fluids
• Nguyen, C. T. and Oberlack, M. (2024c). “Hidden intermittency in turbulent jet flows”.
Under review with Physical Review Research
• Nguyen, C. T. and Oberlack, M. (2024d). “Passive scalar statistics in a turbulent round
jet: symmetry theory and direct numerical simulation”. Under review with Journal of
Fluid Mechanics

1.1. Motivation

Turbulence, characterized by its complex and seemingly unpredictable fluid motion, remains
one of the most challenging and fascinating phenomena in fluid dynamics. Despite decades of
intensive research, turbulence continues to defy complete understanding, making it an essential
area for investigation. The study of turbulence holds profound implications for numerous
practical applications across various scientific and engineering disciplines.
Pioneering works by Kolmogorov et al. (1991), a translation of their work published in 1941,
laid the foundation for the theory of turbulence, introducing the concept of energy cascade
and scaling laws, in specific the 5/3 law, that govern the behavior of turbulent flows at small
scales. The insights of Kolmogorov into the statistical properties of turbulence have guided
generations of researchers in their quest to unravel the mysteries of turbulent motion.
One of the fundamental concepts in turbulence research is self-similarity, which describes the
property of turbulent flows exhibiting similar statistical characteristics across a wide range of
length and time scales. Self-similarity analysis provides a powerful framework for identifying
underlying patterns and structures within turbulent flows, enabling researchers to develop
predictive models and enhance their understanding of complex flow phenomena.

1
Probably the most well-known example of similarity in turbulence is the algebraic decay motion
of the kinetic energy in isotropic turbulence addressed in the ground-breaking paper by Karman
and Howarth (1938). These ideas have been extensively employed in the literature and an
overview including a new exponential decay mode may be taken from George (2009). Other
problems, such as homogeneous or wall-bounded shear turbulence have also been tackled using
similarity methods though to a considerably reduced extend. An overview on various classical
similarity problems in turbulence may be taken from classical text books such as from Pope
(2000). Unfortunately, many of the publications do not make use of any first principle basis
such as the employment of the Navier-Stokes equations (NSEs) or any statistical equations that
can be derived from the former.

A canonical example for classical similarity problems to study fundamental turbulence phe-
nomena involves the turbulent round jet flow, characterized by its axisymmetric nature and
turbulent mixing properties. The dynamics of turbulent jets exhibit remarkable self-similar
behavior, with flow structures and statistical properties remaining consistent across a wide
range of scales. Experimental studies, such as those conducted by Wygnanski and Fiedler
(1969) and Hussein et al. (1994), have provided valuable insights into the structure and evolu-
tion of turbulent round jets. These studies have highlighted the importance of self-similarity
analysis in understanding the underlying dynamics of turbulent flows and have paved the way
for further research in this area.

The relevance of turbulent round jets extends beyond academic curiosity, finding applications
in various engineering fields. In aerospace engineering, turbulent jets play a critical role in
jet propulsion systems, where understanding and optimizing jet dynamics are essential for
enhancing engine performance and fuel efficiency. In environmental engineering, turbulent
jets are central to the dispersion of pollutants in the atmosphere and water bodies, influencing
air quality and ecosystem health.

Furthermore, turbulent round jets serve as a testing ground for validating turbulence models
and simulation techniques, providing benchmarks for assessing the accuracy and reliability of
numerical predictions as seen in the recent dissertation of Klingenberg (2022). By refining our
understanding of turbulent jet dynamics, researchers can improve the predictive capabilities of
turbulence models, leading to more efficient design strategies and safer engineering practices.

In summary, the study of turbulence and self-similarity analysis holds both theoretical signif-
icance and practical importance across a range of scientific and engineering disciplines. By
investigating turbulent round jets and applying self-similarity analysis, researchers can deepen
their understanding of turbulent flows, advance predictive modeling capabilities, and drive
innovation in various fields.

1.2. Symmetry-based turbulence theory and round jet flows

In the theory of turbulent jets, efforts have been dedicated to find scaling parameters using
experimental and numerical data with the aim of accurately describing and modeling real jet
flows. In the theoretical approach, self-similarity analysis has been shown to be very useful in
finding these parameters as mentioned in the previous section.

2
Seminal contributions by Townsend (1956, 1976) have played a large role in this regard,
as they identified classical scaling parameters using self-similarity analysis. Moreover, they
postulated the universality of self-similar solutions, i.e., self-similar solutions are unique and
independent of initial conditions. This implies that for jets, self-similarity is expected with
their scaling laws being a strong attractor.

Early experiments by Bradbury (1965), Heskestad (1965) and Gutmark and Wygnanski (1976)
revealed variations in self-similar profiles of turbulent properties, challenging the notion of
universality. George (1989) subsequently proposed that self-similarity is not universal and,
therefore, depends on initial conditions, contrary to what is conveyed by Townsend (1956).

Further, George introduced equilibrium similarity analysis, a novel approach to determine

similarity scales from statistical equations. Over the past three decades, significant progress in
similarity methods applied to turbulence has been made with the method of equilibrium analysis.
This method has proven highly effective in tackling diverse turbulence problems, including
isotropic turbulence (George, 2009), canonical wall-bounded shear flows like boundary layers
(George and Castillo, 1997), and channel and pipe flow turbulence (Wosnik et al., 2000), among
others. In a subsequent work (George, 2012), he then proposed that the equilibrium similarities
are equally valid for wakes and jets. It gained widespread acceptance in experimental and
numerical studies of jets such as those by Sadeghi et al. (2014), Sadeghi and Pollard (2012),
and Xu and Antonia (2002). However, the reliance of the method on an extended ansatz
approach means it does not fully unravel the underlying mechanisms by using first principles
to construct solutions for turbulent flow statistics. Consequently, this remains an open question,
one which this dissertation aims to address through Lie symmetry theory.

Einstein’s groundbreaking work on special relativity in the early 20th century highlighted
Lie symmetry principles as a fundamental feature of physics, constraining the permissible
laws governing physical systems. As quantum mechanics emerged in the 1920s, symmetries
became established as the foundational framework of physics, serving as the cornerstone for
understanding and mathematically modeling new physical laws yet to be uncovered (see, e.g.,
D. J. Gross, 1996; K. Brading and E. Castellani, 2003).

Today, symmetries stand as the primary guiding principle in the exploration of fundamental
physics, offering insights into the structure of the universe. A global community of researchers
actively extends the principles of Lie symmetries, exploring various possible paths of develop-
ment. These include the investigation of generalized symmetries, approximate symmetries,
non-classical or conditional symmetries, and the establishment of connections to the Painlevé
test and inverse scattering theory, among others. A comprehensive overview of the litera-
ture, particularly concerning the application of these methods, can be found in the works of
Ibragimov (1994, 1995, 1996).

Despite the widespread recognition of Lie symmetries in theoretical physics, their axiomatic im-
portance in statistical turbulence theory has often been overlooked. In many cases, researchers
have implicitly utilized Lie symmetries to construct ansatz-based self-similar solutions, without
explicitly referencing their elementary basis. This gap presents an opportunity for further
exploration, as the application of Lie symmetries in turbulence theory could offer a deeper
understanding of turbulent flows and potentially lead to the development of more rigorous
and comprehensive theoretical frameworks.

3
In the mid-1990s, significant advancements were made in turbulence theory, particularly in the
development of a Lie symmetry based approach utilizing the infinite set of multi-point moment
equations (MPMEs). This work reached its culmination with the habilitation dissertation by
Oberlack (2000a), accompanied by several subsequent publications (Oberlack, 1999, 2001;
Oberlack and Guenther, 2003; Oberlack and Waclawczyk, 2006). These contributions success-
fully unified turbulent scaling laws across various shear flows, including those with rotation,
although the scaling laws were limited to the mean velocity at that time.

In the study of turbulence, the Reynolds decomposition Ui = U i + ui is commonly employed,

where ( · ) denotes ensemble averaging and ui signifies the fluctuation. Within the Reynolds-
averaged Navier-Stokes (RANS) equation, the moment Hij = Ui Uj decomposes into Hij =
U i U j + ui uj , with Rij = ui uj representing the Reynolds stress tensor (RST). The traditional
R-approach, introduced by the Reynolds decomposition, relies on fluctuation correlations,
constituting a non-linear and non-local set of equations coupled between different velocity
moment orders such as in the RANS equation or the Reynolds stress transport equations. In
contrast, the H-approach, where the Reynolds decomposition is omitted, presents an infinite set
of linear equations, allowing for a more convenient treatment in the context of symmetries due
to coupling only between neighboring equations. Further details can be found in Section 3.1.2.
This approach has led to the derivation of new scaling symmetries and an infinite set of new
statistical symmetries in Oberlack and Rosteck (2010).

The first symmetry of this infinite set was probably found by Kraichnan (1965), who called it
random Galilean invariance, which is significant for Kolmogorov’s 5/3-law, in the context of
turbulence modeling. Even the very early turbulence models up to the most recently developed
statistical models are essentially all consistent with this symmetry.

The dissertation by Rosteck (2013) addressed several complex mathematical challenges within
the MPMEs. His key insight was demonstrating the weak coupling between the order of
the equations and therefore of the symmetries. This allowed for a systematic approach to
computing symmetries, starting from lowest order equations and gradually extending to higher
order ones. These insights paved the way for the derivation of additional statistical symmetries,
particularly in wall-bounded turbulent shear flows for not only the mean velocity but also for
the second order correlations.

However, in Wacławczyk et al. (2014) it was shown that the H-approach is not sufficient to
understand statistical symmetries. There, the statistical symmetries, only admitted by the
MPMEs and not the NSEs, where transferred to the Lundgren-Novikov-Monin (LMN) multi-
point probability density function (PDF) equation of turbulence (see Lundgren, 1967; Monin,
1967; Novikov, 1968). The LMN approach revealed that the statistical scaling symmetry
corresponds to a measure of intermittency, while the statistical translation symmetry describes
the non-Gaussianity of the PDF.

The practical implications of these theoretical advancements are showcased by Avsarkisov et al.
(2014) through various large-scale direct numerical simulations (DNS), validating logarithmic
scaling laws in turbulent Poiseuille flow with wall-transpiration. Further, just recently in
Oberlack et al. (2022), turbulent scaling laws of arbitrary H-moments in the log and core
region of a turbulent channel flow were derived and validated using a new DNS at a friction
Reynolds number of Reτ = 104 (Hoyas et al., 2022). Key to the scaling laws is the statistical
scaling symmetry, describing the measure of intermittency.

4
First efforts of the derivation of scaling laws for jet flows have been conducted by Sadeghi
et al. (2018, 2021) for temporally evolving turbulent plane jets. These efforts resulted in the
derivation of scaling laws up to second-order moments of the velocity and a passive scalar,
demonstrating remarkable agreement with DNS data. Alcántara-Ávila et al. (2024) expanded
on this with symmetry-based turbulent scaling laws for streamwise velocity and temperature
moments of arbitrary order which have been validated by DNS at various Reynolds numbers.

Subsequently, this dissertation directs its attention towards spatially evolving turbulent round
jet flows, recognizing their canonical nature and importance in various engineering applications,
as highlighted earlier. In this dissertation, we embark on conducting the two largest DNS of
spatially evolving turbulent round jet flows to date at the two Reynolds numbers 3500 and
7000 based on the orifice diameter and the bulk velocity at the orifice. In addition, we extract
statistical data of velocity and passive scalar moments up to the tenth order. Symmetry-induced
scaling laws are derived up to moments of arbitrary order for this flow, and rigorously validated
against the DNS data. This investigation aims to deepen our understanding of turbulent round
jet flows, potentially leading to improvements in future turbulence models and providing
insights into manipulating flow characteristics for specific engineering objectives.

1.3. Outline of this work

The subsequent section provides an outline of this dissertation.

In Chapter 2 the theoretical groundwork of the present dissertation is established, beginning

with the derivation of the NSEs in Section 2.2. Following this, in Section 2.2.1, the NSEs
are transformed from Cartesian to cylindrical coordinates due to the axisymmetric nature of
turbulent round jet flows, and in Section 2.2.2, they are also non-dimensionalized. Finally, in
Section 2.3, the convection-diffusion equation is derived to describe a passive scalar.

Turbulence and its statistical description is introduced in Chapter 3. In Section 3.1, fundamental
concepts of probability theory are applied to continuous and discrete random variables. The
newly gained knowledge is then consolidated in the example of a Gaussian distribution in
Section 3.1.1. In Section 3.1.2, probability theory is then applied to turbulence. Subsequently,
equations from the R-approach, mainly the RANS equations in Section 3.2 and the Reynolds
stress transport equations in Section 3.3, are derived and compared to the equations from
the H-approach being the MPMEs in Section 3.4, multi-point scalar equations (MPSEs) in
Section 3.5 and multi-point velocity-scalar correlation equations (MPVSCEs) in Section 3.6.

In Chapter 4, the relevant concepts of the Lie symmetry theory for this dissertation are briefly
presented. Section 4.1 focuses on the derivation of symmetries for algebraic equations while
Section 4.2 extends the derivation of symmetries to differential equations.

The subsequent three chapters present the results of this dissertation. The results of the
large-scale DNS at two Reynolds numbers are presented in Chapter 5. The computational
domain of the DNS is specified in Section 5.1, followed by a detailed discussion of various types
of statistical data at both Reynolds numbers. The DNS data includes statistics up to third order
and the PDF for the axial velocity. Notably, and not previously reported, the PDF is plotted
over the radius to investigate self-similarity.

5
Chapter 6 employs the Lie symmetry method to analyze turbulent round jet flows. Invariant
solutions, also called turbulent scaling laws, are derived from Lie symmetries of the MPMEs
in Section 6.1. In Section 6.2 the turbulent scaling laws are then validated using the afore-
mentioned DNS at both Reynolds numbers and compared. Lastly, in Section 6.3, the Gaussian
behavior of the axial H-moments is explored and derived from the symmetries of the reduced
MPMEs using the turbulent scaling laws.
The discoveries in the two preceding chapters are consolidated in Chapter 7 to extend our
knowledge to passive scalars. Keeping the structure of those chapters, the statistical data of first
to third order are presented from the conducted DNS at both Reynolds numbers in Section 7.1.
This also includes mixed velocity and passive scalar moments and is complemented by the
PDF of the passive scalar. In Section 7.2, turbulent scaling laws for a passive scalar are derived
using the Lie symmetry method and generalized for both velocity and passive scalar moments,
alongside mixed moments. Again, the turbulent scaling laws are validated against the DNS
data and the results at both Reynolds numbers are compared. Interestingly, the Gaussian
behavior, discussed in Section 7.2.4, is also occurring in passive scalar moments and in mixed
axial velocity-scalar moments in Section 7.2.6.
Finally, the dissertation concludes in Chapter 8, providing an overview of the findings and
presenting future research directions based on the results.

6
2. Governing Equations

The primary focus of this dissertation is on turbulent round jet flows accompanied by a passive
scalar. As previously outlined in Section 1.1, this canonical flow is frequently investigated in
the context of turbulence characterized by its axisymmetric nature, turbulent mixing properties
and self-similar behavior.
In order to describe turbulent round jet flows, we explore the fundamentals of mass and mo-
mentum conservation which are then brought together to derive the NSEs, initially in Cartesian
coordinates, which are well-known for its description of hydrodynamic flows. However, it is
advantageous to express the governing equations in a cylindrical coordinate system because
of the axisymmetry inherent in round jet flows. For this reason, the NSEs are also shown in
cylindrical coordinates though they shall not be derived in detail since the derivation is cum-
bersome. It is common to write the NSEs in a non-dimensionalized form which is also shown
in this chapter. The transformation renders the governing equations independent of specific
unit systems allowing the transferability of findings across experimental setups. This holds
particular relevance in engineering contexts, where insights learned from model experiments
can enhance the design and optimization of real-world applications, such as in aircraft or
vehicles.
Lastly, the convection-diffusion equation is introduced, derived from a concentration balance.
It is a fundamental equation in fluid dynamics that can also describe the evolution of a scalar
quantity advected by the fluid flow. It combines convection and diffusion processes to model
the transport and dispersion of passive scalars such as temperature, concentration of a chemical
species, or contaminant density, without affecting the flow itself, hence the scalar is “passive”.
In turbulent flows, the passive scalar transport equation plays a crucial role in understanding the
mixing and dispersion of scalar quantities, which are often key factors in various engineering
and environmental processes. It is commonly used in computational fluid dynamics (CFD), to
model the propagation of various passive scalar quantities.
Both the NSEs and the convection-diffusion equation serve as the cornerstone for subsequent
analyses within the framework of Lie symmetry analysis and its application in DNS.

2.1. Conservation laws

The ideal description of a gas or liquid in field theory involves assuming that all field quantities
are defined at every point within the domain through local spatial and statistical averaging at a
molecular level. Recognizing that not every infinitesimally small point in a domain is occupied
by gas molecules raises concerns about defining fluid velocity as an averaged field quantity,
potentially leading to errors. In such cases, it becomes more sensible to describe the statistical

7
behavior of thermodynamic systems, such as through the Boltzmann equation. However, the
corresponding equations and statistical descriptions of fluid behavior on a molecular level can
be overly complex for most flows and fluid applications.

To simplify the analysis, it is more convenient to regard the flow domain as a continuum, where
the aforementioned field quantities are spatially averaged in a continuous field, allowing for
the calculation of local gradients of these fields, provided that significant errors do not arise
across all field points.

The justification for this approach is grounded in the mean free path λ, representing the average
length of particle collisions, being sufficiently small in relation to a problem-specific length
scale Lc , expressed as the Knudsen number

λ
Kn =  1, (2.1)
Lc

where the continuum hypothesis applies when Kn≤0.01 (Rapp, 2017).

For the purpose of explaining the use of the continuum hypothesis, a brief calculation (Pope,
2000) is presented, where the mean free path is λ ≈ 6 · 10−8 m for air under atmospheric
conditions. Considering typical length scales in common fluid applications rarely smaller
than Lc = 10−4 m, the calculated Knudsen number Kn = 0.0006 justifies the assumption of
continuum mechanics.

The adoption of continuum mechanics, based on sufficiently low Knudsen numbers, gives rise
to the foundational concept of conservation laws in field theory. These laws assert that certain
quantities within a closed domain remain conserved over time, only changing through external
production, destruction inside the domain or via flux over the boundaries. In addition to the
discussed conservation laws for mass and momentum, other laws exist for energy, angular
momentum, vorticity, and more, particularly for incompressible fluids. The presentation of
basic equations in continuum mechanics are derived from Spurk and Aksel (2010), with a
detailed derivation of these principles available in that reference.

2.1.1. Continuity equation

The conservation of mass states that in a closed system, i.e., a system with no flux over the
boundaries, the mass of that system is constant over time, since no quantity of mass is added
or removed, meaning
D
ZZZ
ρ dV = 0, (2.2)
Dt
V (t)

where
D ∂ ∂
= + Ui (2.3)
Dt ∂t ∂xi
is the material derivative and t, Ui , xi , ρ and V denote the time, flow velocity, spatial coordinate,
density and volume, respectively. Here, the index i describe the three spatial directions in the
Cartesian index notation according to the Einstein summation convention (Einstein, 1916).

8
By applying Reynolds theorem (Reynolds et al., 1903), Equation (2.2) can be reformulated to

D
ZZZ ZZZ ZZ
∂ρ
ρ dV = + ρ ui ni dS, (2.4)
Dt ∂t
V (t) V Sb

where ni is the the surface normal vector which points away from the boundary surface Sb .

Furthermore, Gauss’ theorem (Arfken and Weber, 2005) allows the transformation of the
surface integral on the right hand side of Equation (2.4), which gives
ZZZ  
∂ρ ∂
+ (ρUi ) dV = 0. (2.5)
∂t ∂xi
V

Applying the chain rule to the second term and considering the definition of the material
derivative (2.3) leads to
Dρ
ZZZ  
∂Ui
+ρ dV = 0. (2.6)
Dt ∂xi
V

Both (2.5) and (2.6) have to hold for any integral bound and thus

∂ρ ∂
+ (ρUi ) = 0, (2.7a)
∂t ∂xi
Dρ ∂Ui
+ρ = 0. (2.7b)
Dt ∂xi

Even though both equations are equivalent, we want to focus on (2.7b) to introduce incom-
pressibility. For Mach numbers M a  1, a fluid is considered incompressible, i.e., the density
ρ is constant. Here, M a is defined as
ul
Ma = , (2.8)
c
with ul being the local velocity and c describing the speed of sound (Acheson, 1990). This
reduces Equation (2.7b) to the commonly known continuity equation

∂Ui
= ∇ · U = 0, (2.9)
∂xi

where the flow velocity vector U in a Cartesian coordinate system is given as

U = (U1 , U2 , U3 )T , (2.10)

while the nabla operator ∇ in a Cartesian coordinate system is defined as

 T
∂ ∂ ∂
∇= , , . (2.11)
∂x1 ∂x2 ∂x3

Therefore the continuity equation (2.9) states that the divergence of the velocity field U is
zero for an incompressible fluid which is assumed throughout this work.

9
2.1.2. Momentum equation

The momentum equation, a basic law of mechanics, states that the sum of forces is equal
to the mass times its acceleration. Hereby, the body forces and surface forces, like pressure
and friction forces, are implied. In order to derive the momentum equation, consider the
momentum Ii of a particle of mass m is given as

Ii = mUi , (2.12)

where Ui describes the velocity of the particle in the direction xi . In continuum mechanics, the
momentum of a fluid body of volume V (t) is considered, meaning
ZZZ
Ii = ρUi dV, (2.13)
V (t)

with Ui = Ui (x1 , x2 , x3 , t) being the local velocity of a field point (x1 , x2 , x3 ).

The momentum (2.13) is only changed through surfaces forces ti on the boundary Sb or
through volume forces ki across the volume V , yielding
DUi
ZZZ ZZ ZZZ
ρ dV = ti dS + ρki dV. (2.14)
Dt
V Sb V

The surface force vector ti can be represented as a linear mapping of the surface normal vector
nj through the Cauchy stress tensor τji as

ti = τji nj . (2.15)

With Gauss’s theorem, the surface integral in (2.14) can be reformulated by inserting (2.15),
resulting in
DUi ∂τji
ZZZ  
ρ − − ρki dV = 0. (2.16)
Dt ∂xj
V

Analogously to Section 2.1.1, for the integral above to be zero independent of the integration
boundaries, the integrand has to be zero, obtaining
DUi 1 ∂τji
= + ki = 0, (2.17)
Dt ρ ∂xj

which is recognized as Cauchy’s momentum equation. When combined with the continuity
equation (2.9), it serves as the foundation for the derivation of the NSEs.

2.2. Navier-Stokes equations

The NSEs constitute a fundamental set of partial differential equations (PDEs) that govern
the motion of fluids. Initially formulated by Navier (1827) and later refined with great rigor
by Stokes (1849), these equations are based on several crucial physical principles. Among
these principles are Newton’s law of viscosity, which establishes a relationship between the

10
shear stress in a fluid and the rate of distortion of fluid elements, the conservation of mass
(Section 2.1.1), which asserts that the mass of an isolated system remains constant over time
and Newton’s second law (Section 2.1.2), which describes the relationship between the forces
acting on a system and its resulting motion. The following derivation is based on Batchelor
(2000) and Spurk and Aksel (2010), which can be consulted for a detailed derivation.
To formulate the NSEs of motion, the stress tensor τij is decomposed into an isotropic part
−P δij , where P denotes the pressure, and the deviatoric stress tensor dij . It is convenient to
express dij in terms of the fluid velocity Ui . For this, Batchelor (2000) provides an expression
for the deviatoric stress in isotropic Newtonian fluids as

dij = 2µeij + µ00 ekk δij , (2.18)

where µ, µ00 and δij are the dynamic viscosity, second coefficient of viscosity and the Kronecker
delta, respectively, while the rate-of-strain tensor eij is given by
 
1 ∂Ui ∂Uj
eij = + . (2.19)
2 ∂xj ∂xi

Substituting τij = −P δij + dij into Equation (2.17), while considering Equations (2.18) and
(2.19), yields the general equation of motion

DUi
    
∂ 00 ∂Uk ∂ ∂Ui ∂Uj
ρ = ρki + −P + µ + µ + . (2.20)
Dt ∂xi ∂xk ∂xj ∂xj ∂xi

Neglecting temperature effects, i.e., spatial and temporal dependence of µ and µ00 , Equation
(2.20) simplifies to

DUi ∂P ∂Uk ∂ 2 Ui
ρ = ρki − + (µ + µ00 ) +µ . (2.21)
Dt ∂xi ∂xk ∂xj ∂xj

Assuming incompressibility through the continuity equation (2.9), the divergence terms vanish,
i.e., ∂Uk /∂xk = 0, leading to the incompressible NSEs of motion

DUi ∂P ∗ ∂ 2 Ui
= ki − +ν , (2.22)
Dt ∂xi ∂xj ∂xj

where P ∗ = P /ρ is the pressure rescaled by the constant density and ν = µ/ρ is the kinematic
viscosity. Equation (2.22), accompanied by the continuity equation (2.9), form a set of four
non-linear PDEs.
It is noted that in common practice the Stokes’ hypothesis is applied for compressible flows
which states that the bulk viscosity κ = µ00 + 2µ/3 is zero (Buresti, 2015). This hypothesis
simplifies the description of compressible flows leading to µ00 = −2µ/3 which links both
coefficients of viscosity. A detailed discussion on the validity of this hypothesis can be found in
the aforementioned work.
In practical scenarios, neglecting the body force ki is justified when gravitational forces are
dominated by inertial and shear forces. This yields the simplified form of the incompressible
NSEs
∂Ui ∂Ui ∂P ∂ 2 Ui
+ Uj =− +ν , (2.23)
∂t ∂xj ∂xi ∂xk ∂xk

11
 ∂Ui
= 0, (2.24)
∂xi
where the asterisk notation in P is omitted hereon after for the sake of readability. The NSEs
can also be written in its coordinate system independent form

∂U
+ (U · ∇) U = −∇P + ν∇2 U , (2.25)
∂t
∇ · U = 0, (2.26)
where different coordinate systems can be easily applied. This is advantageous for specific flow
problems e.g. round jet flows where cylindrical coordinates are commonly deployed.

2.2.1. Navier-Stokes equations in cylindrical coordinates

In this subsection, the NSEs in cylindrical coordinates are showcased as this is a useful repre-
sentation for round jet flows. However, the transformation will not be derived in detail here.
The resulting NSEs are taken from White (2016).
The cylindrical coordinate system uses the two-dimensional polar coordinates (r, ϕ) and
combines them with the third Cartesian coordinate z which leads to the three-dimensional
coordinate system (r, ϕ, z). The cylindrical coordinates are related to the Cartesian coordinates
through

x = r cos ϕ,
y = r sin ϕ, (2.27)
z = z,

while the flow velocity vector in a cylindrical coordinate system is given by

U = (Ur , Uϕ , Uz )T . (2.28)

Using tensor calculus, the momentum equations (2.25) can then be written in their coordinate
system independent form as

∂Ur ∂Ur Uϕ ∂Ur Uϕ2 ∂Ur

r: + Ur + − + Uz =
∂t ∂r r ∂ϕ r ∂z
1 ∂ 2 Ur ∂ 2 Ur
   
∂P Ur 1 ∂ ∂Ur 2 ∂Uϕ
− +ν − 2 + r + 2 + − 2 , (2.29)
∂r r r ∂r ∂r r ∂ϕ2 ∂z 2 r ∂ϕ
∂Uϕ ∂Uϕ Uϕ ∂Uϕ Ur Uϕ ∂Uϕ
ϕ: + Ur + + + Uz =
∂t ∂r r ∂ϕ r ∂z
1 ∂ 2 Uϕ ∂ 2 Uϕ
   
1 ∂P Uϕ 1 ∂ ∂Uϕ 2 ∂Ur
− +ν − 2 + r + 2 + + 2 , (2.30)
r ∂ϕ r r ∂r ∂r r ∂ϕ2 ∂z 2 r ∂ϕ
∂Uz ∂Uz Uϕ ∂Uz ∂Uz
z: + Ur + + Uz =
∂t ∂r r ∂ϕ ∂z
1 ∂ 2 Uz ∂ 2 Uz
   
∂P 1 ∂ ∂Uz
− +ν r + 2 + , (2.31)
∂z r ∂r ∂r r ∂ϕ2 ∂z 2

12
while the continuity equation (2.26) in its coordinate system independent form is expressed by

1 ∂ (rUr ) 1 ∂Uϕ ∂Uz

+ + = 0. (2.32)
r ∂r r ∂ϕ ∂z

The advantage of the NSEs in cylindrical form becomes obvious for axisymmetric flows, since
all dependent variables are independent of the coordinate ϕ and Uϕ = 0, which simplifies the
equations significantly.

2.2.2. Non-dimensional Navier-Stokes equations

An essential concept in fluid dynamics, and in physics in general, is non-dimensionalization, i.e.,

the independence of systems from specific unit systems. The expression of physical quantities in
any units should not affect the underlying physics of the system. A system achieves unambiguous
description when both dependent and independent variables are dimensionless. The core
objective of non-dimensionalization is to identify essential dimensionless numbers that describe
the behavior of the physics inherent to the system. In engineering applications for example,
this concept allows for two independently conducted model experiments to be comparable,
even in the face of differences in size, velocity, or fluid viscosity.
The concept of dimensional analysis in fluid mechanics, initially introduced by Reynolds (1883),
plays a central role in the comparison of fluid flow models, irrespective of their size or velocity.
In this context, the NSEs (2.23) and (2.24) are subjected to non-dimensionalization. This
process involves the non-dimensionalization of independent variables, such as space xi and
time t, as well as dependent variables, such as velocity Ui and pressure P , using characteristic
scales. These scales include the characteristic length scale Lc , which represents the physical
domain, e.g., a pipe diameter D or the length of a flat plate L, and the characteristic velocity
scale Uc , which can be the free-stream velocity U∞ for a boundary layer flow or a bulk velocity
Ub for a pipe or jet flow.
The non-dimensionalized quantities, denoted by an asterisk symbol, are expressed as

xi Ui Uc P
x∗i = , Ui∗ = , t∗ = t , P∗ = . (2.33)
Lc Uc Lc Uc2

Note that P has been already rescaled by ρ in Equation (2.22) and is additionally rescaled by
Uc2 for non-dimensionalization here. Using (2.33), the NSEs (2.23) and (2.24) can then be
rewritten to
∂Ui∗ ∗
∗ ∂Ui ∂P ∗ 1 ∂ 2 Ui∗
+ U j ∗ = − ∗ + , (2.34a)
∂t∗ ∂xj ∂xi Re ∂x∗k ∂x∗k
∂Ui∗
= 0, (2.34b)
∂x∗i

with the Reynolds number Re

Uc Lc
Re = , (2.35)
ν
describing the ratio of the inertial forces to friction forces, symbolizing similarity between
systems and experiments.

13
The comparability of analyses, experiments, and numerical studies is crucial in establishing a
framework where the focus lies on the actual physics rather than the unit system. For further
exploration of dimensional analysis in fluid dynamics applications, readers may refer to Simon
et al. (2017) and Sad Chemloul (2020). For readability purposes, the asterisk symbol is omitted
from this point forward. In symbolic form, the momentum balance equations therefore yield
∂U 1 2
+ (U · ∇) U = −∇P + ∇ U. (2.36)
∂t Re

2.3. Convection-diffusion equation

This section is based on Fermigier M. (2017) and Bird et al. (2007) and are referred to for
interested readers. The dynamics of a passive scalar field Θi are governed by the convection-
diffusion equation. In general, this equation describes a physical phenomenon where physical
quantities, e.g., temperature, passive scalar concentration or mass, are transferred inside a
physical system due to diffusion and convection. It can be derived from a concentration balance
on a fixed control volume V . To be more specific, the variation of time of a given species
i within V is considered, where the rate of change of concentration is equal to the net flux
through Sb = ∂V of a species i
ZZZ ZZ
∂ci
dV = − nc,i · nS dS, (2.37)
∂t
V Sb

where ci , nc,i , nS are a concentration of a species i, the concentration flux of a species i and
the unit vector normal to a surface Sb , respectively. The concentration balance is therefore
given for each species i.
To define the convective and diffusive concentration fluxes, the speeds of motion of each species
i needs to be determined. If uc,i is the global speed of a concentration ci then the flux in a
fixed reference frame is nc,i = ci uc,i . The average speed is then
P
i ci uc,i
U= P . (2.38)
i ci

With Fick’s law (Fick, 1855), an expression for diffusion processes is given with

Ji = −Di ∇ci , (2.39)

where Di is the diffusion coefficient of a species i. Since small changes from thermodynamic
equilibrium are considered, the relationship between the fluxes and concentration gradients
are linear.
The total flux of the species contains the convective concentration flux due to the average flow
ci U and the diffusion (2.39)
nc,i = ci U − Di ∇ci . (2.40)
Inserting the equation above into (2.37) and applying the divergence theorem then leads to
the convection-diffusion equation
∂ci
+ U · ∇ci = Di ∇2 ci . (2.41)
∂t

14
For a passive scalar Θ, the diffusion coefficient is proportional to the inverse of the Péclet
number P e and is defined as the product of Re and the Prandtl number P r
ν
Pr = , (2.42)
α
leading to
Uc Lc
Pe = = Re · P r, (2.43)
α
where α describes the thermal diffusivity. The non-dimensionalized passive scalar transport
equation is then written as
∂Θ 1 2
+ (U · ∇) Θ = ∇ Θ. (2.44)
∂t Pe
This equation can describe a concentration of a fluid transported by a velocity field and diffused
by molecular effects without interacting with the flow dynamics. Another way to think of this
equation can be a dye in a flow field that does not interact with it.
This chapter has laid the groundwork by introducing all the fundamental governing equations
relevant to this work. In the succeeding chapter, we discuss these equations within the context
of turbulence.

15
3. Statistical Description of Turbulence

Turbulent flows are inherently chaotic and unsteady. At its core, turbulence is characterized
by irregular and random fluctuations in fluid motion, making it challenging to analyze using
deterministic approaches. Instead, statistical methods offer a means to capture the statistical
properties and overarching patterns within turbulent flows.
This chapter introduces the characterization and analysis of random variables within turbu-
lence. We begin by establishing a solid foundation in probability theory essentials. Central
are fundamental concepts including PDFs, statistical moments and central moments which
collectively are the basis for the statistical description of turbulent flows.
Subsequently, we examine the Gaussian distribution, renowned for its well-defined statistical
properties. The two central moments, skewness and kurtosis, are introduced as metrics for
quantifying asymmetry and tail behavior relative to a Gaussian distribution, thereby offering
descriptive tools for the probabilistic nature of turbulence.
Next, this foundation is applied to turbulence commencing with the Reynolds decomposition,
which separates velocities into mean and fluctuating parts to derive various statistical equations
from the governing equations in Chapter 2. The well-known closure problem of turbulence
arising from the RANS and Reynolds stress transport equations (fluctuation or R-approach)
are introduced. Subsequently, the statistical momentum equations based on the instantaneous
velocities (instantaneous or H-approach), also called MPMEs, are presented which are more
convenient in the context of Lie symmetry analysis.
Finally, we present the statistical equations of passive scalars and velocity-scalar correlations of
the H-approach with the MPSE and the MPVSCEs, respectively.

3.1. Description of random variables

In this section, the basics of probability theory, drawing upon Pope (2000), are briefly reviewed
before applying it to the context of turbulence, particularly for velocities. To describe turbulence,
concepts of probability theory are applied where continuous random variables are described
by PDFs.
A PDF is defined as a function that characterizes the relative likelihood that an event in a
sample space, i.e., a collection of possible outcomes of a random experiment, occurs. As a
measure, a PDF f (V ) is always positive definite and normalized such that the likelihood of any
event occurring is 100% meaning
Z ∞
f (V ) dV = 1. (3.1)
−∞

17
From this, it follows that
f (±∞) = 0. (3.2)

Any statistical moment of order n for a continuous random variable with f ∈ C 1 , where C 1 is
a class with functions having at least a 1st derivative that is continuous in its domain, can be
calculated from a PDF through
Z ∞
Vn = V n f (V ) dV. (3.3)
−∞

Classical moments emerge from this equation, e.g., when n = 1, the first statistical moment
is received also referred to as the mean or expected value. Moments of order n of discrete
random variables are calculated from
N
1 X n
V n = lim Vi , (3.4)
N →∞ N
i=1

where Equations (3.3) and (3.4) are equivalent if one allows Schwartz distributions for f .

Central moments are moments of a PDF about the random variables mean. They quantify
deviations from the mean rather than from zero. For continuous random variables, these
moments are defined as Z ∞
vn = (V − V )n f (V ) dV, (3.5)
−∞

whereas for discrete random variables, the central moments are expressed as

N
1 X
v n = lim (Vi − V )n . (3.6)
N →∞ N
i=1

Classical
p central moments, derived from v , include the variance v and the standard deviation
n 2

σ = v2.

Standardizing central moments of higher order by the standard deviation allows these moments
to become scale invariant, facilitating the comparison of various PDFs. Two commonly examined
quantities in this context are skewness, given by

v3
S= , (3.7)
σ3
which describes the asymmetry of a PDF. A negative value indicates a tail on the left side of
the distribution, while a positive value suggests a tail on the right side. A skewness of zero
indicates, but not necessitates, symmetry.

Additionally, kurtosis is defined as

v4
K= , (3.8)
σ4
and it serves as a measure of the distribution’s tail weight. A higher kurtosis implies larger
tails, and conversely, a lower kurtosis indicates smaller tails in the distribution.

18
3.1.1. Gaussian distribution

Building upon these foundations, we now possess the tools to comprehend the Gaussian distri-
bution or normal distribution. It is characterized by a bell-shaped curve, which is completely
defined by its mean V and standard deviation σ. The PDF of a Gaussian distribution is given by
 2
1 −1 V −V
f (V ) = √ e 2 σ
. (3.9)
σ 2π
This distribution serves as a benchmark in statistical analysis. Its well-defined shape and
characteristics make it a reference point for understanding the behavior of other distributions.
Skewness and kurtosis values of alternative distributions are often compared to those of the
Gaussian distribution, helping researchers interpret the symmetry and tail behavior of different
datasets.
The Gaussian distribution is characterized by special properties. All higher standard moments
v n of a Gaussian distribution are determined by the variance (Papoulis and Pillai, 2002) with
(
vn 0 if n odd
= (3.10)
σ n
(n − 1)!! if n even,

where ( · )!! is the double factorial (Arfken and Weber, 2005) defined as

k · (k − 2) · (k − 4) · . . . · 2 for even k > 0,



k!! = k · (k − 2) · (k − 4) · . . . · 1 for odd k > 0, (3.11)
for k ∈ {−1, 0}.

1


In the case of a Gaussian distribution, the skewness (n = 3) is therefore S = 0, indicating

symmetry while the kurtosis (n = 4) of the Gaussian distribution is K = 3 indicating a standard
tail behavior.

3.1.2. Statistics in the context of turbulence

With the statistical framework for random variables established, we can now apply it to
turbulence.
Reynolds (1895) suggested that the velocity U and pressure P can be described as the sum of
the time-average and the fluctuation

U = U + u, P = P + p, (3.12)

also known as the Reynolds decomposition, where u = 0 follows by time-averaging (3.12).

This decomposition allows us to relate instantaneous moments, e.g, Ui Uj , to its fluctuating
moments, e.g., ui uj . Multiplying the velocity in Equation (3.12) by itself and averaging yields

Ui Uj = Ui Uj + ui Uj + Ui uj + ui uj = Ui Uj + ui Uj + Ui uj + ui uj = Ui Uj + ui uj . (3.13)

This approach can be extended to higher moments by multiplying (3.12) by itself repeatedly,
although we omit the details here.

19
We define moments based on instantaneous velocities as

(0)
Hij = Ui Uj , (3.14)

while moments based on fluctuation velocities are denoted as

(0)
Rij = ui uj , (3.15)

(0)
where the superscript ( · )(0) indicates the correlation at one point. Rij is also known as the
Reynolds stress tensor which is symmetric and equivalent to the variance of the velocities. In
this dissertation only R-moments at one point are considered. Therefore, the superscript (0)
for R-moments is omitted unless explicitly stated otherwise. The standard deviation,p often
referred to as the root-mean square (rms) velocity or turbulent intensity, is defined as Rij ,
providing a measure of the intensity of turbulent fluctuations in the flow.

The introduction of Hij and Rij distinguishes two distinct approaches: the omission of the
Reynolds decomposition is called the H-approach, while the conventional method is denoted
as the R-approach.

3.2. Reynolds-averaged Navier-Stokes equations

With the theoretical framework established so far, the RANS equations can be derived. This
classical approach based on fluctuation velocities is referred to as the R-approach as mentioned
above. By implementing the decomposition of U and P in (3.12) into the NSEs (2.23) and
(2.24), we receive
∂Ui ∂Ui ∂P ∂ 2 Ui ∂Rij
+ Uj =− +ν − , (3.16)
∂t ∂xj ∂xi ∂xk ∂xk ∂xj

and
∂Ui
= 0, (3.17)
∂xi
respectively. Here, the introduction of the RST introduces 6 unknowns. However, the number
of unknowns exceeds the available equations, resulting in a closure problem also known as the
turbulence closure problem. To address this, additional equations are necessary to describe
the RST, called the Reynolds stress transport equations which will be derived in the following.

An equivalent equation for a passive scalar can be derived from Equation (2.44) by decomposing
the passive scalar with
Θ = Θ + θ, (3.18)

which is not shown here, as the focus of this dissertation lies on the H approach. It is important
to note that all discussed issues in this and the following section also apply to the R-equations
for passive scalars.

20
3.3. Reynolds stress transport equations

In order to derive the Reynolds stress transport equations (Chou, 1945), which describes
the unknown Reynolds stress tensor Rij , the equation based on the fluctuation velocities are
required. The continuity equation of fluctuation velocities can be received by subtracting the
continuity equation (2.24) and the averaged continuity equation (3.17), receiving

∂ui
= 0, (3.19)
∂xi

while the momentum balance equations of the fluctuation velocities is analogously determined
by subtracting momentum balance equations (2.23) and the averaged momentum balance
equations (3.16)

∂ui ∂ui ∂Ui ∂p ∂ 2 ui ∂Rij ∂ui uj

Fi = + Uj + uj + −ν − + = 0. (3.20)
∂t ∂xj ∂xj ∂xi ∂xk ∂xk ∂xj ∂xj

By cross-multiplying Equation (3.20) with the velocity fluctuation u

Fi uj + ui Fj = 0, (3.21)

the Reynolds stress transport equations are

   
∂Rij ∂Rij ∂ ∂ui ∂uj ∂Uj ∂Ui
+ Uk =− (δjk pui + δik puj ) + p + − Rik + Rjk
∂t ∂x
| {z } | {z k} ∂x k ∂xj ∂xi ∂xk ∂xk
| {z } | {z } | {z }
Lij Cij Πdij Πsij Pij

∂ ∂ui ∂uj ∂ 2 Rij

− Rijk − 2ν +ν , (3.22)
∂xk ∂xk ∂xk ∂xk xk
| {z } | {z } | {z }
Tij ij Vij

where the terms Lij , Cij , Πdij , Πsij , Pij , Tij , ij , Vij are the local temporal term, convection,
pressure diffusion, pressure strain, production, turbulent diffusion, dissipation and viscous
diffusion, respectively, and
Rijk = ui uj uk . (3.23)
Further, it is noted that the pressure strain and pressure diffusion term can be simplified to the
velocity-pressure gradient correlation tensor
 
∂p ∂p
Πij = Πdij − Πsij = − ui + uj . (3.24)
∂xj ∂xi

In the Reynolds stress transport equations (3.22), all terms except for Lij , Cij , Pij and Vij are
unclosed which further complicates the closure problem. This issue continues for higher order
R-equations. In practice, these unclosed terms are typically closed using empirical models,
which replace the unclosed terms. However, these models are often specialized and may only
be applicable to certain types of flows. In the next section, we will look at the equations based
on the instantaneous approach or H-approach which maintain the closure problem but gain
properties that are advantageous for Lie symmetry analysis.

21
3.4. Multi-point moment equations

The H-MPMEs are introduced in this section. Due to the complexity of Reynolds stress transport
equations (3.22) based on the R-approach, Oberlack and Rosteck (2010) proposed a novel
approach, suggesting the omission of the Reynolds decomposition (3.12), i.e, they advocate
for handling the moments based on instantaneous velocities (H-approach). The ensuing
definitions and equations are presented in accordance with Oberlack and Rosteck (2010).
Readers interested in further elaborations are encouraged to consult the original work for
comprehensive details.
The H-MPMEs are based on the moments of instantaneous velocities

(3.25)



Hi{n} = Hi(1) i(2) ...i(n) = Ui(1) x(1) , t · . . . · Ui(n) x(n) , t ,

where Ui(n) x(n) , t is the nth velocity in the direction i(n) = 1, 2, 3 at the nth point x(n) in

space. In the following, we omit the inclusion of time t to shorten the notation although it is
implied unless stated otherwise. The correlation at one point, i.e., x = x(n) , is denoted by
(0)
Hi{n} . Further, the connection to the mean velocity is given by Hi{1} = Ui(1) x(1) = Ui , while

the relation between the H-moment of second order with the Reynolds stress tensor is
(0) (0)
Hi{2} = Hij = Hi Hj + Rij = Ui Uj + Rij . (3.26)

Further, the multi-point velocity-pressure correlation tensor is defined as

Iin−1 [l] = Ui(1) x(1) · . . . · Ui(l−1) x(l−1) P x(l) Ui(l+1) x(l+1) · . . . · Ui(n) x(n)
n

Y
(3.27)


= P x(l) Ui(j) x(j) .
j=1,j6=l

With the aforementioned definitions, the H-MPMEs are then derived with the momentum
equations in (2.34) rearranged to

∂Ui (x, t) ∂Ui (x, t) ∂P (x, t) 1 ∂ 2 Ui (x, t)

Mi (x) = + Uk (x, t) + − = 0, (3.28)
∂t ∂xk ∂xi Re ∂xk ∂xk
and
n
X n

Y
(3.29)


Ni{n} = Mi(a) x(a) Ui(b) x(b) = 0,
a=1 b=1,b6=a

receiving the H-MPMEs

 
n
" #
∂Hi{n} X ∂Hi{n+1} i(n) 7→k(l)  x(n) 7→ x(l) ∂Ii{n−1} [l] 2
1 ∂ Hi{n}
Ni{n} = + + − = 0,
∂t ∂xk(l) ∂xi(l) Re ∂xk(l) ∂xk(l)
l=1
(3.30)
where

(3.31)
 



Hi{n+1} i(n) 7→k(l)  x(n) 7→ x(l) = Ui(1) x(1) · . . . · Uk(l) x(l) Ui(n+1) x(n+1) .

The H-MPMEs are considerably more convenient than the R-equations, especially for the
symmetry analysis conducted in this dissertation. They present an infinite linear system of

22
equations with coupling exclusively between orders n and n + 1. Their simplicity enables
generalization to arbitrary moments, a feat not possible with the R-equations. However, despite
this advantage, the interpretability of H-moments is more complicated compared to their R
counterparts. In Lie symmetry analysis, the focus shifts to H-moments, which can nonetheless
be transformed into R-moments due to their mathematical equivalence.

3.5. Multi-point scalar moment equation

Similarly to the H-MPMEs, the MPSE can be derived for Equation (2.44) gaining the advantages
equivalent to the H-MPMEs in (3.30).
The starting point is the rearranged passive scalar equation (2.44) at one point x ∈ R3

∂Θ(x, t) ∂Θ(x, t) 1 ∂ 2 Θ(x, t)

S(x) = + Uk (x, t) − = 0. (3.32)
∂t ∂xk P e ∂xk ∂xk
The two-point passive scalar equation can then be received through statistically averaging the
cross-multiplication of (3.32) and the passive scalar at different points x and y ∈ R3

S2 = S(x)Θ(y) + S(y)Θ(x) = 0, (3.33)

leading to
∂HΘΘ (x, y) ∂HΘΘk (x, y, x) ∂HΘΘk (y, x, y)
S2 = + +
∂t ∂xk ∂y
 2 2
 k
1 ∂ HΘΘ (x, y) ∂ HΘΘ (x, y)
− + = 0, (3.34)
Pe ∂xk ∂xk ∂yk ∂yk
where
HΘΘ (x, y) = Θ(x)Θ(y) (3.35)
and
HΘΘk (x, y, x) = Θ(x)Θ(y)Uk (x). (3.36)
This can be extended up to an arbitrary order m by introducing
m
X m

Y
(3.37)


Sm = S x(a) Θ x(b) = 0,
a=1 b=1,b6=a

where Θ x(i) is a passive scalar at different points x(i) ∈ R3 . Plugging Equations (3.32) into

(3.37) then yields

m
"  
∂HΘ{m} X ∂HΘ{m+1} [Θ(m+1) 7→k(l) ] x(m+1) 7→ x(l)
Sm = +
∂t ∂xk(l)
l=1
#
2
1 ∂ ∂HΘ{m}
− = 0 for m = 1, . . . , ∞, (3.38)
P e ∂xk(l) ∂xk(l)

where
HΘ{m} = Θ(x(1) )Θ(x(2) ) · . . . · Θ(x(m) ) (3.39)

23
and
 
HΘ{m+1} [Θ(l) 7→k(l) ] x(l) 7→ x(p) =
Θ(x(1) ) · . . . · Θ(x(l−1) )Uk(l) (x(p) )Θ(x(l+1) ) · . . . · Θ(x(m+1) ). (3.40)

3.6. Multi-point velocity-scalar correlation equations

In this section, the derivation of the MPVSCEs is reviewed, which forms a correlation between
the velocity and passive scalar moments. This equation also generalizes the MPMEs and the
MPSE to one equation as demonstrated in the following.
Using the passive scalar equation (3.32) and the moment equations (3.28), the two-point
velocity-scalar correlation equations can be derived with

TiΘ = Mi (x)Θ(y) + Ui (x)S(y), (3.41)

which, by inserting the two aforementioned equations, can be extended to

∂HiΘ (x, y) ∂HiΘk (x, y, x) ∂HΘik (y, x, y)

TiΘ = + + (3.42)
∂t ∂xk ∂yk
2 2
 
∂P (x)Θ(y) 1 ∂ HiΘ (x, y) 1 ∂ HiΘ (x, y)
+ − + = 0.
∂xi Re ∂xk ∂xk P e ∂yk ∂yk

Similarly to Equation (3.37), the two-point velocity-scalar correlation equations can be ex-
tended to the MPVSCEs for n velocities at n different points and m passive scalars at m different
points

n
X n
Y n+m
Y
Ti{n} Θ{m} = Mi(a) (x(a) ) Ui(c) (x(c) ) Θ(x(d) )
a=1 c=1,c6=a d=n+1
n+m
X n+m
Y n
Y
+ S(x(a) ) Θ(x(c) ) Ui(d) (x(d) ) = 0, (3.43)
a=n+1 c=n+1,c6=d d=1

for n, m = 1, . . . , ∞. By introducing the moment equations (3.28) and the passive scalar
equation (3.32) into Equation (3.43), we receive the MPVSCEs

∂Hi{n} Θ{m}
Ti{n} Θ{m} = +
∂t  
n
!
X ∂Hi{n+1} Θ{m} i(n+m+1) 7→k(l)  x(n+m+1) 7→ x(l) ∂Ii{n−1} Θ{m} [l] 2
1 ∂ Hi{n} Θ{m}
+ −
∂xk(l) ∂xi(l) Re ∂xk(l) ∂xk(l)
l=1
 
n+m
!
X ∂Hi{n+1} Θ{m} i(n+m+1) 7→k(l) x(n+m+1) 7→ x(l)
 2
1 ∂ Hi{n} Θ{m}
+ − = 0, (3.44)
∂xk(l) P e ∂xk(l) ∂xk(l)
l=n+1

where
Hi{n} Θ{m} = Ui(1) (x(1) ) · . . . · Ui(n) (x(n) )Θ(x(n+1) ) · . . . · Θ(x(n+m) )

24
and the pressure correlation term is defined as

Ii{n−1} Θ{m} [l] = Ui(1) (x(1) ) · . . . · P (x(l) ) · . . . · Ui(n) (x(n) )Θ(x(n+1) ) · . . . · Θ(x(n+m) ). (3.45)

The MPVSCEs (3.44) are therefore a generalization of the MPMEs and the MPSE. The MPMEs
can be received for m = 0, while the MPSE arise for n = 0.
In the next chapter, the Lie symmetry theory is introduced, which is the main method used in
the present dissertation in order to analyze the governing equations that have been described
in this chapter.

25
4. Lie Symmetry Theory

Most people have an intuitive understanding for symmetries as unique properties inherent in
certain geometric objects. An example is rotational symmetry, or a rotational transformation
that leaves an object, such as a cylinder or a sphere, unchanged. This concept extends beyond
geometric shapes to differential equations, and proves to be valuable not only as a mathematical
tool, but also as a conceptual framework to enhance our understanding of physics and, by
extension, turbulence.
In this chapter, we begin by laying the groundwork with a discussion on symmetries of algebraic
equations. Building upon this foundation, we extend this to symmetries of differential equations,
introducing the jet notation, in order to apply symmetry theory to the derivatives of the
dependent variables and hence to the H-MPMEs in Chapter 6 as well as to the MPSE and
MPVSCEs in Chapter 7.
For a more in-depth exploration of this topic, interested readers are encouraged to read the
works of Bluman et al. (2010) and Bluman (2002) on which this chapter is based on. Note that
the following explanation is quite elementary and kept very brief. A deeper explanation requires
considerable understanding of differential geometry for which we refer to Olver (2000).

4.1. Symmetries of algebraic equations

In this section, we first focus on symmetries of algebraic equations before diving into symmetries
of differential equations in the following section.
A symmetry, defined as a variable transformation of a set of variables x, in its global form

T : x∗ = φ(x; a) (4.1)

with an arbitrary constant a ∈ R, referred to as the group parameter, transforms an equation

F (x) = 0 (4.2)

such that
F (x) = 0 ⇐⇒ F (x∗ ) = 0, (4.3)
i.e, such that equation F is invariant under the variable transformation in Equation (4.1).
Within this dissertation, it is assumed that symmetries possess group properties and, further-
more, qualify as simply-connected Lie groups as defined in Hawkins (2000). With this, φ can
be expressed and reparameterized in a manner, where a = 0 serves as the identity element

x∗ = φ(x; a = 0) = x. (4.4)

27
Hereby, the composition operation is

φ(φ(x; a1 ); a2 ) = φ(x; a1 + a2 ). (4.5)

Further, the global form (4.1) can be expanded as a Taylor series around a = 0, leading to

∂φi ∂φi
x∗i = φi (x; a = 0) + a + O(a2 ) = xi + a + O(a2 ), (4.6)
∂a a=0 ∂a a=0

where O is the order of magnitude.

In Lie group theory, Lie’s first theorem is essential for the following step. According to this
theorem, the linear parameter a in Equation (4.6) plays an important role by uniquely deter-
mining the complete action of the transformation (4.1). This insight inherently leads to the
infinitesimal generator
∂
X = ξi , (4.7)
∂xi
where the infinitesimals ξi are given by

∂φi
ξi = , (4.8)
∂a a=0

which defines the invariance of (4.2) under (4.1). With Lie’s theorem, Equation (4.8) is
sufficient to recover (4.1) which is based on various group properties such as Equations (4.4)
and (4.5). With the definition of an infinitesimal generator (4.7), the equation F in (4.2) is
also considered invariant under a symmetry with an infinitesimal generator X, if and only if

XF |F =0 = 0 (4.9)

in a neighborhood of any solution F = 0 of maximal rank (Olver, 2000), which is an alternative

form to the global form of a symmetry in Equation (4.1).

4.2. Symmetries of differential equations

The extension of the symmetry theory from algebraic equations to differential equations is
straightforward. In the following, the PDF system

F (x, y, y (1) , . . . , y (p) ) = 0 (4.10)

is considered, where y are the dependent variables and y (p) refers to the pth derivative of y.
The jet notation is introduced through

∂yi
yi,j =
∂xj
∂yi (4.11)
yi,jk =
∂xj ∂xk
..
.

28
where the derivatives in jet notation should be treated as variables not related to the dependent
variables y or the independent variables x. The relation is only important in view of how the
transformation from (x, y) 7→ (x∗ , y ∗ ) extends to the jet variables, which we will show in the
following. First, we consider the global form

T : x∗ = φ(x, y; a), y ∗ = ψ(x, y; a) (4.12)

for which, analogously to above in Equation (4.6), the global form can be written in an
infinitesimal form through a Taylor expansion in a, receiving

∂φi ∂φi
x∗i = φi (x, y; a = 0) + a + O(a2 ) = xi + a + O(a2 ), (4.13)
∂a a=0 ∂a a=0

∂ψi ∂ψi
yi∗ = ψi (x, y; a = 0) + a + O(a2 ) = yi + a + O(a2 ). (4.14)
∂a a=0 ∂a a=0

The infinitesimal generator can now be extended to

∂ ∂
X = ξi + ηk , (4.15)
∂xi ∂yk

where
∂ψk
ηk = (4.16)
∂a a=0

denotes the infinitesimals of the dependent variables. To include the derivatives, the infinitesi-
mal generator is further extended to a so-called prolonged operator that emerges through the
chain rule of differentiation. The nth prolonged operator is

∂ ∂ ∂ ∂
X (n) = ξi + ηk +ηk;j1 + ηk;j1 j2 + · · ·, (4.17)
∂xi ∂yk ∂yk,j1 ∂yk,j1 j2
| {z }
n terms

where ; denotes that jn does not indicate the derivatives in the jet notation and
n−1
Dηk;j1 ...jn−1 X Dξjm
ηk;j1 ...jn = − ηk;j1 ...jn−1 jm , (4.18)
Dxn Dxn
m=1

with the total differentiation operator is given as

D ∂ ∂ ∂
= + yk,i + yk,ij ··· . (4.19)
Dxi ∂xi ∂yk ∂yk,j

To construct invariant solutions from symmetries of the system (4.10), we may consider the
sought solution
y = f (x) (4.20)
which satisfies the invariant surface condition

X [y − f (x)]|y=f (x) = 0. (4.21)

29
Equation (4.21) results in a hyperbolic PDE system, which may be solved using the method of
characteristics, i.e., the invariant surface condition (4.21) can be rewritten as
dx1 dx2 dxi dy1 dy2 dyk
= = ··· = = = = ··· = . (4.22)
ξ1 (x, y) ξ2 (x, y) ξi (x, y) η1 (x, y) η2 (x, y) ηk (x, y)
Solution to system (4.22) is then determined in terms of the reduced set of invariants x
e(x, y)
and ue (x, y), where x
e denote the new independent variables and ye denote the new dependent
variables. The functional form ye(e x) must then be substituted into the governing equation
(4.10) to find the reduced set of equations.
This theory also introduces an algorithmic approach to calculate symmetries of any given
equation F . This can be done by assuming a general infinitesimal generator X and calculating
the prolongation according to Equation (4.17) as needed for F to cover all the derivatives. After
inserting into Equation (4.9), the resulting system of PDEs can be solved for the infinitesimals
ξi and ηk . These systems are generally tedious to solve, but with the advent of computer
algebra systems, this effort has been reduced significantly. It is noted that no ansatz is generally
required in order to find all simply-connected Lie-point symmetries of F but rely solely on the
algorithmic approach and the local existence of solutions to F (Olver, 2000).

4.2.1. Lie symmetries of a one-dimensional heat equation

In the following, the calculation of the infinitesimals, and hence the Lie symmetries, shall be
exemplified with the one-dimensional heat equation
∂u ∂2u
= . (4.23)
∂t ∂x2
In order to stick with the convention introduced above the variables shall be renamed to
u → y1 , x → x1 , t → x2 , (4.24)
which leads to the PDE
∂y1 ∂ 2 y1
= , (4.25)
∂x2 ∂x21
and is written in jet notation as
F = y1,2 − y1,11 = 0. (4.26)
Since this equation is a PDE of second order, the prolonged operator (4.17) is
∂ ∂ ∂
X (2) =ξ1 (x1 , x2 , y1 ) + ξ2 (x1 , x2 , y1 ) + η1 (x1 , x2 , y1 )
∂x1 ∂x2 ∂y1
∂ ∂
+ η1;2 + η1;11 , (4.27)
∂y1,2 ∂y1,11
where non-relevant terms, such as η1;1 , have already been neglected. The remaining terms can
then be determined with Equation (4.18), receiving
Dη1 Dξ1 Dξ2
η1;1 = − y1,1 − y1,2 ,
Dx1 Dx1 Dx1
Dη1 Dξ1 Dξ2
η1;2 = − y1,1 − y1,2 ,
Dx2 Dx2 Dx2
Dη1;1 Dξ1 Dξ2
η1;11 = − y1,11 − y1,12 , (4.28)
Dx1 Dx1 Dx1

30
where the total differentiation operator (4.19) is

D ∂ ∂ ∂ ∂
= + y1,1 + y1,11 + y1,12 ,
Dx1 ∂x1 ∂y1 ∂y1,1 ∂y1,2
D ∂ ∂ ∂ ∂
= + y1,2 + y1,21 + y1,22 . (4.29)
Dx2 ∂x2 ∂y1 ∂y1,1 ∂y1,2

At this point it becomes obvious that employing computer algebra systems significantly simplifies
the calculation of the prolonged operator. Finally, after inserting (4.28) into the prolonged
operator (4.27) under the consideration of the total differentiation operator (4.29), we can
apply

X (2) F = X (2) (y1,2 − y1,11 ) = 0, (4.30)

which leads to

∂ 2 η1 ∂ 2 η1 ∂ 2 ξ1 ∂ 2 ξ2
   
∂η1 ∂ξ1 ∂ξ1 ∂ξ2
− + y1,1 −2 − + + y1,2 2 − + 2
∂x2 ∂x21 ∂x1 y1 ∂x2 ∂x21 ∂x1 ∂x2 x1
2 2 2
   
2 ∂η1 ∂ ξ1 ∂ξ1 ∂ ξ2 3 ∂ ξ1
+ y1,1 −2 + y1,1 y1,2 2 +2 + y1,1
∂y1 ∂x1 y1 ∂y1 ∂x1 y1 ∂y12
∂ 2 ξ2 ∂ξ2 ∂ξ2
2
+ y1,1 y1,2 2 + y1,12 2 + y1,12 y1,1 2 = 0, (4.31)
∂y1 ∂x1 ∂y1

where the jet variables have been factored out. To solve this equation, the monomials have
to be zero since the jet variables are independent, i.e., the coefficients of the monomials in
Equation (4.31) give a system of nine PDEs which can be solved to receive the infinitesimals

a5 x21
   
a6
η1 = y1 a3 − + x2 − x1 + f2 (x1 , x2 ) (4.32)
2 2 2
ξ1 = a1 + a4 x1 + a5 x1 x2 + a6 x2 (4.33)
ξ2 = a2 + 2a4 x2 + a5 x22 (4.34)

with

∂ 2 f2 (x1 , x2 ) ∂ 2 f2 (x1 , x2 )
= (4.35)
∂x22 ∂x21

31
and hence the symmetries

∂
X1 = :
∂x1
x̃1 = x1 + ; x̃2 = x2 ; ỹ1 = y1 , (4.36)
∂
X2 = :
∂x2
x̃1 = x1 ; x̃2 = x2 + ; ỹ1 = y1 , (4.37)
∂
X3 =y1 :
∂y1
x̃1 = x1 ; x̃2 = x2 ; ỹ1 = e y1 , (4.38)
∂ ∂
X4 =x1 + 2x2 :
∂x1 ∂x2
x̃1 = e x1 ; x̃2 = e x2 ; ỹ1 = y1 , (4.39)
 2 
∂ 2 ∂ x1 x2 ∂
X5 =x1 x2 + x2 − + y1 :
∂x1 ∂x2 4 2 ∂y1
2
x1
x1 x2 y1 −1
x̃1 = ; x̃2 = ; ỹ1 = √ e 4 1−x2 , (4.40)
1 + x2 1 − x2 1 − x2
∂ x1 ∂
X6 =x2 − y1 :
∂x1 2 ∂y1
 
x̃1 = x1 + x2 ; x̃2 = x2 ; ỹ1 = y1 e− 2 (x1 + 2 x2 ) . (4.41)

While the approach relying on infinitesimal generators results in a linear system of PDEs for
computing the invariants, other methods, e.g., the moving frame method, necessitates solving
a system of non-linear algebraic equations. This system is instead fully coupled and introduces
complexity that varies depending on the problem. Given that the infinitesimal-based method
enables a decoupled solution of the system of linear PDEs, it was determined to be more
effective for the current objectives.
With the newly gained knowledge, we are now able to utilize Lie symmetry theory in Chapter 6
and Chapter 7. But first, the subsequent chapter introduces round jet flows and the two DNS
that have been performed to study this flow.

32
5. Direct Numerical Simulations of a Turbulent
Round Jet Flow at two Reynolds Numbers

In this chapter, two DNS of spatially evolving turbulent round jet flows at different Reynolds
numbers conducted in this dissertation are discussed. This chapter is heavily based on the
publication Nguyen and Oberlack (2024a), where the DNS data of the turbulent round jet flow
at Re = 3500 is thoroughly discussed and Nguyen and Oberlack (2024b), which compares the
data to a turbulent round jet flow at Re = 7000.
The chapter begins by providing background on prior research in the field, including experi-
ments of turbulent round jet flow followed by preceding works on DNS of such flows. Following
this, the specifics of the present DNS are discussed, which involves the computational domain,
notably consisting of two parts and therefore two simultaneous running simulations. Finally,
the results obtained from the two DNS are presented and compared with existing literature.
These results comprise statistical moments up to the third order as well as PDFs of the axial
velocity. Additional DNS datasets of moments up to the tenth order and passive scalar moments
are presented and discussed in the subsequent chapters.
Round jet flows describe fluid flows ejected from a circular inlet into a large area. A sketch
of a round jet is depicted in Figure 5.1. At a sufficiently high Reynolds number, this flow

Figure 5.1.: Schematic view of a round jet flow. Fluid is blown through a nozzle with diameter D. The mean axial
velocity of the jet is denoted by U z while the mean axial centerline velocity is denoted by U z,c .

becomes turbulent, an important property that goes beyond academic curiosity. As mentioned
in Section 1.1, turbulent jets play a critical role in jet propulsion systems, where understanding
and optimizing jet dynamics is essential to improving engine performance and fuel efficiency.
In addition, turbulent round jets serve as test beds for validating turbulence models and

33
simulation techniques, providing benchmarks for assessing the accuracy and reliability of
numerical predictions through its unique properties. In particular, turbulent round jet flows
exhibit self-similarity in the far field, i.e., radial profiles of different statistical quantities can
be scaled to converge to a single similarity function.
Universal self-similarity for turbulent round jets has long been discussed, where this state
should be reached asymptotically and independent of initial conditions (Townsend, 1976).
Yet, theoretical analysis show that these self-similar profiles can depend on initial conditions
(George, 1989). Experimental results have indicated that, although a single jet is self-similar,
the profiles differ in various experiments (Hussein et al., 1994; Panchapakesan and Lumley,
1993a; Wygnanski and Fiedler, 1969). First assumptions have been inaccurate measurements
due to disturbances by temperature or uncontrollable turbulence in the ambient region which
could be regulated by numerical simulations. However, since turbulent round jets grow in the
axial direction, numerical simulations are computationally expensive.
Therefore, until the 1980s, the detailed statistical study of turbulent jet flows relied primarily
on experimental methods. Wygnanski and Fiedler (1969) conducted their experiment at a
Reynolds number of Re = 100,000, where the Reynolds number is defined as
Ub D
Re = , (5.1)
ν
with Ub , D, and ν being the inlet bulk velocity, inlet diameter and kinematic viscosity, re-
spectively. They took velocity measurements up to a normalized axial distance of z/D = 100
for first to third order statistics highlighting a substantial disparity in the distance where
self-similarity is attained ranging from z/D = 20 for the first order moment up to z/D = 70 for
the second order moment. This is explained with the energy, where self-similarity is reached
in steps. First, the mean velocity reaches self-similarity which amounts to a certain production
of fluctuation which in turn allows the reach of equilibrium for the second order moments.
Considerably later, Panchapakesan and Lumley (1993a) conducted an experiment at Re =
11,000 taking velocity measurements in range of z/D = 30 − 150 showing self-similarity up
to fourth order moments. In parallel, Hussein et al. (1994) conducted their experiment at
Re = 95,500 taking velocity measurements between z/D = 50 − 122 while the setup was
based on Wygnanski and Fiedler (1969). Both groups were able to describe the turbulent
kinetic energy (TKE) balance but made different assumptions for the dissipation and pressure
diffusion to close that balance. Consequently, the results of both show disagreement due to
these assumptions.
Xu and Antonia (2002) conducted an experiment measuring a jet exiting from a smooth
contraction nozzle and from a pipe with a pipe flow profile at Re = 86,000, taking measurements
in the range z/D = 1 − 75. They show that a contraction jet reaches a state of self-preservation
earlier than a fully developed pipe jet. This contrasts the experiments of Ferdman et al. (2000),
where they state that a pipe jet reaches self-preservation earlier. However, both agree that
for a fully developed pipe profile the far-field decay rate is smaller than in a top-hat velocity
distribution attributing it to a lack of potential core and surrounding mixing layer, which lead
to a slower centerline mixing. Thus, the turbulence intensities on the centerline are smaller
than those of a top-hat profile according to Ferdman et al. (2000).
Recently, Darisse et al. (2015) collected results on a slightly heated round jet at Re = 140,000
where he also generated data on scalar quantities which were gathered at z/D = 30. They

34
measured pure and mixed moments up to the third order, which allowed not only the deter-
mination of the TKE balance, where, again, the pressure diffusion term was modeled and
dissipation is found by closing the balance, but also the passive scalar transport balance, in
which the dissipation is also calculated by closing the balance.
With the advent of supercomputers in the 1980s, DNS for simulating the motion of fluid flow
gained more popularity. Unlike RANS or large-eddy simulation (LES), which rely on turbulence
models derived from experimental or DNS data to approximate the effects of small-scale
turbulence, DNS aims to provide a complete description of the flow field by resolving all
relevant scales of motion. While turbulence models can mitigate the computational cost, they
often do so at the expense of accuracy. Therefore, DNS offers a high-fidelity representation of
turbulent flows without the need for turbulence models, making it particularly attractive for
fundamental studies where accuracy is of high importance.
However, DNS faces several challenges and limitations. Achieving adequate grid resolution to
resolve all relevant scales of motion can be particularly challenging for turbulent flows with
large Reynolds numbers. Additionally, DNS requires small time steps to accurately capture fast
time scales, leading to increased computational demands.
Despite these challenges, DNS has applications across various fields, including fundamental
research, engineering design, and environmental studies, providing insights into fundamental
aspects of fluid dynamics. The ability of DNS to directly simulate turbulent flows at high
fidelity makes it a valuable tool for investigating complex flow phenomena and improving our
understanding of fluid dynamics while simultaneously providing data to optimize turbulence
models.
One of the first DNS of a spatially developing turbulent round jet was by Boersma et al. (1998).
They studied the effect of inflow conditions on the self-similarity scaling up to a box length
of z/D = 45 at Re = 2400, showing that self-similarity are dependent on initial conditions.
However, they were unable to resolve the far-field due to limited computational resources so
the fluctuations have not yet fully reached self-similarity.
Taub et al. (2013) conducted a DNS at Re = 2000 at the same box length where he extracted
statistics up to the third order, including the TKE terms and the terms of the Reynolds stress
transport equations terms, directly. The data is extensively compared to previous studies.
They also conclude that there are inconsistencies in the dissipation profiles due to different
approximations in previous experimental studies, which may also be due to the Reynolds
number dependence of the dissipation (Bogey and Bailly, 2006).
Later, Shin et al. (2017) found self-similarity in a DNS at Re = 7290 after z/D = 15 for the
first and after z/D = 25 for the second moment in a box length of z/D = 60. Only in their
study, the radial profiles of the moments of second order are shown at various distances z but
the collapse into a single similarity function is rather unsatisfactory since the statistics have
been only taken over 80D/Ub units. They show that moments up to the fourth order and the
one-point PDF of mass-weighted stream-age, the product of the jet fluid residence time and
the jet fluid mass fraction, are self-similar.
In contrast, the present study presents high-quality statistics from two new spatially evolving
turbulent round jet DNS at a comparably large Reynolds numbers with a long box. The first
case, conducted at Re = 3500, involves statistics collected over 75,000D/Ub time units where
the velocity profile of a turbulent pipe flow is utilized as an inlet. Additionally, a second new

35
DNS at Re = 7000 was conducted under identical conditions with statistics averaged over
30,000D/Ub . The statistics are compared to the aforementioned studies and the collapse at
various distances from the orifice are shown, which most studies omit. Also, the TKE budget
terms are calculated directly from the simulation and the PDF of the axial velocity is discussed.
Moreover, we use the extensive dataset obtained from these DNS in Chapter 6 to compare
against an extended symmetry theory for computing high moment scaling laws, where we also
address the non-universality of similarity solutions.
The NSEs (2.9) and (2.36) are solved numerically using the CFD code Nek5000 developed by
Fischer et al. (2008). This code is based on a high-order spectral element method (SEM) (see
e.g. Patera, 1984) and is implemented in both Fortran 77 and C. The code is well established in
the scientific turbulence community and is highly scalable over a million ranks using Message
Passing Interface (MPI) for parallelization and has parallel I/O capabilities, utilizing either
MPI-IO or a custom kernel. As a timestepping method, the implicit second order backward
differentiation formula (BDF2) scheme is used and for optimal efficiency the polynomial order
N = 7 has been deployed.
The remarkable scalability of the code enabled simulations at the supercomputer in Munich,
SuperMUC-NG, on 512 nodes, each equipped with 48 cores, for the DNS at Re = 3500 while
the DNS at Re = 7000 employed 2000 nodes which amounts to 24,576 and 96,000 cores,
respectively. Mesh optimization resulted in restart files of approximately 3GB for the Re = 3500
case and about 15GB for the Re = 7000 case.
To extract statistical moments efficiently, a statistics toolbox (Rezaeiravesh et al., 2019) was
employed which computes 44 variables, allowing for the calculation of mean velocities, Reynolds
stress, TKE budgets and Reynolds stress budgets. Initially designed to compute statistical
moments up to the third order on-the-fly, it was extended to accommodate moments up to
the tenth order for the present dissertation. Output for one set of statistical data amounted
to approximately 80GB for the Re = 3500 case and about 415GB for the Re = 7000 case. We
have listed the additional 75 variables that have been extracted in Table 5.1. The passive scalar
and mixed moments are discussed in Chapter 7.
Leveraging the spectral nature of the code, we extracted the statistics at arbitrary points.
Additionally, azimuthal direction averaging was performed during post-processing to account
for the axisymmetric nature of the turbulent jet, yielding statistics in the r and z directions.

36
Table 5.1.: List of 75 variables that have been computed with the extension of the statistics toolbox. Note that these
variables have been computed on top of the 44 variables already calculated by the toolbox (Rezaeiravesh
et al., 2019). The position describes where they are stored in the statistical data output.

Position Variable Position Variable Position Variable

1 Ur 26 Ur2 Θ2 51 Uz6
2 Uϕ 27 Uz2 Θ2 52 Ur7
3 Θ 28 Ur2 Θ3 53 Uz7
4 Uz 29 Uz2 Θ3 54 Ur8
5 P 30 Ur2 Θ4 55 Uz8
6 Ur2 31 Uz2 Θ4 56 Ur9
7 Uϕ2 32 Ur3 57 Uz9
8 Θ2 33 Uz3 58 Ur10
9 Uz2 34 Ur3 Θ 59 Uz10
10 P2 35 Uz3 Θ 60 Θ3
11 Ur Θ 36 Ur3 Θ2 61 Θ4
12 Uz Θ 37 Uz3 Θ2 62 Θ5
13 Ur Uz 38 Ur3 Θ3 63 Θ6
14 Ur Uϕ 39 Uz3 Θ3 64 Θ7
15 Ur Θ2 40 Ur4 65 Θ8
16 Uz Θ2 41 Uz4 66 Θ9
17 Ur Θ3 42 Ur4 Θ 67 Θ10
18 Uz Θ3 43 Uz4 Θ 68 Uϕ3
19 Ur Θ4 44 Ur4 Θ2 69 Uϕ4
20 Uz Θ4 45 Uz4 Θ2 70 Uϕ5
21 Ur Θ5 46 Ur5 71 Uϕ6
22 Uz Θ5 47 Uz5 72 Uϕ7
23 Ur2 Θ 48 Ur5 Θ 73 Uϕ8
24 Uϕ Uz 49 Uz5 Θ 74 Uϕ9
25 Uz2 Θ 50 Ur6 75 Uϕ10

5.1. Computational domain of the direct numerical simulations

We presently adopt a fully developed turbulent pipe flow as an inlet condition, as self-similarity
may occur closer to the inlet. For this, the computational domain has been split into two
domains, i.e., a periodic pipe flow to generate the inlet condition and the main computational
domain to capture the turbulent jet flow. The selection of the two bulk Reynolds numbers,
Re = 3500 and Re = 7000, is deliberate. Experimental observations show that pipe flows are
considered laminar at Re < 2300 and fully turbulent at Re > 2900 (Schlichting and Gersten,
2017). The choice of Re = 3500 falls within the range of fully turbulent pipe flows, ensuring
relevance to established turbulent flow studies. Furthermore, a DNS at this Reynolds number
remains computationally feasible given the available computational resources. On the other
hand, Re = 7000 represents a significantly higher Reynolds number, enabling the identification

37
of observable differences between simulations.

5.1.1. Computational domain of the pipe flow

The computational domain for the pipe in radial, azimuthal and axial direction is r ∈ [0, 0.5D], ϕ ∈
[0, 2π), z ∈ [0, 5D] with 132 cells in the r, ϕ plane and 30 cells in z-direction for the DNS at
Re = 3500 (see Figure 5.2). The Re = 7000 case has 756 cells in the r, ϕ plane and 72 cells
in z-direction. The mesh has the usual refinement towards the near-wall region. Factoring

0.5

r/D 0

−0.5
−0.5 0 0.5
r/D
0.5

r/D 0

−0.5
−5 −4 −3 −2 −1 0
z/D

Figure 5.2.: Cross-sectional view of the computational domain for the pipe at Re = 3500. The N = 7 Gauss-
Lobatto-Legendre (GLL) points has been included in the mesh.

in the N = 7 GLL points, the pipe mesh has around 5.4 million grid points for the Re = 3500
case, while for the Re = 7000 case, the pipe mesh contains around 74 million grid points. The
approach for determining the grid size is given in Appendix A.1.1.

The boundary conditions (BCs) at the wall are no-slip and impermeable wall while at z = −5D
and z = 0 periodicity has been employed to approximate a large infinite system. To generate a
turbulent flow field in the pipe, a small disturbance in all three directions, being sine functions
in the range of 10% of the maximum value of the laminar pipe flow profile, is superposed to

38
Figure 5.3.: A cross-sectional view of the pseudo-color visualized magnitude of the instantaneous velocity field at
Re = 3500 (left) and Re = 7000 (right). The values range from 0 (blue) to 1.4 (red).

the initially laminar pipe flow profile

 
6 1
Uz = 2 2
−r . (5.2)
D 4

To non-dimensionalize all quantities, the pressure gradient and viscosity are set, such that the
bulk velocity is 1 and the Reynolds number is constant. After running the simulation until a
fully developed turbulent pipe flow has been reached (see Figure 5.3), the velocities at the
cross section z = 0 are interpolated onto the inlet of the computational domain of the jet flow
at each timestep which is introduced in the following. In Figure 5.3, it can also be observed
that for a higher Reynolds number, the range of scales increases, although the large scales
dominate in the central region of the flow.

5.1.2. Computational domain of the jet flow

The main computational box of the jet DNS is a cone with a diameter of 4D at the inlet
z/D = 0 and 64D at the outlet z/D = 75 to investigate the near and far-field behavior
of self-preservation in axial direction (see Figure 5.4). This ensures the capture of the jet
spreading while preventing interactions with the lateral far-field boundaries, thus guaranteeing
high-quality statistics.

The computational mesh consists of an inner pipe-type mesh and an outer conical mesh that
expands linearly in the axial direction. As the jet spreads into the far field, the length scales
grow so that the mesh size can be increased gradually in z-direction. Therefore in z-direction,
the computational mesh is geometrically stretched, with a factor of 1.006 resulting in 234 cells
for the Re = 3500 case and with a factor of 1.003 resulting in 275 cells for the Re = 7000 case.
In the outer conical mesh we have the ambient region. Here, a coarser mesh is sufficient which
is achieved by geometrically stretching the cells in r-direction. For the Re = 3500 case, the
stretching factor is 1.06, resulting in 25 cells, while for the Re = 7000 case, the factor is 1.09,
resulting in 36 cells. For the Re = 3500 case, the main computational mesh has around 180,000

39
 Table 5.2.: BCs of the main computational domain
Region Velocity BC

{(r, ϕ, 0)T | r ∈ [0, 0.5], ϕ ∈ [0, 2π)} Interpolated BC

{(r, ϕ, 0)T | r ∈ (0.5, 2), ϕ ∈ [0, 2π)} Dirichlet BC
{(0.4z + 2, ϕ, z)T | ϕ ∈ [0, 2π), z ∈ [0, 75]} Open BC
{(r, ϕ, 75)T | r ∈ [0, 32], ϕ ∈ [0, 2π)} Open BC

cells and, considering the N = 7 GLL points, amounts to 240 million degrees of freedom
(DOFs). For the Re = 7000 case, the mesh includes 900,000 cells yielding 1.24 billion DOFs.
The BC at z/D = 0 outside of the jet inlet has been set as a no-slip wall, while the lateral
boundaries and the jet exit at z/D = 75 are configured as open boundaries following the
approach introduced by Dong et al. (2014). The open boundary is designed to locally prevent
the kinetic energy influx into the domain at the outflows by balancing the energy influx with
the effective stress, addressing a key challenge in DNS of jet flows that leads to instabilities.
Otherwise, it functions as a standard open outflow BC with zero pressure, thereby enabling
entrainment of the fluid into the boundary. The BCs have been summarized in Table 5.2.
Additionally, special attention has been paid to obtaining very accurate statistics, and for all the
statistics in this dissertation for the Re = 3500 case, averages have been taken over 75,000D/Ub
times. Despite the significantly larger box, this corresponds to nearly 200 passes of a particle
through the entire computational domain of the jet and is a factor of 25 longer than the largest
simulations to date which, to the author’s knowledge, was conducted in Sharan and Bellan
(2021). Due to limited computational resources, the statistics for the DNS at Re = 7000 have
only been averaged over 30,000D/Ub , which is still a factor of 10 longer than the DNS in the
previously mentioned work.
For both simulations, the target Courant number has been set to 2.5 with an initial timestep
of ∆t = 5 · 10−3 for the Re = 3500 case and ∆t = 2.5 · 10−3 for the Re = 7000 case which
has been adjusted dynamically to satisfy the Courant number during the simulation. The
simulations have run for a combined 95 million core-hours on SuperMUC-NG.
Figure 5.5 shows a snapshot of isosurfaces of the q-criterion (Hunt et al., 1988) at q = 0.01 for
the DNS at Re = 3500. Here, q serves as a criterion used to capture vortices and is defined as
1
(tr(∇u))2 − tr(∇u · ∇u) , (5.3)

q=
2
where tr(·) is the trace of a square matrix. It can be observed that the transition to turbulence
occurs close to the orifice which can be identified by the spread of the jet compared to earlier
studies (see Figure 5.6). In these studies, the spreading, and hence the decay zone, starts
at a greater distance from the orifice. The virtual origin can be determined by finding the
theoretical source of the spreading, which leads to z0 = 0 in Equation (5.4) below.
In addition to the NSEs, a convection-diffusion equation for a passive scalar with a Prandtl
number P r = 0.71 has been solved on the main computational box. The passive scalar is
introduced at the inlet as a constant with Θ = 1. Similarly to the velocity BC, an energy-stable
open BC (Liu et al., 2020) has been implemented on the lateral boundaries and the jet exit.
This setup and the results are discussed in detail in Chapter 7.

40
Figure 5.4.: Cross-sectional view of the main computational box for the Re = 3500 case at z/D = 0 (left) which is
scaled linearly in z-direction to obtain the whole main computational box (right). The GLL points are
omitted for better visibility. The DNS uses curved elements which are not apparent in this figure.

1.5
velocity magnitude

0.5

0.0
Figure 5.5.: A cross section of q-criterion isosurfaces at q = 0.01 of the conducted jet DNS at Re = 3500 colored
with the velocity magnitude.

5.2. Mean velocity statistics

Classically, the mean velocities of a self-preserving jet are scaled with the centerline velocity
(Boersma et al., 1998)
U z,c (z) Bu D
= , (5.4)
U0 z − z0
where
U 0 = 1.38 (Re = 3500), U 0 = 1.29 (Re = 7000) (5.5)

refer to the mean axial centerline velocity at the orifice for the DNS at the corresponding
Reynolds numbers, Bu to the decay constant and z0 to the virtual origin. The mean axial
centerline velocity at the orifice (5.5) corresponds to the centerline velocity of the turbulent
pipe flow. It is well known, that the velocity gradient near the wall grows with increasing
Reynolds number (Batchelor, 2000). In order to satisfy the constant bulk velocity condition

41
(see. Section 5.1.1), the centerline velocity must decrease, which explains the different values
of (5.5) for the two cases.

In Figure 5.6 the inverse of the mean axial centerline velocity according to Equation (5.4)
is compared with other DNS and experimental data. Presently, the Reynolds decomposition
for the velocity in Equation (3.12) is adopted. Values for the present DNS and those from

14

12

10
U 0 /U z,c

6
1.1
4 1.05

1
2
0.95
2 4 6 8 10
0
0 10 20 30 40 50 60 70
z/D

Figure 5.6.: Inverse of the mean axial centerline velocity according to Equation (5.4) plotted over the distance
from the orifice. Present DNS at Re = 3500 ( ), present DNS at Re = 7000 ( ), Boersma et al.
(1998) at Re = 2400 ( ), Babu and Mahesh (2004) at Re = 2400 ( ), Taub et al. (2013) at
Re = 2000 ( ). The inset figure with the same axis shows a magnification at the near-field, where
the cutoff ( ) marks the end of the potential core. Both present simulations there have an earlier
onset of spreading compared to the comparative studies.

other studies are provided in Table 5.3. The present values of the parameters have been
determined by minimizing the sum of the squares of the residuals relative to the DNS data
on the conservative interval z/D = 15 − 65. We purposely do not include the last part of
the box due to boundary effects. Nevertheless, the average relative residual on the interval
z/D = 15 − 75 is 0.2%. The decay rate Bu of the present DNS is slightly lower compared to
previous findings. Additionally, in the DNS at Re = 3500, the virtual origin is z0 = 0, attributed
to early mixing induced by the fully turbulent inlet. The DNS at Re = 7000 is characterized
by z0 = 2, also leading to a later onset of self-similarity. This delay in self-similarity onset
compared to the DNS at the lower Reynolds number is consistently observed in all the DNS
results throughout this dissertation. The inset figure showcases a magnified view close to
the inlet which provides a cutoff where the so-called potential core ends. A potential core
extends from the orifice up to where the mean axial centerline velocity reaches 95%, which
is followed by the decaying zone (Jambunathan et al., 1992). Further, the size of potential
cores is reduced with increasing turbulence intensities and shear layers. In Figure 5.6, the
potential core extends to U 0 /U z,c = 1/0.95 = 1.053, which is indicated by the dashed line
in the magnified view. It can be observed that the potential core of both present simulations
break down significantly earlier than those of the other studies. This is due to the turbulence
introduced by the turbulent inlet which exhibited later in Figures 5.10 and 5.11.

42
 Reference Re Bu z0 η1/2
Wygnanski and Fiedler (1969) 105 5.7 3.0 0.085
Panchapakesan and Lumley (1993a) 1.1 · 104 6.06 0.0 0.096
Hussein et al. (1994) 95.5 · 103 5.8 4.0 0.094
Boersma et al. (1998) 2400 5.9 4.9 0.091
Taub et al. (2013) 2000 5.4 1.3 0.096
Present DNS 3500 5.15 0.0 0.089
7000 5.25 2.0 0.086

Table 5.3.: Various jet parameters of DNS and experiments.

In classical jet theory, the radial coordinate is non-dimensionalized with

r
η= . (5.6)
z − z0
For the jet half-width η1/2 , representing the position where the velocity reaches half the value
of the centerline velocity, we find η1/2 = 0.089 for Re = 3500 and η1/2 = 0.086 for Re = 7000
(see Figure 5.7). These values are in the range of what has been observed previously, as shown
in Table 5.3. The differences are due to the variations in the inlet conditions utilized in different
studies. In Figure 5.8 the mean radial velocity U r /U z,c reaches a maximum outward velocity

1 1

0.8 0.8
Uz /U z,c

0.6 0.6

0.4 0.4

0.2 0.2

0 0
0 0.05 0.1 0.15 0.2 0.25 0 0.05 0.1 0.15 0.2 0.25
η η

Figure 5.7.: Mean axial velocity profiles normalized with the axial centerline velocity according to Equation (5.7)
as functions of η at different distances from the orifice for Re = 3500 (left) and Re = 7000 (right):
z/D = 25 ( ), 35 ( ), 45 ( ), 55 ( ), 65 ( ).

of U r,,max /U z,c = 0.017 at η = 0.056 for both Reynolds numbers. This velocity value is similar
to the value found by Wygnanski and Fiedler (1969) at U r,max /U z,c = 0.015, Panchapakesan
and Lumley (1993a) at U r,max /U z,c = 0.018 and Taub et al. (2013) at U r,max /U z,c = 0.02.
The maximum inward velocity is at η = 0.22 with U r,min /U z,c = −0.022 for the Re = 3500
case and slightly lower at U r,min /U z,c = −0.021 for the Re = 7000 case at the same position.
It is worth noting that the values of U r on the centerline and of U ϕ in the whole domain are
very small and in order of 10−4 , implying a well converged dataset.

43
The Figures 5.7 and 5.8 show an excellent collapse of the DNS data, indicating self-similarity,
according to the classical scaling

Uei (η)
Ui (r, z) = (5.7)
(z − z0 )

for the mean velocity in the range z/D = 25 − 65, where Uei (η) refers to the self-similar solution,
which in the context of Lie theory is an invariant, i.e, the collapsed profiles only dependent of η.
In Figures 5.7 and 5.8, Ui (r, z) has been normalized by the mean axial centerline velocity U z,c
in (5.4) to obtain the invariant. Subsequently, and in Chapter 6, this is also observed for velocity
moments up to n = 10. Furthermore, Figure 5.9 showcases the invariant U f (η = 0) = B DU
z u 0
according to Equation (5.7) being valid in the far-field in z-direction which further supports
the 1/z decay depicted in Figure 5.6. It can be observed that the statistics of the present DNS
are very well converged up to the end of the domain compared to previous studies. The lower
value compared to the other studies is due to the lower decay rate Bu .

·10−2 ·10−2

2 2

1 1
Ur /U z,c

0 0

−1 −1

−2 −2

0 0.05 0.1 0.15 0.2 0.25 0 0.05 0.1 0.15 0.2 0.25
η η

Figure 5.8.: Mean radial velocity profiles normalized with the axial centerline velocity according to Equation (5.7)
as functions of η at different distances from the orifice for Re = 3500 (left) and Re = 7000 (right):
z/D = 25 ( ), 35 ( ), 45 ( ), 55 ( ), 65 ( ).

44
 8

6
e (η = 0)

4
U z

0
0 10 20 30 40 50 60 70
z/D

Figure 5.9.: The invariant Ue (η = 0) plotted over the distance from the orifice. Present DNS at Re = 3500 (
z ),
present DNS at Re = 7000 ( ), Boersma et al. (1998) ( ), Babu and Mahesh (2004) ( ),
Taub et al. (2013) ( ).

5.3. Turbulent intensities

Figures 5.10 and 5.11 compare the components of the turbulent intensity on the centerline with
the DNS from Taub et al. (2013), LES from Bogey and Bailly (2009) and an experiment from
Panchapakesan and Lumley (1993a). In Figure 5.10, a fast increase of the radial component of
the turbulent intensity is observed starting at z/D = 70 which is due to the boundary effects
which are counteracted by employing the outflow boundary condition by Dong et al. (2014),
preventing an uncontrolled growth in the energy imposed by the boundary. This increase is
transported by the large eddies from the boundary into the domain at the end of the box. This
effect, if ever so slightly, can be also observed in Figure 5.9 close to z/D = 75. In Figure 5.11,
this effect does not become apparent although it can be assumed that it does marginally.

In Taub et al. (2013), it is mentioned that both in theirs and Bogey’s study an almost constant
value above the transition region is attained and that experiments reach an asymptotic value
gradually. This seems, however, only a first approximation, as the variance in the turbulent
intensity of the study by Taub et al. (2013) is still quite high and the LES by Bogey and Bailly
(2009) only covers a small axial distance compared to the present simulations.

There are mainly two methods to shorten the transition region, the first one being a smaller
Reynolds number (Bogey and Bailly, 2009) and the second one being a faster initiation of the
turbulent behavior. The incline of the present DNS in Figure 5.9 starts at z/D = 0 compared
to the other DNS studies showing that the transition begins earlier. This is further supported
by the earlier decay of the potential cores in Figure 5.6. However, it is not apparent that a
constant value is achieved earlier which might be due to the higher Reynolds number of the
present DNS at Re = 3500 compared to the other studies where the range is Re = 2000 − 2400.
Additionally, in Fig. 5.10 and 5.11, we observe a non-zero turbulent intensity at z/D = 0
which has not been achieved by the other studies shown in those figures. Taub et al. (2013) has
superimposed small non-physical perturbations onto a top-hat profile, while Bogey and Bailly
(2009) has used divergence free vortex rings near the inlet. Random perturbations introduced
on the velocity inlet profile seems therefore not sufficient to initiate fully turbulent behavior.
Therefore, we conclude that using a turbulent velocity profile may shorten the transition region.
However, a study to confirm this by comparing different inlet conditions at a constant Re

45
would be beneficial. The behavior of the present DNS is closer to the experimental study

0.3

0.25

0.2
u2r,c /U z,c

0.15
q

0.1

0.05

0
0 10 20 30 40 50 60 70 80
z/D

Figure 5.10.: Radial component of the turbulent intensity on the centerline compared to previous studies: present
DNS at Re = 3500 ( ), present DNS at Re = 7000 ( ), Taub et al. (2013) ( ), Bogey and
Bailly (2009) ( ), Panchapakesan and Lumley (1993a) ( ).

by Panchapakesan and Lumley (1993a) as the turbulent intensities slowly increase in the z
far-field region.

0.3

0.25

0.2
u2z,c /U z,c

0.15
q

0.1

0.05

0
0 10 20 30 40 50 60 70 80
z/D

Figure 5.11.: Axial component of the turbulent intensity on the centerline compared to previous studies: present
DNS at Re = 3500 ( ), present DNS at Re = 7000 ( ), Taub et al. (2013) ( ), Bogey and
Bailly (2009) ( ), Panchapakesan and Lumley (1993a) ( ).

46
5.4. Reynolds stress
2
According to the classical scaling, all Reynolds stresses of the present DNS are scaled with U z,c
in (5.4), i.e.,
ug
i uj (η)
ui uj (r, z) = , (5.8)
(z − z0 )2
and are shown in Figure 5.12 for the Re = 3500 case and in Figure 5.13 for the Re = 7000
case, where ug i uj (η) refers to the invariant or similarity variable similar to Equation (5.7).
The profiles are presented for different axial distances from the orifice exhibiting an excellent
collapse for the Re = 3500 case, while the Re = 7000 case reaches self-similarity at a later
stage. The maximum values are slightly larger in the Re = 7000 case for the normal stresses
but similar for ur uz .
As mentioned by the previous works (Hussein et al., 1994; Panchapakesan and Lumley, 1993a;
Taub et al., 2013), the off-center maxima are due to the strong interaction of the jet with
the ambient fluid in that region. In the present result for uϕ uϕ this maximum is even more
dominant than in Taub et al. (2013) and Panchapakesan and Lumley (1993a).

47
 ·10−2 ·10−2
4
4

3
3
ui uj /U z,c
2

2
2

1 1

0 0
0 0.1 0.2 0.3 0 0.1 0.2 0.3
η η
(a) ij = rr (b) ij = ϕϕ
·10−2 ·10−2
6 2

1.5
4
ui uj /U z,c
2

2
0.5

0 0
0 0.1 0.2 0.3 0 0.1 0.2 0.3
η η
(c) ij = zz (d) ij = rz

2
Figure 5.12.: Reynolds stresses ui uj normalized with the axial centerline velocity U z,c according to Equation (5.8)
at different distances from the orifice for the Re = 3500 case: z/D = 25 ( ), 35 ( ), 45 ( ),
55 ( ), 65 ( ).

48
 ·10−2 ·10−2

4 4

3 3
ui uj /U z,c
2

2 2

1 1

0 0
0 0.1 0.2 0.3 0 0.1 0.2 0.3
η η
(a) ij = rr (b) ij = ϕϕ
·10−2 ·10−2
2

6
1.5
ui uj /U z,c
2

4
1

2
0.5

0 0
0 0.1 0.2 0.3 0 0.1 0.2 0.3
η η
(c) ij = zz (d) ij = rz

2
Figure 5.13.: Reynolds stresses ui uj normalized with the axial centerline velocity U z,c according to Equation (5.8)
at different distances from the orifice for the Re = 7000 case: z/D = 25 ( ), 35 ( ), 45 ( ),
55 ( ), 65 ( ).

49
5.5. Reynolds stress transport budgets

The budget equations of the Reynolds stresses in Equation (3.22) are expressed as

Cij + Πij + Pij + Tij + εij + Vij = 0, (5.9)

|{z} |{z} |{z} |{z} |{z} |{z}
convection velocity-pressure production turbulent diffusion dissipation viscous diffusion
gradient correlation

see, e.g., Hoyas and Jiménez (2008) and Mansour et al. (1988), where i, j = r, ϕ, z and the
local temporal term Lij is zero. Here, the terms denote convection, velocity-pressure gradient
correlation, production, turbulent diffusion, dissipation and viscous diffusion, respectively.
Due to the statistical axisymmetry of the flow all terms ( · )ϕz and ( · )rϕ are zero. Further, the
viscous diffusion is negligible due to a high Reynolds number (Taub et al., 2013).

The five non-zero transport equations from (3.22) in cylindrical coordinates are shown in
Equations (5.10)-(5.13) while the brackets refer to the terms in Equation (5.9), respectively.
The viscous term has already been neglected as the present DNS also show that the term is
close to zero:
ur ur transport equation:
     
∂ ∂ ∂p ∂U r ∂U r
0 = − U z ur ur + U r ur ur − 2ur − 2ur ur + 2uz ur
∂z ∂r ∂r ∂r ∂z
| {z } | {z } | {z }
Crr Πrr Prr
" # " 2 2 2 # (5.10)
∂uz u2r ur u2ϕ
 
1 ∂  3 2 ∂ur ∂ur 1 ∂ur
− + rur − 2 − + + 2
∂z r ∂r r Re ∂z ∂r r ∂ϕ
| {z } | {z }
Trr rr

uϕ uϕ transport equation:
     
∂ ∂ uϕ ∂p Ur
0 = − U z uϕ uϕ + U r uϕ uϕ − 2 − 2 uϕ uϕ
∂z ∂r r ∂ϕ r
| {z } | {z } | {z }
Cϕϕ Πϕϕ Pϕϕ
" # " 2  2 2 # (5.11)
∂uz u2ϕ u u2 
1 ∂   r ϕ 2 ∂u ϕ ∂u ϕ 1 ∂u ϕ
− + rur u2ϕ + 2 − + + 2
∂z r ∂r r Re ∂z ∂r r ∂ϕ
| {z } | {z }
Tϕϕ ϕϕ

uz uz transport equation:
     
∂ ∂ ∂p ∂U z ∂U z
0 = − U z uz uz + U r uz uz − 2uz − 2uz uz + 2uz ur
∂z ∂r ∂z ∂z ∂r
| {z } | {z } | {z }
Czz Πzz Pzz
" # " 2 2 2 # (5.12)
∂u3z
 
1 ∂  2  2 ∂uz ∂uz 1 ∂uz
− + ruz ur − + + 2
∂z r ∂r Re ∂z ∂r r ∂ϕ
| {z } | {z }
Tzz zz

50
ur uz transport equation:
   
∂ ∂ ∂p ∂p
0 = − U z ur uz + U r ur uz − ur + uz
∂z ∂r ∂z ∂r
| {z }
Crz
  " 2  uz u2
#
∂U z ∂U z ∂U r ∂U r ∂uz ur 1 ∂  2 ϕ
− ur uz + ur ur + uz uz + ur uz − + ruz ur −
∂z ∂r ∂z ∂r ∂z r ∂r r
| {z } | {z }
Πrz Trz
 
2 ∂ur ∂uz ∂ur ∂uz 1 ∂ur ∂uz
− + + 2
Re ∂z ∂z ∂r ∂r r ∂ϕ ∂ϕ
| {z }
rz
(5.13)

The velocity-pressure gradient correlation term offers more information when split into a
pressure diffusion and pressure strain term Πij = Πdij − Πsij
     
∂p ∂ ∂ur
− 2ur = − 2 pur + 2p
∂r ∂r ∂r
     
2 ∂p 2 ∂ 2 ∂uϕ
− uϕ = − puϕ + p
r ∂ϕ r ∂ϕ r ∂ϕ
     
∂p ∂ ∂uz
− 2uz = − 2 puz + 2p (5.14)
∂z ∂z ∂z
     
∂p ∂p ∂ ∂ ∂ur ∂uz
− ur + uz = − pur + puz + p +p ,
∂z ∂r ∂z ∂r ∂z ∂r

where the terms on the right side in brackets correspond to the pressure diffusion and pressure
strain, respectively.
Each tensor Eij in (5.9) can be rescaled according to

Eij z
Ẽij = 3 (5.15)
U c,z

using the mean axial centerline velocity in (5.4). Employing the present DNS data, all terms in
(5.9) can be viewed in Figures 5.14 and 5.15 for the Re = 3500 case and the Re = 7000 case,
respectively. The terms are calculated in Cartesian coordinates which are then transformed to
its cylindrical form and averaged in azimuthal direction, subsequently.
For sake of readability, only the curves for z/D = 25 are shown in Figure 5.14, although budgets
also possess an excellent self-similarity. This deliberate choice stems from the observation,
as visible in Figure 5.12, of weak deviations near the jet axis for the second moments of the
fluctuations, despite the extensive statistics conducted in this study.
In Figure 5.15 the curves are shown for a later stage at z/D = 35 since the budget terms also
show self-similarity further away from the orifice as observed in Figure 5.13. The budget terms
of the Re = 7000 case are not as well converged as the budget terms of the Re = 3500 case
due to limited computational resources and the large Reynolds number. While the qualitative

51
 0.1 0.1
( · )ij /(U z,c /z)

0 0
3

−0.1 −0.1

−0.2 −0.2
0 0.05 0.1 0.15 0.2 0.25 0 0.05 0.1 0.15 0.2 0.25
η η
(a) ij = rr (b) ij = ϕϕ
0.4

0.2
0.2
( · )ij /(U z,c /z)
3

0
0

−0.2 −0.2

0 0.05 0.1 0.15 0.2 0.25 0 0.05 0.1 0.15 0.2 0.25
η η
(c) ij = zz (d) ij = rz

2
Figure 5.14.: Budgets of the Reynolds stresses at Re = 3500 scaled with U z,c /z at z/D = 25: production ( ),
dissipation ( ), turbulent diffusion ( ), velocity-pressure gradient correlation ( ), convection
( ), pressure diffusion ( ), pressure strain ( ).

behavior of the budget terms remains mostly similar, quantitative comparison is challenging
due to data variance.

From (5.9), the equation for the TKE can be derived k = 1/2 ui ui
∂k
= C k + Πk,d + P k + T k + εk + V k = 0, (5.16)
∂t
where the TKE budget terms are defined as
1
E k = Eii . (5.17)
2
The TKE budget terms and their sum are shown in Figure 5.16 for both Reynolds numbers. For
the same reasons as for the budgets of the Reynolds stresses, here we also focus only on values

52
 0.2

0.1
0.1
( · )ij /(U z,c /z)

0
3

−0.1
−0.1

−0.2 −0.2
0 0.05 0.1 0.15 0.2 0.25 0 0.05 0.1 0.15 0.2 0.25
η η
(a) ij = rr (b) ij = ϕϕ
0.4

0.2
0.2
( · )ij /(U z,c /z)
3

0
0

−0.2 −0.2

0 0.05 0.1 0.15 0.2 0.25 0 0.05 0.1 0.15 0.2 0.25
η η
(c) ij = zz (d) ij = rz

3
Figure 5.15.: Budgets of the Reynolds stresses at Re = 7000 scaled with U z,c /z at z/D = 35: production ( ),
dissipation ( ), turbulent diffusion ( ), velocity-pressure gradient correlation( ), convection
( ), pressure diffusion ( ), pressure strain ( ).

at one distance from the orifice. All terms have been again calculated directly from the DNS
according to (5.15).

The sum of the budget terms, i.e, the error in the TKE balance, of the Re = 3500 case is below
0.0126 and corresponds to 6% error compared to the maximum value of the dissipation for all
values of η in the same scaling as in Figure 5.16 at Re = 3500, which is contrasted to the error of
0.02 in the DNS of Taub et al. (2013). The error of the TKE balance in Figure 5.16 at Re = 7000
is higher at 0.08 corresponding to a 34% error compared to the maximum dissipation value.
The larger error is due to the reasons mentioned above.

The viscous diffusion is close to zero and is therefore not included in both figures. The
production and dissipation term in Figure 5.16 at Re = 7000 are well converged and are

53
slightly larger than in Figure 5.16 at Re = 3500.

The dissipation on the centerline is about twice as large as the convection term. This agrees
with Hussein et al. (1994), Wygnanski and Fiedler (1969) and Darisse et al. (2015) who
conducted experiments at Reynolds numbers in the order of 105 , whereas Panchapakesan and
Lumley (1993a) and Taub et al. (2013), who examined jet flows at Reynolds numbers in the
order of 103 − 104 , found that the values are about the same. In the works of Bogey and
Bailly (2006) and Taub et al. (2013) it is speculated that this might be due to the difference in
Reynolds number even for high Reynolds numbers.

For the Re = 7000 case, a precise comparison is difficult due to the large error at the centerline.
However, the dissipation term is larger than the convection term rather than similar. Since
the studies above have been using top-hat profiles as the inlet condition, this could mean that
a fully turbulent velocity profile at the inlet of the present DNS influences the dissipation
similarly to a higher Reynolds number, i.e., the similarity profiles is dependent of the inlet
condition.

Further, in Figure 5.17, the dissipation of the present DNS is compared to some previous studies
where the dissipation profiles are similar but reach different magnitudes. Only the studies
of Taub et al. (2013), Bogey and Bailly (2009) and the present study have calculated the
dissipation directly, while Hussein et al. (1994) have approximated the dissipation by assuming
local isotropy. Taub et al. (2013) have also calculated the dissipation using the approximation
by Hussein et al. (1994) where the dissipation has been overestimated. Also, Bogey and Bailly
(2009) have conducted a LES, which in itself already contains models. The dissipation of the
present DNS lies in the range of all studies except for Hussein et al., 1994. It is difficult to draw
conclusions by comparing to the other studies since the conditions differs strongly. However,
for the present DNS setup the dissipation grows with increasing Reynolds number.

0.15

0.1
0.1
0.05
( · )k /(U c,z /z)

0 0
3

−0.05

−0.1 −0.1

−0.15
−0.2
−0.2
0 0.05 0.1 0.15 0.2 0.25 0.3 0 0.05 0.1 0.15 0.2 0.25 0.3
η η

Figure 5.16.: Turbulent kinetic energy budgets for the Re = 3500 case at z/D = 25 (left) and for the Re = 7000
case at z/D = 35 (right): production ( ), dissipation ( ), turbulent diffusion ( ), convection
( ), pressure diffusion ( ), sum of the budget terms ( ).

54
 0.4
Present DNS at Re = 3500
Present DNS at Re = 7000
0.35
Taub at Re = 2000
Panchapakesan at Re = 11,000
0.3
Hussein at Re = 100,000
Bogey at Re = 11,000
0.25
|k |/(U c,z /z)
3

0.2

0.15

0.1

5 · 10−2

0
0 5 · 10−2 0.1 0.15 0.2 0.25
η

Figure 5.17.: The dissipation of the turbulent kinetic energy comparing the present DNS at Re = 3500 and
Re = 7000 to the DNS of Taub et al. (2013) at Re = 2000, the experiments of Panchapakesan and
Lumley (1993a) at Re = 11,000, Hussein et al. (1994) at Re = 95,500 and a LES of Bogey and Bailly
(2009) at Re = 11,000.

The previous results of the dissipation may now be used to evaluate Kolmogorov‘s relation
2
k3
 = C , (5.18)
l
where l is the integral length scale
ZZZ
1
l= Rmm d3 r, (5.19)
k

Rmm is the two-point correlation, and r the distance between two points. Implementing
(5.8) and (5.15) into Kolmogorov‘s relation (5.18), where we extended (5.8) to two-point
correlations, we find that z drops out on both sides and hence C is not depending on the local
Reynolds number as was found for certain other flows (Vassilicos, 2015).
Moreover, the dissipation can serve as a means to ensure that the grid size adequately resolves
the smallest scales, as exemplified for the Re = 3500 case, detailed in Appendix A.1.2.

5.6. Probability density functions

In addition to the classical variables of mean velocity, Reynolds stresses, and budgets, extensive
DNS statistics enable the computation of higher-order statistics, i.e., probability density func-
tions (PDF), at high precision. The PDFs are created by collecting Uz at various points every

55
 101 101

100 100

10−1 10−1
f
10−2 10−2

10−3 10−3

10−4 10−4
0 0.5 1 1.5 2 2.5 0 0.5 1 1.5 2 2.5
Uz /U z,c Uz /U z,c

Figure 5.18.: PDFs of Uz (η = 0, z)/U z,c (z) at Re = 3500 (left) and Re = 7000 (right) for z/D = 15 ( ), 25
( ), 35 ( ), 45 ( ), 55 ( ), 65 ( ) compared to a Gaussian ( ).

0.05D/Ub over the span of 75,000D/Ub at Re = 3500 and every 0.025D/Ub over the span of
30,000D/Ub at Re = 7000. The data is then divided in bins where data of a given r and z can
be collected in the same PDF, which contributes to the smoothness of the resulting PDFs. After
binning the data, the probability density for each bin is calculated. This involves counting
the number of data points falling within each bin. To ensure that the PDFs represent true
probabilities, these densities are normalized so that they sum up to 1.
In Figure 5.18, the PDF f of the normalized axial centerline velocity for both Reynolds numbers
also show an excellence collapse for different distances z/D from the orifice. On the centerline,
the PDFs f are well described by a Gaussian distribution.
Additionally, the PDF development at z/D = 28, 42, 56 for varying η is shown in Figure 5.19.
It can be clearly seen that for increasing η the curves slowly deviate from a Gaussian and
non-Gaussian tails develop. By further increasing η, a delta distribution is approached due
to a laminar constant ambiance velocity. Most important, we observe an excellent PDF data
collapse due to the η scaling.
Finally, in Figure 5.20, the Uz -based skewness S and kurtosis K are extracted from the data of
the Re = 3500 and Re = 7000 cases, respectively. These values also show an almost perfect
collapse in the range of similarity. A skewness of S = 0 and kurtosis of K = 3 correspond to
characteristics of a Gaussian distribution. This indicates that on the jet centerline η = 0, the
PDFs closely resemble a Gaussian distribution which supports the insight gained in Figure 5.18.
As η increases, the PDFs show increasing positive skewness with K = 3 until η = 0.09. Beyond
this point, the PDFs deviate from Gaussian behavior, becoming increasingly non-Gaussian as
they approach a delta distribution in the laminar region of the flow.

56
 100
f
10−3

0.18 0.2
10−6 0.12 0.14 0.16
2 0.08 0.1
0 0 0.02 0.04 0.06
η
Uz /U z,c

100
f
10−3

0.18 0.2
10−6 0.12 0.14 0.16
2 0.08 0.1
0 0 0.02 0.04 0.06
η
Uz /U z,c

Figure 5.19.: PDFs of Uz (η, z)/U z,c (z) at Re = 3500 (above) and Re = 7000 (below) for z/D = 28 ( ), 42
( ), 56 ( ).

15 15

10 10
K, S

5 5

0 0
0 0.05 0.1 0.15 0.2 0 0.05 0.1 0.15 0.2
η η

Figure 5.20.: Kurtosis K (above) and skewness S (below) of Uz (η, z)/U z,c (z) at Re = 3500 (left) and Re = 7000
(right) for z/D = 28 ( ), 42 ( ), 56 ( ). K = 3, S = 0 (dashed) are the Gaussian values.

57
5.7. Conclusive remarks

In conclusion, turbulent round jet flows at Re = 3500 and Re = 7000 have been successfully
simulated. The size of the box extends up to 75D in axial and up to 65D in radial direction
to ensure the caption of the far-field statistics. Statistics are collected over 75,000D/Ub units
for the Re = 3500 case, equivalent to 200 passes of a particle through the domain. For the
Re = 7000 case, statistics are collected over 30,000D/Ub units, equating to 80 passes. Previous
studies have not achieved a DNS at this box size with statistics of comparable quality.
High-quality statistics have been generated for the velocity statistics of first and second order
one-point moments. Furthermore, all the terms of the Reynolds stress equations and the TKE
equation have been calculated directly from the DNS data. Additionally, for the axial velocity,
PDFs have been extracted in radial and axial direction, a novel contribution not reported in
previous literature. The mean axial centerline velocity shows a highly converged 1/z behavior
throughout the entire domain. All the radial velocity profiles of the first and second order
statistics show a remarkable collapse based on classical scaling. Notably, all radial velocity
profiles of first and second order statistics demonstrate remarkable collapse based on classical
scaling. The observed slight axial increase in turbulent intensities aligns with experimental
findings, which previous DNS studies failed to capture due to limitations in box size.
Employing a fully turbulent velocity profile as an inlet reduces the the size of the potential core
significantly leading to an earlier decay. Additionally, this velocity profile may have similar
effects on dissipation an increasing Reynolds number which could be confirmed by comparing
to a jet simulation with a top-hat profile at identical conditions.
From comparing the DNS results at the two Reynolds numbers, we identified a shift in the virtual
origin, with the magnitude of the Reynolds normal stresses increasing with higher Reynolds
numbers, while the Reynolds shear stress remains constant. Analysis of turbulent intensities
indicates that self-similarity shifts to greater distances from the orifice with increasing Reynolds
number, a trend observed across other statistical quantities.
Furthermore, the TKE budgets highlight larger production and dissipation rates for the higher
Reynolds number case. However, it is essential to note that the budget terms of Reynolds stress
are not fully converged due to limited computational resources, evident from the large error in
the TKE budgets for the Re = 7000 case.
In the next chapter, we investigate hidden intermittency in a turbulent round jet based on
symmetry methods and the present DNS data set.

58
6. Lie Symmetry Analysis of a Turbulent Round
Jet Flow

This chapter is heavily based on the publication Nguyen and Oberlack (2024c) and repre-
sents the main results of the Lie symmetry analysis of a turbulent round jet flow for pure
hydrodynamics.
Most likely Karman and Howarth (1938) were first to suggest that turbulent scaling laws are
similarity solutions of two-point moment equations for the decay of isotropic turbulence. For
this, they employed a power-law ansatz though this was limited to second moments. For
various turbulent shear flows, classical scalings have been derived based on self-preservation
analysis by Townsend (1956, 1976).
Both laminar and turbulent jets exhibit a strong tendency to self-similarity, and the respective
scaling laws appear to be a strong attractor. For the laminar round jet flow, an exact solution
has been generated by Schlichting (1933), again based on a similarity ansatz. For a turbulent
round jet, it has been observed that mean axial velocity profiles closely follow a Gaussian
(List, 1982). Mentions of this in plane jets lead back to Bradbury (1965). Investigations of
higher moments have been usually limited to the second and third moments, i.e., the Reynolds
stresses and their transport budgets (Darisse et al., 2015; Hussein et al., 1994; Panchapakesan
and Lumley, 1993a).
In Oberlack (2001), symmetries have been applied for the first time to turbulent shear flows
to generate invariant solutions, equivalent to turbulent scaling laws, for the mean velocity
of stationary parallel turbulent shear flows. However, this approach was limited to the mean
velocity at that time. We shortly recapitulate from Chapter 3 that in the RANS equation,
the moment Hij = Ui Uj is decomposed into Hij = U i U j + ui uj where Rij = ui uj is the
Reynolds stress tensor. With Hij and Rij , constructing moments based on the instantaneous
velocities shall be called the H-approach and the classical approach is called the R-approach.
Instantaneous multi-point moments of arbitrary order n are denoted by Hi{n} and one-point
moments by Uin . The corresponding MPME are an infinite set of linear equations with the only
coupling between “neighboring” equations. Though mathematical fully equivalent (Oberlack
and Rosteck, 2010), this is in contrast to the multi-point R-approach, which is non-linear and
exhibits a strong coupling among various orders of moments.
By rigorously extracting symmetry transformations from the considerably simpler H-MPMEs,
invariant solutions can be derived consistently to be shown below. For temporally evolving
turbulent plane jets, this has been achieved for velocity moments (Sadeghi et al., 2018). Further,
just recently in Oberlack et al. (2022), turbulent scaling laws of arbitrary H-moments for the
log and the core region of a turbulent channel flow were derived and validated using a new
DNS at a friction Reynolds number of Reτ = 10,000. Key ingredients are statistical symmetries

59
identified in Oberlack and Rosteck (2010) and Wacławczyk et al. (2014) as they are only
admitted by the MPMEs and not by the NSEs. One of them defines a measure of intermittency
(Wacławczyk et al., 2014).

6.1. Constructing the turbulent velocity scaling laws

Subsequently, invariant solutions, and hence turbulent velocity scaling laws, are derived
using Lie symmetry analysis as introduced in Chapter 4. To recapitulate, symmetries are
defined as variable transformations that do not
modify the form of the differential equation.
Consider a system of PDEs F x, y, y (1) , ..., y (p) = 0 and a transformation x∗ = φ(x, y; a) and
y ∗ = ψ(x, y; a), where x, y, a are the independent variables, dependent variables, respectively,
and a ∈ R a continuous parameter. Then, that transformation is a Lie symmetry if
 
F x, y, y (1) , ..., y (p) = 0
  (6.1)
⇐⇒ F x∗ , y ∗ , y ∗(1) , ..., y ∗(p) = 0

holds, where y (p) refers to the pth derivative of y.

6.1.1. Symmetries of the multi-point moment equations

The H-MPMEs (3.30) are considerably more simple for a symmetry analysis than the R-
equations, such as the RANS equations (3.16) and the Reynolds stress transport equations
(3.22), due to their linearity and coupling only being between the order n and n + 1.

Presently, the following symmetries of the MPMEs (3.30) are relevant for the turbulent round
jet flow. A translation in space
∗
T x : t∗ = t, x∗i = xi + axi , Ui = Ui , Hi∗{n} = Hi{n} , (6.2)

and the scaling in space and time characterized by their group parameters aSx and aSt ,
respectively,
∗
TS : t∗ = eaSt t, x∗i = eaSx xi , Ui = eaSx −aSt Ui ,
∗
Hi∗{n} = en(aSx −aSt ) Hi{n} , P = e2(aSx −aSt ) P , (6.3)

have their foundation in the NSEs (2.23) and (2.24) in the limit of Re → ∞ and are transferred
to the H-approach.

Note that the NSEs only admit a reduced set of scaling symmetries compared to the Euler
equations, i.e., aSt = 2aSx , in (6.3) to be readily shown by substituting the scaling symmetries
(6.3) into the momentum balance equations (2.23). However, for turbulent flows at larger
Reynolds numbers, viscosity is only dominant in the range of the Kolmogorov length ηk . Based
on this, a boundary-layer type of asymptotic expansion has been conducted in Oberlack (2000b)
showing that for length scales beyond ηk , the scaling is defined by the Euler equations.

60
A third scaling symmetry of the MPMEs (3.30)
∗
T Ss : t∗ = t, x∗i = xi , Ui = eaSs Ui , Hi∗{n} = eaSs Hi{n} (6.4)

defines a measure of intermittency (Wacławczyk et al., 2014) and is of statistical nature only,
i.e., it has no counterpart in the NSEs as pointed out above.
These three symmetries form the basis for the derivation of turbulent scaling laws or invariant
solutions for turbulent round jet flows after consideration of certain invariants to be shown
below.

6.1.2. Derivation of invariant solutions

Subsequently, we focus on one-point statistics, i.e., x(n) = x, and in this case Equation (3.25)
yields
Hi{n} = U[i] n, (6.5)
where ( · )[·] means that no summation of the quantity ( · ) is implied. One of the reasons is the
comparability with the DNS data where only statistics at one point have been extracted.
Also, we transition to a cylindrical coordinate system, specified by the coordinates
i = r, ϕ, z (6.6)
for which the transformation of the variables is straightforward and can be done through
multiplication with a corresponding transformation matrix. This has no effect on the scaling
properties of the moments in the symmetries (6.3) and (6.4) since the group parameters appear
as prefactors of the moments and are therefore independent of this coordinate transformation.
To elaborate on this, a corresponding matrix A with the associated metric terms is considered for
the transformation between the two systems. If A is the Cartesian-to-cylindrical transformation
matrix, then
Hc∗
{n} = AH{n} ,
∗
Hc{n} = AH{n} (6.7)
are the moments in a cylindrical coordinate system. Most important, this has no influence
on the scaling properties of the moments. For the scaling of Hi{n} due to, e.g., the scaling
symmetry (6.3), we have
Hi∗{n} = en(aSx −aSt ) Hi{n} , Hi∗{n} = eaSs Hi{n} . (6.8)

Therefore multiplying (6.8) with A yields

Hic∗
{n}
= en(aSx −aSt ) Hic{n} , Hic∗
{n}
= eaSs Hic{n} , (6.9)

hence, invariant surface conditions (see Equation (4.22)) for the Cartesian and cylindrical
systems are equivalent.
With the symmetries (6.2)-(6.4) the invariant surface condition (4.22) for invariant solutions
of a spatially developing turbulent round jet can be generated
dr dz dUr Uz
= = (6.10)
aSx r aSx z + az [2(aSx − aSt ) + aSs ]Ur Uz
dU[i] dU[i]
n
= = ... = n
.
[(aSx − aSt ) + aSs ]U[i] [n(aSx − aSt ) + aSs ]U[i]

61
Initially, the group parameters aSx , aSt , aSs , az are arbitrary. However, invariant solutions
must satisfy the properties of a given flow, in our case the invariant of a turbulent round jet
flow. In most cases, invariants cause a symmetry breaking, i.e., they impose a constraint on
the group parameters, which helps to understand the significance of these parameters. The
invariant of a turbulent round jet flow is the conservation of momentum, which states that the
momentum introduced by the inlet is constant for any distance from the orifice. This invariant
can be derived from the RANS equations (3.16) in cylindrical coordinates for a axisymmetric
jet
1 ∂(rUr2 ) ∂(Ur Uz ) u2ϕ ∂P
+ − =− , (6.11)
r ∂r ∂z r ∂r
1 ∂(rUr Uz ) ∂(Uz2 ) ∂P
+ =− , (6.12)
r ∂r ∂z ∂z
where the equation of continuity (3.17) has been inserted and the viscousR ∞ terms have been
already omitted. Equation (6.11) can be integrated on both sides with r ( · )dr0 to receive an
explicit expression for the mean pressure P
Z ∞" 2 2
#
(U ) ∂(Ur Uz ) u ϕ
P = P0 − Ur2 + r
+ − dr0 , (6.13)
r r ∂z r

where P0 is the total pressure. In the following, the pressure gradient from Equation (6.12) is
eliminated by employing Equation (6.13). By integrating across the entire jet and considering
that Ur and the turbulence quantities approach zero as r → ∞, a momentum integral invariant
I for turbulent round jets is obtained (see e.g. Hussein et al., 1994)

1 d
Z ∞ Z ∞
 1 ∞ 2
Z
1 2
2
Uz − u + uϕ r dr − I = −
2 Ur Uz r dr +
2
Ur r dz. (6.14)
0 2 r 2 dz 0 2 0

In Hussein et al. (1994), they argue that the first term on the right side of (6.14) is less than
0.02I and the second term on the right side is an order of magnitude smaller than the first
term on the right side. By neglecting these two terms, the invariant from Hussein et al. (1994)
yields Z ∞
1 2 
IG = 2
Uz − u + uϕ r dr.
2 (6.15)
0 2 r
However, by preserving the complete momentum flux in (6.14) and by considering

u2ϕ = Uϕ2 , (6.16)

since Uϕ = 0, we receive
Z ∞
1 2  1 d 
Ur Uz r r dr, (6.17)


IO = Uz2 − 2
Ur + Uϕ +
0 2 2 dz

where IO refers to the invariant used in this dissertation. IO is used because it describes the
complete form of the momentum conservation compared to IG (6.15). Implementing the
symmetries (6.3)-(6.4) into Equation (6.17) through
∗ ∗
r = e−aSx r∗ , Ui2 = e−2(aSx −aSt )−aSs Ui2 , Ur Uz = e−2(aSx −aSt )−aSs Ur Uz (6.18)

62
gives
Z ∞
∗ 1  2∗ ∗
 1 d  
∗ ∗
IO = e−4aSx +2aSt −aSs Uz2 − Ur + Uϕ2 + U U
r z r r∗ dr∗ , (6.19)
0 2 2 dz ∗

where IO has to be invariant and therefore implies the symmetry breaking

1
aSt = 2aSx + aSs . (6.20)
2
On the contrary, the invariant IG (6.15) leads to two symmetry breakings

aSt = 2aSx , aSs = 0. (6.21)

This contrasts with the classical integral invariant (Schlichting and Gersten, 2017)
Z ∞
2
IC = Uz r dr, (6.22)
0

which assumes an axially vanishing pressure gradient and neglects Reynolds stresses. Not
neglecting the Reynolds stress, i.e., u2z 6= 0, results in
Z ∞
IC =
0 Uz2 r dr (6.23)
0

and notably, this invariant leads to an identical symmetry breaking as in IO (6.17). Symmetry
theory reveals no difference between the invariants IC 0 and IO , with IO being less restrictive
by formally allowing a pressure gradient in the axial direction.
Integrating the invariant surface condition (6.10) and considering the symmetry breaking
(6.20) through IO (6.17) gives the invariants, which form the similarity variable
r
η= (6.24)
z − z0

and the turbulent scaling laws

U
e (η)
z
U z (r, z) = 1 ∗ ,
(z − z0 )1− 2 aSs
U]
r Uz (η)
Ur Uz (r, z) = ,
(z − z0 )2
Ufn (η)eci,n n
[i]
n (r, z) =
U[i] 1 ∗ ∗
, (6.25)
(z − z0 )n(1+ 2 aSs )−aSs

where
az aSs
z0 = − , a∗Ss = (6.26)
aSx aSx
and (f· ) indicates the similarity variables or invariants. The constant of integration, i.e.,
invariant η in Equation (6.24) is obtained from integrating the first two terms, while turbulent
scaling laws (6.25) result from integrating the second and remaining terms.

63
The integration constants ci,n represent the invariants related to the moments and get into the
exponent after exponentiation in Equation (6.25). The invariants ci,n are in general functions
of other invariants and presently this is η. Remarkably, the DNS data prove that the ci,n do not
depend on n (see Figure 6.1), representing an invariant with respect to the moment order n.
This allows factoring out the η-dependence as seen in Equation (6.25).

The Equations (6.24) and (6.25) satisfy the symmetry breaking (6.20) of both moment integrals
IO (6.17) and IC 0 (6.23). Additionally, a symmetry reduction of the MPMEs (3.30) (see.
Section 6.3.1), IO (6.17) and IC 0 (6.23) is possible. Further, for n = 2, the H-moments are
unaffected by a∗Ss , representing an invariant regarding a variation of the inflow condition and
therefore generating a new hypothesis for further turbulent jet studies.

The invariant (6.17) induces a one-parameter family of a∗Ss -dependent scaling laws in (6.25),
although DNS data for both Reynolds number show

a∗Ss = aSs = 0. (6.27)

Despite this, aSs , i.e., the symmetry (6.4) denoting intermittency, is crucial for near-wall scaling
laws (Oberlack et al., 2022) and allows variation in inflow parameters and turbulent decay
behavior as suggested by George (1989), and is not influenced by the Reynolds number for
turbulent pipe inlets. The intermittency is presently hidden through (6.27) but reappears in
subsequent reduced equations.

6.2. Validation of the turbulent velocity scaling laws on the centerline

To validate the scaling laws (6.24) and (6.25), we utilize the DNS data set at Re = 3500 and
Re = 7000 extensively discussed in Chapter 5. Further data of velocity moments up to the
tenth order, which have not been covered there, are used for the following validation. For
validation, we first set r = 0 in (6.25)

[i] (η = 0)αi,n e
Ufn ci n
n = U n (r = 0, z) =
U[i],c [i] 1 ∗ ∗
, (6.28)
(z − z0 )n(1+ 2 aSs )−aSs

where αi,n is extracted from the invariant such that

Ufn
[i] (η = 0) = 1 (6.29)

and a∗Ss = 0 as mentioned in Equation (6.27). From Table 5.3, we can extract the virtual
orifices for the two Reynolds numbers.

In the next subsection, the constants αi,n and ci are fitted to the DNS data of the axial velocity
moments at the two Reynolds numbers. Subsequently, they are fitted to DNS data of the radial
and azimuthal velocity moments.

64
6.2.1. Results for the axial velocity scaling law

Figure 6.1 shows the exponential prefactors αi,n eci n for both Reynolds numbers. The data of
the DNS at Re = 3500 follows the exponential prefactors αz ecz n in (6.28) with

cz = 2.123, αz,n = αz = 0.673 (6.30)

quite accurately and are rather similar to the scaling laws of a turbulent channel flow (Oberlack
et al., 2022). The DNS data of the Re = 7000 case gives the constants

cz = 2.127, αz = 0.662 (6.31)

for Equation (6.28). The values of cz and αz are very similar, highlighted by a strong overlap,
which suggests an independence from Reynolds numbers. However, this should be further
investigated by conducting additional simulations with a wide range of Reynolds numbers.

6.2.2. Results for the radial and azimuthal velocity scaling laws

For the radial and azimuthal moments the behavior is quite different. On the axis r = 0, the
radial and azimuthal mean velocities are zero, i.e.,

Hi{n} (r = 0, z) = Ri{n} (r = 0, z). (6.32)

The radial and azimuthal moments can therefore be modeled as a centered Gaussian process
in a very good approximation. This means all the higher moments n are determined by the
second moment (see the definition in (3.10)) with

n
αr,n ecr n = αϕ,n ecϕ n = αr,2 e2cr 2
(n − 1)!!, (6.33)

where
αr,2 e2cr = 1.81 (Re = 3500), αr,2 e2cr = 1.87 (Re = 7000), (6.34)
can be taken from Figure 6.1. ( · )!! refers to the double factorial, described in Equation (3.11),
and αr,n = 0 for odd n. With Stirling’s formula for factorials (Tweddle, 2003) we obtain

n √
cr,n = cr = −0.5, αr,n = αr,2 e2cr 2
2e(n − 1)n/2 . (6.35)

Therefore for all velocity moments of both Reynolds numbers, the scaling is quite similar as
depicted in Figure 6.1.

6.3. Gaussian behavior of the mean axial velocity

The present analysis reveals a Gaussian-type behavior in η not only in the mean U
e (η) but also
z
for all higher instantaneous axial velocity moments. Based on the scaled higher axial moments
in (6.25), we find
Uzn fn (η) = e−γn η2
n
=U z (6.36)
Uz,c

65
 1010
109
108
107
106
αi,n eci n

105
104
103
102
101
100
10−1
0 1 2 3 4 5 6 7 8 9 10
n
Figure 6.1.: The exponential prefactors αi,n eci n of the DNS according to the velocity scaling law on the centerline
(6.28) are shown for each moment up to order n = 10. For Re = 3500, they are presented for each
direction i = r ( ), z ( ) and determined with the DNS data ( ). Likewise, for Re = 7000, they are
displayed for each direction i = r ( ), z ( ) and determined with the DNS data ( ).

being an excellent representation of all Uz moments. At first, this was an empirical observation
but subsequently, it is cast into the group theoretical framework.

The self-preservation of moments up to order ten at five stages between z/D = 25 − 65 can
be taken from Figure 6.2 for Re = 3500 and in Figure 6.3 for Re = 7000. They exhibit an
almost perfect collapse of the data by normalizing with the scaling law (6.28) and become
increasingly narrower as the moment order increases. The remaining four moments for each
of the two Reynolds numbers can be depicted in Appendix A.2. Focusing on the exponent
of Equation (6.36) alone is a very sensitive measure for the data. In Figure 6.4, deviations of
the DNS data at both Reynolds numbers from (6.36) are only visible at the edge of the jet in a
semi-logarithmic representation due to insufficient statistics of very small values. A quadratic
fit of γn

γn = −1.54n2 + 61.95n + 27.10 (Re = 3500), γn = −1.51n2 + 63.39n + 28.21 (Re = 7000)
(6.37)
is shown in Figure 6.5 for both Reynolds numbers.

An interesting analogy arises with the moments of the structure function for small-scale turbu-
lence (Frisch and Kolmogorov, 1995; Sreenivasan and Antonia, 1997) where the experimental
data deviate from pure dimensional scaling, related to Kolmogorov’s linear scaling of the
exponent in n. The non-linearity of the exponent in n leads to heavy-tails of the PDF, i.e.,
rare events in a Gaussian are more likely, which is attributed to internal intermittency (Nelkin,
1994). Similarly, the non-linear behavior of the Gaussian curves (6.37) in n within γn strongly
contrasts with the perfect dimensional scaling of all moments in the scaling laws (6.25) where

66
 1 1

0.8 0.8
Uzn /Uz,c

0.6 0.6
n

0.4 0.4

0.2 0.2

0 0
0 0.05 0.1 0.15 0.2 0.25 0 0.05 0.1 0.15 0.2 0.25
η η
(a) n = 1 (b) n = 2
1 1

0.8 0.8
Uzn /Uz,c

0.6 0.6
n

0.4 0.4

0.2 0.2

0 0
0 0.05 0.1 0.15 0.2 0.25 0 0.05 0.1 0.15 0.2 0.25
η η
(c) n = 4 (d) n = 6
1 1

0.8 0.8
Uzn /Uz,c

0.6 0.6
n

0.4 0.4

0.2 0.2

0 0
0 0.05 0.1 0.15 0.2 0.25 0 0.05 0.1 0.15 0.2 0.25
η η
(e) n = 8 (f) n = 10

Figure 6.2.: The radial profiles of the nth axial moment normalized with the scalings in (6.25) at different distances
from the orifice for the Re = 3500 case: z/D = 25 ( ), 35 ( ), 45 ( ), 55 ( ), 65 ( ).
The black solid lines indicate the Gaussian behavior from (6.36) using γn from (6.37).

67
 1 1

0.8 0.8
Uzn /Uz,c

0.6 0.6
n

0.4 0.4

0.2 0.2

0 0
0 0.05 0.1 0.15 0.2 0.25 0 0.05 0.1 0.15 0.2 0.25
η η
(a) n = 1 (b) n = 2
1 1

0.8 0.8
Uzn /Uz,c

0.6 0.6
n

0.4 0.4

0.2 0.2

0 0
0 0.05 0.1 0.15 0.2 0.25 0 0.05 0.1 0.15 0.2 0.25
η η
(c) n = 4 (d) n = 6
1 1

0.8 0.8
Uzn /Uz,c

0.6 0.6
n

0.4 0.4

0.2 0.2

0 0
0 0.05 0.1 0.15 0.2 0.25 0 0.05 0.1 0.15 0.2 0.25
η η
(e) n = 8 (f) n = 10

Figure 6.3.: The radial profiles of the nth axial moment normalized with the scalings in (6.25) at different distances
from the orifice for the Re = 7000 case: z/D = 25 ( ), 35 ( ), 45 ( ), 55 ( ), 65 ( ).
The black solid lines indicate the Gaussian behavior from (6.36) using γn from (6.37).

68
 100 100

10−3 10−3

10−6 10−6
Uzn
10−9 10−9

10−12 10−12

10−15 10−15
0 0.05 0.1 0.15 0.2 0 0.05 0.1 0.15 0.2
η η

Figure 6.4.: The radial profiles (blue, dashed) of the 1st (top) up to the 10th (bottom) axial moment and the
corresponding curves (black) from (6.36) shown in a semi-logarithmic plot at z/D = 45 for the
Re = 3500 case (left) and the Re = 7000 case (right).

n (r, z) ∝ z −n due to a∗ = 0. Therefore, the non-linearity in n of γ gives rise to intermit-

U[i] Ss n
tency. This is supported by the shapes of the PDF and the curves of skewness and kurtosis of
the fluctuations in Section 5.6 clearly show “large-scale intermittency”.

The distinctly non-linear behavior in (6.37) shows that higher-order moments cannot be
ignored and simultaneously raises the question of their determination based on first principle.
Subsequently, we will see that there is a type of “hidden intermittency” being a combination of
two statistical symmetries resulting from the reduced system in the η-Ufn variables.
[i]

6.3.1. Symmetry reduction of the multi-point moment equations

A symmetry reduction can be applied to the MPMEs (3.30), i.e., the set of variables of (3.30)
can be reduced, so that an ordinary differential equation (ODE) is obtained. This allows us
to extract the Gaussian behavior (6.36). Consider the moment equation in z-direction for a
statistically stationary jet, i.e., the MPMEs (3.30) at one point, neglecting the pressure term, in
polar coordinates independent of t and ϕ

1 ∂(rUr Uzn−1 ) ∂(Uzn )

+ =0 (6.38)
r ∂r ∂z

with n = 1, 2, where n = 1 yields the continuity equation. Transforming (6.38) by inserting

the similarity variable (6.24) and the scaling laws (6.25) gives the ODE

^
^ dUr Uzn−1 dUfn
n−1
Ur Uz + η − η 2 z − nη U
fn = 0, (6.39)
dη dη z

69
 600
Re = 3500
Re = 3500 fit
500 Re = 7000
Re = 7000 fit
400

γn 300

200

100

0
1 2 3 4 5 6 7 8 9 10
n
Figure 6.5.: Constants γn from (6.36) are shown for each moment up to order n = 10 determined by fitting to
the DNS yielding γn = −1.54n2 + 61.95n + 27.10 (Re = 3500) and γn = −1.51n2 + 63.39n + 28.21
(Re = 7000).

where all moments are functions of η only. In turn, from Equation (6.39), we may derive an
^
equation for Ur Uzn−1 by employing the Gaussian behavior. Using (6.36) in (6.38), we receive

^ 1 h 2
i
Ur Uzn−1 = 2γn η 2 − n + 2 e−γn η + (n − 2) , (6.40)


2γn η

^
where the boundary condition limη→0 Ur Uzn−1 = 0 has been used. Comparing (6.40) and the
DNS data in Figure 6.6 for n = 1, 2 at Re = 3500, we observe a remarkable collapse for Ur and
Ur Uz , which typically converge poorly.

6.3.2. Hidden intermittency

For understanding the symmetry, and hence the intermittency behavior of a turbulent jets,
let us recapitulate three generic symmetries of linear equations, particularly of the moment
equation (3.30)
 
n
" #
∂Hi{n} X ∂Hi{n+1} i(n) 7→k(l)  x(n) 7→ x(l) ∂Ii{n−1} [l] 2
1 ∂ Hi{n}
+ + − = 0. (6.41)
∂t ∂xk(l) ∂xi(l) Re ∂xk(l) ∂xk(l)
l=1

However, we will focus mainly on the reduced moment equations (6.39) as an example:
1. Linear equations such as the MPMEs (6.41) always possess the scaling of the dependent
variables as given in the statistical symmetry (6.4).

70
 ·10−2
4

2 n=2
Ur Uzn−1 /Uz,c
n

n=1
−1

−2

−3
0 0.05 0.1 0.15 0.2 0.25 0.3
η

Figure 6.6.: The radial profiles of Ur Uzn−1 for the Re = 3500 at different distances z from the orifice compared
to the solution in (6.40) ( ): z = 25 ( ), z = 35 ( ), z = 45 ( ), z = 55 ( ), z = 65
( ).

2. When given in conservation form such as the MPMEs (6.41), linear equations possess a
shift in the dependent variables, i.e., the addition of a constant (Wacławczyk et al., 2014)
though presently not required.

3. Linear equations always possess a third generic symmetry, i.e., the superposition principle,
which we will present in (6.44) for the reduced moment equations.

The reduced moment equations (6.39) admit the statistical symmetry

∗
T Ss ∗ fn ∗ = eaSs U
: η = η, U fn , U^ n−1 ^
= eaSs Ur Uzn−1 , (6.42)
z z r Uz

which apparently is a reminiscence of the statistical symmetry (6.4) and constitutes a measure
of intermittency (Wacławczyk et al., 2014). The second scaling symmetry has been inherited
from the two-parameter group aSx , aSt from the scaling symmetry (6.3)
∗
fn ∗ = U
T S : η ∗ = eaS η, U fn , U^ n−1 ^
= eaS Ur Uzn−1 (6.43)
z z r Uz

and apparently has been reduced to a one-parameter group aS by the symmetry reduction
through the similarity variable (6.24) and the invariant solutions (6.25). Classical superposition
principle is also a statistical symmetry, as the NSEs admit no correspondence, and for the
reduced Equation (6.39) it reads as follows
∗ 0
TS : η ∗ = η, fn ∗ = U
U fn 0 ,
fn + U ^ ^ ^
Ur Uzn−1 = Ur Uzn−1 + Ur Uzn−1 , (6.44)
z z z

71
where the prime quantities are any additional independent solution of the original linear
system, here of Equation (6.39).
In order to obtain the link to the Gaussian functions (6.36) of the Uz -moments, we first want
to generate an invariant basis solution from the first two symmetries (6.42) and (6.43). This
is analogous to the invariant surface condition (6.10) and reads as follows

^
dη dU
fn dUr Uzn−1
= z
= . (6.45)
as η fn
aSs U ^ n−1
z (as + aSs )Ur Uz
The integration of the first two terms immediately yields a power-law with an initially arbitrary
exponent q = aSs /aS
fn = C η q .
U z n (6.46)
In order to unravel the the link to (6.36), we rewrite it as a converging Taylor series
∞
−γn η 2
X (−γn )k 1
2 1
3
e = η 2k = 1 − γn η 2 + γn η 2 − γn η 2 + ..., (6.47)
k! 2 6
k=0

which proves that, using the superposition symmetry (6.44), the solutions (6.46) are the basis
elements of the Gaussian behavior (6.36) with
∞
2
X
e−γn η = Cnk η qk , (6.48)
k=0

where
(−γn )k
Cnk = , qk = 2k. (6.49)
k!
Likewise, the importance of the intermittency symmetry (6.42) in the exponent qk becomes
clear.
This shows that the intermittency symmetry in the combination of (6.42)-(6.44) is a key to
the understanding the intermittency behavior. This discloses the “hidden intermittency” due to
a certain combination of statistical symmetries, where extensive use of symmetry (6.44) was
made, reflecting the superposition principle.
The above solutions (6.36) and (6.40) are a combination of statistical symmetries, which may
directly be verified by implementing them into (6.39). In other words, (6.36) and (6.40) are
described by the key symmetries (6.42)-(6.44) for the “large-scale” intermittency properties of
a turbulent jet. Simultaneously, it represents the superposition principle of statistical moments
as defined in Equation (6.44). Interestingly, the Gaussian behavior (6.36) and the invariant
solution (6.40) show a quite intricate structure though are composed of rather elementary
symmetries.

6.4. Conclusive remarks

In this chapter, we have shown that symmetry methods yield turbulent scaling laws generated
by statistical turbulence theory that go beyond simple power-laws, logarithmic or exponential

72
laws. These scaling laws have been verified against data that has been generated by the
large-scale simulations in Chapter 5.
Once again, statistical symmetries originating from linear properties of the infinite hierarchy
of moment equations and having no correspondence in the NSEs, form the basis for the scaling
laws, which are the key to a possible variation of the turbulent decay.
Beyond the statistical scaling symmetry, we introduce, for the first time in turbulence, the
repetitive use of the superposition principle, representing a symmetry of all linear differential
equations (Bluman et al., 2010) and serving as an elementary basis here. This methodology
allows the derivation of fundamental scaling laws for various turbulent canonical flows, where
these scaling laws are solutions of the infinite moment hierarchy. While symmetry methods
enable the derivation of more complex scaling laws, determining the scaling parameters remains
an open question, especially regarding γn in (6.37) as a function of n.
Notably, an analogy emerges with small-scale turbulence structure function moments (Frisch
and Kolmogorov, 1995; Sreenivasan and Antonia, 1997), where internal intermittency (Nelkin,
1994) causes deviations from pure dimensional scaling, related to Kolmogorov’s linear scaling
in n.
Moreover, the validation of scaling laws through DNS at two distinct Reynolds numbers has
revealed minimal changes in the scaling law parameters, except for the virtual origin, which
could be attributed to inherent variance. However, to confirm this result, additional simulations
or experiments at higher Reynolds numbers with similar conditions must be performed to
verify this statement, since the parameters may exhibit a small dependence on the Reynolds
number. Furthermore, additional research employing Lie symmetry analysis could be pursued
to analytically discover this dependence.
In the next chapter, the focus is on passive scalars in a turbulent round jet flow, where again
the extensive DNS data set of the passive scalar is discussed and the turbulent scaling laws are
extended to passive scalar and mixed velocity-scalar moments.

73
7. Passive Scalar Statistics in a Turbulent
Round Jet Flow

This chapter is heavily based on the publication Nguyen and Oberlack (2024d) and represents
the main results of the analysis of the statistics of a passive scalar in a turbulent round jet flow
and a subsequent Lie symmetry analysis which unifies and therefore generalizes the turbulent
scaling laws for the velocity and scalar moments.

To validate scaling laws, especially of higher order moments, high-fidelity statistics are impera-
tive, which can be achieved by DNS. However, there have been limited DNS studies for turbulent
round jet flows so far, especially those with a long box and incorporating an additional passive
scalar. Passive scalars are scalar quantities that are not actively involved in the flow physics
such as the temperature or a concentration of a fluid. The Reynolds number in these studies
has often been kept low, and the averaging process has been limited to small time frames due
to computational constraints.

Notably, Boersma et al. (1998) initiated DNS for turbulent round jet flows as mentioned in
Chapter 5. As an extension of this work, Lubbers et al. (2001) examined the self-similarity
of a passive scalar concentration at Re = 2000 and a Schmidt number Sc = 1 in a box with
the length z/D = 40. The statistics have been extracted over 80D/Ub time units where Ub
refers to the bulk velocity at the inlet. The results show that the mean concentration in the far
field is self-similar. However, the root mean square of the concentration fluctuations are not
self-similar.

In Babu and Mahesh (2005) a DNS of a turbulent jet at Re = 2400 and Sc = 1 is performed. The
data are averaged over 1400D/Ub time units. The instantaneous radial profiles of the velocity
and passive scalar exhibit similarities and alternate between “top-hat” and “triangle” profiles,
both spatially and temporally. These profiles are influenced by entrainment from the free-
stream into the jet, resulting in a mean Gaussian profile as a function of r. Diffusion-dominated
regions are observed, occurring closer to the jet center and as “brush-like” regions near the
jet edge. The width of these brush-like regions decreases with increasing Re, suggesting a
transition in mixing behavior.

Simulations by Gilliland et al. (2012) has focused on a DNS at Re = 2400 and a LES at
Re = 68,000 in a box length of z/D = 50 to investigate scalar intermittency in a turbulent
round jet. This simulation is quite recent and the averaging for both simulations has been
performed over about 900D/Ub time units. The results emphasize the importance of external
intermittency in scalar mixing. The study highlights the need for improved subgrid scale
models in LES to enhance the accuracy of predicting external intermittency over a wide range
of Reynolds numbers.

75
However, considering the computational cost of DNS led many researchers to opt for experimen-
tal approaches. In a notable study, Birch et al. (1978) has conducted experiments measuring the
turbulent concentration parameters of a free round methane jet up to the fourth moment. Their
results have revealed a departure from Gaussianity in the mean passive scalar concentration
along the centerline, as indicated by the negative skewness of the PDF.
Dowling and Dimotakis (1990) present an experimental investigation of the turbulent con-
centration field formed when a free turbulent jet at Re = 5000, 16,000, 40,000 mixes with
gas entrained from a quiescent reservoir at a Prandtl numbers between P r = 1 − 1.2. Laser-
Rayleigh scattering measurements taken between z/D = 20 − 90 show a nearly independent
behavior with respect to the Reynolds number near the centerline of the jet for the scaled PDF
of the jet fluid concentration.
Antonia and Mi (1993) have measured the average temperature dissipation using parallel
cold wires at Re = 19,000 and a Péclet number of P e = 83 at a distance of z/D = 30 from
the orifice. The resulting components of the average temperature dissipation, particularly the
radial and azimuthal values, were found to be nearly equal, only slightly larger than the axial
component. The deviation from isotropy of the temperature dissipation was small, especially
when compared to results in other free shear flows.
In the work of Panchapakesan and Lumley (1993b), flying hot-wire measurements of a helium
jet were made up to third-order moments, including mixed moments, measured between
z/D = 50 − 120. The measurements, consistent with earlier studies of helium jets, provided
insights into the mean velocity field, the concentration decay constant, and the radial profiles
of mean velocity and mean concentration.
More recently, Darisse et al. (2015) studied a slightly heated turbulent round air jet at Re =
14 · 103 using temperature as a passive scalar. Laser doppler velocimetry and laser doppler
velocimetry-cold-wire thermometry measurements were performed for the variance, third-order
moments and mixed moments at z/D = 30. This allows the investigation of all passive scalar
transport budget terms except the dissipation term, which is derived by closing the balance.
While previous studies, such as Birch et al. (1978), have ventured into the investigation of
moments up to the fourth order in turbulent round jets, the investigation of moments beyond
the third order has been relatively uncharted territory.
In Chapter 6, we applied Lie symmetry analysis to a spatially evolving turbulent round jet
flow, where we were able to derive turbulent velocity scaling laws up to an arbitrary order and
validate them against DNS data of velocity moments up to the tenth order. Furthermore, we
found that the effects of intermittency are hidden in the symmetries of the governing equations
highlighting the importance of Lie symmetry analysis in improving our understanding of
turbulence.
In this chapter, DNS data on passive scalar moments up to the tenth order and mixed velocity-
scalar moments up to the sixth order shall be presented. Further, the topics discussed in
Chapter 6 shall be expanded onto the passive scalar and mixed moments for a round jet flow.
A detailed description of the DNS methodology can be found in Chapter 5. In addition to that,
a passive scalar has been solved using Equation (2.44) for air at P r = 0.71, and is introduced
at the inlet as a top-hat profile with Θ = 1. At z/D = 0 outside the jet inlet, the passive scalar
is set to zero. To avoid non-physical behavior, the thermal open BC, as presented in Liu et al.
(2020), is implemented for the passive scalar at the lateral boundaries and at the jet outlet at

76
 Table 7.1.: BCs of the main computational domain for a passive scalar
Region Passive Scalar BC

{(r, ϕ, 0)T | r ∈ [0, 0.5], ϕ ∈ [0, 2π)} Dirichlet BC

{(r, ϕ, 0)T | r ∈ (0.5, 2), ϕ ∈ [0, 2π)} Dirichlet BC
{(0.4z + 2, ϕ, z)T | ϕ ∈ [0, 2π), z ∈ [0, 75]} Thermal open BC
{(r, ϕ, 75)T | r ∈ [0, 32], ϕ ∈ [0, 2π)} Thermal open BC

z/D = 75 as detailed in Table 7.1. This BC is energy-stable, and ensures that the contribution
of the open BC will not cause the passive scalar to increase over time.

7.1. Passive scalar statistics

In this section, we investigate the statistical properties associated with a passive scalar in
turbulent round jet flows. Leveraging insights from the DNS, we explore various facets of
the passive scalar field, ranging from mean profiles to higher-order statistical moments and
PDFs. Through rigorous comparisons with previous experimental results, we aim to outline
both agreements and differences, providing a comprehensive view of the dynamics governing
turbulent round jet flows.

7.1.1. Mean passive scalar statistics

The Reynolds decomposition (3.12) is equally applicable to a passive scalar

Θ = Θ + θ. (7.1)

Analogously, the decomposition of a passive scalar extends to

2
HΘΘ = Θ + RΘΘ (7.2)

and mixed moments are denoted by

HiΘ = Ui Θ + RiΘ , (7.3)

where RΘΘ = θ2 and RiΘ = ui θ. The mean centerline passive scalar Θc is scaled by the mean
centerline passive scalar at the orifice Θ0 , which is inversely related to the distance from the
orifice
Θc (z) BΘ D
= . (7.4)
Θ0 z − z0
From the present DNS at Re = 3500, we deduce a decay constant of BΘ = 4.9 with a virtual
origin at z0 = −0.132 which means the virtual origin is close to the virtual origin of the velocity
z0 = 0 from classical mean velocity scaling (5.4). For the DNS at Re = 7000, the decay constant
BΘ = 5.0 at the virtual origin z0 = 2.7 is obtained. The virtual origin shifts to a larger value for
increasing Re just like for pure hydrodynamics (see Table 5.3). The parameters are compared
to Lubbers et al. (2001), where they found z0 = 0.5 for the concentration and z0 = 5.5 for

77
the axial velocity. Figure 7.1 showcases a consistent and smoother linear trend of the inverse
passive scalar up to z/D = 75 for both simulations compared to other experiments and DNS.
In the magnified view in Figure 7.1, the onset of decay can be exhibited. Similar to a potential
core of the mean axial centerline velocity in Section 5.2, we have highlighted the limit where
the passive scalar reaches 95% of its maximum, i.e., where the inverse reaches 1/0.95 = 1.053.
It can be observed that both present simulations reach the decay zone earlier than the other
studies. This has also been observed for the mean inverse axial centerline velocity which is
due to the turbulent pipe velocity profile causing an earlier mixing and therefore an earlier
breakdown of the potential core.
For the mean axial velocity, the decay constant is Bu = 5.15 for the Re = 3500 case and
Bu = 5.25 for the Re = 7000 case as seen in Table 5.3, implying that Θc decays faster than
the mean axial centerline velocity U z,c in both cases. Lubbers et al. (2001) have observed the
same, but their virtual origins for the velocity and the passive scalar were further apart.

15
Θ0 /Θc

10

1.1

1.05
5
1

0.95
2 4 6 8 10

0
0 10 20 30 40 50 60 70
z/D

Figure 7.1.: Mean inverse passive scalar at the centerline over the distance of the orifice: Birch et al. (1978) at
Re = 16,000 ( ), Babu and Mahesh (2005) at Re = 2400 ( ), Lubbers et al. (2001) at Re = 2000
( ), present DNS at Re = 3500 ( ), present DNS at Re = 7000 ( ). In the magnified view,
the potenial core of the passive scalar is marked by a dashed line ( ).

In Figure 7.2, Θc reaches self-similarity at z/D = 15 for the Re = 3500 case while it reaches
self-similarity at z/D = 20 for the Re = 7000 case. For the lower Reynolds number, this
onset occurs earlier than the assumed z/D = 20 in Lubbers et al. (2001). Observing self-
similarity at a later stage for increasing Reynolds number was also seen in the DNS data of
pure hydrodynamics in Chapter 5. Upon comparing the profile of U z,c with that of Θc , the
latter shows a slightly larger spreading rate. The half-width of the mean passive scalar for both
Reynolds numbers, η1/2,Θ = 0.109, exceeds that of the velocity, η1/2,u = 0.089 at Re = 3500
and η1/2,u = 0.086 at Re = 7000, in Section 5.2. Remarkably, η1/2,Θ closely aligns with the
observation in Lubbers et al. (2001), where in their DNS the half-width is η1/2,Θ = 0.108.
In summary, the characteristics of Θc are consistent with those observed by Lubbers et al.
(2001). The profiles exhibit self-similarity and a faster spread compared to the velocity profile.
Moreover, the virtual origins of the scalar and velocity appear to have shifted closer to each

78
 1.2 1.2

1 1

0.8 0.8
Θ/Θc

0.6 0.6

0.4 0.4

0.2 0.2

0 0
0 0.05 0.1 0.15 0.2 0.25 0.3 0 0.05 0.1 0.15 0.2 0.25 0.3
η = r/(z − z0 ) η = r/(z − z0 )

Figure 7.2.: Mean passive scalar Θ scaled with the centerline passive scalar Θc (7.4) for the Re = 3500 case (left)
plotted as a function of the similarity coordinate η at different distances from the orifice: z/D = 15
( ), 25 ( ), 35 ( ), 45 ( ), 55 ( ). For the Re = 7000 case (right), Θ/Θc is plotted at:
z/D = 20 ( ), 30 ( ), 40 ( ), 50 ( ), 60 ( ).

other compared to the DNS by Lubbers et al. (2001), which may be attributed to the turbulent
pipe inlet. The Reynolds number only has a strong influence on the virtual origin out of all the
parameters discussed in this section.

7.1.2. Second order passive scalar and mixed statistics

The variance of the passive scalar fluctuation RΘΘ is shown in Figure 7.3 for different distances
z/D from the orifice. The values at η = 0 being RΘΘ /Θ2c = 0.048 at Re = 3500 and RΘΘ /Θ2c =
0.05 at Re = 7000 agree with the results of Darisse et al. (2015) and Lubbers et al. (2001),
both of which report values around 0.04 and slightly above. However, the off-axis peak being
RΘΘ,max /Θc = 0.067 exceeds the value of RΘΘ,max /Θc = 0.05 measured by Darisse et al.
(2015). Since the present simulations exhibit the same value, this could be due to the different
inlet condition compared to the top-hat velocity profile in the experiment by Darisse et al.
(2015).

Lubbers
√ et al. (2001) notes in Figure 7.4 that the normalized rms of the passive scalar fluctuation
RΘΘ shows an increasing trend rather than a horizontal line, the former indicating non-
self-similarity. This trend can be attributed to the smaller axial length of their computational
box, as the present
√ DNS at Re = 3500 shows a slight increase up to z/D = 35. After this
point, however, RΘΘ /Θ √ c stabilizes around 0.22, supporting the conclusion by Dowling and
Dimotakis (1990), that RΘΘ lies between 0.23 and 0.24, with the present DNS slightly below
this range. Interestingly, the initial maximum close at z/D = 6 of the Re = 3500 case is
smoothed out for the Re = 7000 case and follows the experiment by Birch et al. (1978) at
Re = 16,000 almost perfectly before settling at a similar value of 0.22 as the Re = 3500 case.

79
 ·10−2 ·10−2
7 7

6 6

5 5
RΘΘ /Θ2c

4 4

3 3

2 2

1 1

0 0
0 0.05 0.1 0.15 0.2 0.25 0.3 0 0.05 0.1 0.15 0.2 0.25 0.3
η η

Figure 7.3.: Variance of the passive scalar fluctuations RΘΘ at Re = 3500 (left) scaled with the centerline passive
scalar Θ2c (7.4) plotted as a function of the similarity coordinate η at different distances from the
orifice: z/D = 25 ( ), 35 ( ), 45 ( ), 55 ( ). For the Re = 7000 case (right), RΘΘ /Θ2c is
plotted at: z/D = 30 ( ), 40 ( ), 50 ( ), 60 ( ).

0.4

0.3
RΘΘ /Θc

0.2
√

0.1

0
0 10 20 30 40 50 60 70
z/D
√
Figure 7.4.: Centerline rms of the passive scalar fluctuation RΘΘ scaled with the centerline mean passive scalar
Θc (7.4) over the distance z from the orifice: Darisse et al. (2015) ( ), Birch et al. (1978) ( ), Babu
and Mahesh (2005) ( ), Lubbers et al. (2001) ( ), present DNS at Re = 3500 ( ), present
DNS at Re = 7000 ( ).

The turbulent heat fluxes shown in Figure 7.5 exhibit a close resemblance to the radial profiles
in Darisse et al. (2015). In Darisse et al. (2015), the measured maximum values for the heat
fluxes are RrΘ /(Uz,c,max Θc ) = 2.2 and RzΘ,max /(Uz,c Θc ) = 3, which is consistent with the
present values being RrΘ /(Uz,c,max Θc ) = 2.2 and RzΘ,max /(Uz,c Θc ) = 3.2 observed in Figure
7.5. In addition, the value at η = 0 for RzΘ /(Uz,c Θc ) is about 0.024 in Darisse et al. (2015),
closely matching 0.025 in the present DNS. In Figure 7.6, the heat fluxes of the simulation at

80
 ·10−2 ·10−2
2.5 3.5

3
2
2.5

RzΘ /(Uz,c Θc )
RrΘ /(Uz,c Θc )

1.5 2

1 1.5

1
0.5
0.5

0 0
0 0.05 0.1 0.15 0.2 0.25 0.3 0 0.05 0.1 0.15 0.2 0.25 0.3
η η

Figure 7.5.: Turbulent heat fluxes RrΘ and RzΘ normalized with Uz,c Θc for the Re = 3500 case at different
distances from the orifice: z/D = 25 ( ), 35 ( ), 45 ( ), 55 ( ).

·10−2 ·10−2
2.5 3.5

3
2
2.5
RzΘ /(Uz,c Θc )
RrΘ /(Uz,c Θc )

1.5 2

1 1.5

1
0.5
0.5

0 0
0 0.05 0.1 0.15 0.2 0.25 0.3 0 0.05 0.1 0.15 0.2 0.25 0.3
η η

Figure 7.6.: Turbulent heat fluxes RrΘ and RzΘ normalized with Uz,c Θc for the Re = 7000 case at different
distances from the orifice: z/D = 30 ( ), 40 ( ), 50 ( ), 60 ( ).

81
Re = 7000 are shown which values are equivalent to the present DNS at Re = 3500 except for
the value at η = 0 for RzΘ /(Uz,c Θc ) being slightly larger at about 0.027.
Considering that Darisse et al. (2015) conducted their experiment at Re = 16,000 with a
top-hat velocity profile, a higher Reynolds number and a turbulent velocity inlet might increase
the axial turbulent heat flux on the centerline.

7.1.3. Third order mixed statistics

There is only quite few experimental data on the radial profiles of the third-order mixed statistics
RiΘΘ , due to their sensitivity to disturbances in experimental setups such as measurement
inaccuracies or boundary effects. In contrast, the present DNS at Re = 3500 provides a well
converged data set, allowing a detailed examination of these statistics which can also be
compared to the DNS data at Re = 7000. Comparable experiments are mainly the helium
jet of Panchapakesan and Lumley (1993b) and the heated jet by Darisse et al. (2015) which
provide a basis for validation. Notably, Antonia et al. (1975) also contributed data for RiΘΘ ,
although using a jet with a strong co-flow. Nevertheless, the shape of the data is similar to the
experiments of Darisse et al. (2015), but with different magnitudes.
The radial profiles of the normalized third-order mixed statistics of the present DNS are
showcased in Figures 7.7 and 7.8 for RiΘΘ at both Reynolds numbers. Qualitatively, the

·10−3 ·10−3
3

2
2
RzΘΘ /(Uz,c Θ2c )
RrΘΘ /(Uz,c Θ2c )

1 1

0 0

−1
−1

−2
0 0.05 0.1 0.15 0.2 0.25 0.3 0 0.05 0.1 0.15 0.2 0.25 0.3
η η

Figure 7.7.: Turbulent heat fluxes RrΘΘ and RzΘΘ at Re = 3500 normalized with Uz,c Θ2c at different distances
from the orifice: z/D = 25 ( ), 35 ( ), 45 ( ), 55 ( ).

shapes of the aforementioned experiments are similar to those of the present DNS. While
the maxima of the present DNS at both Reynolds numbers, being at RrΘΘ,max /(Uz,c Θ2c ) =
2.7 · 10−3 and RzΘΘ,max /(Uz,c Θ2c ) = 2 · 10−3 , closely match RrΘΘ,max /(Uz,c Θ2c ) = 2.8 · 10−3
and RzΘΘ,max /(Uz,c Θ2c ) = 2 · 10−3 (Darisse et al., 2015), the minima have slightly larger
magnitudes at RrΘΘ,min /(Uz,c Θ2c ) = −1.6 · 10−3 and RzΘΘ,min /(Uz,c Θ2c ) = −1 · 10−3 compared
to RrΘΘ,min /(Uz,c Θ2c ) = −1 · 10−3 and RzΘΘ,min /(Uz,c Θ2c ) = −0.3 · 10−3 (Darisse et al., 2015).

82
 ·10−3 ·10−3
3

2
2

RzΘΘ /(Uz,c Θ2c )

RrΘΘ /(Uz,c Θ2c )

1 1

0 0

−1
−1

−2
0 0.05 0.1 0.15 0.2 0.25 0.3 0 0.05 0.1 0.15 0.2 0.25 0.3
η η

Figure 7.8.: Turbulent heat fluxes RrΘΘ and RzΘΘ at Re = 7000 normalized with Uz,c Θ2c at different distances
from the orifice: z/D = 30 ( ), 40 ( ), 50 ( ), 60 ( )

In comparison, Panchapakesan and Lumley (1993b) present an even larger minimum at

RzΘΘ,min /(Uz,c Θ2c ) = −1.5 · 10−3 . RzΘΘ has not yet converged at η = 0, but we receive
RzΘΘ /(Uz,c Θ2c ) = −0.5 · 10−3 at Re = 3500 and RzΘΘ /(Uz,c Θ2c ) = −0.6 · 10−3 at Re = 7000
compared to −0.2 · 10−3 (Darisse et al., 2015) and −0.6 · 10−3 (Panchapakesan and Lumley,
1993b).

The same conclusion as for the heat fluxes in Section 7.1.2 can be drawn here, that the increase
of the centerline value of RzΘΘ might be lead back to higher Reynolds numbers and a turbulent
velocity inlet. However, it is noted that the centerline value is not well converged.

7.1.4. Probability density functions

In Figure 7.9, the PDFs of the passive scalar on the centerline at different orifice distances exhibit
excellent collapse, closely resembling a Gaussian distribution. The calculation method of the
axial velocity PDFs have been described in Section 5.6 and have been performed analogously
for passive scalars. The PDFs are slightly negatively skewed which can also be read off later in
Figure 7.11 at η = 0. The negative skewness on the centerline for a passive scalar was also
observed in the methane jet of Birch et al. (1978). The Gaussian curves of both Reynolds
numbers are very similar due to Gaussian curves only being dependent on the variance, i.e.,
RΘΘ , which we already found to be similar in Figure 7.3.

Figure 7.10 depicts the radial PDF evolution at z/D = 28, 42, 56. As η increases, significant
deviations from Gaussian behavior and the emergence of strongly heavy and skewed tails
become apparent. The PDF approaches a delta distribution for large η due to the passive scalar
being zero in the ambient region. Additionally, we can extract that the PDFs collapse due to
the scaling of η.

83
 101 101

100 100

10−1 10−1
f
10−2 10−2

10−3 10−3

10−4 10−4
0 0.5 1 1.5 2 0 0.5 1 1.5 2
Θ/Θc Θ/Θc

Figure 7.9.: PDFs of Θ(η = 0, z)/Θc (z) for the Re = 3500 case (left) and Re = 7000 case (right) at z/D = 15
( ), 25 ( ), 35 ( ), 45 ( ), 55 ( ), 65 ( ) compared to a Gaussian ( ).

102

f 10−2

0.18 0.2
10−6 0.12 0.14 0.16
2 0.08 0.1
1 0 0 0.02 0.04 0.06
η
Θ(η, z)/Θc (z)

102

f 10−2

0.18 0.2
10−6 0.12 0.14 0.16
2 0.08 0.1
1 0 0 0.02 0.04 0.06
η
Θ(η, z)/Θc (z)

Figure 7.10.: PDFs of Θ(η, z)/Θc (z) at Re = 3500 (above) and Re = 7000 (below) for z/D = 28 ( ), 42 ( ),
56 ( ).

84
 10

6
K, S
4

0 0.02 0.04 0.06 0.08 0.1 0.12 0.14 0.16 0.18 0.2
η

Figure 7.11.: Kurtosis K (above) and skewness S (below) of Θ(η, z)/Θc (z) at Re = 3500 for z/D = 28 ( ), 42
( ), 56 ( ). K = 3, S = 0 (dashed) are the Gaussian values.

10

6
K, S
4

0 0.02 0.04 0.06 0.08 0.1 0.12 0.14 0.16 0.18 0.2
η

Figure 7.12.: Kurtosis K (above) and skewness S (below) of Θ(η, z)/Θc (z) at Re = 7000 for z/D = 28 ( ), 42
( ), 56 ( ). K = 3, S = 0 (dashed) are the Gaussian values.

Quantitative insights from Figure 7.11 and 7.12 show the collapse of skewness S and kurtosis
K for z/D = 28, 42, 56. The Gaussian values, S = 0 and K = 3, serve as references. From
both figures we can also see that the PDF becomes more sub-Gaussian as η increases. This is
consistent with the observations in Birch et al. (1978), where the distribution becomes broader
with increasing η. On the centerline, i.e., at η = 0, the kurtosis for the Re = 7000 case is
K ≈ 3, while for Re = 3500 the distribution is slightly super-Gaussian since K > 3.

For larger η, the mean value of Θ(η, z)/Θc (z) approaches zero. As a result, the PDF can no
longer be approximated by a Gaussian because it would require negative values of Θ(η, z)/Θc (z)
to maintain a Gaussian distribution. In Birch et al. (1978), the distribution becomes bimodal

85
 3

2.5

f 1.5

0.5

0
0 0.2 0.4 0.6 0.8 1 1.2 1.4
Θ/Θc

Figure 7.13.: PDFs of Θ(η, z = 28)/Θc (z = 28) for η = 0.139 (solid) and η = 0.149 (dashed).

and they interpret this as intermittency in the flow. This is similar to what is observed exemplary
in Figure 7.13 for the Re = 3500 case at η = 0.139 for z/D = 28. The DNS data at Re = 7000
also shows this bimodal behavior, although a graphical representation has been omitted here.
As the distance from the centerline increases to η = 0.149, the probability of zero passive scalar
quickly becomes dominant which indicates a laminar ambient region with no passive scalar.

7.2. Symmetry-based scaling laws for a passive scalar

In this section, the turbulent scaling laws for a passive scalar will be derived by employing
Lie symmetry analysis to the MPVSCEs (3.44). The proposed scaling laws are then validated
against the DNS data of the turbulent round jet flow at both Reynolds numbers Re = 3500 and
Re = 7000.

7.2.1. Symmetries of the governing equations

The following set of selected symmetries are relevant for the derivation of the scaling laws of a
turbulent round jet, which have their origin in the NSEs (2.23), (2.24) and the passive scalar
equations (2.44). In the limit of Re, P e → ∞ they are transferred to the H-approach which
consist of a translation symmetry in space
∗
Tx : t∗ = t, x∗i = xi + axi , Θ = Θ,
∗
HΘ {m}
= HΘ{m} , Hi∗{n} Θ{m} = Hi{n} Θ{m} (7.5)

86
and a scaling symmetry in space, time and of the passive scalar described by the parameters
aSx , aSt and aSθ , respectively,
∗
TS : t∗ = eaSt t, x∗i = eaSx xi , Θ = eaSθ Θ, (7.6)
∗
HΘ {m}
=e maSθ
HΘ{m} , Hi∗{n} Θ{m} =e maSθ +n(aSx −aSt )
Hi{n} Θ{m} ,

where n and m refer to the moment order of the velocity and passive scalar, respectively. As
mentioned in Section 6.1.1, it is noteworthy that the NSEs (2.23) and (2.24) exhibit a restricted
set of scaling symmetries compared to the Euler equations, in particular aSt = 2aSx in the
symmetry above which can be shown by substituting the scaling symmetry (7.5) into the NSEs
(2.23). For turbulent flows at higher Re, however, the viscosity dominates only on length scales
comparable to the Kolmogorov length scales (Oberlack, 2000a).
Additionally, a statistical symmetry on the basis of the linear MPVSCEs in Equation (3.44) is
considered
∗
T Ss : t∗ = t, x∗i = xi , Θ = eaSs Θ, (7.7)
∗
HΘ {m}
=e aSs
HΘ{m} , Hi∗{n} Θ{m} =e aSs
Hi{n} Θ{m} ,

which is a measure of intermittency (Wacławczyk et al., 2014) and does not appear in either
the NSEs (2.23), (2.24) or the passive scalar equation (2.44). For the linear system of the
MPSE (3.38), there exists also the generic symmetry of superposition, which does not play a
role in the present subsection, but only comes into play in Section 7.2.4.

7.2.2. Derivation of invariant solutions

The combination of the symmetries (7.5)-(7.7) leads to the invariant surface condition for the
invariant solutions, i.e., turbulent scaling laws of the velocity and passive scalar moments of
a spatially evolving turbulent round jet (see Equation (4.22)). Although the system (3.38)
constitutes an arbitrary multi-point moments hierarchy, in the following, we will only consider
one-point statistics and scaling laws, i.e., x(1) = x(2) = . . . = x(n+m) , and therefore
n Θm .
Hi{n} Θ{m} = U[i] (7.8)

The related invariant surface condition reads

dr dz dU[i]
n Θm
= = n Θm
aSx r aSx z + az [n(aSx − aSt ) + maSθ + aSs ]U[i]
dΘ dΘm
= = ... = . (7.9)
[aSθ + aSs ]Θ [maSθ + aSs ]Θm
In the presented methodology, symmetry breaking is introduced through flow-specific invariants
for the turbulent jet, namely the conservation of thermal energy (Sadeghi et al., 2021) given by
Z ∞
IΘ = Uz Θ r dr (7.10)
0

and the full momentum integral (6.17). The group parameters of the symmetries (7.5)-(7.7)
are constrained through these invariants.

87
The symmetry breaking is induced by implementing the symmetries (7.5)-(7.7) in (7.10)
through
∗
Uz Θ = e−(aSθ +aSx −aSt +aSs ) Uz Θ , r = e−aSx r∗ , (7.11)
which gives Z ∞
∗
IΘ = e −(aSθ +3aSx −aSt +aSs )
Uz Θ r∗ dr∗ (7.12)
0
and in the full momentum invariant IO (6.17) for which we obtain Equation (6.19). Inferring
IΘ (7.10) and IO (6.17) being invariants, i.e., constants, result in the symmetry breaking

aSθ + 3aSx − aSt + aSs = 0, 4aSx − 2aSt + aSs = 0 (7.13)

to maintain invariance. When solved for aSt and aSθ , we find

1 1
aSt = 2aSx + aSs , aSθ = −aSx − aSs . (7.14)
2 2
In Chapter 6, we found that the statistical symmetry (7.7), giving a measure of intermittency
(Wacławczyk et al., 2014), allows for a possible variation of the turbulent decay behavior of the
moments in z-direction agreeing with the hypothesis of George (1989). Although, the DNS
data has shown that the intermittency symmetry is broken, i.e., (7.7) aSs = 0, this theory still
holds. Nevertheless, in Chapter 6, the intermittency reappears in the reduced equations of the
MPMEs (3.30) through the scaling laws. Subsequently, it is shown that this also applies to
passive scalars as detailed later in Section 7.2.4.
Under the consideration of the symmetry breaking (7.14), the turbulent scaling laws can be
derived by integrating the invariant surface condition in (7.9), receiving
r
η= , (7.15)
z − z0
U^n Θm (η)eci,nm (n+m)
Uin Θm (r, z) = i
1 ∗ ∗
, (7.16)
(z − z0 )−(n+m)(1+ 2 aSs )+aSs
Θ
g m (η)ecm m
Θm (r, z) = 1 ∗ ∗
, (7.17)
(z − z0 )−m(1+ 2 aSs )+aSs

where z0 = −az /aSx is the virtual origin and a∗Ss = aSs /aSx .
The invariant η (7.15) is received by integrating the first two terms of (7.9), the scaling law
(7.16) results from the integration of the second and third term and finally the scaling law (7.17)
appears after integrating the second and the remaining terms in (7.9). The integration constants
ci,nm and cm emerge from the integration and move into the exponent after exponentiation.
From (7.16), we see not only the turbulent scaling law of the mixed moments but also a
generalization of the velocity and passive scalar scaling laws. For m = 0, the velocity scaling
laws in (6.25) are received, while for n = 0 we obtain the passive scalar scaling law (7.17).

7.2.3. Validation of the scaling laws on the centerline

Just as in Chapter 6, the scaling parameter of the statistical symmetry a∗Ss generates a one-
parameter family of scaling laws which, as conjectured by George (1989), may be induced

88
by a variation of the inflow condition. However, the purely hydrodynamic turbulence DNS in
Chapter 6 has already shown that a∗Ss = aSs = 0 (6.27) and we find this result confirmed by
the coupling of the NSEs (2.9), (2.36) with the scalar equation (2.44). Furthermore, it can
be inferred that a∗Ss cannot be influenced by a passive scalar but is solely defined through the
velocity inflow condition. Physically, it can be explained by the dominance of convection over
the diffusion since we have P e  1 for the passive scalar. Under this assumption, the passive
scalar “follows” the velocity field due to convection.
The validation of the scaling laws (7.15)-(7.17) is first conducted on the centerline, where
r = 0. To facilitate a comparison with the DNS results, the scaling laws presented in (7.16)
and (7.17) are expressed as

Θ,m e
Θ
gm (η = 0)α cm m
Θm (r = 0, z) = , (7.18)
(z − z0 )m

where αΘ,m is extracted so that Θ

gm (η = 0) = 1 and

i Θ (η = 0)αiΘ,nm e
U^
n m ci,nm (n+m)
Uin Θm (r = 0, z) = , (7.19)
(z − z0 )n+m

where we again extract αiΘ,nm so that the invariant U^i Θ (η = 0) = 1. These formulations
n m

allow for a direct comparison between the scaling laws and the DNS data at the centerline.

Results for the passive scalar scaling law

Analyzing the DNS data shown in Figure 7.14, we conclude that for n = 0 at Re = 3500 we
obtain
cm = c = 1.755, αΘ,m = αΘ = 0.722 (7.20)
for the purely passive scalar scaling (7.17). The parameters for the Re = 7000 case for n = 0
are
c = 1.782, αΘ = 0.706. (7.21)
Qualitatively, similar to Figure 6.1 for velocity scaling laws, the figure for the Re = 7000 case
would look the same as Figure 7.14, therefore it is not shown here. Note that both constants
are independent of the moment order m. This implies that c is an invariant with respect to m.
So far, this cannot be derived directly from symmetry theory, but support is essentially from
the DNS data.

Results for the mixed velocity-scalar scaling law

Extending this to the scaling law of the mixed moments for i = z in (7.19), it can be observed
that the constants in Figure 7.14 span a plane. For i = r, ϕ, a “discrete” plane is formed
because the moments are zero for odd n on the centerline in Figure 7.15 for the Re = 3500.
As explained in Chapter 6, this results from the fact that the pure velocity scaling law (m = 0)
for i = r, ϕ is described by a centered Gaussian process. Since the analysis is performed on the
centerline, the scaling law (7.19) gives the same values for i = r, ϕ due to symmetry. Therefore,

89
 108
αzΘ,nm ecz,nm (n+m)

105

102
10
8
10−1 6
0 4
2
4 2 m
6
8
n 10 0

Figure 7.14.: The exponential prefactor αzΘ,nm ecz,nm (n+m) ( ) from Equation (7.19) determined with the DNS
data at Re = 3500 for mixed moments up to n + m = 6 ( ), for Θ moments up to m = 10 ( ) and
Uz moments up to n = 10 ( ) is shown. Additionally, Equation (7.19) is highlighted for m = 0 ( )
and n = 0 ( ).

only i = r will be considered in the following. The prefactors of the scaling law in (7.19) can
be characterized by an arithmetic weighing of the moments with
nαi,n + mαΘ
αiΘ,nm = (7.22)
n+m
and
nci + mc
ci,nm = , (7.23)
n+m
where
√ n
αz,n = αz = 0.673, cz = 2.123, αr,n = 2e [1.81(n − 1)] 2 , cr = −0.5 (Re = 3500)
(7.24)
and
√ n
αz,n = αz = 0.662, cz = 2.127, αr,n = 2e [1.87(n − 1)] 2 ,
(Re = 7000) cr = −0.5
(7.25)
are taken from Chapter 6. Due to the similarity of the values to the Re = 3500 case, the
graphical representation of the Re = 7000 case has been omitted and only Figure 7.15 is shown,
which represents the Re = 3500 case.

90
 108
αrΘ,nm ecr,nm (n+m)

105

102
10
8
10−10 6
2 4
4 2 m
6
8
n 10 0

Figure 7.15.: The exponential prefactor αrΘ,nm ecr,nm (n+m) ( ) from (7.19) determined with the DNS data at
Re = 3500 for mixed moments up to n + m = 6 ( ), for Θ moments up to m = 10 ( ) and Ur
moments up to n = 10 ( ) is shown. Additionally, Equation (7.19) is highlighted for m = 0 ( )
and n = 0 ( ).

7.2.4. Gaussian behavior of the higher passive scalar moments

Interestingly, the passive scalar moments behave similarly to the axial velocity moments.
The higher instantaneous H-moments of the passive scalar exhibit a Gaussian-like behavior
in η, similar to what is observed for the axial velocity moments. By analyzing the scaled
higher moments of the passive scalar according to Equation (7.17), we find that a Gaussian
function provides a highly accurate representation for these moments. The Gaussian function
is expressed as
Θm 2
m (η) = e−γm η ,
m
=Θ g (7.26)
Θc
where η represents the dimensionless radial coordinate as defined in Equation (7.15). This
Gaussian function captures the characteristic distribution of the passive scalar, showing a rapid
decrease in amplitude as the distance from the central region increases. While it may be
m
tempting to assume Θm ≈ Θ , implying that γm is linear and thus trivial due to the strong
dominance of the the mean passive scalar Θ over the fluctuations, our analysis demonstrates
that γm does not behave linearly for larger moments as is already implied by Figure 7.3 for
RΘΘ and the fluctuations are significant.

In Figures 7.16 and 7.17, we observe the self-preservation of six selected moments of the passive
scalar, ranging up to order ten, at multiple stages along the spatial domain z/D = 25 − 55
for both Reynolds numbers. Except for the central region, there is a collapse of the data
when normalized using the scaling law described in Equation (7.18). This collapse suggests
an universal scaling behavior for the moments of the passive scalar, reinforcing the idea of
self-preservation. Additionally, the Gaussian-like curves become progressively narrower as the

91
moment order increases. The remaining four moments for both Reynolds numbers can be
depicted in Appendix A.2.
The DNS data up to the tenth axial moment is presented in semi-logarithmic scaling in Fig-
ure 7.18 for both Reynolds numbers, demonstrating a parabolic trend, suggesting a strong
agreement with Equation (7.26), although the specific value of γm remains undetermined. The
deviation of the DNS data from the parabolic curve at the edge is likely attributed to external
intermittency as seen in Figure 7.10.
A comparison between Figures 7.16, 7.17 and Figure 7.18 reveals that the turbulent scaling
laws hold true within the fully turbulent regime. It is important to note that turbulent scaling
laws are usually only valid within the confines of a fully turbulent regime and become less
applicable as we move into the ambient environment.
The estimation of γm in m for the passive scalar moments can be obtained from Figure 7.19
through a straightforward nonlinear fitting procedure. The resulting fit yields

γm = −1.27m2 + 37.52m + 27.34 (Re = 3500), γm = −1.39m2 + 39.59m + 24.04 (Re = 7000)
(7.27)
m
with a coefficient of determination of R2 = 0.99, clearly showing that Θm 6≈ Θ and therefore
strong intermittency is inferred. This aligns with what is observed in the PDFs in Figure 7.10,
as well as in the higher standardized moments in Figures 7.11 and 7.12, where intermittency
is evident.

92
 1.2

0.8
Θm /Θm
c

0.6

0.4

0.2

0
0 0.1 0.2 0.3 0 0.1 0.2 0.3
η η
(a) m = 1 (b) m = 2

1.2

0.8
Θm /Θm
c

0.6

0.4

0.2

0
0 0.1 0.2 0.3 0 0.1 0.2 0.3
η η
(c) m = 4 (d) m = 6

1.2

0.8
Θm /Θm
c

0.6

0.4

0.2

0
0 0.1 0.2 0.3 0 0.1 0.2 0.3
η η
(e) m = 8 (f) m = 10

Figure 7.16.: The radial profiles of the mth axial moment normalized with the scalings in (7.17) at different distances
from the orifice for the Re = 3500 case: z/D = 25 ( ), 35 ( ), 45 ( ), 55 ( ). The black
solid lines indicate the Gaussian from Equation (7.26) using γm from Equation (7.27).

93
 1.2

0.8
Θm /Θm
c

0.6

0.4

0.2

0
0 0.1 0.2 0.3 0 0.1 0.2 0.3
η η
(a) m = 1 (b) m = 2

1.2

0.8
Θm /Θm
c

0.6

0.4

0.2

0
0 0.1 0.2 0.3 0 0.1 0.2 0.3
η η
(c) m = 4 (d) m = 6

1.2

0.8
Θm /Θm
c

0.6

0.4

0.2

0
0 0.1 0.2 0.3 0 0.1 0.2 0.3
η η
(e) m = 8 (f) m = 10

Figure 7.17.: The radial profiles of the mth axial moment normalized with the scalings in (7.17) at different distances
from the orifice for the Re = 7000 case: z/D = 30 ( ), 40 ( ), 50 ( ), 60 ( ). The black
solid lines indicate the Gaussian from Equation (7.26) using γm from Equation (7.27).

94
 10−1 10−1

10−4 10−4

Θm 10−7 10−7

10−10 10−10

10−13 10−13
0 0.05 0.1 0.15 0.2 0 0.05 0.1 0.15 0.2
η η

Figure 7.18.: The radial profiles (blue, dashed) of the 1st (top) up to the 10th (bottom) axial moment at Re = 3500
(left) and Re = 7000 (right) and the corresponding Gaussian (black) from Equation (7.26) shown in
a semi-logarithmic plot at z/D = 45.

300
Re = 3500
Re = 3500 fit
250 Re = 7000
Re = 7000 fit
200
γm
150

100

50
1 2 3 4 5 6 7 8 9 10
m
Figure 7.19.: Constants γm from Equation (7.26) are shown for each moment up to order n = 10 determined by
fitting to the DNS yielding the following fits: γm = −1.27m2 + 37.52m + 27.34 (Re = 3500) and
γm = −1.39m2 + 39.59m + 24.04 (Re = 7000).

95
7.2.5. Hidden intermittency in passive scalars

To get a first hint towards the understanding of the emergence of Equation (7.26) as a function
of the radial coordinate η, the MPSE (3.38) is transformed into cylindrical coordinates. Note
that unlike in Section 6.3.1, where we derived the reduced equations for up to n = 2 at one point,
here we are able to transform the MPSE (3.38) for arbitrary order moments. Subsequently, the
scaling laws (7.15)-(7.17) are introduced into the cylindrical-transformed MPSE, receiving
m
"  
∂ H̃Θ{m} X 1 ∂η(l) H̃Θ{m+1} [Θ(m+1) 7→r(l) ] x̃(m+1) 7→ x̃(l)
+
∂ t̃ η(l) ∂η(l)
l=1
 #
∂ H̃Θ{m+1} [Θ(m+1) 7→z(1) ] x̃(m+1) 7→ x̃(1)
−x̃k(l)
∂ x̃k(l)
m
"  
X 1 ∂ H̃Θ{m+1} [Θ(m+1) 7→ϕ(l) ] x̃(m+1) 7→ x̃(l)
+
η(l) ∂ ϕ̃(l)
l=2
 #
∂ H̃Θ{m+1} [Θ(m+1) 7→z(l) ] x̃(m+1) 7→ x̃(l)
+
∂ z̃(l)
   
1 ∗
− (m + 1) 1 + aSs − aSs H̃Θ{m+1} [Θ(m+1) 7→z(1) ] x̃(m+1) 7→ x̃(1) = 0 (7.28)
∗
 
2
for m = 1, . . . , ∞ and x̃k(l) 6= ϕ̃(l) with a reference point x̃(1)
T
r(1)


T
x̃(1) = η(1) , 0, 0 = , 0, 0 (7.29)
z(1) − z0

and the points x̃(l)

T
r(l) z(l)


T
x̃(l) = η(l) , ϕ̃(l) , z̃(l) = , ϕ̃(l) , (7.30)
z(1) − z0 z(1) − z0

for l = 2, . . . , ∞ and where t̃ emerges from the scaling symmetry (7.6) as

t
t̃ = 1 ∗ . (7.31)
(z − z0 )2+ 2 aSs
Detailed steps for obtaining Equations (7.28) and (7.31) can be taken from Appendix A.3.
It is important to note that x̃(l) must contain ϕ̃(l) to ensure their unique definition, e.g.,
when looking at correlations on circular arcs around the jet center. Equation (7.28) has been
reduced by one independent variable through the scaling laws. After a symmetry reduction,
the symmetry properties change in relation to the original equation. In particular, scaling
symmetries can disappear, but various symmetries are also inherited.
Presently, symmetries admitted by Equation (7.28) include

T˜ S : t̃∗ = eãSt t̃, ∗

η(·) = eãS η(·) , ϕ̃∗(·) = ϕ̃(·) , ∗
z̃(·) = z̃(·) ,
∗
H̃Θ {m}
= emãSθ −ãS H̃Θ{m} , H̃r∗{1} Θ{m} = emãSθ −ãSt H̃r{1} Θ{m} ,
H̃ϕ∗ {1} Θ{m} = emãSθ −ãSt H̃ϕ{1} Θ{m} , H̃z∗{1} Θ{m} = emãSθ −ãS −ãSt H̃z{1} Θ{m} (7.32)

96
with H̃Θ{m+1} [Θ(m+1) 7→i(l) ] = H̃i{1} Θ{m} and i(l) = r(l) , ϕ(l) , z(l) . However, caution is necessary
regarding the scaling symmetries given here, because although inserting (7.32) into the
cylindrical MPMEs (7.28) seems to prove the symmetry properties, the symmetries (7.32) are
not necessarily symmetries of the complete system (3.44), being the MPVSCEs, in reduced
form. Therefore, we are dealing with a kind of partial invariance, which is not to be confused
with the mathematical concept of partially invariant solutions (see e.g. Meleshko, 2005).
The statistical symmetry with the scaling parameter ãSs is inherited from the statistical sym-
metry (7.7)

T˜ Ss : x̃∗(·) = x̃(·) , H̃Θ

∗
{m}
= eãSs H̃Θ{m} , H̃i∗{1} Θ{m} = eãSs H̃i{1} Θ{m} (7.33)

and the superposition symmetry of linear systems is presently given by

T̃S : t∗ = t, x̃∗(·) = x̃(·) , ∗

H̃Θ {m}
0
= H̃Θ{m} + H̃Θ {m}
,
H̃i∗{1} Θ{m} = H̃i{1} Θ{m} + H̃i0{1} Θ{m} , (7.34)

where the prime quantities are any additional independent solution of reduced MPSE (7.28).
Here, the superposition symmetry (7.34) plays a crucial role, as will be explained later.
With the symmetries (7.32)-(7.33), a invariant surface condition is, once again, derived

dη dH̃Θ{m}
= = ..., (7.35)
ãS η (mãSθ − ãS + ãSs )H̃Θ{m}

where non-relevant terms have been omitted. Integrating the invariant surface condition (7.35)
leads to
H̃Θ{m} = Cm η q , (7.36)
where Cm is the integration constant and
ãSθ ãSs
q=m + − 1. (7.37)
ãS ãS
For the connection of the Gaussian behavior of the passive scalar moments (7.26) to (7.36), the
superposition symmetry (7.34) reveals its importance, allowing for a superposition of (7.36).
The Taylor series of Equation (7.26)
2 1
2 1
3
e−γm η = 1 − γm η 2 + γm η 2 − γm η 2 + . . . (7.38)
2 6
shows that Equation (7.36) is the basis function of (7.26). Therefore, the sum of the basis
functions (7.36) can be written as
∞
2
X
H̃Θ{m} = Cmk η qk = e−γm η , (7.39)
k=0

where
(−γm )k
Cmk = , qk = 2k, (7.40)
k!
showing the importance of the intermittency symmetry and the superposition symmetry which
are both only admitted by the equations in the H-formulation and not by the Euler or NSEs.

97
7.2.6. Gaussian behavior of the higher mixed moments

We have gathered data for higher mixed moments up to n + m = 6, aiming to demonstrate

their composition as a combination of axial velocity and passive scalar moments, similar to
Subsection 7.2.3.
As illustrated in Figures 7.20 and 7.21, we observe that the mixed moments also exhibit
Gaussian behavior at high precision. This characteristic is attributed to both the higher axial
velocity and passive scalar moments displaying similar trends. In Chapter 6 for the axial
velocity moments, we established, that for Gaussian behavior of the axial velocity moments
(6.36), the Gaussian coefficient γn is approximated by

γn = −1.54n2 + 61.95n + 27.10 (Re = 3500),

γn = −1.51n2 + 63.39n + 28.21 (Re = 7000). (7.41)

From Gaussian behaviors (7.26) and (6.36) we may invoke

−γnm η 2
z Θ (η) = e
U^ (7.42)
n m .

The subsequent approximation is a non-linear n-m surface fit to the DNS data as shown in
Figure 7.20 and in Figure 7.21 is

γnm = −1.51n2 − 1.29m2 − 0.03nm + 61.20n + 37.34m + 30.75 (Re = 3500),

2 2
γnm = −1.49n − 1.39m − 0.47nm + 63.08n + 38.97m + 29.17 (Re = 7000), (7.43)

where the coefficient of determination is R2 = 0.99 for both Reynolds numbers and a repre-
sentation of the fit is shown in Figure 7.22 for the Re = 3500 case. Here, the figure of the
Re = 7000 case is omitted due to the strong similarity to Figure 7.22.
The Gaussian coefficient γn and γm in (7.41) and (7.27), respectively, may be compared with
γ in (7.43) by tentatively assuming U
nm
fn Θ
g
z
m ≈U ^n Θm which implies γ + γ ≈ γ
z n m . Com-
nm
paring the Gaussian coefficients of (7.27), (7.41) and (7.43) reveals an approximate halving of
the constant terms and the emergence of a cross-product term, while the remaining terms only
differ by 2%. This observation suggests that the moments are coupled and demonstrate some
degree of cross-correlation. For the Re = 7000 case, the cross-correlation is quite significant at
approximately 30% of the quadratic terms. Nevertheless, the approximation (7.43) reduces
approximately to the special forms γm (7.27) and γn (7.41) in the limit cases n = 0 and m = 0.

98
 10−1

10−3
Uzn Θm

10−5

10−7

10−9

10−11
0 0.05 0.1 0.15 0.2 0 0.05 0.1 0.15 0.2
η η
(a) n = 1 (b) n = 2

10−1

10−3
Uzn Θm

10−5

10−7

10−9

10−11
0 0.05 0.1 0.15 0.2 0 0.05 0.1 0.15 0.2
η η
(c) n = 3 (d) n = 4

10−1

10−3

10−5
Uzn Θm
10−7

10−9

10−11
0 0.05 0.1 0.15 0.2
η
(e) n = 5

Figure 7.20.: The radial profiles (blue, dashed) of m = 1 (top) up to the n + m = 6 (bottom) axial mixed moment
Uzn Θm at Re = 3500 and the corresponding Gaussian (black) from Equation (7.26) shown in a
semi-logarithmic plot at z/D = 45.

99
 10−1

10−3
Uzn Θm

10−5

10−7

10−9

10−11
0 0.05 0.1 0.15 0.2 0 0.05 0.1 0.15 0.2
η η
(a) n = 1 (b) n = 2

10−1

10−3
Uzn Θm

10−5

10−7

10−9

10−11
0 0.05 0.1 0.15 0.2 0 0.05 0.1 0.15 0.2
η η
(c) n = 3 (d) n = 4

10−1

10−3

10−5
Uzn Θm
10−7

10−9

10−11
0 0.05 0.1 0.15 0.2
η
(e) n = 5

Figure 7.21.: The radial profiles (blue, dashed) of m = 1 (top) up to the n + m = 6 (bottom) axial mixed moment
Uzn Θm at Re = 7000 and the corresponding Gaussian (black) from Equation (7.26) shown in a
semi-logarithmic plot at z/D = 45.

100
 γnm
γnm fit

400
γnm

200
10
8
6
0 4
2
4 2 m
6
8
n 10 0
Figure 7.22.: Constants γnm from Equation (7.42) at Re = 3500 are shown for pure moments up to order n, m = 10
and mixed moments up to order n + m = 6 determined by fitting to the DNS yielding the following
fit: γnm = −1.51n2 − 1.29m2 − 0.03nm + 61.2n + 37.34m + 30.75.

7.3. Conclusive remarks

This chapter has two integral components: a comprehensive DNS investigation of a turbulent
round jet flow including a scalar at two Reynolds numbers and the derivation of group invariant
solutions using Lie symmetry analysis applied to the MPVSCEs, validated by the DNS.
The DNS effort yielded new and detailed data on passive scalar moments up to the tenth order
for a turbulent round jet at Re = 3500 and P r = 0.71, where the Prandtl number corresponds
to the value of air. Additionally, previously unreported mixed velocity-passive scalar moments
up to the sixth order were extracted. The highly converged data, averaged over 75,000D/Ub
time units exhibited a remarkable collapse in radial profiles for various statistics in z/D = 10
intervals between z/D = 15 − 55. Further, a DNS at Re = 7000 and P r = 0.71 has been
conducted and the statistical data are compared to the data of the DNS at the lower Reynolds
number.
The turbulent pipe velocity profile has shifted the virtual origins of the velocity and passive
scalar moments closer to each other. The large numerical box reveals self-similarity in the
far-field for the rms of the passive scalar fluctuation which has not been visible by previous
numerical studies due to their small box size. In the near-field, a higher Reynolds number
smooths out the initial maximum visible in the rms of the passive scalar fluctuation for the
lower Reynolds number case. The self-similar profiles of the second and third order correlation
between the passive scalar and the axial velocity show an increase on the centerline which is
attributed to the turbulent velocity inlet and a higher Reynolds number. The passive scalar
PDFs on the centerline are slightly negatively skewed compared to the axial velocity PDFs.
For larger distances from the centerline, the PDFs become bimodal, indicating intermittency,
before approaching a delta distribution in the ambient region.

101
Lie symmetry analysis has been applied to a system of differential equations consisting of the
MPSE and the MPVSCEs, a combination of the MPSE and MPMEs. Both have been derived in
their H-form (instantaneous form) from which, using their symmetries, a set of self-similar
solutions, also called turbulent scaling laws, have been derived. The resulting scaling laws are
then validated against the aforementioned DNS, demonstrating a remarkable agreement for
passive scalar moments up to the tenth order and for mixed moments up to the sixth order,
linking the scaling laws for the velocity and passive scalar moments. In particular, the group
parameter of the statistical scaling symmetry mentioned above, is uniquely determined by
the velocity inflow condition. This means, self-similarity of a passive scalar is universal for a
specific velocity inflow condition.
An intriguing observation emerged from the study of the pure instantaneous passive scalar
moments, revealing a Gaussian-like behavior found in radial direction only in statistical sym-
metries, similar to the pure axial velocity moments. The higher moments were accurately
represented by Gaussian functions linked by a non-linear coefficient γm , challenging the as-
sumption of linear behavior due to the dominance of the mean. This non-linear behavior,
exemplified by the parameter γm , signifies intermittency, which is supported by the findings in
the PDF being heavy-tailed for larger η. Furthermore, the Gaussian behavior extends to mixed
axial velocity-scalar moments, where the coefficient γnm indicates that the axial velocity and
the passive scalar exhibit a clear correlation.

102
8. Conclusion

As we near the end of this dissertation, this final chapter is used to reflect on the insights gained.
In this chapter, we summarize our findings and articulate the significance of our discoveries.
In addition, we envision the potential trajectories and emerging research opportunities that lie
ahead in symmetry-based turbulence theory and turbulent round jet flows.

8.1. Closing remarks

In this dissertation, two large-scale direct numerical simulations (DNS) of turbulent round jet
flows at the Reynolds numbers Re = 3500 and Re = 7000 with a passive scalar at a Prandtl
number of P r = 0.71, describing air, have been performed. The highly converged DNS data
is thoroughly discussed and compared to previous studies. Furthermore, Lie symmetries and
exact statistical descriptions of turbulence are introduced, leading to the derivation of invariant
solutions, known as turbulent scaling laws, for velocity and scalar moments of arbitrary order.
Subsequently, the DNS data is used to validate these turbulent scaling laws.
The size of the numerical box at these Reynolds numbers is unprecedented, extending up
to 75D in the axial direction and 65D in the radial direction to ensure accurate collection
of far-field statistics. The statistical data is averaged over 75,000D/Ub units for the DNS at
Re = 3500, equivalent to 200 passes of a particle through the domain, and over 30,000D/Ub
units for the DNS at Re = 7000, corresponding to 80 passes of a particle. Previous studies have
not achieved a DNS at this box size with statistics of comparable quality.
The high-quality statistics generated include first and second order one-point velocity moment
statistics, Reynolds stress and turbulent kinetic energy budgets, generated directly from the
DNS data. Additionally, novel contributions include the generation of axial velocity probability
density functions (PDFs) at different r/D and z/D. Notably, the all statistics exhibit a highly
converged 1/z n behavior, where n is the moment order, throughout the domain, while the
radial velocity profiles show a remarkable collapse based on classical scaling. Furthermore, the
axial velocity PDFs also exhibits exceptional self-similarity. The observed slight axial increase in
turbulent intensities aligns with experimental findings, previously inaccessible in DNS studies
due to box size limitations.
The use of a fully turbulent velocity profile as an inlet reduces the size of the potential core
causing an earlier onset of decay. Moreover, this velocity profile may have effects on dissipation
similar to an increasing Reynolds number.
Comparison of the DNS results at the two Reynolds numbers reveals a shift in the virtual origin,
with the magnitude of the Reynolds normal stresses increasing at higher Reynolds numbers

103
while the Reynolds shear stress remains constant. Analysis of the turbulent intensities shows
that the self-similarity shifts to greater distances from the orifice with increasing Reynolds
number, which can also be observed across other statistical quantities.

Turbulent scaling laws have been derived through Lie symmetry analysis for a turbulent round
jet flow, using statistical symmetries originating from the linear properties of multi-point
moment equations (MPMEs). These symmetries, which are not found in the Navier-Stokes
equations (NSEs), form the basis of our investigation. The introduction of the H-approach,
rooted in instantaneous moments within the MPMEs, stands in contrast to the fluctuation
moments of the R-approach, such as in the Reynolds-averaged Navier-Stokes (RANS) equations.
Further, the group parameter of the statistical symmetry is key to a possible variation of the
turbulent decay behavior in turbulent round jets. The variability of the turbulent decay was
initially suggested by George (1989), who questioned the classical scaling of 1/z, although he
used a simplified form of the momentum conservation invariant whereas here the full form has
been utilized.

Beyond the statistical scaling symmetry, we introduce, for the first time in turbulence, the
repeated use of the superposition principle, which is a symmetry of all linear differential
equations (Bluman et al., 2010) and serves as an elementary basis here. This principle allows
the derivation of fundamental scaling laws for various turbulent canonical flows, representing
solutions of the infinite moment hierarchy.

Furthermore, validating the scaling laws through DNS at the two distinct Reynolds numbers
has revealed minimal changes in the scaling law parameters, with the exception of the virtual
origin. This discrepancy can be attributed to the inherent variance within the flow.

The DNS effort also yields new and detailed insights into the passive scalar moments up to
the tenth order for a turbulent round jet at Re = 3500 and P r = 0.71. Additionally, previously
unreported mixed velocity-passive scalar moments up to the sixth order are extracted. The
highly converged data exhibit an almost perfect collapse in radial profiles for various statistics
in z/D = 10 intervals between z/D = 15 − 55.

The DNS data of the two simulations reveal, that the turbulent pipe velocity profile moves
the virtual origins of the velocity and passive scalar moments closer together. Due to the
present numerical box, self-similarity is indicated in the far-field for the root-mean square
(rms) of the passive scalar fluctuation, which has not been observable in prior numerical
investigations due to the small box size. In the near-field, a larger Reynolds number smooths
out the initial maximum visible in the rms of the passive scalar fluctuation for the lower Reynolds
number. The self-similar profiles of the second and third order correlation between the passive
scalar and the axial velocity is larger on the centerline compared to earlier studies, which is
ascribed to the turbulent velocity inlet and different Reynolds numbers. The passive scalar
PDFs on the centerline are slightly negative skewed compared to the axial velocity PDFs. As
the distance from the centerline increases, the PDFs become bimodal, indicating intermittency,
and eventually approach a delta distribution in the ambient region.

The Lie symmetry analysis has also been applied to a system of differential equations comprising
the multi-point scalar equation (MPSE) and the multi-point velocity-scalar correlation equations
(MPVSCEs), a combination of the MPSE and MPMEs. The comparison of turbulent scaling
laws with the DNS data shows agreement to passive scalar moments up to the tenth order and
for mixed moments up to the sixth order, linking the scaling laws for the velocity and passive

104
scalar moments. According to Lie symmetry theory, the self-similar profiles of a passive scalar
are universal for a given velocity inflow profile.
Lastly, the examination of the pure instantaneous passive scalar moments has revealed a
Gaussian-like behavior in the radial direction found only in statistical symmetries, similar
to the pure axial velocity moments. The higher moments were accurately represented by
Gaussian functions linked by a non-linear coefficient γm , challenging the assumption of linear
behavior due to the dominance of the mean. This non-linear behavior, described by γm , implies
symmetry and is supported by the heavy-tails of the PDF for a larger distance from the centerline.
Furthermore, the Gaussian behavior extends to mixed axial velocity-scalar moments, where
the coefficient γnm indicates that the axial velocity and the scalar have a clear correlation.

8.2. Outlook

The present dissertation opens many avenues for further investigation. The DNS data serves as a
basis for a further exploration of turbulent round jet flow, such as the investigation of correlation
functions or structure functions. In addition, the high quality DNS data at Re = 3500 can be
used to optimize existing turbulence models, however, additional computational resources are
required to achieve well-converged statistics for the DNS at Re = 7000.
Lie symmetry analysis has suggested possible variations in turbulence decay, although these
hypotheses remain unconfirmed by the DNS data. Experimental or numerical validation is
needed, potentially achieved by varying the inlet conditions using synthetic jets, fractal grids,
or chevrons. Interestingly, the parameter describing this variation does not manifest itself in
second order moments, which could be verified by these experiments or simulations. Manipu-
lation of the turbulence decay is an important insight, which implies a possible manipulation
of the spreading rate of a jet. This, in turn, is relevant to many engineering applications, such
as turbulent mixing to achieve high combustion efficiencies and reduce pollutant emissions in
combustion chambers.
Furthermore, while the scaling parameters show minimal changes between simulations, addi-
tional studies at higher Reynolds numbers are essential to verify these observations. In addition,
exploring the dependence of these parameters on the Reynolds number through further Lie
symmetry analysis could prove this analytically although the scaling symmetry of space and
time would merge into one.
The Gaussian distribution of the axial velocity moment profiles has been derived from the
reduced moment equations only up to the second order moment. For the extension of the
derivation to arbitrary orders, the MPMEs have to be transformed to cylindrical coordinates
and reduced, which is quite complicated due to the pressure correlation term. While symmetry
methods enable the derivation of more complex scaling laws, the determination of the scaling
parameters remains an ongoing challenge, in particular for the Gaussian exponent as a function
of the moment order.
The extension of symmetry theory to the PDF approach based on the Lundgren-Novikov-Monin
(LMN) hierarchy can bridge the gap between classical moment scaling, which has its roots in
engineering-based turbulence research, and PDF scaling, which is based on more physically
oriented approaches. All symmetries derived for the MPMEs have their counterparts in the PDF

105
formulation and these allow for the construction of invariant solutions for the PDF. Conditions,
such as the non-negativity of the PDF and boundary conditions in the sample velocity space,
impose constraints on the various group parameters of the symmetries. These additional
conditions on the scaling law parameters contribute to a deeper understanding of the roots of
these parameters.
Finally, symmetry theory, grounded in first principles, is proving to be a powerful tool in turbu-
lence research. Previous findings using various different approaches have been consolidated
within symmetry theory, such as the logarithmic wall law (Oberlack et al., 2022) and both the
classical scaling as well as the scaling proposed by George (1989) for turbulent round jet flows,
as presently discussed. Symmetry theory has greatly aligned with the results obtained from the
DNS of spatially evolving turbulent round jet flows conducted in this dissertation. Consequently,
the combination of symmetry theory and DNS sets an example for the exploration of additional
flows, holding great promise for advancing our understanding of turbulent phenomena.

106
Bibliography

Acheson, D. J. (1990). Elementary fluid dynamics. Oxford applied mathematics and computing
science series. Oxford and New York: Clarendon Press and Oxford University Press.
Alcántara-Ávila, F., García-Raffi, L. M., Hoyas, S., and Oberlack, M. (2024). “Validation of
symmetry-induced high moment velocity and temperature scaling laws in a turbulent channel
flow: (accepted)”. In: Phys. Rev. E.
Antonia, R. A. and Mi, J. (1993). “Temperature dissipation in a turbulent round jet”. In: J.
Fluid Mech. 250, pp. 531–551.
Antonia, R. A., Prabhu, A., and Stephenson, S. E. (1975). “Conditionally sampled measurements
in a heated turbulent jet”. In: J. Fluid Mech. 72.03, p. 455.
Arfken, G. B. and Weber, H.-J. (2005). Mathematical methods for physicists. 6th ed. Boston:
Elsevier.
Avsarkisov, V., Hoyas, S., Oberlack, M., and García-Galache, J. P. (2014). “Turbulent plane
Couette flow at moderately high Reynolds number”. In: J. Fluid Mech. 751.
Babu, P. and Mahesh, K. (2005). “Direct Numerical Simulation of Passive Scalar Mixing in
Spatially Evolving Turbulent Round Jets”. In: 43rd AIAA Aerospace Sciences Meeting and
Exhibit. AIAA.
Babu, P. C. and Mahesh, K. (2004). “Upstream entrainment in numerical simulations of spatially
evolving round jets”. In: Phys. Fluids 16.10, pp. 3699–3705.
Batchelor, G. K. (2000). An Introduction to Fluid Dynamics. Cambridge mathematical library.
Cambridge: Cambridge University Press.
Birch, A. D., Brown, D. R., Dodson, M. G., and Thomas, J. R. (1978). “The turbulent concentra-
tion field of a methane jet”. In: J. Fluid Mech. 88.3, pp. 431–449.
Bird, R. B., Stewart, W. E., and Lightfoot, E. N. (2007). Transport phenomena. 2. ed. rev. New
York, N.Y.: Wiley.
Bluman, G. W., Cheviakov, A. F., and Anco, S. C. (2010). Application of Symmetry Methods to
Partial Differential Equations. Springer.
Bluman, G. W. (2002). Symmetry methods for differential equations. 2nd ed. Vol. v. 154. Applied
mathematical sciences. New York: Springer.
Boersma, B. J., Brethouwer, G., and Nieuwstadt, F. T. M. (1998). “A numerical investigation on
the effect of the inflow conditions on the self-similar region of a round jet”. In: Phys. Fluids
10.4, pp. 899–909.

107
Bogey, C. and Bailly, C. (2006). “Large eddy simulations of transitional round jets: Influence of
the Reynolds number on flow development and energy dissipation”. In: Phys. Fluids 18.6,
p. 065101.
Bogey, C. and Bailly, C. (2009). “Turbulence and energy budget in a self-preserving round jet:
direct evaluation using large eddy simulation”. In: J. Fluid Mech. 627, pp. 129–160.
Bradbury, L. J. S. (1965). “The structure of a self-preserving turbulent plane jet”. In: Journal of
Fluid Mechanics 23.1, pp. 31–64.
Buresti, G. (2015). “A note on Stokes’ hypothesis”. In: Acta Mechanica 226.10, pp. 3555–3559.
Chou, P. Y. (1945). “On velocity correlations and the solutions of the equations of turbulent
fluctuation”. In: Quarterly of Applied Mathematics 3, pp. 38–54.
D. J. Gross (1996). “The role of symmetry in fundamental physics”. In: Proc. Nat. Acad. Sci. 93,
pp. 14256–14259.
Darisse, A., Lemay, J., and Benaïssa, A. (2015). “Budgets of turbulent kinetic energy, Reynolds
stresses, variance of temperature fluctuations and turbulent heat fluxes in a round jet”. In: J.
Fluid Mech. 774, pp. 95–142.
Dong, S., Karniadakis, G. E., and Chryssostomidis, C. (2014). “A robust and accurate outflow
boundary condition for incompressible flow simulations on severely-truncated unbounded
domains”. In: J. Comput. Phys. 261, pp. 83–105.
Dowling, D. R. and Dimotakis, P. E. (1990). “Similarity of the concentration field of gas-phase
turbulent jets”. In: J. Fluid Mech. 218.-1, p. 109.
Einstein, A. (1916). “Die Grundlage der allgemeinen Relativitätstheorie”. In: Annalen der Physik
354.7, pp. 769–822.
El Khoury, G. K., Schlatter, P., Noorani, A., Fischer, P. F., Brethouwer, G., and Johansson, A. V.
(2013). “Direct Numerical Simulation of Turbulent Pipe Flow at Moderately High Reynolds
Numbers”. In: Flow, Turbulence and Combustion 91.3, pp. 475–495.
Ferdman, E., Otugen, M. V., and Kim, S. (2000). “Effect of Initial Velocity Profile on the
Development of Round Jets”. In: Journal of Propulsion and Power 16.4, pp. 676–686.
Fermigier M. (2017). Transport equations: Mass and heat balances.
Fick, A. (1855). “Ueber Diffusion”. In: Annalen der Physik 170.1, pp. 59–86.
Fischer, P. F., Lottes, J. W., and Kerkemeier, S. G. (2008). nek5000 Web page, https : / /
[Link]/.
Frisch, U. and Kolmogorov, A. N. (1995). Turbulence: The legacy of A.N. Kolmogorov / Uriel
Frisch. Cambridge: Cambridge University Press.
George, W. K. (1989). The self-preservation of turbulent flows and its relation to initial conditions
and coherent structures. Advances in Turbulence.
George, W. K. (2009). “Is There an Asymptotic Effect of Initial and Upstream Conditions on
Turbulence?” In: Proceedings of the ASME Fluids Engineering Division Summer Conference,
2008. New York N.Y.: ASME, pp. 647–672.

108
George, W. K. (2012). “Asymptotic Effect of Initial and Upstream Conditions on Turbulence”.
In: Journal of Fluids Engineering 134.6, pp. 061203–061203-27.
George, W. K. and Castillo, L. (1997). “Zero-Pressure-Gradient Turbulent Boundary Layer”. In:
Applied Mechanics Reviews 50.12, pp. 689–729.
Gilliland, T., Ranga-Dinesh, K. K. J., Fairweather, M., Falle, S. A. E. G., Jenkins, K. W., and
Savill, A. M. (2012). “External Intermittency Simulation in Turbulent Round Jets”. In: Flow,
Turbulence and Combustion 89.3, pp. 385–406.
Gutmark, E. and Wygnanski, I. (1976). “The planar turbulent jet”. In: J. Fluid Mech. 73.3,
pp. 465–495.
Hawkins, T. (2000). Emergence of the theory of Lie groups: An essay in the history of mathematics,
1869-1926. Sources and Studies in the History of Mathematics and Physical Sciences. New
York: Springer.
Heskestad, G. (1965). “Hot-Wire Measurements in a Plane Turbulent Jet”. In: Journal of Applied
Mechanics 32.4, pp. 721–734.
Hoyas, S., Oberlack, M., Alcántara-Ávila, F., Kraheberger, S. V., and Laux, J. (2022). “Wall
turbulence at high friction Reynolds numbers”. In: Phys. Rev. Fluids 7.1, p. 014602.
Hoyas, S. and Jiménez, J. (2008). “Reynolds number effects on the Reynolds-stress budgets in
turbulent channels”. In: Phys. Fluids 20.10, p. 101511.
Hunt, J. C. R., Wray, A. A., and Moin, P. (1988). “Eddies, Streams, and Convergence Zones in
Turbulent Flows”. In: Proceedings of the Summer Program 1988.
Hussein, H. J., Capp, S. P., and George, W. K. (1994). “Velocity measurements in a high-
Reynolds-number, momentum-conserving, axisymmetric, turbulent jet”. In: J. Fluid Mech.
258, pp. 31–75.
Ibragimov, N. H. (1994). CRC Handbook of Lie Group Analysis of Differential Equations, Volume
1: Symmetries, Exact Solutions, and Conservation Laws. CRC Press.
Ibragimov, N. H. (1995). CRC Handbook of Lie Group Analysis of Differential Equations, Volume
2: Applications in engineering and physical sciences. CRC Press.
Ibragimov, N. H. (1996). CRC Handbook of Lie Group Analysis of Differential Equations, Volume
3: New trends in theoretical developments and computational methods. CRC Press.
Jambunathan, K., Lai, E., Moss, M. A., and Button, B. L. (1992). “A review of heat transfer
data for single circular jet impingement”. In: International Journal of Heat and Fluid Flow
13.2, pp. 106–115.
K. Brading and E. Castellani (2003). Symmetries in Physics: Philosophical Reflections. Cambridge
University Press.
Karman, T. von and Howarth, L. (1938). “On the Statistical Theory of Isotropic Turbulence”.
In: Proceedings of the Royal Society of London. Series A - Mathematical and Physical Sciences
164.917, pp. 192–215.

109
Klingenberg, D. (2022). “Development of novel Reynolds-averaged Navier-Stokes turbulence
models based on Lie symmetry constraints”. PhD thesis. Darmstadt: Technical University of
Darmstadt.
Kolmogorov, A. N., Levin, V., Hunt, J. C. R., Phillips, O. M., and Williams, D. (1991). “The local
structure of turbulence in incompressible viscous fluid for very large Reynolds numbers”.
In: Proceedings of the Royal Society of London. Series A: Mathematical and Physical Sciences
434.1890, pp. 9–13.
Kraichnan, R. H. (1965). “Lagrangian-history closure approximaiton for turbulence”. In: Phys.
Fluids 8.4, pp. 575–598.
List, E. J. (1982). “Turbulent Jets and Plumes”. In: Annual Review of Fluid Mechanics 14.1,
pp. 189–212.
Liu, X., Xie, Z., and Dong, S. (2020). “On a simple and effective thermal open boundary
condition for convective heat transfer problems”. In: International Journal of Heat and Mass
Transfer 151, p. 119355.
Lubbers, C. L., Brethouwer, G., and Boersma, B. J. (2001). “Simulation of the mixing of a
passive scalar in a round turbulent jet”. In: Fluid Dynamics Research 28.3, pp. 189–208.
Lundgren, T. S. (1967). “Distribution functions in the statistical theory of turbulence”. In: Phys.
Fluids 10, pp. 969–975.
Mansour, N. N., Kim, J., and Moin, P. (1988). “Reynolds-stress and dissipation-rate budgets in
a turbulent channel flow”. In: J. Fluid Mech. 194.-1, p. 15.
Meleshko, S. V. (2005). Methods for constructing exact solutions of partial differential equations:
Mathematical and analytical techniques with applications to engineering. 2005th. Mathematical
and analytical techniques with applications to engineering. New York: Springer Science &
Business Media.
Monin, A. S. (1967). “Equations of turbulent motion”. In: J. Appl. Math. Mech. 31(6), pp. 1057–
1068.
Navier, M. (1827). “Mémoire sur les lois du mouvement des fluides”. In: Mém. de l’Acad. 6,
pp. 389–416.
Nelkin, M. (1994). “Universality and scaling in fully developed turbulence”. In: Advances in
Physics 43.2, pp. 143–181.
Nguyen, C. T. and Oberlack, M. (2024a). “Analysis of a turbulent round jet based on direct
numerical simulation data at large box and high Reynolds number”. In: Physical Review Fluids
9.7, p. 074608.
Nguyen, C. T. and Oberlack, M. (2024b). “Comparative Study of Turbulent Round Jet Flows
through Direct Numerical Simulation at Medium-High Reynolds Numbers”. In: Under review
with Phys. Fluids.
Nguyen, C. T. and Oberlack, M. (2024c). “Hidden intermittency in turbulent jet flows”. In:
Under review with Phys. Rev. Res.
Nguyen, C. T. and Oberlack, M. (2024d). “Passive scalar statistics in a turbulent round jet:
symmetry theory and direct numerical simulation”. In: Under review with J. Fluid Mech.

110
Novikov, E. A. (1968). “Kinetic equations for a vortex field”. In: Sov. Phys.-Dokl. 12, pp. 1006–
1008.
Oberlack, M. (1999). “Similarity in non-rotating and rotating turbulent pipe flows”. In: J. Fluid
Mech. 379, pp. 1–22.
Oberlack, M. (2000a). “Symmetrie, Invarianz und Selbstähnlichkeit in der Turbulenz”. In:
Habilitation thesis.
Oberlack, M. (2001). “A unified approach for symmetries in plane parallel turbulent shear
flows”. In: J. Fluid Mech. 427, pp. 299–328.
Oberlack, M. and Guenther, S. (2003). “Shear-Free Turbulent Diffusion – Classical and New
Scaling Laws, Fluid Dynamic Research”. In: Fluid Dyn. Res. 33, pp. 453–476.
Oberlack, M., Hoyas, S., Kraheberger, S. V., Alcántara-Ávila, F., and Laux, J. (2022). “Turbulence
Statistics of Arbitrary Moments of Wall-Bounded Shear Flows: A Symmetry Approach”. In:
Phys. Rev. Lett. 128.2, p. 024502.
Oberlack, M. and Rosteck, A. (2010). “New statistical symmetries of the multi-point equations
and its importance for turbulent scaling laws”. In: Discrete & Continuous Dynamical Systems -
S 3.3, pp. 451–471.
Oberlack, M. and Waclawczyk, M. (2006). “On the extension of Lie group analysis to functional
differential equations”. In: Arch. Mech. 58, pp. 597–618.
Oberlack, M. (2000b). “Asymptotic Expansion, Symmetry Groups, and Invariant Solutions of
Laminar and Turbulent Wall-Bounded Flows”. In: ZAMM J. Appl. Math. Mech. 80, pp. 791–
800.
Olver, P. J. (2000). Applications of Lie groups to differential equations. 2nd ed. [1st softcover
print.] Vol. 107, Ed.2000. Graduate Texts in Mathematics. New York: Springer.
Panchapakesan, N. R. and Lumley, J. L. (1993a). “Turbulence measurements in axisymmetric
jets of air and helium. Part 1. Air jet”. In: J. Fluid Mech. 246, pp. 197–223.
Panchapakesan, N. R. and Lumley, J. L. (1993b). “Turbulence measurements in axisymmetric
jets of air and helium. Part 2. Helium jet”. In: J. Fluid Mech. 246, pp. 225–247.
Papoulis, A. and Pillai, S. U. (2002). Probability, Random Variables, and Stochastic Processes.
Fourth. Boston: McGraw Hill.
Patera, A. T. (1984). “A spectral element method for fluid dynamics: Laminar flow in a channel
expansion”. In: J. Comput. Phys. 54.3, pp. 468–488.
Pope, S. B. (2000). Turbulent flows. Cambridge: Cambridge University Press.
Rapp, B. E. (2017). “Chapter 9 - Fluids”. In: Microfluidics: Modelling, Mechanics and Mathematics.
Ed. by Bastian E. Rapp. Micro and Nano Technologies. Oxford: Elsevier, pp. 243–263.
Reynolds, O. (1883). “XXIX. An experimental investigation of the circumstances which deter-
mine whether the motion of water shall be direct or sinuous, and of the law of resistance in
parallel channels”. In: Philosophical Transactions of the Royal Society of London 174, pp. 935–
982.

111
Reynolds, O. (1895). “On the dynamical theory of incompressible viscous fluids and the
determination of the criterion”. In: Philosophical Transactions of the Royal Society of London.
(A.) 186, pp. 123–164.
Reynolds, O., Brightmore, A. W., and Moorby, W. H. (1903). Papers on Mechanical and Physical
Subjects: The sub-mechanics of the universe. Vol. 3. Papers on Mechanical and Physical Subjects.
The University Press.
Rezaeiravesh, S., Vinuesa, R., and Schlatter, P. (2019). A statistics toolbox for turbulent pipe
flow in Nek5000.
Rosteck, A. (2013). “Scaling laws in turbulence - A theoretical approach using Lie-point
symmetries”. In: Doctoral thesis.
Sad Chemloul, N.-E. (2020). Dimensional analysis and similarity in fluid mechanics. 1st edition.
Hoboken: ISTE Ltd / John Wiley and Sons Inc.
Sadeghi, H., Oberlack, M., and Gauding, M. (2018). “On new scaling laws in a temporally
evolving turbulent plane jet using Lie symmetry analysis and direct numerical simulation”.
In: J. Fluid Mech. 854, pp. 233–260.
Sadeghi, H., Oberlack, M., and Gauding, M. (2021). “New symmetry-induced scaling laws of
passive scalar transport in turbulent plane jets”. In: J. Fluid Mech. 919.
Sadeghi, H., Lavoie, P., and Pollard, A. (2014). “The effect of Reynolds number on the scaling
range along the centreline of a round turbulent jet”. In: Journal of Turbulence Volume 15,
pp. 335–349.
Sadeghi, H. and Pollard, A. (2012). “Effects of passive control rings positioned in the shear
layer and potential core of a turbulent round jet”. In: Phys. Fluids 24.
Schlichting, H. (1933). “Laminare Strahlausbreitung”. In: ZAMM - Zeitschrift für Angewandte
Mathematik und Mechanik 13.4, pp. 260–263.
Schlichting, H. and Gersten, K. (2017). Boundary-Layer Theory. 9th ed. 2017. Berlin Heidelberg:
Springer Berlin Heidelberg and Imprint: Springer.
Sharan, N. and Bellan, J. (2021). “Investigation of high-pressure turbulent jets using direct
numerical simulation”. In: J. Fluid Mech. 922.
Shin, D., Sandberg, R. D., and Richardson, E. S. (2017). “Self-similarity of fluid residence time
statistics in a turbulent round jet”. In: J. Fluid Mech. 823, pp. 1–25.
Simon, V., Gomaa, H., and Weigand, B. (2017). Dimensional Analysis for Engineers. 1st ed. 2017.
Mathematical Engineering. Cham: Springer International Publishing and Imprint: Springer.
Spurk, J. and Aksel, N. (2010). Strömungslehre: Einführung in die Theorie der Strömungen. 8.
Aufl. 2010. Springer-Lehrbuch. Berlin, Heidelberg: Springer Berlin Heidelberg.
Sreenivasan, K. R. and Antonia, R. A. (1997). “THE PHENOMENOLOGY OF SMALL-SCALE
TURBULENCE”. In: Annual Review of Fluid Mechanics 29.1, pp. 435–472.
Stokes, G. G. (1849). “On the theories of the internal friction of fluids in motion, and of the
equilibrium and motion of elastic solids”. In: Trans. Cambr. Phil. Soc. 8, pp. 287–319.

112
Taub, G. N., Lee, H., Balachandar, S., and Sherif, S. A. (2013). “A direct numerical simulation
study of higher order statistics in a turbulent round jet”. In: Phys. Fluids 25.11, p. 115102.
Townsend, A. A. (1956). The structure of turbulent shear flow. 1st ed. Cambridge University
Press Cambridge [Eng.] and New York.
Townsend, A. A. (1976). The structure of turbulent shear flow. 2nd ed. Cambridge University
Press Cambridge [Eng.] and New York.
Tweddle, I. (2003). James Stirling’s Methodus Differentialis: An Annotated Translation of Stirling’s
Text. Sources and Studies in the History of Mathematics and Physical Sciences. London:
Springer London.
Vassilicos, J. C. (2015). “Dissipation in Turbulent Flows”. In: Annual Review of Fluid Mechanics
47.1, pp. 95–114.
Wacławczyk, M., Staffolani, N., Oberlack, M., Rosteck, A., Wilczek, M., and Friedrich, R. (2014).
“Statistical symmetries of the Lundgren-Monin-Novikov hierarchy”. In: Phys. Rev. E 90.1,
p. 013022.
White, F. M. (2016). Fluid mechanics. Eighth edition. New York NY: McGraw-Hill Education.
Wosnik, M., Castillo, L., and George, W. K. (2000). “A theory for turbulent pipe and channel
flows”. In: J. Fluid Mech. 421, pp. 115–145.
Wygnanski, I. and Fiedler, H. (1969). “Some measurements in the self-preserving jet”. In: J.
Fluid Mech. 38.3, pp. 577–612.
Xu, G. G. and Antonia, R. A. (2002). “Effect of different initial conditions on a turbulent round
free jet”. In: Exp. Fluids 33, pp. 677–683.

113
A. Appendix

A.1. Estimation of the grid sizes

A.1.1. Grid size of a turbulent pipe flow

To estimate the grid size of a turbulent pipe flow, the first step involves the calculation of the
wall unit for a given Re. This begins with determining the friction factor fD , which can be
computed using Prandtl’s friction law for smooth pipes

1
= 2 log10 ( fD Re) − 0.8. (A.1)
p
√
fD

Subsequently, the friction velocity uτ is calculated using

r
fD
uτ = U , (A.2)
8 b

where Ub is the bulk velocity of the pipe which is typically set to 1 in simulations. The wall
unit is then expressed as
ruτ
r+ = . (A.3)
ν
With this, the grid spacing measured in wall units can then be set such that ∆r+ max ≤
5, ∆Rϕ+ max ≤ 5, ∆z + max ≤ 10 (see e.g. El Khoury et al., 2013). It is crucial that one grid
point approximates the distance of 20
1
δv from the wall, where

δv = ν/uτ (A.4)

represents the viscous length scale, as outlined by Pope (2000). Consequently, the grid spacing
near the wall should be approximately ∆r+ ≈ 0.05.

A.1.2. Grid size of a turbulent round jet flow

To assess the resolution of the smallest scales, we can determine the Kolmogorov length scale.
For instance, in the case of Re = 3500, the Kolmogorov length scale on the centerline can be
expressed as
 3 1/4
ν
ηK,c = = 7 · 10−4 z (A.5)
kc

115
and the dissipation of the kinetic energy scales as Equation (5.15), i.e.,
 3
3
1 k 3 1 k Bu DU 0 Bu DU 0
c = ec U z,c = e
 ≈ kc , (A.6)
z c z − z0 z4
e
z
3
where U z,c is the scaling of the centerline velocity (5.4). The dissipation |e
kc | = 0.21 has been
taken from Figure 5.16 at Re = 3500 on the centerline and the remaining parameters are
extracted from Table 5.3 and Equation (5.5) at Re = 3500. The criterion for a good resolution
of the smallest scales is
∆x
= 2.1 (A.7)
ηk
according to Pope (2000). We can see, e.g., for z = 25D the Kolmogorov length scale is
ηK = 0.02D and the grid spacing on the centerline yields

∆x, ∆y = 0.042D, (A.8)

while at z = 65D, the Kolmogorov length scale is ηK = 0.05D and the grid spacing is

∆x, ∆y = 0.11D. (A.9)

Similarly, the grid spacing for the simulation at Re = 7000 can be calculated. With |e
kc | = 0.23
in Figure 5.16 at Re = 7000 and the parameters of the DNS at Re = 7000 from Table 5.3 and
Equation (5.5), the Kolmogorov length scale on the centerline is

ηK,c = 4.510−4 z. (A.10)

The required grid spacing for z = 25D and z = 65D are ∆x, ∆y = 0.011D and ∆x, ∆y =
0.029D, respectively.

A.2. Self-preservation of additional moments

In this appendix, the remaining four moments are shown at multiple stages along the spatial
domain. The velocity moments for Re = 3500 and for Re = 7000 are showcased in Figures A.1
and A.2, respectively. Additionaly, the passive scalar moments for Re = 3500 and for Re = 7000
are presented in Figures A.3 and A.4 correspondingly. These figures complete the ten statistical
moment orders extracted from the DNS.

116
 1 1

0.8 0.8
Uzn /Uz,c

0.6 0.6
n

0.4 0.4

0.2 0.2

0 0
0 0.05 0.1 0.15 0.2 0.25 0 0.05 0.1 0.15 0.2 0.25
η η
(a) n = 3 (b) n = 5
1 1

0.8 0.8
Uzn /Uz,c

0.6 0.6
n

0.4 0.4

0.2 0.2

0 0
0 0.05 0.1 0.15 0.2 0.25 0 0.05 0.1 0.15 0.2 0.25
η η
(c) n = 7 (d) n = 9

Figure A.1.: The radial profiles of the nth axial moment normalized with the scalings in (6.25) at different distances
from the orifice for the Re = 3500 case: z/D = 25 ( ), 35 ( ), 45 ( ), 55 ( ), 65 ( ).
The black solid lines indicate the Gaussian behavior from (6.36) using γn from (6.37).

117
 1 1

0.8 0.8
Uzn /Uz,c

0.6 0.6
n

0.4 0.4

0.2 0.2

0 0
0 0.05 0.1 0.15 0.2 0.25 0 0.05 0.1 0.15 0.2 0.25
η η
(a) n = 3 (b) n = 5
1 1

0.8 0.8
Uzn /Uz,c

0.6 0.6
n

0.4 0.4

0.2 0.2

0 0
0 0.05 0.1 0.15 0.2 0.25 0 0.05 0.1 0.15 0.2 0.25
η η
(c) n = 7 (d) n = 9

Figure A.2.: The radial profiles of the nth axial moment normalized with the scalings in (6.25) at different distances
from the orifice for the Re = 7000 case: z/D = 25 ( ), 35 ( ), 45 ( ), 55 ( ), 65 ( ).
The black solid lines indicate the Gaussian behavior from (6.36) using γn from (6.37).

118
 1.2

0.8
Θm /Θm
c

0.6

0.4

0.2

0
0 0.1 0.2 0.3 0 0.1 0.2 0.3
η η
(a) m = 3 (b) m = 5

1.2

0.8
Θm /Θm
c

0.6

0.4

0.2

0
0 0.1 0.2 0.3 0 0.1 0.2 0.3
η η
(c) m = 7 (d) m = 9

Figure A.3.: The radial profiles of the mth axial moment normalized with the scalings in (7.17) at different distances
from the orifice for the Re = 3500 case: z/D = 25 ( ), 35 ( ), 45 ( ), 55 ( ). The black
solid lines indicate the Gaussian from Equation (7.26) using γm from Equation (7.27).

119
 1.2

0.8
Θm /Θm
c

0.6

0.4

0.2

0
0 0.1 0.2 0.3 0 0.1 0.2 0.3
η η
(a) m = 3 (b) m = 5

1.2

0.8
Θm /Θm
c

0.6

0.4

0.2

0
0 0.1 0.2 0.3 0 0.1 0.2 0.3
η η
(c) m = 7 (d) m = 9

Figure A.4.: The radial profiles of the mth axial moment normalized with the scalings in (7.17) at different distances
from the orifice for the Re = 7000 case: z/D = 30 ( ), 40 ( ), 50 ( ), 60 ( ). The black
solid lines indicate the Gaussian from Equation (7.26) using γm from Equation (7.27).

120
A.3. Symmetry reduction of the multi-point scalar equation

This appendix shows a step-by-step derivation of the reduced MPSE in cylindrical coordinate
form. For this, the MPSE (3.38) is transformed to cylindrical coordinates while considering
Re−1 , P e−1 = 0 has been employed, receiving
 
∂HΘ{m} 1 ∂r(1) HΘ{m+1} [Θ(m+1) 7→r(1) ] x(m+1) 7→ x(1)
+
∂t r(1) ∂r(1)
 
∂HΘ{m+1} [Θ(m+1) 7→z(1) ] x(m+1) 7→ x(1)
+
∂z(1)
m
"  
X 1 ∂r(l) HΘ{m+1} [Θ(m+1) 7→r(l) ] x(m+1) 7→ x(l)
+
r(l) ∂r(l)
l=2
 
1 ∂HΘ{m+1} [Θ(m+1) 7→ϕ(l) ] x(m+1) 7→ x(l)
+
r(l) ∂ϕ(l)
 #
∂HΘ{m+1} [Θ(m+1) 7→z(l) ] x(m+1) 7→ x(l)
+ = 0, (A.11)
∂z(l)

where the reference point is

T
x(1) = r(1) , 0, z(1) (A.12)

and

T
x(l) = r(l) , ϕ(l) , z(l) (A.13)

for l = 2, . . . , ∞. For the sake of completeness, the scaling of the time t shall be considered.
Using the scaling symmetry (7.6), the invariant surface condition (7.9) extends to

dt dr dz
= = = .... (A.14)
aSt t aSx r aSx z + az

By considering the symmetry breaking (7.14) and integrating the first and third term of (A.14),
the scaling law of the time t yields

t
t̃ = 1 ∗ . (A.15)
(z − z0 )2+ 2 aSs

In the next step, x is reduced by one independent variable by considering the scaling law
(7.15), leading to the reference point
T
r(1)


T
x̃(1) = η(1) , 0, 0 = , 0, 0 (A.16)
z(1) − z0

and the points x̃(l)

T
r(l) z(l)


T
x̃(l) = η(l) , ϕ̃(l) , z̃(l) = , ϕ̃(l) , (A.17)
z(1) − z0 z(1) − z0

121
for l = 2, . . . , ∞. Finally, the scaling laws (A.15) and (7.15)-(7.17) derived through Lie
symmetry analysis are deployed into (A.11), obtaining the reduced MPSE
m
"  
∂ H̃Θ{m} X 1 ∂η(l) H̃Θ{m+1} [Θ(m+1) 7→r(l) ] x̃(m+1) 7→ x̃(l)
+
∂ t̃ η(l) ∂η(l)
l=1
 #
∂ H̃Θ{m+1} [Θ(m+1) 7→z(1) ] x̃(m+1) 7→ x̃(1)
−x̃k(l)
∂ x̃k(l)
m
"  
X 1 ∂ H̃Θ{m+1} [Θ(m+1) 7→ϕ(l) ] x̃(m+1) 7→ x̃(l)
+
η(l) ∂ ϕ̃(l)
l=2
 #
∂ H̃Θ{m+1} [Θ(m+1) 7→z(l) ] x̃(m+1) 7→ x̃(l)
+
∂ z̃(l)
   
1
− (m + 1) 1 + a∗Ss − a∗Ss H̃Θ{m+1} [Θ(m+1) 7→z(1) ] x̃(m+1) 7→ x(1) = 0 (A.18)
 
2

and x̃k(l) 6= ϕ̃(l) .

122