Geometry and dynamics of passive scalar structures in compressible turbulent mixing

A comprehensive understanding of compressibility effects on turbulent mixing processes is of great relevance to advance air-breathing super- and hypersonic propulsion systems.^1,2 Shock wave-induced mixing enhancement has been proposed to accommodate the very short timescale requirements for nonpremixed combustion in these applications.³ Furthermore, compressibility effects can originate from heat release in non-premixed and, especially, in premixed combustion.⁴ Despite the identification of the scalar dissipation rate as a key parameter in non-premixed combustion several decades ago,^5,6 the study of compressibilty effects on scalar mixing has been the focus only in more recent direct numerical simulations (DNS) of homogeneous isotropic turbulence (HIT)^7–12 and also shock-turbulence interaction (STI).^13–15 These numerical studies confirm that the amplification of turbulence kinetic energy (TKE) and transverse vorticity variance in conjunction with decreased turbulence length scales^16–18 result in an increase in the scalar dissipation rate across a shock wave, as suggested by earlier experiments.^19,20

While the TKE amplification saturates²¹ for Mach number, M, larger than 3, the amplification 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.¹⁵ The shock-induced enhancement of scalar mixing is expressed by a reduced scalar Taylor-microscale and an intensified decay of scalar variance in the post-shock region.¹⁴ The interaction of spatial velocity and scalar gradients contributes to the scalar transport.¹¹ The increase in scalar dissipation across the shock is associated with the evolution of the alignment between the scalar gradient, $∇ϕ$⁠, the vorticity, $ω$⁠, and the strain-rate tensor eigenvectors, denoted as $α,β$⁠, and $γ$⁠,^22–24 and ordered by their respective eigenvalues from most extensional to most compressive, $α≥β≥γ$⁠. Recent statistical volumetric analyses of HIT^11,23 and STI^14,15 have shown that the scalar gradient is most aligned with the most compressive eigendirection $γ$ of the strain-rate tensor, followed by $α$ and $β$⁠,^25–27 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.²⁸ 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.^29,30 In contrast, the alignment with the intermediate and most compressive 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 flow regions of high dilatation was recently reported.¹² A strong anti-correlation found between the scalar gradient and the vorticity²² has been interpreted structurally as the passive scalar being wound around vortex tubes. The decrease in scalar gradient alignment with the intermediate eigenvector $β$ and the vorticity $ω$ is thus expected since the alignment between $β$ and $ω$ is strengthened by the shock-interaction.^14,31 The alignment analysis has been utilized in the development of SGS models for the flux of a passive scalar.³² 

In the aforementioned studies, the nature of compressible turbulent mixing is investigated by means of statistical analyses of the volumetric scalar field in HIT or in terms of the streamwise evolution of transverse plane averages in STI. Complementing statistical approaches, structural analyses have provided fundamental insight into the mechanisms of momentum and energy transfer in turbulence. For example, in the buffer layer of near-wall turbulence,^33,34 streaks of streamwise velocity are found to carry most of the TKE, and quasi-streamwise vortices organize dissipation and Reynolds stresses,³⁵ whereas in the logarithmic layer, most of the momentum transfer originates from wall-attached eddies³⁶ and vortices organize into self-similar clusters.³⁷ Hairpin vortices are also frequently identified in wall bounded flows^38–43 as well as in homogeneous shear turbulence^44–46 and turbulent mixing layers,^47–51 where rollers and braids are also found to characterize coherent motions. Structural analyses of turbulent flows have led to the development of structure-based models of turbulence fine scales^52–55 and subgrid-scale models for large eddy simulations.⁵⁶ 

To complement statistical studies of the scalar mixing, the research objective of the present work 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.

Time-resolved numerical^57–59 and experimental^60,61 studies have previously been conducted to analyze the dynamics of coherent flow structures. Recording the evolution of individual finite-sized scalar structures enables to investigate how the scalar mixing is represented at the structural level in the form of the breakup processes and the lifetime of structures until complete diffusion. 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 structures and hence are also examined from a structural viewpoint. For example, it is of interest how the distribution of the alignment between the scalar gradient $∇ϕ$ with the different strain-rate eigenvectors $α,β$⁠, and $γ$ on the isosurface correlates with the local geometry of the surface. Relating evolution of the geometry of the passive scalar structures to the underlying flow physics reveals how the shock-induced structure deformation amplifies 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.

To study the geometrical properties of turbulence structures different description techniques have been suggested in the literature. For example, geometrical descriptors based on Minkowski functionals have been used to characterize enstrophy structures in DHIT.⁶² In this study, we apply the geometrical characterization method introduced by Bermejo-Moreno and Pullin.⁶³ Whereas the aforementioned techniques describe turbulence structures for instantaneous snapshots of a given dataset, the present work focuses on the evolution of those geometrical properties over time and the relation to the temporally or spatially evolving physics of flow field. Therefore, a novel method to track individual flow structures throughout the flow field is introduced in Sec. III adding temporal coherence to the geometrical analysis. The method enables simultaneous tracking of intertwined structures with significant concave regions over large sampling intervals by combining region- and attribute-based algorithms. The tracking attributes incorporate geometric properties of the structures obtained from differential geometry analyses, which are then related to the physical quantities relevant to scalar mixing.

The methods presented in this work are applicable to any kind of flow features defined as closed surfaces embedded in a three-dimensional space. The surfaces could be derived from a numerical simulation as well as originate from experimental data. In this work, the surfaces (structures) are extracted from the volumetric passive scalar fields of direct numerical simulations by isosurfacing. Flow feature extraction methods themselves are extrinsic to the presented methods and not of the subject of this work.

The paper is organized as follows. We describe in Sec. II the conducted DNS of passive scalar mixing in HIT and STI with details on the well-defined initial geometry of the scalar structures. In Sec. III, we summarize the properties used for the structure differential-geometrical analysis and introduce the flow feature tracking method. The results of the temporal geometrical passive scalar structure analysis and their relation to the underlying flow physics are presented in Sec. IV.

To elucidate the effect of shock waves on scalar mixing from a structural viewpoint, we compare temporally decaying HIT and the canonical statistically stationary interaction of a nominally planar shock wave with spatially decaying isotropic turbulence (STI) at two different Mach numbers. We perform shock-capturing DNS following the methodology described in Gao et al.¹⁵ Based on the earlier work of Larsson et al.,^64,65 the numerical methodology combines a fifth-order weighted essentially non-oscillatory (WENO) scheme with HLL flux splitting in regions near discontinuities (shocks or large scalar gradients) and a sixth-order central difference scheme elsewhere. The compressible Navier–Stokes equations for a calorically perfect gas are complemented with additional advection-diffusion equations for the passive scalars,

Here, ρ is the density, $ϕm$ is the mth passive scalar, u_j the jth Cartesian velocity component (j = 1, 2, 3), and D_m the scalar diffusivities. The turbulence is characterized by the Taylor microscale Reynolds number, $Reλ=ρ¯u′rmsλ/μ¯$⁠, and the turbulence Mach number, $Mt=3u′rms/c̃$⁠, where μ is the dynamic viscosity, $λ=[u′2u′2¯/(∂2u2)2¯](1/2)$ is the Taylor microscale, and $u′rms=(u′iu′i¯/3)1/2$ is the root-mean square value of the velocity fluctuations. The dynamic viscosity μ obeys the power law $μ=μref(T/Tref)3/4$⁠, where $μref$ and $Tref$ are the reference viscosity and temperature, respectively. For the STI, the mean flow Mach number is $M=ũ1,u/c̃$⁠. In these definitions, Reynolds and Favre averages of a flow quantity, f, are indicated as $f¯$ and $f̃$⁠, respectively. The corresponding fluctuations are $f′=f−f¯$ and $f″=f−f̃$⁠.

The STI simulations are conducted with $Reλ=40, Mt=0.2$ and $M={1.5,3.0}$⁠, where $Reλ$ and M_t values are taken just upstream of the shock location. With this choice of flow parameters, the STI is in the wrinkled-shock regime.⁶⁴ The initial values of $Reλ$ and M_t in the corresponding DHIT configuration are chosen as $(Reλ,Mt)=(46,0.23)$ for the STI with M = 1.5 (denoted STI-M1.5) and $(Reλ,Mt)=(42,0.21)$ for the STI with M = 3.0 (denoted STI-M3) to account for the decay of $Reλ$ and M_t from the inlet to the shock location in the STI. The Schmidt number, defined as the ratio of the kinematic viscosity $ν=μ/ρ$ of the fluid and the diffusivity D_m of the mth passive scalar, is chosen to be Sc = 1 in all presented DHIT and STI simulations, characteristic of supersonic mixing of gases. For the STI flow configuration, a resolution of $1280×384×384$ for the domain of size $4π×2π×2π$ was found to be grid converged for the range of flow parameters considered here (see Ref. 15). The DHIT has a cubic domain of edge length $2π$⁠, and the grid resolution in each coordinate direction of the DHIT configuration matches the transverse resolution of the STI simulation.

Precursor DHIT simulations are initialized following the work of Ristorcelli and Blaisdell,⁶⁶ with the initial spectra for the velocity field set to $E(k)∝k4 exp (−2k2/k02)$⁠, where $k0=4$ is the wavenumber at which the energy peaks. The DHIT simulations are performed up to a time instant in which the turbulence is fully developed, and the Taylor microscale Reynolds number and the turbulence Mach number match their target values of $Reλ=46$ (42) and $Mt=0.23(0.21)$⁠, respectively. No shocklets are observed at the turbulence Mach numbers considered in the DHIT.

At that time instant, and for each simulation, three passive scalar fields, 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. To avoid discontinuities, the scalar fields are initialized using a smooth radial profile that transitions between zero $(ϕmin)$ outside and unity $(ϕmax)$ inside the targeted spherical shapes as follows:

where the subscript $α={1,2,4}$ indicates the scalar field with $Nα$ initial spheres of radius $αλ$ uniformly distributed (a distance $3αλ$ apart) in each coordinate direction. For each sphere, $Rα$ is the inner radius where the passive scalar is unity, and $ΔRα$ is the radial distance over which the scalar transitions to zero away from the inner radius. $rαk=x′ix′i$ is the radial coordinate relative to the center $x0i(k)$ of kth sphere, and $x′i=xi−x0i(k)$⁠.

The three initial configurations of passive scalar structures studied in this work are summarized in Table I. The same sampling interval length of $Δt=τ0/12.7$ is used for the tracking for all structures in the STI with M = 1.5 and corresponding DHIT, where $τ0=λ/urms$ is the initial large-eddy turnover time of the DHIT configuration. To account for the higher upstream velocity, the sampling frequency is doubled in the STI with M = 3.0 and in the corresponding DHIT.

To generate the STI inflow database, three individual DHIT precursor simulations without passive scalars are performed until reaching the target values $Reλ=46$ (42) and $Mt=0.23(0.21)$⁠. From each of those three time instants, two additional DHIT realizations are obtained by applying rotations with respect to the y and z coordinate axes. The resulting nine DHIT realizations are blended,⁶⁷ obtaining a $18π×2π×2π$ cuboid grid in which the scalar structures are initialized as previously described for the DHIT (Fig. 1). The same initial geometry is used for the DHIT and STI study to compare the passive scalar mixing in both configurations (Table I). The outlet (back) pressure⁶⁴ is adjusted to fix the averaged shock location at a distance $Δx≈2$ downstream of the inlet. This averaged shock location will be taken as the origin of the x axis. A sponge layer⁶⁴ of length π is introduced in front the outlet to avoid acoustic reflections from the outflow boundary into the subsonic region of the domain. Therefore, the effective streamwise length of the STI configuration relevant to the tracking is reduced to the region starting at the position downstream of the inlet where the structures first become closed surfaces, which depends on the initial size of the structures, to the starting position of the sponge layer.

The shapes and the number of individual scalar structures extracted from the underlying scalar field depend on the contour value $ϕiso$ used for iso-surfacing. In the STI configuration, increasing $ϕiso$ from $ϕmin$ breaks larger connected structures into disconnected smaller ones and the number of extracted structures is increased. For all three scalar types increasing $ϕiso$ beyond 0.3, the gain of structures due to this breakup process is counterbalanced by structures going below the threshold used for iso-surfacing, resulting in an overall decay of the number of extracted structures. Due to the structure deformation induced by the shock, the immediate post-shock region is prone to the formation of widely connected isosurfaces and therefore prescribes the percolation limit⁶⁸ of the system. The formation of largely connected structures jeopardizes the goal of investigating individual structure evolution. For the given initial spacing, an isovalue of $ϕiso=0.6$ is chosen to define the passive scalar structure in all presented configurations, which decreases the initial size of the structures insignificantly compared to the nominal radius of the spheres, but prevents the formation of large, interconnected isosurfaces in the post-shock region. An isovalue greater than 0.6 would simplify the structure tracking by reducing the number of educed structures and their interactions, but causes an undesirable loss of information. 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 tracking to ensure a reliable calculation of the differential-geometry properties from the surface triangulation (see Sec. III).

To analyze the dynamics and evolving characteristics of flow features (e.g., passive scalar isosurfaces) in physical space, we propose a local, hybrid attribute- and region-based, explicit tracking methodology (Fig. 2). Local methods assume a limited motion of features between consecutive frames in contrast to global methods.^69–71 Explicit tracking requires the independent extraction of flow features from a series of time instants (frames), followed by the matching of extracted features for each pair of consecutive frames (solving the correspondence problem). Region-based algorithms rely on the spatial overlap of features in consecutive frames, requiring a fairly high sampling frequency for datasets that contain small, fast-moving features.^72–75 Instead, attribute-based algorithms solve the correspondence problem using a set of abstract attributes to describe the features.^76–79 Our proposed hybrid method does not impose any requirement of spatial overlap of corresponding structures between frames and is therefore less restrictive in terms of the sampling frequency by combining the spatial location and non-local geometric attributes of flow features and introducing a set of physically realizable constraints.

In the matching process, flow structures in consecutive time frames t_n and $tn+1$ will be referred to as “source” and “target” structures, respectively. Interactions of source and target structures within a sampling interval are described by the following types of events.⁷⁶ In a “creation” event, a new structure appears in the target frame. In a “continuation” event, a structure continues from one frame to the next without interacting with another structure. In a “split” event, a source structure breaks into two or more target structures. In a “merge” event, two or more source structures coalesce into one target structure. In a “disappearance” event, a source structure vanishes. The complexity of events can depend on the sampling interval. For example, as the sample interval increases, a sequence of two first-order splits would eventually be detected as a single second-order split. Large sampling intervals create a need to additionally consider “compound” events, in which multiple source structures merge and split into multiple target structures.⁸⁰ 

The attributes used to find correspondences include the dimensionless center of the oriented bounding box (OBB) ${xc,yc,zc}$ of a structure and three dimensionless geometric attributes: the global compactness parameter $λ̂$⁠, the feature shape index $Ŝ$⁠, and curvedness $Ĉ$⁠. The compactness is defined in terms of the structure volume V and area A by

where the normalization factor is chosen such that $λ̂=1$ for spheres. The non-local quantities $Ŝ$ and $Ĉ$ are obtained from first- and second-order moments of the area-based joint probability density function of two pointwise differential geometry properties:⁶³ the absolute value of the shape index, S, and the dimensionless curvedness, C, defined in terms of the principal curvatures ${κ1,κ2}$ as⁸¹ 

For brevity, the dimensionless curvedness and the absolute shape index are referred to as curvedness and shape index. Larger curvedness implies a shorter radius of curvature and, hence, a more curved surface. The shape index classifies surface points as hyperbolic S < 0.5, parabolic S = 0.5, or (concave or convex) elliptical $0.5<S<1$⁠. S = 0 and S = 1 correspond to symmetrical saddle and umbilical points, respectively. $S∈[0,1]$ is undefined for a locally planar surface, for which C = 0 (note $C∈R+$⁠). Regions of the space of geometric attributes ${Ŝ,Ĉ,λ̂}$ are associated with particular surface shapes: blob-like structures concentrate near the $(Ŝ,Ĉ,λ̂)=(1,1,1)$ point, tube-like structures are near the axis $(1/2,1,λ̂)$⁠, and sheet-like structures arise for lower values of $Ĉ$ and $λ̂$⁠.

Hence, each structure will be represented as a point in a six-dimensional attribute space ${xc,yc,zc;Ŝ,Ĉ,λ̂}$ that combines spatial and geometric shape attributes. The spatial attributes are nondimensionalized with the domain size in each direction. All structures in a frame form a cloud of points in this attribute space. The proposed tracking method is applicable to any other choice of feature attributes.

A search for correspondences between each pair of source and target point clouds is conducted utilizing first the full set of attributes and, second, only its spatial subset, $x={xc,yc,zc}$⁠. The search for correspondences is performed in multiple stages in which both the source and the target clouds participate as query and searched clouds. An estimate $x′$ of the spatial attributes in the searched cloud is obtained by a linear projection based on the surface-averaged velocity, $vavg$⁠, and the tracking step $Δt$⁠,

The projection is performed forward (+) or backward (–) in time from the source and target frames, respectively.

Structures involved in continuation events do not interact with other structures and, thus, are generally characterized by small changes in geometric shape compared with splits and merges. Therefore, for a continuing source structure, a corresponding target structure is expected in its neighborhood in the full attribute space, and candidate targets can be detected by a nearest-neighbor (NN) search^82,83 in this space. The NN search is conducted bidirectionally: the nearest neighbor in the target cloud is found for each source point and vice versa. The bidirectional, unrestricted NN search guarantees that all structures have at least one corresponding structure, which can lead to falsely detected correspondences for disappearing source structures or newly created target structures. These false correspondences will be rejected by the constraining process, as explained in Sec. III B.

The NN search for each source (target) point predominantly detects targets involved in potential continuation or merge (split) events.

Since topology-changing (merge, split, and compound) events affect the geometric attributes more abruptly than the spatial attributes, an additional correspondence search in the spatial subspace is conducted excluding the geometrical attributes.

This search is conducted as a bidirectional range, rather than a NN search, and thus requires the definition of a search radius. An individual search radius is assigned to each query point based on the half-diagonal of its OBB, with a provided tolerance $εR$ (Fig. 3). Since the search radius is associated with the query structure size, splits (mergers) are more likely detected when a source (target) structure is used as a query. Hence, prior to the radius search using the source (target) points as query, the target (source) backward (forward) projection of spatial attributes is conducted.

A confidence value $ci∈[0,1]$ is defined for each found correspondence and search stage. For the bidirectional NN search, the confidence value is computed as

where $ωN N$ is a specified weight, $cN N<1$ is a minimum NN confidence value, and $dN N,i$ is the Euclidean distance between query and match points in the attribute space for the ith found correspondence, whereas $d¯N N$ and $σN N$ are the mean and standard deviation of the distances of all found correspondences in this search stage. For the radius search, the confidence is defined as $cR,i=ωR(1−dR,i/ri)$ where ω_R is a specified weight, $dR,i$ is the center distance between the matched structures and r_i is the search radius of the query. Correspondences between the same pair of source and target structures found more than once in this multi-stage search are merged into a single correspondence with a cumulative confidence value.

Correspondence constraints reject single correspondence candidates deemed not physically realizable, avoiding the formation of unphysical, complex events.

A sequence of region-based constraints of increasing accuracy and computational cost is applied to each candidate correspondence and a confidence value is computed in each constraint.

In a first constraint, the minimal spatial distance between the OBBs of source and target structures involved in a correspondence serves as a measure for the realizability (Fig. 4). The distance between the source and target OBBs (which are convex sets) is calculated using Gilbert–Johnson–Keerthi's (GJK) algorithm^84,85 based on Minkowski's separating axis theorem (SAT).⁸⁶ The OBB distance is compared to a measure defined as the maximum of $μ̂=3V/A$ of the source and target structures. If the minimal distance of the OBBs is larger than this measure (up to a specified tolerance), the correspondence is rejected. A confidence value is assigned to passing correspondences as

where $dmin,i$ is the minimum distance between the OBBs, and $ωO B B$ and $εO B B$ are specified confidence weight and tolerance values, respectively. The measure of the source and target structures is indicated as $μ̂s,i$ and $μ̂t,i$⁠, respectively.

Since OBBs are inaccurate representations of highly convoluted structures (Fig. 5), a second constraint is applied to the remaining correspondences. This more accurate and computationally expensive constraint solves distance queries between the actual isosurfaces, rather than between their OBBs, by calculating the minimum of point-to-point distances⁸⁷ between the surfaces. Again, the distance between structures is compared to the maximum of the $μ̂$ parameter of source and target structures to determine the confidence value of passing correspondences. The confidence value follows Eq. (7), replacing $cO B B,i$⁠, $ωO B B$⁠, and $εO B B$ with the confidence, weight, and tolerance of the proximity constraint, respectively.

An event groups correspondences with common source or target structures. By inspection of the set of correspondences, a list of candidate events describing the structure interactions between source and target frames is formed.

A sequence of constraints is applied to test the physical realizability of the events and to examine if an event can be simplified into subevents, defined as subsets of the correspondences of the whole event (i.e., each accepted subevent satisfies the event constraint).

To avoid false merge events resulting from small disappearing structures proximal to other structures, a first constraint is applied that removes merging source structures with a number of points of the surface triangulation below a minimum threshold. Second, a conservation constraint is applied based on an attribute assumed to be conserved in the (sub)event within a specified tolerance. The conserved attribute must differ from those utilized in the correspondence search. For the tracking of passive scalar structures, their mass (estimated as the product of the surface-averaged density and the volume enclosed by the structure) is chosen as the conserved attribute used in the event constraining process. The realizability of the (ith sub)event based on a conserved attribute α is assessed by comparing the sum of the conserved attribute between the n source and m target structures of the (sub)event, for which a score $si,c$ is defined as

where $αi,js$ and $αi,jt$ are the values of the conserved attribute of the event's sources and targets, respectively. A (sub)event is considered realizable if its score is above a specified threshold. 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 subevents to be tested is given by $c=∑r=1nn!/(n−r)!r!$⁠. Similarly, detected compound events are tested for decomposition into simpler events based on the conserved quantity.

We define structures as closed surfaces (note that the geometric attributes $Ĉ$ and $λ̂$ include the volume in their definition). Surfaces intersecting periodic boundaries (thus open) can be closed by reconnection,⁶³ but they must also be duplicated on all crossing boundaries to ensure a correct search for correspondences [Fig. 6(a)]. To detect potentially boundary-crossing structures in the source (target) frame prior to the correspondence search, their location in the target (source) frame is estimated by the forward (backward) linear projection. The subsequent correspondence search does not distinguish between duplicated and original structures. Therefore, multiple correspondences can result that connect the same source and target structures via their periodic duplicates, but represent the same fundamental correspondence [Fig. 6(b)]. Duplicated correspondences are thus removed after the search [Fig. 6(c)].

Graph data structures represent pairwise relationships between objects and are therefore ideal to organize the temporal evolution of flow features by mapping accepted events between pairs of frames. Each vertex of the graph represents a flow structure at a certain instant in time, while each edge indicates a correspondence. Since vertices are ordered in temporal sequence, the graph is directed and acyclic. An implicit vertex clustering results in a two-layer hierarchy of subgraphs, defined as tree and branch (Fig. 7). Each branch represents the lifetime of one individual structure. Depending on the structure interaction, some trees in the graph may be disjoint from all other trees. The union of all trees describes the evolution of all structures over all frames of the dataset.

The source (target) vertex of a split (merge) has an out (in) degree greater than one, i.e., more than one correspondence with target (source) structures. Multiple vertices of degree higher than one exist in a compound event. The vertex clustering approach aims to continue branches in the building of these event types. Therefore, the dominant target (source) of each split (merger) must be selected based on a set of criteria. For a compound event, the process is applied to the sub-split and sub-merge events to select dominant sources and targets.

Allowing n criteria, the total score of the ith candidate correspondence is calculated as a weighted sum of the scores of each criterion,

where ω_j is the specified weight for the jth criterion and s_ij is the score of the ith candidate in the jth criterion.

The first criterion to decide on the dominant edge(s) of an event is based on the comparison between a specified conserved property C of the candidate and the reference, such as the volume or (estimated) mass of the structure. For a (sub-)split or (sub-)merger, the reference is the splitting source or the merged target, respectively. The score of the ith candidate in this conservation-based criterion is determined as

The second criterion considers the confidence value of each correspondence as the score to prioritize high confidence in the selection of the dominant edge.

For (sub-)merge events, a third criterion based on the lifetime of the source structure before the event is used. The lifetime-based score of the ith source candidate is defined as

where $ti,0$ is the time of the frame in which the ith candidate is first present as an individual structure. Finally, the correspondence candidate that maximizes the total score, $si,total$⁠, is selected as the dominant source or target of the (sub-)merge or (sub-)split event.

Branches are classified based on the connectivity of their start and end vertices.⁵⁹ The start vertex of a primary branch results from a creation event (with no predecessor), whereas its end vertex is a disappearing structure (with no successor). Branches that do not satisfy both criteria are considered secondary (Fig. 7).

We apply the tracking algorithm introduced in Sec. III to analyze the temporal evolution of passive scalar structures extracted from the DHIT and STI numerical datasets discussed in Sec. II. Figure 8 shows the initialized structures and their temporal evolution in the DHIT at three different time instants. The tracking parameters regarding the correspondence search, event formation, constraints, and graph building stages, listed in Tables II–IV, are chosen by inspection of the detected events and applied to all simulation cases. For small structures, defined by having a number of surface triangulation points smaller than 1000, the correspondence constraining, Sec. III B, applies a relaxed fraction $εrel=2ε$ to provide a higher tolerance accounting for rapid dissipation on small scales. Similarly, the conservation-based constraint skips splitting and continuing source structures with less than 1000 surface triangulation points.

We first compare the temporal (DHIT) and spatial (STI) distributions of different event types in Sec. IV A, and the departures from volume preservation for each event type formed by non-material scalar structures in Sec. IV B. We then focus on geometric changes of the structures undergoing split and merge events in Sec. IV C, the correlation of the structure geometry with physical quantities relevant to mixing in Sec. IV D, and the temporal trajectories in the space of geometric parameters in Sec. IV E.

To compare the structure evolution in DHIT and STI, the time coordinate of the DHIT simulations is translated into an equivalent streamwise coordinate of the STI as

where x_s is the shock location, and u_u and u_d are the mean streamwise flow velocities upstream and downstream of the shock.

The equivalent-streamwise distribution of different event types in STI and DHIT configurations is presented in Fig. 9, where the streamwise coordinate is normalized by the k₀ wavenumber and the initial Taylor micro-scale, λ, of the DHIT configurations. Expected creations correspond to structures entering the STI domain and to the initial structures in the DHIT. Scalar dissipation eventually reduces the scalar value below the isosurfacing threshold, resulting in the disappearance of the structures. The peak location of creations is shifted downstream as the structures increase in initial size from $ϕλ$ to $ϕ4λ$ and become closed surfaces (triggering the creation) at a farther downstream location (Fig. 9). In the post-shock region, a small number of structures are falsely detected as newly created. Those structures are small in volume and correspond to targets of continuations or (sub-)splits. Due to passive scalar dissipation, once a structure vanishes below the iso-surfacing threshold, it will not reemerge (triggering new creation events) since the scalar value will remain below the threshold.

The distribution of disappearances is linked to that of splits as most (early) splits have one dominant structure similar in size to the parent structure (see Sec. IV E) and a number of smaller targets that disappear in close proximity to the split event. The normalized number of splits and disappearances is larger in the STI than in the DHIT configuration, indicative of shock-induced mixing enhancement, and further increases with the shock Mach number and the initial size of the structures. For all initial sizes, a higher shock Mach number induces more splits in closer proximity to the shock. While disappearances and splits peak approximately at the same streamwise location for $ϕλ$ and $ϕ2λ$ in STI and DHIT, they peak closer to the shock for the $ϕ4λ$ structures, which indicates increased benefit of mixing of those larger structures by the shock interaction. Disappearance events extend farther downstream than split events, corresponding to scalar diffusion leading to continuation and disappearance without further splitting.

The number of compound events is larger in STI than in DHIT, increasing with Mach number, and peaks closely after the shock. The sub-splits hidden in compound events contribute to the number of disappearing structures in Fig. 9 closely behind the shock, especially for the $ϕ4λ$ structures, as the number of compounds per incoming structure is larger for increasing initial size of the structure (Fig. 9). The normalized number of mergers is small compared to the number of splits. Similar to sub-splits, there are a considerable number of sub-merge events hidden in compounds. While the number of pure splits is clearly larger than the number of compounds, the latter is comparable or even larger than the number of pure mergers, especially in the STI configuration, where more compounds, and hence more sub-mergers, occur. While there is typically one sub-merge per compound event, multiple sub-splits of potentially higher order occur per compound (see Fig. 10 for representative examples). Compounds are of importance to capture all cases of structure interactions in the scalar mixing process. Therefore, it is essential for the tracking to account for compound events rather than neglecting them as in most previous methods.

Volume distributions of disappearing $ϕ4λ$ structures in STI-M3 and DHIT are shown in Fig. 11(a). 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 induced by the shock. While the largest disappearing structures in DHIT and STI are similar in size (right tails), the smallest disappearing structure is significantly larger in the DHIT than in the STI configuration (left tails). This is mainly attributed to an increased streamwise grid resolution of the STI immediately downstream of the shock while using the same minimum number of points for surface triangulations to be considered in the tracking.

A similar grid-resolution effect is observed in the joint source and target volume histograms of continuing $ϕλ$ structures in STI-M3 and DHIT [Fig. 11(b)]. A deviation from volume-preserving continuation events is observed on small scales, due to scalar dissipation, and for larger shock-crossing structures due to the compression of the structure and the scalar dissipation amplification by the shock. Consistent with the distribution of splits (Fig. 9), the number of splits per initial/incoming structure is larger for the initially bigger spheres in $ϕ4λ$ (Fig 12). The splitting source volume, V_s, can be larger than the initial structure volume due to prior merge events. While the largest splitting volumes are naturally found for the $ϕ4λ$ structures, the smallest splitting structures have similar volumes for both scalar fields as the dissipation equally affects both fields on the same scales, owing to the same Schmidt number.

For both DHIT and STI, the smallest splitting volumes are significantly larger than the smallest disappearing volumes. Hence, 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 dissipation. Structures of small volumes, near the end of their lifetime, exhibit a relatively high compactness, $λ̂$⁠, which decreases the probability of a split (see Sec. IV E). Figure 12 shows that large splitting structures tend to conserve the volume in split event. In contrast, for small splitting structures, with a volume below 10% of that of a sphere of radius equal to the Taylor microscale, the sum of the target volumes is less than the source volume due to increased scalar dissipation at smaller scales.

The geometric signature ${Ŝ,Ĉ,λ̂}$ of source and target structures of split and merge events is investigated, focusing on the effects of the initial structure size and the shock with Mach number M = 1.5 (similar conclusions were found for the higher Mach number M = 3.0). The joint histograms $(λ̂,Ĉ)$ and $(Ŝ,Ĉ)$ of the geometrical attributes of source and target structures of split events in DHIT and STI, shown in Fig. 13, indicate that splitting source structures are characterized by lower values of the compactness, $λ̂$⁠, and the shape index, $Ŝ$⁠, than the resulting target structures, while the curvedness, $Ĉ$⁠, is comparable between sources and targets. Structures characterized by low compactness, $λ̂$⁠, are more stretched and likely to contain a ligament-like portion (e.g., see Fig. 10), which increases the probability of these structures becoming the source of a split.

The compactness $λ̂$ of splitting sources in DHIT (Fig. 13) 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 translates into a higher curvedness, $Ĉ$⁠, and an increased amount of hyperbolic surface points, which lowers the shape index (⁠$Ŝ<0.5$⁠). In contrast, elliptical points (⁠$Ŝ>0.5$⁠) are dominant in splitting $ϕλ$ structures.

For $ϕ2λ$ and $ϕ4λ$ structures, the (⁠$Ŝ, Ĉ$⁠) joint PDF of source structures is essentially enclosed by that of the target structures (Fig. 13). As the analysis of geometric trajectories presented in Sec. IV E will show, split events early in the structure lifetime typically contain a dominant target structure 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. Compared to the DHIT, the compression of the passive scalar structures by the shock in the STI decreases the compactness, $λ̂$⁠, and shape index, $Ŝ$⁠, of splitting source structures (Fig. 13).

The $(λ̂,Ĉ)$ and $(Ŝ,Ĉ)$ joint PDFs of the geometrical attributes of source and target structures of merge events in DHIT and STI (not shown) present opposite conclusions than for splits. Merged target structures are characterized by a lower compactness. The difference between merging sources and merged targets in terms of the shape index is less pronounced than for split events. The dependence of the geometric signature on the initial structure size observed for splitting source structures also applies to merged target structures.

The local geometry of a surface can be characterized by the pointwise values of curvedness C and the shape index S on the surface [Figs. 14(a) and 14(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 [Fig. 14(c)]. To elucidate correlations of local physical quantities on a passive scalar isosurface with its local geometry, the average alignment of $∇ϕ$ and $γ$ conditioned on S and C accounting for the area-based joint PDF is obtained [Fig. 14(d)]. For the sample structure shown in Fig. 14, a clear correlation between flat regions (low C) and high alignment of $∇ϕ$ and $γ$ is observed on the surface renderings [Figs. 14(a)–14(c)] and the conditioned joint PDF [Fig. 14(d)].

Figure 15 shows surface distributions, conditioned on S and C, of alignments between $∇ϕ$⁠, the strain-rate eigenvectors and vorticity, and the scalar gradient magnitude obtained from $ϕ4λ$ isosurfaces centered around four different streamwise intervals in the STI-M1.5 simulation. Positive correlations with C < 5.0 are observed for the alignment of $∇ϕ$ (which is normal to the isosurface at every point) with $α$ and $β$⁠, indicating that regions of higher curvature present better alignment between the scalar gradient and the most extensional and intermediate strain-rate eigenvectors. The opposite trend is observed for the alignment with the most compressive strain-rate eigenvector, $γ$⁠, indicating that flatter regions of the isosurfaces present larger alignment between $∇ϕ$ and $γ$⁠. For C < 5.0, these correlations are strengthened as the isosurfaces are transported downstream. For highly curved regions (C > 5.0) of structures in the post-shock domain, alignment with $α$ and $β$ slightly decreases, whereas alignment with $γ$ increases. Across the shock, the scalar gradient alignment with $α$ (⁠$γ$⁠) increases (decreases) for the entire range of the shape index, S. A decreased alignment with $β$ downstream of the shock is more pronounced at parabolic points (S = 0.5) of the isosurfaces.

For the alignment of scalar gradient with the most extensive strain-rate eigendirection, $α$⁠, the correlation with S is less pronounced than with C: at the farther 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 eigendirection, $β$⁠, farther downstream. Contrary to $α$ and $β$⁠, alignment of $∇ϕ$ with the most compressive strain-rate eigendirection, $γ$⁠, shows a negative correlation with the shape index for hyperbolic points. In the elliptical range of shape index (S > 0.5), however, the correlation is positive, but becomes less pronounced farther downstream. The alignment between scalar gradient and vorticity, $ω$⁠, shows a positive correlation when conditioned on C and with S in its hyperbolic range. Hence, the anti-correlation between $∇ϕ$ and $ω$ is more pronounced on flatter surface regions.

The correlation between the magnitude of $∇ϕ$ and the local geometry exhibits similar trends to those found for the alignment between $∇ϕ$ and $γ$⁠. This implies that most of the dissipation after the transient and its enhancement after crossing the shock occurs on flatter regions (lower C). The ordering of the ensemble means by streamwise location agrees with the streamwise evolution of the surface-averaged scalar gradient magnitude described in the  Appendix. Similar correlations are found in the DHIT and STI-M3 configurations and for all initial structure sizes (not shown).

The deformation of passive scalar structures is analyzed by means of their trajectories in the geometrical feature space ${Ŝ,Ĉ,λ̂}$⁠. To focus on the temporal evolution of primary structures, we consider only graph branches whose start vertex results from a creation event and end vertex is a disappearing structure.

Initially, all $ϕ4λ$ structures in the DHIT are spheres, and thus their trajectories have a common starting point in the feature space at $(Ŝ,Ĉ,λ̂)=(1,1,1)$⁠, as shown in Fig. 16. These $ϕ4λ$ isosurfaces become highly convoluted with time (see also Fig. 8), causing a rapid increase in the curvedness, $Ĉ$⁠, and a decrease in the shape index, $Ŝ$⁠, due to the occurrence of hyperbolic surface points (S < 0.5). The deviation from the spherical shape also lowers the compactness, $λ̂$⁠, of the structures as they are stretched. After reaching the maximum curvedness, all structures follow a nearly linear trajectory in the $Ĉ$-$λ̂$ plane toward the approximate $(0.7,0.25)$ location, lowering both their curvedness and compactness. The shape index, $Ŝ$⁠, follows a rapid initial transition from 1 (spherical shape) to 0.5, reached as the maximum $Ĉ$ is attained. The subsequent decrease in $Ĉ$ and $λ̂$ has a lesser effect on $Ŝ$⁠, indicating a balance between the surface area coverage of hyperbolic and elliptic local geometries on the convoluted surfaces.

Until the occurrence of the first non-continuation events, the change in shape from one frame to the next (given by the distance between consecutive trajectory points in Fig. 16) monotonically decreases due to the viscous decay of the underlying turbulent background field. As the stretching increases (for $λ̂<0.5$⁠), the first splits occur in the evolution of the structures. Highly stretched structures tend to split over a sequence of frames. However, despite this cascade of split events, the change in shape of the primary structure itself is comparable to that of a continuation event. Splits occurring in this early stage of the evolution tend to be dominated by one target structure (similar in shape to the splitting source structure) and other, much smaller targets. As the dominant target continues the branch of the primary structure, its geometrical signature is only moderately altered by the split. This geometric similarity between the source and dominant target structure is also observed in the joint PDFs of geometrical attributes (Fig. 13).

After these split cascades, the structures have reached their minimum compactness, and subsequent split events tend to increase the compactness and the shape index of the primary structure while preserving its curvedness, consistent with the joint PDF of splits in the geometrical feature space discussed in Sec. IV C. These later splits do not form a cascade even though the branch section between splits characterized by a decreasing curvedness and compactness remains short in some cases. The final stage of primary structure lifetime is characterized by continuation events with increasing amount of elliptical surface points and increasing compactness, making further splitting of the structure less likely. This trend is in accordance with the streamwise distribution of splits compared to disappearance events, discussed in Sec. IV A.

From the collection of individual geometrical trajectories, ensemble statistics are calculated to obtain the mean trajectory and its standard deviation for primary structures of each scalar field (⁠$ϕλ, ϕ2λ$⁠, and $ϕ4λ$⁠) in DHIT and STI (M = 1.5, and 3.0). Ensemble trajectories, shown in Fig. 17, enable a characterization of the effects of the initial size and the shock compression on the geometric evolution of the structures.

While all three mean trajectories have the spherical starting point $(λ̂,Ŝ,Ĉ)=(1,1,1)$⁠, the initial evolution of the curvedness varies with initial size. Whereas the ensemble curvedness $Ĉ$ initially increases for medium (⁠$ϕ2λ$⁠) and, especially, large (⁠$ϕ4λ$⁠) structures, it decreases for smaller (⁠$ϕλ$⁠) structures as they are less convoluted by the background turbulence eddies. With increasing initial structure size, the ensemble mean trajectory reaches lower values of the shape index, $Ŝ$⁠, favoring hyperbolic points, and also lower compactness, $λ̂$⁠, indicative of an increased convolution and stretching. At the end of their lifetime, structures of all the scalar fields considered tend toward the same geometrical attributes: moderate compactness and curvedness with predominantly elliptical points.

The dependence of surface convolution and stretching on initial structure size is reflected in the mean area evolution of primary structures in DHIT (Fig. 18), suggesting a degree of self-similarity. Exponential area growth has been reported for the evolution of Lagrangian structures in statistically stationary turbulence.^88–92 In the present study, despite the temporally decaying turbulence and diffusion acting on the (non-Lagrangian) passive scalar structures, an exponential time dependence of the form $A(t/τ0)/A0=a exp (b t/τ0)$ accurately models the early, stretching-dominated surface area evolution, with a and b constant. A larger initial area growth is observed for initially larger structures and Taylor-Reynolds number (⁠$Reλ=46$ compared with 42, corresponding to the inflow turbulence used for STI cases with M = 1.5 and 3.0, respectively). After the 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 diffusion-dominated regime, the area decrease also follows an exponential decay. 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 significantly, occur in this period of exponential surface area decay.

In the STI configurations, the structures are already deformed when they become closed surfaces upon entering the computational domain. Therefore, their mean trajectories in the geometrical feature space do not have a common spherical starting point and differ based on initial structure size (see center and bottom rows of Fig. 17). In the M = 3.0 case, consecutive $ϕ4λ$ structures in the streamwise direction, decelerated by the shock, merge with each other and pile up within a distance from the shock below the initial structure radius, generating structures of large streamwise extent downstream of the shock. The number of primary $ϕ4λ$ structures with very long lifetimes in the STI configurations (M = 1.5 and 3.0) is too small to be statistically significant, so the final steps of their lifetime are excluded from the mean ensemble trajectory. Therefore, the final values of the geometrical properties attained at the end of the mean trajectories of $ϕ4λ$ structures differ from those of the smaller $ϕλ$ and $ϕ2λ$ structures.

Shock-enhanced scalar dissipation lowers the lifetime of the primary structures in comparison to the corresponding DHIT as the structures are diffused faster below the threshold used for iso-surfacing (Table V). Compared to the DHIT, structure compression induced by the shock in the STI cases lowers the curvedness $Ĉ$ (flattening the surface in directions normal to the compression) and the compactness $λ̂$ (resulting in higher stretching, thus favoring splits). The deformation induced across the shock, particularly for the higher 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-M3 case, the curvedness of the $ϕλ$ and $ϕ2λ$ structures reaches a minimum after the structures are completely processed by the shock and increases afterward due to shock-induced turbulence amplification.

The lifetime reduction indicates that the mixing of larger structures benefits more from the shock-interaction: for the STI-M1.5 case, the mean disappearance location of $ϕλ$ structures is reduced by 24%, while for $ϕ2λ$ structures, the reduction is 47%. For the higher Mach number case, STI-M3, the lifetime 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 presented in the  Appendix (Fig. 22).

Turbulence amplification across the shock causes a surface area re-growth of $ϕ2λ$ structures in the post-shock region that is stronger for higher Mach number (Fig. 18). This period of area re-growth correlates with the decreasing compactness parameter, $λ̂$⁠, in the feature-space trajectories (Fig. 17) after the structures have been processed by the shock. While the amplification of the surface-averaged vorticity magnitude is similar for all initial structure sizes (Fig. 22), the area re-growth is larger for intermediate (⁠$ϕ2λ$⁠) than for small (⁠$ϕλ$⁠) 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 (shock-amplified) 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 significantly.

Ensemble trajectory statistics are extended to include the evolution of two physical quantities averaged on each surface: the magnitude of the scalar gradient (indicative of the scalar dissipation on the surface) and the vorticity magnitude. These ensemble trajectories in the parameter space that combines geometrical features and averaged physical quantities are shown in Fig. 19 for the DHIT case. For $ϕ4λ$ structures, the scalar gradient magnitude exhibits a rapid initial increase due to stretching of the initial scalar transitional region. The peak in curvedness is reached significantly earlier than the peak in the scalar gradient magnitude. The geometry of $ϕ4λ$ structures, especially their curvedness and shape index, develops 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 splits, not only is the geometrical signature of the primary structure marginally affected, 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, reaching lower values on the larger $ϕ4λ$ structures due to their longer lifetime.

The individual trajectory closest to the ensemble mean in terms of geometry and scalar gradient magnitude (and thus, scalar dissipation) is also shown in Fig. 19 for $ϕ4λ$ and $ϕλ$ structures. Figure 21 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 onto the surface. The area-based joint and marginal probability density functions of S and C for the structures at the last highlighted step are also shown, indicating that for both structure sizes the largest curvedness on the surfaces is found at parabolic points (S = 0.5) that occur on the ridges of the convoluted surfaces. Previous studies have reported highly curved surface elements of propagating surfaces in isotropic turbulence to be of cylindrical shape,⁹⁰ consistent with such parabolic points found presently. The correlation, found in Sec. IV D, between flatter 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. This correlation can be observed in the joint and marginal distribution of averaged alignment conditioned on C shown in Fig. 21.

Mean trajectories of primary structures in STI with M = 1.5 differ 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 [Fig. 20(a)]. The exclusion of the final steps of the lifetime for lack of statistical significance prevents the mean $ϕ4λ$ trajectory from reaching values of the scalar gradient and vorticity magnitudes as low as in the trajectories of $ϕλ$ and $ϕ2λ$ structures. Whereas the increase in scalar gradient magnitude depends on the size of the primary structures, the increase in the vorticity magnitude is the same for all structure sizes.

The stronger compression and quicker closure of incoming structures in the STI-M3 configuration change the ensemble mean trajectories in terms of geometrical and physical structure properties compared to the lower Mach number case [Fig. 20(b)]. The smaller the initial size of the structures, the lower their curvedness is at the point of highest scalar dissipation. Compared to the DHIT (Fig. 19), $ϕ4λ$ structures reach significantly higher average scalar dissipation, and at lower curvedness resulting from the flattening of the surfaces 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 pileup 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 final stage of their ensemble mean trajectories.

The average alignment between scalar gradient and the most compressive strain-rate eigenvector conditioned on S and C at highlighted steps of the closest-to-mean trajectories (Fig. 21) emphasizes again the found correlation between flatter surface regions and high $∇ϕ$-$γ$ alignment for both initial structure sizes.

A hybrid approach of attribute- and region-based tracking has been introduced to study the temporal evolution of flow features, defined as surfaces of any kind in a turbulent flow field. Spatial and geometrical feature attributes are utilized to generate candidates of temporal feature relation. A sequence of region-based correspondence and attribute-based event realizability constraints has been proposed to test the physical plausibility of detected correspondence candidates. A graph data structure is built to organize the evolution and record the interactions of individually tracked flow features over time. The graph data structure is analyzed and queried to gain insight into the temporal flow feature evolution.

The tracking methodology was applied to investigate the dynamics of passive scalar structures in temporally decaying HIT and statistically stationary STI at Mach numbers 1.5 and 3.0. The flow features of interest were defined as isosurfaces of the underlying passive scalar field. The comparison of both flow configurations allowed us to elucidate the effect of shock waves on the scalar mixing of finite-sized structure through their interactions and evolution of surface-mapped physical quantities relevant to mixing.

Shock-induced deformation facilitates the structural breakup process increasing the occurrence of splits and disappearances, shifting those events closer to the shock for the increasing Mach number. As expected, initially larger structures take longer to diffuse, but their mixing benefits 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 significant than for smaller structures.

Small-structure interactions deviate from volume-preservation as a consequence of increased scalar dissipation at fine scales, whereas events involving larger structures show a high degree of volume preservation. Structures that have shrunk to a certain volume (⁠$10−2$ and $10−3$ times the volume of the Taylor-scale sphere for DHIT and STI, respectively), 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 diffusion. The geometric signature of structures involved in splits and mergers shows 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 configuration.

Larger scalar gradient (hence, dissipation), its shock-induced amplification, and alignment with the most compressive strain-rate eigendirection were found to correlate with flatter surface regions. These correlations emerged shortly after the early 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 flow). 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 configurations, lowering the compactness and curvedness of the structures. Less compact (i.e., more stretched) structures are more likely to split, while a lower curvedness increases the proportion of flatter surface regions where scalar dissipation dominates. The sudden shock interaction flattens the structures and enhances mixing more for larger structures. The post-shock amplification of turbulence increases the curvedness again, inducing a period of re-growth of the surface area and stretching of the structures that increases the probability of splits and enhances scalar mixing.

Financial support was provided by the Army Research Office (Grant No. W911NF2010096, Program Manager M. Munson) and the 2018 Summer Program at the Center for Turbulence Research (Stanford University). Computation for the work described in this paper was supported by the University of Southern California's Center for High-Performance Computing and an award of computer time provided by the Innovative and Novel Computational Impact on Theory and Experiment (INCITE) program. This research used resources of the Argonne Leadership Computing Facility, which is a DOE Office of Science User Facility supported under Contract No. DE-AC02-06CH11357. The authors are grateful to Professor Johan Larsson for providing the Hybrid numerical flow solver upon which extensions for passive scalar transport were implemented. Fruitful discussions with Profs. Sanjiva Lele and Lin Fu during the 2018 Summer Program at the Center for Turbulence Research, Stanford University are gratefully acknowledged.

Ensemble statistics of the streamwise evolution of surface-averaged physical quantities of the flow field computed on individually tracked isosurfaces of $ϕλ, ϕ2λ$⁠, and $ϕ4λ$ passive scalars in STI and DHIT are presented in Figs. 22 and 23. Since the effect of the shock on these quantities is more pronounced for a higher Mach number, only results of the STI configuration with M = 3.0 are included. The streamwise evolution of the ensemble statistics confirms findings from statistical analyses of the volumetric scalar field as summarized in Sec. I and is presented here for validation purposes.

Whereas the scalar gradient magnitude of $ϕλ$ structures decays immediately after the initialization, stretching of the (larger) non-zero scalar gradient region of the $ϕ4λ$ structures causes an initial increase in the scalar gradient magnitude. After the initial increase, the magnitude of scalar gradient monotonically decays in the DHIT configuration as the scalar mixing progresses.

In contrast, the surface-averaged scalar gradient magnitude is amplified as the structures pass through the shock wave in the STI configuration. This amplification, indicative of scalar dissipation enhancement, is stronger for a larger initial size of the structures. After the scalar gradient magnitude has reached its maximum, its decay is more pronounced than in the corresponding DHIT as the scalar dissipation rate is amplified by the shock, consistent with prior statistical analyses.^14,15

The authors have no conflicts to disclose.

The data that support the findings of this study are available from the corresponding author upon reasonable request.