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

The flow topologies of compressible large-scale distorted flames are studied bymeans 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 the Q-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 -QW 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.


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 (, , and ) were defined and some useful properties of the -- 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  is zero, and hence the majority of studies have concentrated on the two invariants  and .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 - 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 wherein  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 - 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][17][18][19] mainly focused on the variation of heat release rate, pressure, temperature, vorticity, shock waves, and so forth.Considering these, the properties of compressible shockflame 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.

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  [14].The eigenvalues ( 1 ,  2 and  3 ) of  satisfy the following characteristic equation: The quantities , , and  appearing in (1) are the first, second, and third invariants of , respectively, and the first invariant  reflects the compressibility of a fluid element; that is,  > 0 signifies fluid element compression,  < 0 signifies fluid element expansion, and  = 0 indicates that fluid element is incompressible [4,5].Furthermore,  may be broken up into a symmetric rate-of-strain tensor  and an antisymmetric rate-of-rotation tensor .The invariants of  can be written as   ,   , and   , and the invariants of  can be written as   ,   , and   .And the relations among the above invariants can be proved as follows: It should be noted that   = 0 by the invariant derivation of , and a new variable   is used in (2).In the above formulas, nonzero   represents the enstrophy density and   signifies the viscous dissipation of a fluid element.Hence, in view of  =   +   ,  > 0 indicates the fluid element rotation is dominant and  < 0 indicates the fluid element strain is dominant.Moreover,   represents the dissipation production rate, and   represents the enstrophy production rate.As a result,  reflects the competition between   and   .
The flow topology can be inferred using the threedimensional -- space, which can be partitioned into complicated spatial regions by a set of surfaces.The detailed surfaces can be expressed as  1 - 3 [4]: The above  1 - 3 surfaces split the -- space into different partitions, and each of these partitions corresponds to a particular topology.To simplify the analysis, the  value is fixed first, and we could study the flow topology in the visualized two-dimensional - diagram.Figure 1 shows three representative two-dimensional - diagrams at discrete different  values, wherein 1-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 - diagrams might be found in [14,15].

Numerical Method and Physical Model
where  is the density;   is the th velocity component ( = 1, 2, 3);  is the pressure;  is the total energy density;  = /( − 1) + 0.5 ∑ 3 =1  2  + , where  is the adiabatic index,  is the chemical energy release, and  is the mass fraction of the reactant.  is the viscous stress tensor: where   is the Kronecker function.The kinematic viscosity , diffusion , and heat conduction  have a similar exponential dependence on temperature [16]: where,  0 ,  0 , and  0 are constants,   is the specific heat at constant pressure, and  = 0.7.In this paper, ω is the chemical reaction rate, and a one-step reaction model is used: where  is the preexponential factor,   is the activation energy, and  is the universal gas constant.All the related parameters can be found in [18,19].

Numerical Scheme.
Equations ( 4) can be expressed as follows: The splitting algorithm is applied to solve these equations.Specifically, the inviscid part of (, , ) 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.

Numerical Validations.
The three-dimensional shockflame 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 C 2 H 4 /3O 2 /4N 2 with initial temperature  0 = 293 K at initial pressure  0 = 13.3 kPa.The initial spherical flame radius is 19 mm, and the incident shock wave with  = 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  = 0 m and  = 0 m and the right boundary of the  direction.The planes at  = 19 mm and  = 38 mm are planes of symmetry, and the left boundary condition in the  direction is the zerogradient 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 C 2 H 4 /3O 2 /4N 2 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 3 gives the comparisons between the experimental Schlieren images [20] and the computational Schlieren images.Figures 3(a To better understand the flame distortion process, Figure 4 presents the instantaneous three-dimensional flame results after the incident shock wave ( = 189 s) and after the reflected shock wave ( = 501 s), respectively; the isosurface with red color corresponds to the flame surface ( = 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].

Results and Discussion
Considering the above detailed topological methodology and numerical results, Figure 5 presents the distribution of the dimensionless  value at two typical times corresponding to Figure 4, wherein ⟨  ⟩ 1/2 is used for the nondimensionalization of the  value, and ⟨  ⟩ signifies the average value of the enstrophy   in the flame zone ( < 0.99).It is found that the distribution peak of the  value locates near zero, which indicates the distorted flame behind the shock wave mainly behaves as an incompressible flow.However, the  value has a larger scope and a more gentle change behind the reflected shock wave ( = 501 s) than that behind the incident shock wave ( = 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.When a  value is selected, the complicated threedimensional -- space could be simplified to a twodimensional - 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  and  conditioned upon selected dimensionless  values in the distorted flame zone.As discussed above, different  values signify the various levels of dilatations; namely,   2. It is worth noting that the probability of every selected  value must be moderate (to ensure enough fluid elements can be used for the probability statistics, as shown in Figure 5).4.1.Joint p.d.f. of - Invariants.Figures 6 and 7 present the joint p.d.f. of  and  invariants, wherein the dashed lines correspond to the bounds of different topology partitions in Figure 1.Clearly, the joint p.d.f. of the - 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 - 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  < 0 and  decreasing or the case of  > 0 and  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.

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 of - is  made in Figure 8.It can be found that foci (corresponding to 1 and 4 in Table 1) dominate the flame zone when  = 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., 4 in Table 1), which is in agreement with the results in [6].
In the flame zone with negative dilatation (i.e.,  > 0), the fluid element should stay in a state of stable compression; therefore, the probability of unstable 1 and stretching 4 decreases while the probability of stable and compressing 5 increases gradually.Further, the foci pattern corresponding to stable 4 and 5 can promote the rise of fluid element vorticity in the compressed state.
While, in the flame zone with positive dilatation (i.e.,  < 0), the fluid element should stay in a state of unstable stretching, the probability of stable 4 and compressing 1 Advances in Materials Science and Engineering decreases while probability of unstable and stretching 7 increases gradually.Further, the foci pattern corresponding to unstable 1 and 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 2 and stable 3 with node/saddle/saddle topology pattern both show monotonic variation with different  values.The probability of stable 3 increases gradually while the  value increases, and unstable 2 is obviously dominant in the expanding region ( < 0) which means the flame expansion makes more fluid elements unstable.Figure 8 also shows the probabilities of stable 6 and unstable 8 are tiny, and these indicate that the pure node topology patterns are basically not affected by the different dilatation levels.
According to Table 1, 1-8 can be divided into foci patterns (including 1, 4, 5, and 7) and node patterns (including 2, 3, 6, and 8), where foci represent rotation and node represents strain, or can be divided into stable patterns (including 3, 4, 5, and 6) and unstable patterns (including 1, 2, 7, and 8).Figure 9 shows the probability evolution of rotation () and strain () with the dimensionless  value at  = 501 s, wherein the probability of rotation is obtained from the probability summation of 1, 4, 5, and 7, and the probability of strain is obtained from the probability summation of 2, 3, 6, and 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 ( = 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 10 presents the probability evolution of stable and unstable flow topology patterns with the dimensionless  value at the same times as Figure 9, wherein the probability of stable topology is obtained from the probability summation of 3, 4, 5, and 6, and the probability of unstable topology is obtained from the probability summation of 1, 2, 7, and 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  values.That is, in the flame zone with negative dilatation ( > 0), a fluid element tends towards stability and the stable topology probability increases gradually, while, in the flame zone with positive dilatation ( < 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. .In the literature [5], it has been reported that a vortex tube is dominant when the probability distribution is located near the   axis; irrotational dissipation is dominant when the probability distribution is located near the − *  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  *  -  is inclined towards the   axis in the flame zone with zero dilatation (as shown in Figure 11(d)), which indicates that rotation is dominant when  = 0, while, in the region with dilatation level increases (positive and negative), the joint p.d.f. is inclined towards the − *  axis, which implies that irrotational dissipation increases in the compression and expansion regions.
After 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. of  *  -  spreads along the dashed line   = − *  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 − *  axis, which also implies that the irrotational dissipation is enhanced in the compression and expansion regions.
On the basis of these observations, Figures 11 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 largescale compressible flame region is higher than that when  = 0, and this conclusion agrees with the work by Wang and Lu [15].

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 invariant  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 the - 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 invariant  *  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.
) and3(b)  show the scenarios of the interaction between the incident shock wave and the spherical flame, while Figures3(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 4 :DistributionFigure 5 :
Figure 4: Different flame structures at selected typical times (isosurface with red color corresponds to flame surface; yellow color represents pressure).

Figure 6 :
Figure 6: Joint p.d.f. of - invariants in the flame zone at  = 189 s.

Figure 9 :Figure 10 :
Figure 9: Probability evolution of rotation and strain in the flame zone ( = 501 s).

4. 3 .
Joint p.d.f. of  *  -  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 invariant  *  of the rate-of-strain tensor S * (S = S * + I  /3, where I is the identity matrix, and   is the principal diagonal element of tensor S).Figures 11 and 12 show the joint p.d.f. of  *  -  in the flame zone corresponding to Figures 6 and 7, respectively, where  *  is negative and both  *  and   are nondimensionalized by ⟨  ⟩.The dashed lines in these figures mean   = − *