- Email: [email protected]
Direct numerical simulation of fine flow structures of subsonic-supersonic mixing layer
JID:AESCTE AID:105431 /FLA
[m5G; v1.261; Prn:3/10/2019; 11:20] P.1 (1-18)
Aerospace Science and Technology ••• (••••) ••••••
1
Contents lists available at ScienceDirect
67 68
2 3
Aerospace Science and Technology
4
69 70 71
5
72
6
www.elsevier.com/locate/aescte
7
73
8
74
9
75
10
76
11 12 13
Direct numerical simulation of fine flow structures of subsonic-supersonic mixing layer
77 78 79 80
14 15 16 17 18 19 20
Liu Yang
a,∗
a
, Ma Dong , Fu Benshuai
b,c
a
, Li Qiang , Zhang Chenxi
d
a
Science and Technology on Combustion, Internal Flow and Thermo-Structure Laboratory, Northwestern Polytechnical University, Xi’an, China b No. 713 Institute, China Shipbuilding Industry Co., Zhengzhou 450015, PR China c Henan Key Laboratory of Underwater Intelligence Equipment, Zhengzhou 450015, PR China d School of Astronautics, Northwestern Polytechnical University, Xi’an, China
a r t i c l e
i n f o
25 26 27 28 29 30 31 32
83 84 85 86 88
a b s t r a c t
89
23 24
82
87
21 22
81
Article history: Received 29 July 2019 Received in revised form 4 September 2019 Accepted 25 September 2019 Available online xxxx Keywords: Subsonic-supersonic mixing layer Compressibility Direct numerical simulation Dynamic modal decomposition
33 34 35 36 37 38 39 40 41 42 43 44 45 46
Subsonic-supersonic mixing flow is one of the important research fields in turbulence research. For example, in the combustion chamber of the rocket ramjet combination engine, the rocket jet flow and the air inflow is a typical shear mixing flow, which has the characteristics of high convection Mach number and large flow parameter gradient. It is of great significance to explore the development law and flow structure of subsonic-supersonic mixing layer to enhance the mixing and thus improve the performance of rocket ramjet engine. In this paper, in order to further study the evolution process and obtain the evolution mechanism and fine flow structure of subsonic-supersonic mixing layer, three groups of high-order format direct numerical simulations of typical working conditions are carried out. The convective Mach numbers are 0.69, 0.92 and 1.27, respectively. The direct numerical simulation results show that: (1) the large-scale structure in the subsonic-supersonic shear mixing flow formed the small-scale structure earlier, and the flow was dominated by the small-scale and broken eddy structure; (2) the development process of the subsonic-supersonic shear mixing flow was mainly divided into three stages: laminar flow zone, transition zone and development zone. The transition zone was a transition zone from two-dimensional to three-dimensional, which began to show three-dimensional flow characteristics. The vorticity shows a twisted-winding structure. The large-scale Coherent structure, the Λ eddy and hairpin eddy, can be observed, and the Λ eddy structure is gradually elongated with the increase of convective Mach number; (3) the shocklets structure mainly occurs in the stage when the mixing layer has not fully developed. The starting point of the shocklets structure is near the vortex core. With the increase of convective Mach number, the position of the shocklets generation moves to the front of the vortex core; (4) the dynamic mode decomposition results show that the eigenvalues of each mode are basically located in the unit circle, which indicates that the calculated modes are stable; in addition, it is found that there is a certain dominant frequency in the flow structure, which can provide a reference for the active mixing enhancement method. © 2019 Elsevier Masson SAS. All rights reserved.
90 91 92 93 94 95 96 97 98 99 100 101 102 103 104 105 106 107 108 109 110 111 112
47
113
48
114
49
115 116
50 51
1. Introduction
52 53 54 55 56 57 58 59
In rocket ramjet engines, the rocket gas with high-temperature and high-speed and air flow with low-temperature and low-speed will form a Subsonic-supersonic shear flow with a large parameter gradient. Due to the high compressibility of the shear flow, the growth rate of the shear layer is small, the turbulence and the Reynolds stress are small, the stability of the shear layer is enhanced, and it is not easy to form a large-scale quasi-ordered
60 61 62 63 64 65 66
*
Corresponding author. E-mail address: [email protected] (L. Yang).
https://doi.org/10.1016/j.ast.2019.105431 1270-9638/© 2019 Elsevier Masson SAS. All rights reserved.
structure, leading to difficulties in the transfer of energy and momentum in the two flows, which directly affects the engine’s blending and combustion efficiency. The remarkable feature of the subsonic-supersonic shear mixed flow is its strong compressibility. The shear mixing flow was first studied on the low-speed incompressible flow blending problem and then gradually extended to high-speed compressible flow mixing. Bogdanoff et al. [1] proposed the concept of convective Mach number (Mc), which describes the effect of compressibility on the mixing layer by convection Mach number. It has been studied that the compressibility of the flow can reduce the growth rate of the shear layer and enhance the stability of flow [2,3].
117 118 119 120 121 122 123 124 125 126 127 128 129 130 131 132
JID:AESCTE
AID:105431 /FLA
2
1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48
[m5G; v1.261; Prn:3/10/2019; 11:20] P.2 (1-18)
L. Yang et al. / Aerospace Science and Technology ••• (••••) ••••••
In the free shear layer, due to Kelvin-Helmholtz (K-H) instability, the large-scale vortex structure is rolled up, paired, and fused as the flow field develops. These large-scale vortex structures still exist under supersonic flow conditions. First, Brown and Roshko et al. [4] demonstrated the existence of large-scale reverse-order structures of mixed layers in experiments at low-speed flows. Papamoschou and Roshko et al. [5] demonstrated the existence of large-scale structures in the compressible mixing layer. Later, many researchers began to capture the flow field structure by means of schlieren, shadow, plane laser Mie scattering (PLMS), particle image velocimetry (PIV) and other measurements [6–8]. Olsen et al. [9] used the schlieren technique to capture the shear flow field, and the outline of the shear layer was observed from the schlieren. Clemens and Petullo et al. [10] used optical measurements to observe the mixing layer. It was found that the vortex structure of the compressible mixing layer under the same convective Mach number is different from that of the incompressible mixing layer. The pulsation velocity obtained in the experiment revealed that the compressible shear mixed layer vortex structure has threedimensional characteristics. Experimental research on the fine structure of the shear layer has been difficult in recent decades. Numerical simulation method has become a very important means to study the fine flow of turbulent flow in shear layer. Numerical simulation method has become a very important means to study the precise flow of shear layer turbulence. Considerable recent progress has been made in the study of DES series of methods [11–13]. Doris et al. [14] studied the three-dimensional compressive shear mixing flow (Mc = 0.64) using the LES method. The results show that both the streamwise and the spanning pulsation velocity have selfsimilarity. Fu Song and Li Qibing et al. [15] researched the supersonic shear flow by the LES method and found that small shock wave structures appeared in the flow field, and the compressibility reduced the parameters such as Reynolds stress. Two-dimensional and three-dimensional simulations of a supersonic mixing layer are performed using an in-house hybrid LES code by Konark et al. [16]. They found that three-dimensional predictions of growth rate and statistical moments in the self similar regime are consistent with past experimental and direct numerical simulation studies. Recently, some researchers have studied the effects of shock and compressibility on the shear layer [17,18]. At present, researchers have conducted a lot of research on shear mixing flow and made great progress. However, further research is needed in the fine flow structure of the subsonicsupersonic shear mixed layer. Therefore, in this paper, the relevant parameters were changed to study the fine flow structure of the shear mixed layer, and the flow structure under different working conditions was compared.
49 50
2. Direct numerical simulation method
51 52 53 54 55 56 57 58
Through theoretical analysis, it is found that decreasing density and increasing temperature can reduce Reynolds number of the flow, thus greatly reducing calculated amount of direct numerical simulation. Therefore, this paper conducts direct numerical simulation research on the subsonic-supersonic shear mixing flow under the condition of low pressure and high temperature, and obtains the fine flow structure of the shear mixing flow.
59 60
2.1. Numerical computation method
61 62 63 64 65 66
2.1.1. Numerical computation platform In the direct numerical simulation study of this paper, the OpenCFD program developed by Li Xinliang is adopted [19]. The core of the program is the compressible solver of N-S equation. The numerical method of the program is the finite difference method,
67 68 69 70 71 72 73 74 75 76 77 78 79 80 81 82 83 84 85
Fig. 1. Turbulence energy spectrum analysis diagram.
86 87
which includes a variety of high-precision difference schemes currently popular, such as central difference scheme, upwind difference scheme and WENO difference scheme, etc., and also includes the difference scheme with its own characteristics, such as upwind compact scheme, group speed control difference scheme, etc. In the direct numerical simulation, the fourth order RungeKutta scheme is adopted for the time term, the sixth order central difference method is adopted for the viscous term, and the seventh order WENO method is adopted for the non-viscous term. The entire simulation process was completed on the Tianhe-2 supercomputer of the national supercomputing center in Guangzhou, using 15 nodes and 480 cores for calculation, and the calculation time of a single example was about 150 hours.
88 89 90 91 92 93 94 95 96 97 98 99 100 101
2.1.2. Numerical computation method validation More detailed experiments and direct numerical simulation studies have been carried out on the supersonic-supersonic shear mixing flow. In the study of this paper, the turbulence energy spectrum of direct numerical simulation results is analyzed, and the validity of grid scale is judged by the turbulence energy spectrum. Based on isotropic turbulence theory proposed by Kolmogorov, the slope of the turbulence energy spectrum in the inertial range is −5/3, and the absolute value in the dissipative zone is greater than −5/3. In the direct numerical simulation performed in this study, an intercepted mixed flow fully develops a section of the area. The position of the spectrum is selected in the streamwise range of 400–600. The resulting turbulent energy spectrum is shown in Fig. 1. The slope of the energy spectrum in the low-wavenumber region is −5/3, which satisfies Kolmogorov’s turbulence theory. In the high wavenumber region, a −5 slope indicates the minimum grid size of the direct numerical simulation observed in the dissipative region, thereby satisfying the grid requirements of the direct numerical simulation.
102 103 104 105 106 107 108 109 110 111 112 113 114 115 116 117 118 119 120 121 122
2.2. Numerical calculation parameters
123 124
2.2.1. Numerical calculation conditions All parameters are normalized by primary stream flow parameters. Given the density (d), velocity (U ), temperature (T ) and other parameters of the incoming flow, dimensionless flow parameters are shown in Table 1, and flow parameters reflecting the actual flow are shown in Table 2. In this chapter, three groups of direct numerical simulation calculations were carried out, and the convective Mach Numbers were 0.69, 0.92 and 1.19, respectively. The
125 126 127 128 129 130 131 132
JID:AESCTE AID:105431 /FLA
[m5G; v1.261; Prn:3/10/2019; 11:20] P.3 (1-18)
L. Yang et al. / Aerospace Science and Technology ••• (••••) ••••••
1 2
Table 1 Dimensionless flow parameters.
3
Primary stream
67 68
6 7 8
U T d Mc
69
Secondary stream
4 5
3
1 1 1 –
D1
D2
D3
0.199 1.768 0.565 0.69
0.133 0.786 1.273 0.92
0.111 1.222 0.818 1.19
70 71 72 73 74 75
9 10 11
76
Table 2 Actual flow parameters.
77 78
12
Primary stream
13 14 15 16 17 18 19
U , m/s T *, K T, K P , kPa Ma Mc
Secondary stream
79
D1
D2
D3
D1
80
845.2 800.0 444.4 10 2.0 0.69
1267.9 1800.0 1000.0 10 2.0 0.92
1524.9 1800.0 642.9 10 3.0 1.19
168.6 800.0 785.9 10 0.3 –
81 82 83 84 85
20
86
21
87
22
Fig. 3. Mach number distribution contours of D1–D3.
24 25 26 27 28 29
Fig. 2. Schematic diagram of computational area.
30 31 32
corresponding dimensionless Reynolds Numbers were 1346, 523 and 1309, respectively.
33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60 61 62 63 64 65 66
88 89
23
2.2.2. Mesh and boundary conditions (1) Mesh The finite-difference method is difficult to deal with complex geometry, so the classical cubic configuration in the plane shear mixing flow is adopted here, as shown in Fig. 2. The reference length is 5 × 10−4 m. The dimensionless length and the number of grid nodes in X , Y and Z directions are 600, 120, 60 and 3500, 500, 300, respectively, and the total amount of grid is 525 million. In the direction of X , the first 90% grid nodes are equally spaced, while the later grid nodes are distributed according to the parabola distribution law, and the grid node spacing gradually increases. In the Y direction, the grid in the central region is dense, and the upper and lower walls are sparse. The shock waves reflected from the wall surface are numerically absorbed and dissipated to a certain extent. In the Z direction, uniform mesh is adopted. (2) Boundary conditions 1) Inlet: Pressure (p), density (ρ ) and velocity (u , v , w) were given at the inlet of primary flow and secondary flow. Random disturbances of a certain frequency are added to the five basic parameters at the entrance to simulate a fully developed turbulence in a limited computational domain.
ϕi = A sin(ωt + ξ ) where ϕi is the variable, A the disturbance amplitude, ω the disturbance frequency, and ξ is the random phase within [0, 2π ]. Perturbation frequency is determined on the basis of the spectral information in the mixing layer obtained from a previous calculation [31]. The primary and secondary flow averages use a tangent function to simulate the velocity gradient in the mixing layer given the boundary layer formed at the separator.
2) Outlet: At the exit, the supersonic region and subsonic region are judged according to the local Mach number. Zero gradient boundary conditions are adopted for the supersonic region. The density of the subsonic region is given by the backpressure, and the remaining parameters are given by zero gradient boundary conditions. 3) The boundary conditions on the top and bottom surface: Each parameter is given a specific value, for example, v , w set to 0. 4) The boundary conditions in the spanwise direction: Periodic boundary conditions are adopted as follows:
90 91 92 93 94 95 96 97 98 99 100
ϕz=0 = ϕz=Lz
101
where ϕ is a variable, representing the pressure, velocity and other parameters in the calculation.
103
102 104 105
3. Calculation result analysis
106 107
3.1. Parameter field distribution
108 109
In the calculation results, at a certain moment, the calculation domain is cut off along the X Y plane, and the Mach number distribution contours of three working conditions from D1 to D3, as shown in Fig. 3, is obtained. It can be clearly found that the large-scale coherent structure in the shear mixing flow forms the small-scale vortex structure early, and the flow is dominated by the small-scale and broken vortex structure, which is quite different from the incompressible flow and the compressible flow with small convective Mach number. In addition, as the convection Mach number increases in turn, flow transition delay, and due to the large dynamic pressure difference between primary and secondary streams, the shear flow direction has a certain deflection in the secondary flow area. The density gradient contours can be obtained by calculating the gradient of the density field and taking the model, as shown in Fig. 4. The change of flow structure in the flow field can be clearly observed by using dark color display in the area where the density gradient varies greatly. As can be seen from the above two groups of images, the development process of compressible shear mixing flow is mainly divided into three stages: laminar flow zone, transition zone and full development zone. The fluid in laminar flow zone develops stably, and there is a linear interface between the two flows,
110 111 112 113 114 115 116 117 118 119 120 121 122 123 124 125 126 127 128 129 130 131 132
JID:AESCTE
4
AID:105431 /FLA
[m5G; v1.261; Prn:3/10/2019; 11:20] P.4 (1-18)
L. Yang et al. / Aerospace Science and Technology ••• (••••) ••••••
1
67
2
68
3
69
4
70
5
71
6
72
7
73
8
74
9
75
10
76
11
77
12
78
13
79
14
80
15
81
16
82
17
83
18
84
19
85
20
86
21
87
22
88
23
89
24
90
25
91
26
92
27
93
28
94
29
95
30
96
31
97
32
98
33
99
34
100
35
101
36
102
37
103 104
38 39
Fig. 4. Density gradient contours of D1–D3. (For interpretation of the colors in the figure(s), the reader is referred to the web version of this article.)
41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60 61 62 63 64 65 66
105 106
40
which presents a two-dimensional characteristic; Transition zone is a transition zone from two-dimensional to three-dimensional. The coherent structure with a slightly larger scale can be observed and begins to present three-dimensional flow characteristics; The full development zone presents the characteristics of full threedimensional flow, without large-scale vortex structure, and the flow is dominated by small-scale and broken vortex structure. The length of each stage of the mixing layer under each working condition is marked by a red dotted line in Fig. 4. It can be found that with the increase of the convection Mach number, the transition position of the mixing layer begins to be delayed, and the transition zone length is shortened. In order to observe the development characteristics of shear mixing flow in each stage more clearly, each zone was magnified and displayed respectively. Since the structure of laminar flow zone is relatively simple, it is no longer shown here, and only the detailed flow structure of transition zone and full development zone is shown here. (1) Transition zone The Mach number distribution contours and density gradient contours of the mixing layer in the transition zone are shown in Fig. 5. After transition from laminar flow state of mixing layer, coherent structure similar to that caused by K-H instability is formed. After a relatively short distance development, large-scale coherent structure is gradually deformed and broken, and relatively chaotic small-scale vortex structure appears. By comparing the images un-
der different working conditions, we can also find that as the convection Mach number increases, the transition area decreases, and the coherent structure becomes less and less obvious. By example D3, large-scale coherent structure can hardly be seen. (2) Full development zone The Mach number distribution contours and density gradient contours of the mixing layer in the full development zone are shown in Fig. 6. In this region, the flow structure of the mixing layer becomes more chaotic, and large-scale coherent structure is no longer seen, only small-scale vortex structure close to dissipative state can be observed, and the flow structure appears chaotic and disorderly.
107 108 109 110 111 112 113 114 115 116 117 118 119
3.2. Mixing layer thickness
120 121
Using the average velocity thicknessδ = y U 2−0.1 U − y U 2+0.9 U as the characterization method of the mixing layer thickness, the distribution curve of the thickness of the shear layer along the streamwise is obtained, which is shown in Fig. 7(a). At the beginning of the flow, the shear layer is laminar, the thickness of the shear layer grows very slowly in the flow direction, and the development of the shear layer thickness under three different conditions is very close to 75 mm. At the streamwise of about 100 mm, the shear layer under three working conditions began to turn, and the growth rate of the shear layer began to increase. Especially, the larger the convective Mach number, the later the transition
122 123 124 125 126 127 128 129 130 131 132
JID:AESCTE AID:105431 /FLA
[m5G; v1.261; Prn:3/10/2019; 11:20] P.5 (1-18)
L. Yang et al. / Aerospace Science and Technology ••• (••••) ••••••
5
1
67
2
68
3
69
4
70
5
71
6
72
7
73
8
74
9
75
10
76
11
77
12
78
13
79
14
80
15
81
16
82
17
83
18
84
19
85
20
86
21
87
22
88
23
89
24
90
25
91
26
92
27
93
28
94
29
95
30
96
31
97
32
Fig. 5. Mach number distribution contours and density gradient contours of transition zone.
98
33
99
34
100
35
101
36
102
37
103
38
104
39
105
40
106
41
107
42
108
43
109
44
110
45
111
46
112
47
113
48
114
49
115
50
116
51
117
52
118
53
119
54
120
55
121
56
122
57
123
58
124
59
125
60
126
61
127
62
128
63
129
64
130 131
65 66
Fig. 6. Mach number distribution contours and density gradient contours of full development zone.
132
JID:AESCTE
6
AID:105431 /FLA
[m5G; v1.261; Prn:3/10/2019; 11:20] P.6 (1-18)
L. Yang et al. / Aerospace Science and Technology ••• (••••) ••••••
1
67
2
68
3
69
4
70
5
71
6
72
7
73
8
74
9
75
10
76
11
77
12
78
13
79
14
80
15
81
16
82
17
83
18
84
19
85
20
86
21
87
22
88
23
89
24
90
25
91
26
92
27
93
28
94
29
95
30
96
31
97
32
98 99
33 34
Fig. 7. Thickness and thickness growth rate of mixing layer of D1–D3.
100
35
101
36
102
37
103
38
104
39
105
40
106
41
107
42
108
43
109
44
110
45
111
46
112
47
113
48
114
49
115
50
116
51
117
52
118
53
119
54
120
55
121
56
122
57
123
58
124
59
125
60
126
61
127
62
128
63
129
64
130 131
65 66
Fig. 8. Isocontours of ∇ · U = −0.01 and pressure distribution contours.
132
JID:AESCTE AID:105431 /FLA
[m5G; v1.261; Prn:3/10/2019; 11:20] P.7 (1-18)
L. Yang et al. / Aerospace Science and Technology ••• (••••) ••••••
7
1
67
2
68
3
69
4
70
5
71
6
72
7
73
8
74
9
75
10
76
11
77
12
78
13
Fig. 9. Isocontours of ∇ · U = −0.01 and locally enlarged contours of pressure distribution.
14
79 80
15
81
16
82
17
83
18
84
19
85
20
86
21
87
22
88
23
89
24
90
25
91
26
92
27
93
28
94
29
95
30
96
31
97
32
98
33
99
34
100
35
101
36
102
37
103
38
104
39
105
40
106 107
41
Fig. 10. Q-criterion iso-surface contours colored by Ma of D1–D3.
42 44 45 46 47 48 49 50 51 52 53 54 55
position. The figure shows that the lengths of the shear layer transitions under different convective Mach numbers are different. The thickness variation of the transition zone can be regarded as a linear growth trend. There is some irregularity in the thickness of the shear layer in the development zone, but Overall, the subsonicsupersonic shear layer thickness still shows an approximate linear change after the transition zone. Therefore, in the shear layer transition zone and the development zone, a plurality of streamwise shear layer thicknesses are extracted for linear fitting and the shear layer thickness growth rate is calculated as shown in Fig. 7(b, c). The calculated thickness growth of the dimensionless shear layer decreases as the convective Mach number increases.
56 57
3.3. Shocklets structure
58 59 60 61 62 63 64 65 66
108 109
43
By analyzing the density gradient contour of the shear mixing flow, the shocklets structure is found in the shear mixing flow. With the development of the mixing layer, the high-speed rotation of the large-scale vortex structure further reduces the pressure in the center of the vortex core, so that the tangential velocity of the outer edge of the vortex is relative to the velocity of the supersonic flow in the upper layer, and a shocklets is formed near the outer edge of the vortex in the supersonic flow. Kida [20] first discovered
the shocklets structure in the numerical simulation of turbulent flows, and There is no unified and effective method to display the structure of shocklets now. Here, the identification method of velocity divergence proposed by Freund et al. [21,30] is adopted to find out where the divergence value in the flow field is less than a certain negative value. Because the fluid particle by shock wave, its density has a leap, accordingly, it should be negative infinity divergence theory, and numerical calculation for a larger negative, here take ∇ · U = −0.01, as shown in Fig. 8 black strip curve. It can be found from Fig. 8 that the shocklets structure extracted with ∇ · U as the criterion mainly occurs at the stage when the mixing layer has not been fully developed, and is generated by the coherent structure of the mixing layer and propagates to the upper boundary of the flow. In this region, the mixing layer has not been fully developed and there are few disturbances to the primary and secondary streams, so the shocklets structure can be observed clearly. In the zone of full development, shocklets structure still exists in the primary stream, but the shocklets structure at this time is different from the shocklets structure in the transition zone, and there are many disturbances in the primary and secondary streams, so it is difficult to produce a clear shocklets structure.
110 111 112 113 114 115 116 117 118 119 120 121 122 123 124 125 126 127 128 129 130 131 132
JID:AESCTE
8
AID:105431 /FLA
[m5G; v1.261; Prn:3/10/2019; 11:20] P.8 (1-18)
L. Yang et al. / Aerospace Science and Technology ••• (••••) ••••••
1
67
2
68
3
69
4
70
5
71
6
72
7
73
8
74
9
75
10
76
11
77
12
78
13
79
14
80
15
81
16
82
17
83
18
84
19
85
20
86
21
87
22
88
23
89
24
90
25
91
26
92
27
93
28
94
29
95
30
96
31
97
32
98
33
99
34
100
35
101
36
102
37
103
38
104
39
105
40
106
41
107
42
108
43
109
44
110
45
111
46
112
47
113
48
114
49
115
50
116
51
117
52
118
53
119 120
54 55
Fig. 11. Vorticity contours in the downstream plane and the
ωz contours in the spanwise plane of D1–D3.
57 58 59 60 61 62 63 64 65 66
121 122
56
The starting point of the shocklets structure is located near the vortex core, which presents as a local minimum in the pressure distribution. In the study conducted by Zhou [22] (Mc = 0.7), shocklets were generated after the vortex core, while in the study conducted by Rossmann [23] (Mc = 1.7), shocklets were generated before the vortex core. This was because the convection Mach number was higher in Rossmann’s study, and the strong disturbance caused the early generation of shocklets structure. It can be seen from the ∇ · U iso-surface contours and the locally enlarged contours of the pressure distribution contours in
Fig. 9, with the increase of the convective Mach number, the location of the shocklets move before the vortex core, which is consistent with the above research results.
123 124 125 126
3.4. Fine flow structures
127 128
The coherent structure of turbulence is a large-scale vortex structure which is relatively organized. Its scale and structure have certain universality and repeatability for a certain type of flow. Research on the evolution process of the coherent structure can help
129 130 131 132
JID:AESCTE AID:105431 /FLA
[m5G; v1.261; Prn:3/10/2019; 11:20] P.9 (1-18)
L. Yang et al. / Aerospace Science and Technology ••• (••••) ••••••
9
1
67
2
68
3
69
4
70
5
71
6
72
7
73
8
74
9
75
10
76
11
77
12
78
13
79
14
80
15
81
16
82
17
83
18
84
19
85
20
86
21
87
22
88
23
89 90
24
Fig. 11. (continued)
25
91
26
92
27
93
28
94
29
95
30
96
31
97
32
98
33
99
34
100
35
101
36
102
37
103
38
104
39
105
40
106
41
107
42
108
43
109
44
110
45
111
46
112
47
113
48
114
49
115
50
116
51
117
52
118
53
119
54
120
55
121
56
122
57
123
58
124
59
125
60
126
61
127
62
128
63
129
64
130 131
65 66
Fig. 12.
ωz contours in the spanwise plane of D1–D3.
132
JID:AESCTE
10
AID:105431 /FLA
[m5G; v1.261; Prn:3/10/2019; 11:20] P.10 (1-18)
L. Yang et al. / Aerospace Science and Technology ••• (••••) ••••••
1
67
2
68
3
69
4
70
5
71
6
72
7
73
8
74
9
75
10
76
11
77
12
78
13
79
14
80
15
81
16
82
17
83
18
84
19
85
20
86
21
87
22
88
23
89
24
90
25
91
26
92
27
93
28
94
29
95
30
96
31
97
32
98
33
99
34
100
35
101
36
102
37
103
38
104
39
105
40
106
41
107
42
108
43
109
44
110
45
111 112
46 47
Fig. 13. Top view and oblique view Λ vortex structure of D1–D3.
49 50 51 52 53 54 55 56 57 58 59 60 61 62 63 64 65 66
113 114
48
us have a deeper understanding of the development process of transition from layer to turbulent state and the flow mechanism of transition and turbulence, so as to provide a basis for the research on turbulence control, turbulence model and large eddy simulation method. Similar to the incompressible and compressible supersonicsupersonic mixing layer, the compressible subsonic-supersonic mixing layer also presents a linear growth state, as shown in Fig. 10, the Q-criterion iso-surface contours with Mach number staining. It can be clearly found that with the increase of convective Mach number, the compressibility increases, the thickness of mixing layer becomes thinner, the transition position is delayed, and the large-scale vortex structure decreases in the flow, while the small-scale vortex structure increases. Fig. 11 shows the vorticity contours in the downstream plane and the ωz contours in the spanwise plane. It can be found that when two flows meet, the first is laminar flow state. At this time,
there is only viscous shear action between two flows, but no exchange of momentum and energy. After a certain distance, the two flows become unstable. At this time, the expression of the flow is similar to the coherent structure caused by K-H instability in the incompressible shear mixing flow, presenting a flow form dominated by two-dimensional state. With the gradual development of the coherent structure, the size increases gradually, and then the large-scale vortex structure begins to break up, forming small-scale vortex structure, until it completely breaks down and reaches the fully mixing flow state. Fig. 12 shows the ωz contours in the plane with different downstream directions. For the working conditions of different convection Mach Numbers, the transition position of the mixing layer is different. Before the mixing layer begins transition, it is generally believed that the flow at this time is dominated by the two-dimensional state. It can be found from Fig. 12(a) that the flow structure in this part is not a strict two-dimensional state and
115 116 117 118 119 120 121 122 123 124 125 126 127 128 129 130 131 132
JID:AESCTE AID:105431 /FLA
[m5G; v1.261; Prn:3/10/2019; 11:20] P.11 (1-18)
L. Yang et al. / Aerospace Science and Technology ••• (••••) ••••••
11
Mach number increases, the Λ vortex structure is gradually elongated. With vortex structure further development, after Λ vortex structure, the vortex structure becomes more clutter. A U-shaped vortex structure, commonly known as hairpin vortex, can be observed in this zone, as shown in Fig. 14.
1 2 3 4 5 6
3.5. Analysis of proper orthogonal decomposition
8
The amount of data in the DNS calculation is huge. Only 49 snapshots are taken here for POD modal decomposition. The POD modal analysis is performed on the velocity U in the X direction parallel to the XY plane, and the eigenvalue λ of each mode is used to represent the energy ratio of each mode, and the POD modal energy and the cumulative energy curve of each order are obtained. The POD modal energy and the modal cumulative energy curve are obtained as shown in Fig. 15. It can be found from Fig. 15 that the first-order modal energy accounts for more than 99% of the total energy, and the proportion of the first-order modal energy in the three working conditions is 99.67%, 99.43%, and 99.56%, respectively. This indicates that the main flow in the compressible mixing flow is characterized by an average flow state, and the high-order mode has a small energy content, which can be used to describe its unsteady flow field characteristics more accurately. By performing FFT analysis on the modal coefficients of each order, the frequencies of the modes of each order can be obtained, as shown in Fig. 16 for the 2–5th mode frequency curves of D1–D3. It can be seen from the figure that the higher frequency of the modal energy is distributed between 0 and 30 kHz, but the dominant frequencies of the modes are not unique. Due to the small amount of calculated data selected, there may be large deviations. However, this provides a new idea for the enhanced blending of shear mixed flow, that is, the main frequency affecting the shear mixing flow can be obtained by FFT analysis of the modal coefficient, and then the corresponding frequency is applied to the velocity U in the X direction at the inlet boundary. The interference can make the mixing flow turn into a turbulent state as soon as possible.
10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38
4. Analysis of Dynamic Mode Decomposition (DMD)
40 42 43 44 45 46 47 48 49 50 51 52 53
Fig. 14. Hairpin vortex structure after Λ vortex structure of the D1–D3.
55 56
60 61 62 63 64 65 66
71 72 74 76 77 78 79 80 81 82 83 84 85 86 87 88 89 90 91 92 93 94 95 96 97 98 99 100 101 102 103 104 106 107
41
59
70
105
39
58
69
75
9
57
68
73
7
54
67
has been differentiated in the spanwise direction. In Fig. 12(b), an obvious spanwise structure has already appeared. By Fig. 12(c), the flow start to become chaotic, and by Fig. 12(d), the flow structure has become completely chaotic. In order to display the Λ vortex structure in the development process of mixing layer, the Q-criterion contours stained with Mach number was locally enlarged, as shown in Fig. 13. As can be seen from Fig. 13, when shear mixing flow begins to transition, a two-dimensional coherent structure similar to that caused by K-H instability is developed, i.e., Λ vortex structure. As the convective
It is very important for shear mixing flow to accurately describe and understand its flow structure change. However, it is very difficult to accurately and quantitatively extract the complex flow structure, and it requires inverse and eigenvalue calculation of the large matrix. DMD method is a flow field decomposition method developed on the basis of linear Koopman mapping. The DMD method can obtain the main structure of the flow field, directly obtain the corresponding frequency of each mode, and judge the stability of each mode. The DMD method is only based on the snapshot of the flow field and is not limited by the flow type. The DMD method is only based on the snapshot of the flow field and is not limited by the flow type. The decomposition of the flow field is based on the dynamic characteristics of the flow field. So the extracted bases are orthogonal to each other in spatiotemporal evolution, and the main spatiotemporal characteristics of the flow field can be obtained. The DMD method was first proposed by Schmid et al. [24], and then the DMD method was used to analyze the flow field data obtained by numerical method and experiment respectively. Rowley et al. [25] used DMD method to analyze the nonlinear flow in the transverse jet flow, and the results showed that this method could extract the main flow structure and frequency. In recent years, a large number of scholars have continuously improved and optimized the DMD method, such as Optimized DMD (Opt-DMD) and Optimal MD (Opt-MD) [26]. In order to more
108 109 110 111 112 113 114 115 116 117 118 119 120 121 122 123 124 125 126 127 128 129 130 131 132
JID:AESCTE
AID:105431 /FLA
12
[m5G; v1.261; Prn:3/10/2019; 11:20] P.12 (1-18)
L. Yang et al. / Aerospace Science and Technology ••• (••••) ••••••
1
67
2
68
3
69
4
70
5
71
6
72
7
73
8
74
9
75
10
76
11
77
12
78
13
79
14
80
15
81
16
82
17
83
18
84
19
85
20
86
21
87
22
88
23
89
24
90
25
91
26
92
27
93
28
94
29
95
30
96
31
97
32
98
33
99
34
100
35
101
36
102
37
103
38
104
39
105
40
106
41
107
42
108
43
109
44
110
45
111
46
112
47
113
48
114 115
49
Fig. 15. POD modal energy curve.
50 52 53 54 55 56 57 58 59 60 61 62 63 64
accurately extract the characteristic modes that have significant influence on the development of the flow field, Jovanovic et al. [27] combined with the optimization theory formed the sparse improved dynamic modal decomposition (DMDSP) method. Currently, DMD and its improved methods have been widely used in complex flow problems such as shock wave boundary layer interference [28], boundary layer transition [29] and so on. Snapshots of n moments obtained by experiment or numerical simulation can be written as a snapshot sequence matrix X and Y, and the time interval between any two snapshots is t, that is
Assuming that the flow field v i +1 can be represented by linear mapping A ∈ R m×m and the flow field v i , namely
Y = [ A v 1 A v 2 · · · A v n −1 ] = A X
128
X = [ v 1 v 2 · · · v n −1 ]
For the matrix X with rank r, DMD method need to find matrix B ∈ C r ×r to approximate high dimension matrix A, then matrix B is the optimal low-dimensional estimation matrix of matrix A.
130
65 66
116 117
51
Y = [v 2 v 3 · · · vn]
118 119 120
v i +1 = A v i
121
The above assumption is true for any snapshot at any time. If the dynamic system itself is nonlinear, the process is a linear estimation process, that is, the nonlinear estimation is realized through the linear hypothesis. There are
123
122 124 125 126 127 129 131 132
JID:AESCTE AID:105431 /FLA
[m5G; v1.261; Prn:3/10/2019; 11:20] P.13 (1-18)
L. Yang et al. / Aerospace Science and Technology ••• (••••) ••••••
13
1
67
2
68
3
69
4
70
5
71
6
72
7
73
8
74
9
75
10
76
11
77
12
78
13
79
14
80
15
81
16
82
17
83
18
84
19
85
20
86
21
87
22
88
23
89
24
90
25
91
26
92
27
93
28
94
29
95
30
96
31
97
32 34 35 36 37 38
DMD decomposition results can be obtained by solving the eigenvalues and eigenvectors of B. The matrix B can be obtained by the singular value decomposition of the matrix X , i.e. H
39
X = UΣV
40
A ≈ U BU H
41 42 43 44 45 46 47 48 49 50 51 52
Carry out singular value decomposition for matrix X , then: U ∈ C , U H U = I , U H is the complex conjugate transpose matrix of the matrix U , I is the unit matrix; V ∈ C r × N ; Σ is diagonal matrix, and diagonal elements from top to bottom are singular values from large to small. The calculation of matrix B can be summarized as the following minimization problem: M ×r
2 min B Y − U B Σ V H F where, · 2F is Frobenius norm. Then, the optimal form of matrix B is:
53 54 55 56 57 58 59 60 61 62 63 64 65 66
98
Fig. 16. D1–D3 2–5th mode frequency.
33
B D M D = U H Y V Σ −1
99
It is assumed that the matrix B is an accurate mapping matrix in the r-dimensional subspace, rather than an estimation of the matrix A, that is:
100 101 102 103 104
xk+1 = Bxk
105
Then the matrix U can be seen as the approximate mapping from xk to v k , namely:
106
v k ≈ U xk
109
In order to show the contribution of modal to flow field, the modal amplitude is defined as αi = z iH x0 . If α j is bigger, so that the corresponding DMD modal for the effect of initial condition selected is bigger, is:
111
107 108 110 112 113 114 115
v k ≈ U xk =
r
ϕ i μki α j
i =1
For modal energy ϕ i , the following matrix norm is often used for characterization:
Solve the eigenvalues μi and eigenvector z i of the matrix B D M D , then the ith dynamic mode is:
ϕ i =
ϕ i = U zi
4.1. DMD analysis of X direction velocity U in XY plane
Modal magnification g i = Re(lg μi / t ) and modal frequency ωi = Im(lg μi / t ) can be obtained by logarithmic form of eigenvalues μi . When judging the stability, if the magnification is positive, the corresponding mode is unstable; if the magnification is negative, the corresponding mode is stable; if the magnification is zero, the corresponding mode is periodic; if the eigenvalue falls within the unit circle, it represents the stable mode, and vice versa.
Because of the large scale of direct numerical simulation results, only 49 snapshots were taken for DMD modal analysis. The DMD modal analysis of the X -direction velocity U in the X Y plane obtained the eigenvalue distribution as well as the modal energy and frequency relation as shown in Fig. 17, and part of the DMD modal contours of D3 as shown in Fig. 16.
116 117 118 119 120 121
1
122
V H Σ −1 Z
123 124 125 126 127 128 129 130 131 132
JID:AESCTE
14
AID:105431 /FLA
[m5G; v1.261; Prn:3/10/2019; 11:20] P.14 (1-18)
L. Yang et al. / Aerospace Science and Technology ••• (••••) ••••••
1
67
2
68
3
69
4
70
5
71
6
72
7
73
8
74
9
75
10
76
11
77
12
78
13
79
14
80
15
81
16
82
17
83
18
84
19
85
20
86
21
87
22
88
23
89
24
90
25
91
26
92
27
93
28
94
29
95
30
96
31
97
32
98
33
99
34
100
35
101
36
102
37
103
38
104
39
105
40
106
41
107
42
108
43
109
44
110
45
111
46
112
47
113
48
114
49
115
50
116
51
117
52
118 119
53 54
Fig. 17. Modal eigenvalue distribution and modal energy versus frequency of D1–D3.
56 57 58 59 60 61 62 63 64 65 66
120 121
55
Fig. 17(a) shows the eigenvalue distribution of each order of DMD mode, where the horizontal axis is the real part of the modal eigenvalue, and the vertical axis corresponds to the imaginary part. As shown in the figure, the eigenvalues are basically near the unit circle, and individual modes are located in the unit circle, which indicates that the calculated modes are all stable. The 1 mode energy values of the three operating conditions are 5921.5, 5349.8, and 5680.9, respectively, and their mode frequencies are all zero. As can be seen from Fig. 17(b), the energy of 1 order modal is at least 3 orders of magnitude larger than that of other
modal, indicating that 1 order modal is the main modal of the flow field. Fig. 18(a) is 1 order modal contour, which is consistent with the average velocity field of velocity U in the XY plane, indicating that the DMD method can capture the main flow characteristics of the mixing layer. As can be seen from Fig. 18, 2 order modal contour is the same as 3 order modal contour, which is because the modes obtained by DMD decomposition are conjugate modes, that is, the two modes corresponding to the two eigenvalues of conjugate to each other are actually one mode.
122 123 124 125 126 127 128 129 130 131 132
JID:AESCTE AID:105431 /FLA
[m5G; v1.261; Prn:3/10/2019; 11:20] P.15 (1-18)
L. Yang et al. / Aerospace Science and Technology ••• (••••) ••••••
15
1
67
2
68
3
69
4
70
5
71
6
72
7
73
8
74
9
75
10
76
11
77
12
78
13
79
14
80
15
81
16
82
17
83
18
84
19
85
20
86
21
87
22
88
23
89
24
90
25
91
26
92
27
93
28
94
29
95
30
96
31
97
32
98
33
99
34
100
35
101
36
102
37
103
38
104
39
105
40
106
41
107
42
108
43
109
44
110
45
111 112
46 47
Fig. 18. Part of the DMD modal contours of D3.
113
48
114
49
115
50
4.2. DMD analysis of Y direction velocity V in XY plane
51
that the DMD method can capture the main flow characteristics of the mixing layer.
116
5. Conclusion
119
52 53 54 55 56 57 58 59 60 61 62 63 64 65 66
The DMD modal analysis of the Y -direction velocity V in the X Y plane obtained the eigenvalue distribution as well as the modal energy and frequency relation as shown in Fig. 19, and part of the DMD modal contours of D3 as shown in Fig. 18. Fig. 19(a) shows the eigenvalue distribution of each order of DMD mode. As shown in the figure, the eigenvalues are basically near the unit circle, and individual modes are located in the unit circle, which indicates that the calculated modes are all stable. As can be seen from Fig. 19(b), the energy of 1 order modal is only equal to the order of other modal. All modes have higher energy, and the frequency distribution of higher modal energy is between 0–20 kHz. Fig. 20(a) is 1 order modal contour, which is consistent with the average velocity field of velocity V in the X Y plane, indicating
117 118 120
In this paper, the flow structure of the subsonic-supersonic mixing layer is obtained by direct numerical simulation with high order accuracy. The results show that the flow structure of the subsonic-supersonic mixing layer has the following characteristics: (1) The large-scale structure in the subsonic-supersonic shear mixing flow forms the small-scale structure early, and the flow is dominated by small-scale and broken vortex structure, which is quite different from the incompressible flow and the compressible flow with small convective Mach number. (2) With the increase of convection Mach number, the growth rate of dimensionless mixing layer thickness decreases, transition delay, and transition zone length is shortened. Moreover, due to
121 122 123 124 125 126 127 128 129 130 131 132
JID:AESCTE
16
AID:105431 /FLA
[m5G; v1.261; Prn:3/10/2019; 11:20] P.16 (1-18)
L. Yang et al. / Aerospace Science and Technology ••• (••••) ••••••
1
67
2
68
3
69
4
70
5
71
6
72
7
73
8
74
9
75
10
76
11
77
12
78
13
79
14
80
15
81
16
82
17
83
18
84
19
85
20
86
21
87
22
88
23
89
24
90
25
91
26
92
27
93
28
94
29
95
30
96
31
97
32
98
33
99
34
100
35
101
36
102
37
103
38
104
39
105
40
106
41
107
42
108
43
109
44
110
45
111
46
112
47
113
48
114
49
115
50
116 117
51 52
Fig. 19. Modal eigenvalue distribution and modal energy versus frequency of D1–D3.
54 55 56 57 58 59 60 61 62 63 64 65 66
118 119
53
the large dynamic pressure difference between primary and secondary streams, the shear flow direction has a certain deflection in the secondary stream zone. (3) The development process of compressible shear mixing flow is mainly divided into three stages: laminar flow zone, transition zone and full development zone. The fluid in laminar flow zone develops stably, and there is a linear interface between the two flows, which presents a two-dimensional characteristic; Transition zone is a transition zone from two-dimensional to three-dimensional. The coherent structure with a slightly larger scale can be observed and begins to present three-dimensional flow characteristics; The full development zone presents the characteristics of full three-
dimensional flow, without large-scale vortex structure, and the flow is dominated by small-scale and broken vortex structure. (4) Shocklets structure mainly occurs at the stage when the mixing layer has not been fully developed. It is generated by the coherent structure of the mixing layer and propagates to the upper boundary of the flow. The starting point of shocklets structure is located near the vortex core. As the convective Mach number increases, the location of the small shock wave moves towards the vortex core. (5) DMD analysis results show: The eigenvalues of DMD modes of each order are basically located in the unit circle, which indicates that the calculated modes of each order are stable; Only the
120 121 122 123 124 125 126 127 128 129 130 131 132
JID:AESCTE AID:105431 /FLA
[m5G; v1.261; Prn:3/10/2019; 11:20] P.17 (1-18)
L. Yang et al. / Aerospace Science and Technology ••• (••••) ••••••
17
1
67
2
68
3
69
4
70
5
71
6
72
7
73
8
74
9
75
10
76
11
77
12
78
13
79
14
80
15
81
16
82
17
83
18
84
19
85
20
86
21
87
22
88
23
89
24
90
25
91
26
92
27
93
28
94
29
95
30
96
31
97
32
98
33
99
34
100
35
101
36
102
37
103
38
104
39
105
40
106
41
107 108
42 43
Fig. 20. Part of the DMD modal contours of D3.
45 46 47 48 49 50 51 52 53 54 55 56 57 58
modal energy of the velocity U in the X direction is at least 3 orders of magnitude greater than that of other modes, which is the main mode of the flow field. The magnitude of the mode energy of the velocity in other directions is not different, and there is no major mode; The mode obtained by DMD decomposition is conjugate mode, that is, the two modes corresponding to the two eigenvalues of conjugate to each other are actually one mode, and the corresponding modal contour is the same; In addition, it is found that there is a certain dominant frequency in the flow structure, which can provide a reference for the active mixing enhancement method. Declaration of competing interest
59 60
None declared.
61 62
Acknowledgements
63 64 65 66
109 110
44
This research was supported by the National Natural Science Foundation of China (Grant No. 51506178) and the Aeronautical Science Foundation of China (Grant No. 2015ZA53009).
References
111 112
[1] D.W. Bogdanoff, Compressibility effects in turbulent shear layers, AIAA J. 21 (6) (1983) 926–927. [2] J.L. Hall, P.E. Dimotakis, H. Rosemann, Experiments in nonreacting compressible shear layers, AIAA J. 31 (12) (1993) 2247–2254. [3] C. Pantano, S. Sarkar, A study of compressible effects in the high-speed turbulent shear layer using direct simulation, J. Fluid Mech. 451 (2002) 329–371. [4] G.L. Brown, A. Roshko, On density effects and large structure in turbulent mixing layers, J. Fluid Mech. 64 (1974) 775–814. [5] D. Papamoschou, A. Roshko, Observations of supersonic free shear layers, ¯ Sadhan a¯ Acad. Proc. Eng. Sci. 12 (1988) 1–14. [6] T. Rossmann, An Experimental Investigation of High-Compressibility Mixing Layers, Stanford University, USA, 2002. [7] T.J. Wang, J. Chen, X.T. Shi, et al., Experimental evidence for non-linear growth in compressible mixing layer, Sci. China, Phys. Mech. Astron. 57 (5) (2014) 963–970. [8] C.W.H. Van Doorne, J. Westerweel, Measurement of laminar, transitional and turbulent pipe flow using stereoscopic-PIV, Exp. Fluids 42 (2) (2007) 259–279. [9] M.G. Olsen, J.C. Dutton, Planar velocity measurements in a weakly compressible mixing layer, J. Fluid Mech. 431 (2003) 189–222. [10] N.T. Clemens, S.P. Petullo, D.S. Dolling, Large-scale structure evolution in supersonic interacting shear layers, AIAA J. 34 (10) (1996) 2062–2070. [11] L. Zhou, R. Zhao, W. Yuan, An investigation of interface conditions inherent in detached-eddy simulation methods, Aerosp. Sci. Technol. 74 (2018) 46–55.
113 114 115 116 117 118 119 120 121 122 123 124 125 126 127 128 129 130 131 132
JID:AESCTE
18
1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21
AID:105431 /FLA
[m5G; v1.261; Prn:3/10/2019; 11:20] P.18 (1-18)
L. Yang et al. / Aerospace Science and Technology ••• (••••) ••••••
[12] L. Wang, S. Fu, Detached-eddy simulation of flow past a pitching NACA 0015 airfoil with pulsed actuation, Aerosp. Sci. Technol. 69 (2017) 123–135. [13] Y. Han, Y. He, Y. Tian, F. Zhong, J. Le, Assessment of a hybrid RANS/LES simulation method and URANS method in depicting the unsteady motions of flow structures in a scramjet combustor, Aerosp. Sci. Technol. 72 (2018) 114–122. [14] L. Doris, C. Teuaud, L.T. Phnoc, LES of spatially developing 3D compressible mixing layer, Comput. Fluid Solid Mech. 328 (7) (2000) 567–573. [15] S. Fu, Q. Li, Numerical simulation of compressible mixing layers, Int. J. Heat Fluid Flow 27 (5) (2006) 895–901. [16] K. Arora, K. Chakravarthy, Large eddy simulation of supersonic, compressible, turbulent mixing layers, Aerosp. Sci. Technol. 86 (2019) 592–598. [17] Sandham, Effects of compressibility and shock-wave interactions on turbulent shear flows, Flow Turbul. Combust. 97 (1) (2016) 1–25. [18] X. Fang, C. Shen, M. Sun, et al., Effects of oblique shock waves on turbulent structures and statistics of supersonic mixing layers, Phys. Fluids 30 (11) (2018). [19] Xinliang Li, Dexun Fu, Yanwen Ma, Direct numerical simulation of hypersonic boundary layer transition over a blunt cone, AIAA J. 46 (11) (2008) 2899–2913. [20] S. Kida, S.A. Orszag, Enstrophy budget in decaying compressible turbulence, J. Sci. Comput. 5 (1) (1990) 1–34. [21] J.B. Freund, P. Moin, S.K. Lele, Compressibility effects in a turbulent annular mixing layer. Part 2. Mixing of a passive scalar, J. Fluid Mech. 421 (2000) 269–292.
[22] Q. Zhou, F. He, M.Y. Shen, Direct numerical simulation of a spatially developing compressible plane mixing layer: flow structures and mean flow properties, J. Fluid Mech. 711 (2012) 437–468. [23] T. Rossmann, M.G. Mungal, R.K. Hanson, Evolution and growth of large-scale structures in high compressibility mixing layers, J. Turbul. 3 (2002) N9. [24] P.J. Schmid, Dynamic mode decomposition of numerical and experimental data, J. Fluid Mech. 656 (2010) 5–28. ´ S. Bagheri, et al., Spectral analysis of nonlinear flows, J. [25] C.W. Rowley, I. Mezic, Fluid Mech. 641 (2009) 115–127. [26] A. Wynn, D.S. Pearson, et al., Optimal mode decomposition for unsteady flows, J. Fluid Mech. 733 (2) (2013) 473–503. ´ P.J. Schmid, J.W. Nichols, Sparsity-promoting dynamic mode de[27] M.R. Jovanovic, composition, Phys. Fluids 26 (2) (2014) 024103. [28] M. Grilli, P.J. Schmid, S. Hickel, N.A. Adams, Analysis of unsteady behaviour in shockwave turbulent boundary layer interaction, J. Fluid Mech. 700 (6) (2012) 16–28. [29] T. Sayadi, P.J. Schmid, J.W. Nichols, P. Moin, Reduced-order representation of near-wall structures in the late transitional boundary layer, J. Fluid Mech. 748 (2) (2014) 278–301. [30] J.B. Freund, S.K. Lele, P. Moin, Compressibility effects in a turbulent annular mixing layer. Part 1. Turbulence and growth rate, J. Fluid Mech. 421 (2000) 229–267. [31] Zhang Liu, et al., Direct numerical simulation of subsonic-supersonic mixing layer, Acta Astronaut. 153 (2018) 50–59.
67 68 69 70 71 72 73 74 75 76 77 78 79 80 81 82 83 84 85 86 87
22
88
23
89
24
90
25
91
26
92
27
93
28
94
29
95
30
96
31
97
32
98
33
99
34
100
35
101
36
102
37
103
38
104
39
105
40
106
41
107
42
108
43
109
44
110
45
111
46
112
47
113
48
114
49
115
50
116
51
117
52
118
53
119
54
120
55
121
56
122
57
123
58
124
59
125
60
126
61
127
62
128
63
129
64
130
65
131
66
132
[m5G; v1.261; Prn:3/10/2019; 11:20] P.1 (1-18)
Aerospace Science and Technology ••• (••••) ••••••
1
Contents lists available at ScienceDirect
67 68
2 3
Aerospace Science and Technology
4
69 70 71
5
72
6
www.elsevier.com/locate/aescte
7
73
8
74
9
75
10
76
11 12 13
Direct numerical simulation of fine flow structures of subsonic-supersonic mixing layer
77 78 79 80
14 15 16 17 18 19 20
Liu Yang
a,∗
a
, Ma Dong , Fu Benshuai
b,c
a
, Li Qiang , Zhang Chenxi
d
a
Science and Technology on Combustion, Internal Flow and Thermo-Structure Laboratory, Northwestern Polytechnical University, Xi’an, China b No. 713 Institute, China Shipbuilding Industry Co., Zhengzhou 450015, PR China c Henan Key Laboratory of Underwater Intelligence Equipment, Zhengzhou 450015, PR China d School of Astronautics, Northwestern Polytechnical University, Xi’an, China
a r t i c l e
i n f o
25 26 27 28 29 30 31 32
83 84 85 86 88
a b s t r a c t
89
23 24
82
87
21 22
81
Article history: Received 29 July 2019 Received in revised form 4 September 2019 Accepted 25 September 2019 Available online xxxx Keywords: Subsonic-supersonic mixing layer Compressibility Direct numerical simulation Dynamic modal decomposition
33 34 35 36 37 38 39 40 41 42 43 44 45 46
Subsonic-supersonic mixing flow is one of the important research fields in turbulence research. For example, in the combustion chamber of the rocket ramjet combination engine, the rocket jet flow and the air inflow is a typical shear mixing flow, which has the characteristics of high convection Mach number and large flow parameter gradient. It is of great significance to explore the development law and flow structure of subsonic-supersonic mixing layer to enhance the mixing and thus improve the performance of rocket ramjet engine. In this paper, in order to further study the evolution process and obtain the evolution mechanism and fine flow structure of subsonic-supersonic mixing layer, three groups of high-order format direct numerical simulations of typical working conditions are carried out. The convective Mach numbers are 0.69, 0.92 and 1.27, respectively. The direct numerical simulation results show that: (1) the large-scale structure in the subsonic-supersonic shear mixing flow formed the small-scale structure earlier, and the flow was dominated by the small-scale and broken eddy structure; (2) the development process of the subsonic-supersonic shear mixing flow was mainly divided into three stages: laminar flow zone, transition zone and development zone. The transition zone was a transition zone from two-dimensional to three-dimensional, which began to show three-dimensional flow characteristics. The vorticity shows a twisted-winding structure. The large-scale Coherent structure, the Λ eddy and hairpin eddy, can be observed, and the Λ eddy structure is gradually elongated with the increase of convective Mach number; (3) the shocklets structure mainly occurs in the stage when the mixing layer has not fully developed. The starting point of the shocklets structure is near the vortex core. With the increase of convective Mach number, the position of the shocklets generation moves to the front of the vortex core; (4) the dynamic mode decomposition results show that the eigenvalues of each mode are basically located in the unit circle, which indicates that the calculated modes are stable; in addition, it is found that there is a certain dominant frequency in the flow structure, which can provide a reference for the active mixing enhancement method. © 2019 Elsevier Masson SAS. All rights reserved.
90 91 92 93 94 95 96 97 98 99 100 101 102 103 104 105 106 107 108 109 110 111 112
47
113
48
114
49
115 116
50 51
1. Introduction
52 53 54 55 56 57 58 59
In rocket ramjet engines, the rocket gas with high-temperature and high-speed and air flow with low-temperature and low-speed will form a Subsonic-supersonic shear flow with a large parameter gradient. Due to the high compressibility of the shear flow, the growth rate of the shear layer is small, the turbulence and the Reynolds stress are small, the stability of the shear layer is enhanced, and it is not easy to form a large-scale quasi-ordered
60 61 62 63 64 65 66
*
Corresponding author. E-mail address: [email protected] (L. Yang).
https://doi.org/10.1016/j.ast.2019.105431 1270-9638/© 2019 Elsevier Masson SAS. All rights reserved.
structure, leading to difficulties in the transfer of energy and momentum in the two flows, which directly affects the engine’s blending and combustion efficiency. The remarkable feature of the subsonic-supersonic shear mixed flow is its strong compressibility. The shear mixing flow was first studied on the low-speed incompressible flow blending problem and then gradually extended to high-speed compressible flow mixing. Bogdanoff et al. [1] proposed the concept of convective Mach number (Mc), which describes the effect of compressibility on the mixing layer by convection Mach number. It has been studied that the compressibility of the flow can reduce the growth rate of the shear layer and enhance the stability of flow [2,3].
117 118 119 120 121 122 123 124 125 126 127 128 129 130 131 132
JID:AESCTE
AID:105431 /FLA
2
1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48
[m5G; v1.261; Prn:3/10/2019; 11:20] P.2 (1-18)
L. Yang et al. / Aerospace Science and Technology ••• (••••) ••••••
In the free shear layer, due to Kelvin-Helmholtz (K-H) instability, the large-scale vortex structure is rolled up, paired, and fused as the flow field develops. These large-scale vortex structures still exist under supersonic flow conditions. First, Brown and Roshko et al. [4] demonstrated the existence of large-scale reverse-order structures of mixed layers in experiments at low-speed flows. Papamoschou and Roshko et al. [5] demonstrated the existence of large-scale structures in the compressible mixing layer. Later, many researchers began to capture the flow field structure by means of schlieren, shadow, plane laser Mie scattering (PLMS), particle image velocimetry (PIV) and other measurements [6–8]. Olsen et al. [9] used the schlieren technique to capture the shear flow field, and the outline of the shear layer was observed from the schlieren. Clemens and Petullo et al. [10] used optical measurements to observe the mixing layer. It was found that the vortex structure of the compressible mixing layer under the same convective Mach number is different from that of the incompressible mixing layer. The pulsation velocity obtained in the experiment revealed that the compressible shear mixed layer vortex structure has threedimensional characteristics. Experimental research on the fine structure of the shear layer has been difficult in recent decades. Numerical simulation method has become a very important means to study the fine flow of turbulent flow in shear layer. Numerical simulation method has become a very important means to study the precise flow of shear layer turbulence. Considerable recent progress has been made in the study of DES series of methods [11–13]. Doris et al. [14] studied the three-dimensional compressive shear mixing flow (Mc = 0.64) using the LES method. The results show that both the streamwise and the spanning pulsation velocity have selfsimilarity. Fu Song and Li Qibing et al. [15] researched the supersonic shear flow by the LES method and found that small shock wave structures appeared in the flow field, and the compressibility reduced the parameters such as Reynolds stress. Two-dimensional and three-dimensional simulations of a supersonic mixing layer are performed using an in-house hybrid LES code by Konark et al. [16]. They found that three-dimensional predictions of growth rate and statistical moments in the self similar regime are consistent with past experimental and direct numerical simulation studies. Recently, some researchers have studied the effects of shock and compressibility on the shear layer [17,18]. At present, researchers have conducted a lot of research on shear mixing flow and made great progress. However, further research is needed in the fine flow structure of the subsonicsupersonic shear mixed layer. Therefore, in this paper, the relevant parameters were changed to study the fine flow structure of the shear mixed layer, and the flow structure under different working conditions was compared.
49 50
2. Direct numerical simulation method
51 52 53 54 55 56 57 58
Through theoretical analysis, it is found that decreasing density and increasing temperature can reduce Reynolds number of the flow, thus greatly reducing calculated amount of direct numerical simulation. Therefore, this paper conducts direct numerical simulation research on the subsonic-supersonic shear mixing flow under the condition of low pressure and high temperature, and obtains the fine flow structure of the shear mixing flow.
59 60
2.1. Numerical computation method
61 62 63 64 65 66
2.1.1. Numerical computation platform In the direct numerical simulation study of this paper, the OpenCFD program developed by Li Xinliang is adopted [19]. The core of the program is the compressible solver of N-S equation. The numerical method of the program is the finite difference method,
67 68 69 70 71 72 73 74 75 76 77 78 79 80 81 82 83 84 85
Fig. 1. Turbulence energy spectrum analysis diagram.
86 87
which includes a variety of high-precision difference schemes currently popular, such as central difference scheme, upwind difference scheme and WENO difference scheme, etc., and also includes the difference scheme with its own characteristics, such as upwind compact scheme, group speed control difference scheme, etc. In the direct numerical simulation, the fourth order RungeKutta scheme is adopted for the time term, the sixth order central difference method is adopted for the viscous term, and the seventh order WENO method is adopted for the non-viscous term. The entire simulation process was completed on the Tianhe-2 supercomputer of the national supercomputing center in Guangzhou, using 15 nodes and 480 cores for calculation, and the calculation time of a single example was about 150 hours.
88 89 90 91 92 93 94 95 96 97 98 99 100 101
2.1.2. Numerical computation method validation More detailed experiments and direct numerical simulation studies have been carried out on the supersonic-supersonic shear mixing flow. In the study of this paper, the turbulence energy spectrum of direct numerical simulation results is analyzed, and the validity of grid scale is judged by the turbulence energy spectrum. Based on isotropic turbulence theory proposed by Kolmogorov, the slope of the turbulence energy spectrum in the inertial range is −5/3, and the absolute value in the dissipative zone is greater than −5/3. In the direct numerical simulation performed in this study, an intercepted mixed flow fully develops a section of the area. The position of the spectrum is selected in the streamwise range of 400–600. The resulting turbulent energy spectrum is shown in Fig. 1. The slope of the energy spectrum in the low-wavenumber region is −5/3, which satisfies Kolmogorov’s turbulence theory. In the high wavenumber region, a −5 slope indicates the minimum grid size of the direct numerical simulation observed in the dissipative region, thereby satisfying the grid requirements of the direct numerical simulation.
102 103 104 105 106 107 108 109 110 111 112 113 114 115 116 117 118 119 120 121 122
2.2. Numerical calculation parameters
123 124
2.2.1. Numerical calculation conditions All parameters are normalized by primary stream flow parameters. Given the density (d), velocity (U ), temperature (T ) and other parameters of the incoming flow, dimensionless flow parameters are shown in Table 1, and flow parameters reflecting the actual flow are shown in Table 2. In this chapter, three groups of direct numerical simulation calculations were carried out, and the convective Mach Numbers were 0.69, 0.92 and 1.19, respectively. The
125 126 127 128 129 130 131 132
JID:AESCTE AID:105431 /FLA
[m5G; v1.261; Prn:3/10/2019; 11:20] P.3 (1-18)
L. Yang et al. / Aerospace Science and Technology ••• (••••) ••••••
1 2
Table 1 Dimensionless flow parameters.
3
Primary stream
67 68
6 7 8
U T d Mc
69
Secondary stream
4 5
3
1 1 1 –
D1
D2
D3
0.199 1.768 0.565 0.69
0.133 0.786 1.273 0.92
0.111 1.222 0.818 1.19
70 71 72 73 74 75
9 10 11
76
Table 2 Actual flow parameters.
77 78
12
Primary stream
13 14 15 16 17 18 19
U , m/s T *, K T, K P , kPa Ma Mc
Secondary stream
79
D1
D2
D3
D1
80
845.2 800.0 444.4 10 2.0 0.69
1267.9 1800.0 1000.0 10 2.0 0.92
1524.9 1800.0 642.9 10 3.0 1.19
168.6 800.0 785.9 10 0.3 –
81 82 83 84 85
20
86
21
87
22
Fig. 3. Mach number distribution contours of D1–D3.
24 25 26 27 28 29
Fig. 2. Schematic diagram of computational area.
30 31 32
corresponding dimensionless Reynolds Numbers were 1346, 523 and 1309, respectively.
33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60 61 62 63 64 65 66
88 89
23
2.2.2. Mesh and boundary conditions (1) Mesh The finite-difference method is difficult to deal with complex geometry, so the classical cubic configuration in the plane shear mixing flow is adopted here, as shown in Fig. 2. The reference length is 5 × 10−4 m. The dimensionless length and the number of grid nodes in X , Y and Z directions are 600, 120, 60 and 3500, 500, 300, respectively, and the total amount of grid is 525 million. In the direction of X , the first 90% grid nodes are equally spaced, while the later grid nodes are distributed according to the parabola distribution law, and the grid node spacing gradually increases. In the Y direction, the grid in the central region is dense, and the upper and lower walls are sparse. The shock waves reflected from the wall surface are numerically absorbed and dissipated to a certain extent. In the Z direction, uniform mesh is adopted. (2) Boundary conditions 1) Inlet: Pressure (p), density (ρ ) and velocity (u , v , w) were given at the inlet of primary flow and secondary flow. Random disturbances of a certain frequency are added to the five basic parameters at the entrance to simulate a fully developed turbulence in a limited computational domain.
ϕi = A sin(ωt + ξ ) where ϕi is the variable, A the disturbance amplitude, ω the disturbance frequency, and ξ is the random phase within [0, 2π ]. Perturbation frequency is determined on the basis of the spectral information in the mixing layer obtained from a previous calculation [31]. The primary and secondary flow averages use a tangent function to simulate the velocity gradient in the mixing layer given the boundary layer formed at the separator.
2) Outlet: At the exit, the supersonic region and subsonic region are judged according to the local Mach number. Zero gradient boundary conditions are adopted for the supersonic region. The density of the subsonic region is given by the backpressure, and the remaining parameters are given by zero gradient boundary conditions. 3) The boundary conditions on the top and bottom surface: Each parameter is given a specific value, for example, v , w set to 0. 4) The boundary conditions in the spanwise direction: Periodic boundary conditions are adopted as follows:
90 91 92 93 94 95 96 97 98 99 100
ϕz=0 = ϕz=Lz
101
where ϕ is a variable, representing the pressure, velocity and other parameters in the calculation.
103
102 104 105
3. Calculation result analysis
106 107
3.1. Parameter field distribution
108 109
In the calculation results, at a certain moment, the calculation domain is cut off along the X Y plane, and the Mach number distribution contours of three working conditions from D1 to D3, as shown in Fig. 3, is obtained. It can be clearly found that the large-scale coherent structure in the shear mixing flow forms the small-scale vortex structure early, and the flow is dominated by the small-scale and broken vortex structure, which is quite different from the incompressible flow and the compressible flow with small convective Mach number. In addition, as the convection Mach number increases in turn, flow transition delay, and due to the large dynamic pressure difference between primary and secondary streams, the shear flow direction has a certain deflection in the secondary flow area. The density gradient contours can be obtained by calculating the gradient of the density field and taking the model, as shown in Fig. 4. The change of flow structure in the flow field can be clearly observed by using dark color display in the area where the density gradient varies greatly. As can be seen from the above two groups of images, the development process of compressible shear mixing flow is mainly divided into three stages: laminar flow zone, transition zone and full development zone. The fluid in laminar flow zone develops stably, and there is a linear interface between the two flows,
110 111 112 113 114 115 116 117 118 119 120 121 122 123 124 125 126 127 128 129 130 131 132
JID:AESCTE
4
AID:105431 /FLA
[m5G; v1.261; Prn:3/10/2019; 11:20] P.4 (1-18)
L. Yang et al. / Aerospace Science and Technology ••• (••••) ••••••
1
67
2
68
3
69
4
70
5
71
6
72
7
73
8
74
9
75
10
76
11
77
12
78
13
79
14
80
15
81
16
82
17
83
18
84
19
85
20
86
21
87
22
88
23
89
24
90
25
91
26
92
27
93
28
94
29
95
30
96
31
97
32
98
33
99
34
100
35
101
36
102
37
103 104
38 39
Fig. 4. Density gradient contours of D1–D3. (For interpretation of the colors in the figure(s), the reader is referred to the web version of this article.)
41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60 61 62 63 64 65 66
105 106
40
which presents a two-dimensional characteristic; Transition zone is a transition zone from two-dimensional to three-dimensional. The coherent structure with a slightly larger scale can be observed and begins to present three-dimensional flow characteristics; The full development zone presents the characteristics of full threedimensional flow, without large-scale vortex structure, and the flow is dominated by small-scale and broken vortex structure. The length of each stage of the mixing layer under each working condition is marked by a red dotted line in Fig. 4. It can be found that with the increase of the convection Mach number, the transition position of the mixing layer begins to be delayed, and the transition zone length is shortened. In order to observe the development characteristics of shear mixing flow in each stage more clearly, each zone was magnified and displayed respectively. Since the structure of laminar flow zone is relatively simple, it is no longer shown here, and only the detailed flow structure of transition zone and full development zone is shown here. (1) Transition zone The Mach number distribution contours and density gradient contours of the mixing layer in the transition zone are shown in Fig. 5. After transition from laminar flow state of mixing layer, coherent structure similar to that caused by K-H instability is formed. After a relatively short distance development, large-scale coherent structure is gradually deformed and broken, and relatively chaotic small-scale vortex structure appears. By comparing the images un-
der different working conditions, we can also find that as the convection Mach number increases, the transition area decreases, and the coherent structure becomes less and less obvious. By example D3, large-scale coherent structure can hardly be seen. (2) Full development zone The Mach number distribution contours and density gradient contours of the mixing layer in the full development zone are shown in Fig. 6. In this region, the flow structure of the mixing layer becomes more chaotic, and large-scale coherent structure is no longer seen, only small-scale vortex structure close to dissipative state can be observed, and the flow structure appears chaotic and disorderly.
107 108 109 110 111 112 113 114 115 116 117 118 119
3.2. Mixing layer thickness
120 121
Using the average velocity thicknessδ = y U 2−0.1 U − y U 2+0.9 U as the characterization method of the mixing layer thickness, the distribution curve of the thickness of the shear layer along the streamwise is obtained, which is shown in Fig. 7(a). At the beginning of the flow, the shear layer is laminar, the thickness of the shear layer grows very slowly in the flow direction, and the development of the shear layer thickness under three different conditions is very close to 75 mm. At the streamwise of about 100 mm, the shear layer under three working conditions began to turn, and the growth rate of the shear layer began to increase. Especially, the larger the convective Mach number, the later the transition
122 123 124 125 126 127 128 129 130 131 132
JID:AESCTE AID:105431 /FLA
[m5G; v1.261; Prn:3/10/2019; 11:20] P.5 (1-18)
L. Yang et al. / Aerospace Science and Technology ••• (••••) ••••••
5
1
67
2
68
3
69
4
70
5
71
6
72
7
73
8
74
9
75
10
76
11
77
12
78
13
79
14
80
15
81
16
82
17
83
18
84
19
85
20
86
21
87
22
88
23
89
24
90
25
91
26
92
27
93
28
94
29
95
30
96
31
97
32
Fig. 5. Mach number distribution contours and density gradient contours of transition zone.
98
33
99
34
100
35
101
36
102
37
103
38
104
39
105
40
106
41
107
42
108
43
109
44
110
45
111
46
112
47
113
48
114
49
115
50
116
51
117
52
118
53
119
54
120
55
121
56
122
57
123
58
124
59
125
60
126
61
127
62
128
63
129
64
130 131
65 66
Fig. 6. Mach number distribution contours and density gradient contours of full development zone.
132
JID:AESCTE
6
AID:105431 /FLA
[m5G; v1.261; Prn:3/10/2019; 11:20] P.6 (1-18)
L. Yang et al. / Aerospace Science and Technology ••• (••••) ••••••
1
67
2
68
3
69
4
70
5
71
6
72
7
73
8
74
9
75
10
76
11
77
12
78
13
79
14
80
15
81
16
82
17
83
18
84
19
85
20
86
21
87
22
88
23
89
24
90
25
91
26
92
27
93
28
94
29
95
30
96
31
97
32
98 99
33 34
Fig. 7. Thickness and thickness growth rate of mixing layer of D1–D3.
100
35
101
36
102
37
103
38
104
39
105
40
106
41
107
42
108
43
109
44
110
45
111
46
112
47
113
48
114
49
115
50
116
51
117
52
118
53
119
54
120
55
121
56
122
57
123
58
124
59
125
60
126
61
127
62
128
63
129
64
130 131
65 66
Fig. 8. Isocontours of ∇ · U = −0.01 and pressure distribution contours.
132
JID:AESCTE AID:105431 /FLA
[m5G; v1.261; Prn:3/10/2019; 11:20] P.7 (1-18)
L. Yang et al. / Aerospace Science and Technology ••• (••••) ••••••
7
1
67
2
68
3
69
4
70
5
71
6
72
7
73
8
74
9
75
10
76
11
77
12
78
13
Fig. 9. Isocontours of ∇ · U = −0.01 and locally enlarged contours of pressure distribution.
14
79 80
15
81
16
82
17
83
18
84
19
85
20
86
21
87
22
88
23
89
24
90
25
91
26
92
27
93
28
94
29
95
30
96
31
97
32
98
33
99
34
100
35
101
36
102
37
103
38
104
39
105
40
106 107
41
Fig. 10. Q-criterion iso-surface contours colored by Ma of D1–D3.
42 44 45 46 47 48 49 50 51 52 53 54 55
position. The figure shows that the lengths of the shear layer transitions under different convective Mach numbers are different. The thickness variation of the transition zone can be regarded as a linear growth trend. There is some irregularity in the thickness of the shear layer in the development zone, but Overall, the subsonicsupersonic shear layer thickness still shows an approximate linear change after the transition zone. Therefore, in the shear layer transition zone and the development zone, a plurality of streamwise shear layer thicknesses are extracted for linear fitting and the shear layer thickness growth rate is calculated as shown in Fig. 7(b, c). The calculated thickness growth of the dimensionless shear layer decreases as the convective Mach number increases.
56 57
3.3. Shocklets structure
58 59 60 61 62 63 64 65 66
108 109
43
By analyzing the density gradient contour of the shear mixing flow, the shocklets structure is found in the shear mixing flow. With the development of the mixing layer, the high-speed rotation of the large-scale vortex structure further reduces the pressure in the center of the vortex core, so that the tangential velocity of the outer edge of the vortex is relative to the velocity of the supersonic flow in the upper layer, and a shocklets is formed near the outer edge of the vortex in the supersonic flow. Kida [20] first discovered
the shocklets structure in the numerical simulation of turbulent flows, and There is no unified and effective method to display the structure of shocklets now. Here, the identification method of velocity divergence proposed by Freund et al. [21,30] is adopted to find out where the divergence value in the flow field is less than a certain negative value. Because the fluid particle by shock wave, its density has a leap, accordingly, it should be negative infinity divergence theory, and numerical calculation for a larger negative, here take ∇ · U = −0.01, as shown in Fig. 8 black strip curve. It can be found from Fig. 8 that the shocklets structure extracted with ∇ · U as the criterion mainly occurs at the stage when the mixing layer has not been fully developed, and is generated by the coherent structure of the mixing layer and propagates to the upper boundary of the flow. In this region, the mixing layer has not been fully developed and there are few disturbances to the primary and secondary streams, so the shocklets structure can be observed clearly. In the zone of full development, shocklets structure still exists in the primary stream, but the shocklets structure at this time is different from the shocklets structure in the transition zone, and there are many disturbances in the primary and secondary streams, so it is difficult to produce a clear shocklets structure.
110 111 112 113 114 115 116 117 118 119 120 121 122 123 124 125 126 127 128 129 130 131 132
JID:AESCTE
8
AID:105431 /FLA
[m5G; v1.261; Prn:3/10/2019; 11:20] P.8 (1-18)
L. Yang et al. / Aerospace Science and Technology ••• (••••) ••••••
1
67
2
68
3
69
4
70
5
71
6
72
7
73
8
74
9
75
10
76
11
77
12
78
13
79
14
80
15
81
16
82
17
83
18
84
19
85
20
86
21
87
22
88
23
89
24
90
25
91
26
92
27
93
28
94
29
95
30
96
31
97
32
98
33
99
34
100
35
101
36
102
37
103
38
104
39
105
40
106
41
107
42
108
43
109
44
110
45
111
46
112
47
113
48
114
49
115
50
116
51
117
52
118
53
119 120
54 55
Fig. 11. Vorticity contours in the downstream plane and the
ωz contours in the spanwise plane of D1–D3.
57 58 59 60 61 62 63 64 65 66
121 122
56
The starting point of the shocklets structure is located near the vortex core, which presents as a local minimum in the pressure distribution. In the study conducted by Zhou [22] (Mc = 0.7), shocklets were generated after the vortex core, while in the study conducted by Rossmann [23] (Mc = 1.7), shocklets were generated before the vortex core. This was because the convection Mach number was higher in Rossmann’s study, and the strong disturbance caused the early generation of shocklets structure. It can be seen from the ∇ · U iso-surface contours and the locally enlarged contours of the pressure distribution contours in
Fig. 9, with the increase of the convective Mach number, the location of the shocklets move before the vortex core, which is consistent with the above research results.
123 124 125 126
3.4. Fine flow structures
127 128
The coherent structure of turbulence is a large-scale vortex structure which is relatively organized. Its scale and structure have certain universality and repeatability for a certain type of flow. Research on the evolution process of the coherent structure can help
129 130 131 132
JID:AESCTE AID:105431 /FLA
[m5G; v1.261; Prn:3/10/2019; 11:20] P.9 (1-18)
L. Yang et al. / Aerospace Science and Technology ••• (••••) ••••••
9
1
67
2
68
3
69
4
70
5
71
6
72
7
73
8
74
9
75
10
76
11
77
12
78
13
79
14
80
15
81
16
82
17
83
18
84
19
85
20
86
21
87
22
88
23
89 90
24
Fig. 11. (continued)
25
91
26
92
27
93
28
94
29
95
30
96
31
97
32
98
33
99
34
100
35
101
36
102
37
103
38
104
39
105
40
106
41
107
42
108
43
109
44
110
45
111
46
112
47
113
48
114
49
115
50
116
51
117
52
118
53
119
54
120
55
121
56
122
57
123
58
124
59
125
60
126
61
127
62
128
63
129
64
130 131
65 66
Fig. 12.
ωz contours in the spanwise plane of D1–D3.
132
JID:AESCTE
10
AID:105431 /FLA
[m5G; v1.261; Prn:3/10/2019; 11:20] P.10 (1-18)
L. Yang et al. / Aerospace Science and Technology ••• (••••) ••••••
1
67
2
68
3
69
4
70
5
71
6
72
7
73
8
74
9
75
10
76
11
77
12
78
13
79
14
80
15
81
16
82
17
83
18
84
19
85
20
86
21
87
22
88
23
89
24
90
25
91
26
92
27
93
28
94
29
95
30
96
31
97
32
98
33
99
34
100
35
101
36
102
37
103
38
104
39
105
40
106
41
107
42
108
43
109
44
110
45
111 112
46 47
Fig. 13. Top view and oblique view Λ vortex structure of D1–D3.
49 50 51 52 53 54 55 56 57 58 59 60 61 62 63 64 65 66
113 114
48
us have a deeper understanding of the development process of transition from layer to turbulent state and the flow mechanism of transition and turbulence, so as to provide a basis for the research on turbulence control, turbulence model and large eddy simulation method. Similar to the incompressible and compressible supersonicsupersonic mixing layer, the compressible subsonic-supersonic mixing layer also presents a linear growth state, as shown in Fig. 10, the Q-criterion iso-surface contours with Mach number staining. It can be clearly found that with the increase of convective Mach number, the compressibility increases, the thickness of mixing layer becomes thinner, the transition position is delayed, and the large-scale vortex structure decreases in the flow, while the small-scale vortex structure increases. Fig. 11 shows the vorticity contours in the downstream plane and the ωz contours in the spanwise plane. It can be found that when two flows meet, the first is laminar flow state. At this time,
there is only viscous shear action between two flows, but no exchange of momentum and energy. After a certain distance, the two flows become unstable. At this time, the expression of the flow is similar to the coherent structure caused by K-H instability in the incompressible shear mixing flow, presenting a flow form dominated by two-dimensional state. With the gradual development of the coherent structure, the size increases gradually, and then the large-scale vortex structure begins to break up, forming small-scale vortex structure, until it completely breaks down and reaches the fully mixing flow state. Fig. 12 shows the ωz contours in the plane with different downstream directions. For the working conditions of different convection Mach Numbers, the transition position of the mixing layer is different. Before the mixing layer begins transition, it is generally believed that the flow at this time is dominated by the two-dimensional state. It can be found from Fig. 12(a) that the flow structure in this part is not a strict two-dimensional state and
115 116 117 118 119 120 121 122 123 124 125 126 127 128 129 130 131 132
JID:AESCTE AID:105431 /FLA
[m5G; v1.261; Prn:3/10/2019; 11:20] P.11 (1-18)
L. Yang et al. / Aerospace Science and Technology ••• (••••) ••••••
11
Mach number increases, the Λ vortex structure is gradually elongated. With vortex structure further development, after Λ vortex structure, the vortex structure becomes more clutter. A U-shaped vortex structure, commonly known as hairpin vortex, can be observed in this zone, as shown in Fig. 14.
1 2 3 4 5 6
3.5. Analysis of proper orthogonal decomposition
8
The amount of data in the DNS calculation is huge. Only 49 snapshots are taken here for POD modal decomposition. The POD modal analysis is performed on the velocity U in the X direction parallel to the XY plane, and the eigenvalue λ of each mode is used to represent the energy ratio of each mode, and the POD modal energy and the cumulative energy curve of each order are obtained. The POD modal energy and the modal cumulative energy curve are obtained as shown in Fig. 15. It can be found from Fig. 15 that the first-order modal energy accounts for more than 99% of the total energy, and the proportion of the first-order modal energy in the three working conditions is 99.67%, 99.43%, and 99.56%, respectively. This indicates that the main flow in the compressible mixing flow is characterized by an average flow state, and the high-order mode has a small energy content, which can be used to describe its unsteady flow field characteristics more accurately. By performing FFT analysis on the modal coefficients of each order, the frequencies of the modes of each order can be obtained, as shown in Fig. 16 for the 2–5th mode frequency curves of D1–D3. It can be seen from the figure that the higher frequency of the modal energy is distributed between 0 and 30 kHz, but the dominant frequencies of the modes are not unique. Due to the small amount of calculated data selected, there may be large deviations. However, this provides a new idea for the enhanced blending of shear mixed flow, that is, the main frequency affecting the shear mixing flow can be obtained by FFT analysis of the modal coefficient, and then the corresponding frequency is applied to the velocity U in the X direction at the inlet boundary. The interference can make the mixing flow turn into a turbulent state as soon as possible.
10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38
4. Analysis of Dynamic Mode Decomposition (DMD)
40 42 43 44 45 46 47 48 49 50 51 52 53
Fig. 14. Hairpin vortex structure after Λ vortex structure of the D1–D3.
55 56
60 61 62 63 64 65 66
71 72 74 76 77 78 79 80 81 82 83 84 85 86 87 88 89 90 91 92 93 94 95 96 97 98 99 100 101 102 103 104 106 107
41
59
70
105
39
58
69
75
9
57
68
73
7
54
67
has been differentiated in the spanwise direction. In Fig. 12(b), an obvious spanwise structure has already appeared. By Fig. 12(c), the flow start to become chaotic, and by Fig. 12(d), the flow structure has become completely chaotic. In order to display the Λ vortex structure in the development process of mixing layer, the Q-criterion contours stained with Mach number was locally enlarged, as shown in Fig. 13. As can be seen from Fig. 13, when shear mixing flow begins to transition, a two-dimensional coherent structure similar to that caused by K-H instability is developed, i.e., Λ vortex structure. As the convective
It is very important for shear mixing flow to accurately describe and understand its flow structure change. However, it is very difficult to accurately and quantitatively extract the complex flow structure, and it requires inverse and eigenvalue calculation of the large matrix. DMD method is a flow field decomposition method developed on the basis of linear Koopman mapping. The DMD method can obtain the main structure of the flow field, directly obtain the corresponding frequency of each mode, and judge the stability of each mode. The DMD method is only based on the snapshot of the flow field and is not limited by the flow type. The DMD method is only based on the snapshot of the flow field and is not limited by the flow type. The decomposition of the flow field is based on the dynamic characteristics of the flow field. So the extracted bases are orthogonal to each other in spatiotemporal evolution, and the main spatiotemporal characteristics of the flow field can be obtained. The DMD method was first proposed by Schmid et al. [24], and then the DMD method was used to analyze the flow field data obtained by numerical method and experiment respectively. Rowley et al. [25] used DMD method to analyze the nonlinear flow in the transverse jet flow, and the results showed that this method could extract the main flow structure and frequency. In recent years, a large number of scholars have continuously improved and optimized the DMD method, such as Optimized DMD (Opt-DMD) and Optimal MD (Opt-MD) [26]. In order to more
108 109 110 111 112 113 114 115 116 117 118 119 120 121 122 123 124 125 126 127 128 129 130 131 132
JID:AESCTE
AID:105431 /FLA
12
[m5G; v1.261; Prn:3/10/2019; 11:20] P.12 (1-18)
L. Yang et al. / Aerospace Science and Technology ••• (••••) ••••••
1
67
2
68
3
69
4
70
5
71
6
72
7
73
8
74
9
75
10
76
11
77
12
78
13
79
14
80
15
81
16
82
17
83
18
84
19
85
20
86
21
87
22
88
23
89
24
90
25
91
26
92
27
93
28
94
29
95
30
96
31
97
32
98
33
99
34
100
35
101
36
102
37
103
38
104
39
105
40
106
41
107
42
108
43
109
44
110
45
111
46
112
47
113
48
114 115
49
Fig. 15. POD modal energy curve.
50 52 53 54 55 56 57 58 59 60 61 62 63 64
accurately extract the characteristic modes that have significant influence on the development of the flow field, Jovanovic et al. [27] combined with the optimization theory formed the sparse improved dynamic modal decomposition (DMDSP) method. Currently, DMD and its improved methods have been widely used in complex flow problems such as shock wave boundary layer interference [28], boundary layer transition [29] and so on. Snapshots of n moments obtained by experiment or numerical simulation can be written as a snapshot sequence matrix X and Y, and the time interval between any two snapshots is t, that is
Assuming that the flow field v i +1 can be represented by linear mapping A ∈ R m×m and the flow field v i , namely
Y = [ A v 1 A v 2 · · · A v n −1 ] = A X
128
X = [ v 1 v 2 · · · v n −1 ]
For the matrix X with rank r, DMD method need to find matrix B ∈ C r ×r to approximate high dimension matrix A, then matrix B is the optimal low-dimensional estimation matrix of matrix A.
130
65 66
116 117
51
Y = [v 2 v 3 · · · vn]
118 119 120
v i +1 = A v i
121
The above assumption is true for any snapshot at any time. If the dynamic system itself is nonlinear, the process is a linear estimation process, that is, the nonlinear estimation is realized through the linear hypothesis. There are
123
122 124 125 126 127 129 131 132
JID:AESCTE AID:105431 /FLA
[m5G; v1.261; Prn:3/10/2019; 11:20] P.13 (1-18)
L. Yang et al. / Aerospace Science and Technology ••• (••••) ••••••
13
1
67
2
68
3
69
4
70
5
71
6
72
7
73
8
74
9
75
10
76
11
77
12
78
13
79
14
80
15
81
16
82
17
83
18
84
19
85
20
86
21
87
22
88
23
89
24
90
25
91
26
92
27
93
28
94
29
95
30
96
31
97
32 34 35 36 37 38
DMD decomposition results can be obtained by solving the eigenvalues and eigenvectors of B. The matrix B can be obtained by the singular value decomposition of the matrix X , i.e. H
39
X = UΣV
40
A ≈ U BU H
41 42 43 44 45 46 47 48 49 50 51 52
Carry out singular value decomposition for matrix X , then: U ∈ C , U H U = I , U H is the complex conjugate transpose matrix of the matrix U , I is the unit matrix; V ∈ C r × N ; Σ is diagonal matrix, and diagonal elements from top to bottom are singular values from large to small. The calculation of matrix B can be summarized as the following minimization problem: M ×r
2 min B Y − U B Σ V H F where, · 2F is Frobenius norm. Then, the optimal form of matrix B is:
53 54 55 56 57 58 59 60 61 62 63 64 65 66
98
Fig. 16. D1–D3 2–5th mode frequency.
33
B D M D = U H Y V Σ −1
99
It is assumed that the matrix B is an accurate mapping matrix in the r-dimensional subspace, rather than an estimation of the matrix A, that is:
100 101 102 103 104
xk+1 = Bxk
105
Then the matrix U can be seen as the approximate mapping from xk to v k , namely:
106
v k ≈ U xk
109
In order to show the contribution of modal to flow field, the modal amplitude is defined as αi = z iH x0 . If α j is bigger, so that the corresponding DMD modal for the effect of initial condition selected is bigger, is:
111
107 108 110 112 113 114 115
v k ≈ U xk =
r
ϕ i μki α j
i =1
For modal energy ϕ i , the following matrix norm is often used for characterization:
Solve the eigenvalues μi and eigenvector z i of the matrix B D M D , then the ith dynamic mode is:
ϕ i =
ϕ i = U zi
4.1. DMD analysis of X direction velocity U in XY plane
Modal magnification g i = Re(lg μi / t ) and modal frequency ωi = Im(lg μi / t ) can be obtained by logarithmic form of eigenvalues μi . When judging the stability, if the magnification is positive, the corresponding mode is unstable; if the magnification is negative, the corresponding mode is stable; if the magnification is zero, the corresponding mode is periodic; if the eigenvalue falls within the unit circle, it represents the stable mode, and vice versa.
Because of the large scale of direct numerical simulation results, only 49 snapshots were taken for DMD modal analysis. The DMD modal analysis of the X -direction velocity U in the X Y plane obtained the eigenvalue distribution as well as the modal energy and frequency relation as shown in Fig. 17, and part of the DMD modal contours of D3 as shown in Fig. 16.
116 117 118 119 120 121
1
122
V H Σ −1 Z
123 124 125 126 127 128 129 130 131 132
JID:AESCTE
14
AID:105431 /FLA
[m5G; v1.261; Prn:3/10/2019; 11:20] P.14 (1-18)
L. Yang et al. / Aerospace Science and Technology ••• (••••) ••••••
1
67
2
68
3
69
4
70
5
71
6
72
7
73
8
74
9
75
10
76
11
77
12
78
13
79
14
80
15
81
16
82
17
83
18
84
19
85
20
86
21
87
22
88
23
89
24
90
25
91
26
92
27
93
28
94
29
95
30
96
31
97
32
98
33
99
34
100
35
101
36
102
37
103
38
104
39
105
40
106
41
107
42
108
43
109
44
110
45
111
46
112
47
113
48
114
49
115
50
116
51
117
52
118 119
53 54
Fig. 17. Modal eigenvalue distribution and modal energy versus frequency of D1–D3.
56 57 58 59 60 61 62 63 64 65 66
120 121
55
Fig. 17(a) shows the eigenvalue distribution of each order of DMD mode, where the horizontal axis is the real part of the modal eigenvalue, and the vertical axis corresponds to the imaginary part. As shown in the figure, the eigenvalues are basically near the unit circle, and individual modes are located in the unit circle, which indicates that the calculated modes are all stable. The 1 mode energy values of the three operating conditions are 5921.5, 5349.8, and 5680.9, respectively, and their mode frequencies are all zero. As can be seen from Fig. 17(b), the energy of 1 order modal is at least 3 orders of magnitude larger than that of other
modal, indicating that 1 order modal is the main modal of the flow field. Fig. 18(a) is 1 order modal contour, which is consistent with the average velocity field of velocity U in the XY plane, indicating that the DMD method can capture the main flow characteristics of the mixing layer. As can be seen from Fig. 18, 2 order modal contour is the same as 3 order modal contour, which is because the modes obtained by DMD decomposition are conjugate modes, that is, the two modes corresponding to the two eigenvalues of conjugate to each other are actually one mode.
122 123 124 125 126 127 128 129 130 131 132
JID:AESCTE AID:105431 /FLA
[m5G; v1.261; Prn:3/10/2019; 11:20] P.15 (1-18)
L. Yang et al. / Aerospace Science and Technology ••• (••••) ••••••
15
1
67
2
68
3
69
4
70
5
71
6
72
7
73
8
74
9
75
10
76
11
77
12
78
13
79
14
80
15
81
16
82
17
83
18
84
19
85
20
86
21
87
22
88
23
89
24
90
25
91
26
92
27
93
28
94
29
95
30
96
31
97
32
98
33
99
34
100
35
101
36
102
37
103
38
104
39
105
40
106
41
107
42
108
43
109
44
110
45
111 112
46 47
Fig. 18. Part of the DMD modal contours of D3.
113
48
114
49
115
50
4.2. DMD analysis of Y direction velocity V in XY plane
51
that the DMD method can capture the main flow characteristics of the mixing layer.
116
5. Conclusion
119
52 53 54 55 56 57 58 59 60 61 62 63 64 65 66
The DMD modal analysis of the Y -direction velocity V in the X Y plane obtained the eigenvalue distribution as well as the modal energy and frequency relation as shown in Fig. 19, and part of the DMD modal contours of D3 as shown in Fig. 18. Fig. 19(a) shows the eigenvalue distribution of each order of DMD mode. As shown in the figure, the eigenvalues are basically near the unit circle, and individual modes are located in the unit circle, which indicates that the calculated modes are all stable. As can be seen from Fig. 19(b), the energy of 1 order modal is only equal to the order of other modal. All modes have higher energy, and the frequency distribution of higher modal energy is between 0–20 kHz. Fig. 20(a) is 1 order modal contour, which is consistent with the average velocity field of velocity V in the X Y plane, indicating
117 118 120
In this paper, the flow structure of the subsonic-supersonic mixing layer is obtained by direct numerical simulation with high order accuracy. The results show that the flow structure of the subsonic-supersonic mixing layer has the following characteristics: (1) The large-scale structure in the subsonic-supersonic shear mixing flow forms the small-scale structure early, and the flow is dominated by small-scale and broken vortex structure, which is quite different from the incompressible flow and the compressible flow with small convective Mach number. (2) With the increase of convection Mach number, the growth rate of dimensionless mixing layer thickness decreases, transition delay, and transition zone length is shortened. Moreover, due to
121 122 123 124 125 126 127 128 129 130 131 132
JID:AESCTE
16
AID:105431 /FLA
[m5G; v1.261; Prn:3/10/2019; 11:20] P.16 (1-18)
L. Yang et al. / Aerospace Science and Technology ••• (••••) ••••••
1
67
2
68
3
69
4
70
5
71
6
72
7
73
8
74
9
75
10
76
11
77
12
78
13
79
14
80
15
81
16
82
17
83
18
84
19
85
20
86
21
87
22
88
23
89
24
90
25
91
26
92
27
93
28
94
29
95
30
96
31
97
32
98
33
99
34
100
35
101
36
102
37
103
38
104
39
105
40
106
41
107
42
108
43
109
44
110
45
111
46
112
47
113
48
114
49
115
50
116 117
51 52
Fig. 19. Modal eigenvalue distribution and modal energy versus frequency of D1–D3.
54 55 56 57 58 59 60 61 62 63 64 65 66
118 119
53
the large dynamic pressure difference between primary and secondary streams, the shear flow direction has a certain deflection in the secondary stream zone. (3) The development process of compressible shear mixing flow is mainly divided into three stages: laminar flow zone, transition zone and full development zone. The fluid in laminar flow zone develops stably, and there is a linear interface between the two flows, which presents a two-dimensional characteristic; Transition zone is a transition zone from two-dimensional to three-dimensional. The coherent structure with a slightly larger scale can be observed and begins to present three-dimensional flow characteristics; The full development zone presents the characteristics of full three-
dimensional flow, without large-scale vortex structure, and the flow is dominated by small-scale and broken vortex structure. (4) Shocklets structure mainly occurs at the stage when the mixing layer has not been fully developed. It is generated by the coherent structure of the mixing layer and propagates to the upper boundary of the flow. The starting point of shocklets structure is located near the vortex core. As the convective Mach number increases, the location of the small shock wave moves towards the vortex core. (5) DMD analysis results show: The eigenvalues of DMD modes of each order are basically located in the unit circle, which indicates that the calculated modes of each order are stable; Only the
120 121 122 123 124 125 126 127 128 129 130 131 132
JID:AESCTE AID:105431 /FLA
[m5G; v1.261; Prn:3/10/2019; 11:20] P.17 (1-18)
L. Yang et al. / Aerospace Science and Technology ••• (••••) ••••••
17
1
67
2
68
3
69
4
70
5
71
6
72
7
73
8
74
9
75
10
76
11
77
12
78
13
79
14
80
15
81
16
82
17
83
18
84
19
85
20
86
21
87
22
88
23
89
24
90
25
91
26
92
27
93
28
94
29
95
30
96
31
97
32
98
33
99
34
100
35
101
36
102
37
103
38
104
39
105
40
106
41
107 108
42 43
Fig. 20. Part of the DMD modal contours of D3.
45 46 47 48 49 50 51 52 53 54 55 56 57 58
modal energy of the velocity U in the X direction is at least 3 orders of magnitude greater than that of other modes, which is the main mode of the flow field. The magnitude of the mode energy of the velocity in other directions is not different, and there is no major mode; The mode obtained by DMD decomposition is conjugate mode, that is, the two modes corresponding to the two eigenvalues of conjugate to each other are actually one mode, and the corresponding modal contour is the same; In addition, it is found that there is a certain dominant frequency in the flow structure, which can provide a reference for the active mixing enhancement method. Declaration of competing interest
59 60
None declared.
61 62
Acknowledgements
63 64 65 66
109 110
44
This research was supported by the National Natural Science Foundation of China (Grant No. 51506178) and the Aeronautical Science Foundation of China (Grant No. 2015ZA53009).
References
111 112
[1] D.W. Bogdanoff, Compressibility effects in turbulent shear layers, AIAA J. 21 (6) (1983) 926–927. [2] J.L. Hall, P.E. Dimotakis, H. Rosemann, Experiments in nonreacting compressible shear layers, AIAA J. 31 (12) (1993) 2247–2254. [3] C. Pantano, S. Sarkar, A study of compressible effects in the high-speed turbulent shear layer using direct simulation, J. Fluid Mech. 451 (2002) 329–371. [4] G.L. Brown, A. Roshko, On density effects and large structure in turbulent mixing layers, J. Fluid Mech. 64 (1974) 775–814. [5] D. Papamoschou, A. Roshko, Observations of supersonic free shear layers, ¯ Sadhan a¯ Acad. Proc. Eng. Sci. 12 (1988) 1–14. [6] T. Rossmann, An Experimental Investigation of High-Compressibility Mixing Layers, Stanford University, USA, 2002. [7] T.J. Wang, J. Chen, X.T. Shi, et al., Experimental evidence for non-linear growth in compressible mixing layer, Sci. China, Phys. Mech. Astron. 57 (5) (2014) 963–970. [8] C.W.H. Van Doorne, J. Westerweel, Measurement of laminar, transitional and turbulent pipe flow using stereoscopic-PIV, Exp. Fluids 42 (2) (2007) 259–279. [9] M.G. Olsen, J.C. Dutton, Planar velocity measurements in a weakly compressible mixing layer, J. Fluid Mech. 431 (2003) 189–222. [10] N.T. Clemens, S.P. Petullo, D.S. Dolling, Large-scale structure evolution in supersonic interacting shear layers, AIAA J. 34 (10) (1996) 2062–2070. [11] L. Zhou, R. Zhao, W. Yuan, An investigation of interface conditions inherent in detached-eddy simulation methods, Aerosp. Sci. Technol. 74 (2018) 46–55.
113 114 115 116 117 118 119 120 121 122 123 124 125 126 127 128 129 130 131 132
JID:AESCTE
18
1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21
AID:105431 /FLA
[m5G; v1.261; Prn:3/10/2019; 11:20] P.18 (1-18)
L. Yang et al. / Aerospace Science and Technology ••• (••••) ••••••
[12] L. Wang, S. Fu, Detached-eddy simulation of flow past a pitching NACA 0015 airfoil with pulsed actuation, Aerosp. Sci. Technol. 69 (2017) 123–135. [13] Y. Han, Y. He, Y. Tian, F. Zhong, J. Le, Assessment of a hybrid RANS/LES simulation method and URANS method in depicting the unsteady motions of flow structures in a scramjet combustor, Aerosp. Sci. Technol. 72 (2018) 114–122. [14] L. Doris, C. Teuaud, L.T. Phnoc, LES of spatially developing 3D compressible mixing layer, Comput. Fluid Solid Mech. 328 (7) (2000) 567–573. [15] S. Fu, Q. Li, Numerical simulation of compressible mixing layers, Int. J. Heat Fluid Flow 27 (5) (2006) 895–901. [16] K. Arora, K. Chakravarthy, Large eddy simulation of supersonic, compressible, turbulent mixing layers, Aerosp. Sci. Technol. 86 (2019) 592–598. [17] Sandham, Effects of compressibility and shock-wave interactions on turbulent shear flows, Flow Turbul. Combust. 97 (1) (2016) 1–25. [18] X. Fang, C. Shen, M. Sun, et al., Effects of oblique shock waves on turbulent structures and statistics of supersonic mixing layers, Phys. Fluids 30 (11) (2018). [19] Xinliang Li, Dexun Fu, Yanwen Ma, Direct numerical simulation of hypersonic boundary layer transition over a blunt cone, AIAA J. 46 (11) (2008) 2899–2913. [20] S. Kida, S.A. Orszag, Enstrophy budget in decaying compressible turbulence, J. Sci. Comput. 5 (1) (1990) 1–34. [21] J.B. Freund, P. Moin, S.K. Lele, Compressibility effects in a turbulent annular mixing layer. Part 2. Mixing of a passive scalar, J. Fluid Mech. 421 (2000) 269–292.
[22] Q. Zhou, F. He, M.Y. Shen, Direct numerical simulation of a spatially developing compressible plane mixing layer: flow structures and mean flow properties, J. Fluid Mech. 711 (2012) 437–468. [23] T. Rossmann, M.G. Mungal, R.K. Hanson, Evolution and growth of large-scale structures in high compressibility mixing layers, J. Turbul. 3 (2002) N9. [24] P.J. Schmid, Dynamic mode decomposition of numerical and experimental data, J. Fluid Mech. 656 (2010) 5–28. ´ S. Bagheri, et al., Spectral analysis of nonlinear flows, J. [25] C.W. Rowley, I. Mezic, Fluid Mech. 641 (2009) 115–127. [26] A. Wynn, D.S. Pearson, et al., Optimal mode decomposition for unsteady flows, J. Fluid Mech. 733 (2) (2013) 473–503. ´ P.J. Schmid, J.W. Nichols, Sparsity-promoting dynamic mode de[27] M.R. Jovanovic, composition, Phys. Fluids 26 (2) (2014) 024103. [28] M. Grilli, P.J. Schmid, S. Hickel, N.A. Adams, Analysis of unsteady behaviour in shockwave turbulent boundary layer interaction, J. Fluid Mech. 700 (6) (2012) 16–28. [29] T. Sayadi, P.J. Schmid, J.W. Nichols, P. Moin, Reduced-order representation of near-wall structures in the late transitional boundary layer, J. Fluid Mech. 748 (2) (2014) 278–301. [30] J.B. Freund, S.K. Lele, P. Moin, Compressibility effects in a turbulent annular mixing layer. Part 1. Turbulence and growth rate, J. Fluid Mech. 421 (2000) 229–267. [31] Zhang Liu, et al., Direct numerical simulation of subsonic-supersonic mixing layer, Acta Astronaut. 153 (2018) 50–59.
67 68 69 70 71 72 73 74 75 76 77 78 79 80 81 82 83 84 85 86 87
22
88
23
89
24
90
25
91
26
92
27
93
28
94
29
95
30
96
31
97
32
98
33
99
34
100
35
101
36
102
37
103
38
104
39
105
40
106
41
107
42
108
43
109
44
110
45
111
46
112
47
113
48
114
49
115
50
116
51
117
52
118
53
119
54
120
55
121
56
122
57
123
58
124
59
125
60
126
61
127
62
128
63
129
64
130
65
131
66
132