# Flow Topology of Three-Dimensional Spherical Flame in Shock Accelerated Flows

Authors: Yuejin Zhu; Lei Yu; Gang Dong; Jianfeng Pan; Zhenhua Pan
Journal: Advances in Materials Science and Engineering
 (2016)

---

## Abstract

The flow topologies of compressible large-scale distorted flames are studied by means of the analysis of the invariants of the velocity gradient tensor (VGT). The results indicate that compressibility plays a minor role in the distorted flame zone. And the joint probability density function (p.d.f.) of theQ
-
R diagram appears as a teardrop shape, which is a universal feature of turbulence. Therefore, the distorted flame exhibits the characteristic of large-scale turbulence combustion, especially behind the reflected shock wave, while the p.d.f. of the Q
S
⁎
-
Q
W diagram implies that the dissipation is enhanced in the compression and expansion regions, where it is higher than that when P
=
0. Furthermore, we identify that the flame evolution is dominated by rotation by means of a quantitative statistical study, and the SFS topology is the predominant flow pattern. Not surprisingly, negative dilatation could suppress the unstable topologies, whereas positive dilatation could suppress the stable topologies.


---

## Body

## 1. Introduction

To date, great effort has been made to study the flow topology of turbulent flows. The topological study is based on an analysis of the invariants of the velocity gradient tensor (VGT), which is inherently interesting and significant in the flow dynamics and mechanisms. Knowledge of the invariants is particularly important to visualize the different flow patterns [1] and to deduce the vortical structures [2]. The tensor invariants also can be employed to reveal the important turbulence processes such as dissipation, vortex stretching, and scalar mixing [3].A general classification approach is to analyse the flow topology, as set out by the work of Chong et al. [4], where the three invariants (P, Q, and R) were defined and some useful properties of the P
-
Q
-
R space were illustrated. Subsequently, Perry and Chong [5] clearly pointed out that this topology classification method could be applied to every point in the flow field, and they used this topology classification to describe the large-scale and fine scale eddying motions. For incompressible flows, the first invariant P is zero, and hence the majority of studies have concentrated on the two invariants Q and R. Chacín et al. [6] studied the velocity field of incompressible turbulent boundary layers and found that the joint probability density functions (p.d.f.s) of the Q
-
R diagram showed a trend towards a teardrop shape. Furthermore, the teardrop shape appears to be common to a large variety of flows, such as boundary layers [7], grid turbulence [8], and free shear flow [9].The previous works mainly focused on the flow topology in incompressible flows whereinP is zero, and many meaningful results were obtained. Compared to the above research, some flow topology studies on compressible flow were performed as well. Chen et al. [10] examined the flow topology in turbulent compressible mixing layer flows. Maekawa et al. [11] applied the topology method in compressible isotropic turbulent flows and found that fluid motions corresponding to high dissipation rates were characterized by unstable node/saddle/saddle type. Pirozzoli and Grasso [12] studied the influence of initial compressibility on flow topology in compressible isotropic turbulence and found that the p.d.f.s of the Q
-
R diagram exhibited a universal teardrop shape as in the incompressible case. Lee et al. [13] also investigated the influence of compressibility on velocity gradients in compressible turbulence, and they observed that the dilatation level increased progressively with the turbulent Mach number. Moreover, Suman and Girimaji [14] found that the velocity gradient structure and flow topology were similar to incompressible turbulence when conditioned upon zero dilatation. In addition, Wang and Lu [15] identified the preference topologies in compressible turbulent boundary layers and found that the locally compressed region was more stable than the locally expanding region. To the best of our knowledge, however, a relevant study of flow topology has never been performed for a typical compressible shock-flame interaction problem, where the relevant previous shock-flame interaction studies [16–19] mainly focused on the variation of heat release rate, pressure, temperature, vorticity, shock waves, and so forth.Considering these, the properties of compressible shock-flame interaction are investigated by means of a detailed numerical simulation. The characteristics of the invariants and flow topologies are analysed using the VGT in this paper. The purpose is to reveal the underlying mechanism of a fluid element rotation and strain in the complicated large-scale flame deformation processes, which can help us to achieve an improved understanding of the complicated phenomena in compressible reactive flows.


## 2. Invariants and Flow Topology

Perry and Chong [5] pointed out that the three-invariant topological classification can be applied to an arbitrary point in the flow field. Further, the topology at an arbitrary point can be deduced with knowledge of the eigenvalues of the three-dimensional VGT A [14]. The eigenvalues (λ
1, λ
2 and λ
3) of A satisfy the following characteristic equation:(1)
λ
3
+
P
λ
2
+
Q
λ
+
R
=
0
.The quantitiesP, Q, and R appearing in (1) are the first, second, and third invariants of A, respectively, and the first invariant P reflects the compressibility of a fluid element; that is, P
>
0 signifies fluid element compression, P
<
0 signifies fluid element expansion, and P
=
0 indicates that fluid element is incompressible [4, 5]. Furthermore, A may be broken up into a symmetric rate-of-strain tensor S and an antisymmetric rate-of-rotation tensor W. The invariants of S can be written as P
S, Q
S, and R
S, and the invariants of W can be written as P
W, Q
W, and R
W. And the relations among the above invariants can be proved as follows:(2)
P
=
P
S
+
P
W
;
Q
=
Q
S
+
Q
W
;
R
-
R
S
=
-
R
P
.It should be noted thatR
W
=
0 by the invariant derivation of W, and a new variable R
P is used in (2). In the above formulas, nonzero Q
W represents the enstrophy density and Q
S signifies the viscous dissipation of a fluid element. Hence, in view of Q
=
Q
S
+
Q
W, Q
>
0 indicates the fluid element rotation is dominant and Q
<
0 indicates the fluid element strain is dominant. Moreover, R
S represents the dissipation production rate, and R
P represents the enstrophy production rate. As a result, R reflects the competition between R
S and R
P.The flow topology can be inferred using the three-dimensionalP
-
Q
-
R space, which can be partitioned into complicated spatial regions by a set of surfaces. The detailed surfaces can be expressed as r
1
–
r
3 [4]:(3)
r
1
:
1
3
P
Q
-
2
9
P
2
-
2
27
-
3
Q
+
P
2
3
/
2
-
R
=
0
,
r
2
:
1
3
P
Q
-
2
9
P
2
+
2
27
-
3
Q
+
P
2
3
/
2
-
R
=
0
,
r
3
:
P
Q
-
R
=
0
.The abover
1
–
r
3 surfaces split the P
-
Q
-
R space into different partitions, and each of these partitions corresponds to a particular topology. To simplify the analysis, the P value is fixed first, and we could study the flow topology in the visualized two-dimensional Q
-
R diagram. Figure 1 shows three representative two-dimensional Q
-
R diagrams at discrete different P values, wherein A
1
–
A
8 represent different topology patterns, respectively, and the detailed descriptions and corresponding acronyms are listed in Table 1. Moreover, the detailed explanations of different flow topologies on the Q
-
R diagrams might be found in [14, 15].Table 1
Acronyms and descriptions of various flow partitions and topologies [15].
Partition
Acronym
Descriptions
A
1
UFC
Unstable focus compressing
A
2
UN/S/S
Unstable node/saddle/saddle
A
3
SN/S/S
Stable node/saddle/saddle
A
4
SFS
Stable focus stretching
A
5
SFC
Stable focus compressing
A
6
SN/SN/SN
Stable node/node/node
A
7
UFS
Unstable focus stretching
A
8
UN/UN/UN
Unstable node/node/nodeFigure 1
Partitions of the two-dimensionalQ
-
R diagrams at different P values [14, 15].
(a)
P
=
0
(b)
P
>
0
(c)
P
<
0


## 3. Numerical Method and Physical Model

### 3.1. Governing Equation

The three-dimensional, compressible, and reactive Navier-Stokes equations are expressed as follows:(4)
∂
ρ
∂
t
+
∂
∂
x
j
ρ
u
j
=
0
,
∂
ρ
u
i
∂
t
+
∂
∂
x
j
ρ
u
i
u
j
=
-
∂
p
∂
x
i
+
∂
τ
i
j
∂
x
j
,
∂
E
∂
t
+
∂
∂
x
j
u
j
E
+
p
u
j
=
-
∂
q
j
∂
x
j
k
∂
T
∂
x
j
+
∂
∂
x
j
u
i
τ
i
j
,
∂
ρ
Y
∂
t
+
∂
∂
x
j
ρ
Y
u
j
=
∂
∂
x
j
D
ρ
∂
Y
∂
x
j
-
ρ
ω
˙
,where ρ is the density; u
i is the ith velocity component (i
=
1,2
,
3); p is the pressure; E is the total energy density; E
=
p
/
(
γ
-
1
)
+
0.5
ρ
∑
i
=
1
3
u
i
2
+
ρ
q
Y, where γ is the adiabatic index, q is the chemical energy release, and Y is the mass fraction of the reactant. τ
i
j is the viscous stress tensor: (5)
τ
i
j
=
ρ
υ
∂
u
i
∂
x
j
+
∂
u
j
∂
x
i
-
2
3
δ
i
j
∂
u
k
∂
x
k
,where δ
i
j is the Kronecker function. The kinematic viscosity υ, diffusion D, and heat conduction k have a similar exponential dependence on temperature [16]:(6)
υ
=
υ
0
T
n
ρ
,
D
=
D
0
T
n
ρ
,
k
ρ
c
p
=
k
0
T
n
ρ
,where, υ
0, D
0, and k
0 are constants, c
p is the specific heat at constant pressure, and n
=
0.7. In this paper, ω
˙ is the chemical reaction rate, and a one-step reaction model is used:(7)
ω
˙
=
A
ρ
Y
exp
⁡
-
E
a
R
T
,where A is the preexponential factor, E
a is the activation energy, and R is the universal gas constant. All the related parameters can be found in [18, 19].


### 3.2. Numerical Scheme

Equations (4) can be expressed as follows:(8)
∂
U
¯
∂
t
+
∂
F
¯
∂
x
+
∂
G
¯
∂
y
+
∂
H
¯
∂
z
=
0
.The splitting algorithm is applied to solve these equations. Specifically, the inviscid part of (F
¯
,
G
¯
,
H
¯) is computed by a high accuracy five-point TVD scheme [21], the viscous part is solved by the second-order central difference scheme, and a second-order Runge-Kutta method is adopted for unsteady time marching.


### 3.3. Numerical Validations

The three-dimensional shock-flame computational setup is determined by the experiment of Thomas et al. [20]. Figure 2 gives the schematic of this computation, where the computational domain (gray box) is 170 mm in length, 19 mm in width, and 38 mm in height; such a quadrant domain of a shock tube has been applied in the studies of inert shock-bubble interaction [22] and reactive shock-flame interaction [23]. The premixed combustible gas is C2H4/3O2/4N2 with initial temperature T
0
=
293 K at initial pressure p
0
=
13.3 kPa. The initial spherical flame radius is 19 mm, and the incident shock wave with M
a
=
1.7 initially moves towards the right and interacts with the spherical flame; it reflects from the end wall and interacts with the distorted flame again. The adiabatic and no-slip wall condition is enforced on the planes of y
=
0 m and z
=
0 m and the right boundary of the x direction. The planes at z
=
19 mm and y
=
38 mm are planes of symmetry, and the left boundary condition in the x direction is the zero-gradient condition. Uniform hexahedral meshes are adopted and the mesh size is 0.38 mm, which is 2 times smaller than the reaction zone of laminar C2H4/3O2/4N2 flame under the same initial conditions [17]; thus it is enough for describing the large-scale flame evolution and distortion correctly. The calculation time step is 0.0684 µs; hence the Courant number is smaller than 0.2, which satisfies the stability criteria of time marching.Figure 2
Schematic of three-dimensional computational setup.Figure3 gives the comparisons between the experimental Schlieren images [20] and the computational Schlieren images. Figures 3(a) and 3(b) show the scenarios of the interaction between the incident shock wave and the spherical flame, while Figures 3(c) and 3(d) correspond to the interaction between the reflected shock wave and the distorted flame. Comparisons between experimental and computational Schlieren images imply the reliability of the numerical method and chemical reaction model.Figure 3
Comparisons between experimental (top) [20] and computational (bottom) Schlieren images: (a) t
=
99 µs; (b) t
=
189 µs; (c) t
=
382 µs; (d) t
=
501 µs.
(a)
(b)
(c)
(d)To better understand the flame distortion process, Figure4 presents the instantaneous three-dimensional flame results after the incident shock wave (t
=
189 μs) and after the reflected shock wave (t
=
501 μs), respectively; the isosurface with red color corresponds to the flame surface (Y
=
0.99) and yellow color represents pressure (15 kPa for Figure 4(a), 80 kPa for Figure 4(b)). For all the details of this computation, please refer to our previous work in [18, 19].Figure 4
Different flame structures at selected typical times (isosurface with red color corresponds to flame surface; yellow color represents pressure).
(a)
(b)


## 3.1. Governing Equation

The three-dimensional, compressible, and reactive Navier-Stokes equations are expressed as follows:(4)
∂
ρ
∂
t
+
∂
∂
x
j
ρ
u
j
=
0
,
∂
ρ
u
i
∂
t
+
∂
∂
x
j
ρ
u
i
u
j
=
-
∂
p
∂
x
i
+
∂
τ
i
j
∂
x
j
,
∂
E
∂
t
+
∂
∂
x
j
u
j
E
+
p
u
j
=
-
∂
q
j
∂
x
j
k
∂
T
∂
x
j
+
∂
∂
x
j
u
i
τ
i
j
,
∂
ρ
Y
∂
t
+
∂
∂
x
j
ρ
Y
u
j
=
∂
∂
x
j
D
ρ
∂
Y
∂
x
j
-
ρ
ω
˙
,where ρ is the density; u
i is the ith velocity component (i
=
1,2
,
3); p is the pressure; E is the total energy density; E
=
p
/
(
γ
-
1
)
+
0.5
ρ
∑
i
=
1
3
u
i
2
+
ρ
q
Y, where γ is the adiabatic index, q is the chemical energy release, and Y is the mass fraction of the reactant. τ
i
j is the viscous stress tensor: (5)
τ
i
j
=
ρ
υ
∂
u
i
∂
x
j
+
∂
u
j
∂
x
i
-
2
3
δ
i
j
∂
u
k
∂
x
k
,where δ
i
j is the Kronecker function. The kinematic viscosity υ, diffusion D, and heat conduction k have a similar exponential dependence on temperature [16]:(6)
υ
=
υ
0
T
n
ρ
,
D
=
D
0
T
n
ρ
,
k
ρ
c
p
=
k
0
T
n
ρ
,where, υ
0, D
0, and k
0 are constants, c
p is the specific heat at constant pressure, and n
=
0.7. In this paper, ω
˙ is the chemical reaction rate, and a one-step reaction model is used:(7)
ω
˙
=
A
ρ
Y
exp
⁡
-
E
a
R
T
,where A is the preexponential factor, E
a is the activation energy, and R is the universal gas constant. All the related parameters can be found in [18, 19].


## 3.2. Numerical Scheme

Equations (4) can be expressed as follows:(8)
∂
U
¯
∂
t
+
∂
F
¯
∂
x
+
∂
G
¯
∂
y
+
∂
H
¯
∂
z
=
0
.The splitting algorithm is applied to solve these equations. Specifically, the inviscid part of (F
¯
,
G
¯
,
H
¯) is computed by a high accuracy five-point TVD scheme [21], the viscous part is solved by the second-order central difference scheme, and a second-order Runge-Kutta method is adopted for unsteady time marching.


## 3.3. Numerical Validations

The three-dimensional shock-flame computational setup is determined by the experiment of Thomas et al. [20]. Figure 2 gives the schematic of this computation, where the computational domain (gray box) is 170 mm in length, 19 mm in width, and 38 mm in height; such a quadrant domain of a shock tube has been applied in the studies of inert shock-bubble interaction [22] and reactive shock-flame interaction [23]. The premixed combustible gas is C2H4/3O2/4N2 with initial temperature T
0
=
293 K at initial pressure p
0
=
13.3 kPa. The initial spherical flame radius is 19 mm, and the incident shock wave with M
a
=
1.7 initially moves towards the right and interacts with the spherical flame; it reflects from the end wall and interacts with the distorted flame again. The adiabatic and no-slip wall condition is enforced on the planes of y
=
0 m and z
=
0 m and the right boundary of the x direction. The planes at z
=
19 mm and y
=
38 mm are planes of symmetry, and the left boundary condition in the x direction is the zero-gradient condition. Uniform hexahedral meshes are adopted and the mesh size is 0.38 mm, which is 2 times smaller than the reaction zone of laminar C2H4/3O2/4N2 flame under the same initial conditions [17]; thus it is enough for describing the large-scale flame evolution and distortion correctly. The calculation time step is 0.0684 µs; hence the Courant number is smaller than 0.2, which satisfies the stability criteria of time marching.Figure 2
Schematic of three-dimensional computational setup.Figure3 gives the comparisons between the experimental Schlieren images [20] and the computational Schlieren images. Figures 3(a) and 3(b) show the scenarios of the interaction between the incident shock wave and the spherical flame, while Figures 3(c) and 3(d) correspond to the interaction between the reflected shock wave and the distorted flame. Comparisons between experimental and computational Schlieren images imply the reliability of the numerical method and chemical reaction model.Figure 3
Comparisons between experimental (top) [20] and computational (bottom) Schlieren images: (a) t
=
99 µs; (b) t
=
189 µs; (c) t
=
382 µs; (d) t
=
501 µs.
(a)
(b)
(c)
(d)To better understand the flame distortion process, Figure4 presents the instantaneous three-dimensional flame results after the incident shock wave (t
=
189 μs) and after the reflected shock wave (t
=
501 μs), respectively; the isosurface with red color corresponds to the flame surface (Y
=
0.99) and yellow color represents pressure (15 kPa for Figure 4(a), 80 kPa for Figure 4(b)). For all the details of this computation, please refer to our previous work in [18, 19].Figure 4
Different flame structures at selected typical times (isosurface with red color corresponds to flame surface; yellow color represents pressure).
(a)
(b)


## 4. Results and Discussion

Considering the above detailed topological methodology and numerical results, Figure5 presents the distribution of the dimensionless P value at two typical times corresponding to Figure 4, wherein Q
W
1
/
2 is used for the nondimensionalization of the P value, and Q
W signifies the average value of the enstrophy Q
W in the flame zone (Y
<
0.99). It is found that the distribution peak of the P value locates near zero, which indicates the distorted flame behind the shock wave mainly behaves as an incompressible flow. However, the P value has a larger scope and a more gentle change behind the reflected shock wave (t
=
501 µs) than that behind the incident shock wave (t
=
189 µs). This tendency indicates the strengthening of compressibility and dilatation of the distorted flame zone behind the reflected shock wave, and the compressibility characteristic becomes more obvious.Figure 5
Probability distribution of dimensionlessP value.When aP value is selected, the complicated three-dimensional P
-
Q
-
R space could be simplified to a two-dimensional Q
-
R diagram, as shown in Figure 1. In this way, complexity of the flow topology could be reduced, and we could investigate the influence of various levels of dilatations on the flow topology. Hence, in this paper, we analyse the joint probability density function (p.d.f.) of the dimensionless Q and R conditioned upon selected dimensionless P values in the distorted flame zone. As discussed above, different P values signify the various levels of dilatations; namely, P
<
0 (positive dilatation) means fluid element expansion and P
>
0 (negative dilatation) means fluid element compression. The specific selected dimensionless P values can be seen in Table 2. It is worth noting that the probability of every selected P value must be moderate (to ensure enough fluid elements can be used for the probability statistics, as shown in Figure 5).Table 2
Selected dimensionlessP values.
Times
P
/
Q
W
1
/
2
t = 189 µs
—
-
0.4
-
0.2
-
0.1
0
0.1
0.2
—
t = 501 µs
-
0.8
-
0.4
-
0.2
-
0.1
0
0.1
0.2
0.4


### 4.1. Joint p.d.f. ofQ-R Invariants

Figures6 and 7 present the joint p.d.f. of Q and R invariants, wherein the dashed lines correspond to the bounds of different topology partitions in Figure 1. Clearly, the joint p.d.f. of the Q
-
R diagram in the flame zone appears as a teardrop shape (as shown in Figures 6(d) and 7(e)). And the statistical characteristics about incompressible turbulence [6, 24] and compressible turbulence [14, 15] indicate that the teardrop shape of Q
-
R is a universal feature of turbulence. Therefore, Figures 6 and 7 imply that the distorted flame exhibits the characteristic of large-scale turbulence combustion, especially for the flame behind the reflected shock wave, which appears to have a more standard teardrop shape, and this indicates that the small scales of the turbulence have become more developed in the flame zone behind the reflected shock wave. For the case of P
<
0 and P decreasing or the case of P
>
0 and P increasing, the corresponding numbers of fluid elements both decrease, as shown in Figure 5; hence, the joint p.d.f. gets smaller and smaller gradually. However, the basic probability distribution shape always keeps the teardrop shape, from which it can be concluded that the compressed and expanding regions both demonstrate turbulence combustion characteristics.Figure 6
Joint p.d.f. ofQ
-
R invariants in the flame zone at t
=
189 µs.
(a)
P
/
Q
W
1
/
2
=
-
0.4
(b)
P
/
Q
W
1
/
2
=
-
0.2
(c)
P
/
Q
W
1
/
2
=
-
0.1
(d)
P
/
Q
W
1
/
2
=
0
(e)
P
/
Q
W
1
/
2
=
0.1
(f)
P
/
Q
W
1
/
2
=
0.2Figure 7
Joint p.d.f. ofQ
-
R invariants in the flame zone at t
=
501 µs.
(a)
P
/
Q
W
1
/
2
=
-
0.8
(b)
P
/
Q
W
1
/
2
=
-
0.4
(c)
P
/
Q
W
1
/
2
=
-
0.2
(d)
P
/
Q
W
1
/
2
=
-
0.1
(e)
P
/
Q
W
1
/
2
=
0
(f)
P
/
Q
W
1
/
2
=
0.1
(g)
P
/
Q
W
1
/
2
=
0.2
(h)
P
/
Q
W
1
/
2
=
0.4


### 4.2. Quantitative Analysis of Various Topologies

Due to the small scales of the turbulence that have become more developed behind the reflected shock wave, a quantitative study of the occurrence probabilities of various topologies ofQ
-
R is made in Figure 8. It can be found that foci (corresponding to A
1 and A
4 in Table 1) dominate the flame zone when P
=
0, which means rotation dominates the flow field with zero dilatation. Considering that the incompressible characteristics play an important role in the whole flame zone (as seen in Figure 5), rotation is the representative flow pattern in the distorted flame zone. Moreover, the invariants of the VGT show that a preference exists for the stable focus stretching (SFS) topology (i.e., A
4 in Table 1), which is in agreement with the results in [6].Figure 8
Occurrence probabilities of various topologies in the flame zone (t
=
501 µs).In the flame zone with negative dilatation (i.e.,P
>
0), the fluid element should stay in a state of stable compression; therefore, the probability of unstable A
1 and stretching A
4 decreases while the probability of stable and compressing A
5 increases gradually. Further, the foci pattern corresponding to stable A
4 and A
5 can promote the rise of fluid element vorticity in the compressed state.While, in the flame zone with positive dilatation (i.e.,P
<
0), the fluid element should stay in a state of unstable stretching, the probability of stable A
4 and compressing A
1 decreases while probability of unstable and stretching A
7 increases gradually. Further, the foci pattern corresponding to unstable A
1 and A
7 can promote the decline of fluid element vorticity in the expanding state, and this trend is in agreement with the fact that vorticity decreases when the flame zone expands [18]. In addition, the probabilities of unstable A
2 and stable A
3 with node/saddle/saddle topology pattern both show monotonic variation with different P values. The probability of stable A
3 increases gradually while the P value increases, and unstable A
2 is obviously dominant in the expanding region (P
<
0) which means the flame expansion makes more fluid elements unstable. Figure 8 also shows the probabilities of stable A
6 and unstable A
8 are tiny, and these indicate that the pure node topology patterns are basically not affected by the different dilatation levels.According to Table1, A
1
–
A
8 can be divided into foci patterns (including A
1, A
4, A
5, and A
7) and node patterns (including A
2, A
3, A
6, and A
8), where foci represent rotation and node represents strain, or can be divided into stable patterns (including A
3, A
4, A
5, and A
6) and unstable patterns (including A
1, A
2, A
7, and A
8). Figure 9 shows the probability evolution of rotation (F) and strain (N) with the dimensionless P value at t
=
501 µs, wherein the probability of rotation is obtained from the probability summation of A
1, A
4, A
5, and A
7, and the probability of strain is obtained from the probability summation of A
2, A
3, A
6, and A
8 in Figure 8. From this figure, we observe that the rotation probability decreases and the strain probability increases with dilatation level increases (positive and negative), but, as a whole, the percentage of rotation is slightly larger than that of strain, particularly in the flame zone dominated by incompressible characteristic (P
=
0), and this implies that the rotation is dominant in the flame zone again, which is in accordance with the results of Figure 8.Figure 9
Probability evolution of rotation and strain in the flame zone (t
=
501 µs).Figure10 presents the probability evolution of stable and unstable flow topology patterns with the dimensionless P value at the same times as Figure 9, wherein the probability of stable topology is obtained from the probability summation of A
3, A
4, A
5, and A
6, and the probability of unstable topology is obtained from the probability summation of A
1, A
2, A
7, and A
8 in Figure 8. From this figure, we could find that the probabilities of stable and unstable topology patterns both exhibit a monotonic variation with the different P values. That is, in the flame zone with negative dilatation (P
>
0), a fluid element tends towards stability and the stable topology probability increases gradually, while, in the flame zone with positive dilatation (P
<
0), a fluid element tends towards instability and the unstable topology dominates, in accordance with the conclusions of Suman and Girimaji [14] on compressible isotropic turbulence and Wang and Lu [15] on compressible turbulent boundary layers.Figure 10
Probability evolution of stable and unstable topologies in the flame zone (t
=
501 µs).


### 4.3. Joint p.d.f. ofQ
S
∗-Q
W Invariants

In addition, the relationship between dissipation and enstrophy in the flame zone is also studied. For comparison with the literature, dissipation in the flame zone is represented by the second invariantQ
S
∗ of the rate-of-strain tensor S
∗ (S
=
S
∗
+
I
S
i
i
/
3, where I is the identity matrix, and S
i
i is the principal diagonal element of tensor S). Figures 11 and 12 show the joint p.d.f. of Q
S
∗
-
Q
W in the flame zone corresponding to Figures 6 and 7, respectively, where Q
S
∗ is negative and both Q
S
∗ and Q
W are nondimensionalized by Q
W. The dashed lines in these figures mean Q
W
=
-
Q
S
∗. In the literature [5], it has been reported that a vortex tube is dominant when the probability distribution is located near the Q
W axis; irrotational dissipation is dominant when the probability distribution is located near the -
Q
S
∗ axis; and a vortex sheet is dominant when the probability distribution is located near the dashed line. In Figure 11, we find that the joint p.d.f. of Q
S
∗
-
Q
W is inclined towards the Q
W axis in the flame zone with zero dilatation (as shown in Figure 11(d)), which indicates that rotation is dominant when P
=
0, while, in the region with dilatation level increases (positive and negative), the joint p.d.f. is inclined towards the -
Q
S
∗ axis, which implies that irrotational dissipation increases in the compression and expansion regions.Figure 11
Joint p.d.f. ofQ
S
∗
-
Q
W invariants in the flame zone at t
=
189 µs.
(a)
P
/
Q
W
1
/
2
=
-
0.4
(b)
P
/
Q
W
1
/
2
=
-
0.2
(c)
P
/
Q
W
1
/
2
=
-
0.1
(d)
P
/
Q
W
1
/
2
=
0
(e)
P
/
Q
W
1
/
2
=
0.1
(f)
P
/
Q
W
1
/
2
=
0.2Figure 12
Joint p.d.f. ofQ
S
∗
-
Q
W invariants in the flame zone at t
=
501 µs.
(a)
P
/
Q
W
1
/
2
=
-
0.8
(b)
P
/
Q
W
1
/
2
=
-
0.4
(c)
P
/
Q
W
1
/
2
=
-
0.2
(d)
P
/
Q
W
1
/
2
=
-
0.1
(e)
P
/
Q
W
1
/
2
=
0
(f)
P
/
Q
W
1
/
2
=
0.1
(g)
P
/
Q
W
1
/
2
=
0.2
(h)
P
/
Q
W
1
/
2
=
0.4After the reflected shock wave interacts with the distorted flame, the flame zone is in a high turbulent and reactive state. The joint p.d.f. ofQ
S
∗
-
Q
W spreads along the dashed line Q
W
=
-
Q
S
∗ in Figure 12, which indicates that there exist simultaneously high dissipation and enstrophy and that a vortex sheet is dominant in the flame zone with zero dilatation (as shown in Figure 12(e)). However, in the region with dilatation level increases (positive and negative), the joint p.d.f. is inclined towards the -
Q
S
∗ axis, which also implies that the irrotational dissipation is enhanced in the compression and expansion regions.On the basis of these observations, Figures11 and 12 reflect that the percentage of rotation (vorticity) decreases gradually in the region with the dilatation level increases (positive and negative), which is in accordance with Figure 9. In addition, it is worth noting that dissipation in the large-scale compressible flame region is higher than that when P
=
0, and this conclusion agrees with the work by Wang and Lu [15].


## 4.1. Joint p.d.f. ofQ-R Invariants

Figures6 and 7 present the joint p.d.f. of Q and R invariants, wherein the dashed lines correspond to the bounds of different topology partitions in Figure 1. Clearly, the joint p.d.f. of the Q
-
R diagram in the flame zone appears as a teardrop shape (as shown in Figures 6(d) and 7(e)). And the statistical characteristics about incompressible turbulence [6, 24] and compressible turbulence [14, 15] indicate that the teardrop shape of Q
-
R is a universal feature of turbulence. Therefore, Figures 6 and 7 imply that the distorted flame exhibits the characteristic of large-scale turbulence combustion, especially for the flame behind the reflected shock wave, which appears to have a more standard teardrop shape, and this indicates that the small scales of the turbulence have become more developed in the flame zone behind the reflected shock wave. For the case of P
<
0 and P decreasing or the case of P
>
0 and P increasing, the corresponding numbers of fluid elements both decrease, as shown in Figure 5; hence, the joint p.d.f. gets smaller and smaller gradually. However, the basic probability distribution shape always keeps the teardrop shape, from which it can be concluded that the compressed and expanding regions both demonstrate turbulence combustion characteristics.Figure 6
Joint p.d.f. ofQ
-
R invariants in the flame zone at t
=
189 µs.
(a)
P
/
Q
W
1
/
2
=
-
0.4
(b)
P
/
Q
W
1
/
2
=
-
0.2
(c)
P
/
Q
W
1
/
2
=
-
0.1
(d)
P
/
Q
W
1
/
2
=
0
(e)
P
/
Q
W
1
/
2
=
0.1
(f)
P
/
Q
W
1
/
2
=
0.2Figure 7
Joint p.d.f. ofQ
-
R invariants in the flame zone at t
=
501 µs.
(a)
P
/
Q
W
1
/
2
=
-
0.8
(b)
P
/
Q
W
1
/
2
=
-
0.4
(c)
P
/
Q
W
1
/
2
=
-
0.2
(d)
P
/
Q
W
1
/
2
=
-
0.1
(e)
P
/
Q
W
1
/
2
=
0
(f)
P
/
Q
W
1
/
2
=
0.1
(g)
P
/
Q
W
1
/
2
=
0.2
(h)
P
/
Q
W
1
/
2
=
0.4


## 4.2. Quantitative Analysis of Various Topologies

Due to the small scales of the turbulence that have become more developed behind the reflected shock wave, a quantitative study of the occurrence probabilities of various topologies ofQ
-
R is made in Figure 8. It can be found that foci (corresponding to A
1 and A
4 in Table 1) dominate the flame zone when P
=
0, which means rotation dominates the flow field with zero dilatation. Considering that the incompressible characteristics play an important role in the whole flame zone (as seen in Figure 5), rotation is the representative flow pattern in the distorted flame zone. Moreover, the invariants of the VGT show that a preference exists for the stable focus stretching (SFS) topology (i.e., A
4 in Table 1), which is in agreement with the results in [6].Figure 8
Occurrence probabilities of various topologies in the flame zone (t
=
501 µs).In the flame zone with negative dilatation (i.e.,P
>
0), the fluid element should stay in a state of stable compression; therefore, the probability of unstable A
1 and stretching A
4 decreases while the probability of stable and compressing A
5 increases gradually. Further, the foci pattern corresponding to stable A
4 and A
5 can promote the rise of fluid element vorticity in the compressed state.While, in the flame zone with positive dilatation (i.e.,P
<
0), the fluid element should stay in a state of unstable stretching, the probability of stable A
4 and compressing A
1 decreases while probability of unstable and stretching A
7 increases gradually. Further, the foci pattern corresponding to unstable A
1 and A
7 can promote the decline of fluid element vorticity in the expanding state, and this trend is in agreement with the fact that vorticity decreases when the flame zone expands [18]. In addition, the probabilities of unstable A
2 and stable A
3 with node/saddle/saddle topology pattern both show monotonic variation with different P values. The probability of stable A
3 increases gradually while the P value increases, and unstable A
2 is obviously dominant in the expanding region (P
<
0) which means the flame expansion makes more fluid elements unstable. Figure 8 also shows the probabilities of stable A
6 and unstable A
8 are tiny, and these indicate that the pure node topology patterns are basically not affected by the different dilatation levels.According to Table1, A
1
–
A
8 can be divided into foci patterns (including A
1, A
4, A
5, and A
7) and node patterns (including A
2, A
3, A
6, and A
8), where foci represent rotation and node represents strain, or can be divided into stable patterns (including A
3, A
4, A
5, and A
6) and unstable patterns (including A
1, A
2, A
7, and A
8). Figure 9 shows the probability evolution of rotation (F) and strain (N) with the dimensionless P value at t
=
501 µs, wherein the probability of rotation is obtained from the probability summation of A
1, A
4, A
5, and A
7, and the probability of strain is obtained from the probability summation of A
2, A
3, A
6, and A
8 in Figure 8. From this figure, we observe that the rotation probability decreases and the strain probability increases with dilatation level increases (positive and negative), but, as a whole, the percentage of rotation is slightly larger than that of strain, particularly in the flame zone dominated by incompressible characteristic (P
=
0), and this implies that the rotation is dominant in the flame zone again, which is in accordance with the results of Figure 8.Figure 9
Probability evolution of rotation and strain in the flame zone (t
=
501 µs).Figure10 presents the probability evolution of stable and unstable flow topology patterns with the dimensionless P value at the same times as Figure 9, wherein the probability of stable topology is obtained from the probability summation of A
3, A
4, A
5, and A
6, and the probability of unstable topology is obtained from the probability summation of A
1, A
2, A
7, and A
8 in Figure 8. From this figure, we could find that the probabilities of stable and unstable topology patterns both exhibit a monotonic variation with the different P values. That is, in the flame zone with negative dilatation (P
>
0), a fluid element tends towards stability and the stable topology probability increases gradually, while, in the flame zone with positive dilatation (P
<
0), a fluid element tends towards instability and the unstable topology dominates, in accordance with the conclusions of Suman and Girimaji [14] on compressible isotropic turbulence and Wang and Lu [15] on compressible turbulent boundary layers.Figure 10
Probability evolution of stable and unstable topologies in the flame zone (t
=
501 µs).


## 4.3. Joint p.d.f. ofQ
S
∗-Q
W Invariants

In addition, the relationship between dissipation and enstrophy in the flame zone is also studied. For comparison with the literature, dissipation in the flame zone is represented by the second invariantQ
S
∗ of the rate-of-strain tensor S
∗ (S
=
S
∗
+
I
S
i
i
/
3, where I is the identity matrix, and S
i
i is the principal diagonal element of tensor S). Figures 11 and 12 show the joint p.d.f. of Q
S
∗
-
Q
W in the flame zone corresponding to Figures 6 and 7, respectively, where Q
S
∗ is negative and both Q
S
∗ and Q
W are nondimensionalized by Q
W. The dashed lines in these figures mean Q
W
=
-
Q
S
∗. In the literature [5], it has been reported that a vortex tube is dominant when the probability distribution is located near the Q
W axis; irrotational dissipation is dominant when the probability distribution is located near the -
Q
S
∗ axis; and a vortex sheet is dominant when the probability distribution is located near the dashed line. In Figure 11, we find that the joint p.d.f. of Q
S
∗
-
Q
W is inclined towards the Q
W axis in the flame zone with zero dilatation (as shown in Figure 11(d)), which indicates that rotation is dominant when P
=
0, while, in the region with dilatation level increases (positive and negative), the joint p.d.f. is inclined towards the -
Q
S
∗ axis, which implies that irrotational dissipation increases in the compression and expansion regions.Figure 11
Joint p.d.f. ofQ
S
∗
-
Q
W invariants in the flame zone at t
=
189 µs.
(a)
P
/
Q
W
1
/
2
=
-
0.4
(b)
P
/
Q
W
1
/
2
=
-
0.2
(c)
P
/
Q
W
1
/
2
=
-
0.1
(d)
P
/
Q
W
1
/
2
=
0
(e)
P
/
Q
W
1
/
2
=
0.1
(f)
P
/
Q
W
1
/
2
=
0.2Figure 12
Joint p.d.f. ofQ
S
∗
-
Q
W invariants in the flame zone at t
=
501 µs.
(a)
P
/
Q
W
1
/
2
=
-
0.8
(b)
P
/
Q
W
1
/
2
=
-
0.4
(c)
P
/
Q
W
1
/
2
=
-
0.2
(d)
P
/
Q
W
1
/
2
=
-
0.1
(e)
P
/
Q
W
1
/
2
=
0
(f)
P
/
Q
W
1
/
2
=
0.1
(g)
P
/
Q
W
1
/
2
=
0.2
(h)
P
/
Q
W
1
/
2
=
0.4After the reflected shock wave interacts with the distorted flame, the flame zone is in a high turbulent and reactive state. The joint p.d.f. ofQ
S
∗
-
Q
W spreads along the dashed line Q
W
=
-
Q
S
∗ in Figure 12, which indicates that there exist simultaneously high dissipation and enstrophy and that a vortex sheet is dominant in the flame zone with zero dilatation (as shown in Figure 12(e)). However, in the region with dilatation level increases (positive and negative), the joint p.d.f. is inclined towards the -
Q
S
∗ axis, which also implies that the irrotational dissipation is enhanced in the compression and expansion regions.On the basis of these observations, Figures11 and 12 reflect that the percentage of rotation (vorticity) decreases gradually in the region with the dilatation level increases (positive and negative), which is in accordance with Figure 9. In addition, it is worth noting that dissipation in the large-scale compressible flame region is higher than that when P
=
0, and this conclusion agrees with the work by Wang and Lu [15].


## 5. Conclusions

A three-dimensional numerical simulation of flame evolution in shock accelerated flows is conducted. The flow topologies and statistical properties of compressible large-scale distorted flame are studied using the invariants of the VGT, which is helpful to understand the intrinsic flow dynamics. The main conclusions are summarized as follows:(1)
The probability distribution of the first invariantP in the flame zone shows that compressibility plays a minor role in the large-scale distorted flame in the present study.
(2)
The joint p.d.f. of theQ
-
R diagram appears as a teardrop shape, which means the large-scale distorted flame behind the reflected shock wave is in a high turbulent and reactive state. And different dilatation levels (positive and negative) can affect the topology characteristics.
(3)
The occurrence probabilities of various topologies are analysed. We find that the rotation is the representative flow pattern in the distorted flame zone, and the SFS topology is favored. Further, negative dilatation can promote an increase of the percentage of stable topology, while positive dilatation can promote an increase of the percentage of unstable topology.
(4)
Additionally, dissipation (represented by the second invariantQ
S
∗ of the rate-of-strain tensor) in the compressible large-scale distorted flame zone is enhanced with the dilatation level increases, which is higher than that in the incompressible zone.



---

*Source: 3158091-2016-11-17.xml*