Document [original]

AIP Advances 11, 085322 (2021); https://doi.org/10.1063/5.0062145 11, 085322

Spatial and temporal dynamics of a

supersonic mixing layer with a blunt base

Cite as: AIP Advances 11, 085322 (2021); https://doi.org/10.1063/5.0062145

Submitted: 02 July 2021 • Accepted: 29 July 2021 • Published Online: 17 August 2021

Lantian Li (李蓝天) and Hao Li (李浩)

ARTICLES YOU MAY BE INTERESTED IN

Monotonic variation in carbon-related defects with Fermi level in different conductive

types of GaN

AIP Advances 11, 085321 (2021); https://doi.org/10.1063/5.0054483

Simulation of convective combustion reactions in PBX based on DEM–CPM

AIP Advances 11, 085326 (2021); https://doi.org/10.1063/5.0062549

Abnormal threshold voltage shifts in p-channel low temperature polycrystalline silicon TFTs

under deep UV irradiation

AIP Advances 11, 085328 (2021); https://doi.org/10.1063/5.0060553

AIP Advances ARTICLE scitation.org/journal/adv

Spatial and temporal dynamics of a supersonic

mixing layer with a blunt base

Cite as: AIP Advances 11, 085322 (2021); doi: 10.1063/5.0062145

Submitted: 2 July 2021 •Accepted: 29 July 2021 •

Published Online: 17 August 2021

Lantian Li (李蓝天)1and Hao Li (李浩)2,a)

AFFILIATIONS

1Science and Technology on Scramjet Laboratory, College of Aerospace Science and Engineering, National University

of Defense Technology, Changsha, Hunan 410073, People’s Republic of China

2Institut für Strömungsmechanik und Technische Akustik (ISTA), Technische Universität Berlin,

Müller-Breslau-Straße 8, D-10623 Berlin, Germany

a)Author to whom correspondence should be addressed: [email protected]

ABSTRACT

A supersonic mixing layer with a blunt base is of practical significance to engineering. Two flow configurations with splitter thicknesses of

1 mm (TN) and 5 mm (TK) are simulated using large eddy simulation. The cluster-based network model (CNM) projects the supersonic

mixing layer into a ten-centroid based low-dimensional dynamical system. The CNM’s outputs of TN and TK cases are compared in order to

better understand the spatial and temporal physics. The given baseline case (TN) demonstrates a quasi-steady dynamics with a periodic visit

between ten centroids. Each cluster occupies a nearly uniform space region and is also populated with equal probability. The CNM identifies

ten centroids associated with these two flow regimes observed in the TK case: Kelvin–Helmholtz vortex and vortex pairing. According to the

resolved centroids, increasing the thickness of the splitter plate complicates the flow structures and expands the high-dimensional state space.

The CNM presents probable state transitions, revealing that the temporal dynamics in the whole field exhibits highly intermittent behaviors,

with large shape modifications but small fluctuations in turbulent kinetic energy. In the near-wake field, the reattachment point and shock

wave behave similarly that they move downstream and upstream alternatively. The blunt base supersonic mixing layer, in aggregate, increases

the turbulent kinetic energy by 20.5%.

(http://creativecommons.org/licenses/by/4.0/). https://doi.org/10.1063/5.0062145

I. INTRODUCTION

A mixing layer is a canonical flow observed in many prac-

tical engineering issues, such as noise reduction,1,2 drag reduc-

tion,3aero-optic effect,4and mixing enhancement.5–7 Over the past

few decades, numerous studies on mixing layers have been car-

ried out. Starting with the incompressible regime, the mixing layer

has been extensively investigated from the perspectives of veloc-

ity and pressure fields,9large-scale coherent structures,10,11 and

their merging/pairing in the entrainment process.12 Subsequently,

Bell and Mehta13 discovered that the incompressible mixing lay-

ers are able to maintain the self-similar state with the averaged

profiles of the first and second-order flow quantities. Following

the deep insights into the incompressible one, research interest

in compressible mixing layers was highly encouraged owing to

their wide engineering applications. Ortwerth and Shine11 exper-

imentally confirmed the presence of large-scale structures in the

supersonic mixing layer. However, due to the limitations of exper-

imental optical devices, fine structures cannot be captured. Goebel

et al.14,15 first discovered the self-similar state in the supersonic

mixing layer, which was later confirmed by the following simula-

tions16 and experiments.17,18 Bogdanoff19 first proposed the convec-

tive Mach number (Mc) to quantify the compressible effects, and

subsequent researchers have continued to recognize it. Apart from

the heat release,20,21 increasing Mc reduces the expansion of the

compressible mixing layer.15,22 These studies, like the ones men-

tioned above, are mostly concerned with flow visualization and

statistical analysis of flow quantities, with minimal emphasis on

dynamics.

The military value of hypersonic aircraft has recently attracted

much attention. The development of combined cycle engines is

critical as they are the most promising propulsion system for the

single-stage-to-orbit air-breathing launch vehicle and the reusable

recce/strike airplane platform.23 The engineering development of

AIP Advances 11, 085322 (2021); doi: 10.1063/5.0062145 11, 085322-1

AIP Advances ARTICLE scitation.org/journal/adv

combined cycle engines has reignited the research interest in the

fundamental studies24–27 and mixing enhancement techniques28,29

inthesupersonicmixinglayer.However,thesplitterplatewitharel-

ativelythintrailingedgeisgenerallyincludedintheexperiments30–32

for the study convenience, especially in the supersonic regimes since

the flow structures are pure at this moment. The velocity inlet in

simulations is assigned as a tanh-type profile without taking into

account a splitter plate.24–26 These methods are appropriate to con-

centrateonthecoreissuesofinterest,i.e.,asimplemixinglayer,nev-

ertheless,disjointfromtheengineeringpractice.Inacombinedcycle

engine, the fuel supply and protective cooling system are imperative

constituent parts, which incur the inevitable increase in the thick-

ness of the splitter plate. Thus, the supersonic mixing layer with a

blunt base deserves deep exploration.

Some pioneers have illustrated the characteristics of the super-

sonic mixing layer downstream a splitter plate via experimental

means33 and numerical methods.34,35 Flow structures and statis-

tic analysis are frequently entailed, but dynamical analysis is rarely

concerned. The nonlinear nature in the dynamical system and a

large amount of data from both experiments and numerical simu-

lations pose great challenges to understanding the dynamics. The

reduced-order model (ROM) is a key enabler for physical under-

standing. Laizet et al.35 simulated incompressible mixing layers hav-

ingvariousgeometriesofthicksplittersbymeansofdirectnumerical

simulation. The influence of the trailing edge shape is carefully con-

sidered on the spatial development of the mixing layer. In addi-

tion, the proper orthogonal decomposition (POD) is employed to

understand the transition from the wake to the mixing layer. Low-

dimensional models were also built using POD to understand the

temporal and spatial evolution of the mixing layer.36,37 However,

the POD is conducted by solving the eigenvalue problem on the

correlation matrix of any two fields. It tends to mix the tempo-

ral and spatial frequencies. The POD modes are closely associated

with the eigenvectors instead of the practical physics, which are

hard to be interpreted. Dynamic mode decomposition (DMD) is

an exciting tool to overcome these drawbacks. Song et al.38 con-

ducted DMD analysis to explore the sound generation mecha-

nisms in mixing layers. DMD successfully captures both the far-

field acoustics and the near-field dynamics with a few modes. Soni

and De39 executed DMD to segregate the coherent structures on

the basis of the pure frequency in the supersonic mixing layer.

The dominant modes capture the breakdown, pairing/merging, and

stretching of vortices to aid in understanding the contribution to

supersonic mixing. However, DMD modes cannot be ranked

according to their importance for flow physics. Due to the lin-

earization of the governing equations in the algorithms,40 the DMD

generally shows weak robustness to the nonlinear dynamics with

complicated dynamics, especially with highly intermittent dynam-

ics.41 Given the success and drawbacks of standard DMD, a few

variants of DMD have been proposed, such as recursive DMD,42

spatial–temporal DMD,43,44 and optimized DMD,45 just to name

a few.

Most drawbacks mentioned in the POD and DMD can

be avoided by means of the cluster-based reduced-order model

(CROM). Burkardt et al.46,47 first introduced an approach of cen-

troidal Voronoi tessellation (CVT) to produce a ROM of the

Navier–Stokes equation. Subsequently, Kaiser et al.48 proposed a

cluster-based Markov model (CMM) that integrates the cluster

analysis with a Markov model to identify the transitions between

different spatial modes. Then, the CMM is employed to analyze

the wake of a high-speed train,49 twisted cylinder flow,50 super-

sonic mixing layer,51 and turbine wake.52 Shahbazi and Esfaha-

nian53 built a ROM using cluster analysis and orthogonal clus-

ter analysis to analyze the lead-acid battery, which continuously

enrichedthe family memberofthe CROM. IntheCMM, the cluster-

ing is employed to coarse-grain the data into representative states,

and the temporal evolution is modeled as a probabilistic Markov

model. In this model, the state vector of cluster probabilities con-

verges to a fixed value representing the post-transient attractor,

which leads to the disappearance of the dynamics. To surmount

the drawback of the CMM in modeling dynamic evolution, we

proposed a cluster-based network model (CNM) in Refs. 54 and

55, which is a universal method for data-driven modeling of non-

linear dynamics from time-resolved snapshots without any prior

knowledge. The CNM, a branch of CROM, combines a cluster

analysis of an ensemble observations from simulations or experi-

ments and a network model for transitions between the different

flow characteristics educed from the cluster analysis. The snap-

shots are coarse-grained into a small number of clusters with the

representative states called centroids. The dynamics is described

by a network model with “constant velocity flights” between the

centroids.

This study seeks to understand the flow dynamics of a blunt

basesupersonicmixinglayerusing theCNM insteadofthestatistical

analysis.Consideringvariousconfigurations andinfluence factorsin

this kind of flow, this paper is studied based on two given examples.

This paper is organized as follows: Sec. II briefly describes the phys-

ical flow model and large eddy simulation method. The grid inde-

pendence analysis and code verification are carefully carried out to

ensure the accuracy of numerical results. Section III introduces the

algorithmofthecluster-basednetworkmodelstep bystep, including

the pre- and post-processing for the clustering results as well as the

method for creating a network model based on them. To give read-

ers an overall expression, the fundamental flow features are briefly

describedinSec. IV A. SectionIVBcomparestheCNM results from

the two cases in order to provide an in-depth insight into the blunt

base supersonic mixing layer. Finally, Sec. Vsummarizes some key

findings.

II. NUMERICAL SIMULATION SETUP

A. Physical model

The physical model in the present simulation is the supersonic

mixing layer downstream a splitter plate. We include the splitter

plate in the computational domain to reproduce realistic situation

as in our previous experiment.57 Dual parallel streams encounter in

the trailing edge of the splitter plate. Interactions happen and grad-

ually generate a mixing layer whose visualization thickness increases

along with the streamwise orientation. Figure 1 illustrates the two-

dimensional numerical sketch. The flow is described in a Cartesian

coordinate where the origin is located in the geometrical center of

the starting edge of the splitter plate. The xaxis is aligned with

the streamwise direction and yaxis with the transverse direction.

U,P, and Trepresent the inlet velocity, static pressure, and static

temperature, respectively. The subscripts “1” and “2” denote the

upper and lower streams, respectively. In this study, two splitter

AIP Advances 11, 085322 (2021); doi: 10.1063/5.0062145 11, 085322-2

AIP Advances ARTICLE scitation.org/journal/adv

FIG.1. The two-dimensionalsketchof the numericalsimulation setup. The numer-

ical domain is rectangular with the dimension excluding a splitter plate. The letters

“A,” “B,” and “C” are used to denote the different boundaries. Two-dimensional

simulation is acceptable at Mc =0.22.

plates are considered in the numerical simulations: the thin split-

ter plate (“denoted by TN”) and thick splitter plate (denoted by

“TK”). Note that only two-dimensional simulations are conducted

since the two-dimensional vortex structures that originated from

the Kelvin–Helmholtz (K–H) instability play the absolutely domi-

nant role in the weakly compressible regime with Mc <0.4.25,56 The

mass transfer and combustion in the mixing layer are not involved

in this simulation. The two splitter plates (TN and TK) have the

same dimension in the length Ls=80 mm but different thicknesses:

hTN =1 mm and hTK =5 mm. The convective Mach number is

expressed by

Mc =u1−u2

a1+a2(1)

with the corresponding convective velocity

Uc=U1a2+U2a1

a1+a2. (2)

Note that ais the local Mach number. The convective velocity Ucis

the averaged moving downstream velocity of the characteristic flow

structures. The Reynolds number is based on the convective veloc-

ity Ucand thickness of the splitter plate h. They read ReTN =3.46

×104and ReTK =1.73 ×105. The computational domain extends a

rectangular area excluding the splitter plate,

Ω=Lx×H−Ls×h. (3)

Many previous literature studies have demonstrated that the

splitter plate is of great significance to the vortex dynamics of the

mixing layer.33–35 Herein, we stress the differences between this

studyand other numerical simulations of the supersonic mixing lay-

ers. A tanh-type velocity profile is generally employed as the inlet

condition for a rectangular domain. This estimation is convenient

for the fundamental analysis on the mixing layer vortex but stay

away from the engineering practice.

The dimensions of the computational domain and flow param-

eters in this numerical study are inspired from our previous numer-

ical6,7 and experimental study.57 Notably, the TK case is the same as

our experimental study in the supersonic wind tunnel.6The TN case

is benchmarked as a baseline. The parameters of dual streams are

thoroughly listed in Table I.

TABLE I. The parameters of dual incoming streams for numerical simulation.

Upper Lower

Convective Mach number Mc 0.22

Mach number a1,a22.12 3.18

Convective velocity Uc(m/s) 576.9

Density ρ(kg/m3) 0.0556 0.0834

Density ratio s=ρ2/ρ11.5

Velocity ratio r=U2/U11.2

Streamwise velocity U(m/s) 519 623

Static pressure Ps(Pa) 2 640 2 520

Total pressure P0(Pa) 21 360 101 325

Boundary layer thickness δ(mm) 0.5 0.5

We have to state that many vital variables will affect the super-

sonic mixing layer’s evolution in the downstream. These variables

include the thickness and geometrical shape of the splitter plate, the

thickness of the initial boundary layer, the Reynolds number, and

theconvectiveMachnumber.This paperdoesnotdiscussthesevari-

ables and focuses on the dynamics using the cluster-based network

model introduced in Sec. III.

B. Numerical method

In our present simulation, the compressible mixing layer is

computed using an in-house large eddy simulation (LES) solver

based on the finite volume method (FVM). In this solver, the mes-

sage passing interface (MPI) based on parallel architecture using

domain decomposition and distributed memory is employed to

gain high calculation efficiency. This solver has been successfully

applied to the simulations of the supersonic mixing layer and cavity

flow.6,7,51

The two-dimensional Navier–Stokes equations are filtered

using the Deardorff box filtering58 and resulting equations are as

follows:

Continuity equation,

∂¯

∂t+∂

∂xj(¯

ρ˜

uj)=0. (4)

Momentum equation,

∂

∂t(¯

ρ˜

ui)+∂

∂xj(¯

ρ˜

ui˜

uj+¯

pδij −˜

τij +τsgs

ij )=0. (5)

Energy equation,

∂

∂t(¯

ρE)+∂

∂xj[(¯

ρE+¯

p)˜

uj+˜

qj−˜

ui˜

τij +Hsgs

j+σsgs

j]=0. (6)

Here, ρdenotes the density. ui(i=1,2)is the velocity vector in

Cartesian coordinates and pis the pressure. The total energy of the

systemisthesumoftheinternalenergyandkineticenergy.Then,the

filtered total energy is determined by the sum of the filtered internal

energy, the resolved kinetic energy, and the subgrid kinetic energy.

The subgrid-scale stress resulting from the filtering operations is

modeled by the subgrid-scale turbulences based on the Boussinesq

AIP Advances 11, 085322 (2021); doi: 10.1063/5.0062145 11, 085322-3

AIP Advances ARTICLE scitation.org/journal/adv

hypothesis59 and can be computed by

τsgs

ij −1

3τkkδij =−2μtSij, (7)

where μtis the subgrid-scale turbulent viscosity. The Sij is the rate-

of-strain tensor for the resolved scale defined by

Sij =1

2(∂ui

∂xj+∂uj

∂xi). (8)

Note that the isotropic part of the subgrid-scale stress τkk is added

to the filtered pressure term. The compressible form of the subgrid

stress tensor is defined as

τsgs

ij =¯

ρ˜

uiuj−¯

ρ˜

ui˜

uj. (9)

The term can be split into its isotropic and deviatoric parts,

τsgs

ij =τsgs

ij −1

3τkkδij

´¹¹¹¹¹¹¹¹¹¹¹¹¹¹¹¹¹¹¹¹¹¹¹¹¹¹¸¹¹¹¹¹¹¹¹¹¹¹¹¹¹¹¹¹¹¹¹¹¹¹¹¹¹¶

deviatoric +1

3τkkδij

´¹¹¹¹¹¹¸¹¹¹¹¹¹¶

isotropic

. (10)

The deviatoric part of the subgrid-scale stress tensor is modeled

using the compressible form of the Smagorinsky model,60

τsgs

ij −1

3τkkδij =−2μt(Sij −1

3Skkδij), (11)

where the SGS eddy viscosity is

μt=ρL2

S∣¯

S∣. (12)

In the above equation,

∣¯

S∣≡√2¯

Sij ¯

Sij. (13)

The mixing length for subgrid scales Lsreads

Ls=min(κd,CsΔ). (14)

In Eq. (14),κis Karman’s constant, dis the distance to the clos-

est wall, and Csis the Smagorinsky constant. Δis the scale of the

local grid scale and can be computed according to the volume of the

computational cell using

Δ=V1/3, (15)

and ∣¯

S∣is the modulus of rate-of-strain for the resolved scales.

Following Ref. 6,Csis set as 0.1 for the reliable results, which

is also found to yield the best results for a wide range of

turbulence.

The time-marching is performed by means of a second-

order, implicit dual time step scheme to improve the calculat-

ing efficiency. A Lower-Upper Symmetric Gauss–Seidel (LU-SGS)

method is employed to achieve the inner product. The gradients

are calculated by the Green–Gauss cell-based method. The second-

order spatially accurate upwind scheme with the low-diffusion Roe

Flux Difference Splitting (low Roe-FDS) formulation is utilized in

this solver to reduce the dissipation during calculations, which

is an updated version of the Roe-FDS. The time step is fixed at

Δt=5×10−8s to facilitate the clustering snapshots at convenience.

The Courant–Friedrichs–Lewy (CFL) number is less than 0.5. Note

that the physical time step Δtis limited only by the level of desired

temporal accuracy. The pseudo-time step Δτin the dual time

step scheme is determined by the CFL condition of the implicit

time-marching scheme. The calculation case in this paper simu-

lates nearly 2 ×105time steps to ensure accuracy of the statistical

results.

To simulate the actual situation of the boundary layer, the

velocity conditions of two inlets with a Blasius flat-plate boundary

layer spanning the inlet were determined. The two inlets (denoted

by “A” in Fig. 1) are tested by imposing random noise on a mean.

The mean profile has a 1/7th power law implemented, and a tur-

bulent kinetic energy of 0.4% is employed. The non-slip boundary

and constant temperature conditions (denoted by “B” in Fig. 1)

are adopted in the boundaries of the upper wall, bottom side-

wall, and splitter plate. The outflow condition (denoted by “C” in

Fig. 1) is set as the non-reflective condition61 for calculation sta-

bility. The reflection waves from the outflow boundary are well

controlled.

The supersonic mixing layer with the weak compressibility

effect (Mc ≤0.4, Mc =0.22 in this paper) will mainly exhibit two-

dimensional vortex structures.8In addition, the mixing layer with

these conditions will exhibit three-dimensional behaviors in a 3D

view, but the two-dimensional vortex tubes play the dominant role.

Readers are referred to our published results with nearly the same

conditions in Ref. 7.

C. Grid independence analysis and code verification

To ensure the simulation’s validity and reliability, the grid

independence analysis is conducted using three distinct grid sys-

tems, namely, coarse, moderate, and refined grids. Table II con-

tains the detailed grid information, including the minimum scale of

Δxand Δyrequired to satisfy the constraints of x+≤1 and y+≤1.

The calculation area is discretized into a structured non-uniform

grid.By comparing the statistical profiles at the same spatial posi-

tion, the effect of grid scales on the numerical results is inves-

tigated. As illustrated in Fig. 2, the dimensionless time-averaged

temperature and pressure for the TK case at x/H=3 are com-

pared using coarse, moderate, and fine grids, respectively. Three

curves for both quantities basically coincide, indicating that no

obvious difference exists between numerical results from three

grid scales. Additionally, the results from moderate and refined

TABLE II. Three sets of meshes for the 2D numerical simulation of the supersonic

mixing layer.

Grid Coarse Moderate Refined

nx×ny838 ×260 1001 ×301 1201 ×351

Cell number 193980 272301 375361

Δymin (m) 2.00 ×10−51.00 ×10−51.00 ×10−5

Δymax (m) 5.00 ×10−45.00 ×10−45.00 ×10−4

Δxmin (m) 2.00 ×10−51.00 ×10−51.00 ×10−5

Δxmax (m) 5.00 ×10−45.00 ×10−42.00 ×10−4

AIP Advances 11, 085322 (2021); doi: 10.1063/5.0062145 11, 085322-4

AIP Advances ARTICLE scitation.org/journal/adv

FIG. 2. Comparison of time-averaged profiles in three grid systems at x/H=3

for the TK case. (a) Time-averaged static pressure and (b) time-averaged static

temperature.

meshes are highly comparable. As a result, selecting a moderate

grid scale for the computational grid can result in cost and time

savings.

We conduct code verification to ensure the solver’s reliabil-

ity. As shown in Fig. 3, the statistical results obtained using the

LES solver for the TK case agree reasonably well with the exper-

imental data. Nota bene, the experimental data are taken from

Ref. 6.

In summary, this LES solver performs admirably well when

used to simulate the supersonic mixing layer numerically.

FIG. 3. Streamwise velocity profiles obtained from the numerical results (black

solid line) and experimental data (red dots) at (a) x/H=2.5 and (b) x/H=3.5 for

the TK case.

III. CLUSTER-BASED NETWORK MODEL

The cluster-based network model paves a new avenue to model

the high-dimensional nonlinear dynamics in a data-driven manner.

Figure 4 shows the technical sketch of the cluster-based network

model exemplified by the incompressible mixing layer. Cluster anal-

ysis groups the high-dimensional state space into a given number

of subspaces with the representative states called centroids in an

unsupervised manner. Obviously, cluster analysis is indeed one of

the dimension reduction techniques. Then, the dynamics is repro-

duced by the CNM on a directed network, where the nodes are the

AIP Advances 11, 085322 (2021); doi: 10.1063/5.0062145 11, 085322-5

AIP Advances ARTICLE scitation.org/journal/adv

FIG. 4. The methodology sketch of the cluster-based network model. The snapshots u(x,t)are collected at constant sampling frequency. These snapshots are grouped

into Kclusters in terms of their similarity meanwhile obtaining centroids ckand cluster affiliations. The representative states, i.e., centroids, are considered as the network

nodes. The direct transition probabilities between network nodes are identified from the training data, which describe all possible transitions from the probabilistic point of

view. Finally, the complicated dynamics in the high-dimensional state space is converted into a simple dynamics between centroids.

centroids, i.e., the representative states for the subspaces of the sys-

tem. The transition probability and transition time between nodes

can be inferred from the original data. Generally, there are three

steps for the CNM analyzing a dynamical system. The data should

be collected at first. These data will be coarse-grained into a given

number of clusters to get the centroids and time series of labels for

each observation. The network is built temporally by means of the

timeseries oflabelsand spatiallybymeans ofthecentroids. Network

science becomes increasingly popular in nearly all fields of mathe-

matical modeling, including computer sciences, biology, and fluid

dynamics.62–64

A. Data collection

The CNM is a purely data-driven model to resolve a dynamic

system without any prior knowledge. The starting point of the CNM

is the data collection of Mtime-discrete N-dimensional state of

the system. Theoretically, any variables can be fed into the CNM

due to the great robustness to the nonlinear system. Recommended

by the model developers, the interested variables, especially reflect-

ing the property of a dynamic system, should be first picked up.

The state represented by the collected data can comprise the full

state, the low-dimensional representation of the full state, or an

observation. Note that these data can be from the numerical sim-

ulations or experiments. In addition, the model quality is directly

related to its training data. In particular, the sampling frequency for

snapshots must be selected that the relevant dynamics are resolved.

The total time range covering these snapshots should ensure the

statistic convergence of the dynamics. Therefore, the CNM needs

a larger amount of training data than other data-driven models

like DMD. In the present study, we collect Mtime-resolved veloc-

ity snapshots in a steady domain Ωas the observation for the

CNM. These snapshots are obtained from the numerical simula-

tion described in Sec. II. The velocity field is collected at an instant

time step Δtcand denoted as um(x)∶=u(x,tm),m=1,...,M. The

time step Δtcshould be short enough to capture the character-

istic flow features. We emphasize the difference between Δtin

LES and Δtcthat two variables should satisfy the requirement

of Δt≤Δtc.

B. Cluster analysis

1. Data preprocessing: Compression using POD

We consider the POD method as a data compression tool to

reduce the computational loads, but this step is not mandatory for

the CNM. The data compression is highly recommended for the

clusteringoflarge high-dimensional datasets sinceitispretty expen-

sive.Lumley65 firstintroducedPODinto theturbulencecommunity.

This method aims to find optimal basis functions by maximizing

the mean square projection of the system using a spatial two-point

correlation in an eigenvalue problem. The inner product of two

square-integrable velocity fields u1and u2∈L2(Ω)is given as

(u1,u2)Ω∶=∫Ωu1(x)⋅u2(x)dx. (16)

The corresponding norm reads ∥u∥Ω∶=√(u,u)Ω. All collected

snapshots are averaged in time to get the mean flow um. The

snapshot-based POD is applied to the velocity fluctuations,

vm(x)=u(x,tm)−um(x). (17)

The two-point correlation matrix C=(Cmn)∈RM×Mof the fluctua-

tions is first determined,

Cmn =1

M(vm,vn)Ω. (18)

The eigenvalue problem for the correlation matrix Cis solved to

determine the POD modes and time coefficients (i.e., mode ampli-

tudes),

Cei=λiei,i=1,...,M. (19)

Each POD mode is a linear combination of the snapshot

fluctuations,

ui=1

√Mλi

∑

m=1em

ivm,i=1,...,N. (20)

AIP Advances 11, 085322 (2021); doi: 10.1063/5.0062145 11, 085322-6

AIP Advances ARTICLE scitation.org/journal/adv

Note that any two POD modes are mutually orthonormal, i.e.,

(ui,uj)Ω=δij,∀i,j∈{1,...,M}. The POD mode amplitudes can be

expressed by

i∶=√λiMei

m,i=1,...,N. (21)

Finally, the POD method decomposes each snapshot into a lin-

ear combination of the time coefficients ai(t)and spatial modes ui,

u(x,t)=um(x)+N

∑

i=1ai(t)ui(x), (22)

where N=M−1 means the lossless POD decomposition.

Cluster analysis identifies the flow pattern via the similarity

metric between any two snapshots. The similarity is measured by

the Euclidean distance, expressed by

d(um,un)=√∫Ωdx∥um−un∥2. (23)

In this equation, the whole computational domain Ωis used for the

integration. In the cluster algorithm, many repetitions are required

to gain a better result, which leads to extremely expensive calcula-

tion. Note that any two POD modes are mutually orthonormal. The

distance between any two velocity fields can be converted into the

distance of the corresponding POD mode amplitudes, given as

d(um,un)=∥um−un∥Ω=∥am−an∥. (24)

In this way, a large amount of the area/volume integrals in the

computation of Euclidean distance for cluster analysis are avoided

since the most expensive process has been finished by comput-

ing the two-point correlation matrix C18 in the POD method.

Therefore, each snapshot is represented by a time coefficient vector

am=[a1(tm),a2(tm),...,aM−1(tm)]. The dataset of velocity snap-

shots is converted into a large matrix Acontaining a series of POD

time coefficient vectors,

A=⎛

⎜

⎝

a1(t1)⋅⋅⋅ aM−1(t1)

⋮ ⋮ ⋮

a1(tM)⋅⋅⋅ aM−1(tM)⎞

⎟

⎠. (25)

The subsequent clustering and analysis in this paper are con-

ducted based on this matrix. In Sec. III B 2, the cluster analysis

algorithm is described based on the POD time coefficients.

2. Clustering analysis to identify characteristic

flow patterns

The cluster analysis is a dimension reduction tool that dis-

cretizes the high-dimensional state space in an unsupervised man-

ner. To be specific, it partitions the observations/snapshots that are

mutually similar to Krelatively homogeneous clusters. The homo-

geneity of each cluster is quantified by the cluster diameter Dk.

The representative state of each cluster is described by the centroid

ck,k=1,...,K. The similarity of any two velocity fields umand un

is measured using Euclidean distance D. As mentioned in Eq. (24),

the distance of two velocity snapshots umand unbased on the norm

associated with the Hilbert space ℒ2(Ω)of square-integrable func-

tions is converted to the distance between two corresponding POD

time coefficients amand an. Therefore, the cluster analysis is applied

to a sequence of POD time coefficients am=am(tm),m=1,...,M.

In this study, cluster analysis is achieved using the unsupervised k-

means++algorithm66 that minimizes the inner-cluster variance and

maximizes the inter-cluster variance.

In the k-means++algorithm, a set of centroids ckis ran-

domly initialized and their performance can be evaluated by the

total inner-cluster variance Jof snapshots amwith respect to their

corresponding centroids,

J(c1,...,ck)=K

∑

k=1∑

am∈𝒞k

d(am,ck). (26)

Thek-means++algorithmsolvesthisoptimizationproblemthrough

iterations to pursue the minimum inner-cluster variance. The itera-

tions will stop until the convergence is reached, and meanwhile, a

set of optimal centroids copt,a

k(k=1,...,K)associated with the POD

time coefficients are determined,

(copt,a

1,...,copt,a

K)=argmin

c1,...,cK

V(ca

1,...,ca

K). (27)

Thecluster analysis assignsacluster indextoeach snapshot thatcor-

responds to the nearest centroid, resulting in the cluster affiliation

function, k(am)=argmin

i∥am−ca

i∥. (28)

This function defines Kcluster regions as Voronoi cells 𝒞ksur-

rounding centroids. In order to communicate with other following

variables conveniently, a characteristic function is proposed as an

alternative of Eq. (28) and reads

k∶=⎧

⎪

⎨

⎪

⎩

1 ifam∈𝒞k,

0 otherwise. (29)

The centroid is defined as the mean flow of all snapshots

belonging to a cluster,

ck=1

nk∑

am∈𝒞k

am=1

∑

m=1Tm

kam, (30)

where nkis the number of snapshots in cluster 𝒞kand reads

nk=M

∑

m=1Tm

k. (31)

3. Data postprocessing after cluster analysis

To visualize the centroids, the flow structures are reconstructed

via the inverse process of the POD by computing the linear combi-

nationof the PODcoefficientvector of centroidac

i,kandPOD modes

ui, i.e.,

ck(x)=um(x)+N

∑

i=1ac

i,kui(x). (32)

Similar to POD and DMD, the transverse fluctuation velocity v′is

visualized to show the flow patterns rather than including the mean

flow.

AIP Advances 11, 085322 (2021); doi: 10.1063/5.0062145 11, 085322-7

AIP Advances ARTICLE scitation.org/journal/adv

The performance of cluster analysis is evaluated by the degree

of homogeneity of the resolved subspaces. The diameter of a clus-

ter is the distance between the cluster’s two most distant observa-

tions. The cluster diameter is the maximum Euclidean distance [see

Eq. (23)] between any two velocity fields amand anwithin a cluster,

and it is expressed as

Dk=max

m,n(Dmn :am,an∈𝒞k). (33)

The standard deviation is a measure of how closely data are

clustered around the mean. Basically, low standard deviation means

that data are clustered around the mean, whereas a high standard

deviation indicates dispersed data. The standard deviation is defined

Rk=¿

À1

nk∑

am∈𝒞k∥am−ac

k∥2. (34)

The dissimilarity between any two representative states can be

quantified by the Euclidean distance between their two correspond-

ing centroids in the steady domain Ω. This leads to a K×Ksym-

metric distance matrix d. Each entry dij of this matrix is defined as

dij =d(ci,cj)=√∫Ωdx∥ci−cj∥2. (35)

The clusters resolve the regions of the dynamical trajectory

associated with various turbulent kinetic energies (TKEs). Herein,

we estimate the mean TKE of each cluster using the snapshots in

this section,

TKEk=∑M

m=1Tm

k∥am−a∥2

2∑M

m=1Tm

. (36)

The probabilistic representation of a dynamical system is

defined as

qk∶=nk

M=1

∑

m=1Tm

k. (37)

A probability qkis approximated in this equation by the ratio of

snapshots in cluster kto the total snapshots M. This parameter

contains data about the cluster population.

C. Network model for state transitions

In Sec. III B 2, the k-means++algorithm assigns a cluster index

to each snapshot from the training dataset, which leads to a discrete-

time cluster affiliation function k(am). From a probabilistic stand-

point, the link between any two clusters is described using a derived

network model. The cluster affiliation function is used to infer the

direct transition matrix Qand its related time scale matrix Tin

the network. The network model, unlike the cluster-based Markov

model,48,51 ignores inner-cluster residency and only considers inter-

cluster switches. The element of Qij and Tij describes an event in

which the cluster jtakes time Tij to transit to cluster iwith the prob-

ability Qij. The direct transition probability Qij is inferred from the

training data as

Qij =nij

nj,i,j=1,...,K, (38)

where nij is the number of transitions from 𝒞jto 𝒞iand njis the

number of transitions departing from 𝒞jregardless of the terminal

points. Note that the key diagonals of the matrix Tare all 0 as the

lingering in the same cluster is omitted.

A time-delay gained from the training data is embedded into

the low-dimensional dynamical system, in contrast to the constant

time step of the cluster-based Markov model.48,51 Assume tn(tn+1)

is the time of the first (last) snapshot to enter (leave) the cluster 𝒞i.

The residence time in cluster 𝒞iis computed using the following

formula: τi=tn+1−tn. (39)

The spatial mode in the low-dimensional dynamical system is

represented by the centroid ci, which represents the spatial mode in

the low-dimensional dynamical system, supposes to stay at the cen-

ter of this time range (tn+tn+1)/2. The individual transition time

from cluster jto cluster iis then defined as half the residence time of

both clusters and reads

τij =τi+τj

2=tn+2−tn

2. (40)

Because of the dynamical system’s complexity, the transition tra-

jectory between two clusters may change, resulting in distinct time

scales. To characterize time scale, we average the transition dura-

tions between clusters, imitating the definition of a centroid, which

is the mean of all snapshots belonging to this cluster. The entry of

the averaged transition time matrix reads

Tij =⟨τij⟩. (41)

The cluster-based network model can also predict the future

evolution of a dynamical system with a predetermined initial state.

Thevalidationofthemodelisofsignificancetogainaccurateperfor-

mance. We do not introduce the detailed algorithm in this paper to

focus on the physical understanding of the supersonic mixing layer.

Readers are referred to Ref. 54 for the details.

IV. RESULTS AND DISCUSSION

In this section, the TN and TK cases are simulated using

the large eddy simulation. The dynamics is thoroughly investi-

gated by comparing CNM output for both TN and TK cases.

Section IV A provides an overview of the flow characteristics

in order to gain a fundamental understanding. The supersonic

mixing layer is characterized by a large range of time and spa-

tial scales, intrinsic high dimensionality, and nonlinear dynamics.

Section IV B presents the clustering-derived representative spatial

modes (i.e., centroids). Cluster analysis is an entirely data-driven

dimension reduction technique for identifying flow patterns. The

supersonic mixing layer’s high-dimensional state space is projected

into K=10 centroids. In Sec. IV C, the cluster-based network

model is used to reproduce the high-dimensional nonlinear com-

plex dynamical system of the supersonic mixing layer. To compre-

hend the dynamics, the state temporal transitions are thoroughly

analyzed.

A. Flow features

Numerous publications attest to the significance of the role of

the splitter plate in the development of the supersonic mixing layer.

The fundamental flow characteristics of the blunt-base supersonic

mixing layer are visualized in Fig. 5 using the numerical schlieren

AIP Advances 11, 085322 (2021); doi: 10.1063/5.0062145 11, 085322-8

AIP Advances ARTICLE scitation.org/journal/adv

FIG. 5. The fundamental features of the supersonic mixing layer visualized via numerical schlieren67 in the TK case.

method.67 This sketch clearly illustrates the layout of wake-mixing

layer and shock-mixing layer interactions, which differs from the

layout observed in the majority of simulations using a tanh-type

velocity profile. Along with the expansion waves, the supersonic

streams begin to separate at the trailing edge of the splitter plate.

Instead of the immediate confrontation, two streams are separated

at the rear of the splitter plate by a low-speed wake. When two shear

layers approach the splitter plate’s centerline, compression occurs

as the flow attempts to realign with the streamwise direction. The

compressionwavesoriginated intheshearlayer’ssupersonicregions

graduallycoalesce into therecompressionoblique shockwave.Their

reflected waves from the upper or lower wall interact with the rede-

veloping mixing layer and exert influence on the growth process

via the modification of compressibility effects along the streamwise

direction. Following the disappearance of the wake, the two main-

streams merge into a redeveloping mixing layer. As illustrated in

Fig. 5, roll-up typical structures, i.e., K–H vortices, develop down-

stream of the flow field due to the free mixing layer’s inherent

K–H instability. We assume that the wake instability has a signif-

icant effect on the expansion of the mixing layer, which will be

analyzed in Sec. IV B. Maintaining K–H instability results in the

pairing and merging process of adjacent vortices, which leads to

the growth of the mixing layer. The curvature shocklets are cap-

tured astonishingly well in the vicinity of the large-scale structures

in Fig. 5.

FIG. 6. Instantaneous vorticity contour for the two cases.

The instantaneous vorticity fields for these two cases are com-

pared in Fig. 6 to demonstrate the dynamic difference. Consider-

ing the intermittent dynamics of the wake-mixing layer, these two

instantaneous images can only provide a primary understanding of

flow features. These two examples demonstrate significant differ-

ences in wake size, vortex shedding position, and vortex size. As for

the TN case, the wake-mixing layer exhibits a rapid roll-up of K–H

vortices following a small wake. This indicates that the wake has a

negligible effect on the mixing layer dynamics. In comparison, the

shedding positions in the TK case are significantly delayed down-

stream, but larger-scale structures emerge at the start of the redevel-

oping mixing layer. Note that the occurrence position of shedding

is not fixed and generally shows intermittent dynamics. This is fur-

ther confirmed in the centroids captured in Sec. IV B by the cluster

analysis.

B. Spatial modes: Characteristic flow states

To implement the cluster analysis, M=2000 post-transient

velocity snapshots are collected with an instant time step of Δτ

=5Δt=2.5 ×10−7s. The sampling frequency 1/Δτis sufficiently

small to capture the relevant dynamics, and the time range of

M=2000 snapshots is satisfactory for the accurate statistics. The

optimal number of clusters Kis a trade-off between the resolu-

tion of dynamics and modeling time horizon. These snapshots are

first compressed using the POD method to withstand the compu-

tational loads generated by the k-means++algorithm’s iterations.

These snapshots, denoted by POD time coefficient vectors, are par-

titioned into K=10 clusters in this paper. The number of clusters

K=10 is confirmed to be reasonable for the current dataset, consis-

tentwith ourpreviouspublications inRefs.7,51, and54.We remark

on the critical nature of model validation for the sake of model

accuracy. The validation of the CNM via the correlation function

is omitted here in order to focus on the mechanism understanding

and also for the writing conciseness. Readers are referred to Refs. 7,

51, and 54 for more details.

The homogeneity of each cluster is quantified using the cluster

diameter Dkand standard deviation Rk, as shown in Fig. 7. Similar

diameter and standard deviation distributions in a case demonstrate

the k-means++algorithm’s superior performance, as each clus-

ter is relatively homogeneous in size. The M=2000 snapshots are

AIP Advances 11, 085322 (2021); doi: 10.1063/5.0062145 11, 085322-9

AIP Advances ARTICLE scitation.org/journal/adv

FIG. 7. Cluster diameter Dk(a) and standard deviation Rk(b) of K=10 clusters for TK (gray) and TN (white).

sufficient to fill the high-dimensional state space of the supersonic

mixing layer. The clustering process discretizes high-dimensional

state space into a set of Ksubspaces. Clustering projects the high-

dimensional state space into K-dimensional subspace. The diameter

quantifiesthe sizeofa subspace andstandarddeviations describethe

distributeddensity of observations(snapshots)surroundingthe cor-

responding centroids. Each cluster’s Dkand Rkare nearly identical

in the TN case, indicating a relatively homogeneous partition of the

high-dimensional state space and also implying periodic dynamics.

Almost all diameters of TK are more than twice as large as those of

the TN cases. The standard deviations of TK (at least k=1,...,8)

are over three times larger than that of the TK case. This implies

that,abstractlyspeaking,theTKcasespansasignificantlylargerstate

space than the TN case with dispersedly distributed snapshots.

The K=10 centroids can be used to illustrate the representa-

tiveflowpatternsof clusters, as shown in Figs. 8and9.Centroidsare

FIG. 8. K=10 centroids in the whole computational domain for the TN case. Only the transverse fluctuation velocity v′is visualized.

AIP Advances 11, 085322 (2021); doi: 10.1063/5.0062145 11, 085322-10

AIP Advances ARTICLE scitation.org/journal/adv

FIG. 9. K=10 centroids in the whole computational domain for the TK case. Only the transverse fluctuation velocity v′is visualized.

reconstructed using the POD method’s reverse process, as described

inEq. (32). Each centroid has two components in this study: stream-

wise velocity and transverse velocity. We visualize the transverse

fluctuation velocity v′in this section as it is critical for mass mixing

and approaches the vortex structures. We will refer to the transverse

fluctuation velocity v′as centroids in the following. These two cases

exhibit radically different flow characteristics. In the TN case, the

redeveloping mixing layer forms immediately after the trailing edge

of the splitter plate. The wake is small enough that its effect on the

development of the mixing layer can be completely ignored. This

resultsina significantly greater advance ofthereattachmentpointin

theTNcasethanintheTKcaseasillustratedinFig.5.Theexpansion

angle βbetween the upper and lower expansion waves, as shown in

Fig. 8, is pretty small. The spatial structure in the TN case is approx-

imately uniform structures in the streamwise direction. Starting at

x/H=1, the area occupied by the vortex structures is amplified and

gradually decreases until x/H=1.5. Another similar process occurs

betweenx/H=1.5andx/H=2.25andalsobetweenx/H=2.25and

x/H=3.

Increasing the thickness of the splitter plate greatly affects

the spatial and temporal dynamics of the supersonic mixing layer.

The TK case contains a variety of structures of varying sizes and

shapes. The wake behind the trailing edge has a very low transverse

fluctuation velocity v′where the streamwise motions dominate the

dynamics of this region. The shapes of the K=10 centroids reveal

two distinct subsets: Kelvin–Helmholtz vortices (K–H k=1...8)

andvortexpairing (VP, k=9,10). It is widely acknowledgedthatthe

VP behavior significantly generates more energy than K–H behav-

ior. In the K–H subset, the wake effects cease to exist at around

x/H=1.5, followed by the fast expansion of the mixing area to

approximately x/H=2 with the peak. Similarly, after the peak, the

mixing area decreases to around x/H=2.5. Subsequently, another

increase begins. Various flow behaviors, such as vortex pairing,

merging, and tearing, are well captured in the VP subset. In addi-

tion, rare events, such as quadruple-vortex pairing/merging, are also

clearlyresolvedincentroidsk=9,10,asdenoted by “QP”and“QM”

(QP′and QM′). Due to the wide value range of the legend, the

wake flow near the trailing edge of the splitter plate is deliberately

eliminated in Fig. 9. This area has been magnified to show the fine

flow structures of the wake, as illustrated in Fig. 10. The absolute

v′amplitudes are only equal to 1 as the maximum ΔUapproaches

100. The reattachment points in these centroids are not fixed, which

are related to the redeveloping mixing layer’s diverse behaviors.

The shock-wave/mixing layer interactions are well captured in the

majority of centroids, denoted by “S”inFig. 9. In comparison to

the VP subset, the reattachment points of the K–H subset are gen-

erally delayed, except for centroid k=1. The unstable position of

the reattachment point strengthens the initial turbulence of the

AIP Advances 11, 085322 (2021); doi: 10.1063/5.0062145 11, 085322-11

AIP Advances ARTICLE scitation.org/journal/adv

FIG. 10. K=10 centroids near the trailing edge of the splitter plate for the TK case. Only the transverse fluctuation velocity v′is visualized. These figures corresponds to

Fig. 9 with a new legend. The star symbol (⋆) and solid lines denote the reattachment points and reattachment shock waves, respectively. The intersection location of two

shock waves approximates the reattachment point. We emphasize that the denoted positions are roughly determined according to the visualization.

redeveloping mixing layer, implying the presence of a variety of

dynamical behaviors downstream.

The geometrical difference of the flow structures in the same

regime can be overall measured using the Euclidean distance

d(ci,cj), which results in the symmetry distance matrix D∈R10×10

[see Eq. (35)]. Figure 11 depicts the symmetrical distance matrix

Dfor TN and TK, respectively. The cluster distance matrix can be

usedto estimatethetrajectory lengthforcluster transitionsaswell as

the geometric relationship between the centroids. The TN’s distance

matrix reveals the close association of two adjacent centroids, for

example,centroidcnismutuallysimilartocentroidcn+1.Thisindeed

implies that the flow structures in the convective process move in a

periodic manner. The adjacent centroids exhibit a small phase dif-

ferencealmost entirelyindependent ofstructureshape modification.

The TK’s distance matrix is more informative than the TN’s since it

further highlights the difference between two regimes as mentioned

above: the K–H and VP.

C. Dynamics analysis: State temporal transitions

The CNM extracts a sequence of events from a probabilistic

point of view, providing an in-depth insight into the flow dynamics

of the supersonic mixing layer, especially the intermittent behavior.

After cluster analysis, the k-means algorithm assigns an affiliation

FIG. 11. Cluster distance matrix D[see Eq. (35)]. (a) TN and (b) TK. Distances are linearly scaled with all zero elements in the diagonal.

AIP Advances 11, 085322 (2021); doi: 10.1063/5.0062145 11, 085322-12

AIP Advances ARTICLE scitation.org/journal/adv

FIG. 12. Cluster affiliation function (a) TK and (b) TN.

to each snapshot, which is the index of the corresponding near-

est centroid of the snapshot. As illustrated in Fig. 12, the cluster

affiliation function depicts the cluster temporal transition among

training data, i.e., these M=2000 snapshots. We establish a rule

that the first cluster to appear is designated cluster k=1. Then, the

subsequent newly appearing cluster is considered to be k+1. These

cluster transitions are summarized using a direct transition matrix

(DTM) Qwith the averaged time scale, i.e., averaged transition time

matrix T, as shown in Figs. 13 and 14, respectively. The TN case

exhibits strictly periodic dynamics with a cycled state visit at cen-

troid 1 →2,...,9 →10. As shown in Fig. 16, there is no fluctuation

fortheturbulentkinetic energyduring theperiodicdynamics,which

indicates the same contributions of each cluster to the flow mixing.

The resolved centroids move downstream continuously with a short

distance between adjacent centroids. As shown in Fig. 8, a vortex

denoted as “B” convects downstream from centroid k=1 until this

structure propagates out of the computational domain in centroid

k=10. The vortex’s geometrical modification is almost unheeded

except the phase difference in this cycled behavior. In addition, the

expansion angle denoted by “β” in Fig. 8 formed by the two expan-

sion waves is not changed in the periodic motion. The observations

belonging to each cluster are equally distributed in population with

a probability qTN

k≈0.1,k=1,...,10, as shown by white rectangles

in Fig. 15.

We must emphasize that the TN case displays pretty spe-

cial dynamics only with K–H vortices in the community of the

supersonic mixing layer. The disappearance of the vortex pair-

ing could result from the limited length of the computational

domain. In addition, the length of the splitter plate Ls, inflow con-

ditions of the numerical simulation, and Reynolds number could

also significantly affect the dynamics. We focus on the dynamics

of the given supersonic mixing layer instead of discussing these

factors.

Increasing the thickness of the splitter plate entangles the flow

dynamics of the TK case with a highly random and chaotic state. As

shown in Fig. 16, the energy of all TK’s clusters is much higher than

thatof TN.However,the populationexhibitsmaldistribution, which

isassembledinthefirstseveralclusters.TheVPregimek=9,10with

largerTKEispopulatedataprettysmallprobability.TheK–Hsubset

containsthe majority ofstatetransitionsamong TK’s clusters,which

are accompanied by a small fluctuation in turbulent kinetic energy

(seeFig.16).Thephenomenonis explained by the probability distri-

bution qk, which shows that the K–H clusters are overall populated

at 96.45%. The direct transition matrix Qcan be used to investigate

the dynamics of state transitions. The cluster k=1 is assumed to be

the initial state as this cluster is populated at the largest probabil-

ity. In the dynamics, only the most probabilistic transition route is

considered. The cluster k=1 prefers cluster k=3 as the next state.

Then, within the K–H regime, there is one cyclic motion with high

FIG. 13. Cluster-based network model for TN case. (a) Direct transition matrix Q. (b) Averaged transition time matrix T.Qand Thave the same geometry. White areas in

two figures denote the zero elements. All elements of Qare 1.

AIP Advances 11, 085322 (2021); doi: 10.1063/5.0062145 11, 085322-13

AIP Advances ARTICLE scitation.org/journal/adv

FIG. 14. Cluster-based network model for the TK case. (a) Direct transition matrix Q. (b) Averaged transition time matrix T.

probability: 3 →5→6→3. The flow starts with an advanced reat-

tachment point (denoted by a ⋆in Fig. 10), revealing fast develop-

ment of the mixing layer of dual mainstreams at around x/H=1.24

in cluster k=1 (see Fig. 10). We emphasize that the reattachment

point and reattachment shock are determined roughly according

to the visualization instead of the accurate detection method. The

numerical schlieren67,68 could be a good choice. This state takes

longertotransitintotheclusterk=3[seeFig.14(b)],wherethereat-

tachment point and redeveloping mixing layer are both delayed to

around x/H=1.30. The expansion angle is also found to decrease.

The cluster k=5 follows the trend of the two characteristic flow fea-

tures being further postponed to x/H=1.38 and then advanced to

x/H=1.3 in cluster k=3. During the intermittent evolution, the

size of the wake oscillates in a similar way. In contrast to the loop

motions in the TN case, the wake effects on the downstream mixing

layer are greatly amplified, affecting the downstream redeveloping

mixing layer. The advanced (delayed) reattachment point, as shown

in Figs. 9 and 10, indicates the complicated (simple) vortex struc-

tures in the near field but simple (complicated) vortex structures in

the far field. The turbulent kinetic energy carried by each cluster in

the cyclic motion is remarkably similar, with only minor variations.

FIG. 15. Cluster population probability qkfor the TK (gray) and TN (white) case.

Beginning with cluster k=1, the carried turbulent kinetic energy

gradually climbs to the peak with cluster k=5. In the near field, the

centroid c5exhibits a fast expansion of the K–H vortex, while in the

far field, vortex pairing/merging and vortex tearing can be observed.

With a smaller probability, the cluster k=1 could select the

cluster k=2 as the next state. Then, the following transition is

→8→5→6→3. The cluster k=2 passes through cluster k=8

to enter the route the same as that with the highest probabil-

ity. The reattachment point displays a slower downstream motion

in the cyclic motion. Meanwhile, the TKE fluctuation exhibits a

similar trend that climbs into the peak at cluster k=5 and then

decreases. In these two dynamical routes, the cluster k=5 is the

only node that must be passed. The cluster k=3 serves as the

terminal state. In the K–H regime, a delayed reattachment point

reveals more complicated structures in the far field with the larger

TKE.The cluster k=5 acts as a bifurcation point in both dynami-

cal routes. Apart from returning to cluster k=3 via cluster k=6,

the cluster k=5 has four additional evolutions, including a direct

FIG. 16. The turbulent kinetic energy of each cluster for the TK (gray) and TN

(white).

AIP Advances 11, 085322 (2021); doi: 10.1063/5.0062145 11, 085322-14

AIP Advances ARTICLE scitation.org/journal/adv

return to cluster k=1 forming a closed-loop, passing through

cluster k=4 with a lower probability, transiting to cluster k=2

followed by a new evolution, and immediately transiting to cluster

k=7. Throughout these evolutions, the reattachment point oscil-

lates forward and downward in response to the TKE fluctuation

and the stretching of the flow structures. The VP regime k=9,10

can be merely accessed by the cluster k=8 and then reverted to

the K–H regime via cluster k=5. The VP regime only accounts for

a negligible portion with 3.55%, as shown in Fig. 15. This implies

that the VP remains a rare event in the supersonic mixing layer

with a 5 mm splitter plate despite the fact that the individual cen-

troid of the VP regime contains more turbulent kinetic energy than

the K–H regime, as shown in Fig. 16. This indicates that the K–H

regime is found difficult to reach a high-energy level. The state tran-

sition from cluster k=9 to cluster k=10 is unidirectional in the

VP subset. The flow structures retain their original shapes during

this evolution, as illustrated in in Fig. 9. With rotational deforma-

tion,theQPandQMstructuresmigratedownstreamevolveintoQP′

and QM′. To quantify the overall energy fluctuation in a particular

dynamical system, we define the weight-average cluster TKE as the

overall captured energy E=∑K

k=1qkTKEk. In the TK case, the over-

all fluctuation energy (E=1.76) is 20.5% larger than that of the TN

case (E=1.46).

V. CONCLUSION

In this article, the flow dynamics of a blunt base supersonic

mixing layer is analyzed using the cluster-based network model

(CNM). The mechanism of the wake/mixing layer attractor is com-

prehensively understood in comparison to a given benchmark with

the eliminated wake effect. With only ten physical centroids, the

CNM resolves the dynamical system of the broadband supersonic

mixing layer. The direct transition matrix Qand averaged transition

time matrix Tdefine the spatial and time scales of the ten-centroid

based low-dimensional model, respectively. In the given baseline

case, the CNM identifies only one flow pattern existing in the

dynamical system: the K–H regime. The given baseline case exhibits

strictly periodic transitions between these ten K–H clusters with

a nearly constant turbulent kinetic energy. Each cluster is equally

populated with a probability qTN

k≈0.1,k=1,...,10. Increasing the

thickness of the splitter plate undergoing a wake results in highly

complicated flow behaviors from spatial and temporal points of

view. In the supersonic blunt base mixing layer, the CNM identi-

fies two distinct flow regimes: K–H and VP. The cluster diameter

Dkand standard deviation Rkreveal that the TK case spans signif-

icantly larger state space than the TN case. Each cluster of the TK

case in the given regime brings much larger TKE than that of the TN

case. In contrast to the TN case, the TK case results in vortex pair-

ing/merging and the formation of more complex flow structures. In

termsoftemporaldynamics,theflowstructuresexhibit highly inter-

mittent behaviors with small amplitude fluctuation. The CNM has

identified some probable routes for state transitions and the related

bifurcation point, which aids in the design of control strategy. The

controldesignusingtheCNMisviachangeddirecttransitionproba-

bilityQandchangedtransitiontimeT.AppendixC4inRef.54gives

a framework of the flow control strategy.

The present study gives an example to analyze the given

blunt base supersonic mixing layer using the CNM. The resolved

instantaneous behaviors are described in two-dimensional flow

fields as graphical means instead of one-dimensional signals. Addi-

tionally, the CNM provides a set of useful methods for quantifying

thesimilarity of representative states, carried energies, mode switch-

ing time scales, and metrics in high-dimensional state space. In

general, the CNM opens a novel automatable avenue for nonlinear

dynamical modeling.

Whenthesplitter plate thickness is small or when the inlet tanh

profile is embedded, the dominant instability inclines to be con-

vective in nature. The wake behind the splitter plate’s trailing edge

behaves global instability, which exacerbates the initial instability of

the redeveloping mixing layer. The asymmetric conditions in the

upward wake flow further complicate the flow dynamics. However,

increasing the thickness from 1 to 5 mm does not significantly mag-

nify the turbulent kinetics. The TKE is only raised by 20.5%, which

is much less than the cavity-actuated supersonic mixing layer.7In

practice, the thickness of the splitter plate is invariably increased

and the flow parameters are also temporally changed during engine

operation. The comparison of direct transition matrices reveals the

dynamical difference in regimes with various initial conditions, i.e.,

the thickness of the splitter plate in this study. This paper inspires

the critical nature of geometrical dimensions for splitter plate while

developing a supersonic mixer. There should exist a set of geometri-

cal parameters that are optimal for given incoming flow parameters.

Indeed, the CNM provides a means of optimizing the flow system to

meettheneedsofresearchers.Thismeritsfurtherexploration,which

the authors intend to do in the future.

ACKNOWLEDGMENTS

H.L. appreciates the National Natural Science Foundation of

China (Grant No. 91441121) and Graduate Student Research Inno-

vation Project of Hunan Province (Grant No. CX2018B027). He

gratefullyacknowledgesthe support of the ChinaScholarshipCoun-

cil (CSC) (Grant No. CSC201803170267) during his study in Tech-

nischeUniversitätBerlinandtheexcellentworkingconditionsofthe

Hermann-Föttinger-Institut. He greatly appreciates the computing

sources and technical support provided by the National Supercom-

puter Center in GuangZhou. We appreciate valuable stimulating

discussions with Professor J. G. Tan and Professor Bernd Noack.

Theauthors declare thereisnoconflict of interestregardingthe

publication of this paper.

DATA AVAILABILITY

The data that support the findings of this study are available

from the corresponding author upon reasonable request.

REFERENCES

1A. V. G. Cavalieri, P. Jordan, Y. Gervais, M. Wei, and J. B. Freund, “Intermittent

sound generation and its control in a free-shear flow,” Phys. Fluids 22, 115113

(2010).

2D. Casalino, F. Diozzi, R. Sannino, and A. Paonessa, “Aircraft noise reduction

technologies: A bibliographic review,” Aerosp. Sci. Technol. 12, 1–17 (2008).

3R. Li, B. R. Noack, L. Cordier, J. Borée, and F. Harambat, “Drag reduction of a

car model by linear genetic programming control,” Exp. Fluids 58, 103 (2017).

4X.-W. Sun, W. Liu, and Z.-X. Chai, “Method of investigation for numerical

simulation on aero-optical effect based on WCNS-E-5,” AIAA J. 57(5), 1–13

(2019).

AIP Advances 11, 085322 (2021); doi: 10.1063/5.0062145 11, 085322-15

AIP Advances ARTICLE scitation.org/journal/adv

5T. C. Island, W. D. Urban, and M. G. Mungal, “Mixing enhancement in com-

pressible shear layers via sub-boundary layer disturbances,” Phys. Fluids 10,

1008–1020 (1998).

6H.Li,J.Tan,D.Zhang,andJ.Hou,“Investigationsonflowandgrowthproperties

of a cavity-actuated supersonic planar mixing layer downstream a thick splitter

plate,” Int. J. Heat Fluid Flow 74, 209–220 (2018).

7J. Tan, H. Li, and B. R. Noack, “On the cavity-actuated supersonic mixing layer

downstream a thick splitter plate,” Phys. Fluids 32, 096102 (2020).

8A. W. Vreman, N. D. Sandham, and K. H. Luo, “Compressible mixing

layer growth rate and turbulence characteristics,” J. Fluid Mech. 320, 235–258

(1996).

9B. Jones, “Statistical investigation of pressure and velocity fields in the turbulent

two-stream mixing layer,” in 4th Fluid and Plasma Dynamics Conference, 1971.

10G. L. Brown and A. Roshko, “On density effects and large structure in turbulent

mixing layers,” J. Fluid Mech. 64, 775–816 (1974).

11P. Ortwerthand A.Shine, “Onthe scaling ofplane turbulentshear layers,”Laser

Digest, Report No. AFWL-TR-77-118, 1977, p. 115.

12P. E. Dimotakis and G. L. Brown, “The mixing layer at high Reynolds number:

Large-structure dynamics and entrainment,” J. Fluid Mech. 78, 535–560 (1976).

13J. H. Bell and R. D. Mehta, “Development of a two-stream mixing layer from

tripped and untripped boundary layers,” AIAA J. 28, 2034–2042 (1990).

14S. G. Goebel, J. C. Dutton, H. Krier, and J. P. Renie, “Mean and turbulent

velocitymeasurements of supersonic mixing layers,” Miner. Deposita29, 263–272

(1994).

15S. G. Goebel and J. C. Dutton, “Experimental study of compressible turbulent

mixing layers,” AIAA J. 29, 538–546 (1991).

16C. Pantano and S. Sarkar, “A study of compressibility effects in the high-speed

turbulentshear layerusing directsimulation,” J. FluidMech. 451,329–371 (2002).

17W.D.UrbanandM.G.Mungal,“Planarvelocitymeasurementsincompressible

mixing layers,” J. Fluid Mech. 431, 189–222 (2001).

18M. G. Olsen and J. C. Dutton, “Planar velocity measurements in a weakly

compressible mixing layer,” J. Fluid Mech. 486, 51–77 (2003).

19D. W. Bogdanoff, “Compressibility effects in turbulent shear layers,” AIAA J.

21, 926–927 (1983).

20S. Ghosh and R. Friedrich, “Effects of radiative heat transfer on the turbulence

structure in inert and reacting mixing layers,” Phys. Fluids 27, 055107 (2015).

21M. Ferrer, P. José, G. Lehnasch, and A. Mura, “Compressibility and heat

release effects in high-speed reactive mixing layers I: Growth rates and turbulence

characteristics,” Combust. Flame 180, 284–303 (2017).

22B.Thurow,M.Samimy,andW.Lempert,“Compressibilityeffectsonturbulence

structures of axisymmetric mixing layers,” Phys. Fluids 15, 1755–1765 (2003).

23T. Kanda, K. Tani, and K. Kudo, “Conceptual study of a rocket-ramjet

combined-cycle engine for an aerospace plane,” J. Propul. Power 23, 301–309

(2007).

24Q. Dai, T. Jin, K. Luo, and J. Fan, “Direct numerical simulation of particle dis-

persion in a three-dimensional spatially developing compressible mixing layer,”

Phys. Fluids 30, 113301 (2018).

25Q. Dai, T. Jin, K. Luo, and J. Fan, “Direct numerical simulation of a three-

dimensional spatially evolving compressible mixing layer laden with particles. I.

Turbulent structures and asymmetric properties,” Phys. Fluids 31, 083302 (2019).

26D. Chong, Y. Bai, Q. Zhao, W. Chen, J. Yan, and Y. Hong, “Direct numerical

simulation of vortex structures during the late stage of the transition process in a

compressible mixing layer,” Phys. Fluids 33, 054108 (2021).

27T. Matsushima, K. Nagata, and T. Watanabe, “Wavelet analysis of shearless

turbulent mixing layer,” Phys. Fluids 33, 025109 (2021).

28E. J. Gutmark, K. C. Schadow, and K. H. Yu, “Mixing enhancement in super-

sonic free shear flows,” Annu. Rev. Fluid Mech. 27, 375–417 (1995).

29J. Tan, D. Zhang, and L. Lv, “A review on enhanced mixing methods in

supersonic mixing layer flows,” Acta Astronaut. 152, 310–324 (2018).

30X.-x. Fang, C.-b. Shen, M.-b. Sun, R. D. Sandberg, and P. Wang, “Flow struc-

turesof alobed mixer andeffects ofstreamwise vortices onmixing enhancement,”

Phys. Fluids 31, 066102 (2019).

31X.-x. Fang, C.-b. Shen, M.-b. Sun, H.-b. Wang, and P. Wang, “Turbulent struc-

tures and mixing enhancement with lobed mixers in a supersonic mixing layer,”

Phys. Fluids 32, 041701 (2020).

32S. Yadala, N. Benard, M. Kotsonis, and E. Moreau, “Effect of dielectric barrier

discharge plasma actuators on vortical structures in a mixing layer,” Phys. Fluids

32, 124111 (2020).

33R. D. Mehta, “Effect of velocity ratio on plane mixing layer development:

Influence of the splitter plate wake,” Exp. Fluids 10, 194–204 (1991).

34N. D. Sandham and R. D. Sandberg, “Direct numerical simulation of the early

developmentof a turbulentmixing layerdownstreamof asplitter plate,” J.Turbul.

10, N1 (2009).

35S. Laizet, S. Lardeau, and E. Lamballais, “Direct numerical simulation of

a mixing layer downstream a thick splitter plate,” Phys. Fluids 22, 015104

(2010).

36C. Braud, D. Heitz, P. Braud, G. Arroyo, and J. Delville, “Analysis of the wake

mixing-layer interaction using multiple plane PIV and 3D classical POD,” Exp.

Fluids 37, 95–104 (2004).

37M. Rajaee, S. K. F. Karlsson, and L. Sirovich, “Low-dimensional description

of free-shear-flow coherent structures and their dynamical behaviour,” J. Fluid

Mech. 258, 1–29 (1994).

38G. Song, F. Alizard, J.-C. Robinet, and X. Gloerfelt, “Global and Koopman

modes analysis of sound generation in mixing layers,” Phys. Fluids 25, 124101

(2013).

39R. K. Soni and A. De, “Investigation of mixing characteristics in strut injectors

using modal decomposition,” Phys. Fluids 30, 016108 (2018).

40P. J. Schmid, “Dynamic mode decomposition of numerical and experimental

data,” J. Fluid Mech. 656, 5–28 (2010).

41K. Taira, S. L. Brunton, S. T. M. Dawson, C. W. Rowley, T. Colonius, B. J.

McKeon, O. T. Schmidt, S. Gordeyev, V. Theofilis, and L. S. Ukeiley, “Modal

analysis of fluid flows: An overview,” AIAA J. 55, 4013–4041 (2017).

42B. R. Noack, W. Stankiewicz, M. Morzy´

nski, and P. J. Schmid, “Recursive

dynamicmodedecompositionoftransientandpost-transientwakeflows,” J. Fluid

Mech. 809, 843–872 (2016).

43S. Le Clainche and J. M. Vega, “Spatio-temporal Koopman decomposition,” J.

Nonlinear Sci. 28, 1793–1842 (2018).

44S. Le Clainche, J. M. Pérez, and J. M. Vega, “Spatio-temporal flow structures

in the three-dimensional wake of a circular cylinder,” Fluid Dyn. Res. 50, 051406

(2018).

45K. K. Chen, J. H. Tu, and C. W. Rowley, “Variants of dynamic mode decompo-

sition: Boundary condition, Koopman, and Fourier analyses,” J. Nonlinear Sci. 22,

887–915 (2012).

46J. Burkardt, M. Gunzburger, and H.-C. Lee, “Centroidal Voronoi tessellation-

based reduced-order modeling of complex systems,” SIAM J. Sci. Comput. 28,

459–484 (2006).

47J. Burkardt, M. Gunzburger, and H.-C. Lee, “POD and CVT-based reduced-

ordermodelingofNavier–Stokesflows,”Comput.MethodsAppl.Mech.Eng.196,

337–355 (2006).

48E. Kaiser, B. R. Noack, L. Cordier, A. Spohn, M. Segond, M. Abel, G. Daviller, J.

Östh, S. Krajnovi´

c, R. K. Niven et al., “Cluster-based reduced-order modelling of

a mixing layer,” J. Fluid Mech. 754, 365–414 (2014).

49J. Östh, E. Kaiser, S. Krajnovi´

c, and B. R. Noack, “Cluster-based reduced-order

modelling of the flow in the wake of a high speed train,” J. Wind Eng. Ind.

Aerodyn. 145, 327–338 (2015).

50Z. Wei, Z. Yang, C. Xia, and Q. Li, “Cluster-based reduced-order modeling of

the wake stabilization mechanism behind a twisted cylinder,” J. Wind Eng. Ind.

Aerodyn. 171, 288–303 (2017).

51H. Li and J. Tan, “Cluster-based Markov model to understand the transition

dynamics of a supersonic mixing layer,” Phys. Fluids 32, 056104 (2020).

52N. Ali, M. Calaf, and R. B. Cal, “Cluster-based probabilistic structure dynamical

model of wind turbine wake,” J. Turbul. 22, 497–516 (2021).

53A. A. Shahbazi and V. Esfahanian, “Reduced-order modeling of lead-acid bat-

tery using cluster analysis and orthogonal cluster analysis method,” Int. J. Energy

Res. 43, 6779–6798 (2019).

54H. Li, D. Fernex, R. Semaan, J. Tan, M. Morzy´

nski, and B. R. Noack, “Cluster-

based network model,” J. Fluid Mech. 906, A21 (2021).

55D. Fernex, B. R. Noack, and R. Semaan, “Cluster-based network

modeling—From snapshots to complex dynamical systems,” Sci. Adv. 7,

eabf5006 (2021).

AIP Advances 11, 085322 (2021); doi: 10.1063/5.0062145 11, 085322-16

AIP Advances ARTICLE scitation.org/journal/adv

56N. T.Clemens and M.G. Mungal,“Large-scalestructure andentrainment in the

supersonic mixing layer,” J. Fluid Mech. 284, 171–216 (1995).

57H.Li,J.-g.Tan,J.-w.Hou,andD.-d.Zhang,“Investigationsofself-excitedvibra-

tion in splitter plate with a cavity in the supersonic mixing layer,” Aerosp. Sci.

Technol. 74, 120–131 (2018).

58M. Germano, “Turbulence: The filtering approach,” J.Fluid Mech. 238,325–336

(1992).

59J. HinzeandJ.Hinze, Turbulence,McGraw-HillClassicTextbook ReissueSeries

(McGraw-Hill, 1975).

60J. Smagorinsky, “General circulation experiments with the primitive equations:

I. The basic experiment,” Mon. Weather Rev. 91, 99–164 (1963).

61T. J. Poinsot and S. K. Lelef, “Boundary conditions for direct sim-

ulations of compressible viscous flows,” J. Comput. Phys. 101, 104–129

(1992).

62K. Taira, A. G. Nair, and S. L. Brunton, “Network structure of two-dimensional

decaying isotropic turbulence,” J. Fluid Mech. 795, R2 (2016).

63A. Guha and K. Pradhan, “Secondary motion in three-dimensional branching

networks,” Phys. Fluids 29, 063602 (2017).

64K. Pradhan and A. Guha, “Fluid dynamics of oscillatory flow in three-

dimensional branching networks,” Phys. Fluids 31, 063601 (2019).

65J. L.Lumley, “Thestructure ofinhomogeneous turbulent flows,”in Atmospheric

Turbulence and Radio Wave Propagation (University of California Press, 1967)

pp. 166–177.

66J. MacQueen, “Some methods for classification and analysis of multivariate

observations,” in Proceedings of the Fifth Berkeley Symposium on Mathemat-

ical Statistics and Probability (University of California Press, 1967), Vol. 1,

pp. 281–297.

67A. Hadjadj and A. Kudryavtsev, “Computation and flow visualization in high-

speed aerodynamics,” J. Turbul. 6(16), 33–81 (2005).

68C. M. Stack, R. Bhaskaran, and M. C. Adler, and D. V. Gaitonde “Influence of

the splitter plate thickness on the near-wake dynamics of compressible turbulent

mixing layers,” in AIAA SciTech 2019 Forum, 2019.

AIP Advances 11, 085322 (2021); doi: 10.1063/5.0062145 11, 085322-17