University of Southern California Dissertations and Theses

Tracking and evolution of compressible turbulent flow structures

(USC Thesis Other)

Tracking and evolution of compressible turbulent flow structures

PDF

Download a page range

Download transcript

Copy asset link

Request this asset

Transcript (if available)

Content Tracking and evolution of compressible turbulent ow structures by Jonas Buchmeier A Dissertation Presented to the FACULTY OF THE GRADUATE SCHOOL UNIVERSITY OF SOUTHERN CALIFORNIA In Partial Fulllment of the Requirements for the Degree DOCTOR OF PHILOSOPHY (Aerospace Engineering) May 2022 Copyright 2022 Jonas Buchmeier Acknowledgements First of all, I would like thank my advisor, Prof. Iv an Bermejo-Moreno, for all his guidance and patience throughout my PhD journey. I felt in good hands while working with Iv an. My PhD work builds up on his own PhD on the geometry of turbulence and would not have been possible without it. Also I want to thank him for taking me to the Summer Research Program of the Center for Turbulence Research at Stanford University in 2018, which was quite a special experience and honor. I would like to thank Profs. Julian Domaradzki, Carlos Pantano-Rubino, Aiichiro Nakano, Paul Newton and Mihailo Jovanovic for agreeing to be in my qualifying exam and defense committee and showing interest in this research. I also want to express my thanks to Prof. Stefan Hickel for his advice during my Master thesis at the Technical University of Munich and Dr. Bernhard Eisfeld of the German Aerospace Center for awakening my passion for uid dynamics during my Bachelor thesis project. My deepest gratefulness belongs to my parents, Harald and Ulla, for their endless support and encour- agement throughout this journey. I would like to thank my lab mates Xiangyu, Jon and Naili and also my oce mates Arturo, Raye, Ruiyang and Keyvan for all the fun we had in and outside the oce. Special thanks goes to Alex for the great collaboration during his stay at USC and his friendship. This work would not have been possible without nancial support received from the Army Research Oce (Grant W911NF2010096) and computational resources provided by USC's Center for Advanced Research Computing and an award of computer time provided by the Innovative and Novel Computational Impact on Theory and Experiment (INCITE) program of the Argonne Leadership Computing Facility, which is a DOE Oce of Science User Facility supported under Contract DE-AC02-06CH11357. ii Table of Contents Acknowledgements ii List of Tables v List of Figures vi Abstract ix Chapter 1: Introduction 1 1.1 Motivation . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 1 1.2 Survey of ow feature tracking methods . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 3 1.2.1 Implicit feature extraction and tracking . . . . . . . . . . . . . . . . . . . . . . . . . . 3 1.2.2 Explicit feature extraction and correspondence based tracking . . . . . . . . . . . . . . 5 1.3 Outline of the present work . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 8 Chapter 2: Flow feature tracking methodology 10 2.1 Introduction . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 10 2.2 Event denitions in correspondence search . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 12 2.3 Structure attributes for event detection and rejection . . . . . . . . . . . . . . . . . . . . . . . 14 2.4 Detection of correspondences . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 16 2.4.1 Projection of spatial attributes . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 17 2.4.2 Attribute scaling . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 18 2.4.3 Nearest-neighbor search in the attribute space . . . . . . . . . . . . . . . . . . . . . . 18 2.4.4 Radius search in the spatial subspace . . . . . . . . . . . . . . . . . . . . . . . . . . . 21 2.4.5 Detection of unmatched structures . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 24 2.4.6 Collecting correspondences form the multi-stage search . . . . . . . . . . . . . . . . . . 24 2.5 Constraining of correspondences . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 24 2.5.1 Bounding volume constraints . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 25 2.5.2 Surface proximity constraints . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 29 2.6 Detection of events . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 30 2.7 Constraining and optimization of events . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 30 2.7.1 Event-driven correspondence constraining . . . . . . . . . . . . . . . . . . . . . . . . . 31 2.7.2 Combinatorial extraction and conservation-based constraining of subevents . . . . . . 32 2.7.3 Sequence of event realizability constraints . . . . . . . . . . . . . . . . . . . . . . . . . 34 2.8 Treatment of periodic boundaries . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 34 2.9 Graph representation of temporal evolution of ow features . . . . . . . . . . . . . . . . . . . 36 2.10 Graph mining for pattern recognition in ow feature evolution . . . . . . . . . . . . . . . . . 42 2.10.1 Criteria for pattern evaluation . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 46 2.10.2 Application to a synthetic tracking graph . . . . . . . . . . . . . . . . . . . . . . . . . 48 iii Chapter 3: Passive scalar structures in isotropic turbulence and shock-turbulence interac- tion 54 3.1 Introduction . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 54 3.2 Numerical setup and initialization of passive scalar structures . . . . . . . . . . . . . . . . . . 55 3.3 Downstream distribution of events . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 63 3.4 Structure volumes of event types . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 67 3.5 Split and merge events in geometrical feature space . . . . . . . . . . . . . . . . . . . . . . . . 69 3.6 Downstream evolution surface-averaged physical quantities . . . . . . . . . . . . . . . . . . . 72 3.7 Correlations between local surface geometry and physical quantities on the surface . . . . . . 77 3.8 Trajectories of primary structures in the feature space . . . . . . . . . . . . . . . . . . . . . . 80 3.8.1 Geometrical trajectories and surface area evolution . . . . . . . . . . . . . . . . . . . . 80 3.8.2 Trajectories with scalar dissipation and vorticity . . . . . . . . . . . . . . . . . . . . . 86 Chapter 4: Vortical structures in temporally developing compressible mixing layers 92 4.1 Introduction . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 92 4.2 Problem formulation and numerical setup . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 93 4.2.1 Setup of the ow conguration . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 93 4.2.2 Setup of the tracking algorithm . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 95 4.3 Validation of self-similarity and growth rate of the mixing layer . . . . . . . . . . . . . . . . . 97 4.4 Identication, extraction and general statistics of vortical structures . . . . . . . . . . . . . . 102 4.5 Creations neighboring existing vortices . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 106 4.6 Detection of hairpin-like vortical structures . . . . . . . . . . . . . . . . . . . . . . . . . . . . 108 4.7 Hairpins in the geometrical feature space . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 113 4.8 Temporal distribution and formation of hairpin-like vortices . . . . . . . . . . . . . . . . . . . 115 4.9 Cross-stream and spanwise growth of hairpin vortices . . . . . . . . . . . . . . . . . . . . . . . 123 Chapter 5: Conclusion 128 Bibliography 131 iv List of Tables 3.1 Initial conguration of passive scalar structures used in the DHIT and STI conguration . . . 57 3.2 Settings for NN and radius search and correspondence constraining in passive scalar application 62 3.3 Settings of event constraining in passive scalar application . . . . . . . . . . . . . . . . . . . . 63 3.4 Settings for graph building in passive scalar application . . . . . . . . . . . . . . . . . . . . . 63 3.5 Mean disappearance location, start location and spatial length of split cascade of primary structures . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 84 4.1 Grid properties for mixing layer congurations. . . . . . . . . . . . . . . . . . . . . . . . . . . 95 4.2 Settings for the correspondence search and constraining. . . . . . . . . . . . . . . . . . . . . . 97 4.3 Settings for event constraining and event building in the graph. . . . . . . . . . . . . . . . . . 97 4.4 Begin, end and growth rate of the self-similar regimes. . . . . . . . . . . . . . . . . . . . . . . 97 4.5 Input parameter of the hairpin lter algorithm . . . . . . . . . . . . . . . . . . . . . . . . . . 111 v List of Figures 1.1 Classication of ow feature tracking methods . . . . . . . . . . . . . . . . . . . . . . . . . . . 4 2.1 Overview of the tracking algorithm . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 11 2.2 Generic shape enclosed by AABB and OBB and illustration of radius search . . . . . . . . . . 22 2.3 Illustration of AABB volume and distance constraint . . . . . . . . . . . . . . . . . . . . . . . 26 2.4 Illustration of separating axis theorem and local surface proximity . . . . . . . . . . . . . . . 28 2.5 Duplication of boundary intersecting or crossing structures . . . . . . . . . . . . . . . . . . . 35 2.6 Event reconnection under periodic boundary conditions . . . . . . . . . . . . . . . . . . . . . 36 2.7 Graph representation of structure evolution . . . . . . . . . . . . . . . . . . . . . . . . . . . . 37 2.8 Vertex clustering strategies for graph building . . . . . . . . . . . . . . . . . . . . . . . . . . . 39 2.9 Classication of branches in graph . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 41 2.10 Synthetic graph for graph mining demonstration . . . . . . . . . . . . . . . . . . . . . . . . . 43 2.11 Labeled synthetic input graph . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 45 2.12 Illustration of graph compression and graph isomorphism . . . . . . . . . . . . . . . . . . . . 47 2.13 Illustration of graph compression and generation of parallel edges . . . . . . . . . . . . . . . . 47 2.14 Compressed synthetic graph after the rst Subdue iteration . . . . . . . . . . . . . . . . . . . 50 2.15 Compressed input graph after the 8-th iteration . . . . . . . . . . . . . . . . . . . . . . . . . . 51 2.16 Best pattern extracted from a synthetic input graph . . . . . . . . . . . . . . . . . . . . . . . 52 2.17 Embedding of pattern instances . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 52 2.18 Embedding and decompression of pattern instances . . . . . . . . . . . . . . . . . . . . . . . . 53 2.19 Embedding and decompression of pattern instances of varying compressed nodes . . . . . . . 53 3.1 Initial radial scalar prole and regions of scalar discontinuities . . . . . . . . . . . . . . . . . . 57 3.2 Renderings of passive scalar structure evolution in DHIT . . . . . . . . . . . . . . . . . . . . . 58 3.3 Plane cuts of transitional regions of scalar gradient . . . . . . . . . . . . . . . . . . . . . . . . 60 3.4 Problem setup of passive scalar mixing in STI and blended DHIT precursor simulations . . . 61 3.5 Number of extracted iso-surfaces depending on iso-value . . . . . . . . . . . . . . . . . . . . . 62 3.6 Hierarchical layout of the graph of the 4 scalar in DHIT . . . . . . . . . . . . . . . . . . . . 64 3.7 Hierarchical layout of the graph of the 2 scalar in DHIT . . . . . . . . . . . . . . . . . . . . 64 3.8 Downstream distribution of dierent event types . . . . . . . . . . . . . . . . . . . . . . . . . 65 3.9 Representative compound events of 2 structures in the post-shock region . . . . . . . . . . 66 vi 3.10 Volume histogram of disappearing (a) and 4 (b) structures in DHIT and STI M = 1:5 . 68 3.11 Joint histogram of volumes V s of source and V t targets of continuation events . . . . . . . . . 69 3.12 Joint histogram of volumes V s of source and V t of summed targets of split events . . . . . . . 70 3.13 Joint histogram of volumes of summed source and summed targets of compound events . . . 71 3.14 Joint histograms of geometrical attributes of source and target structures of split events . . . 71 3.15 Absolute shape index mapped on 2 structures in DHIT during split events . . . . . . . . . 72 3.16 Joint histograms of geometrical attributes of source and target structures of merge events . . 73 3.17 Ensemble mean of the downstream evolution of the scalar gradient magnitude, vorticity mag- nitude and the alignment between the scalar gradient and vorticity averaged on isosurfaces. . 74 3.18 Ensemble mean of the downstream evolution of the scalar gradient alignment with strain-rate eigenvectors averaged on scalar iso-surfaces . . . . . . . . . . . . . . . . . . . . . . . . . . . . 76 3.19 Dimensionless curvedness, absolute shape index and the alignment between the eigenvector of the strain-rate tensor and the scalar gradientr 4 mapped on an isosurface of the 4 . . 78 3.20 Alignment analysis conditioned on local surface geometry . . . . . . . . . . . . . . . . . . . . 79 3.21 Trajectories in the geometric feature space of individually tracked primary 4 structures in DHIT . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 82 3.22 Early and later higher-order split of 4 structures in DHIT . . . . . . . . . . . . . . . . . . . 82 3.23 Ensemble mean trajectories of primary structures in geometric space for DHIT and STI . . . 83 3.24 Temporal evolution of the normalized surface area of primary structures in DHIT. . . . . . . 84 3.25 Streamwise evolution of the surface area of primary and 2 structures in DHIT and STI normalized by its initial DHIT value . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 85 3.26 Ensemble mean trajectories of primary structures in geometric space with scalar gradient and vorticity magnitude for DHIT . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 87 3.27 Ensemble mean trajectories of primary structures in geometrical feature space with scalar gradient and vorticity magnitude (STI M = 1:5). . . . . . . . . . . . . . . . . . . . . . . . . 88 3.28 Ensemble mean trajectories of primary structures in geometrical feature space with scalar gradient and vorticity magnitude (STI M = 3:0). . . . . . . . . . . . . . . . . . . . . . . . . . 89 3.29 Renderings of primary 4 and structures in DHIT . . . . . . . . . . . . . . . . . . . . . . 90 3.30 Renderings of 4 and in STI M = 3:0 . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 91 4.1 Illustration of the initial condition of the temporally developing mixing layer . . . . . . . . . 94 4.2 Grid specication of the mixing layer congurations . . . . . . . . . . . . . . . . . . . . . . . 96 4.3 OBB tree presentation of a hairpin vortex . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 96 4.4 Density xy-planes of the M c = 0:7 mixing layer at dierent times . . . . . . . . . . . . . . . . 98 4.5 Momentum thickness evolution and growth rate . . . . . . . . . . . . . . . . . . . . . . . . . . 99 4.6 Self-similar streamwise velocity proles . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 100 4.7 Individual streamwise turbulence intensity proles in the cross-stream direction . . . . . . . . 101 4.8 Ensemble mean prole of streamwise, cross-stream and spanwise turbulence intensities and Reynolds shear stress . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 101 4.9 Iso-value evolution and number of extracted structures from the mixing layer . . . . . . . . . 103 vii 4.10 Detected event types of Q structures in the mixing layer over time . . . . . . . . . . . . . . . 103 4.11 Temporal evolution of the geometrical signature in the early stage of the M c = 0:3 mixing layer104 4.12 Selection of -vortices . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 105 4.13 Alignment between Q-surface normal and vorticity. Cross-stream locations of creation events 106 4.14 Geometrical variety of the extracted Q-structures in the self-similar regimes . . . . . . . . . . 107 4.15 Newly created structures and their rst event in geometric feature space . . . . . . . . . . . . 107 4.16 Creation-merge pattern in the self-similar regime of the M c = 0:3 mixing layer . . . . . . . . 108 4.17 Creation and subsequent merge bending a tube-like structures . . . . . . . . . . . . . . . . . . 108 4.18 Creation and subsequent merge forming legs on tube-like structures . . . . . . . . . . . . . . 109 4.19 Velocity eld surrounding an individual hairpin vortex . . . . . . . . . . . . . . . . . . . . . . 110 4.20 Large connected hairpins and shape index on individual hairpin . . . . . . . . . . . . . . . . . 110 4.21 Elements of the hairpin lter algorithm . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 112 4.22 Structures geometrically similar, but physically dissimilar to a hairpin vortex. . . . . . . . . . 113 4.23 Leg vorticity component - cross coordinate clustering rejections during hairpin ltering. . . . 114 4.24 Examples of detected hooks and arches. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 115 4.25 Dimensionless curvedness mapped on hairpins . . . . . . . . . . . . . . . . . . . . . . . . . . . 115 4.26 Hairpin vortices in the geometrical feature space . . . . . . . . . . . . . . . . . . . . . . . . . 116 4.27 Geometrical signatures of hairpin-like vortices grouped by mixing layer stage . . . . . . . . . 117 4.28 Temporal distribution of hairpin detections . . . . . . . . . . . . . . . . . . . . . . . . . . . . 118 4.29 Connected structures breaking into individual hairpin-like structures . . . . . . . . . . . . . . 120 4.30 Hooks developing hairpin shape . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 120 4.31 Short-lived hairpin in the M c = 0:3 mixing layer . . . . . . . . . . . . . . . . . . . . . . . . . 121 4.32 Long-lived hairpin in the M c = 0:3 mixing layer . . . . . . . . . . . . . . . . . . . . . . . . . . 122 4.33 Hairpin in geometrical feature space enhanced by surface-averaged streamwise vorticity com- ponent . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 123 4.34 Lifetime of a structure developing temporarily hairpin-like shape ( = 280 - 376) . . . . . . . 124 4.35 Lifetime of a structure developing temporarily hairpin-like shape ( = 96 - 172) . . . . . . . . 125 4.36 Mean cross stream locations of hairpins . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 126 4.37 Cross-stream extent of center plane crossing hairpins . . . . . . . . . . . . . . . . . . . . . . . 126 4.38 Cross-stream and spanwise growth of hairpin vortices . . . . . . . . . . . . . . . . . . . . . . . 127 viii Abstract An algorithm is introduced to track the evolution of individual structures in turbulent ows. Flow structures are dened as closed surfaces of any kind in numerical or experimental data sets. To track evolving structures over time and to capture their interactions with each other, correspondences between structures in consecutive time steps are found in a higher-dimensional feature space, which is constructed from dierential geometric and spatial structure attributes. After detected correspondences are tested for physical realizability, structure interactions are mapped into a graph data structure which represents the evolution and interactions of all tracked structures over time. The hybrid attribute- and region-based algorithm allows large sampling intervals in the database by capturing complex compound interactions resulting from fast feature dynamics. The tracking graph is queried to retrieve individual and statistical information on the evolution of the structures. Graph data mining techniques are adapted to recognize common patterns in the ow feature evolution. In a rst application, the tracking algorithm is utilized to conduct a structure-based numerical analysis of passive scalar mixing in decaying homogeneous isotropic turbulence (DHIT) and shock-turbulence inter- action (STI) canonical congurations. Structures of a scalar eld are dened as disconnected isosurfaces of the particular underlying scalar eld in shock-capturing direct numerical simulations (DNS) of DHIT, with Taylor microscale Reynolds number Re = 40 and turbulence Mach number M t = 0:2, and STI cases in which the turbulence is processed by a shock wave at Mach numbers M = 1:5 and 3:0. The analysis focuses on the temporal evolution of ensembles of passive scalar structures, initialized as collections of uniformly spaced spheres of varying scales, commensurate with the initial Taylor microscale, , of the underlying tur- bulence. Changes in the feature surface geometry are related to the underlying physical processes relevant to the mixing (turbulent transport, scalar dissipation, and shock compression). Temporal surface convolution increases for initially larger structures, resulting in a higher probability of locally hyperbolic geometries where break-up into smaller structures occurs. Shock-induced deformation of the structures amplies breakup pro- cesses, enhancing mixing, particularly for larger structures. Mixing enhancement by the shock is manifested as an amplication of the surface-averaged scalar gradient, which increases for initially larger structures. The alignment between the scalar gradient and the most extensional strain-rate eigendirection on the scalar isosurfaces also increases across the shock. Larger magnitudes of the scalar gradient and its alignment with ix the most compressive strain-rate eigendirection correlate with atter surface regions. Shock-induced struc- ture compression increases the area coverage of at regions, where the amplication of scalar gradient is localized. In a second application of the tracking algorithm, the inhibiting eect of compressibilty on the turbulence intensities in temporarily developing mixing layers is analyzed from a structural viewpoint. In the congu- rations of initial momentum thickness Reynolds number Re ;0 = 160 and varying convective Mach numbers M c =f0:3; 0:7; 1:1g the dynamics and geometrical properties of vortical structures dened by iso-surfacing the Q-criterion eld are investigated. Using a time-dependent iso-value based on statistics of the Q-eld, less vortices are extracted with increasing convective Mach number, expressing the inhibiting eect of compress- ibilty on the structural level. Higher compressibilty lowers the curvedness of the extracted Q-structures. The creation of vortices in the neighborhood of existing vortices and the subsequent merging of both vortices is recognized as a prominent pattern of the vortex dynamics. Other patterns of statistical signicance are split-dominated or represent the short-lived evolution of isolated structures which disappear shortly after their creation without interacting with others. To evaluate the occurrence, persistence and signicance of hairpin-like vortices in the mixing layers, a detection method for such special vortices is introduced. For all convective Mach numbers, the number of hairpin detections peaks at Re 400 in the pre-self-similar stages of the mixing layers and decays afterwards, suggesting that hairpins are not a dominant ow feature once the turbulent background is developed. The time span in which vortices preserve their hairpin-like shape varies, but most hairpins are rather short-lived as they cross the center planes of the mixing layers and are partially pulled into the opposite free streams. The majority of individual hairpins is formed by the structural breakup process of largely connected hairpins which originate from vortices which grow in the early stage of the mixing layer through merging of tube-like vortices. A number of hairpins have a tilted orientation with respect to the center plane, especially in the stages after the pre-self-similar regime. For hairpins that have quasi-streamwise orientation cross-stream and spanwise growth of equal rates occurs simultaneously. x Chapter 1 Introduction 1.1 Motivation Despite its chaotic, multi-scale character, organized motions and turbulence structures have been observed in turbulent uid ow for a long time. Hussain (1986) links the motivation to study coherent structures to a desire to nd order in apparent disorder. Although ambiguous denitions of what a coherent structure is can be found in the literature, there is an intuitive understanding that coherent structures are associated with motions that have some sort of spatial organization and temporal coherence leading to a certain lifetime of the coherent structure in which its appearance actively or passively evolves. Complementing statistical approaches to study the nature of turbulent ow, structural analyses have provided fundamental insight into the mechanisms of momentum and energy transfer in turbulence. For example, in the buer layer of near-wall turbulence streaks of streamwise velocity are found to carry most of the turbulent kinetic energy (TKE) (Kline et al., 1967; Smith & Metzler, 1983) and quasi-streamwise vortices organize dissipation and Reynolds stresses (Jim enez, 2013), whereas in the logarithmic layer, most of the momentum transfer origi- nates from wall-attached eddies (Townsend, 1961) and vortices organize into self-similar clusters (Del Alamo et al., 2006). Structural analyses of turbulent ows have led to the development of structure-based models of turbulence ne scales (Townsend, 1951; Tennekes, 1968; Lundgren, 1982; Pullin & Saman, 1993) and subgrid-scale models for large eddy simulations (Misra & Pullin, 1997). To gain further insight into the structural mechanisms of turbulent ow, dierent description techniques for the geometrical properties of turbulence structures have been suggested in the literature (Bermejo-Moreno & Pullin, 2008; Leung et al., 2012). By utilizing these techniques the geometry of an individual turbulence structures can be related to the underlying ow physics. The structural approach can be expanded from turbulence structures to other ow features such as droplets in multiphase ows, transition zones or shock waves. Silver (1995) and Post et al. (2002) summarize dierent techniques to dene, extract and visualize features from a ow eld. Since ow features are an abstraction of the original data and typically occupy a relatively small region of the computational domain, their extraction in general leads to a reduction of data, as emphasized, for example, by Van Walsum et al. (1996). Extracted features can be isolated from the 1 overall data set, which simplies their individual exploration compared to visualizing the entire data set at once. Time-resolving computational and experimental techniques to study uid motion, such as direct numeri- cal simulation (DNS) or time-resolved three-dimensional particle image velocimetry (3D-PIV), have granted access to instantaneous, three-dimensional ow elds and hence to the dynamics of turbulence structures. Whereas the aforementioned geometrical techniques describe turbulence structures for instantaneous snap- shots of a given data set, the present work focuses on the evolution of those geometrical properties over time and the relation to the temporally or spatially evolving physics of the ow eld. Therefore, a novel method to track individual ow structures throughout the ow eld is introduced in §2 adding temporal coherence to the geometrical analysis. The method enables simultaneous tracking of intertwined structures with signicant concave regions over large sampling intervals, by combining region- and attribute-based algo- rithms. The tracking attributes incorporate geometric properties of the structures obtained from dierential geometry analyses. In our present tracking methodology, ow features are dened as closed surfaces of any kind embedded in a three-dimensional space, but could be generalized to other dimensions (e.g., ow features in two-dimensional turbulence). In principle, the surfaces could be derived from a numerical simulation, as well as originate from experimental data. In this work, the surfaces are extracted from volumetric scalar elds by isosurfacing. However, we emphasize that this ow feature extraction is extrinsic to the tracking itself, as the tracking is explicit, see §1.2.2. The algorithm tracks individual, disconnected isosurface components of the chosen eld. The polygonal mesh discretization of the isosurface is obtained by applying the marching cube algo- rithm (Lorensen & Cline, 1987) followed by a disconnection step to isolate individual isosurfaces. Features, structures and isosurfaces are used synonymously in this work. The shape of the structures, characterized from their dierential geometry properties, plays a key role in the presented tracking methodology. In addi- tion to its contribution to the tracking eort via geometrical attributes used in the correspondence search, we aim to relate the temporal development of the geometry of the structures to the underlying ow physics. The method is designed to track all structures extracted from the ow simultaneously over large sampling intervals, which pose challenges to existing tracking methods due to the increased inter-frame evolution. To study the evolution of individual structures and their mutual interactions, the tracking results are encoded in a graph data structure, which is a representation of the evolution of all tracked features. In a rst application presented in this work the proposed tracking method is used to investigate passive scalar mixing in compressible turbulence from a structural point of view. Adding to ndings based on statistical analyses, the tracking of individual nite-sized scalar structures enables to investigate how the scalar mixing is represented at the structural level in the form of the structural breakup process and the structure lifetime until complete diusion. Physical quantities relevant to the mixing, such as the magnitude of the scalar gradient and its alignment with eigendirections of the strain-rate tensor, can be mapped on the tracked surfaces and hence are also examined from a structural viewpoint. The tracking application and 2 geometrical analysis allows to study how the shock-induced structure deformation amplies the structural mixing process. In particular, the temporal geometrical analysis of the scalar structures permits to identify the predominant shapes of increased scalar dissipation as the structures pass through the shock wave. Another ow conguration addressed in this work is the mixing layer formed between two parallel streams. Statistical analyses show an inhibiting eect of compressibilty on the turbulence intensities in the mixing layer (Papamoschou & Lele, 1993). By utilizing the proposed tracking algorithm the eect of compressibility investigated in regard to the dynamics and geometrical properties of individual vortical structures in the mixing layer. Despite of their proposal in the 1950s, the occurrence, persistence and relevance of hairpin vortices in shear ow congurations (such as the mixing layer) is still subject to ongoing debate. Therefore, a lter is proposed to detect hairpin-like structures during the development of the mixing layer. After identifying structures which are temporarily of hairpin-like shape, their individual evolution and geometry is examined. The tracking allows to study how hairpin vortices are formed and how they persist in the developing turbulent background of the mixing layer. 1.2 Survey of ow feature tracking methods The features of a time-depending ow eld may have evolving characteristics such as their spatial location and shape. Analyzing the dynamics of ow features in physical space is subject to feature tracking. Dierent feature tracking techniques have been introduced in the literature that dier on how the features are dened and extracted from the ow eld. A classication of feature tracking algorithms is presented in gure 1.1. In explicit feature tracking the extraction of features is conducted as a pre-step to the tracking itself. Because the extraction is performed independently for each time instant (frame) and the temporal coherency of each feature is not utilized in the extraction process, the correspondence problem between features in successive frames arises. Solving the correspondence problem means nding the relation between feature in successive frames. For example, if two features extracted from successive frames are matched in a one-to-one corre- spondence, the two features represent the same feature at the two dierent time instants. Without this correspondence information the individual feature evolution cannot be analyzed. In general, each feature can have multiple forward and backward correspondences which adds signicant complexity to the corre- spondence problem. Avoiding the correspondence analysis between frames is one motivation for implicit feature tracking methods, which either apply a spatio-temporal feature extraction rather than a pure spatial extraction (as in explicit tracking) or use a prediction in the extraction process based on features extracted in previous frames, and hence also include temporal coherence in the feature extraction process. 1.2.1 Implicit feature extraction and tracking Weigle & Banks (1998) and Ji et al. (2003) present implicit algorithms to track time-varying isosurfaces and interval volumes by isocontouring in higher dimensions across the entire spatio-temporal domain. Assuming 3 Feature tracking methods Implicit extraction & tracking Explicit extraction Correspondence based tracking Spatiotemporal isocountouring & feature ow elds Path predictive Physics based Attribute based Region based Physics based Hybrid Figure 1.1: Classication of ow feature tracking methods and proposed hybrid approach of attribute and region based tracking. that correspondingR 3 isosurfaces from two consecutive time steps overlap with each other, isocontouring in R 4 is used to nd correspondences implicitly. Given aR 3 isosurface at time t n the overlapping isosurface at timet n+1 is found by computing the connected isosurface inR 4 that intersects with the givenR 3 isosurface and then slicing it along the time dimension to get the overlapping R 3 isosurface at t n+1 . Hence the correspondence problem for isosurfaces in R 3 is implicitly solved. To reduce the computational cost of isocontouring inR 3 andR 4 during the tracking, in a follow-up work, Ji & Shen (2004) used a precomputed correspondence table between isosurface components for any possible isovalue. Then, in the actual tracking, the correspondence search reduces to look-up operations from this pre-computed table. The construction of the look-up table is based on the observation that the correspondence relationship between R 3 isosurfaces in consecutive time steps can only change at critical isovalues. At critical isovalues the number of isosurface components inR 4 changes, since the temporal correspondences change. The tracking method proposed by Bajaj et al. (2002) is based on generating the isosurfaces in consecutive time steps progressively from the previous time step by performing a temporal contour propagation of existing isosurfaces and seed generation of new components, instead of extracting isosurfaces independently from dierent time steps. Sohn & Bajaj (2005) construct the topology of time-varying isosurfaces resulting in a time-dependent contour tree. Correspondences are computed between two instants of the contour tree. Topological events, such as the split of one isosurface into multiple ones or the merge of multiple isosurfaces into one, are represented in a topology-change graph. All of the aforementioned implicit spatio-temporal tracking methods are based on a feature-overlapping criterion, equivalently found in a certain class of explicit methods, as described later. Theisel & Seidel (2003) introduced an alternative approach to represent the time-varying behavior of features not as a higher- dimensional isosurface but as stream lines or surfaces of higher-dimensional vector elds. This approach of \feature ow elds" is applicable to local ow feature and was tested, in particular, to track critical points of 2D and 3D vector elds. Even though in feature ow elds no explicit correspondence matching is performed, features may not be connected over time by the streamlines of the feature ow eld, if the spatial displacement between the features in consecutive time steps is too large and, therefore, an overlap criterion is somewhat implied. Theisel et al. (2005) apply the concept of `feature ow elds' to extract and 4 track vortex core lines in a formulation based on `parallel vectors' as introduced by Peikert & Roth (1999). The concept of parallel vectors derives two vector elds from the original vector eld, such as itself and its curl, and a vortex core is located where these two vector elds are parallel. Streamline/surface integration of the feature ow eld is used to extract/track the parallel vector lines. Another implicit tracking method for line-type features was proposed by Bauer & Peikert (2002) to track vortex cores through dierent spatial and temporal scales. Rather than utilizing a spatio-temporal extraction process, Banks & Singer (1995) apply a vorticity- predictor pressure-corrector scheme to reconstruct and implicitly track elongated vortex tubes in turbulent ow. Muelder & Ma (2009) also propose an implicit tracking method in a prediction-correction manner. Assuming the features have fairly predictable paths, prediction methods of dierent order are used to estimate the feature region in a frame based on the previous frames. In a correction step the estimated feature region is either shrunk or grown to extract the actual feature in the current frame. Hence, the method assumes an overlap between the estimated and actual feature region. In general, low temporal resolution in the ow eld constitutes a challenge to tracking its features. To track volumetric ow features, especially over large sampling intervals, Sauer et al. (2014) proposed a method using trajectories of corresponding Lagrangian particle data to predict the location of the feature in the next available time instant. Another particle-based, path-predictive tracking method was applied by Schafhitzel et al. (2008) on segmented vortex structures to study perturbation events on arbitrary temporarily tracked vortex structures. The vortex core line is tracked in a predictor-corrector manner by placing particles on the core line. The particles are integrated along a given unsteady vector eld and tested against overlapping with core lines in the consecutive time step. 1.2.2 Explicit feature extraction and correspondence based tracking In contrast to the described implicit tracking methods, the feature extraction step is a pre-step to the actual tracking eort in explicit tracking. As mentioned earlier, in explicit tracking the correspondence problem arises due to the exclusion of temporal coherency in the feature extraction. Explicit tracking is therefore also called correspondence-based tracking and related methods can be divided in three categories based on how the search for correspondences is conducted: physics-, region- and attribute-based correspondence search. Another way of classication is distinguishing between local and global tracking techniques as described later on. Since the features are extracted from the ow eld, their dynamics are in principle governed by the same equations as the uid motion. Clyne et al. (2012) proposed a feature path prediction method which is based on the physics of the underlying solutions to the equation of uid motion. To start the path prediction, for each feature, the point of minimum dissipation is found and advected by the velocity eld to the consecutive frame. A streamline passing through this advected point is integrated in forward and backward directions. 5 Possible candidates for correspondence with the original feature are features in the consecutive frame which intersect the streamline. The best match is found by minimizing the relative volume dierence between candidate and original features. Early work on region-based correspondence search for explicit tracking was conducted by Silver & Wang (1996, 1997, 1998). The main motivation of the region-based search is based on the observation that, if the sampling frequency is high enough, corresponding features in consecutive time steps usually have a spatial overlap. To track one feature, rst, all overlapping features in the consecutive frame are found. To select the best match from those candidates a degree of overlap is computed in terms of a normalized volume- dierence test. This approach has been extended by Chen et al. (2003) to be applied to feature tracking in data sets obtained by application of adaptive mesh renement (AMR) techniques, whereby each feature can span multiple renement levels. A \feature tree" is obtained representing the feature in multiple levels of renement. A tracking method based on region-based correspondence was applied to track vortical structures in forced and decaying isotropic turbulence by Kareem et al. (2007) and Abdel-Kareem (2011). Lozano-Dur an & Jim enez (2014) use region-based correspondence search to study the evolution of coherent structures in turbulent channels, organizing the found correspondences over time into a graph that represented temporal structure interactions. Another region-based tracking method was proposed by Saikia & Weinkauf (2017), where feature regions are dened by merge trees that, in general, segment a scalar eld into dierent regions. The similarity between two subtree regions in consecutive frames is measured by considering the volume overlap and the 2 distance between voxel intensity histograms of two regions. A graph data structure is built over all time steps, where the nodes are the subtree regions in each time step and the edges represent all possible connections between regions. The edge weight is dened by the combined distance from the volume overlap and the histogram distance. Since the graph incorporates connection information for all regions over all time steps, it can be used to make tracking decisions globally instead of considering two consecutive time steps. By applying Djikstra's shortest path algorithm, the tracking path for each feature is determined. A similar approach, using a graph to obtain tracking paths, was used by Widanagamaachchi et al. (2012), which also utilizes merge trees for the denition of regions. In this method the features contained in a track (branch) of the graph have similar intensity levels, whereas in Saikia & Weinkauf (2017) the intensity level of a region of the merge tree can change over time. The same method was later adapted by Widanagamaachchi et al. (2015) to track features in embedded surfaces with application to turbulent combustion. As mentioned earlier, region-based methods rely on the spatial overlap of feature in consecutive frames. Hence, they require a fairly high sampling frequency, especially when small, fast moving features are present in the data. In contrast to the region-based feature tracking methods, attribute-based methods try to solve the correspondence problem by considering only a small set of abstract attributes describing the feature. Therefore attribute-based matching costs signicantly less memory and a correspondence search in multiple stages becomes aordable. The feature attributes are utilized to uniquely represent each feature in the search 6 for correspondences. Therefore, the choice of attributes used in this representation is fundamental and may depend on the application. To solve the correspondence problem based on an attribute representation of the features, matching criteria need to be established for the attributes of source and target feature. Samtaney et al. (1994) introduced a list of evolutionary events in the lifetime of a feature: continuation, creation, disappearance, split and merge, see §2.2. To detect those events an attribute-based method which uses the feature centroid as the primary attribute in the search is proposed. In a rst step, continuing features are determined by nding the closest feature in the consecutive frame in terms of the centroids. The matched candidates are tested with respect to additional attributes, such as second order moments of inertia, volume and mass. In a second step remaining unmatched features are tested for split and merge events. Feature evolution over multiple frames and interactions between features in the form of the detected events are represented in a directed acyclic graph. A similar tracking method based on centroid location and volume as key attributes was adapted by Chan et al. (2018) to track bubbles in breaking waves. Reinders et al. (1999) and Reinders et al. (2001) propose an attribute-based method in a prediction- verication scheme. Their attribute set contains centroid position, volume, mass and the best tting ellipsoid (in terms of position, scaling, and orientation) of each feature. The prediction utilizes a linear extrapolation using the previous two frames and is applied in forward and backward directions. Hence, multiple passes are executed for each frame. The previously described tracking methods make the common assumption that the feature motion in between consecutive frames is limited and fall into the category of local tracking techniques. This assumption is fairly obvious in terms of the region-based methods, since they rely on the spatial overlap. The attribute- based method of Samtaney et al. (1994) also tests the closest centroids in consecutive frames for continuation events. The prediction-verication scheme in Reinders et al. (2001) makes the assumption that features evolve predictably. Hence, it also implies certain restriction on the feature evolution between frames. Contrary, global tracking techniques do not impose any restrictions on feature motion between frames. In other words, features may move anywhere in the domain in between frames. Hence, these methods do not perform the tracking of each feature independently with local search approaches. In a global tracking algorithm all the features from two consecutive frames are considered together and a globally best match is pursued. The technique introduced by Ji & Shen (2006) seeks the globally best match by using the Earth Mover's Distance (EMD), which is also known as Wasserstein metric, as a matching cost measure and propose a \branch-and- bound" method to nd the global minimal cost eciently. EMD is a statistical tool used to measure the distance between two probability distributions over a region. However, although in principle no restrictions are applied, \the cost algorithm is weighted heavily in favor of features that are in close proximity to the tracked feature" as pointed out by Clyne et al. (2012). The Wasserstein metric was also used by Soler et al. (2018) to track topological features in time-varying scalar data. Caban et al. (2007) introduce a non-local tracking method which is attribute-based with no bias toward nearby features. It is based on the observation that volume segments in time-varying data have textural 7 features that are characteristic to that particular subvolume. Here a feature is associated with a subvolume of the domain. The texture of a subvolume is dened by 26 textural attributes. Correspondences are found in terms of the best texture match in subsequent time steps. First, a tracking window is estimated based on the dynamics between frames in the data set. For data sets where little is known about those dynamics, the window can be large and, in principle, span the entire time domain. Then, similarity measurements in terms of texture are conducted for candidate subvolumes. Since correspondences are found based on textural properties, the best match is not necessarily found in the features spatial neighborhood. 1.3 Outline of the present work The ow feature tracking algorithm proposed as a local, hybrid approach combining attribute- and region- based principles is described in §2. After outlining the algorithm in §2.1, the dierent types of ow structure interactions which need to be captured by the algorithm are dened in §2.2. §2.3 presents the structure attributes used to abstract ow features in the methodology. The multi-stage search for structure corre- spondence candidates is described in §2.4. The constraining process of obtained correspondence candidates is addressed in §2.5. §2.6 describes how candidates for complete events are formed from the constrained structure-to-structure correspondences, while the constraining and optimization procedure for event can- didates is detailed in §2.7. A strategy to treat periodic domain boundaries in the context of ow feature tracking is proposed in §2.8. The design of the graph data structure storing the accepted ow feature inter- actions over time is explained in §2.9. §2.10 introduces a graph mining approach to extract common patterns of ow feature interactions from the tracking graph. The introduced algorithm is applied to track ow features extracted from scalar (i.e., single-component) elds governed by dierent physical mechanisms. §3 investigates the temporal evolution of passive scalar structures with well-dened initial geometry in canonical homogeneous isotropic turbulence and shock- turbulence interaction ow congurations. The motivation and numerical setup of this ow congurations is given in §3.1 and §3.2. The temporal (DHIT) and downstream (STI) distribution of detected events is presented in §3.3. Structure volumes associated with dierent event types are shown in §3.4. Split and merge events are investigated in geometrical feature space in §3.5. §3.6 presents the downstream evolution of surface-averaged physical quantities relevant for the mixing process. Correlations between the local sur- face geometry and those physical quantities are highlighted in §3.7. The evolution of primary structures is presented in §3.8.1 focusing on geometrical feature space trajectories, surface area evolution and relation to scalar dissipation and vorticity on the surface. In a second application, motivated in §4.1, the ow feature tracking methodology the evolution of vortical structures is studied in temporarily evolving compressible mixing layers, for which the problem formulation and numerical setup is described in §4.2. The mixing layer simulations are validated in §4.3. §4.4 addresses how vortices are extracted from the mixing layer and presents their general statistics. By utilizing the 8 graph mining technique introduced in §2.10, §4.5 analyses merge events of structures newly created in close proximity of existing vortices. A detection method for individual hairpin-like vortices is presented in §4.6. The geometrical signature of detected hairpin-like structures is investigated in §4.7. The temporal distribution, and hence the relation to the Reynolds number, of hairpin-like structure detections is discussed in §4.8, which also investigates the formation of individual hairpins. The analysis of the dynamics of vortical structures in mixing layer concludes with a discussion of the cross-stream and spanwise growth of hairpin-like structures in the pre-self-similar stage of the mixing layer. A conclusion re ects the thesis in §5. 9 Chapter 2 Flow feature tracking methodology 2.1 Introduction We propose a tracking algorithm that follows the idea of explicit, correspondence-based tracking. Hence, the extraction of structures from each given time instant of the ow eld and the generation of attributes describing those features are performed as a pre-step to the actual tracking. This extraction and character- ization step is conducted for each time frame independently and described in §2.3 in detail. The exclusion of temporal coherence in the feature extraction leads to the correspondence problem, which is the starting point of the tracking eort. The proposed tracking algorithm can be organized in three main blocks: 1) correspondence search & constraining; 2) event formation & constraining; and 3) dynamic graph building. While we seek to solve the correspondence problem in rst block of the algorithm, the event formation and, especially, graph building go beyond the solution of that correspondence problem. The three blocks are car- ried out in sequence for two consecutive frames and repeated in a time loop for which the sampling interval length does not have to be constant, as illustrated in gure 2.1(a). After the structures have been extracted and their attributes are generated, a frame-to-frame correspon- dence search is conducted in multiple stages utilizing dierent subsets of the structure attribute set. The search process involves a higher dimensional neighborhood in the attribute space and, optionally, a radius search in the spatial attribute subspace. Each search is conducted in the forward direction from source to target frame and also backwards from target to source as described in detail in §2.4. The outcome of the search is a list of potential structure-to-structure correspondences. The proposed search method has its key counterpart in a constraining process by which the found correspondences are tested for physical realizability on a structure-to-structure correspondence level (see §2.5). The algorithm allows for exibility in terms of the types and order of the correspondence constraints. In general, the correspondence constraining can be based on structure attributes, regions, and/or physical laws. However, all correspondence constraints here discussed are region-based, adding the hybrid character to our algorithm. Contrary to existing tracking ap- proaches, the region-based constraining presently introduced does not rely on spatial overlap between source 10 Structure extraction & characterization Correspondence search Correspondence constraining Event formation Event constraining Dynamic graph building Graph analysis & further post-processing Attributes Correspondence Candidates Time advancing Regions Accepted Correspondences Event Candidates Attributes Accepted Events Event history Graph Tracking (a) Source & target point cloud Projection of spatial attributes Attribute scaling Bidirectional NN search in attribute space & radius search in spatial subset Duplicate correspondence removal & Condence collection Creations & disappearances detection Found correspondences (b) Figure 2.1: (a) Overview of the tracking work ow identifying three main blocks: correspondence search & constraining, event formation & constraining, and dynamic graph building. (b) Overview of the correspon- dence search process with the bidirectional nearest-neighbor search in the full attribute space and the radius search in the spatial subspace as its central element. and target structures. The structure region information is used to constrain what has already been found and is not meant to be utilized in the search. As described in 2.2, only continuation events consist of a single structure-to-structure correspondence. Multi-correspondence events, such as split and merge, are formed based on the structure-to-structure corre- spondences that pass the constraining process, 2.6. The formation of potential events is adding information that can be used to further test for physical realizability. Again, in principle the event constraining is open to be performed on a attribute, region or physics basis. In contrast to the implemented correspondences constraints, all implemented event constraints are attribute-based, §2.7. Similar to the correspondence con- straints, the algorithm allows exibility regarding the type and order of the event constraints. Depending on the application, in a sequence of event constraints the same event constraint type may be applied to dierent event types with dierent tolerances. The event constraining process often leads to signicant event 11 simplications, emphasizing again that the proposed search and constraining process needs to be seen as a unity. The accepted events are used to dynamically extend a graph data structure at the end of each tracking step (§2.9). The graph is directed and acyclic indicating the progressing tracking time, its vertices represent structures at a certain time instant and its edges display the found structure to structure correspondences between successive time steps. The building of multi-correspondence events requires a decision making process in terms of the continuation, creation and termination of branches which is based on the overall condence that was assigned to correspondence in the search and structure properties, such as lifetime and volume. The entire graph represents the lifetime and interactions of all structures present in the ow eld over time. Each branch of the graph represents the evolution of one individual structure. Hence, the graph can be queried to further post process the individually tracked structures. Even though the main purpose of the graph is to store the found structure evolution and interactions over time, the event history information contained in the graph can also be utilized in the event constraining of subsequent tracking steps. Periodic boundaries are commonly found in numerical datasets from which structures are extracted (e.g., instantaneous volumetric snapshots / frames obtained from simulations of turbulent ows). An adequate handling of domain periodicity poses additional challenges to the development of the tracking algorithm, as it aects not only the structure extraction pre-step, but also both the correspondence and event components. To provide an overview of the tracking process, the treatment of periodic boundaries in the tracking is not shown in gure 2.1(a). This aspect of the algorithm is discussed in section 2.8. The modular character of the algorithm makes certain steps optional and depending on the application and purpose of the tracking. In addition, new types of realizability constraints which may be needed for new application can be added easily and without modication of the architecture. The presented algorithm requires user input in form of search parameter and constraint settings. To tune these settings for a specic application, the algorithm oers several forms of visualization which provide feedback about the accuracy of the detected feature interactions to the user. Visualizing the feature inter- action can be seen as a rst verication step that delivers useful insight into how the search and constrain parameter need to be adjusted to obtain reasonable tracking results for the specic application. 2.2 Event denitions in correspondence search To track an evolving structure over time the structure needs to be identied in a series of consecutive time instants (frames). The sampling interval (or tracking time step) between consecutive frames may not be constant, but is assumed to be suciently small to match structures in consecutive frames and capture interactions between them. In the matching process, structures in a frame at time t n are referred to as `source' structures, while structures in the consecutive frame t n+1 are called `target' structures. Samtaney 12 et al. (1994) introduced a classication of events that describe interactions of structures within a sampling interval as follows: 1. Continuation: A structure continues from one frame to the next without interacting with another structure. The source structure has only one forward correspondence. 2. Creation: A new structure appears in the target frame. The target structure has no backward corre- spondence. 3. Split: A source structure fractures into two or more target structures. The source structure has multiple forward correspondences. 4. Merge: Two or more source structures merge into one target structure. The target structure has multiple backward correspondences. 5. Disappearance: A source structure disappears. The source structure has no forward correspondence. While a continuation is a one-to-one correspondence, split and merge events may be seen as a group of one-to-one correspondences. A generalized event can be dened as a collection of one-to-one correspondences. The order of a split or merge is dened as the dierence between the number of target and source structures. The complexity of events is related to the sampling interval. For example, as the sample interval length is increased, a cascade of two rst-order splits would eventually be detected as a single second-order split. A reduced sampling frequency decreases the amount of data which needs to be stored for the tracking application. Since our tracking algorithm aims to be applicable to datasets of large sampling intervals, we add a generalized event type, that combines splits and merges of multiple order, to the list of events: 6. Compound event: Multiple source structures merge and split into multiple target structures. While compound events have been neglected by most tracking algorithms found in the literature, Guo et al. (2015) reported the observation of compound events during the tracking of magnetic ux vortices. A higher temporal resolution is needed to decompose compound events into regular events. The goal of the tracking is to generate a list of those events for each individual structure when the matching process is repeated over a sequence of frames. This way the temporal evolution of each structure and its interactions with others are recorded. In addition to these events from one frame to the next, events occurring over several frames as combinations of the previously dened frame-to-frame events may be dened. For example, a merge of a structure with another structure from which it was previously split may be seen as a `reconnection event' and a sequence of split events can be seen as a `cascade split event.' However, the occurrence of those multi-frame events highly depends on the sampling frequency. For example, a reconnection event (split and merge) may be seen as a continuation if the sampling frequency is lower. In contrast to the frame-to-frame events that need to be detected by the tracking algorithm, multi-frame events do not have to be taken into account in solving the correspondence problem and are subject to the graph analysis following the tracking. 13 2.3 Structure attributes for event detection and rejection In attribute-based correspondence search, correspondences between structures in successive time steps are based on the comparison of attributes of the structures. An attribute of a structure is a quantitative property, such as its size or orientation. A group of attributes is a ngerprint used to identify a structure, in other words, to distinguish a structure from other structures in the same time instant. In general, the ngerprint of a structure may change from one time instant to the next. Dierent event types aect the attributes of a structure from one time instant to the next in dierent ways. For example, in most cases, a split event may aect the size of the structure more signicantly than a continuation event. Since the unique description of a structure at each time instant is an essential part of the tracking, the choice of attributes to form the ngerprint of the structure is a key element. Since a good ngerprint separates structures, including the spatial location of a structure in the group of attributes is a natural choice. The location of a structure can be represented by a number of attributes. The full information about the spatial region occupied by a structure is contained in the spatial coordinates of all points forming the isosurface. In the present algorithm, this spatial information is reduced to three attributes that are reserved to represent the location of the structure. Depending on the application dierent choices for those three attributes may be made. Candidates to represent the spatial information of a structure in three attributes can be based on solely the geometry of the structure, such as the three coordinates of the center of the bounding box of the structure, or including other properties such as its centroid. In all applications of tracking algorithm presented in this work, we use the center of either the axis-aligned or oriented bounding box of a structure to dene its location. Note that, although attributes like the center of the bounding box are meant to describe the position of the structure, the location dened by those three attributes is not always within the actual boundaries of the structure. In addition to the spatial location, the shape of a structure can be used to distinguish the structure from other structures in the time instant. Similar to the three attributes representing the location of a structure, we are looking for a reduced set of geometrical attributes that classies the geometry of isosurfaces. Dierent geometry description techniques have been suggested in the literature. Geometrical descriptors based on Minkowski functional have been used by Leung et al. (2012). In this study, we apply the geometrical characterization method introduced by Bermejo-Moreno & Pullin (2008). From the three main steps of this methodology, extraction, characterization and classication, we do not make use of the latter one. Starting from a three dimensional scalar eld at a time instant the extraction of structures begins with a multi-scale decomposition of the given elds into component elds representing dierent scales of the original eld. This rst step is not applied in our study and we reduce the extraction method to an isocontouring of the original eld. Disconnected regions of the isosurfaces extracted from a volumetric scalar eld dene the individual structures to be tracked over time. Thus, the structures subject to the tracking are disconnected closed isosurfaces of the original eld. Since the geometrical characterization and hence the tracking is limited to closed surfaces, it is essential to reconnect and close isosurfaces that intersect periodic domain boundaries, if 14 those are present in the simulation. Structures intersecting non-periodic domain boundaries are not closed, however, may still be included in the geometrical analysis and tracking attempt by closing them either with simply using the boundary that they intersect or by closing them using their mirrored extension, Demanet & Ying (2007). The geometrical characterization of each individual isosurface starts with the computation of two dierential- geometry properties at each point of the surface. Following Koenderink & Van Doorn (1992) the absolute value of the shape indexS and the dimensionless curvednessC are locally determined based on the principal curvaturesf 1 ; 2 g as S = 2 arctan 1 + 2 1 2 C = 3 V A r 2 1 + 2 2 2 (2.1) , where the structure volume V and areaA are used for a non-dimensionalization of the curvedness. Higher values of the curvedness correspond to lower values of the radius of curvature and hence a stronger locally curved surface. The absolute shape index classies the local shape of the surface as hyperbolic S < 0:5, parabolic S = 0:5 or (concave or convex) elliptical 0:5 < S < 1. Extreme values of S are associated with a symmetrical saddle point S = 0 and an umbilical point S = 1, while the shape index is undened for a planar surface. Since the attribute representation of each structure is based on attribute of dierent types, the value range of each attribute needs to be considered and a scaling operation into a common range becomes necessary, see section 2.4.2. While the absolute shape index is bounded as 0 S 1, the dimensionless curvedness only has the planar surface as a lower bound, thus C2R + . However typical values of C range in [0; 3]. The area-based joint probability density function P (S;C) of the pointwise quantities S and C is obtained provide a non-local description of the geometry of the surface. Based on the rst and second moments of the joint probability density function P (S;C) a so-called feature center ^ S; ^ C is dened which lies closer to region of high density than the mean center (asymmetry of the joint pdf). In addition upper and lower distances d S u ;d S l ;d C u ;d C l of the joint pdf in each of its variable are obtained. The set of parameter which span the so called feature space of parameters is completed by the global compactness parameter, which is dened as ^ = 3 p 36 V 2=3 A (2.2) and has a range ^ 2 [0; 1]. The compactness parameter is also known as the sphericity since it measures how close the shape of an object is to a perfect sphere ( ^ = 1) (Wadell, 1935). As a result of the geometrical characterization, each isosurface is represented as a point in the 7-dimensional feature space of parameters f ^ S; ^ C; ^ ;d S u ;d S l ;d C u ;d C l g. The feature space of parameters is designed to geometrically classify surfaces of dierent shape, e.g. identify clusters in the high-dimensional space. However, combing the feature space of parameters with the three attributes describing the location of each surface would construct a 10-dimensional space in which the search for correspondences needs to be conducted. Searches in high dimensional space are subject to the \curse of dimensionality" problem as observed by Bellman (2015). Furthermore with k being the number of dimensions, in kd-tree data structures which are suitable for searches in space of low 15 to medium dimensionality the number of data points N should be N 2 k in order to make the search signicantly more ecient than an brute-force search. Thus, the number of dimensions of the space in which the correspondence search in conducted should be large enough. to distinguish structures but not bigger than that. Consequently we combine the location information of the structure not with the full feature space of geometrical parameters, but with the subsetf ^ S; ^ C; ^ g. Typical regions in this reduced parameter space are near the (f ^ S; ^ C; ^ g) = (1; 1; 1) point for blob-like structures and (1=2; 1; ^ ) for tube-like structures, which become more stretched with a lower value of ^ . Sheet-like structures are characterized by low values of ^ C and ^ . The proposed tracking method is not limited to a particular choice of feature attributes. Depending on the application, dierent attribute sets may be chosen. However, in the present work, the set of attributes describing each structure in any given time frame is formed by the three spatial attributesfx c ;y c ;z c g and the three geometrical attributesf ^ S; ^ C; ^ g. Hence each structure is represented as a point in a six dimensional attribute space and all structures in a frame form a cloud of points in this attribute space. As described in section §2.4.2, a scaling and weighting operation is applied to the group of attributes. Dierent events aect the structure attributes in dierent ways. While the geometrical shape of the structure is marginally evolving in a continuation, splits and merges can signicantly aect the geometrical signature. Thus, to detect all types of events dierent subsets of the group of attributes are used in the search of correspondences (see §2.4). In addition to the attributes that are used for the correspondence search, another set of attributes can be used to test and reject found correspondences and events based on physical realizability. To test for realizability, the algorithm follows region-based approaches to constrain correspondences and attribute-based approaches to constrain events. In ow congurations with a dominant mean ow direction, the direction of the integrated velocity of the structures may be used as an attribute to test directionality of the found event. Depending on the application a structure attribute may be identied which can be assumed to be conserved within the event within a user-dened tolerance. Examples for such attributes could be the overall volume or mass which may be assumed to conserved from all source to all target structures of an event (see §2.7.2). 2.4 Detection of correspondences After all structures from the source and target frame have been mapped into the attribute space, a multi- stage search for correspondences between these two clouds of points is conducted utilizing rst the full set of attributes and then its spatial subset only. Matching source points with target points by comparing the structure attributes translates into minimizing the distance between matched points in a space which is spanned by those attributes or a subset of them. The search for correspondences is conducted in multiple stages in which both the source and the target clouds participate as query and searched clouds. An overview is presented in gure 2.1(b). In general, the diculty of frame-to-frame correspondence detection and the 16 complexity of the resulting events highly depend on the sampling frequency of the dataset. The proposed combined method of searching and constraining aims to be applicable to datasets with low sampling fre- quency where corresponding structures do not necessarily spatially overlap in consecutive frames. However, a suciently high sampling rate is assumed. The limiting value for the sampling frequency is application de- pendent. In a dataset with relatively slow feature dynamics, in terms of the attributes spanning the attribute space, a lower sampling frequency may be sucient than in applications characterized by quick feature dy- namics. The minimal required sampling frequency also depends on the scale of the structures subject to the tracking in relation to the structure velocity. For most tracking methods (especially region-based) small, fast moving features result in the most restrictive sampling frequency. However, in the attribute space, structures are represented in terms of spatial location and geometric shape. In other words, structures that move fast relative to their size but are characterized by minor changes in shape from frame to frame may result in a comparable minimum sampling frequency than structures that move more slowly relative to their size but deform faster.Thus, measures to evaluate the sampling frequency based on location and shape are needed. 2.4.1 Projection of spatial attributes As a pre-step to the correspondence search, the velocity of the structures can be utilized to compute an estimated position used to improve the spatial subsetfx c ;y c ;z c g of the structure attributes. The goal of the projection is to further minimize the distance of matching points in the attribute space before the actual search in performed. Either the source structures can be projected forward or the target structures backward. The direction in which the spatial projection is conducted prior to each stage of the correspondence search is chosen based on the event type whose detection is the main objective of that particular stage, as detailed in §2.4.3-2.4.4. In contrast to Muelder & Ma (2009), who predict the feature region in the target frame as an extrapolation using the feature's positions in the previous frames, we compute the estimate x 0 = (x 0 c ;y 0 c ;z 0 c ) for the spatial subset of attributes in the previous or successive frame from the current subset x = (x c ;y c ;z c ) as a linear projection x 0 = x tv avg (plus / minus sign for forward / backward projection) based on the average structure velocity dened as v avg = 1 A Z A vdA (2.3) In general,the quality of the projection depends on the length of the sampling interval. A higher-order projection could be obtained by utilizing the velocity information, not only from the current source or target frame but also from previous or subsequent frames. However, when evaluating the relevance of the spatial projection in general, and especially of including higher order contributions, it needs to be considered that topology-changing events, such as split and merge, can change the structure's center location from source to target frame more signicantly than predicted by the velocity-based displacement, where no change in 17 topology is assumed. In addition, application-specic phenomena that accelerate or decelerate structures (such as shock waves in a compressible ow) further diminish the benets of the projection. To account for these eects the spatial projection may be zonally adapted, for example turned o across a shock wave. The idea of a spatial projection will be further utilized in the constraining process of found correspon- dences (discussed in §2.5), where instead of the structure's spatial attributes, the structure's bounding volumes or the structure itself will be projected. In contrast to the projection of the spatial attributes fx c ;y c ;z c g, this spatial projection of the structure or its bounding volume may benet from including a rotational contribution, in addition to the translational component. 2.4.2 Attribute scaling Since the spatial attributes describe the location of a structure in physical space, the domain size of the application limits their range. In calculating the Euclidean distance with coordinates of dierent scales the coordinate with larger range would dominate the distance. Hence, a rescaling (standardization) of the attributes is necessary. We rescale an attribute as ~ = ( 0 )!, where 0 is the attribute origin and ! is the weight of the attribute. Dierent weights can be assigned to the geometric attributes, depending what attributes are to be emphasized. For the tracking application presented in this work, a weight dened as the inverse domain length is used for the spatial attributes to rescale them into a range [0; 1] and balance the contribution from spatial and geometrical attributes to the distance calculation. 2.4.3 Nearest-neighbor search in the attribute space Regarding the inter-frame evolution of a structure that continues from a source to a target frame and has no interactions with other structures neither in the source nor the target frame, the following observations can be made: 1. The spatial displacement of the structure is limited. Thus, the corresponding target structure may not be the spatially closest to the source structure, but can be found somewhere in its spatial neighborhood. 2. A continuation event can be characterized by small topology changes compared to split and merge events. The geometrical deformation of the structure is limited. Therefore, for a continuation event, a corresponding target structure can be expected in the close neighborhood of the source structure in the attribute space. Hence, a nearest-neighbor (NN) search in the full six-dimensional attribute space is conducted to detect potential continuation events. The outcome of a NN search depends on the search direction, e.g. if the nearest neighbor in the target cloud is found for each source point or vice versa. Hence, we conduct the NN search in both directions, starting with source-to-target as the rst step of the rst stage of the correspondence search. Dierent data structures and algorithms have been proposed in the literature to conduct a NN search in higher dimensional space, as for example summarized by Yianilos (1993). We conduct the NN search in the attribute space based on the construction of a kd-tree data structure as introduced by Bentley (1975); Friedman et al. (1977). In general, a kd-tree is a space-partitioning data 18 structure that hierarchically organizes the data points based on their position in the k-dimensional space aiming to allow ecient access to the organized points by positionally querying. To organize the data, the space is recursively split into half-spaces by axis-aligned hyperplanes. Since half-spaces are constructed in the k-dimensional space, the kd-tree is a binary tree. Each node of the tree is not only associated with a half-space but also with one of thek dimensions that is used as the splitting direction to generate the node's children. Although each nodes denes the direction of the splitting hyperplane, the construction of the kd-tree is not unique, because the plane can in principle be placed anywhere in the half-space. A standard rule is to split before the median point along the longest dimension of the half-space. Hence, assuming the overall space has comparable depth in each dimension, one cycles through the dimensions when moving down the tree. Once the kd-tree is constructed, several ways exist to nd the nearest neighbor to a query point. We note that always choosing the \good side" (i.e. in direction of the query point) of the binary tree when traversing down the tree can have a high probability of failure, especially when many data points are located close to half-space boundaries. Detailed strategies can be found in Liu et al. (2005), for example. We refer to the cloud which is formed by all query points in the NN search as the query cloud and the searched cloud which is partitioned in form of a kd-tree as the match cloud. In the search for the nearest target neighbor to each source point in our attribute space, the source points serve as the query cloud and the target points form the match cloud. Opposite roles are taken when the NN search is conducted to nd the nearest source neighbor for each target point. At each stage of the correspondence search, we dene a condence value for each found correspondence in the particular stage of the search, allowing the assignment of dierent weights on specic stages of the search. At the end of the search a total condence is computed (see §2.4.6). The condence evaluation is continued in the constraining process. For the bidirectional NN search, we determine the condence value depending on the user preference either as a general value for the NN search c i =c =! NN c NN or as a correspondence-specic value c i =! NN " c NN + (1c NN ) 1 d i d + 3 # (2.4) where ! NN and c NN < 1 are the user-provided weight and minimum value for the NN condence. d i is the Euclidean distance between query and match points in the attribute space for thei-th found correspondence, and d and are the mean and standard deviation of the distances of all found correspondences in this stage. Note that, before weighting, the condence value is bounded between c NN and unity, where a unitary condence is associated with a vanishing distance in the attribute space and hence maximal condence in the correspondence. Depending on the relative location of point in the source and target clouds in the attribute space, a correspondence may be found in both search directions and a condence value will be assigned to the same correspondence in both directions. If the mean and the standard deviation, d and , are included 19 in the condence determination (equation 2.4), then the condence value will not be the same in both search directions, in general. In contrast to a range search, the NN is not restricted to a maximal distance. Events other than continuation lead to a more pronounced separation of source and target points in the attribute space because of the associated topology changes. However, since the NN search is not rejecting correspondences by comparing to a certain maximal allowed distance, correspondences are still going to be found for the source and target points involved in these non-continuation events. The validity of these correspondences depends signicantly on the number of source and target structures and their location in the attribute space. The NN search for each source point has a tendency to detect either potential continuation or merge events, while the NN search for each target point detects predominantly continuation and split candidates. To detect merging structures, a forward projection of the source spatial attributes re ects the converging character of the merge event, while for a split event, a backward target projection captures the diverging targets. Hence, prior to the NN search from source to target, a forward projection of the source spatial attributes is conducted. Likewise, prior to the NN search from target to source, a backward projection of the target spatial attributes is performed, as described in §2.4.1. However, the validity of split and merge events found in the NN search stages of the search is uncertain due to the unpredictable changes of the geometrical attributesf ^ S; ^ C; ^ g that occur in a merge or split event, however certain tendencies exist, see §3.5. This problematic can be further illustrated by considering the behavior of a NN search in terms of disappearance and creation events. Consider a (synthetic) frame pair where a structure continues from the source to target frame and another small source structure disappears, e.g. has no true forward correspondence. However, visible to the NN search are two structures in the source frame and one in the target frame. Hence, the unrestricted NN search for both source structures will result in matching with the same target structure, while the NN search for the target structure will match the target structure with only one of the two source structures. Combing all found correspondences of the two NN searches a false merge event is detected. Analogously, a false split event is detected when in reality only a single continuation and creation takes place. These simple thought experiments emphasize a need to test the validity all correspondences found in the NN search, especially the non-continuation events, and a need to have an alternative search method to nd candidates for split and merge events. However, applying the NN searches without restriction still has a justication and certain advantages. The hypothetical example with two source structures and one target structure could also be used to illustrate an event where the source which was disappearing in the previous example is instead a truly merging / incoming structure which is moving fast and is small (hence is causing only minor topology changes in the merge). The dominant correspondence of this merge is the correspondence found by the NN search in both directions. However, a false disappearance and continuation instead of a merge may have been detected if the NN search would have been conducted with a maximum allowed distance. The proposed bidirectional NN search needs to be seen as one part of a multi-stage search and in combination with the 20 constraining process. The purpose of the search is to nd all possible correspondence candidates, which will be tested for realizability in the constraining stage. 2.4.4 Radius search in the spatial subspace The geometry of a continuing structure evolves smoothly between the source and target frame. The structure is locally deformed, but the target structure is characterized by a global geometric signature f ^ S; ^ C; ^ g which is comparable to its counterpart in the source frame. Contrary to the continuation, events involving interactions with other structures, such as split and merge, cause signicant changes in the shape of the structure. Consequently, corresponding structures are separated in the full attribute space. The change in shape not only modies the geometrical attributes of the structure, but also relocates its center, leading to a larger spatial displacement, in terms of structure centers, compared to continuation events. However, topology-changing events aect the geometrical attribute subspace much more abruptly than the spatial subspace: Despite the spatial displacement, corresponding structures are still expected to be found in their spatial neighborhood. Therefore, a correspondence search in the spatial subspace is conducted to exclude the geometrical attributes from the search. Although excluding the geometry from the search is motivated by jumps in the geometric signature caused by split and merge events, the search in the spatial subspace is still going to detect potential continuation candidates, which are not only neighboring each other in the full attribute space, but also in the spatial subspace. In contrast to the correspondence search employing the full spatial and geometrical attribute set, the search in the spatial subset is conducted in form of a radius search rather than a nearest-neighbor search as illustrated in g. 2.2(b). A k-nearest-neighbor search in the spatial subset would determine the k-nearest neighbors of the query point without restriction to their distance. In contrast, a radius search around a query point nds all match points within a xed distance from the query point. Conducting a k-nearest-neighbor search is not feasible, because the number k of nearest neighbors to be found cannot be specied, since no assumption of the order of the split or merge event can be made at this stage. Instead, searching for all candidates involved in the event comes with the requirement of having to dene a search radius. For the NN search of both stages an attribute scaling and weighting as described in section 2.4.2 is necessary because of the dierent scales of spatial and geometrical attributes. To conduct the radius search in the spatial subset a scaling is not applied, since all spatial attributes are used with equal weights. An individual search radius is assigned to each query point based on the size of the structure. Dierent options to measure the structure size can be considered. The spatial extent of objects is often measured in form of their minimum bounding box. A minimum bounding box of three dimensional object is dened as the rectangular cuboid with minimal volume which encloses all points of the object. Choosing a rectangular box as a bounding rather than any other shape such as a bounding sphere is computationally convenient when a domain is described using a Cartesian coordinate system. Depending on the constraints that need to 21 y x y x C 1 C 2 C 3 C 4 Figure 2.2: (a) Generic two dimensional shape enclosed by its axis-aligned bounding box (red) and oriented bounding box (blue).(b) Illustration of a continuation, rst order merge and disappearance event of generic shapes. OBB of source structures in blue, OBB of target structures in red. Search radius based on OBB half-diagonal without tolerance. Events detected by the radius search are a true continuation and a false second order merge, which again emphasizes the need for a constraining process based physical realizability. be satised by the box dierent types of bounding boxes exists. Figure 2.2(a) shows an isosurface enclosed by its axis-aligned bounding box (AABB) and oriented bounding box (OBB). While the AABB is dened by the minimum and maximum position of the structure in each coordinate direction, the OBB can be characterized by a center point indicating the location, an orthonormal set of three basis vectors indicating the orientation and three scalars indicating the symmetric extent of the box in the direction of its principal axes. The box's three principal axes are referred to as major, medium and minor axis and the corresponding scalars describing the extent are referenced as half-width, half-height and half-depth. In case of the oriented bounding box no restriction apply to the orientation of the box and the orientation is chosen to minimize the box's volume. The axis-aligned bounding box also minimizes its volume, but under the constraint that its principal axes have to be aligned with the axes of the coordinate system. Hence, the axis-aligned bounding box is simply dened by the structures minimum ans maximum x;y;z coordinate. Once the axis-align or oriented bounding box has been found for each query structure in the radius search, it is used to dene the the structures search radius. Possible choices for the search radius are the half-diagonal of the box, its longest half-dimension or its rms-value of all three dimensions. The denition of the structure's search radius is completed by specifying a tolerance factor (> 1). Tracking application presented in this work apply a search radius based on the OBB half-diagonal with a tolerance factor of 1.3. An alternative to forming the search radius based on the structure's bounding box is to use the structure's parameter, = 3V=A, which has been previously used to non-dimensionalize the curvature of the structure, see section 2.3. 22 Similar to the nearest-neighbor search in the full attribute set, the radius search is conducted in both directions. In sequence of the multi-stage correspondence search, we apply the radius search where the source points in the spatial attribute subspace serve as query points and the target cloud is the searched match cloud after the NN search is conducted with the source points as query points. To organize the multi-stage search, we dene a rst stage where the source points are the query points and the target points are the match points. The rst step of this rst stage is the NN search in the full attribute space from source to target cloud, followed by the radius search in the spatial subspace. The second stage of the search consists of the same two steps, but conducting the NN search in opposite direction and using the target points as query points and the source points as match points for the radius search. To perform the radius searches eciently, we organize the respective three dimensional match clouds in a kd-tree data structure. While the query clouds always contain all points of the source or target clouds, the user can optionally lter out points in the match cloud that have been already matched in previous search steps aiming to minimize the time spent on the kd-tree construction and search. Since the search radius is associated with the query structure size, split events (or split parts of a compound event) are more likely detected when a source structure is used as a query, while merge events (or merge parts of a compound event) have a higher likelihood to be found when a target structure is used as query. Hence, prior to the radius search using the source points as query, the target backward projection of spatial attributes is conducted. Likewise, the source cloud forward projection is carried out prior to the radius search using target points as query. In spite of this spatial projection, the detection of split and merge events still benets from the non- directionality (isotropy) of the radius search. For example, a highly convoluted, elongated source structure could split in a way where the line connecting the center of the splitting source structure and the center of the splintered target structure is largely deviating from the direction of the projection. Similar to the NN search, a condence value is calculated for each found correspondence in the radius search. The condence in each correspondence is dened as c i =! R 1 d i r i (2.5) where! R is a specied weight,d i is the center distance between the matched structures and r i is the search radius of the query. Note that, before weighting, the condence value is bounded in [0; 1], where 1 indicates maximal possible condence in the found correspondence. The distance d i of the correspondence is the distance between the projected, rather than the original, source (target) point and the target (source) point. 23 2.4.5 Detection of unmatched structures Based on the assumption that the source and target frames contain structures, carrying out the bidirectional and unrestricted NN search guarantees that all structures have at least one corresponding structure. Conse- quently, no creation and disappearance events can be found when the NN searches are performed. However, structures which are newly created in the target frame or have their nal occurrence in the source frame should not be part of any correspondence. The falsely detected correspondences involving these structures do not pass the constraining process. If the bidirectional NN search is conducted, the constraining process is the only mechanism to nd creation and disappearance events. In contrast, if only the spatial radius search is conducted, structures may be present in the source or target frame that have not been matched with a corresponding structure. Thus, a search is performed for unmatched source and target structures. The found structures participate in a creation or disappearance event. However, the pool of creation and disappearance events is enhanced during the constraining process, see §2.5 and §2.7. For example, a detected continuation that does not satisfy one of the applied constraints generates a disappearance and a creation event. 2.4.6 Collecting correspondences form the multi-stage search At the end of the multi-stage correspondence search a list of detected correspondences is obtained. A correspondence identies the involved source and target structure as well as the condence assigned to the correspondence. In the multi-stage search, the same pair of source and target structures can be matched with each other multiple times. For example, a correspondence may have been established in both direction of the NN search or the radius search generates a correspondence that has previously been found by the NN search. Those correspondences share the information of the involved structures. However, they are not complete duplicates, because the condence values dier based on the distances in the dierent spaces and the specied weight for each search step. The group of correspondences with common source and target structures can be merged into a single correspondence, storing the duplicated source and target information and summarizing the search record of this group in the form of a cumulative condence. 2.5 Constraining of correspondences In contrast to event constraints (§2.7) that test an event in its entirety for physical realizability, corre- spondence constraints evaluate the plausibility of single correspondence candidates found in the search as described in S2.4. Correspondences deemed not physically realizable on a structure to structure basis are rejected as a pre-step to the correspondence-based formation of events, to avoid the formation of complex, non-physical events in the rst place. Although attributes of the involved features can be applied to con- strain correspondences, all constraints proposed in this section are region-based, indicating the hybrid nature of the tracking algorithm. The presented region-based correspondence constraints can be classied in two groups working either with the feature's bounding volume, such as axis-aligned or oriented bounding box, or 24 directly with the polygonal mesh of the structure. Hence, the dierent types of correspondence constraints are characterized by dierent degrees of accuracy and computational cost. In particular, the axis-aligned bounding box representation of a structure, and thus all constraints based on it, are associated with the low- est accuracy, while the constraints invoking the oriented bounding box and polygonal mesh of the structure result in intermediate and high accuracy, respectively. To minimize the overall cost of the correspondence constraining, a sequence of constraints ordered by as- cending accuracy and cost is applied to each correspondence. Assuming the tolerance setting to be consistent throughout the constraints, it can be concluded that a correspondence that is rejected by a constraint of low accuracy would also be rejected by all subsequent constraints of higher accuracy. Hence, the correspondence does not need to be tested in the subsequent constraints. Especially, the more computationally expensive surface proximity constraint (§2.5.2), that works with the polygonal meshes, is only applied to a reduced number of correspondences. To continue the condence evaluation for each correspondence found in the correspondence search, a condence value is computed in each constraint. A found correspondence is rejected, if the calculated condence falls below the specied threshold for any correspondence constraint in the sequence. However, a source will be marked as disappearing only if all its correspondences are rejected. Analogously, for the involved target structure, a creation event will be generated only if there is no other correspondence which passes the constraining and has the same target. 2.5.1 Bounding volume constraints A computationally inexpensive, region-based evaluation of the physical realizability of a structure to structure correspondence can be conducted by utilizing a bounding volume representation of the correspondence's source and target structure. The correspondence constraints proposed in this section use either the AABB or OBB of a structure as its bounding volume, and examine the spatial relation between source and target bounding volumes. Similar to the projection of spatial attributes in the correspondence search, all bounding volume con- straints optionally utilize the surface-averaged structure velocity to either project the bounding volume of the source forward or the bounding volume of the target backward, depending on which projection minimizes the distance between the bounding volumes (dened as the distance between the bounding volume centers). The axis-aligned bounding box volume constraint tests the physical realizability of a correspondence by inspecting the spatial (volumetric) overlap between the resized AABB of the source and target structure in- volved in this correspondence. If the resized AABBs do not overlap, the correspondence is deemed physically unrealizable. The AABB of the source and target structures are resized by a common specied tolerance factor preserving the AABBs aspect ratios, as illustrated in gure 2.3(a). 25 y x y x Figure 2.3: (a)Two dimensional illustration of the axis-aligned bounding box volume constraint. While the solid boxes are the AABBs, the dashed boxes indicate the resized AABB of each structure. The constraint evaluates correspondence as previously found in a radius search (even without tolerance) based on the overlap region between a resized source AABB (blue) and three resized target AABB (red). Here, the correspondence between the source and the tube-like target structure would be rejected by the constraint. (b) Two dimensional illustration of the axis-aligned bounding box distance constraint. Source AABB in blue, while target AABB drawn in red. Displacement between box centers shown in green. Box radius in direction of the displacement highlighted. The volume V o of the overlapping region between the two AABBs can be computed as V o = 3 Y i=1 max min (S max i ;T max i ) max S min i ;T min i ; 0 (2.6) whereS max i ;S min i ;T max i andT min i are the maximal and minimal coordinates of the source and target AABB in i-th dimension. If the resized source and target AABB do not overlap (V o = 0), the correspondence is rejected by the constraint. Overlapping resized AABBs allows the correspondence to pass the constraint and a condence value is assigned as c =! V o min(V S ;V T ) (2.7) where ! is a specied weight and V S and V T are the volumes of the resized source and target AABB. The computed condence is added to the value generated in the search for each correspondence. Alternatively to the volumetric overlap test of the resized source and target AABBs, a correspondence can evaluated by examining the spatial extent or radius of both boxes in the direction given by the displacement between the centers of the AABBs. A two dimensional illustration of the axis-aligned bounding box distance constraint is presented in gure 2.3(b). 26 The radii r s!t and r t!s of the source AABB in direction of the target AABB and vice versa can be computed from the axis-aligned extents of both boxes. If the Euclidean norm of the center displacement is denoted askdk, the correspondence passes the constraint ifkdk=(r s!t +r t!s ) < ", where " is a specied threshold. The assigned condence value of a passing correspondence is determined relative to the threshold " and dened as c =! 1 kdk " (r s!t +r t!s ) (2.8) which is, for a passing correspondence, bounded between [0; 1] before weighting with the specied factor !. Depending on the orientation of the structure, the OBB can enclose signicantly less volume compared to that of the structure's AABB. Analogously to the volumetric overlap test of source and target AABB, the oriented bounding box volume constraint conducts a volumetric overlap test in terms of the resized OBBs of the structures involved in the correspondence. A brute-force algorithm conducts 2 12 6 = 144 edge-face checks to determine if two OBBs overlap or not. In contrast, a fast overlap test for two OBBs can be conducted by applying Minkowski's separating axis theorem (SAT) as suggested in Gottschalk et al. (1996). In general, the SAT can be used to test if two convex sets are disjoint sets in n-d Euclidean space. The SAT examines the disjointedness of two convex objects by conducting an axial projection of both objects onto an axis on which two intervals are formed by the projections. If those intervals do not overlap, then the axis is a separating axis and a separating plane orthogonal to the separating axis exists that separates the two objects, which are consequently proven to be disjoint. However, if there is an overlap between the intervals, the objects may still be disjoint. Further testing using another axial projection onto a dierent axis is needed. For any two points of a (object oriented bounding) box, the box contains the line segment which joins the two points. Hence an OBB is a convex set and the SAT can be used to test the overlap status of two OBBs. In addition, an OBB is also a polytope, which simplies the selection of axes for the axial projection of the SAT. If two convex polytopes are disjoint and there exists a separating plane between them, then there also exists a separating plane which is either parallel to a face of either polytope or parallel to an edge from each polytope. Since the separating axis is orthogonal to the separating plane, the axis can be found either orthogonal to a face or edge of either polytope. Since each OBB has three face orientations and there are nine pairwise combinations of edges, there are 15 cases of potential separating axes. A two dimensional illustration of the SAT applied to two disjoint boxes is presented in gure 2.4(a). The 15 test cases dened by the SAT are sucient to determine if two OBBs overlap or not. The testing process may be terminated early, as soon as one of the tests reveals that a separating axis exists. Following a similar approach to allow for some tolerance in the constraining as in the axis-aligned bound- ing box volume constraint, the SAT can be applied to resized OBBs of the source and target structure. Alternatively to testing the resized OBBs for spatial overlap, a certain distance between the original OBBs can be allowed. The Gilbert-Johnson-Keerthi (GJK) algorithm, Gilbert et al. (1988), iteratively computes the distance between two convex shapes and is therefore applicable to determine the distance between the 27 y x Separating axis Separating line Figure 2.4: (a)Illustration of the separating axis theorem in two dimensions applied to rectangular shapes. (b) Illustration of the minimal proximity between two generic shapes with overlapping oriented bounding boxes. The target (red) structure is enclosed in the concave region of the source (blue) structure. In spite of overlapping OBBs, the structures themselves do not intersect each other, but show a signicant distance (green) between each other compared to their size. source and target OBBs in the OBB distance constraint. To solve the distance query between two convex shapes A and B, the GJK elegantly recasts the problem by considering their Minkowski dierence A B =fa a ab b bka a a2A;b b b2Bg (2.9) The motivation to considerA B rather than the physical space becomes clear, when rephrasing the question about the spatial overlap of A and B to a question about common points of A and B. If A and B intersect, the two objects have at least one common point. Consequently, their Minkowski dierence A B contains the origin. In general, the distance between two convex shapes is the minimum distance between the origin and the Minkowski dierence of the objects. Hence, the GJK iteratively searches for the point on the Minkowski dierence closest to the origin, which is the point of minimum norm of the Minkowski dierence. Instead of computing A B explicitly, its point of minimum norm is approximated by iteratively building a series of simplices, which move closer to the origin in each iteration. A simplex is a convex hull of a set of anely independent points inR n . As the simplex moves closer to the origin, its point of minimum norm becomes a better approximation of the point of minimum norm of the Minkowski dierence. To update the simplex in each iteration, points of the simplex which do not contribute to the closest point to the origin computation, are removed from the simplex, while a vertex from the Minkowski dierence which is closer to the origin is added to the simplex. This vertex is selected by evaluating a support function on the Minkowski dierence in the direction of the current minimal distance vector approximation. Hence, the vertex of A B is selected which is farthest in the direction of the current distance vector. Further implementation details 28 of the GJK can be found in Bergen (1999) and Montanari et al. (2017). In collision detection application the GJK is typically followed by an application of the Expanding Polytope Algorithm, Van Den Bergen (2001), to obtain the penetration depth of intersecting convex shapes. However, in the presented correspondence constraining, a correspondence passes the OBB distance constraint if the source and target OBBs intersect, and the penetration depth is not further considered. If the OBBs are not intersecting, their distance needs to be evaluated to decide whether to accept or reject this correspondence. The computed OBB distance is compared to a constraint measure characterizing the size of the source and target structure. Depending on the application dierent choices of constraint measure r can be made, such as the maximum of the half diagonal of source and target OBBs or the maximum of = 3V=A of the source and target structures, which is the constraint measure used in the tracking application presented in this work. A condence value is assigned to passing correspondences as c =! 1 kd min k " max (r s ;r t ) (2.10) where d min is the minimum distance between the OOBs and ! and " are specied condence weight and tolerance values, respectively. The constraint measure of the source and target structure is indicated as r s and r t . 2.5.2 Surface proximity constraints Although their exible orientation makes OBBs more accurate representations of the ow feature's spatial extent than AABBs, both bounding volume types, being convex objects, are inaccurate representations if the structure is highly convoluted and has signicant concave regions. Figure 2.4 (b) illustrates a conguration where the OBBs of source and target structure overlap, while the actual structures enclose, rather than intersect, each other. Therefore, the narrow phase of the correspondence constraining solves distance queries between the ow features' actual isosurfaces rather than their bounding volumes. Distance queries for concave objects typically require increased computational eort compared to convex shapes. As an optional pre-step to the distance computation the isosurfaces of a correspondence are spatially projected according to the projection in the bounding volume constraining of the same correspondence. The distance between the surfaces is determined by means of a point-to-point distance calculation where the surfaces are interpreted as point clouds. Therefore, the distance query translates into a three dimensional nearest-neighbor (NN) search that determines the nearest neighbor in the points of the target structure for each point in the source structure. The minimum distance of the surfaces is obtained as the minimum of theses point-to-point distances. Similar to the nearest-neighbor search in the correspondence search, a kd-tree data structure is used to conduct the nearest-neighbor searches. Here, the NN search is conducted in lower dimensionality compared to the NN search used to nd correspondences, and a kd-tree specialized for three dimension can be applied (Blanco & Rai, 2014). Furthermore, the surface distance computation 29 can be terminated early, if the current point-to-point distance is smaller than a threshold, based on the grid size, below which the structure are considered to be intersecting. Correspondences of intersecting source and target surfaces are accepted and pass the constraining process based on local surface proximity. Analogously to the OBB distance constraint the distance between non-intersecting structures is compared to a constraint measure, taken as the maximum of source and target structure, to decide if the correspondence should be accepted or rejected. A condence value is assigned to a passing correspondence as in equation 2.10. 2.6 Detection of events At the end of the correspondence constraining process described in §2.5, a unique list of accepted one-to-one correspondences between source and target structures is obtained. By inspection of this list of corresponding structures a list of events as dened in §2.2 can be formed. An event groups of correspondences that have common source or target structures, as a structure in the source or target frame can correspond to zero, one or multiple structures in its complementary frame. Disappearing and newly created structures are detected among all source and target structures as those which do not have any forward or backward correspondence, respectively. A continuation event is identied as a correspondence in which both the query and the match indices are unique and neither the query nor the match index is part of any other correspondence. A group of correspondences with a shared source structure but distinct targets forms a pure split event, in which the duplicated source is the splitting structure. Analogously, a pure merge event of multiple order is build from multiple correspondences, which share the same match index. In the most general case, compound event is generated form a group of correspondences that not only has duplicates in the involved query indices, but also has common match indices. This event formation process generates a list of events describing the multi-structure interactions between source and target frames. Each structure from the source and target frames participates in one single event and, hence, appears only once in the obtained list of events. 2.7 Constraining and optimization of events In contrast to the constraining of isolated structure-to-structure correspondences presented in §2.5, the event constraining process evaluates the physical realizability of single correspondences in the context of the entire event to which they belong and tests the whole event for plausibility. Additionally, the event constraining evaluates if an event can be simplied into subevents, which consist of subsets of the correspondences of the whole event, under satisfaction of the constraint, i.e. each subevent extracted from the original event satises the event constraint individually. 30 2.7.1 Event-driven correspondence constraining The correspondence constraining described in §2.5 assesses the realizability of structure-to-structure corre- spondences based on information contained only in this particular correspondence itself. For example, at the stage of correspondence constraining the type of event an investigated correspondence potentially is involved in is undetermined. Contrary, the event-driven correspondence constraining also tests single structure-to- structure correspondences, but utilizes information from the entire event that the correspondence is part of. In the tracking framework a number of event-driven correspondence constraints are implemented, while their applicability depends on the ow conguration and also on the type of feature of the ow to be tracked. For example, in the feature tracking of passive scalar structures in compressible turbulence presented in §3, small structures discretized as triangulations with a number of points may disappear below specied minimum due to diusion, which should be accounted as a disappearance event. If a small structure disappears in close proximity to another continuing structure, the presented correspondence search could detect a false merge event rather than a continuation and a disappearance. To correct those false merge events or merge- parts of compound events a simple constraint based on the number of points of a structure is introduced. If the number of points of discretized source structures involved in a merge event (or the merge part of a compound event) is below the specied minimum number of points for a structure to be considered in the tracking increased by a user-provided oset, the constraint removes this source structure from the merge event, while the remaining part of the event is kept. A disappearance event is generated for the removed source. If the original event was a pure merge event, removing the source either reduces the order of the merge event or transforms the merge into a continuation event. Analogously, removing a source from a compound event reduces the complexity of the compound or transforms it into a pure split event. Since the number-of-points constraint is only applied to merging sources, the event type needs to be known prior to its application and the constraint cannot be used at the correspondence constraining stage. For ows characterized by a dominant mean ow direction, a constraint evaluating the directionality of an event can be useful to detect falsely formed events. For events other than creation and disappearance, measuring the plausibility of the directionality can be achieved by computing the angle between the surface- averaged velocity of the source structure and the distance vector between the structure centers (spatial attributes) of source and target for each correspondence in the event. The surface-averaged velocity indicates the estimated direction of motion of a source structure. The deviation from this estimation of the actual direction of the distance vector of the detected correspondence should be bounded by specied limits, if there is a dominant direction in the ow. Since topology changing events are expected to lead to larger deviations from the estimated direction of motion of a structure, dierent bounds are suitable for dierent event types. 31 2.7.2 Combinatorial extraction and conservation-based constraining of subevents In addition to evaluating single correspondences in the context of an event, a constraining process is applied to test the physical realizability of all correspondences of an event at once and examines the plausibility of simplied subevents extracted from the original event. This subevent constraining can be applied based on one or multiple attributes, which are assumed to be conserved in total during the (sub)event within a specied tolerance and have not been utilized in the correspondence search. For example, the structure volume, mass or another physical quantity can be chosen. Since the attribute is assumed to be conserved, the realizability of the (sub)event can be tested by measuring the similarity between the source and target sum of the attribute(s). Subevents are extract from an event by erasing one or multiple correspondences from the original event. Based on the similarity between the source and target sum of the conserved attribute a score can be assigned to each subevent, indicating the realizability of the subevent. The score s i;c of thei-th subevent with n sources and m targets based on a conserved attribute is dened as s i;conv = min i;s ; i;t max i;s ; i;t i;s = n X j=0 s i;j i;t = m X j=0 t i;j (2.11) where s i;j and t i;j are the values of the conserved attribute of the event's sources and targets, respectively. A (sub)event is considered to be realizable if its score is bigger than a specied threshold. Depending of the event type under consideration dierent types of subevents can be extracted from the original event. The only subevents that can be extracted from a continuation event are one disappearance and one creation. Hence, if a continuation event does not satisfy the conservative event constraint, there is no alternative other than decoupling the continuation into a disappearance and a creation. The extraction of subevents aims to avoid creation and disappearances as subevents, if possible. Original split and merge events that satisfy the constraint are therefore accepted without considering their potential subevents. A split event can fail the conservation constraint because of two reasons: (1) not all targets of the source are found in the search or meaningful correspondences have been rejected by previously applied constraints. Since at this point in the algorithm it is impossible to add correspondences to an event, there is no mechanism to recover from this rejection by the constraint. The split may be either kept or is completely decoupled into one disappearance and a creation event for each target of the original split. (2) the split is rejected because too many targets have been assigned to the source and passed the previous constraining process. Analogously, a merge can be rejected either because not enough sources are present in the event to fulll the constraint or too many source are connected to the merging target. A mechanism to recover from a split (or merge) for which too many targets (or sources) are present in the original event is to extract subevents by removing one or multiple targets (or sources) from the event and testing the realizability of the resulting subevents. For a split (merge) event with n targets (sources), the number of subevent to be tested is given by c = P n r=1 n!=(nr)!r!. 32 The subevent with the maximum score may be a split (merge) of reduced order or a continuation event. If the optimal score satises the constraint, this subevent is accepted and a creation (disappearance) event is generated for each of the remaining targets (sources). If the maximum score which can be achieved with the given pool of targets (sources) still does not satisfy a specied tolerance of the constraint, the optimized event may still be kept and either a creation (disappearance) is only generated for the remaining targets (sources) or the original event is completely decoupled. Under certain circumstances, the proposed correspondence search can generate compound events in which a pair of sources exist that has a pair of common targets, called double connections in the following. The correspondences between these sources and targets typically dier signicantly in their condence values. To simplify the event and especially to enable the extraction of subevents, the correspondence between these subgroup of source and targets within the compound event with the lowest condence may be removed if its condence is lower than a specied threshold. A compound event which has no double connections or from which double connections have been success- fully removed previously can be decomposed into simpler subevents by classifying its targets based on their backward connections. In a compound without double connections three types of targets can be identied: regular targets have a single backward correspondence; merge targets have multiple backward correspon- dences and are the only target to one of their sources; split targets have multiple backward correspondences and all of their sources have multiple targets. Due to its single backward correspondence, a compound event cannot be decomposed into subevents at a regular target. In contrast, both merge and split targets allow the decomposition of the compound event. However, to avoid disappearances, the way subevents are extracted diers between decompositions based on a split or merge target. Because their sources have multiple forward correspondences, split targets are assigned to all candidate subevents that are generated by decomposing the compound event at the split target. After the score of all candidate subevents is determined the split target is ultimately assigned to the subevent with the maximum score, while the other reduced subevents still need to satisfy the constraint after the split target has been removed. In contrast, since a merge target has one source which has only the merge target as a target, the merge target can only be assigned to the subevent in which this source participates, to avoid that this source becomes a disappearing source. Since the subevents extracted from compound events can themselves be compound events of reduced complexity, the compound decomposition needs to be applied in a recursive manner. Furthermore, the compound decomposition can be applied to the lowest possible level, i.e. aiming for subevents which are as simple as possible and still satisfy the constraint, or stopped at an intermediate level. In the latter case, the decomposition is stopped where further decomposing the subevents would lead to simpler events, which, on the other hand, have lower scores. 33 2.7.3 Sequence of event realizability constraints If an event is constrained by one of the described constraints, the resulting sub-events are generated imme- diately. Therefore the order matters in a series of these constraints, as the events a constraint is acting on depends on the position of the constraint in the series. A general guideline to order event constraints is to go from complexity to simplicity, i.e. constraints acting on complex compound events should be applied rst as they may generate simplied regular sub-events which need to be tested later on in the series in constraints acting on these simpler event types. 2.8 Treatment of periodic boundaries The applied geometrical characterization and classication of disconnected isosurfaces of the original eld is limited to closed surfaces. As the tracking depends on the characterization as a pre-step from which the feature attributes are obtained, the tracking itself is also restricted to closed surfaces. However, this limitation could be overcome if a dierent set of attributes would be used. Because the features are required to be closed surfaces, any interaction of a structure with a domain boundary that leads to open surfaces requires a strategy to close the structures and, consequently, enable their tracking. Structures intersecting with non-periodic domain boundaries may be included in the geometrical analysis and tracking attempt by closing them either by simply using the boundary that they intersect or by closing them using their mirrored extension. Contrary, structures intersecting periodic boundaries are continued at the opposite side of the domain and, hence, these structures can be closed by reconnecting their pieces across the periodic boundary. To reconnect a structure across a periodic domain boundary Bermejo-Moreno & Pullin (2008) apply a matching algorithm for the intersection curves between the boundary and the structure on both sides of the domain. After the structure pieces are reconnected in the periodic direction, the structure is present as one closed surface at one arbitrarily chosen side of the domain used in the reconnection process. Note that if all three dimension are periodic, the periodic reconnection needs to be applied in all three dimensions simultaneously to close a structure that intersects all three boundaries. Periodic reconnection enables the geometrical characterization of previously open surfaces intersecting periodic domain boundaries. However, additional steps are required to apply the tracking under periodic boundary conditions. In order to use the spatial attributes of a feature in the search for correspondences under periodic conditions, boundary crossing structures need to be detected and duplicated according to the periodicity of the domain. The spatial attributes of the duplicated structure are obtained from those of the original structure by applying a translation as specied by the domain size in the periodic direction, while the geometrical attributes are unaected by the duplication. Furthermore, to solve the correspondence problem between source and target frames under periodic conditions, boundary crossing structures must be detected in both frames prior to the correspondence search. To detect potentially boundary-crossing structures, their location in the subsequent (source) or previous (target) frame is estimated by a linear 34 Figure 2.5: Duplication of boundary intersecting or crossing structures under periodic boundary conditions for one frame. Straight arrows indicate the estimated location of the structure in the successive or previous frame due to forward source or backward target projection. Bent arrows show the estimated location of the duplicate structure at the opposite side of the domain. Bent solid arrows illustrate duplication due to periodic structure reconnection, while bent dashed arrows show duplication due to boundary crossing structures. projection, as described in §2.4.1. Structures that have not been reconnected prior to their geometrical characterization are duplicated if the radius search region or the bounding volume (whichever is bigger) of the projected structure intersects or crosses the boundary. Structures that have been reconnected across periodic boundaries are duplicated according to the original location of the structure fragments. However, the potential combination of boundary crossing and intersecting structure also needs to be taken into account (see blue structure in gure 2.5 for an illustrative example). The need to duplicate both sources and targets occurs because a source may correspond to a target which does not meet the criteria to be duplicated (source duplication needed) or vice versa. The search for correspondences does not distinguish between duplicated and original structures. Cor- respondences are individually established for each duplicate (gure 2.6a). However, each original structure and its duplicates are fundamentally the same structure, therefore, all correspondences found need to be col- lected and assigned to one structure after the correspondence search (gure 2.6b). For example, if both the source and the target of a correspondence are duplicated, the correspondence will also be duplicated, i.e. two correspondences are established at dierent sides of a periodic domain, with common source and target and are duplicated correspondence between the same structures. Hence, during the collection of correspondence, duplicated correspondences are removed. 35 S 1 T 0 3 S 0 1 T 3 S 2 T 2 S 0 2 S 3 T 1 S 0 3 T 6 S 4 S 5 T 4 S 0 5 T 0 4 T 5 (a) S 0 1 T 3 S 2 T 2 T 1 S 0 3 T 6 S 4 S 0 5 T 0 4 T 5 (b) Figure 2.6: Correspondence search and event reconnection under periodic boundary conditions. (a) individu- ally found correspondences (green) between duplicated and original source (blue) and target (red) structures. (b) resulting reconnected events. Primes indicate duplicated structures. Each structure is present in only one of the reconnected events in form of the original structure or one of its duplicates. 2.9 Graph representation of temporal evolution of ow features Graph data structures represent pairwise relationships between objects and are therefore ideal to organize the temporal evolution of ow features by mapping the accepted events solving the correspondence problem between pairs of frames onto a graph at the end of each tracking step. Objects are abstracted to form the set of vertices of the graph and each related pair of vertices is an edge of the graph. In the present tracking methodology of turbulent ow structures each vertex of the graph represents a ow structure at a certain instant in time, while each edge indicates the correspondences between source and target structures. Each vertex and edge carries properties associated to the represented structure or correspondence, respectively. In particular, the structure identication and the time instant in which the structure exists are essential vertex properties to interrogate the graph for post-processing purposes. An edge can, for example, store the condence of the correspondence associated with it. Here, edges connect from source to target vertex, i.e. have a direction associated with them and the resulting graph is directed. The source vertex of an edge always represents a structure at one time instant earlier than the structure represented by the target vertex of the edge. At the end of each tracking step, the accepted events are mapped onto the graph and extend it by one frame, resulting in an ordering of vertices into a temporal sequence. Therefore the graph is characterized by a topological ordering, which forbids directed cycles. The graph is a directed acyclic graph. After the list of accepted events is obtained for a tracking step, the graph needs to be updated according to this list of events. How dierent event types are incorporated in the graph building process is explained in the following. Intrinsic to the graph building process is an implicit vertex clustering which results into two- layer hierarchy of subgraphs, dened as tree and branch, of the overall graph data structure as illustrated in gure 2.7. By clustering the vertices and introducing the subgraph hierarchy, the graph re ects the life time 36 t 0 t 1 t 2 t 3 t 4 t 5 t 6 (a) (b) Figure 2.7: (a) Illustrative representation of the temporal evolution (from left to right) and interaction of structures of arbitrary shapes, where each color indicates an individual structure (undergoing merges, splits, creation and disappearances). (b) Representation of the structure evolution and interactions in (a) as a directed acyclic graph, with the same color scheme as the individual structures. Dominant edges in events indicated as solid arrows , while non-dominant edges are shown as dashed arrow . Node and structure colors identify the same structure in multiple frames. Node shapes indicate the event type from which the structure results: creation, continuation, split, merge, and compound events. Trees and branches indicated as solid and dashed wrapping. Branches classied as primary ( , , ), incoming ( , ), outgoing ( , , , ) and reconnecting ( ). Note that a branch can be outgoing from or incoming to multiple branches through compound events. Not shown are branches of connecting type which are outgoing from a branch and incoming to another branch. Branches that have connecting and re-connecting character are considered to be of mixed type, see also g. 2.9. 37 and relationships of the tracked features. Each tree graph is a subgraph of the overall graph data structure, while each branch is a subgraph of a tree. Each branch is a representation of the lifetime of one individual structure, i.e. the rst vertex of a branch is associated with the frame in which this particular structure rst appears, either resulting from a creation or split event. On the other end, the last vertex of a branch represents the nal presence of the structure either due to a following disappearance or a merge event. A tree is a composition of branches. All branches of a tree constitute the lifetime of structures that originated from the same parent structure, which is the root vertex of the tree. From this unique root vertex, all other vertices of the tree are reachable, which requires a minimum of connectivity between the branches of a tree. If structures originating from dierent root structures interact with each other , then their branches within the dierent trees will be connected. The only events that can connect trees are a merge or a compound. Due to the orientation of the edges, for a root vertex of one tree, there are unreachable vertices in the other tree and the trees are weakly connected. On the other hand, there may exist trees in the graph which are disjoint from all other trees. The structures of such an isolated tree only interact with structures that have the same root structure. Finally, the union of all weakly connected or disjoint trees describes the evolution of all structures present over all frames of the dataset. Each structure in the rst source frame of the tracking application is represented as a root vertex, which seeds a directed tree with currently one branch. After a root vertex has been generated for each structure in the rst source frame, the initial graph data structure consists of a number of trees equal to the number of structures in the frame. Each of these trees has one branch containing one vertex, which is the root vertex of the tree. Each root vertex is uniquely associated with a source of an event of the rst tracking step. As the tracking progresses over time, the source structures of all events can always be mapped to the set of added vertices in the previous tracking step. On the other hand, to extend the graph according to newly found events, the target structures of the events are added as new vertices to existing trees, but not necessarily to already existing branches. The only exceptions are target structures that are involved in creation events. By denition, a newly created target does not have an antecessor and is therefore represented as a new root vertex in the graph, similar to the structures in the rst source frame. A creation event always seeds a new tree with one active branch. Contrary, source structures that disappear do not have a successor and, consequently, a disappearance event terminates the branch to which the vertex representing the source structure has been added in the previous tracking step. Depending on how many active branches the tree has, the termination of the branch may also terminate the tree. The degree of a vertex can be used to distinguish dierent types of vertices. In a directed graph, the number of incoming edges of a vertex is its indegree, while the number of outgoing edges is the outdegree. Consequently, a vertex associated with a target of a creation event has an indegree of zero, while a disappearing source is a vertex of outdegree zero. A continuation event extends an existing branch of an existing tree. Since a continuation event is a single correspondence, the source vertex of the continuation has an outdegree of one, while its target vertex has an indegree of one. 38 (a) (b) Figure 2.8: Dierent vertex clustering strategies to build a rst order split in the graph. In (a) the branch of the splitting vertex is continued by one of its target vertices (continuation of the parent branch). Contrary, each of the targets in (b) starts a new branch, while the branch of the splitting vertex is terminated (out branching of the child branch). In both building variants all involved branches are part of the same existing tree. Solid arrows indicate a dominant edge, while dashed arrows represent a non-dominant edge. Because the vertices involved in creation, disappearance or continuation events have degrees of either zero or one, the vertex clustering into branches to build these events in the graph is unique. In contrast, the source vertex of split and the target vertex of a merge have degrees higher than one. Multiple vertices of degree higher than one exist in a compound event. As illustrated in gure 2.8, dierent strategies for the vertex clustering can be applied to build these event types in the graph. For example, the targets of these events could always start a new branch. This strategy applied to a split event does not need further decision making, as the newly created branches are part of the tree of the splitting vertex. However, the same approach used to build a merge event between trees already requires a process that decides to which of the merging trees the new branch should belong, as a branch cannot exist outside of a tree. To avoid this decision making process an alternative would be to not only branch out the target, but seed a new tree all together. An alternative vertex clustering approach aims to continue branches in the building of these event types. In case of a split event, the target vertex that continues the branch of the splitting source vertex would still be interpreted as the same structure after the split rather than a new structure which starts its lifetime after the split. Depending on the split order, multiple target candidates exist to continue the branch of the splitting vertex. Therefore, a strategy is needed to select the dominant correspondence (edge) of the split whose target is continuing the parent branch. The decision can be based on multiple criteria, however the result should consider that this selected target will be interpreted as the source structure after the split. Analogously, for a merge event, a dominant source needs to be selected, so that the merged target continues the branch of the dominant source. For a compound event, the decision-making process is applied to the sub-split and sub-merge events to select dominant sources and targets. 39 To allow a decision based on multiple criteria, a total score is calculated for each of the candidate correspondences as a weighted sum of scores of this correspondence for each criterion. Hence, the total score of the i-th candidate correspondence based on n criteria is dened as s total;i = n X j=0 ! j s ij ; n X j=0 ! j = 1; s ij 2 [0; 1] (2.12) where ! j is the specied weight for the j-th criterion and s ij is the score of the i-th candidate in the j-th criterion, which is bounded by [0; 1]. Since the dominant target of the (sub)split event will be seen as the successor of the splitting source, and the dominant source of a (sub)merge event will be interpreted as the predecessor of the merged target, a criterion to decide on the dominant edge(s) of an event can be based on the comparison between the volume of candidate and the reference volume. For a (sub)split, the reference is the splitting source, whereas for a (sub)merge, the reference is the merged target. The score of the i-th candidate in this volume-based criterion is determined as s i;volume = min (V i ;V ref ) max (V i ;V ref ) (2.13) Other criteria based on any physical quantity of the structures could be introduced in the same way as the volume-based criterion. In addition to the volume-based criterion, the condence value of each correspondence in the event can be used as an score for the candidate as well. This way, a dominant edge is linked to a correspondence of high condence. The condence includes a contribution from the geometrical signature due to the NN search. However a criteria explicitly based on shape could be added. In contrast to dominant target candidates for (sub)splits, the history of source candidates for (sub)merge events is available and can be utilized to form another criterion based on the lifetime of the source structure before the (sub)merge event. To assign a lifetime-based score bounded in [0; 1], the score of the i-th source candidate is dened as s i;lifetime = 1e 1 2 (titi;0) (2.14) , where t i;0 is the time of the frame in which the i-th source candidate is rst present as a individual structure. If two source candidates dier by t in their lifetime, their dierence in the lifetime score depends on the lifetime of the structures due to the exponential slope of the score. Two long-living source structures diering in lifetime by a small fraction of both structures lifetime t obtain a similar lifetime score and the decision will be dominate by another criterion. After the scores of each criterion have been determined for each candidate, the candidate that maximizes the total score s i;total is selected as the dominant source or target of the (sub)merge or (sub)split event. While the dominant target of a (sub)split continues the parent branch, each non-dominant target generates a new branch within the same tree. The outdegree of a splitting source vertex equals the number of targets. Contrary to a split, a merge can connect branches of the same or dierent trees. For a merge between 40 Figure 2.9: Classication of branches based on how they emerge and disappear. Start and end vertices of primary branches ( ) do not have an antecessor or successor, respectively. The start vertex of an incoming branch ( ) has no antecessor, but its end vertex has at least one successor, while the opposite scenario denes an outgoing branch ( ). The start vertices of connecting ( ) and re-connecting ( ) branches have at least one antecessor and their end vertices have at least one successor. The end vertex of a reconnecting branch is connected to a vertex of a branch from which the reconnecting branch originated, while the end vertex of a connecting branch is connected to a separate branch. indicates branch connecting edges. Figure adapted from Lozano-Dur an & Jim enez (2014). trees selecting the dominant source, not only species the continued branch but also answers the question to which tree the merged target is assigned since the continued branch is uniquely associated with a tree. The branches of the non-dominant sources of a merge event are terminated. The indegree of the merged target equals the number of sources of the merge. Finally, a hybrid vertex clustering approach can be used to outbranch targets of split or merge events if their total score or specied subscore is below a threshold. The hybrid approach especially useful for higher order events, where the targets may be interpreted as individual structures due to the multiple interactions of the event. For example, if a source generates multiple approximately equally sized targets in a higher order split, the targets may be outbranched. While post-processing the graph data structure resulting from the application of the tracking algorithm, the branches can be classied based on the connectivity of their start and end vertices as suggested by Lozano-Dur an & Jim enez (2014). The start vertex of a primary branch results from a creation event, i.e. it has no antecessor, whereas its end vertex is a disappearing structure, i.e. it has no successor. Branches that do not satisfy both criteria are considered secondary branches. More specically, incoming branches have start vertices without an antecessor, but their end vertex has a successor in another branch. Their counterpart, outgoing branches, are characterized by having a start vertex with antecessor, but having an end vertex without successor. In contrast, connecting and re-connecting branches are described by start and end vertices that are connected to other branches, as illustrated in gure 2.9. An end vertex of a re-connecting branch is connected to the same branch as its start vertex, whereas the end vertex of a connecting branch is connected to a dierent branch than its start vertex. Note that the connections between branches need not be regular split or merge events, but can also be more complex compound events. For example, an incoming branch could also end in a compound event rather than a pure merge. Therefore, the branches are classied 41 by the indegree of their start vertex and the outdegree of their end vertex, rather than the event types these vertices participate in. 2.10 Graph mining for pattern recognition in ow feature evolution The application of the ow feature tracking algorithm described in the previous sections results in a graph data structure representing the temporal evolution and interactions of all individually tracked structures. To extract common patterns in the structure evolution, data mining techniques can be applied to its graph representation. Mined patterns are subgraphs of the tracking graph which occur frequently and hence are of signicance for the ow feature evolution. While the input graph obtained from the tracking may have disjoint components, patterns are required to be connected subgraphs. In contrast to the frame-to-frame events detected by the correspondence search, patterns span multiple frames and therefore have a certain temporal extent. The actual occurrences of a pattern in the input graph are denoted as the instances of the pattern, while the pattern itself refers to the general denition of the subgraph. One of the established graph mining approaches is the Subdue graph-based relational learning system (Cook & Holder, 1993; Holder et al., 1994; Ketkar et al., 2005). The following describes how the Subdue algorithm is adapted in the present work to recognize common patterns in the temporal interactions of ow features. One of the common challenges of graph mining approaches is the immense number of potentially interesting subgraphs extracted from large graphs. In contrast to mining techniques focusing on frequent subgraph discovery, such as the gSpan algorithm (Yan & Han, 2002), Subdue, in its unsupervised form, searches for patterns which best compress the graph dataset based on the minimum description length (MDL) principle (Rissanen, 1998) and typically reports fewer, but ideally more relevant patterns (Cook & Holder, 2006). To further reduce the number of extracted patterns by allowing slight variations between instances of the same pattern, Subdue applies an inexact graph matching between a pattern and its instances (approximate graph isomorphism). During the ow feature evolution, continuation events typically occur much more frequently than all other event types. Consequently, most subgraphs extracted from a graph representing the feature evolution are also dominated by continuation events. These subgraphs have a high value in terms of the overall graph compression. Therefore the majority of extracted patterns by Subdue in its unsupervised form are continuation-dominated and less interesting in terms of interactions between structures. Hence, a number of measures are introduced to guide the Subdue pattern search towards patterns which allow a certain amount of continuations, but focus on structure interactions and events other than continuations. As a pre-step to the Subdue pattern search, the tracking graph is processed by labeling functions which assign additional attributes to the nodes and edges of the graph. In a labeled graph the pattern mining not only takes the connectivity (structure) of the extracted subgraphs into account, but also the attributes of the 42 Figure 2.10: Synthetic graph used to demonstrate the graph labeling process and serving as input for the Subdue graph miner. Boxes indicate the tree-branch subgraph hierarchy of the graph. Node shapes symbolize event types as detailed in g. 2.7. involved nodes and edges. To illustrate the graph labeling process, a synthetic tracking graph is considered in g. 2.10. The rst two node attributes to be used in the Subdue pattern mining are the in- and outdegree of the node in the tracking graph. The node degrees signicantly constrain the instance matching. For example, without these attributes, a rst-order split could be matched with a sub-split of a second-order split. The unequal outdegree of the splitting node prevents matching those two subgraphs to the same pattern. A third node attribute evaluates the forward and backward distance of the node to another node involved in an event other than continuation. The attribute, denoted as continuation attribute in the following, has a value of 0 if the node is closer than a user-specied limit to a source or target node of a non-continuation event, otherwise the value of this node attribute is 1. The purpose of this continuation node attribute is to guide the pattern search towards patterns including events other than continuations as explained below in more detail. 43 In addition to node labels, attributes are also assigned to the edges of the tracking graph to be accounted for during the pattern mining. For each edge, the outdegree of the source node of the edge and the indegree of the target node of the edge are used as two edge attributes. An additional attribute stores the information whether or not the edge is a dominant edge in the tracking graph. This edge attribute preserves the tree- branch subgraph hierarchy of the tracking graph in the pattern mining process. Assigning the described attributes to the nodes and edges of the synthetic tracking graph of g. 2.10 results in the labeled graph of g. 2.11. Contrary to the previously described attributes, which are based on the properties of the tracking graph itself, attributes used in the pattern mining process can also be based on information extrinsic to the graph. For instance, a node attribute can describe the shape of the structure that the node represents. Introducing such an attribute places a geometric constraint on the pattern mining process. Rather than using the geometric signature of a structure directly as node attribute, the attribute can be obtained from clustering all tracked structures in the geometric feature space. The ID of the cluster to which the structure belongs to is then assigned to the geometric node attribute. Attributes extrinsic to the graph are not utilized in the demonstration of the graph mining using the synthetic graph, but will be employed in mining the graph representing vortical structures in mixing layers in chapter §4. Node and edge attributes other than those introduced by the labeling functions, such as the structure ID or time attribute of a node, are ignored when the tracking graph is fed into the Subdue algorithm. The initial state of the pattern search is the set of all uniquely labeled nodes of the tracking graph. Subdue applies a variant of beam search to expand those pattern candidates: pattern candidates are extended in all possible ways by a single edge and a vertex. The obtained pattern candidates are evaluated and ordered in terms of graph compression (or domain-specic rules as described below). Only the top beam patterns are expanded in the next iteration of the pattern candidate extension, where the beam width is a user- specied input value. This pattern growth process terminates when the search space is exhausted and no further growth can be realized. To focus the pattern candidates on some sort of event activity other than solely continuations, a barrier is introduced into Subdue's pattern growth approach by utilizing the previously introduced continuation node attribute: Instead of extending the pattern candidate in all possible ways, only those extensions are allowed which add nodes to the pattern which have the same continuation attribute value as the nodes already in the pattern. 44 Figure 2.11: Synthetic input graph labeled with the node and edge attributes used in the graph mining. The node labels show the indegree (rst number), outdegree (second number) and continuation (third number) attribute of the nodes. The user-specied node distance used for determination of the continuation attribute is set to 2 in this demonstration. The optional structure geometry-based node attribute is not used in this example. Edge labels are formed from the outdegree attribute of the source vertex of the edge, the indegree of the target vertex of the edge, and the dominant attribute of the edge, separated by the slash symbol. The dominant attribute of an edge has a value of 1 for a dominant edge, and 0 otherwise. Instances of the best pattern of the rst iteration of the Subdue algorithm highlighted in purple. 45 Once the pattern search terminates, a graph compression procedure is invoked which substitutes all instances of the best pattern (based on the evaluation method, see below) in the tracking graph by single nodes, illustrated in g. 2.12 and 2.13. Those compressed nodes represent the pattern in the compressed graph. During the compression procedure, edges, incoming or outgoing from best pattern instances, are redirected to be incoming to or outgoing from the compressed node. The graph compression operation also aects the node and edge attributes of the graph. In contrast to regular nodes, the compressed nodes do not carry the two node degree attributes, but contain an attribute indicative of the the pattern this compressed node is an instance of. The continuation node attribute value is inherited from the nodes the instance compressed. As a consequence, compressed nodes can only be matched with other compressed instances of the same pattern. The source outdegree and target indegree attribute are removed from the edges incoming to or outgoing from compressed nodes, respectively, but the dominant edge attribute is kept to maintain the tree-branch hierarchy of the graph. After compressing the graph, the Subdue pattern search is carried out in another iteration on the compressed graph. This process is repeated, while a list of patterns ordered by value is reported after each iteration and always the best pattern of the current iteration is used for compression.(The Subdue algorithm is applied iteratively on the successively compressed graphs, using the best pattern for compression in each iteration.) Due to the compression the graph size is iteratively reduced and Subdue has a higher chance to nd larger pattern as pattern can involve compressed nodes. The pattern search of each iteration can be constrained by a minimum and maximum number of edges a pattern needs to consist of. A higher minimum number of edges can ensure a certain amount of complexity in the pattern, especially if only a few iterations of the Subdue algorithm are carried out. However, smaller patterns found in early iterations can serve as building blocks of complex pattern of later iterations. It's worth noting, that graph compression can disrupt the temporal ordering of the tracking graph and also generate parallel edges between nodes, which makes the graph a directed multi-graph as in the example shown in g. 2.13b. In contrast to the original Subdue version, the creation of self-loops (edges where the source and target is the same node) is prevented during the graph compression to avoid paths of innite length. 2.10.1 Criteria for pattern evaluation The pattern evaluation based on the graph compression using the MDL principle is replaced by a set of factors measuring the value of a pattern in terms of the application of the mining to ow feature evolution. First, a value of 0 is assigned to patterns consisting solely of nodes having a continuation attribute value of 1. Therefore those pattern candidates are never considered for expansion during the beam search. The coverage and compactness of a pattern are two measures suggested as domain-independent (Holder et al., 1994). The coverage is dened as the ratio of the number of edges in the pattern to the total number of edges in the graph (at the current state of compression) multiplied by the number of instances of the pattern candidate (minus 1). Despite of diering from the MDL principle, the coverage value of a pattern still values 46 (a) (b) Figure 2.12: Illustrative split-dominated graphs labeled with in- and outdegree node attributes and in- , outdegree and dominant edge attributes. The two input graphs in the top row of (a) are similar but not isomorphic. Compression of the subgraphs highlighted in red leads to two isomorphic graphs in the bottom row. Compression of a dierent subgraph selection of the same two input graphs also generates two isomorphic graphs in (b). Figure continued in 2.13 (a) . (a) (b) Figure 2.13: (a) Iterative graph compression leading to isomorphism.(b) Generation of parallel edges through partial compression of reconnecting branches. The example illustrates the preservation of the dominant edge attribute during compression. the graph compression potential of the pattern. The compactness value of a pattern is dened as the ratio between edges and vertices of the pattern. To favor patterns with structure interactions (rather than just continuation events), the product of all in- and outdegrees of the pattern's nodes is considered as another factor of the pattern evaluation. Here, the degrees of the nodes in the subgraph extracted from the graph at the current state of compression are used. The subgraph may contain compressed nodes which are not present in the original graph. For nodes having a degree of zero in the subgraph it is inspected if this degree value coincides with the degree value of the node in the original graph, if the node exists in the original graph. If this is the case, those nodes contribute with a factor of 2 to the product of node degrees to value creation and disappearance events, otherwise the node does not contribute to the node degree product. Note that a node may contribute with its outdegree 47 but not with its indegree and vice versa. The nal degree product of the subgraph is decreased by one before multiplication with the other factors of the pattern evaluation. To favor pattern candidate of high temporal length, the length of the longest path in the subgraph is considered. A path in a directed graph is a sequence of connected, directed edges. For the longest path computation edges are weighted in terms of event activity they are part of: For edges which are not incoming to and outgoing from compressed nodes the sum of the indegree and outdegree edge attribute decreased by one is used as an edge weight. This way, edges of continuations have a weight of 1 and edges of other event types get higher weights. For edges incoming to or outgoing from compressed nodes, the number of edges the node compressed is used as a weight to add attraction to those nodes. For nested compressed nodes the cumulative number of compressed edges is used as a weight. If the edge connects compressed nodes, the sum of the edges compressed by both nodes is used as a weight. The overall value of the pattern candidate is calculated as the product of its value in coverage and compression and its node degree-based and longest path-based values. 2.10.2 Application to a synthetic tracking graph The Subdue algorithm with the described modications and pattern evaluation rules is applied to mine patterns from the labeled synthetic tracking graph of g. 2.11. The minimum number of edges constraint for a pattern is set to 2 and the node distance limited to determine the value of the node continuation attribute is also set to 2 in this application. The synthetic graph compressed by the two instances of the best pattern of the rst iteration is drawn in g. 2.14. With described settings, the best patterns of the Subdue iterations after the 8-th iteration have only one instance in the graph and the mining is terminated at this stage. The compressed tracking graph after the 8-th iteration is drawn in g. 2.15, where the colored compressed nodes correspond to the best pattern of the particular iteration plotted in g.2.16. The node and edge attributes of the best patterns show their values based on the original graph( g. 2.11) rather than based on the connectivity of the pattern itself. When inspecting the best pattern of iterations 6 and 7, it is evident that the subgraphs only partially express the connectivity of the graph region from which they are extracted. For example, the merge node of the best pattern of iteration 7 has a indegree of 3 in the original graph, but in the extracted pattern the node has only 2 incoming edges. It is therefore useful to embed the instances of a pattern into their surroundings in the original graph to gain context information. Naturally, the embedding diers between the instances of the same pattern. The two instances of the best pattern of iteration 7 are embedded in g. 2.10.2. Nodes and edges drawn in solid indicate a pattern instance, while dashed nodes and edges show the connectivity of the instance to the surrounding input graph. One layer of predecessors and successors of the instance nodes is considered to represent the embedding. From the context information of the embedding it is clear that one of the source nodes of the second-order merge diers in in its attributes from the other node and the extracted pattern only partially covers the merge. 48 If a pattern instance contains compressed nodes, those are recursively decompressed before embedding the instance into the surroundings of the the tracking graph. Two scenarios can lead to pattern instances involving compressed nodes: First, a compressed node can be added in the pattern growth during the pattern search. In this case, all instances of the pattern will have a compressed node representing an instance of the same best pattern of a previous iteration. The second way a compressed node can be integrated into a pattern instance is in cases where a pattern instance overlaps with an instances of a best pattern. This best pattern can be the best pattern of a previous, the same as the current or even a subsequent iteration. Here, the compressed node is only added to the instance which overlaps with the other best pattern instance. Consequently, the variation between the pattern instances is increased. The second scenario is shown in g. 2.18, where only the rst instance of the pattern carries a compressed node. Depending on the overlap, it is also possible that dierent instances of the same pattern include compressed nodes representing dierent best patterns as illustrated in g. 2.19. 49 Figure 2.14: Compressed synthetic graph after the rst iteration. All instances of the best pattern of iteration 1, highlighted as purple subgraphs in g. 2.11, have been replaced by compressed nodes indicated in purple. Compressed nodes are label as PATTERN-i-j where i indicates the iteration number and j the instances number. Since always the best pattern of an iteration is used for compression, the pattern number is always 1 and hence not included in the label of the compressed node. 50 Figure 2.15: Compressed input graph after the 8-th Subdue iteration. The graph has been compressed by the best pattern of each iteration. Each best pattern is highlighted with a color corresponding to the pattern denitions in g. 2.16 51 Figure 2.16: Best pattern of each iteration of the Subdue algorithm applied to the synthetic input graph. Colors corresponding to compressed nodes in g, 2.15. Since node shapes can vary among the instances of a pattern, the double circle shape is used as a generic shape. Node and edge attributes show their values based on the original graph (g. 2.11). (a) (b) Figure 2.17: Embedding of the two instances of the best pattern of iteration 7. The nodes and edges drawn in solid indicate the actual pattern instances, while the dashed nodes and edges are the embedding of the instances in the original graph. 52 Figure 2.18: Embedding and decompression of the two instances of pattern 11 of iteration 4. The rst instances overlaps with an instance of the best pattern of iteration 4 in two nodes, highlighted in half-white half-green. Therefore the overlapping instance of the best pattern of iteration 4 is added to the rst instance of pattern 11 of iteration 4. The denition of the best pattern of iteration 4 is given in 2.16. Inserting overlapping instances of best patterns into instances of another pattern increases the variation among these instances. Figure 2.19: Embedding and decompression of the two instances of pattern 6 of iteration 5. The denitions of the involved best patterns are given in 2.16. Both instances contain compressed nodes of dierent pattern types, leading to quite a variation between the instances, when they are decompressed and embedded. 53 Chapter 3 Passive scalar structures in isotropic turbulence and shock-turbulence interaction 3.1 Introduction A comprehensive understanding of compressibility eects on turbulent mixing processes is of great relevance to advance air-breathing super- and hypersonic propulsion systems (Andreopoulos et al., 2000; Urzay, 2018). Shock wave-induced mixing enhancement has been proposed to accommodate the very short time scale requirements for nonpremixed combustion in these applications (Budzinski et al., 1992). Furthermore, com- pressibilty eects can originate from heat release in non-premixed and, especially, in premixed combustion (Bilger, 2004). Despite the identication of the scalar dissipation rate as a key parameter in non-premixed combustion several decades ago(Bilger, 1976; Pitsch & Steiner, 2000), the study of compressibilty eects on scalar mixing has been the focus only in more recent direct numerical simulations (DNS) of homogeneous isotropic turbulence (HIT) (Pan & Scannapieco, 2010, 2011; Ni, 2015, 2016; Danish et al., 2016; Panick- acheril John et al., 2020) and also shock-turbulence interaction (STI)(Tian et al., 2017; Boukharfane et al., 2018; Gao et al., 2020). These numerical studies conrm that the amplication of turbulence kinetic energy (TKE) and transverse vorticity variance in conjunction with decreased turbulence length scales (Lee et al., 1993, 1994; Ryu & Livescu, 2014) result in an increase of the scalar dissipation rate across a shock wave, as suggested by earlier experiments (Hesselink & Sturtevant, 1988; Agui et al., 2005). While the TKE amplication saturates(Lee et al., 1997) for Mach number, M, larger than 3, the am- plication of the scalar dissipation rate does not saturate at least up to M = 5, but is lowered for larger Taylor-Reynolds number, Re , and turbulence Mach number, M t (Gao et al., 2020). The shock-induced enhancement of scalar mixing is expressed by a reduced scalar Taylor-microscale and an intensied decay of scalar variance in the post-shock region (Boukharfane et al., 2018). The interaction of spatial velocity and scalar gradients contributes to the scalar transport (Danish et al., 2016). The increase of scalar dissipation across the shock is associated with the evolution of the alignment between the scalar gradient, r, the vorticity,!, and the strain-rate tensor eigenvectors, denoted as; and (Kerr, 1985; Ashurst et al., 1987; 54 Parashar et al., 2017), and ordered by their respective eigenvalues from most extensional to most compres- sive, . Recent statistical volumetric analyses of HIT (Ashurst et al., 1987; Danish et al., 2016) and STI (Gao et al., 2020; Boukharfane et al., 2018) have shown that the scalar gradient is most aligned with the most compressive eigendirection of the strain-rate tensor, followed by and (Batchelor, 1959; Vedula et al., 2001; Donzis et al., 2010), consistent with Batchelor's turbulence theory. The scalar gradient alignment with the eigenvectors of the solenoidal component of the strain-rate tensor shows very similar behavior to that of incompressible turbulence (John et al., 2019). Scalar gradient alignment with the most extensional eigenvector,, increases across the shock. For reactive scalars in premixed combustion, the alignment is also known to increase due to dilatational motion induced by heat release in the reaction zone (Swaminathan & Grout, 2006; Kim & Pitsch, 2007). In contrast, the alignment with the intermediate and most compres- sive eigenvectors decreases due to enhanced scalar dissipation immediately after the shock. A weakening of alignment between the scalar gradient and the most compressive strain-rate eigenvector in ow regions of high dilatation was recently reported (Panickacheril John et al., 2020). A strong anti-correlation found between the scalar gradient and the vorticity (Kerr, 1985) has been interpreted structurally as the passive scalar being wound around vortex tubes. The decrease of scalar gradient alignment with the intermediate eigenvector and the vorticity ! is thus expected since the alignment between and ! is strengthened by the shock-interaction (Livescu & Ryu, 2016; Boukharfane et al., 2018). The alignment analysis has been utilized in the development of SGS models for the ux of a passive scalar (Pullin, 2000). To complement statistical studies of the scalar mixing, the research objective this chapter is to investigate the dynamics of this process in decaying HIT and its enhancement by shock-interaction in STI from a structural point of view. 3.2 Numerical setup and initialization of passive scalar structures The presented tracking algorithm is applied to study compressible turbulent mixing of passive scalars from a structural rather than volumetric view. To elucidate the eect of shock waves on the geometry of the scalar structures and mixing enhancement, we compare temporally decaying homogeneous isotropic turbulence (DHIT), and the canonical statistically stationary interaction of a nominally planar shock wave with spa- tially decaying isotropic turbulence (STI) at two dierent Mach numbers. Shock-capturing direct numerical simulations of these two ow conguration are performed based on a numerical method described in detail in Larsson & Lele (2009) and Larsson et al. (2013). To include passive scalars in these ow simulations the compressible Navier-Stokes equations for a calorically perfect gas are complemented with additional advection-diusion equations @ t ( m ) +@ j ( m u j D m @ j m ) = 0 (3.1) for the passive scalar transport as described in detail in Gao et al. (2020). Here, is the density, m are the scalars, u j the j-th velocity component (j = 1; 2; 3), and D m the scalar diusivities. The most relevant 55 dimensionless parameters for the DHIT are the Taylor microscale Reynolds number, Re = u 0 rms = , where is the dynamic viscosity, = [u 0 2 u 0 2 =(@ 2 u 2 ) 2 ] (1=2) is the Taylor microscale and u 0 rms = (u 0 i u 0 i =3) 1=2 is the r.m.s value of the velocity uctuation, and the turbulent Mach number, M t = p 3u 0 rms =~ c and additionally the mean ow Mach number,M = ~ u 1;u =~ c, for the STI simulation. In these denitions, Reynolds averages of a ow quantity, f, are indicated as f, while ~ f denotes a Favre averaged ow quantity. The corresponding uctuations are f 0 =f f and f 00 =f ~ f. In the STI ow conguration, these dimensionless parameters are chosen as Re = 40, M t = 0:2 and M =f1:5; 3:0g , where Re and M t are taken just upstream of the shock location. The initial values of Re and M t in the corresponding DHIT conguration are chosen as (Re ;M t ) = (46; 0:23) for the M = 1:5 STI and (Re ;M t ) = (42; 0:21) for the M = 3:0 STI to account for the decay from the inlet to the shock location in the STI. The Schmidt number, dened as the ratio of the kinematic viscosity of the uid and the diusivity D m of them-th passive scalar, is chosen to be Sc = 1 in all presented DHIT and STI simulations, such as typically for supersonic mixing of gases. For the STI ow conguration, a resolution of 1280 384 384 was found to be grid converged, seeGao et al. (2020). The grid resolution in each coordinate direction of the DHIT conguration matches the transverse resolution the STI simulation. In a volumetric study of passive scalar mixing in (compressible) decaying isotropic turbulence the passive scalars are typically initialized to have dened (initial) spectra, such as E(k)/ k 4 e 2k 2 =k 2 0 with k 0 = 4. In contrast, to study the mixing from a structural point of view, we initialize the passive scalar elds with a well-dened initial structure such that their isosurfaces present a specic initial shape. Precursor DHIT simulations are performed and the initial passive scalar elds are added to a time instant in which the turbulence is fully developed (velocity derivative skewness approximately -0.5, temporally increasing dissipation-based length scale L " =K 3=2 =";K =u 00 i u 00 i : =2;" = ij S ij = ) and the Taylor microscale Reynolds number and the turbulent Mach number match their target values of Re = 46 (42)and M t = 0:23 (0:21), respectively. This time instant of fully developed turbulence complemented with initialized passive scalar structures serves as the initial eld of the DHIT simulations and also as the in ow database for the STI conguration, as described at the end of this section. To study geometric changes of canonical shapes in a turbulent background, for each simulation the passive scalar elds, denoted as , 2 , and 4 , are initialized as collections of uniformly spaced spheres of uniform radius equal to one, two and four times, respectively, the initial Taylor microscale, , of the underlying decaying turbulence. The scalar elds are initialized using a smooth radial prole that transitions between zero outside and unity inside the targeted spherical shapes (see gure 3.1). The smooth transition aims to avoid discontinuities in the scalar elds that can result in spurious numerical instabilities. The transition prole describes the scalar value (r) as a function of the distance r from the center of the ellipsoid as (r) = max + min 2 + max min 2 erf 2 4(rr in ) r out r in (3.2) 56 (a) (b) Figure 3.1: (a) Radial scalar prole at initialization. (b) Detected regions of scalar discontinuities (red) during the development from the initial condition in DHIT in a slice of the 4 eld. Symbol Radius Transition Spacing Number 5 3 12 3 2 2 10 6 6 3 4 4 20 12 3 3 Table 3.1: Initial conguration of passive scalar structures used in the DHIT and STI conguration. Initial- ized number of structures per DHIT domain block. Radial transition length from min to max . , where r in and r out are the inner and outer radius, i.e. (r =r in ) = max and (r =r our ) = min . The three initial congurations of passive scalar structures whose temporal evolutions are studied in this work are summarized in table 3.1. The congurations of each structure type in the initial DHIT frame and in the corresponding STI in ow database are identical. The same sampling interval length of t = 0 =12:7 is used for the tracking for all structures in the STI M = 1:5 and corresponding DHIT ow congurations, where 0 = =u rms is the initial large-eddy turnover time of the DHIT conguration. To account for the higher upstream velocity, the sampling frequency is doubled in the M = 3:0 STI and in the corresponding DHIT. Figure 3.2 shows the initialized 4 structures and their temporal evolution in the DHIT conguration at three dierent time instants. The application of the reconnection of surfaces intersecting periodic domain boundaries is demonstrated in 3.2 (c). Despite of initializing the passive scalar structures using the smooth radial prole 3.2, large gradients are detected in the passive scalar elds by applying a jump-based scalar discontinuity sensor, gure 3.1. In domain regions marked by the scalar discontinuity sensor, a fth-order weighted essentially non-oscillatory 57 (a) (b) (c) (d) (e) (f) (g) (h) (i) Figure 3.2: Evolution of structures educed from 4 (a-c), 2 (d-f) and (g-i) scalars in DHIT captured in dierent frames from left to right: t= 0 = 0:0 (a,d,g), t= 0 1:65 (b,e,h), t= 0 = 4:5 (c), t= 0 = 3:5 (f), t= 0 = 2:4 (i). Coloring indicates the identication of the individually tracked structures. Thus, color changes are associated with split and merge events. Periodic reconnection has been applied to structures crossing boundaries. For instance, in (c) the top part of the green structure is located outside the computational domain as it has been periodically reconnected to its lower part continuation inside the domain, leaving the void region on the opposite side (see also g. 2.5 for an illustration of the reconnecting algorithm). The rendering shows the periodically reconnected structure at one of its possible locations, all of which are considered in the correspondence search. 58 (WENO) scheme with HLL ux splitting, Ducros et al. (2000), is applied, while a sixth-order central dier- ence scheme is applied elsewhere. Discontinuities in the hydrodynamic ow eld near the shock (or shocklets, for suciently large M t ) as detected by a dilatation-vorticity based sensor (> 1:2! k ! k ) are analogously treated numerically. Since the radial length r out r in in which the scalar initially transitions from min to max is a x fraction of the nominal initial radius of the passive scalar structures, the region of non-zero initial scalar gradient thickens with increasing nominal radius (gure 3.3). As the structures deform after the initialization, stretching of the transitional region induces local increases of the scalar gradient magnitude, which are more pronounced the longer the initial transitional length is. The increased stretching of larger transitional regions enhances the mixing of the passive scalar eld (Villermaux, 2019). To generate the in ow database used for the STI simulation, three individual DHIT precursor simulations without passive scalars are performed and from each the time instant in which the Taylor-Reynolds number and turbulent Mach number reach their target values Re = 46(42) and M t = 0:23(0:21) to start the DHIT simulation including the passive scalar elds. From each of those three time instants two additional eld are obtained by applying rotations with respect to the y and z coordinate axes. The nine obtained boxes are blended together following a technique introduced in Larsson (2009). After the blending operation the passive scalar structures are initialized on the obtained cuboid grid in the same way as previously described for the DHIT ow conguration, gure 3.4.The same initial geometry is used for the DHIT and STI study to make comparison between the passive scalar mixing in both ow congurations, table 3.1. The nominally planar shock wave is located stationarily a distance 2 downstream of the inlet, which will be taken as the origin of the x axis (k 0 x 0). To avoid acoustic re ections from the out ow boundary into the subsonic region of the domain, a sponge layer of length is introduced before the outlet boundary (Larsson & Lele, 2009). Since damping terms are added to the Navier-Stokes equation in the sponge layer region, the ow eld in the sponge layer is considered to be not physically meaningful. Therefore the eective downstream length of the STI conguration relevant to the tracking is reduced to the region starting at the position downstream of the inlet where the structures rst become closed surfaces, which depends on the initial size of the structures, to the the starting position of the sponge layer. To obtain individual passive scalar structures from the underlying eld by iso-surfacing in principle any vale iso 2 [ min = 0:0; max = 1:0] can be used as an iso-value. Not only the shape of the passive scalar isosurfaces extracted at each frame of the simulation depend on the chosen contour value, but also the total number of educed structures in each step. For the STI application, gure 3.5a reports the total number of extracted structures in a certain frame in the statistically stationary state for contour values ranging from min to max for the considered passive scalar elds ; 2 and 4 . The number of extracted structures peaks around iso = 0:3 and decreases for larger iso-values. The downstream length of the region in which structures are present is signicantly aected by the choice of the iso-value. In general, a lower contour value leads to an increased lifetime of the scalar structures and hence an increased number of extracted structures. 59 (a) Figure 3.3: Scalar gradient of 4 (top row) and 2 (bottom row) visualized in plane cuts of the DHIT M = 1:5 domain: Initial transitional region 4 (top,left), Transitional region 4 att= 0 = 0:98 (top,right), Initial transitional region 2 (bottom,left), Transitional region 2 at t= 0 = 0:98 (bottom,right). 60 shock wave x y z blended region blended region 2 2 2 2 2 HIT in ow database in ow 2 4, STI domain sponge layer Figure 3.4: Problem setup of passive scalar mixing in STI and blended DHIT precursor simulations. Length of in ow database reduced in gure. Adapted from Gao et al. (2020). For iso < 0:3 very large connected regions are obtained and the number of extracted structures is reduced. Due to the structure deformation induced by the shock, the post-shock region is prone to the formation of widely connected isosurfaces as visualized in gure 3.5b. Large isosurfaces in the post-shock region are unfavorable in terms of tracking individual structures. Additionally, surfaces which span across the domain and connect to themselves via periodic boundaries have innite extent, remaining open and therefore lead to a failure of their geometrical characterization. The size of the isosurface of largest extent is found to be quite sensitive to the chosen iso-value iso and also depends on the choice of the initial spacing of the passive scalar structures. For the given spacing, an isovalue of iso = 0:6 is chosen to dene the passive scalar structure in all presented congurations, which decreases the initial size of the structures insignicantly compared to the nominal radius but prevents the formation of large isosurfaces in the post-shock region. From all isosurfaces extracted for iso = 0:6, only those having a polygonal mesh consisting of at least 50 points will be considered in the geometrical characterization (and, hence, in the tracking), to ensure a reliable calculation of the dierential-geometric properties. The tracking methodology is applied to the analysis of structures of passive scalar mixing with the settings presented in tables 3.2-3.4. In the NN search, the spatial attributes are scaled using by the inverse of the transverse domain length, such that they are balanced with the geometrical attributes (table 3.2). This scaling gives more emphasis to the mean shock-normal dominant ow direction in the STI conguration. In the radius search, a large tolerance " R is applied to the OBB half-diagonal to form a search radius that captures all possible correspondence candidates. The same settings are used for both applied correspondence constraints. The obtained distances are compared to a fraction" of the maximum of the values of source and target structures of the tested correspondence. For small structures, dened by having a number of surface 61 (a) (b) Figure 3.5: (a) Number isosurfaces N s of each scalar eld ; 2 ; 4 for dierent isovalues iso extracted from a frame in the statistically stationary state of the STI. (b) Formation of a isosurface of largest extent (green) for iso = 0:3 in the post-shock region. Shock visualized as a dilatation isosurface (gray). Nearest-neighbor search Radius search Correspondence constraints f x f y f z f ^ S f ^ C f ^ c NN ! NN " R ! R l ref " " rel ! 1=2 1=2 1=2 1 1 1 0:5 0:5 1:5 0:1 0:2 0:4 0:1 Table 3.2: Scaling factors, condence value, and weight for the NN search. Tolerance and weight utilized in the radius search. Structure reference length, tolerance, relaxed tolerance and weight used for both OBB distance and surface proximity correspondence constraints, as dened in §2.4.2-2.5. triangulation points smaller than 1000, a relaxed fraction" rel = 2" is used to provide a higher tolerance. Table 3.3 lists the parameter settings for the event constraints ordered according to their sequence of application. First, a conservation-based constraint is applied to compound events where each extracted sub-event needs to satisfy a constraint threshold " in terms of the estimated mass. A relaxed threshold " rel is used if the extracted sub-event is a continuation. Correspondences of the compound event with less condence thanc min may be removed to resolve double connections (see §2.7). After a surface point constraint is deployed to test merge events and potentially simplify remaining compound events, another conservation-based constraint is applied to compound events with a less restrictive setting. Subsequently, merge, split and continuation events are tested by a conservation-based constraint that skips splitting and continuing source structures with less than 1000 surface triangulation points. The weights used for the dominant target/source selection for the construction of split, merge or compound events in the graph data structure are listed in table 3.4. Here, the conservation-based criteria use again the estimated mass of the structures involved in the event to build. 62 Conservation Surface points Conservation Compound Merge, Sub-merge Compound Merge Split Continuation " " rel c min N pts " " rel c min " " " 0:5 0:25 0:75 250 0:4 0:15 0:8 0:6 0:3 0:2 Table 3.3: Settings of the event constraining. Type of the constraint and event type the constraint acts on. Event constraints are applied in sequence from left to right. (Sub-) Split (Sub-) Merge Conservation Condence Lifetime Conservation Condence 0:8 0:2 0:5 0:3 0:2 Table 3.4: Weights used for the dominant target/source selection to build (sub-)split/merge events in the graph data structure. 3.3 Downstream distribution of events The graph representation of the temporal evolution of all tracked 4 structures in the DHITM = 1:5, gure 3.6, indicates that the majority of detected events are continuations of structures which evolve from source to target frames without interacting with others. Second to continuing structures, the graph is dominated by split and, consequently, disappearance events. The reduced lifetime of the 2 structures compared to the 4 structures of larger initial size is apparent in the reduced length of the branches of the graph representing the 2 evolution in gure 3.7, as the sampling frequency is the same in both applications of the tracking. The number of outgoing branches relative to the number of primary branches is reduced in the graph representation of the 2 structure evolution. These ndings motivated by the appearance of the graphs are quantied in the following by investigating the temporal or equivalently spatial (downstream) distribution of dierent event types. Since the STI conguration is statistically stationary, the temporal evolution of individual passive scalar structures is investigated by means of their downstream evolution. Since the structure velocity in x is always positive the mapping from time to downstream location is monotonic, t n < t n+1 ) x(t n ) < x(t n+1 ). To compare the structure evolutions in the temporally decaying DHIT conguration to those of the STI, the time coordinate of the DHIT conguration is translated into a equivalent downstream coordinate as x hit = 8 > > < > > : u u t for t<x s =u u x s +u d (tx s =u u ) else (3.3) , where x s is the shock location and u u and u d are the mean downstream ow velocities upstream and downstream of the shock. The downstream distribution of dierent event types during structure evolution in the STI and DHIT conguration is presented in gure 3.8, where the downstream coordinate is normalized by thek 0 wavenumber and the initial Taylor micro-scale, , of the DHIT conguration. Here, the histogram 63 Figure 3.6: Hierarchical layout of the graph of the 4 scalar in DHIT (evolving from left to right). Colors highlight communities as detected by edge betweenness clustering, Girvan & Newman (2002). Figure 3.7: Hierarchical layout of the graph of the 2 scalar in DHIT (evolving from top to bottom). Colors highlight communities as detected by edge betweenness clustering, Girvan & Newman (2002). 64 (a) (b) Figure 3.8: Downstream distribution of dierent event types for the ; 2 and 4 passive scalar structures in DHIT and STI congurations for M = 1:5 (a) and M = 3:0 (b). Relative (normalized by bin size) his- togram counts are normalized by its bin size and the number of initial (DHIT) or incoming (STI) structures: N c : number of creations, N d : number of disappearances, N s : number of splits, N sub;s : number of sub-splits of compound events. Shock location at x= = 0. count of number of events per downstream interval is normalized by the interval length and the number of incoming (STI) or initial (DHIT) structures. Due to scalar dissipation, no creation of new structures are expected to occur from the underlying passive scalar eld, as the scalar value stays below the threshold used for iso-surfacing. The only expected source of creation events are incoming structures of the in ow database of the STI congurations. Figure 3.8 shows the expected sharp peak of creation events close to the inlet boundary. The peak location is shifted downstream as the structures increase in initial size from to 4 and become closed surfaces at a later downstream location. A small number of creation events is detected in the post-shock region, especially for the 4 structures. Those structures falsely detected as being newly created are small in volume and are small continuing structures or targets of split events in most cases. The distributions of the other event types are also shifted downstream with increasing initial size leading to practically no disappearance, split or merge events for 2 and 4 structures in the pre-shock region. The number and distribution of disappearance events are linked to the number and distribution of split events as the majority of (early) splits have one dominant structure similar in size of the parent structure, see §3.8.1, and a number of smaller targets which disappear in close proximity to the split event. The peak of disappearance events is closer to the shock locatio in STI than in the corresponding DHIT, which is a consequence of the shock-induced mixing enhancement. This dierence between STI and DHIT becomes larger for the higher MAch number M = 3:0. Note that gure 3.8 is not distinguishing between split or 65 (a) (b) Figure 3.9: Representative compound events of 2 structures in the post-shock region of the STI cong- uration. For each event (a)-(b) target structures are articially displaced for clarity. The diagram below each rendering represents the structures as nodes (with colors matching those of the associated rendered structures) and the correspondences as arrows, directed from source to target structures. Structure colors indicate branches as illustrated in g. 2.7. merge events of dierent order. While disappearance and split events start to occur approximately at the same downstream location, the downstream range of disappearance events is larger than for split events as structures (which have been continuing for a while) disappear without prior splitting. For all structure types ( ; 2 ; 4 ), the number of split and, consequently, disappearance events per initial (DHIT) or incoming (STI) structure is larger in the STI than in the DHIT conguration due to the added structure deformation induced by the shock. This dierence between STI and DHIT further increases with the shock Mach number and becomes larger with increasing initial size of the structures. While the peak location is approximately at the same downstream location for and 2 , the disappearances and splits peak further upstream (closer to the shock) in the STI than in the DHIT for the 4 structures, which indicates again the dependence of the structure deformation on the (initial) size of the structures. For both, DHIT and STI, the number of disappearance and splits per initial and incoming structure, respectively, increased with initial size, which can also be seen in the graphs, gure 3.6 and 3.7. The higher shock Mach number M = 3:0 induces more splits in closer proximity to the shock. The number of compound events is larger in the STI than in the DHIT, increasing with higher Mach number and peaks closely after the shock location. The sub-split events hidden in the compound events explain the reduced number of pure split events in the STI and, also, contribute to the number of disappearing structures closely behind the shock, especially for the 4 structures, as the number of compound events per 66 incoming structure is larger for increasing initial size of the structure, gure 3.8.Note that, despite the higher sampling frequency, which reduces the number of sub-splits in the DHIT used for the STI with M = 3:0 compared to the DHIT used for the STI with M = 1:5, the number of sub-splits increases for the higher Mach number in STI. In both DHIT and STI ow congurations and for all structure types the number of merge events per initial/incoming structure is small compared to the number of splits. Similar to sub-split events, there is a considerable number of sub-merge events hidden in compound events. While the number of pure split events is clearly larger than the number of compound events, the number of compound events is comparable or even larger than the number of pure merge events, especially in the STI conguration where more compound events, and hence more sub-merge events, occur. While there is typically one sub-merge per compound event, multiple sub-split events of potentially higher order occur per compound event, see gure 3.9 for some representative examples of compound events. Compound events are of importance to capture all cases of structure interactions in the scalar mixing process. Therefore it is essential for the tracking algorithm to account for compound events rather than neglecting them as in previous feature tracking methodologies. 3.4 Structure volumes of event types In this section we investigate the volumetric size of source and target structures of dierent event types. In addition to providing physical insight, this section also serves as a verication of the detected events. The volume histograms of disappearing and 4 structures in STI M = 1:5 and corresponding DHIT are shown in 3.10. As already seen in the downstream distribution of disappearance events presented in the previous section, the normalized number of disappearing structures per initial/incoming structures is larger for the 4 structures of larger initial size than the structures and for the same initial size more structures per incoming / initial structure disappear in the STI than in the DHIT, due to the increased number of splits per structure. As expected, the volume of disappearing structures of both structure types and in both ow conguration is small compared to the volume of a sphere of radius equal to the initial Taylor micro-scale in the DHIT simulation. The volume range of disappearing structures is approximately the same for both structure types. While the largest disappearing structure in the DHIT conguration is similar in volume to the largest disappearing structure in STI, the smallest disappearing structure in the DHIT is signicantly larger than the smallest disappearing structure in the STI conguration, because of the increased grid resolution of the STI in downstream direction. As both congurations apply the same minimum number of surface points criterion for a structure to be considered in the tracking, a smaller volume can be realized with this minimum number of points in the STI conguration due to the higher resolution. A similar grid resolution eect can be observed in the joint source and target volume histograms of continuing and 4 structures in STI M = 3:0 and DHIT, shown in gure 3.11. A deviation from volume-preserving continuation events is observed on small scales due to scalar dissipation and for larger 67 (a) (b) Figure 3.10: Volume histogram of disappearing (a) and 4 (b) structures in DHIT and STI M = 1:5 normalized by the number of initial and incoming structures. Volume V s of source structures relative to a sphere of radius of the initial Taylor-microscale . shock-crossing structures due to the compression of the structure and the scalar dissipation amplication by the shock. The joint histograms of splitting and 4 structures in STI and DHIT for both considered Mach numbers are formed by the volume of the source structure and the sum of the volume of the target structures of the split events, gure 3.12, where the grid resolution eect can be observed again. In accordance with the downstream distribution of split events, gure 3.8, the number of splitting structures per initial/incoming structure is larger for the initially bigger sphere 4 . Note that the splitting of a structure can occur after the structure participated in a merge event, which allows the ratio of the splitting volume to the reference volume, V s =(4=3) 3 , to be larger than the ratio of the volume of the initialized structures to this reference volume. While the largest splitting volumes are naturally found for the 4 structures, the smallest splitting structures have similar volumes for both scalar elds as the dissipation aects both elds on the same scales. Similar to the disappearing volumes, the largest splitting volume in the DHIT and STI are comparable, while the smallest splitting structure have less volume in the STI than in the DHIT due to the increased grid resolution. For both, DHIT and STI, the smallest splitting volume are signicantly larger than the smallest disappearing volumes. Hence, structures that have shrunk to a certain volume due to scalar dissipation or previous splits, do not split further and simply continue until they reach the end of their lifetime and disappear by dissipation. In agreement, structures of small volumes, near the end of their lifetime, are characterized by a relatively high compactness, ^ , which decreases the probability of a split (see §3.8.1). Figure 3.12 also indicates that the sum of the volume of the target structures of a split approximately conserves the volume of large source structures, while the sum of the target volumes is reduced compared to the source volume for small source volumes due to increased scalar dissipation on smaller scales. Figure 3.13 expresses an increased number of compound events in the STI compared to the DHIT for both structure types, and 4 . The number of occurring compound events per initial/incoming structure 68 (a) (b) Figure 3.11: Joint histogram of volumes V s of source and V t targets of continuation events in the evolution of (left) and 4 (right) structures in DHIT and STI forM = 3:0 normalized by the number of initial and incoming structures. Volumes of the source and target structures are normalized by the volume of sphere of radius of the initial Taylor-microscale. Green reference line indicates volume preserving splits. Marginal histograms included as projections. increases from to 4 structures for both DHIT and STI. Compared to pure split events, the smallest structures participating in compound events are larger as, by denition, a compound event has at least one sub-merge event, and merge events typically occur between structures of relatively large volume (as can be inferred from the volume histogram for merge events, not shown). Similar to split events, for 4 structures the portion of compound events for which the volume is approximately conserved is higher than for the smaller structures as the deviation occurs for smaller volumes. 3.5 Split and merge events in geometrical feature space The geometrical signaturef ^ S; ^ C; ^ g of source and target structures of split and merge events are investigated in this section, focusing on the eects of the initial structure size and the shock. The joint histograms ( ^ ; ^ C) and ( ^ S; ^ C) of the geometrical attributes of source and target structures of split events in DHIT (gure 3.14 (a)) and STI (gure 3.14 (b)), indicate that splitting source structures are characterized by lower values of the compactness, ^ , and the absolute shape index than their targets, while the dimensionless curvedness ^ C is comparable between sources and targets. Structures characterized by low compactness are more stretched and likely to contain a ligament portion, which increases the probability of these structures becoming the source of a split. Figure 3.15 suggests that once the highly stretched ligament section of a source breaks, the produced target structures exhibit a dome-like region in the region previously connected to the former ligament section, thus increasing their absolute shape index. The compactness ^ of splitting sources in DHIT (gure 3.14 (a)) decreases with larger initial size of the structures from to 4 , as the structures become more convoluted and stretched before they split. The higher convolution of the source structures for bigger initial size also 69 (a) (b) (c) (d) Figure 3.12: Joint histogram of volumes V s of source and V t of summed targets of split events in the evolution of (left) and 4 (right) structures in DHIT and STI for M = 1:5 (top) andM = 3:0 (bottom) normalized by the number of initial and incoming structures. Volumes of the source and target structures are normalized by the volume of sphere of radius of the initial Taylor-microscale. Green reference line indicates volume preserving splits. Marginal histograms included as projections. shows in a higher dimensionless curvedness ^ C and an increased amount of hyperbolic surface points, which lowers the absolute shape index ( ^ S < 0:5). In contrast, splitting structures tend to be dominated by elliptical points ( ^ S > 0:5). This dependence of the geometric signature of splitting sources on the initial size is linked to geometrical trajectories of the structures, which have a common starting point in the geometric feature space as all of them are spheres but develop very dierently depending on the initial structure size (see §3.8.1). Furthermore, for 2 and 4 structures, the ^ S- ^ C joint PDF of source structures is essentially enclosed by that of the target structures. As the geometrical trajectories indicate, split events early in the structure lifetime typically contain a dominant target structure which is similar in shape to the splitting source structure. Hence, this target is located in close proximity to the splitting source in the geometrical feature space. The compression of the passive scalar structures by the shock wave leads to lower values of compactness, ^ , and absolute shape index, ^ S, of splitting source structure in the STI, gure 3.14 (b), compared to the DHIT. This eect is most pronounced for 4 structures, which causes the increase of in the number of split events from the DHIT to the STI to be the highest for the 4 structures, as already observed in §3.3. The joint histograms ( ^ ; ^ C) and ( ^ S; ^ C) of the geometrical attributes of source and target structures of merge events in DHIT and STI presented in gures 3.16(a) and (b), respectively, show opposite relations 70 (a) (b) Figure 3.13: Joint histogram of volumes of summed source and summed targets of compound events in the evolution of (left) and 4 (right) structures in DHIT and STI M = 1:5 normalized by the number of initial and incoming structures. Volumes of the source and target structures are normalized by the volume of sphere of radius of the initial Taylor-microscale. Green reference line indicates volume preserving compound events. (a) (b) Figure 3.14: Joint histograms of geometrical attributes of source and target structures of split events in DHIT (a) and M = 1:5 STI (b). Top: , Center: 2 , Bottom: 4 . Contour line: , : 0.2, , : 0.4, , : 0.6, , : 0.8. 71 (a) (b) Figure 3.15: Absolute shape index mapped on 2 structures in DHIT during split events of order one (a) and order two (b). Source structures articially displaced with source on the top and target on the bottom. between source and target geometrical signature than previously observed for split events. Here, the merged target structure is characterized by a lower compactness. The dierence between merging sources and merged targets in terms of the absolute shape index is less pronounced than previously observed for split events. The dependencies of the geometric signature on the initial structure size previously observed for splitting source structures apply to merged target structures. 3.6 Downstream evolution surface-averaged physical quantities Ensemble statistics of the downstream evolution of surface-averaged physical quantities of the ow eld computed as q avg = 1 A Z A qdA (3.4) on individually tracked isosurfaces of and 4 passive scalar in STI and corresponding DHIT congura- tions are investigated in this section. The eect of the shock on these surface-averaged quantities is more pronounced for a higher Mach number M = 3:0 in the STI conguration. The downstream evolution of the ensemble mean and standard deviation of the magnitude of the scalar gradientjrj normalized by its initial mean value in the DHIT is presented in gure 3.17. In both ow congurations, the in uence of the size of the transitional length of the initialized structures (see gure 3.3) on the early development of the scalar gradient magnitude is apparent. While the scalar gradient magnitude of structures immediately decays after the initialization, stretching of the non-zero scalar gradient region of the 4 structures cause an initial increase of the scalar gradient magnitude. After the initial increase the magnitude of scalar gradient monotonically decays in the DHIT conguration as the scalar mixing progresses. In contrast, the scalar gradient magnitude is amplied while the structures pass through the shock wave in the STI congurations. This amplication is larger for a larger initial size of the structures which indicates the scalar dissipation enhancement is larger for those structures. After the scalar gradient magnitude has reached its maximum, 72 (a) (b) Figure 3.16: Joint histograms of geometrical attributes of source and target structures of merge events in DHIT (a) M = 1:5 STI (b). Top: , Center: 2 , Bottom: 4 . Contour line: , : 0.2, , : 0.4, , : 0.6, , : 0.8. its decay is more pronounced than in the corresponding DHIT as the scalar dissipation rate is amplied by the shock (Gao et al., 2020; Boukharfane et al., 2018). Contrary to the downstream evolution of the scalar gradient magnitude, the vorticity magnitude does not exhibit a transient development as the turbulence is already fully developed at the initialization of the structures. The downstream evolution of the ensemble statistics of the surface-averaged vorticity magnitude shown in gure 3.17 is normalized by the initial mean value in the DHIT. As the turbulence decays in the DHIT, the surface-averaged vorticity magnitude monotonically decreases over time. In contrast, the vorticity magnitude is amplied as the structures pass through the shock wave. Unlike the scalar gradient magnitude, the increase of the vorticity magnitude is similar for all initial structure sizes. Similar to the evolution of the scalar gradient magnitude, the vorticity magnitude decays more rapidly in the post-shock region than in the DHIT. The scalar gradient and the vorticity show no preferential alignment at initialization since the turbulence is homogeneous and isotropic. After initialization a transient develops in which the alignment decreases, gure 3.17.A strong anti-correlation between the scalar gradient and the vorticity was found by Kerr (1985), interpreted as the scalar being wound around vortex tubes. In the present STI conguration, the alignment is further decreased when the structures pass through the shock wave. The minimal alignment during transient becomes smaller with increasing initial structure size. Again, the downstream evolution in the STI is shifted 73 Figure 3.17: Ensemble mean of the downstream evolution of the scalar gradient magnitude (left), vorticity magnitude (center) and the alignment between the scalar gradient and vorticity averaged on (row 1 for M = 1:5 and 3 M = 3:0) and 4 (row 2 M = 1:5 and 4 M = 3:0) isosurfaces. Color lled region indicate one standard deviation. Structure location dened by its most downstream surface point. 74 against the DHIT evolution, especially in the initial rapid decay of the alignment. The alignment is increasing for all initial structure size in the post-shock region. While the alignment in the STI has a lower minimum than in the DHIT and stays below the alignment in the DHIT for structures, the minimum is bigger in the STI than in the DHIT for 4 and the alignment reaches similar values in both congurations further downstream in the post-shock region. The amplication of the passive scalar dissipation rate across the shock is associated with evolution of the alignment between the eigenvectors of the strain-rate tensor and the scalar gradient (Kerr, 1985; Ashurst et al., 1987; Danish et al., 2016; Parashar et al., 2017).Alignment analysis of these vectors has been utilized in the development of SGS models for the ux of a passive scalar (Pullin, 2000). Figure 3.18 show the downstream evolution of the ensemble mean and standard deviation of the surface-averaged alignments betweenr and the strain-rate eigenvectors, denoted as; and with their respective eigenvalues ordered from most extensional to most compressive, . Due to the ad-hoc initialization of the passive scalar structures and the isotropic, homogeneous back- ground turbuence, the upstream starting point in the DHIT corresponds to no preferential alignment (all cosines of those angles equal approximately 0:5). After the initialization, a transient development is observed for all alignments on all considered isosurfaces, which lasts longer for larger initial size of the structures. According to Batchelor's theory, the scalar gradient is most aligned with the most compressive eigendirec- tion of the strain-rate tensor Batchelor (1959); Vedula et al. (2001); Donzis et al. (2010). In agreement, the surface-averaged mean alignment with becomes the largest, followed by and. These results, obtained for ensembles of initially spherical structures from surface-averaged data, are consistent with previous volu- metric analyses of HIT (Ashurst et al., 1987; Danish et al., 2016) and STI (Boukharfane et al., 2018; Gao et al., 2020). In addition, the scalar gradient alignment with the eigenvectors of the solenoidal component of the strain-rate tensor shows very similar behavior to that of incompressible turbulence (John et al., 2019). After the transient, the alignment with becomes nearly constant, whereas alignment with () slightly decreases (increases) in a similar way for all considered initial structure sizes. The alignments in the STI exhibit a similar transient development than in the DHIT with a temporal (downstream) shift which becomes more signicant with increasing initial structure size. While the transient development aects the entire structures immediately after their initialization in the DHIT, only the part of an incoming structure in the STI conguration which is already inside of the domain experiences the transient development, while the part of the structure outside of the domain remains in the initial conditions of the in ow database. Once the incoming structures become closed surfaces, the transient development is more advanced on the further downstream region of the surface. Hence, the surface-averaged transient development is shifted compared to the DHIT structures and the shift is more signicant for larger incoming structures, as they become closed surfaces at a more downstream location. A change of slopes is observed as the structures traverse the shock, which is more pronounced for the small structures due to the shorter streamwise distance of interaction and for the higher Mach numer M = 3:0. In particular, the scalar gradient alignment with the 75 Figure 3.18: Ensemble mean of the downstream evolution of the scalar gradient alignment with the (left), (center) and (right) eigenvector of the strain-rate tensor averaged on (row 1 for M = 1:5 and 3 M = 3:0) and 4 (row 2 M = 1:5 and 4 M = 3:0) isosurfaces. Color lled region indicate one standard deviation. Structure location dened by its most downstream surface point. 76 most extensional eigenvector increases across the shock. For reactive scalars in premixed combustion, the alignment is also known to increase due to dilatational motion induced by heat release in the reaction zone Swaminathan & Grout (2006); Kim & Pitsch (2007). The alignment with the intermediate and most compressive eigenvector decreases due to enhanced scalar dissipation immediately after the shock. The decrease in alignment with the intermediate eigenvector and the vorticity ! is thus expected since the alignment between and ! is strengthened by the shock-interaction (Livescu & Ryu, 2016; Boukharfane et al., 2018). A weakening of alignment between the scalar gradient and the most compressive strain-rate eigenvector in ow regions of high dilatation was recently reported by (Panickacheril John et al., 2020). In the post-shock region, the alignments in the STI reach similar values than in the DHIT.Especially for the 2 and 4 structures the variance of the ensemble statistics increases over time in the DHIT and STI conguration. 3.7 Correlations between local surface geometry and physical quantities on the surface The local geometry of a surface can be characterized by the pointwise values of dimensionless curvedness C and the absolute shape indexS on the surface (see gure 3.19a,b). Additionally, pointwise physical quantities can be mapped on the surface, such as the alignment between the scalar gradient and the most compressive eigenvector of the strain-rate tensor (gure 3.19c). To eludicate correlations of local physical quantities on a passive scalar isosurface with its local geometry, the average alignment of scalar gradient and conditioned on S and C accounting for the area-based joint PDF is obtained (gure 3.19e). For the sample structure shown in gure 3.19, a clear correlation between at regions (low C) and high alignment ofr and is observed on the surface renderings (a,c) and the conditioned joint PDF (e). Figure 3.20 shows ensemble mean and standard deviation of the suraface distributions, conditioned on S and C, of alignments betweenr and strain-rate eigenvectors and the vorticity, and the scalar gradient magnitude obtained on 4 isosurfaces, for collections of isosurfaces centered around four dierent stream- wise intervals at any given time of the simulation. Positive correlations with C < 5:0 are observed for the alignment ofr (which is parallel to the normal to the isosurface at every point, n) with and, indicating that regions of higher curvature present better alignment between the scalar gradient and the most exten- sional and intermediate strain-rate eigenvectors. The opposite trend is observed for the alignment with the most compressive strain-rate eigenvector, , indicating that atter regions of the isosurfaces present larger alignment betweenr (and thus n) and . ForC < 5:0 these correlations are strengthened as the isosurfaces are transported downstream. For regions of very large curvedness (C > 5:0) the expected alignment with and slightly decreases, whereas the opposite trend is observed for the alignment with . Across the shock, the scalar gradient alignment with ( ) increases (decreases) for the entire range of the absolute shape 77 (a) (b) (c) (d) (e) Figure 3.19: Local values of the dimensionless curvedness (a), absolute shape index (b) and the alignment between the most compressive eigenvector of the strain-rate tensor and the scalar gradientr 4 (c) mapped on a isosurface of the 4 at t= 0 = 4:3. (d) Area-based joint and marginal probability density functions of S and C. (e) Average alignment of scalar gradient and conditioned on S and C (hue) accounting for the area-based joint PDF (saturation, from white for zero area coverage to the full color for the maximum area coverage), with area-based marginal PDFs ofS andC colored by the average alignment conditioned on each variable. Equispaced contour lines for 0.2 (dotted), 0.4 (dash-dotted), 0.6 (dashed), and 0.8 (solid) of the peak joint PDF value. 78 Figure 3.20: Surface distributions of the alignment between the scalar gradient and the eigenvectors;; of the strain-rate tensor and the vorticity! and the magnitude of the scalar gradient conditioned on local dimensionless curvedness (left) and absolute shape index (right) of 4 in STI M = 1:5 at dierent down- stream locations: k 0 x2 [2:2;1:8] ( ), k 0 x2 [4:0; 4:4] ( ) k 0 x2 [8:0; 8:4]( ). Lines represent ensemble mean of the conditional pdf, while colored region enclose two standard deviations. index, S. A decreased alignment with downstream of the shock is more pronounced at parabolic points (S = 0.5) of the isosurfaces. For the scalar gradient alignment with the most extensive strain-rate eigendirection , the correlation with S are less pronounced than with C: at the further downstream locations there is a slightly positive correlation with the shape index in its hyperbolic range (S < 0:5). The same correlation is found for the alignment with the intermediate eogendirection farther downstream. The ordering of the ensemble means by the three considered locations is according to the downstream evolution of their surface-averaged values, §3.6. Contrary to the and alignment withr, the alignment with the most compressive eigendirection shows a negative correlation with the absolute shape index for hyperbolic points. In the elliptical range of S, however, the correlation is positive, but becomes less pronounced further downstream. The alignment betweenr and ! shows a positive correlation when conditioned on C and with S in its hyperbolic range. Hence, the anti-correlation betweenr and! is more pronounced on atter surface 79 regions. The correlation between the magnitude ofr and the local geometry exhibits similar trend to found for the alignment betweenr and . This implies that most of the dissipation after the transient and its enhancement while crossing the shock occurs on atter regions (lower ^ C). Again, the ordering of the ensemble means by location is in agreement with the downstream evolution of the surface-averaged scalar gradient magnitude described in 3.6. Similar correlations are found in the DHIT and STI M = 3:0 congurations and for all initial structure sizes (not shown). The correlation between the Q-criterion on the passive scalar isosurface and the surface geometry shows similarities to the ndings for the vorticity magnitude (not shown). At the most upstream location Q is positively correlated to C over its entire range. At locations further downstream this correlation becomes more pronounced for C < 2:5. A slightly negative correlation is established for highly curved regions, but the expected value of the Q-criterion is still signicantly higher than for at regions. When conditioned on S Q shows a positive correlation in the hyperbolic range of S and for the more downstream locations peaks at the parabolic value ofS = 0:5, which is also observed for the vorticity magnitude. Note that the ordering of the ensemble means conditioned onS by downstream location is opposite for the Q-criterion than for the vorticity magnitude. 3.8 Trajectories of primary structures in the feature space 3.8.1 Geometrical trajectories and surface area evolution The deformation of passive scalar structures during their temporal evolution can be quantied by investi- gating the structure's trajectory (or path) in the geometrical feature spacef ^ S; ^ C; ^ g. Here, we focus on the temporal evolution of primary structures. Hence, from the graph representing the evolution of all tracked structures only branches whose start vertex results from a creation event and end vertex is a disappearing structure are considered. In the DHIT conguration only the initialized structures and in the STI only the incoming structures are candidates for primary structures, but depending on the structure interaction only a subset of them are truly primary structures. At initialization, all 4 passive scalar structures in the DHIT are spheres and thus their trajectories have a common starting point in the feature space atf ^ S; ^ C; ^ g = (1; 1; 1) (see gure 3.21). As visualized in gure 3.2 (a)-(c), the isosurfaces become highly convoluted shortly after the initialization, causing a rapid increase in the dimensionless curvedness ^ C and a decrease in the absolute shape index ^ S, due to the occurrence of hyperbolic surface points where S < 0:5. The deviation form the spherical shape also lowers the compactness, ^ , of the structures, as they are stretched. After reaching the maximum dimensionless curvedness, all structures follow a nearly linear trajectory in the ^ C- ^ plane towards the (0:7; 0:25) location, thus lowering both their dimensionless curvedness and compactness. The absolute shape index, ^ S, follows a rapid initial transition from 1:0 (spherical shape) to 0:5, reached as the maximum ^ C is attained. The subsequent decrease of ^ C and ^ has a lesser eect on ^ S, indicating a balance between the area coverage of hyperbolic and elliptic geometries on the convoluted surfaces. 80 From the start of the trajectories to the occurrence of the rst non-continuation event the change in the geometric signature from one frame to the next (given by the distance between consecutive trajectory points corresponding to tracking steps) monotonically decreases due to the viscous decay of the underlying turbulent background eld. As the stretching increases ( ^ < 0:5) rst split events occur in the structure evolution. The highly stretched structures tend to split over a sequence of frames. However, despite this cascade of split events, the change in geometrical signature of the structure is comparable to that of a continuation event, implying that these early splits do not alter the geometry signicantly. As visualized in gure 3.22(a), the split events that occur in this early stage of the structure evolution tend to be dominated by one target structure, similar in shape to the splitting source structure, and other smaller targets. As the dominant target continues the branch of the primary structure, its geometrical signature is moderately changed by the split event. This geometrical similarity between the source and the dominant target structure is also observed in the joint histograms of geometrical attributes (see gure 3.14 (a)). After these cascades of split events, the structures have reached their minimum compactness (maximum stretching), and the jump in length between consecutive tracking steps in the trajectory in the geometric feature space is increased as subsequent split events tend to aect the geometrical signature more signicantly (gure 3.22 (b)). The later generations of splits do not form a cascade, even though the branch section between splits remain short in some cases. The structure evolution in between splits is again characterized by a decreasing dimensionless curvedness and compactness. The later generations of splits aect the geometric signature of the primary structure in accordance with the ndings in the joint histograms of split events in the feature space (see §3.5): while the dimensionless curvedness remains relatively constant, the compactness and by trend the absolute shape index increase during the split event. The nal stage of the lifetime of primary structures is characterized by continuation events with increasing area coverage of elliptical surface points and increasing compactness, which makes a further splitting of the structure less likely. This trend is in accordance with the downstream distribution of split compared to disappearance events, discussed in §3.3. From the collection of individual geometrical trajectories, we calculate ensemble statistics resulting in the mean trajectory and its standard deviation (at each tracking step) for structures of each scalar eld ( ; 2 ; 4 ) in DHIT and STI (M = 1:5; 3:0) congurations. These ensemble trajectories, shown in gure 3.23, enable a characterization of the eects of the initial size and the shock compression on the geometric evolution of the structures. Whereas all structures start from a spherical shape ( ^ ; ^ S; ^ C) = (1; 1; 1) in DHIT, the initial increase in curvedness is signicantly less pronounced for the 2 than for the 4 structures. Furthermore, the mean ensemble trajectory of structures does not exhibit an initial increase in ^ C at all, since the initially smaller structures are less convoluted by the background turbulence eddies. Due to the dierences in surface convolution and stretching, the mean area evolution of primary structures in DHIT, shown in gure 3.24, also depends on their initial size, but suggests a certain degree of self-similarity. Exponential area growth has been proposed in the literture for the evolution of Lagrangian structures in statistically stationary turbulence (Etemadi, 1990; Girimaji & Pope, 1990, 1992; Goto & Kida, 2007; Yang 81 Figure 3.21: Trajectories in the geometric feature space of individually tracked primary 4 structures in DHIT. Markers indicate start (lled) and end (empty) points of sections of the branch where only continuation events occur: denotes a section start resulting from a split event, indicates a splitting section end. Half- lled markers are used if the section consists only of one frame. Analogously, and are used for merge and compound events, respectively. (a) (b) Figure 3.22: Higher-order split of 4 structures (translucent) in DHIT with one dominant target structure (blue) occurring at a earlier stage of a structure's evolution (a) and a higher-order split with two dominant target structures (in red and green) occurring at a later stage of a structure's evolution (b). 82 Figure 3.23: Ensemble mean trajectories of primary structures of (green), 2 (blue), and 4 (red), in the ( ^ S; ^ C; ^ ) geometric space, for DHIT (top), STI with M = 1:5 (center), and STI with M = 3:0 (bottom) cases. Arrows at each tracking time step indicate the direction of time. Error bars indicate one standard deviation of ensemble trajectories in each geometric parameter, plotted every eighth tracking step for clarity. The yellow circle marks the mean shock-crossing tracking step for the STI cases. Each trajectory considers tracking steps relative to its own creation event. 83 (a) (b) Figure 3.24: Temporal evolution of the normalized surface area (solid) of primary (green), 2 (blue) and 4 (red) structures in DHIT corresponding to STI M = 1:5 (a) and STI M = 3:0 (b). Errorbars include one standard deviation, plotted every 8-th data point. Dashed line represent tted curves of exponential growth and decay. Yellow marker indicate the range used for the curve tting. STI with M = 1:5 STI with M = 3:0 DHIT STI DHIT STI 2 4 2 2 4 2 ( x d x s )= 10.17 28.99 54.2 7.78 15.33 12.9 31.5 64.6 5.9 11.94 ( x c;s x s )= 3.61 9.93 15.71 3.64 7.34 6.5 12.9 23.9 3.11 5.5 ( x c;e x c;s )= 1.09 2.63 7.40 0.81 3.07 0.51 0.94 2.1 0.39 1.6 Table 3.5: Normalized mean disappearance location (top row), mean start location of split cascade (center row) and mean spatial length of split cascade (bottom row) for primary structures of type ; 2 and 4 in STI for M = 1:5 and M = 3:0 and corresponding DHIT ow congurations. et al., 2010). In the present study, despite the temporally decaying turbulence and diusion acting on the (non-Lagrangian) passive scalar structures, an exponential time dependence of the form A(t= 0 )=A 0 = a exp(bt= 0 ) can be tted to the early, stretching-dominated surface area evolution, with a and b constant. A larger initial area growth is observed for initially larger structures. The area growth rate is larger for the DHIT with a slightly higher initial Taylor-Reynolds number (Re = 46 compared with 42, corresponding to in ow turbulence used for STI cases with M = 1:5 and 3:0, respectively). After the initially exponential growth, the area growth rate decreases as the dissipation takes over the surface stretching and the surface area reaches its maximum. In the following diusion-dominated regime, a exponential decay can be tted to the area decrease. The area decay rate is faster for more stretched surfaces with a larger surface area (i.e., 4 structures). Cascade splits, which potentially shrink the surface area signicantly, occur in this period of exponential surface area decay. As the initial structure size increases (from to 4 ), the ensemble mean trajectory reaches lower values of the shape index, ^ S, favoring hyperbolic points, and compactness, ^ , indicative of an increased convolution. 84 (a) (b) Figure 3.25: Streamwise evolution of the surface area of primary (green) and 2 (blue) structures in DHIT (solid) and STI (dashed) normalized by its initial DHIT value. STI M = 1:5 and corresponding DHIT shown in (a), while STI M = 3:0 and corresponding DHIT in (b). Errorbars include one standard deviation, plotted every 8-th data point. Time coordinate in DHIT is transformed into spatial coordinate using equation 3.3. With increasing initial size, also the mean structure lifetime increases, as it takes longer for the structure to be diused below the threshold used for iso-surfacing. Increased lifetime translates into a mean disappearance location which is farther downstream after translating time to streamwise coordinate (equation 3.3). The lifetime reduction indicates that the mixing of larger structures benets more from the shock-interaction: for the shock interaction with M = 1:5, the mean disappearance location of structures is reduced by 24%, while for 2 structures the reduction is 47%. For the higher Mach number case, the reductions are 54% and 62%, respectively. The dependence of the lifetime reduction on the structure size correlates with that of the surface-averaged scalar gradient magnitude on the structure size. Furthermore, the cascade splits start at a later streamwise location and are spread over a larger streamwise length (see table 3.5). As the pre- and post-shock mean streamwise velocity diers between the two considered STIs, the streamwise locations of the corresponding DHITs dier in table 3.5. At the end of their lifetime, structures of all scalar elds considered tend towards the same geometrical attributes: moderate compactness and dimensionlesss curvedness with predominantly elliptical points. In the STI ow congurations, the structures have already been deformed individually at the time when they become closed surface. Therefore, their trajectories in the feature space do not have a common starting point and dier based on initial structure size (see center and bottom rows of gure 3.23). The more upstream half of the structure is less convoluted than the more downstream half when the structures become closed as the downstream half already spent more time within the domain. This discrepancy between front and rear half is more pronounced for bigger structures (as it takes longer for them to be closed) and for the lower Mach number case M = 1:5 due to the lower pre-shock velocity. As a consequence, especially for the 4 structures the maximum curvedness is less than in the DHIT as the more upstream part of the structure is less convoluted than the more downstream part and thus overall the maximum curvedness is less than in the 85 DHIT. The shock interaction lowers the lifetime of the primary structures in comparison to the corresponding DHIT as the mixing is enhanced (table 3.5). Additionally, the shock induces an earlier starting point of the cascade split events. Compression of the structures by the shock wave lowers the dimensionless curvedness ^ C ( attening the surface normal to the compression) and the compactness ^ (higher stretching, thus favouring split events) compared to the DHIT. The deformation induced across the shock, particularly for M = 3 is more pronounced for the smaller structures, as a larger portion of the structure is compressed by the shock at once. For the STI M = 3 case, the curvedness of the and 2 structures reaches a minimum after the structures are completely processed by the shock, and increases afterwards due to shock-induced turbulence amplication. The amplicaton of turbulence across the shock also leads to a surface area re-growth of 2 structures in the post-shock region that is stronger for the higher Mach number (gure 3.25). This period of area re-growth correlates with the decreasing behavior of the compactness parameter, ^ , in the feature space trajectories (gure 3.23) after the structures have been processed by the shock. While the amplication of the surface-averaged vorticity magnitude is similar for all initial structure sizes (gure 3.17), the area re-growth is larger for 2 than for structures. At this post-shock location the size of the structures is smaller than the Taylor microscale and, therefore, those structures are mainly advected by the (amplied) eddies. In contrast, 2 structures are still large enough to be also stretched by the (smaller) turbulence eddies. Prior to the surface area re-growth in the post-shock region, the higher Mach number compression decreases the surface area more signicantly. Structures of 4 exhibit a similar behavior up to a certain post-shock location where consecutive structures in the streamwise direction, decelerated by the shock, merge with each other and pile up, generating structures of large streamwise extent. A bigger initial spacing of the structures is needed to avoid this pile up process. The cascades of split events start at slightly higher values of dimensionless curvedness and end at slightly lower values of the compactness. Even though the cascade (sequence of consecutive splits) ends there, the structures are still highly stretched, thus further splits are likely to occur at later times of their evolution. In contrast to the DHIT, a number of splits (or compound events) occur for the STI cases, which aect the geometric signature of the primary structures early in their evolution. The number of primary 4 structures in the STI congurations with very long lifetime is too small to be statistically signicant and therefore the nal steps of their lifetime are excluded from the mean ensemble trajectory. Therefore those mean trajectories do not reach similar values of the geometrical properties as the other two mean trajectories at the end of the trajectory. 3.8.2 Trajectories with scalar dissipation and vorticity The study of ensemble trajectory statistics is extended to include the evolution of two physical quantities averaged on the surfaces: the magnitude of the scalar gradient (indicative of the scalar dissipation on the surface) and the vorticity magnitude, shown in gure 3.26 for the DHIT case. For 4 structures, the scalar gradient magnitude exhibits a rapid initial increase, due to stretching, of the initial scalar transitional 86 Figure 3.26: Ensemble mean trajectories (solid) of primary structures of (green), 2 (blue), and 4 (red) for DHIT in the ( ^ S; ^ C; ^ ) geometric space with scalar gradient (top) and vorticity magnitude (bottom). Physical quantities normalized by their initial value in DHIT. Error bars include one standard deviation of ensemble trajectories in each geometric and physical parameter, plotted every eighth tracking step for clarity. For and 4 structures the individual closest-to-mean trajectory (dotted) is plotted with highlighted steps (yellow markers) corresponding to the structures rendered in gure 3.29. region, as previously observed in the streamwise evolution (gure 3.17). The peak in curvedness is reached signicantly earlier than the peak in the scalar gradient magnitude: while the curvedness already decreases, the scalar gradient magnitude still increases. The geometry of 4 structures, especially the curvedness and the absolute shape index, develops much faster than both the scalar gradient and the vorticity magnitude at early stages of the trajectory. In contrast, for the smaller structures, the decrease in scalar gradient magnitude is initially more pronounced than the geometric changes. During the cascade split events, not only is the geometrical signature of the primary structure marginally aected, as previously explained, but also both surface-averaged physical quantities. Since the turbulence is already fully developed at structure initialization, the surface-averaged vorticity magnitude exhibits a monotonic decay in DHIT. Due to their longer lifetime, the vorticity magnitude reaches lower values on the larger 4 structures. The individual trajectory closest to the ensemble mean in terms of geometry and scalar gradient magni- tude (and thus, scalar dissipation) is also shown in gure 3.26 for 4 and structures. Figure 3.29 shows the structures corresponding to four highlighted steps of the closest-to-mean trajectory rendered with the alignment between the scalar gradient and the eigenvector of the strain-rate mapped on the surface. The area-based joint and marginal probability density functions ofS andC for the structures at each highlighted step are also shown, indicating that for both structure sizes (and for all four highlighted steps) the largest curvedness on the surfaces is found at parabolic points (S = 0:5) that occur on the ridges of the convoluted 87 Figure 3.27: Ensemble mean trajectories of primary structures of (green), 2 (blue), and 4 (red) for STI with M = 1:5 in the ( ^ S; ^ C; ^ ) geometric space with normalized scalar gradient (top) and vorticity magnitude (bottom).The yellow star marks the mean shock-crossing tracking step. Other plot elements are the same as in gure 3.26. surfaces. This is in agreement with Girimaji & Pope (1992) who reported that highly curved surface ele- ments of propagating surfaces in isotropic turbulence are of cylindrical shape. The correlation, found in §3.7, between atter surface regions and close alignment between the scalar gradient and the most compressive strain-rate eigenvector, , emerges shortly after the ad-hoc spherical initialization, weakening later in the structure lifetime. The ensemble mean trajectories of primary structures in STI with M = 1:5 dier from those in DHIT due to the open structure evolution close to the inlet and the impact of the shock on the geometrical and physical surface properties (see gure 3.27). As previously noted, the number of primary 4 structures with very long lifetime is too small to be statistically signicant and therefore the nal steps of their lifetime are excluded from the mean ensemble trajectory. Hence, the mean 4 trajectory does not reach scalar gradient and vorticity magnitudes as low as in the mean trajectories of and 2 structures. While the increase in scalar gradient magnitude depends on the size of the primary structures, the increase of the vorticity magnitude is the same for all structure sizes. These ndings are consistent with gure 3.17, where all structure types, not only primary ones, have been investigated. The stronger compression and quicker closure of incoming structures in the STI M = 3:0 conguration change the ensemble mean trajectories in terms of geometrical and physical structure properties compared to the lower Mach number case (gure 3.28). The smaller the initial size of the structures, the lower their curvedness is at the point of highest scalar dissipation. Compared to the DHIT (gure 3.26), 4 structures reach signicantly higher average scalar dissipation, and at lower curvedness resulting from the attening 88 Figure 3.28: Ensemble mean trajectories (solid) of primary structures of (green), 2 (blue), and 4 (red) for STI with M = 3:0 in the ( ^ S; ^ C; ^ ) geometric space with normalized scalar gradient (top) and vorticity magnitude (bottom).For the and 4 structures the individual closest-to-mean trajectory (dotted) is plotted with highlighted steps (yellow markers) corresponding to the structures rendered in gure 3.30. Other plot elements are the same as in gure 3.27. normal to the shock. Considering the evolution after the peak in scalar gradient magnitude, the curvedness of the and 2 structures decays in DHIT, is nearly constant for STI with M = 1:5, and has a period of increase for the STI with M = 3:0. Since the pile-up process of 4 structures farther downstream of the shock generates structures of large streamwise extent, the scalar gradient and vorticity magnitude do not evolve toward the low values observed for smaller ( , 2 ) structures in the nal stage of their ensemble mean trajectories. The average alignment between scalar gradient and the most compressive strain-rate eigenvector condi- tioned on S and C at highlighted steps of the closest-to-mean trajectory (gure 3.30) conrms again the correlation between atter regions and high alignment for both initial structure sizes. 89 (a) (b) (c) (d) (e) (f) (g) (h) (i) (j) (k) (l) (m) (n) (o) (p) Figure 3.29: Renderings of primary 4 (a-d) and (i-l) structures in DHIT with the alignment between scalar gradient and the eigenvector of the strain-rate mapped on the surface, corresponding to the high- lighted steps of the closest-to-mean trajectory in gure 3.26 (left to right, not to scale). The corresponding distribution of average alignment of scalar gradient and conditioned onS andC (hue), accounting for the area-based joint PDF is shown underneath each rendered structure (e-h,m-p). The horizontal direction in the renderings corresponds to the average shock-normal direction. 90 (a) (b) (c) (d) (e) (f) (g) (h) (i) (j) (k) (l) (m) (n) (o) (p) Figure 3.30: Renderings of 4 (a-d) and (i-l) in STI M = 3:0 and corresponding distribution of average alignment of scalar gradient and conditioned on S and C (hue) accounting for the area-based joint PDF (e-h,m-p) corresponding to the highlighted steps of the closest-to-mean trajectories in gure 3.28. Other plot elements are the same as in gure 3.29. 91 Chapter 4 Vortical structures in temporally developing compressible mixing layers 4.1 Introduction The mixing of two parallel streams and, hence, the growth rate of the mixing / shear layer forming between them is substantially inhibited at high Mach number compared to its incompressible counterpart. Statistical analyses have shown that compressibilty has a stabilizing eect which reduces the turbulence levels and growth rate of the mixing layer (Papamoschou & Lele, 1993; Cambon et al., 1992; Friedrich & Bertolotti, 1996). To measure compressibilty a set of dimensionless parameters is used. Denoting the free stream velocity, density and speed of sound in the two free streams as U i ; i ;c i , for i = 1; 2, the convective Mach number, M c = (U 2 U 1 )=(c 2 +c 1 ), is a measure of the compressiblity of the mixing layer, where equal specic heat capacities are assumed (Bogdano, 1983). Other compressibility measures are the turbulent Mach number, M t =u rms =c, whereu rms is the root mean square velocity, and the gradient Mach number,M g =Sl=c, where S denotes the mean velocity gradient and l is a characteristic length scale in the direction of the gradient. The gradient Mach number was found especially useful to distinguish the stronger compressibility eects in shear layers from those in boundary layers (Sarkar, 1995). In the literature, the canonical ow conguration of a mixing layer has been studied in terms of a temporally or spatially developing mixing layer, focusing on the dynamics of large, organized structures (Brown & Roshko, 1974), process of transition to turbulence and the self-similar state of the shear layer including passive and reactive scalar transport (Clemens & Mungal, 1995; Karasso & Mungal, 1996). From perturbed laminar initial (or inlet) conditions the ow develops to a turbulent, self-similar state via Kelvin- Helmholtz instability (Moser & Rogers, 1993; Rogers & Moser, 1994). During the transition to turbulence the occurrence of dierent structures, such as rollers and rib vortices (Moser & Rogers, 1991; Rogers & Moser, 1992) and -vortices, harpin-vortices and \ ower" structures (Zhou et al., 2012) has been observed. The dynamics of these structures have been found to be sensitive to initial (inlet) conditions. Statistical features of the turbulence in its self-similar state, such as the longitudinal velocity structure functions and scaling 92 exponents, follow the observations for homogeneous isotropic turbulence for small scales, but dier from those for larger scales (Attili & Bisetti, 2012). Linear stability analysis revealed that the mixing layer becomes truly three dimensional above a convective Mach number of M c = 0:6 (Sandham & Reynolds, 1990). The formation of shocklets in three-dimensional compressible mixing layer was reported for a convective Mach number of M c = 1:2 (Vreman et al., 1995; Fu & Li, 2006). However, Zhou et al. (2012) observed shocklets in a spatial developing mixing layer at a lower convective Mach number of M c = 0:7. It is of interest how the phenomenology of a decreased growth rate of mixing layers with increased compressibilty can be explained from a standpoint of turbulence structures. To gain insight into these structural mechanisms, the geometrical analysis and the ow feature tracking, §2, are applied to vortical structures extracted as isosurfaces of the Q-criterion from the mixing layer. Despite of the introduction of the hairpin vortex concept by Theodorsen in the 1950s, the occurrence, persistence and signicance of such vortices in shear ow congurations is still subject of debate. Therefore, special interest exists to study the evolution of hairpin-lie vortices in the mixing layer. 4.2 Problem formulation and numerical setup 4.2.1 Setup of the ow conguration Two measures of the thickness of the mixing layer are commonly used in the literature: the vorticity thickness ! (t) = U j@~ u=@yj max (4.1) and the momentum thickness (t) = 1 0 U 2 Z 1 1 U 2 ~ u 1 U 2 + ~ u 1 dy (4.2) where the Favre-avearged streamwise velocity component is ~ u 1 = u 1 = and the free stream velocity dierence is denoted as U =U 2 U 1 . Those two thickness measures can be used as characteristic length scales to form the vorticity thickness Reynolds number Re ! = U ! = and the momentum thickness Reynolds number Re = U =. Besides the convective Mach number, M c , the density ratio between the free streams 2 = 1 was also found to aect the growth rate of the mixing layer. This density ratio may be aected by the gas composition and heat release, such as in high-speed propulsion applications. In high-speed ow, a reduced (momentum thickness) growth rate for unequal free stream densities has been linked to a shift of the dividing streamline towards the low-density side (Pantano & Sarkar, 2002). In the present application the equal-density case 2 = 1 = 1 is considered. Following the conguration setup of a temporally developing mixing layer by Pantano & Sarkar (2002) and Foysi & Sarkar (2010), the mixing of two equally but opposite free streams,U 1 = U 2 = U=2, is 93 L x L y L z Top sponge layer Bottom sponge layer U 1 = U 2 0 ;p 0 U 2 = U 2 0 ;p 0 (a) (b) Figure 4.1: (a)Illustration of the initial condition of the temporally developing mixing layer. Mean streamwise velocity prole and sponge boundary damping indicated in red.(b) The M c = 1:1 mixing layer at the beginning of its self-similar regime ( = 650). Planes with reduced opacity show the density eld (white: = 0:7 , black: = 1:1) where the center plane is articially displaced to the bottom for visibility.Vortical structures are identied as iso-surfaces of the Q-criterion and colored by the streamwise velocity component (blue: u 1 =1:0,red: u 1 = 1:0). presently studied by conducting direct numerical simulations at higher spatial resolution than in previous studies, followed by a geometric characterization and temporal tracking of the resulting turbulent structures. Figure 4.1(a) presents the simulation domain and setup. The ow is initialized such that the mean streamwise velocity component u 1 is prescribed by a hyperbolic tangent prole, while the cross-stream and spanwise mean velocity components are set to zero: u 1 = U 2 tanh y 2 ;0 ! ; u 2 = 0; u 3 = 0 (4.3) Pressure and density are initially uniform across the domain. The vorticity and momentum thickness Reynolds number have initial values of Re !;0 = 640 and Re ;0 = 160. The initial vorticity thickness is specied to be !;0 = 1, while the momentum thickness is ;0 = 0:25 initially. To trigger the transition to turbulence broadband uctuations, u 0 i , are added to the mean velocity components. Those uctuations are obtained from a random eld on which a isotropic turbulence spectrum (Ristorcelli & Blaisdell, 1997) is imposed as E(k)/k 4 e 2k 2 =k 2 0 (4.4) with k 0 = 48. Applying this initialization technique leads to a nonzero compressible velocity component derived based on a small M t expansion of the Navier{Stokes equations. Additionally, (incompressible) uctuations of the pressure eld and linearly related density and temperature uctuations are also initially 94 M c ~ L x ~ L y ~ L z N x N y N z ~ L y;core N y;core 0:3 345 560 86 1024 1127 256 299:7 908 0:7 345 480 86 1024 973 256 253:5 768 1:1 345 400 86 1024 931 256 253:5 768 Table 4.1: Grid properties of the mixing layer conguration for dierent convective Mach number M c simulation cases. The domain extent is normalized by the initial momentum layer thickness. present. To limit the velocity uctuation to the center region of the domain a damping cross-stream function, f(y), is applied in the form of a hyperbolic secant shape following Li & Fu (2003): f(y) =C sech 2 y 2 ;0 ! (4.5) with C = 0:25. The velocity components are initialized as u i = u i +f(y)u 0 i . Since the size of turbulent structures in mixing layers decreases with increasing convective Mach number (Lesieur et al., 2005), three convective Mach numbers of M c = 0:3; 0:7 and 1:1 are considered. Whereas periodic boundary conditions are utilized in streamwise and spanwise directions, non-re ective sponge layers are added to the bottom and top boundaries of the domain to avoid spurious acoustic re ections (Mani, 2012). Dierent sponge lengths of normalized cross-stream extent of ~ L s = 153:25; 113:25 and 73:25 are chosen depending on the convective Mach number of the mixing layer M c = 0:3; 0:7 and 1:1, respectively. All grids share the same spatial resolution in the uniformly-spaced core region of the domain (see table 4.1). The grids have a uniform spacing in stream- and spanwise direction. For the cross-stream core region of the domain,j~ yj < 150 (M c = 0:3) orj~ yj < 126 (M c = 0:7; 1:1), the grid spacing in the y direction is uniform, whereas a grid stretching based on the geometric series is applied towards the bottom and top boundary (gure 4.2). For theM c = 0:7 andM c = 1:1 congurations the beginning of the sponge layers coincides with the beginning of the grid stretching, whereas for the M c = 0:3 case the sponges start in the outer region of the core of the grid. 4.2.2 Setup of the tracking algorithm Vortical structures in the mixing layer, identied and extracted as detailed in §4.4, are individually tracked by applying the ow feature tracking methodology as described in §2. To increase the performance of the tracking algorithm, in the OBB distance constraint, see §2.5, a OBB tree representation of varying depth, g. 4.3, is used rather than always using the root node of the tree to represent the structures in the GJK- distance calculation. Beneting from an improved surface representation, the distances between the leaf nodes of a source-target correspondence pair are calculated and more false correspondences are rejected in the OBB distance constraint, which otherwise would be rejected by the computationally more expansive surface proximity constraint. If a distance between leaf nodes satises the constraint, computation of all other (not yet determined) distances is skipped. The level of the OBB tree used to represent a structure is 95 Figure 4.2: Cross-stream grid coordinates (a) and spacing (b) of the mixing layer domain for M c = 1:1 as a function of the grid point index, j in that cross-stream direction. Figure 4.3: A hairpin vortex (blue) extracted from the M c = 1:1 mixing layer at = 196 and its OBB tree representation (red) of increasing depth from left to right. To capture the bent shape of the structure, the level of the tree needs to be suciently high. based on the compactness parameter ^ and the absolute shape index ^ S: Larger stretching and an increased number of hyperbolic surface points are used as indicators that a higher level of the tree is needed to capture the bent surface shape correctly. The settings of the tracking algorithm are summarized in tables 4.2 and 4.3. The main dierence in the setting compared to the tracking of the passive scalar structures, see §3, are: (i) Since the Q-structures are not convoluted compared to the passive scalar structures, the diagonal of the structure OBB is used as a reference length for the radius correspondence search and also for the spatial constraints rather than the parameter. (ii) The conservation-based events constraints are based on the volume rather than on the estimated mass of the structure. (iii) The weights used for the dominant target/source selection for the construction of split, merge or compound events in the graph data structure are based on conservation (volume) and condence, i.e. the lifetime weight is not used for (sub-)merge events. 96 Nearest-neighbor search Radius search Correspondence constraints f x f y f z f ^ S f ^ C f ^ c NN ! NN l ref " R " R;rel ! R l ref " " rel ! OBB ! local 1=L x 1=L x 1=L x 1 1 1 0:5 0:5 d OBB 1:2 1:4 0:1 d OBB 0:075 0:15 0:1 0:3 Table 4.2: Scaling factors, condence value, and weight for the NN search. Tolerance and weight utilized in the radius search. Structure reference length, tolerance, relaxed tolerance and weight used for both OBB distance and surface proximity correspondence constraints, as dened in §2.4.2-2.5. Conservation Surface points Conservation (Sub-)Split/Merge Compound Merge, Sub-merge Compound Merge Split Cont. Cons. Conf. " " rel c min N pts " " rel c min " " " 0:8 0:7 0:5 350 0:6 0:6 0:75 0:5 0:5 0:25 0:8 0:2 Table 4.3: Settings of the event constraining and event building in the graph. Type of the constraint and event type the constraint acts on. Event constraints are applied in sequence from left to right. Weights (conservation and condence) used for the dominant target/source selection to build (sub-)split/merge events in the graph data structure. 4.3 Validation of self-similarity and growth rate of the mixing layer The temporal development of the mixing layer conguration ofM c = 0:7 is visualized in g. 4.4 by means of the density distribution on the z = 0 plane spanning in streamwise (x) and cross-stream (y) direction. Due to friction between the two streams the temperature rises in the mixing layer, which results in a decrease of density as the overall pressure is nearly constant. Local pressure minima associated with vortices and acoustic waves emanating from the mixing layer are hinted in the density visualization. Self-similar solutions can be derived for the equations and boundary conditions governing a turbulent mixing layer of suciently high Reynolds number (Townsend, 1980). The time necessary to the reach self- similar state is in uenced by the choice of method and its perturbation strength used to generate initial random uctuation in the ow eld. The requirements for the mixing layer to be in the self-similar state are linear thickness growth and the independence of the shapes of the mean velocity and Reynolds stress proles of time (or downstream distance in spatial developing layers) when scaled by the velocity dierence and local mixing layer thickness (Mehta, 1991). These three indications of a self-similar period in the mixing layer evolution are validated in the following for the three considered simulation cases. M c ss;s ss;e _ ;ss 0:3 350 700 0.01308 0:7 500 950 0.009572 1:1 650 1050 0.008739 Table 4.4: Start ss;s and end ss;e time of the self-similar regime for dierent convective Mach number M c . _ ;ss denotes the slope of the linear regression conducted in the self-similar regime. 97 (a) (b) (c) (d) (e) (f) Figure 4.4: Density xy-planes of the M c = 0:7 mixing layer conguration at dierent times: (a) = 4, (b) = 200, (c) = 400, (d) = 600, (e) = 1000, (f) = 1800. Frames (c)-(e) roughly fall within the self-similar regime of the mixing layer = 500 960. The gure is limited to the transversal core region of the domain, see table 4.1. 98 (a) (b) Figure 4.5: (a) Temporal evolution of the momentum thickness for dierent convective Mach numbers M c . Markers indicate evolution interval used for the linear regression in the self-similar period. (b) Mixing layer growth rate in self-similar period for present simulations ( ) and references: `Langley curve' by Birch & Eggers (1972) ( ), Barre et al. (1996) ( ), Freund et al. (2000) ( ), Bonnet et al. (1993) ( ) Fu & Li (2006) ( ), Pantano & Sarkar (2002) ( ), Zhou et al. (2012)( ). The growth of the mixing layer, expressed in terms of momentum thickness, is shown in Fig. 4.5a. A period of linear growth can be found in all three momentum thickness evolutions as indicated by the linear regression of high correlation coecient in this period. With increasing convective Mach number a longer time of transition is observed until the mixing layer reaches its period of linear growth. Again indicative of turbulence inhibitions by higher M c , the momentum thickness growth rate _ = 1 U d dt = 1 ;0 d d (4.6) in the self-similar period, _ ;ss , is lower for higher M c , table 4.4. Following Pantano & Sarkar (2002), the lower convective Mach number case, M c = 0:3, is taken as quasi-incompressible in g. 4.5b to determine the in uence of the Mach number on the mixing layer growth rate in the self-similar regime. For the M c = 0:3 case, the present DNS has a growth rate of _ ;ss = 0:01308, while Pantano & Sarkar (2002) report a slightly higher value of 0.016. While a good agreement with the reference data is found for M c = 0:7, the decrease in _ ;ss relative to the quasi-incompressible growth rate is less pronounced in the present DNS compared to other simulations. Better agreement could be reached by slightly modifying the tting-ranges of the M c = 0:3 and M c = 1:1 cases. However, the time periods leading to the highest amount of self-similarity in the velocity and Reynolds stress proles, discussed in the following, were chosen as the periods of self-similarity. The streamwise velocity proles averaged over homogeneous ow directions are considered in g. 4.6. The velocity is normalized by the free-stream velocity dierence and cross-stream coordinate is scaled by the momentum or vorticity thickness of the time the prole is extracted from. The individual proles taken from a collection of time instants throughout the self-similar regime show indeed a high degree of collapse. 99 Figure 4.6: Self-similar streamwise velocity proles for M c = 0:3 (top),M c = 0:7 (center) and M c = 1:1 (bottom). Time interval between proles 40 in order: , , , , , , , , , , , . Reference data by Bell & Mehta (1990) ( ), Pantano & Sarkar (2002) ( ) and Rogers & Moser (1994) ( ) for the (quasi-)incompressible mixing layer. 100 Figure 4.7: Instantaneous cross-stream (y) proles of streamwise turbulence intensity for the M c = 0:3 (right), M c = 0:7 (center), and M c = 1:1 (left) cases just after initialization ( ) and in the self-similar periods: = 360 680; = 521 922; = 680 1040 forM c = 0:3; 0:7; 1:1. Time interval between proles 40 in order: , , , , , , , , , , . Each prole is spatially averaged in the homogenous directions (x and z). Figure 4.8: Ensemble mean prole of streamwise, cross-stream and spanwise turbulence intensities and Reynolds shear stress forM c = 0:3: Present DNS ( ), Pantano & Sarkar (2002, ), Rogers & Moser (1994, ), Bell & Mehta (1990, ) and Spencer & Jones (1971, ). 101 A collapse of Reynolds stress proles R ij = u 00 i u 00 j : = u i u j : ~ u i ~ u j is considered to be a more sensitive evidence for self-similarity than a collapse of mean velocity proles (Rogers & Moser, 1994). Therefore instantaneous cross-stream proles of the streamwise turbulence intensity are examined for self-similarity in g.4.7. Instead of individual, ensemble mean proles of the streamwise, cross-stream and spanwise turbulence intensities and Reynolds shear stress prole are examined in g. 4.8. The present DNS shows good agreement with previous simulations and experiments. Peak turbulence intensities are found in decreasing order along the streamwise, spanwise, and cross-stream directions, respectively, in agreement with Pantano & Sarkar (2002). 4.4 Identication, extraction and general statistics of vortical structures To investigated the temporal mixing layer conguration, described in §4.2.1, from a structural viewpoint, the present analysis focuses on the inhibiting eect of compressibility, expressed as a set of convective Mach numbers M c =f0:3; 0:7; 1:1g, on the evolution and geometrical properties of vortical structures in the mixing layer. As the vortex identication method, the compressible version of the Q-criterion originally introduced by Hunt et al. (1988) is utilized. Like the -criterion by Dallmann (1983); Vollmers et al. (1983); Chong et al. (1990), 2 -criterion by Jeong & Hussain (1995) or the ci -swirl-strength criterion by Berdahl & Thompson (1993); Zhou et al. (1999), the Q-criterion identies vortices based on the velocity gradient tensor A ij =@V i =@x j . Comprehensive reviews of vortex identication methods can be found in Epps (2017) and Zhang et al. (2018). To summarize the motivation of the Q-criterion, the characteristic polynomial of the velocity gradient tensor det A ij I = 0, which can be formulated as 3 + P 2 + Q + R = 0 (4.7) is considered. Here I is the identity matrix and are the eigenvalues of the velocity gradient tensor. The tensor's invariants P;Q and R can be expressed as (Suman & Girimaji, 2010; Epps, 2017) P =tr[A] =S ii (4.8) Q = 1 2 P 2 tr[A 2 ] = 1 2 P 2 S ij S ji ij ji (4.9) R = 1 3 P 3 + 3PQtr[A 3 ] = 1 3 P 3 + 3PQ S ij S jk S ki 3 ij jk S ki (4.10) where ij = 1=2(A ij A ji ) is the rotation-rate tensor and S ij = 1=2(A ij + A ji ) is the strain-rate tensor. These two tensorial quantities are dened as the antisymmetric and symmetric part of the velocity gradient tensor A ij . The Q-criterion denes a vortex as a region where Q > Q iso , where nominally Q iso > 0. In 102 (a) (b) Figure 4.9: (a) Evolution of the iso-contour value Q iso used to extract individual Q structures from the mixing layers of varying convective Mach number M c . (b) Number of extracted structures based on the iso-value evolution. The yellow marker indicate the begin and end of the self-similar regimes of the mixing layers. (a) (b) Figure 4.10: Smoothed number of detected events of dierent types over time: (a) Creations (solid), disap- pearances (dashed), (b) splits (solid) and merger (dashed) of any order. incompressible ow, this condition is satised where the vorticity magnitude is larger than magnitude of rate of strain. The Q-structures of interest are obtained by iso-contouring the Q scalar eld. Due to the threshold- based character of the Q-criterion, the choice of the iso-contour value Q iso is crucial for the shape, spatial extent and number of Q-structures extracted from the ow eld at an instant. A low iso-value generates large, widely ramied structures not suitable for the tracking of individual structures, while a high iso-value shrinks the extracted Q-structures reducing their complexity and lifetime. To account for the temporally evolving character of the ow eld, a time-dependent iso-value is chosen: From the mean,hQi, and standard deviation, p hQ 2 ihQi 2 , of the scalar Q eld of each frame the iso-value used for extraction is formed as Q iso =hQi+n p hQ 2 ihQi 2 . Here, the standard deviation is found to be the dominant contribution toQ iso . For the lower convective Mach number case (M c = 0:3), a choice ofn = 6:0 was found to extract a satisfactory number of individual Q structures throughout the simulation and mainly structures of tube-like shape in the 103 Figure 4.11: Temporal evolution of the geometrical signature of all structures extracted during the early, pre-self-similar stage ( = 0 - 240) of the M c = 0:3 mixing layer conguration. Each sample is weighted by the inverse of the total number of structures of the frame from which the sample is extracted. self-similar regime and later stages of the simulation. Due to the inhibiting eect of compressibility on the turbulence intensities in the mixing layer, the lower Mach number case is the most prone to formation of large, ramied Q structures. Therefore, a lower choice forn would still avoid largely connected structures for the higher Mach number cases. However, for the sake of consistency and to limit the number of parameter in uencing the study, the same choice of n = 6:0 is made for M c = 0:7; 1:1. Note that, despite the same choice of n, at a given time instant the iso-value Q iso tends to be lower for the higher Mach number, as the standard deviation of the Q eld is lower, see g. 4.9a. After an initial transient, which is longer for larger M c , the iso-value reaches its maximum before the beginning of the self-similar regime for all considered Mach number M c and decays afterwards. Indicative of the expression of turbulence inhibition on the structural level, a larger number of structures is extracted throughout the simulation for lower Mach number, g. 4.9b. Similar to the iso-value evolution, a transient is present in the number of extracted structures and the temporal distribution of detected event types, g. 4.10. Based on the applied iso-contouring and extraction method, the number of structures decays in the self-similar regime for M c =f0:7; 1:1g and only in the M c = 0:3 case the number of extracted structures remains constant throughout the self-similar regime of the mixing layer. The peak event activity coincides with the peak in extracted structures (and iso-value) and is ordered by Mach number. The number of merger exceeds the number of split events only in the very early stage of the simulation. The initialized velocity uctuations, see §4.2.1, generate small, rather compact ( ^ > 0:6) Q-structures in the very rst frames of the simulations. While these structures grow and merge with others, they become more stretched shortly after initialization ( < 40 for M c = 0:3), see g. 4.11. The continued merging process ( = 40 - 120 for M c = 0:3) leads to the formation of large structures characterized by a low value of compactness and more elliptical than hyperbolical surface points ( ^ S > 0:5). These largely connected 104 (a) ! x =(U= !;0 ) -1 -0.5 0 0.5 1 (b) Figure 4.12: Selection of largely connected vortical structures extracted at = 40 (a) and = 80 (b) from the M c = 0:3 mixing layer conguration colored by streamwise vorticity and individual structure color, respectively. Periodic reconnection of domain boundary intersecting structures is applied in streamwise and spanwise direction. The highly stretched structures span across the center plane of the mixing layer. The shape of the vortices is reminiscent of -vortices observed near the inlet of spatially developing mixing layers (Zhou et al., 2012). (ramied) structures, visualized in g. 4.12, are reminiscent of -vortices observed near the inlet of spatially developing mixing layers (Zhou et al., 2012). As the structures span across the center plane of the mixing layer, they split into smaller, less ramied structures shortly after, when parts of the structures are pulled into the oppositely directed free streams. Another indicator of the inhibiting eect of compressibilty on the turbulent dynamics of the mixing layer can be seen in the statistics of the (anti-)alignment between the normal vector of the Q-isosurfaces and the vorticity vector, g. 4.13(a), as it takes longer for the vorticity to become approximately tangential to the surfaces after initialization. Consistent with the Reynolds stress proles, see gure 4.7 and 4.8, the majority of new vortices is created near the center plane of the mixing layer, g. 4.13(b), where the cross-stream location of a structure is dened as the mean of the maximal and minimal cross-stream coordinate of the structure's surface points. Less new vortices are created per frame of the self-similar regime with increasing convective Mach number M c . The geometrical variety of the Q-structures extracted during the self-similar regimes of the mixing layers is shown in g. 4.14. The majority of the structures is characterized by a tube-like shape of decreasing 105 (a) (b) Figure 4.13: (a) Temporal evolution of the ensemble mean of of the surface-averaged angle between the surface normal and the vorticity on the surface. Filled area encloses two standard deviations. After an initial transient, which takes longer for increasingM c , the vorticity becomes approximately tangential to the Q-isosurfaces. (b) Distribution of cross-stream locations of creation events during the self-similar periods of the mixing layers. Counts normalized by number of frames in the self-similar regime. curvedness with increasing convective Mach number M c . The tube-like structures show a large variation in stretching, which decreases with decreasing curvedness. 4.5 Creations neighboring existing vortices Newly created Q-structures in the self-similar regime of the M c = 0:3 mixing layer are geometrically charac- terized as generally of tube-like shape with a relatively higher compactness ^ , g. 4.15. A large portion of the newly created structures disappears below the iso-valueQ iso used for surface extraction in the particular frame within a few tracking time steps and without any interaction with other structures. Contrary to these short-lived, isolated structures, other structures participate in structure interactions shortly after their creations. Therefore, this section focuses on the rst event (other than continuation and disappearance) in the lifetime of newly created structures. Compound events rarely occur as such (at least in the self-similar regime) and created structures mainly have merge or split events as their rst structure interaction. Pre- sented in g. 4.15, created structures that are the source of a split in their rst event typically live longer than merging newly created structures and become more stretched tube-like before the split occurs. If the rst event in the lifetime of a newly created structure is a merge, the structure either merges into another structure (i.e., the newly created structure is a non-dominant source of the merge) or the new structure continues its life after another structure merged into it. While both scenarios occur, the rst one is found to be (much) more frequent. In this case, the structure is created in close proximity to an existing structure and merges into such shortly after its creation. 106 (a) (b) S Figure 4.14: (a) Geometrical variety of all extracted Q-structures with more than 100 surface points in the self-similar regimes of the mixing layers for dierent convective Mach numbers. Contour lines of the PDF indicate 0:1 (solid), 0:4 (dashed) and 0:7 (dotted) of the maximum count. (b) Exemplary structures extracted from theM c = 0:3 at = 344 case corresponding to the letters in the feature space plot. Structures colored by absolute shape index S. Structure renderings are not to scale. Figure 4.15: Newly created structures during the self-similar regime of the M c = 0:3 mixing layer and their geometric signature while participating as a source in the rst merge (red), split (blue) or compound(orange) event after their creation. Newly created structures which do not have any structure interaction during their (short-lived) existence are excluded from the gure. The size of the event source marker is indicative of the lifespan from the creation to the rst event. For merge events, the shows a merge into another existing structure, while a marker is used if another structure merges into the newly created structure. 107 (a) (b) (c) (d) (e) (f) Figure 4.16: Examples of patterns of a newly created structure merging with an existing structure extracted from the self-similar regime of the M c = 0:3 mixing layer. Number of instances of each pattern: a) 50, b) 4, c) 122, d) 3, e) 12, f) 11. Dashed nodes do not belong to the pattern instance itself but to its embedding of the entire tracking graph. Node shapes signal event types following g. 2.7 Figure 4.17: Creation of a new vortex of rather compact shape in the neighborhood of an existing straight, tube-like vortical structure. The merging of both structures yields a tube-like structure of bent shape. Temporal order from left to right. Green arrows indicate vorticity on the structures colored by tracking ID. To further support these ndings based on g. 4.15, the graph mining method, presented in §2.10 is applied. Due to the size of the tracking graph, this graph mining eort is conducted only on the graph representing all structures evolutions in the self-similar regime of the M c = 0:3 mixing layer rather than searching the graph of the entire data set for frequent subgraphs. While the graph miner detects many split-dominated patterns and also the previously mentioned short-lived structure evolution only consisting of a creation, a sequence of continuation varying in length and a disappearance, this section focuses on found patterns involving creation and merge events. Examples of detected creation-merge patterns and their statistical signicance are shown in g. 4.16. Consistent with g. 4.15, the patterns indeed exhibit merging structures shortly after their creation. Depending on the relative position and orientation of the merging structures, the interaction of an existing straight tube-like structure with a newly created structure can yield a bent shape of the resulting structure, g. 4.17. The existing, dominant structure often develops a additional branch (\nger") in the course of merging with the incoming (newly created) structure. During the merge process the incoming, merging structure undergoes a stretching, as seen in g. 4.18(b). The new branch of the merge target structure can be of straight or rather twisted shape depending on the stretching during the merging process. 4.6 Detection of hairpin-like vortical structures Considering the geometrial variety of the extracted vortical structures, see g. 4.14, the conceptual model of a hairpin vortex (Theodorsen, 1952, 1955; Robinson, 1991) is proposing a vortical structure of distinctive shape. While observed in other shear ow congurations, such as homogeneous shear turbulence (Rogers & Moin, 1987; Vanderwel & Tavoularis, 2011; Dong et al., 2017) and mixing layers (Brown & Roshko, 1974; 108 (a) (b) Figure 4.18: The incoming (newly created) structure often develops an additional \nger" on the existing structures. These addition are often found to be twisted as in (b) Rogers & Moser, 1992; Moser & Rogers, 1993; Zhang et al., 2019; Zhou et al., 2012), hairpin vortices are frequently reported in wall-bounded turbulence (Acarlar & Smith, 1987; Chac n et al., 1996; Zhou et al., 1999; Tomkins & Adrian, 2003; Adrian, 2007; Wu & Moin, 2009; Schlatter et al., 2014; Eitel-Amor et al., 2015). The characteristic elements of a hairpin vortex are the spanwise oriented head connecting a pair of quasi-streamwise, counterrotating vortices referred to as the legs of the hairpin (Adrian et al., 2000). The growth of hairpin vortices is related to Townsend's attached eddy hypothesis (Townsend, 1980), which states that the size of eddies are proportional to their distance from the wall. In wall turbulence the hairpin concept oers a mechanism to generate Reynolds shear stress, low-speed streaks and to transport vorticity of the mean shear at the wall outward (Adrian, 2007). Conceptually, the hairpin head undergoes a (small) upward (wall-normal) motion and lies further away from the wall as the leg pair. In agreement with observations of hairpins in boundary layers, the velocity eld (perturbed by smaller-scale turbulence) surrounding an individual hairpin extracted from the M c = 0:3 mixing layer displays an ejection in between the legs of the hairpin, see g. 4.19. The inner part of the head of such individual hairpin is dominated by hyperbolic surface points, while the outer head part and the tips of the legs are elliptical in shape, g. 4.20(a). To investigate the temporal evolution of hairpin-like vortical structures in mixing layers, an algorithm to lter hairpin-like structures from all vortical structures extracted from a turbulent ow eld is proposed in the following. The goal of the algorithm is to lter individual hairpin-like structures, whereas (depending on the iso-value used for iso-surfacing ) largely connected hairpins, g. 4.20(b), are extracted from the ow eld but rejected by the lter. While there is no exact mathematical denition of a hairpin vortex in the literature, hairpin-like structures are ltered by imposing a sequence of geometrical and physical constraints on each extracted vortical structure. A set of input parameters required by the algorithm is listed in table 4.5 with 109 u 2 =U 0.15 0 -0.15 ! x =(U= !;0 ) 1.4 0 -1.4 Figure 4.19: Plane cuts of the velocity eld surrounding an individual hairpin extracted at = 144 from the M c = 0:3 mixing layer crossing its center plane. The velocity vector is colored by the cross-stream velocity component, while the streamwise vorticity component is mapped on the surface. (a) 0.2 0.4 0.6 0.8 S (b) ! x =(U= !;0 ) 1.5 0 -1.5 Figure 4.20: (a) Absolute shape index mapped on the surface of an individual hairpin vortex extracted at = 144 from the M c = 0:3 case. (b) Large, connected hairpin vortices extracted at = 128 from the M c = 0:3 case. Streamwise vorticity component mapped on surface, while green arrows represent the vorticity at sampled surface points. Renderings not to scale. parameter values used for all applications in this work. Processing one structure at a time, the algorithm starts with an application of spectral clustering to the surface points of the structure. To identify candidates for the head and leg pair of the potentially hairpin-like structure, this spatial clustering is conducted with a prescribed number of three clusters. No scaling operation is performed prior to the spectral clustering. Structures having a spatial clustering where any of the three clusters has less than a user-dened minimum size s min are rejected. The three spatial clusters are assigned as a head candidate and two leg candidates based on the number of openings each surface cluster has: The head candidate has two openings, whereas each leg has one opening. If the clustering leads to any other combination of openings, the structure is rejected. Next, the OBB tree of each surface cluster is obtained, where the tree depth is one for a leg and two for the head to better capture its bent shape, see g. 4.21. The two OBB leaf nodes of each leg are used to dene the direction of the leg as the line connecting the centers of the leaf nodes and pointing away from the root node of the head candidate. The angle between the leg directions is computed and compared against a user-provided maximum angle leg . It is observed that a local coordinate systemf^ x OBB 1 ; ^ x OBB 2 ; ^ x OBB 3 g based on the major, medium and 110 Required Optional s min leg " dip d dip 250 =4 0:5 150 0:05 10 Table 4.5: Input parameter of the hairpin lter algorithm: s min minimum size of spatial clusters, leg maximum angle between the leg pair," maximum distance between two samples to be considered in the neighborhood of each other in the leg vorticity - cross coordinate clustering, minimum cluster size in the leg vorticity - cross coordinate clustering. Parameter used for the optional bimodality test of the leg cross coordinate distribution: dip signicance threshold andd dip interval merge distance of the UniDip algorithm. minor direction of the OBB of the entire surface may lead to a poor capturing of the direction of the legs if the structure is rather asymmetric due to unevenly long legs for example, see g. 4.21 (a). To correct this tilt a local coordinate system is constructed as follows: From the two leg directions a mean leg direction ~ l is obtained and projected onto the major-medium plane of the structure OBB as ~ l proj = ~ l ( ~ l ^ x OBB 3 )^ x OBB 3 . A orthnormal basisf~ x 1 ; ~ x 2 ; ~ x 3 g is formed using the unit projected mean leg direction ~ x 1 = ^ l proj , the minor direction ^ x OBB 3 of the structure OBB and their cross product ~ x 3 . Using this local coordinate system it is required for the head candidate to have at least one leaf node in between the legs in terms of the ~ x 3 coordinate. A number of extracted structures approximate a hairpin-like shape, but their leg pairs do not exhibit counter-rotating character which is essential for hairpin vortices, g. 4.22. Therefore, simple statistical constraints are imposed on the distribution of the ~ x 1 component of the vorticity on the candidate surface: (i) The mean of the distribution on the surface cluster is required to have opposite signs for the two legs due to the counter-rotation. (ii) The absolute value the mean is expected larger than the standard deviation. To include tolerance, the coecient of variation is required to be smaller than two. (iii) Since the vorticity changes orientation along the head, the absolute value of the mean on the head cluster is expected smaller than on both legs. A hairpin with straight, approximately parallel legs exhibits a bimodal distribution of ~ x 3 surface point coordinates of the leg pair, i.e., the distribution excludes the head cluster, g. 4.21 (a). To detect such bimodality the probability density function is obtained from the ~ x 3 coordinates of the legs using a kernel density estimation (KDE) where the bandwidth is optimized by utilizing a cross-validation approach. After adding noise to this estimated PDF, random samples are obtained from the distribution and fed into the UniDip algorithm (Maurus & Plant, 2016) to determine the modality of the distribution. However, hairpins are often distorted by surrounding eddies. For example, if a leg is bent inwards, the distribution of the ~ x 3 coordinates of the legs may not be bimodal. To detect those more asymmetric hairpins the bimodality is neither a necessary nor a sucient criterion. It is kept in the algorithm to optionally lter the subset of hairpins that satisfy this criterion, but not applied to obtain the results in the present work. Complementing the ~ x 3 coordinates, the ~ x 1 component of vorticity, ~ ! 1 , is considered as another feature of the leg pair. The z-scores of those two features of the leg pair are obtained independently and the standardized features are clustered by applying the DBSCAN algorithm (Ester et al., 1996; Schubert et al., 2017). In 111 (a) (b) Figure 4.21: Elements of the hairpin detection method from left to right: Leg component of the vorticity mapped on the extracted surface. Green arrows represent the vorticity at sampled surface points. The three clusters resulting from the spatial clustering of the surface are indicated in red, green and blue. The OBB trees of each surface cluster are colored according to the cluster they enclose. A tree depth of one is used for the legs, while the depth is two for the head cluster. Red and green arrows show the direction of each leg, while black arrows are the mean leg direction projected onto the major-medium plane of the global OBB of the surface (black cuboid). Modes in the distribution of the cross coordinates of the legs are highlighted as density intervals. Probability density function obtained from the KDE of the leg cross coordinates shown in black. Histogram obtained from the samples fed into the UniDip algorithm with added noise. Convex hulls enclose the clusters obtained from the clustering of the leg component vorticity - cross coordinate clustering of the legs. Cluster core points plotted with increased point size. Scatter points colored by the leg component of the vorticity. Range of color bar based on leg vorticity component of entire surface to also apply to the surface mapping (left). (a) showing a bimodal distribution of cross coordinates of the legs, while a unimodal distribution was detected for (b). Structures extracted from the M c = 0:3 case at = 148 and = 232, respectively. 112 Figure 4.22: Representative structures approximating a hairpin-like structure in shape, while diering in vorticity distribution on their surface from a hairpin vortex. Leg component of the vorticity mapped on the surfaces, while green arrows indicate vorticity. The structures are detected as non-hairpin-like based on their physical dissimilarity to a hairpin vortex. contrast to the earlier used spectral clustering, the number of clusters is not prescribed but estimated by the DBSCAN algorithm. Requiring a certain amount of compactness in ~ x 3 and ~ ! 1 , it is expected that each leg contributes one cluster to the data, i.e., the DBSCAN algorithm estimates two clusters in total. This approach is robust against smaller eddies which merged into the hairpin impairing not only its shape (~ x 3 ) but also its vorticity distribution due to the noise detection of the DBSCAN algorithm. Fig. 4.23 illustartes dierent types of rejections of vortical structures which passed all other prior constraints. As mentioned earlier, hairpins are often lopsided and the (quasi-streamwise) legs are only rarely counter- rotating pairs of equal strength leading to asymmetry (Robinson, 1991). The proposed hairpin lter detects those imperfect candidates, referred to as hooks and arches in g. 4.24, as hairpin-like. The number of detections of these imperfect hairpins can be greatly reduced by requiring the candidate head surface cluster to have an OBB root node center which has the utmost (smallest or largest) major coordinate ^ x 1 of the entire surface OBB compared to the OBB root node centers of the other two surface clusters. 4.7 Hairpins in the geometrical feature space The proposed detection method, §4.6, is applied to lter hairpin-like structures from all extracted structures of the compressible mixing layers considering all frames from initialization until the end of the simulations in the post-self-similar regime. The lter settings are listed in table 4.5. The detected hairpin-like structures are mapped into the geometrical feature space in g. 4.26. The hairpin vortices are characterized as stretched tube-like structures. While there is a number of detections showing a balance of hyperbolic and elliptic surface points ( ^ S 0:5), elliptic surface points slightly dominate for most cases due to the shape of the outer head region and the tips of the legs, g. 4.20 (a). For the lower convective Mach number (M c = 0:3) the range in geometric signature is in agreement with hairpins recently reported in stably stratied atmospheric boundary layers (Harikrishnan et al., 2021). Similar to the observations regarding structures of all shapes in the self-similar regimes, see g. 4.14, the curvedness of the detected hairpins decreases with increasing convective Mach number. The decrease in curvedness is mainly localized on the hairpin legs rather than 113 Figure 4.23: Representation of dierent types of rejections in the leg vorticity component - cross coordinate clustering step of the hairpin detection method. Convex hulls enclose detected clusters. Core and non-core cluster members are colored by the leg vorticity component, where non-core members are reduced in size. Noise points indicated in black. 114 Figure 4.24: Examples of hooks (top row) and arches (bottom row) detected as hairpin-like in the pre-self- similar stage of the M c = 0:3 mixing layer case. Other plot elements as in g. 4.22. (a) (b) (c) (d) C Figure 4.25: Dimensionless curvedness mapped on hairpins extracted from the M c = 0:3 (a, = 232; ^ C = 0:78), (b, = 360; ^ C = 0:82) and M c = 1:1 (c, = 464; ^ C = 0:63), (d, = 424; ^ C = 0:58) mixing layer. For the higher convective Mach number atter regions are found on the legs of the hairpin. on the head of the structure. Highly stretched hairpin detections are typically characterized by a medium curvedness compared to other detections of the same convective Mach number. Similar geometrical signatures are obtained for hairpins detected at dierent stages of the mixing layer evolution (pre-/post-/ self-similar), see g. 4.27. 4.8 Temporal distribution and formation of hairpin-like vortices The temporal distribution of hairpin-like structure detections during the mixing layer evolution is presented in g. 4.28. If an individual structure is detected to be of hairpin-like shape in multiple frames within the its lifetime, each detection contributes to the distribution. Throughout their lifetime Q-structures can develop and lose hairpin-like shape multiple times with or without interacting with other structures. Therefore the lifetime of a structure temporarily in hairpin-like shape can be split into sequences of detections in consecutive frames of the simulation. The sequences of consecutive detections vary in temporal length, however, the majority is found to be rather short-lived ( < 20, or even shorter). 115 (a) (b) (c) Figure 4.26: Hairpin vortices in the geometrical feature space. Hairpins detected in the M c = 0:3(a), M c = 0:7(b) and M c = 1:1(c) mixing layer. Markers of the same color indicate the same structure detected as hairpin-like in dierent (not necessarily successive) frames. Marker types are indicative of the time of extraction: Pre-self-similar , self-similar , post-self-similar . The lines , , correspond to 0.25, 0.5 and 0.75 contour level respectively. 116 Figure 4.27: Geometrical signatures of hairpin-like vortices detected in the pre-self-similar (top), self-similar (center) and post-self-similar stage of the M c = 0:3 mixing layer. Marker color indicative of tracking ID. 117 (a) (b) (c) Figure 4.28: Temporal distribution of hairpin detections in mixing layers of convective Mach number M c = 0:3 (top),M c = 0:7 (center) andM c = 1:1 (bottom). Number of detections normalized by the total number of extracted structures with at least 750 surface points of the time step. While the solid distributions include all detections, the dashed distributions are based on the rst in sequence detections. Markers on the temporal axis show self-similar periods. The event types forming the hairpin-like structures are included as stacked contributions to the rst in sequence distribution. 118 For all considered convective Mach numbers, the number of hairpin detections peaks in the pre-self-similar stage and decays afterwards. This concentrated occurrence is related to the breakup process of the large, connected structures, see g. 4.12 and g. 4.20(b), which were formed through (sub-)merge events in the very early stage of the mixing layer development. Consequently, split and compound event, i.e., sub-split events, contribute the most to the peak occurrence of individual, disconnected hairpins. Examples of this structural breakup process are rendered in g. 4.29 for the M c = 0:7 mixing layer. The breakup process is observed for connected structures within which the hairpin heads are located in the same free stream side as well as in opposite sides. The breakup of large, connected hairpins is in uenced by the choice of the iso-value used for extracting the Q-structures, i.e., a structure which is observed to breakup for the used iso-value of the frame, may stay connected for a lower iso-value. As discussed in §4.4, the M c = 0:3 case exhibits the most restrictive requirement for the choice of n in Q iso =hQi +n p hQ 2 ihQi 2 to prevent the formation of very large connected structures in the early mixing layer development. While a choice of n = 6 indeed prevents the formation of domain-spanning structures, a signicant portion of hairpin-like structure may still be connected for this choice of n and is therefore not detected as individual hairpin-like structures, which decreases the count of detections in g. 4.29(a). The number of connected hairpins (not accounted for in this analysis) need to be kept in mind when evaluating the peak count of detections for the dierent convective Mach numbers. For the given choice of n = 6 the normalized peak count decreases with increasing convective Mach number fromM c = 0:7 toM c = 1:1. With increasing convective Mach number the occurrence of the peak in the temporal distribution of the hairpin-like detections is delayed in time, but occurs approximately at the same momentum thickness Reynolds number of Re 400. In the self-similar regimes of the mixing layers, which are less prone to the formation of largely connected structures compared to the early stage of the mixing layer development, the normalized number of detections may slightly decrease with increasing convective Mach number. The reduced detection of hairpin-like structures for increasing momentum thickness Reynolds number Re supports observations of hairpins in turbulent boundary layers for which the hairpin regeneration process was reported to be not sustained once the turbulent background is developed (Eitel-Amor et al., 2015). The same authors nd no evidence of hairpin vortices reaching into the wall region for Re > 400, which matches with the momentum thickness Reynolds number of the peaks in the temporal distribution of the hairpin-like detections in g. 4.28. In the subsequent stages of the mixing layer evolution the detected, individual hairpin-like structures contribute less than 5% to the total number of structures with more than 750 surface points of the frame and are therefore not considered as a dominant ow feature. However, a number of hairpins connected to other (not necessarily hairpin-like) structures may still exist, but is not accounted for. Fig. 4.30 demonstrates how hook-like structures with quite unevenly pronounced legs can develop arch- or hairpin-like shape by either splitting the dominated leg or merging with another tube-like structure which serves as the second leg of the hairpin. Note, however, that hook-like structures can also develop more 119 (a) (b) (c) (d) Figure 4.29: Connected, partially hairpin-like structures, see also g. 4.12(b) and g. 4.20(b), breaking into individual hairpins during the pre-self-similar stage of the M c = 0:7 mixing layer. The breakup process is observed for connected structures within which the hairpin heads are located in the same free stream side (b, = 164), opposite sides (c, = 172) and (d, = 192) as well as often in mixed cases (a, = 104). Target structures articially displaced for visibility, renderings are not to scale. (a) (b) Figure 4.30: Hook-like structures developing in arch-like and hairpin-like shape by splitting the dominant leg (a) ( = 156) and merging with a tube-like vortex (b) ( = 252) in the M c = 0:7 mixing layer, not to scale. 120 Figure 4.31: Side view (xy) of a series of merge events generating a short-lived hairpin vortex in theM c = 0:3 mixing layer and subsequent breakup from = 304 (top,left) to = 340 (bottom,right). Structure color represents structure ID, while green arrows are the down sampled vorticity on the structures. Note the tilted orientation of the hairpin with respect to the center plane of the layer indicated in gray. hairpin-like shape in a sequence of continuation events without structure interactions by gradually growing or losing a leg. The formation of a hairpin vortex in a series of merge events is shown in g. 4.31. First, the merging of two tube-like structures generates a hook which develops a bent shape and eventually merges with another structure to form a hairpin. Note, that both, the hook and the more tube-like structure, grow towards each other in this merging process. Interestingly, the hairpin has a tilted orientation with respect to the center plane of the mixing layer. The resulting hairpin is short-lived as it loses its rst leg only one frame after its formation. 121 Figure 4.32: Section ( = 204 - 296) of a lifespan of a vortical structure developing hairpin-like shape in the M c = 0:3 mixing layer, from left to right and top to bottom. Interacting structures included as source or target structures of the respective event. The structure is mostly located above the center plane of the mixing layer indicated in translucent black. Several structure interactions contribute to the persistence of the hairpin in the rendered section and subsequent frames (not shown). 122 Figure 4.33: Trajectory of a vortical structure temporarily developing hairpin-like shape in geometrical fea- ture space enhanced with surface-averaged streamwise vorticity component. The trajectory section between points A and B corresponds to the renderings in g. 4.32. Detected haipin-like structures are highlighted as in the left subplot. In the center subplot markers , and indicate if a step of the trajectory is the target of a split, merge or compound event respectively. Fig. 4.33 presents the full trajectory of a structure which preserves its hairpin shape for a longer se- quence of frames in the geometrical feature space enhanced with the surface-averaged streamwise vorticity component. The temporal evolution of the structure is rendered for a subsection of its trajectory in the cor- responding g. 4.32. The structure originates from a creation event at = 80 and develops tube-like shape, while being orientated mainly in streamwise direction and intersecting with the center plane of the mixing layer. As the tube-like structure becomes more stretched, its streamwise vorticity component intensies. The structure develops hook-like shape as a result of a complex compound event and eventually becomes hairpin-like over a series of frames within which one leg intensies while the other leg is shortened in splits. At the time the structure is rst detected as a hairpin the streamwise vorticity component is almost balanced on the structure. The geometric signature of the hairpin detections is in agreement with those presented in g. 4.26. The sequence of consecutive hairpin detection ends in a split of a leg, which leads to an imbalance in the surface-averaged streamwise vorticity. However, the structure retrieves its hairpin shape in a merge with another structure serving as the new leg. 4.9 Cross-stream and spanwise growth of hairpin vortices The entire lifespan of two structures which temporarily develop hairpin-like shape in the M c = 0:3 mixing layer are shown in g. 4.34 and g. 4.35. Both structures remain close to the center plane of the mixing layer throughout their lifetime, often intersecting with it. Therefore, the streamwise component of the surface- averaged velocity changes its sign multiple times during the structure evolution. Resulting from a split, the 123 Figure 4.34: Lifetime ( = 280 - 376) of a structure temporarily developing hairpin-like shape in theM c = 0:3 mixing layer. Streamwise vorticity mapped on the structure and vorticity sampled as green arrows. Due to the proximity to the center plane of the mixing layer, the surface-averaged streamwise velocity of the surface often changes its sign. Articial structure displacement added for visibility based on the sign of this surface- averaged velocity: Starting from level A a transition upwards indicates a change towards an (on average) left moving structure, while a downwards transition indicates a shift toward a right moving structure. Size of the black arrow indicates the temporal order of these transitions. The yellow plane indicates the center plane of the mixing layer at each level of the gure. The 3D renderings on the right show the structures of each level of the left side of the gure (parallel projection). While the rendering are not to scale, the pink line serves as a reference scale (10 times the initial momentum layer thickness). The structures from which the hairpin-like structure originally split (source of (sub-)split) and eventually merges (target of (sub-)merge) are added with reduced opacity. Note that the head of the hairpin becomes quite tilted at the end of the lifetime. structure shown in g. 4.35 has a hairpin-like shape with the head located in the upper free stream. However, multiple structure interactions later, the structure has developed into a larger hairpin with the head located in the lower free stream. The distribution of the cross-stream location of all detected hairpin-like structures and the rst-in- sequence subset is presented in g. 4.36 for the three considered convective Mach numbers. Here, the cross-stream location of a structure is dened as the mean of its minimal and maximal cross-stream position. Following the Reynolds stress proles, g. 4.7 and g. 4.8, the distributions peak at the center plane and are about symmetrical. The distribution including only the rst-in-sequence detections shows slightly less variance and more detections occur closer to the center plane. Fig. 4.37 investigates the subset of the center plane crossing detected structures and their extent on each free stream side of the mixing layer relative to the momentum layer thickness at the time. The majority of the structures shows an uneven top and bottom cross stream extent, which is also supported by the previous renderings of detected structures. For hairpins whose legs are oriented in the streamwise direction and the 124 Figure 4.35: Lifetime ( = 96 - 172) of a structure developing temporarily hairpin-like shape while crossing the center plane of the mixing layer. As frequently observed, the hairpin loses one leg over time and becomes a hook-like structure. Other plot elements as in g. 4.34. head in spanwise direction, the dominant cross-stream side is the side on which the head is located. Fig. 4.37 includes structures which have a tilted orientation with respect to the center plane of the mixing layer. In turbulent boundary layers the cross-stream and spanwise growth rate are reported to be similar to each other and is associated with Townsend's attached eddy hypothesis (Adrian, 2007). Therefore, the hairpin growth behavior is investigated in the context of the mixing layer conguration, in g. 4.38. The minimal and maximal surface point coordinate in the cross-stream and spanwise direction is used to dene the extent of a hairpin structure in the particular direction, y and z. The evolution of each structure temporarily in hairpin shape is inspected visually to exclude structures that show a signicant tilted orientation and hairpins that are distorted by connecting to other structures at any point in their evolution. While no exact measure is used to exclude those structures, it is ensured that the spanwise extent is representative for distance of the outsides of the legs and the cross-stream extent is an indication of the distance between the head and the legs. Due to these requirements mainly hairpins occurring in the pre-self-similar regime are included in the analysis. While the growth behavior diers between individual hairpins, the overall trend for the M c = 0:3 case is in rough agreement with the ndings reported for hairpins in wall-bounded turbulence. The cross-stream and spanwise structure extents are similar to each other. The results are less clear for the M c = 1:1 mixing layer as a number of hairpins is found to grow in cross-stream direction while approximately preserving their spanwise extent. However, the overall ranges of cross-stream and spanwise extents are similar for both convective Mach numbers. 125 Figure 4.36: Distribution of the cross-stream location of all detected hairpins (solid) and the rst-in-sequence subset of hairpins (dashed) normalized by the momentum thickness. Compare to the Reynolds stress proles in g. 4.7 and g. 4.8 as well as the cross-stream distributions of newly created structures of all shapes in the self-similar regimes in g. 4.13 (b). Figure 4.37: The cross-stream extent of center plane crossing hairpins. Most hairpin-like structures emerge unevenly into the top and bottom free stream. On the basis of prior renderings, the dominant cross-stream side is the side where the head is located. 126 (a) (b) Figure 4.38: Growth of the cross-stream y and spanwise z extent of selected hairpin-like vortices in the M c = 0:3 (a) and M c = 1:1 (b) mixing layers. Colors indicate individual tracking structure ID, while temporal order is shown with arrows. For hairpin detections in non-consecutive frames dashed arrows are used. Note that the hairpins may be aected by structure interactions at any time, which in uence the extent of the structure. Hairpin detection in the pre-self-similar self-similar and post-self-similar stage plotted as , and , respectively. 127 Chapter 5 Conclusion A hybrid approach of attribute- and region-based local tracking has been introduced to study the temporal evolution of ow features, dened as closed surfaces of any kind in a turbulent ow eld. As a pre-step of the explicit tracking, ow feature are extracted and characterized for each frame of the given data set independently. The resulting spatial and geometrical feature attributes are utilized in the search for corre- spondences between feature in consecutive frames aiming to generate candidates of temporal feature relation. A sequence of region-based correspondence and attribute-based event realizability constraints has been pro- posed to test the physical plausibility of the detected correspondence candidates. A graph data structure is built to organize the evolution and record the interactions of individually tracked ow features over time. The graph data structure is analyzed and queried to gain insight of the ow feature dynamics by means of individual structure evolution as well as ensemble statistical evolution. Graph data mining techniques have been adapted to search for common pattern in the interactions between ow features. The tracking algorithm was applied to investigate the evolution of passive scalar structures in temporally decaying homogeneous isotropic turbulence and statistically stationary shock-turbulence interaction at Mach numbers 1:5 and 3:0. Here, the ow features were dened as iso-surfaces of the underlying passive scalar eld. The tracking application and the comparison of both ow congurations allowed to elucidate the eect of shock waves on the scalar mixing on nite-sized structure evolution and surface-mapped physical quantities relevant to mixing. Shock-induced deformation increases the occurrence of split and subsequent disappearance events, shifting those events closer to the shock for increasing Mach number. The structural breakup process facilitates the scalar mixing. As expected, initially larger structures take longer to diuse, but their mixing benets more from the shock-interaction as the number of splits and disappearances per incoming structure increases more (compared to the DHIT) and the upstream shift of event activity is more signicant than for smaller structures. Small-structure interactions were found to deviate from volume-preservation as a consequence of increased scalar dissipation at ne scales, whereas events involving larger structures show a high degree of volume 128 preservation (apart from larger shock-crossing structures which are compressed and aected by the shock- amplied scalar dissipation). Structures that have shrunk to a certain volume due to scalar dissipation or previous splits, tend not to split further, and simply continue until they reach the end of their lifetime and disappear by diusion. The geometric signature of structures involved in split and merge events show opposite trends: splitting (merging) structures are more stretched (compact) and show a larger (smaller) area of hyperbolic points than the resulting split (merged) structures. Compared to DHIT, shock-induced compression of the passive scalar structures increases their stretching and, in turn, the probability of the structure to split in the STI conguration. The enhancement of scalar dissipation by the shock is expressed at the structural level as an amplication of the surface-averaged scalar gradient. This amplication increases with initial structure size, indicating a better mixing enhancement by the shock-interaction for larger structures. Surface-averaged alignments of the scalar gradient with the most extensional strain-rate eigendirection increases across the shock, while the alignment with the intermediate and most compressive eigenvector decreases. This result, previously reported in volumetric statistical analysis of fully-developed scalar mixing, is here conrmed also at the structure level from scalar elds initialized with uniform patterns of well-dened shapes. Larger scalar gradient (hence, dissipation), its shock-induced amplication, and alignment with the most compressive strain-rate eigendirection were found to correlate with atter surface regions. These correlations emerged shortly after the initial convolution of the initially spherical structures by the background turbulence. The geometrical trajectories of the passive scalar structures depend on their initial sizes (relative to the eddy sizes of the background turbulent ow). The amount of surface convolution early in the structure's lifetime increases with initial size and leads to larger curvedness and area-coverage of hyperbolic points, resulting in an increased number of splits per structure. Shock compression alters the geometrical trajectories in the STI congurations, lowers the compactness and curvedness of the structures. Less compact (i.e. more stretched) structures are more likely to split, while a decreased curvedness shows increased area-coverage of at surface regions where scalar dissipation dominates. The sudden shock interaction attens the structures and enhances mixing more for larger structures. The post-shock amplication of turbulence increases the curvedness again, inducing a period of re-growth of the surface area and stretching of the structures, which increases the probability of splits and enhances scalar mixing. In a second application, the ow feature tracking was applied to study the dynamics and geometrical properties of vortical structures in temporarily developing mixing layers of varying convective Mach numbers M c =f0:3; 0:7; 1:1g to complement existing statistical analyses on the inhibiting eect of compressibilty on the turbulence intensities in such ow conguration. The compressible Q-criterion was used to identify vortices in the ow eld. The tracking algorithm was successfully applied with the same types of realizabilty constraints as in the passive scalar mixing. In contrast to the passive scalar application, a time-dependent iso-value was used for the structure extraction to account for the temporarily developing character of the 129 ow. With increasing convective Mach number less vortices were extracted which pronounces the inhibiting eect of compressibilty on the structural level. Higher compressibilty was observed to lower the curvedness of the extracted Q-structures. The application of the graph data mining approach focused on the creation of vortices in the neighborhood of existing vortices. The subsequent merging of both vortices is recognized as a prominent pattern of the vortex dynamics. Other patterns of statistical signicance are split-dominated or highlight the short-lived evolution of isolated structures without structure interactions in their lifetime. A detection method for hairpin-like vortices was introduced based on geometrical and physical constraints. In combination with the tracking, the hairpin lter was used to evaluate the occurrence, persistence and signicance of hairpin-like vortices in the mixing layers. For all convective Mach numbers, the number of hairpin detections peaks atRe 400 in the pre-self-similar stages of the mixing layers and decays afterwards, suggesting that hairpins are not a dominant ow feature once the turbulent background is developed. These ndings support previous reports of the hairpin regeneration process, which was found not to be sustained at higher Reynolds number in wall-bounded turbulence. It would be of great interest to evaluate what in uence dierent choices for the iso-value have on the number of detected hairpins, especially for the lower Mach number case. The time span in which vortices preserve their hairpin-like shape varies, but most hairpins are rather short-lived as they cross the center planes of the mixing layers and are partially pulled into the opposite free streams. This mechanism could explain the reduced lifetime of hairpins in the mixing of two opposite free streams compared to those found in boundary layer ow. The majority of individual hairpins is formed by the structural breakup process of largely connected hairpins which originate from vortices that grow in the early stage of the mixing layer through merging of tube-like vortices. A number of hairpins have a tilted orientation with respect to the center plane, especially in the later stages of the mixing layer. Hairpins with quasi-streamwise orientation exhibit equal rates of simultaneous growth in cross-stream and spanwise direction. Future extensions of the presented ow feature tracking methodology could target to track features not of only one type, but of multiple types simultaneously. This would allow to study relations between the dynamics of ow feature of dierent types. The graph mining technique would greatly benet from a pattern classication method to summarize the pattern extracted from large graphs better. To increase the performance of the algorithm, parallelism could be added to the implementation, especially to the local surface distance computation that is currently the computational bottleneck of the method. Another mode can be introduced to the algorithm in which it does not operate as a pure post-processing tool of a given data set, but in which the tracking is coupled with the solver generating the data set. Thereby an adaptive sampling frequency can be determined on the y by evaluating the complexity of the detected events. 130 Bibliography Abdel-Kareem, W. 2011 Tracking of vortical structures in three-dimensional decaying homogeneous isotropic turbulence. Int J Mod Phys C 22 (12), 1373{1391. Acarlar, M. S. & Smith, C. R. 1987 A study of hairpin vortices in a laminar boundary layer. Part 1. hairpin vortices generated by a hemisphere protuberance. J. Fluid Mech. 175, 1{41. Adrian, R. J. 2007 Hairpin vortex organization in wall turbulence. Phys. Fluids 19 (4), 041301. Adrian, R. J., Meinhart, C. D. & Tomkins, C. D. 2000 Vortex organization in the outer region of the turbulent boundary layer. J. Fluid Mech. 422, 1{54. Agui, J. H., Briassulis, G. & Andreopoulos, Y. 2005 Studies of interactions of a propagating shock wave with decaying grid turbulence: velocity and vorticity elds. J. Fluid Mech. 524, 143. Andreopoulos, Y., Agui, J. H. & Briassulis, G. 2000 Shock wave|turbulence interactions. Ann. Rev. Fluid Mech. 32 (1), 309{345. Ashurst, W. T., Kerstein, A., Kerr, R. & Gibson, C. 1987 Alignment of vorticity and scalar gradient with strain rate in simulated navier{stokes turbulence. Phys. Fluids 30 (8), 2343{2353. Attili, A. & Bisetti, F. 2012 Statistics and scaling of turbulence in a spatially developing mixing layer at re= 250. Phys. Fluids 24 (3), 035109. Bajaj, C., Shamir, A. & Sohn, B.-S. 2002 Progressive tracking of isosurfaces in time-varying scalar elds . Banks, D. C. & Singer, B. A. 1995 A predictor-corrector technique for visualizing unsteady ow. IEEE Trans Vis Comput Graph 1 (2), 151{163. Barre, S., Alem, D. & Bonnet, J. 1996 Experimental study of a normal shock/homogeneous turbulence interaction. AIAA J. 34 (5), 968{974. Batchelor, G. K. 1959 Small-scale variation of convected quantities like temperature in turbulent uid part 1. general discussion and the case of small conductivity. J. Fluid Mech. 5 (1), 113{133. 131 Bauer, D. & Peikert, R. 2002 Vortex tracking in scale-space. In Proceedings of the symposium on Data Visualisation 2002 , pp. 233{. Eurographics Association. Bell, J. H. & Mehta, R. D. 1990 Development of a two-stream mixing layer from tripped and untripped boundary layers. AIAA J. 28 (12), 2034{2042. Bellman, R. E. 2015 Adaptive control processes: a guided tour. Princeton university press. Bentley, J. L. 1975 Multidimensional binary search trees used for associative searching. Commun. ACM 18 (9), 509{517. Berdahl, C. & Thompson, D. 1993 Eduction of swirling structure using the velocity gradient tensor. AIAA J. 31 (1), 97{103. Bergen, G. v. d. 1999 A fast and robust gjk implementation for collision detection of convex objects. J. graph. tools 4 (2), 7{25. Bermejo-Moreno, I. & Pullin, D. 2008 On the non-local geometry of turbulence. J. Fluid Mech. 603, 101{135. Bilger, R. W. 1976 The structure of diusion ames. Combust. Sci. Technol. 13 (1-6), 155{170. Bilger, R. W. 2004 Some aspects of scalar dissipation. Flow, Turb. and Combust. 72 (2), 93{114. Birch, S. F. & Eggers, J. 1972 A critical reviewof the experimental data for developedfree turbulent shear layers. Free Turbulent Shear Flows, VolumeI-ConferenceProceedings, NASASP-321 pp. 11{40. Blanco, J. L. & Rai, P. K. 2014 nano ann: a C++ header-only fork of FLANN, a library for nearest neighbor (NN) with kd-trees. https://github.com/jlblancoc/nanoflann. Bogdanoff, D. 1983 Compressibility eects in turbulent shear layers. AIAA J. 21 (6), 926{927. Bonnet, J., Debisschop, J. & Chambres, O. 1993 Experimental studies of the turbulent structure of supersonic mixinglayers. In 31st Aerospace Sciences Meeting, p. 217. Boukharfane, R., Bouali, Z. & Mura, A. 2018 Evolution of scalar and velocity dynamics in planar shock-turbulence interaction. Shock Waves 28 (6), 1117{1141. Brown, G. L. & Roshko, A. 1974 On density eects and large structure in turbulent mixing layers. J. Fluid Mech. 64 (4), 775{816. Budzinski, J., Zukoski, E. & Marble, F. 1992 Rayleigh scattering measurements of shock enhanced mixing. In 28th Joint Propulsion Conference and Exhibit, p. 3546. Caban, J., Joshi, A. & Rheingans, P. 2007 Texture-based feature tracking for eective time-varying data visualization. IEEE Trans Vis Comput Graph 13 (6), 1472{1479. 132 Cambon, C., Coleman, G. & Mansour, N. 1992 Rapid distortion analysis and direct simulation of compressible homogeneous turbulence at nite mach number . Chac n, J. M., Cantwell, B. J. & Kline, S. J. 1996 Study of turbulent boundary layer structure using the invariants of the velocity gradient tensor. Exp. Therm. Fluid Sci. 13 (4), 308{317. Chan, W., Dodd, M., Johnson, P., Urzay, J. & Moin, P. 2018 Formation and dynamics of bubbles in breaking waves: Part i. algorithms for the identication of bubbles and breakup/coalescence events . Chen, J., Silver, D. & Jiang, L. 2003 The feature tree: Visualizing feature tracking in distributed amr datasets. In IEEE Symposium on Parallel and Large-Data Visualization and Graphics, 2003. PVG 2003., pp. 103{110. IEEE. Chong, M. S., Perry, A. E. & Cantwell, B. J. 1990 A general classication of three-dimensional ow elds. Phys Fluid Fluid Dynam 2 (5), 765{777. Clemens, N. & Mungal, M. 1995 Large-scale structure and entrainment in the supersonic mixing layer. J. Fluid Mech. 284, 171{216. Clyne, J., Mininni, P. & Norton, A. 2012 Physically-based feature tracking for cfd data. IEEE Trans Vis Comput Graph 19 (6), 1020{1033. Cook, D. J. & Holder, L. B. 1993 Substructure discovery using minimum description length and back- ground knowledge. J Artif Intell Res 1, 231{255. Cook, D. J. & Holder, L. B. 2006 Mining graph data. John Wiley & Sons. Dallmann, U. 1983 Topological structures of three-dimensional vortex ow separation. In 16th Fluid and Plasmadynamics Conference, p. 1735. Danish, M., Suman, S. & Girimaji, S. S. 2016 In uence of ow topology and dilatation on scalar mixing in compressible turbulence. J. Fluid Mech. 793, 633{655. Del Alamo, J. C., Jimenez, J., Zandonade, P. & Moser, R. D. 2006 Self-similar vortex clusters in the turbulent logarithmic region. J. Fluid Mech. 561, 329{358. Demanet, L. & Ying, L. 2007 Curvelets and wave atoms for mirror-extended images. In Wavelets XII , , vol. 6701, p. 67010J. International Society for Optics and Photonics. Dong, S., Lozano-Dur an, A., Sekimoto, A. & Jim enez, J. 2017 Coherent structures in statistically stationary homogeneous shear turbulence. J. Fluid Mech. 816, 167{208. Donzis, D. A., Sreenivasan, K. & Yeung, P. 2010 The batchelor spectrum for mixing of passive scalars in isotropic turbulence. Flow Turb. Combust. 85 (3), 549{566. 133 Ducros, F., Laporte, F., Soul eres, T., Guinot, V., Moinat, P. & Caruelle, B. 2000 High-order uxes for conservative skew-symmetric-like schemes in structured meshes: application to compressible ows. J. Comp. Phys. 161 (1), 114{139. Eitel-Amor, G., Orl u, R., Schlatter, P. & Flores, O. 2015 Hairpin vortices in turbulent boundary layers. Phys. Fluids 27 (2), 025108. Epps, B. 2017 Review of vortex identication methods. In 55th AIAA aerospace sciences meeting, p. 0989. Ester, M., Kriegel, H.-P., Sander, J., Xu, X. & others 1996 A density-based algorithm for discov- ering clusters in large spatial databases with noise. In kdd, , vol. 96, pp. 226{231. Etemadi, N. 1990 On curve and surface stretching in isotropic turbulent ow. J. Fluid Mech. 221, 685{692. Foysi, H. & Sarkar, S. 2010 The compressible mixing layer: an les study. Theor Comput Fluid Dyn 24 (6), 565{588. Freund, J. B., Lele, S. K. & Moin, P. 2000 Compressibility eects in a turbulent annular mixing layer. part 1. turbulence and growth rate. J. Fluid Mech. 421, 229{267. Friedman, J. H., Bentley, J. L. & Finkel, R. A. 1977 An algorithm for nding best matches in logarithmic expected time. ACM Trans Math Softw (TOMS) 3 (3), 209{226. Friedrich, R. & Bertolotti, F. 1996 Compressibility eects due to turbulent uctuations. Applied scientic research 57 (3-4), 165{194. Fu, S. & Li, Q. 2006 Numerical simulation of compressible mixing layers. Int. J. Heat Fluid Flow 27 (5), 895{901. Gao, X., Bermejo-Moreno, I. & Larsson, J. 2020 Parametric numerical study of passive scalar mixing in shock turbulence interaction. J. Fluid Mech. 895. Gilbert, E. G., Johnson, D. W. & Keerthi, S. S. 1988 A fast procedure for computing the distance between complex objects in three-dimensional space. IEEE Robot Autom Mag 4 (2), 193{203. Girimaji, S. & Pope, S. 1990 Material-element deformation in isotropic turbulence. J. Fluid Mech. 220, 427{458. Girimaji, S. & Pope, S. 1992 Propagating surfaces in isotropic turbulence. J. Fluid Mech. 234, 247{277. Girvan, M. & Newman, M. E. 2002 Community structure in social and biological networks. Proc. Natl. Acad. Sci. U.S.A 99 (12), 7821{7826. Goto, S. & Kida, S. 2007 Reynolds-number dependence of line and surface stretching in turbulence: folding eects. J. Fluid Mech. 586, 59. 134 Gottschalk, S., Lin, M. C. & Manocha, D. 1996 Obbtree: A hierarchical structure for rapid interference detection. In Proceedings of the 23rd annual conference on Computer graphics and interactive techniques, pp. 171{180. Guo, H., Phillips, C. L., Peterka, T., Karpeyev, D. & Glatz, A. 2015 Extracting, tracking, and visualizing magnetic ux vortices in 3d complex-valued superconductor simulation data. IEEE Trans Vis Comput Graph 22 (1), 827{836. Harikrishnan, A., Ansorge, C., Klein, R. & Vercauteren, N. 2021 Geometry and organization of coherent structures in stably stratied atmospheric boundary layers. arXiv preprint arXiv:2110.02253 . Hesselink, L. & Sturtevant, B. 1988 Propagation of weak shocks through a random medium. J. Fluid Mech. 196, 513{553. Holder, L. B., Cook, D. J., Djoko, S. & others 1994 Substucture discovery in the subdue system. In KDD workshop, pp. 169{180. Washington, DC, USA. Hunt, J., Wray, A. & Moin, P. 1988 Eddies, stream, and convergence zones in turbulent ows. Proceed- ings of the Summer Program (Center for Turbulence Research) pp. 193{208. Hussain, A. F. 1986 Coherent structures and turbulence. J. Fluid Mech. 173, 303{356. Jeong, J. & Hussain, F. 1995 On the identication of a vortex. J. Fluid Mech. 285 (-1), 69. Ji, G. & Shen, H.-W. 2004 Ecient isosurface tracking using precomputed correspondence table. In Proceedings of the Sixth Joint Eurographics-IEEE TCVG conference on Visualization, pp. 283{292. Euro- graphics Association. Ji, G. & Shen, H.-W. 2006 Feature tracking using earth mover's distance and global optimization. In Pacic graphics, , vol. 2. Ji, G., Shen, H.-W. & Wenger, R. 2003 Volume tracking using higher dimensional isosurfacing. In IEEE Visualization, 2003. VIS 2003., pp. 209{216. IEEE. Jim enez, J. 2013 Near-wall turbulence. Phys. Fluids 25 (10), 101302. John, J. P., Donzis, D. A. & Sreenivasan, K. R. 2019 Solenoidal scaling laws for compressible mixing. Phys. Rev. Lett. 123 (22), 224501. Karasso, P. & Mungal, M. 1996 Scalar mixing and reaction in plane liquid shear layers. J. Fluid Mech. 323, 23{63. Kareem, W. A., Izawa, S., Xiong, A. & Fukunishi, Y. 2007 Extraction and tracking of multi-scaled vortices from a homogeneous isotropic turbulent eld. J. Turb. (8), N12. 135 Kerr, R. M. 1985 Higher-order derivative correlations and the alignment of small-scale structures in isotropic numerical turbulence. J. Fluid Mech. 153, 31{58. Ketkar, N. S., Holder, L. B. & Cook, D. J. 2005 Subdue: Compression-based frequent pattern discovery in graph data. In Proceedings of the 1st international workshop on open source data mining: frequent pattern mining implementations, pp. 71{76. Kim, S. H. & Pitsch, H. 2007 Scalar gradient and small-scale structure in turbulent premixed combustion. Phys. Fluids 19 (11), 115104. Kline, S. J., Reynolds, W. C., Schraub, F. A. & Runstadler, P. W. 1967 The structure of turbulent boundary layers. J. Fluid Mech. 30 (4), 741{773. Koenderink, J. J. & Van Doorn, A. J. 1992 Surface shape and curvature scales. Image Vis Comput 10 (8), 557{564. Larsson, J. 2009 Blending technique for compressible in ow turbulence: algorithm localization and accu- racy assessment. J. Comput. Phys. 228 (4), 933{937. Larsson, J., Bermejo-Moreno, I. & Lele, S. K. 2013 Reynolds-and mach-number eects in canonical shock-turbulence interaction. J. Fluid Mech. 717, 293. Larsson, J. & Lele, S. K. 2009 Direct numerical simulation of canonical shock/turbulence interaction. Phys. Fluids 21 (12), 126101. Lee, S., Lele, S. K. & Moin, P. 1993 Direct numerical simulation of isotropic turbulence interacting with a weak shock wave. J. Fluid Mech. 251, 533{562. Lee, S., Lele, S. K. & Moin, P. 1994 Corrigendum: direct numerical simulation of isotropic turbulence interacting with a weak shock wave. J. Fluid Mech. 264, 373{374. Lee, S., Lele, S. K. & Moin, P. 1997 Interaction of isotropic turbulence with shock waves: eect of shock strength. J. Fluid Mech. 340, 225{247. Lesieur, M., M etais, O., Comte, P. & others 2005 Large-eddy simulations of turbulence. Cambridge university press. Leung, T., Swaminathan, N. & Davidson, P. 2012 Geometry and interaction of structures in homoge- neous isotropic turbulence. J. Fluid Mech. 710, 453{481. Li, Q. & Fu, S. 2003 Numerical simulation of high-speed planar mixing layer. Comput. Fluids 32 (10), 1357{1377. Liu, T., Moore, A. W., Yang, K. & Gray, A. G. 2005 An investigation of practical approximate nearest neighbor algorithms. In Advances in neural information processing systems, pp. 825{832. 136 Livescu, D. & Ryu, J. 2016 Vorticity dynamics after the shock{turbulence interaction. Shock Waves 26 (3), 241{251. Lorensen, W. E. & Cline, H. E. 1987 Marching cubes: A high resolution 3d surface construction algorithm. ACM siggraph comput. graph. 21 (4), 163{169. Lozano-Dur an, A. & Jim enez, J. 2014 Time-resolved evolution of coherent structures in turbulent chan- nels: characterization of eddies and cascades. J. Fluid Mech. 759, 432{471. Lundgren, T. S. 1982 Strained spiral vortex model for turbulent ne structure. Phys. Fluids 25 (12), 2193{2203. Mani, A. 2012 Analysis and optimization of numerical sponge layers as a nonre ective boundary treatment. J. Comput. Phys. 231 (2), 704{716. Maurus, S. & Plant, C. 2016 Skinny-dip: clustering in a sea of noise. In Proceedings of the 22nd ACM SIGKDD international conference on Knowledge discovery and data mining, pp. 1055{1064. Mehta, R. 1991 Eect of velocity ratio on plane mixing layer development: In uence of the splitter plate wake. Exp. Fluids. 10 (4), 194{204. Misra, A. & Pullin, D. I. 1997 A vortex-based subgrid stress model for large-eddy simulation. Phys. Fluids 9 (8), 2443{2454. Montanari, M., Petrinic, N. & Barbieri, E. 2017 Improving the gjk algorithm for faster and more reliable distance queries between convex objects. ACM Trans. Graph. (TOG) 36 (3), 1{17. Moser, R. & Rogers, M. 1991 Mixing transition and the cascade to small scales in a plane mixing layer. Phys Fluid Fluid Dynam 3 (5), 1128{1134. Moser, R. D. & Rogers, M. M. 1993 The three-dimensional evolution of a plane mixing layer: pairing and transition to turbulence. J. Fluid Mech. 247, 275{320. Muelder, C. & Ma, K.-L. 2009 Interactive feature extraction and tracking by utilizing region coherency. In 2009 IEEE Pacic Visualization Symposium, pp. 17{24. IEEE. Ni, Q. 2015 Compressible turbulent mixing: Eects of schmidt number. Phys. Rev. E 91 (5), 053020. Ni, Q. 2016 Compressible turbulent mixing: Eects of compressibility. Phys. Rev. E 93 (4), 043116. Pan, L. & Scannapieco, E. 2010 Mixing in supersonic turbulence. Astrophys. J. 721 (2), 1765. Pan, L. & Scannapieco, E. 2011 Passive scalar structures in supersonic turbulence. Phys. Rev. E 83 (4), 045302. 137 Panickacheril John, J., Donzis, D. A. & Sreenivasan, K. R. 2020 Compressibility eects on the scalar dissipation rate. Combust. Sci. Technol. 192 (7), 1320{1333. Pantano, C. & Sarkar, S. 2002 A study of compressibility eects in the high-speed turbulent shear layer using direct simulation. J. Fluid Mech. 451, 329. Papamoschou, D. & Lele, S. K. 1993 Vortex-induced disturbance eld in a compressible shear layer. Phys Fluid Fluid Dynam 5 (6), 1412{1419. Parashar, N., Sinha, S. S., Danish, M. & Srinivasan, B. 2017 Lagrangian investigations of vorticity dynamics in compressible turbulence. Phys. Fluids 29 (10), 105110. Peikert, R. & Roth, M. 1999 The" parallel vectors" operator-a vector eld visualization primitive. In Proceedings Visualization'99 (Cat. No. 99CB37067), pp. 263{532. Ieee. Pitsch, H. & Steiner, H. 2000 Scalar mixing and dissipation rate in large-eddy simulations of non- premixed turbulent combustion. Proc. Combust. Inst. 28 (1), 41{49. Post, F. H., Vrolijk, B., Hauser, H., Laramee, R. S. & Doleisch, H. 2002 Feature extraction and visualisation of ow elds. In Eurographics (STARs). Pullin, D. I. 2000 A vortex-based model for the subgrid ux of a passive scalar. Phys. Fluids 12 (9), 2311{2319. Pullin, D. I. & Saffman, P. G. 1993 On the Lundgren{Townsend model of turbulent ne scales. Phys. Fluids A: Fluid Dyn. 5 (1), 126{145. Reinders, F., Post, F. H. & Spoelder, H. J. 1999 Attribute-based feature tracking. In Data visualiza- tion'99, pp. 63{72. Springer. Reinders, K. F., Frederik, K. & Reinders, J. 2001 Feature-based visualization of time-dependent data . Rissanen, J. 1998 Stochastic complexity in statistical inquiry, , vol. 15. World scientic. Ristorcelli, J. R. & Blaisdell, G. A. 1997 Consistent initial conditions for the dns of compressible turbulence. Phys. Fluids 9 (1), 4{6. Robinson, S. K. 1991 Coherent motions in the turbulent boundary layer. Annu. Rev. Fluid Mech. 23 (1), 601{639. Rogers, M. M. & Moin, P. 1987 The structure of the vorticity eld in homogeneous turbulent ows. J. Fluid Mech 176, 33{66. 138 Rogers, M. M. & Moser, R. D. 1992 The three-dimensional evolution of a plane mixing layer: the kelvin{helmholtz rollup. J. Fluid Mech. 243, 183{226. Rogers, M. M. & Moser, R. D. 1994 Direct simulation of a self-similar turbulent mixing layer. Phys. Fluids 6 (2), 903{923. Ryu, J. & Livescu, D. 2014 Turbulence structure behind the shock in canonical shock{vortical turbulence interaction. J. Fluid Mech. 756. Saikia, H. & Weinkauf, T. 2017 Global feature tracking and similarity estimation in time-dependent scalar elds. In Computer Graphics Forum, , vol. 36, pp. 1{11. Wiley Online Library. Samtaney, R., Silver, D., Zabusky, N. & Cao, J. 1994 Visualizing features and tracking their evolution. Computer 27 (7), 20{27. Sandham, N. & Reynolds, W. 1990 Compressible mixing layer-linear theory and direct simulation. AIAA J. 28 (4), 618{624. Sarkar, S. 1995 The stabilizing eect of compressibility in turbulent shear ow. J. Fluid Mech. 282, 163{186. Sauer, F., Yu, H. & Ma, K.-L. 2014 Trajectory-based ow feature tracking in joint particle/volume datasets. IEEE Trans Vis Comput Graph 20 (12), 2565{2574. Schafhitzel, T., Baysal, K., Rist, U., Weiskopf, D. & Ertl, T. 2008 Particle-based vortex core line tracking taking into account vortex dynamics. In Proceedings International Symposium on Flow Visual- ization, , vol. 8. Schlatter, P., Li, Q., Orl u, R., Hussain, F. & Henningson, D. S. 2014 On the near-wall vortical structures at moderate Reynolds numbers. Eur. J. Mech.-B/Fluids 48, 75{93. Schubert, E., Sander, J., Ester, M., Kriegel, H. P. & Xu, X. 2017 Dbscan revisited, revisited: why and how you should (still) use dbscan. ACM Trans. Database Syst. (TODS) 42 (3), 1{21. Silver, D. 1995 Object-oriented visualization. IEEE Comput Graph Appl 15 (3), 54{62. Silver, D. & Wang, X. 1996 Volume tracking. In Proceedings of Seventh Annual IEEE Visualization'96 , pp. 157{164. IEEE. Silver, D. & Wang, X. 1997 Tracking and visualizing turbulent 3d features. IEEE Trans Vis Comput Graph 3 (2), 129{141. Silver, D. & Wang, X. 1998 Tracking scalar features in unstructured data sets. In Proceedings Visualiza- tion'98 (Cat. No. 98CB36276), pp. 79{86. IEEE. 139 Smith, C. R. & Metzler, S. P. 1983 The characteristics of low-speed streaks in the near-wall region of a turbulent boundary layer. J. Fluid Mech. 129, 27{54. Sohn, B.-S. & Bajaj, C. 2005 Time-varying contour topology. IEEE Trans Vis Comput Graph 12 (1), 14{25. Soler, M., Plainchault, M., Conche, B. & Tierny, J. 2018 Lifted wasserstein matcher for fast and robust topology tracking. arXiv preprint arXiv:1808.05870 . Suman, S. & Girimaji, S. S. 2010 Velocity gradient invariants and local ow-eld topology in compressible turbulence. J. Turbul. (11), N2. Swaminathan, N. & Grout, R. 2006 Interaction of turbulence and scalar elds in premixed ames. Phys. Fluids 18 (4), 045102. Tennekes, H. 1968 Simple model for the small-scale structure of turbulence. Phys. Fluids 11 (3), 669{671. Theisel, H., Sahner, J., Weinkauf, T., Hege, H.-C. & Seidel, H.-P. 2005 Extraction of parallel vector surfaces in 3d time-dependent elds and application to vortex core line tracking. In VIS 05. IEEE Visualization, 2005., pp. 631{638. IEEE. Theisel, H. & Seidel, H.-P. 2003 Feature ow elds. In Proceedings of the symposium on Data visualisa- tion 2003, pp. 141{148. Eurographics Association. Theodorsen, T. 1952 Mechanisms of turbulence. In Proceedings of the 2^< nd> Midwestern Conference on Fluid Mechanics, 1952 . Theodorsen, T. 1955 The structure of turbulence. In 50 Jahre Grenzschichtforschung, pp. 55{62. Springer. Tian, Y., Jaberi, F. A., Li, Z. & Livescu, D. 2017 Numerical study of variable density turbulence interaction with a normal shock wave. J. Fluid Mech. 829, 551{588. Tomkins, C. D. & Adrian, R. J. 2003 Spanwise structure and scale growth in turbulent boundary layers. J. Fluid Mech. 490, 37{74. Townsend, A. 1980 The structure of turbulent shear ow. Cambridge university press. Townsend, A. A. 1951 On the ne-scale structure of turbulence. Proc. Math. Phys. Sci. 208 (1095), 534{542. Townsend, A. A. 1961 Equilibrium layers and wall turbulence. J. Fluid Mech. 11 (1), 97{120. Urzay, J. 2018 Supersonic combustion in air-breathing propulsion systems for hypersonic ight. Annu. Rev. Fluid Mech. 50, 593{627. 140 Van Den Bergen, G. 2001 Proximity queries and penetration depth computation on 3d game objects. In Game developers conference, , vol. 170. Van Walsum, T., Post, F. H., Silver, D. & Post, F. J. 1996 Feature extraction and iconic visualization. IEEE Trans Vis Comput Graph 2 (2), 111{119. Vanderwel, C. & Tavoularis, S. 2011 Coherent structures in uniformly sheared turbulent ow. J. Fluid Mech. 689, 434{464. Vedula, P., Yeung, P. & Fox, R. O. 2001 Dynamics of scalar dissipation in isotropic turbulence: a numerical and modelling study. J. Fluid Mech. 433, 29. Villermaux, E. 2019 Mixing versus stirring. Ann. Rev. Fluid Mech. 51, 245{273. Vollmers, H., Kreplin, H. & Meier, H. 1983 Separation and vortical-type ow around a prolate spheroid-evaluation of relevant parameters. Tech. Rep.. DEUTSCHE FORSCHUNGS-UND VERSUCH- SANSTALT FUER LUFT-UND RAUMFAHRT EV . . . . Vreman, B., Kuerten, H. & Geurts, B. 1995 Shocks in direct numerical simulation of the conned three-dimensional mixing layer. Phys. Fluids 7 (9), 2105{2107. Wadell, H. 1935 Volume, shape, and roundness of quartz particles. J. Geol. 43 (3), 250{280. Weigle, C. & Banks, D. C. 1998 Extracting iso-valued features in 4-dimensional scalar elds. In IEEE Symposium on Volume Visualization (Cat. No. 989EX300), pp. 103{110. IEEE. Widanagamaachchi, W., Chen, J., Klacansky, P., Pascucci, V., Kolla, H., Bhagatwala, A. & Bremer, P.-T. 2015 Tracking features in embedded surfaces: Understanding extinction in turbulent combustion. In 2015 IEEE 5th Symposium on Large Data Analysis and Visualization (LDAV), pp. 9{16. IEEE. Widanagamaachchi, W., Christensen, C., Pascucci, V. & Bremer, P.-T. 2012 Interactive explo- ration of large-scale time-varying data using dynamic tracking graphs. In IEEE symposium on large data analysis and visualization (LDAV), pp. 9{17. IEEE. Wu, X. & Moin, P. 2009 Direct numerical simulation of turbulence in a nominally zero-pressure-gradient at-plate boundary layer. J. Fluid Mech. 630, 5{41. Yan, X. & Han, J. 2002 gspan: Graph-based substructure pattern mining. In 2002 IEEE International Conference on Data Mining, 2002. Proceedings., pp. 721{724. IEEE. Yang, Y., Pullin, D. & Bermejo-Moreno, I. 2010 Multi-scale geometric analysis of lagrangian struc- tures in isotropic turbulence. J. Fluid Mech. 654, 233{270. 141 Yianilos, P. N. 1993 Data structures and algorithms for nearest neighbor search in general metric spaces. Zhang, D., Tan, J. & Yao, X. 2019 Direct numerical simulation of spatially developing highly compressible mixing layer: Structural evolution and turbulent statistics. Phys. Fluids 31 (3), 036102. Zhang, Y., Liu, K., Xian, H. & Du, X. 2018 A review of methods for vortex identication in hydrotur- bines. Renew. Sust. Energ. Rev. 81, 1269{1285. Zhou, J., Adrian, R. J., Balachandar, S. & Kendall, T. 1999 Mechanisms for generating coherent packets of hairpin vortices in channel ow. J. Fluid Mech. 387, 353{396. Zhou, Q., He, F. & Shen, M. 2012 Direct numerical simulation of a spatially developing compressible plane mixing layer: ow structures and mean ow properties. J. Fluid Mech. 711, 437. 142

Abstract (if available)

Linked assets

University of Southern California Dissertations and Theses

Conceptually similar

PDF

Direct numerical simulation of mixing and combustion under canonical shock turbulence interaction

PDF

Numerical study of shock-wave/turbulent boundary layer interactions on flexible and rigid panels with wall-modeled large-eddy simulations

PDF

Supersonic flow study of a capsule geometry using large-eddy simulation

PDF

A Boltzmann model for tracking aerosols in turbulent flows

PDF

Modeling and analysis of parallel and spatially-evolving wall-bounded shear flows

PDF

Hydrodynamic stability of two fluid columns of different densities and viscosities subject to gravity

PDF

Large eddy simulations of turbulent flows without use of the eddy viscosity concept

PDF

The projection immersed boundary method for compressible flow and its application in fluid-structure interaction simulations of parachute systems

PDF

Model based design of porous and patterned surfaces for passive turbulence control

PDF

An experimental study of shock wave attenuation

PDF

Numerical study of flow characteristics of controlled vortex induced vibrations in cylinders

PDF

Development and applications of a body-force propulsor model for high-fidelity CFD

PDF

Pattern generation in stratified wakes

PDF

Re-assessing local structures of turbulent flames via vortex-flame interaction

PDF

On the simulation of stratified turbulent flows

PDF

Aerodynamics at low Re: separation, reattachment, and control

PDF

Numerical study of focusing effects generated by the coalescence of multiple shock waves

PDF

Tsunami-induced turbulent coherent structures

PDF

Large eddy simulations of laminar separation bubble flows

PDF

Boundary layer and separation control on wings at low Reynolds numbers

Asset Metadata

Creator Buchmeier, Jonas (author)

Core Title Tracking and evolution of compressible turbulent flow structures

School Viterbi School of Engineering

Degree Doctor of Philosophy

Degree Program Aerospace Engineering

Degree Conferral Date 2022-05

Publication Date 02/08/2022

Defense Date 01/24/2022

Publisher University of Southern California (original), University of Southern California. Libraries (digital)

Tag flow feature tracking,fluid dynamics,geometry of turbulence,mixing layers,OAI-PMH Harvest,shock compression,turbulent scalar mixing,vortex dynamics

Format application/pdf (imt)

Language English

Contributor Electronically uploaded by the author (provenance)

Advisor Bermejo-Moreno, Ivan (committee chair), Domaradzki, Julian (committee member), Nakano, Aiichiro (committee member), Pantano-Rubino, Carlos (committee member)

Creator Email jbuchmei@usc.edu,jonas.buchmeier@aol.de

Permanent Link (DOI) https://doi.org/10.25549/usctheses-oUC110702279

Unique identifier UC110702279

Legacy Identifier etd-BuchmeierJ-10383

Document Type Dissertation

Format application/pdf (imt)

Rights Buchmeier, Jonas

Type texts

Source 20220214-usctheses-batch-912 (batch), University of Southern California (contributing entity), University of Southern California Dissertations and Theses (collection)

Access Conditions The author retains rights to his/her dissertation, thesis or other graduate work according to U.S. copyright law. Electronic access is being provided by the USC Libraries in agreement with the author, as the original true and official version of the work, but does not grant the reader permission to use the work if the desired use is covered by copyright. It is the author, as rights holder, who must provide use permission if such use is covered by copyright. The original signature page accompanying the original submission of the work to the USC Libraries is retained by the USC Libraries and a copy of it may be obtained by authorized requesters contacting the repository e-mail address given.

Repository Name University of Southern California Digital Library

Repository Location USC Digital Library, University of Southern California, University Park Campus MC 2810, 3434 South Grand Avenue, 2nd Floor, Los Angeles, California 90089-2810, USA

Repository Email cisadmin@lib.usc.edu