Large Eddy Simulation and Dynamic Mode Decomposition of Turbulent Mixing Layers

: Turbulent mixing layers are canonical ﬂow in nature and engineering, and deserve compre-hensive studies under various conditions using different methods. In this paper, turbulent mixing layers are investigated using large eddy simulation and dynamic mode decomposition. The accuracy of the computations is veriﬁed and validated. Standard dynamic mode decomposition is utilized to ﬂow decomposition, reconstruction and prediction. It was found that the dominant-mode selection criterion based on mode amplitude is more suitable for turbulent mixing layer ﬂow compared with the other three criteria based on singular value, modal energy and integral modal amplitude, respectively. For the mixing layer with random disturbance, the standard dynamic mode decomposition method could accurately reconstruct and predict the region before instability happens, but is not qualiﬁed in the regions after that, which implies that improved dynamic mode decomposition methods need to be utilized or developed for the future dynamic mode decomposition of turbulent mixing layers.


Introduction
Turbulent mixing layers are among the most important fundamental turbulent flows in nature and engineering [1][2][3]. In recent years, with the rapid development of computational technologies, large eddy simulation (LES) and direct numerical simulation (DNS) have become the widely used methods to study turbulent mixing layers. For example, McMullan and Garrett [4] implemented LES to investigate the influence of inlet disturbance on the large-scale spanwise and streamwise structure of turbulent mixing layers; Tan et al. [5] used LES to study the influence of splitter plate cavity on the dynamic development of turbulent mixing layer; Zhang et al. [6] investigated the effect of multiple ring-like vortices on mixing in highly compressible turbulent mixing layer via DNS; Baltzer and Livescu [7] studied the asymmetry of two different density fluids in turbulent mixing layers by DNS; Ren et al. [8] analyzed the interactions of vortex, shock-wave and reaction in droplet laden supersonic turbulent mixing layer with DNS; Chen and Wang [9] found the effects of combustion mode on growth of reacting supersonic turbulent mixing layers through DNS. In these studies, high-precision numerical simulations have revealed many characteristics of turbulent mixing layers, thus continuously advancing the research frontier of turbulent mixing layers.
In the process of high-precision numerical simulation, massive spatial and temporal flow field data will be generated. Extracting valuable information from these flow field data has become a research hotspot in fluid mechanics [10]. Existing methods of extracting flow information include proper orthogonal decomposition (POD) [11], dynamic mode decomposition (DMD) [12] and many others. Among these methods, DMD is a relatively new one. It obtains the main characteristics of unsteady flow by analyzing the flow field data acquired from numerical simulation or experimental measurement. DMD is only based on The turbulent mixing layer adopts the mixing layer with different free flow velocities on upper and lower sides and the same temperature and composition. The freestream temperature T ∞ is 298 K and the freestream composition is air. The upper air freestream velocity U 1 equals to 173 m/s, and the lower air freestream velocity U 2 equals to 86.5 m/s. The initial vorticity thickness is given by The Reynolds number based on the initial vorticity thickness and the freestream velocity difference across the mixing layer Appl. Sci. 2021, 11, 12127 3 of 26 where v is the kinematic viscosity of air at T ∞ . The velocity ratio of the mixing layer, The convective velocity The convective Mach number where c ∞ is the speed of sound of air at T ∞ .

Large Eddy Simulation Method 2.2.1. Computational Domain and the Original Grid
The large eddy simulation is conducted in both 2D and 3D computational domain. The 2D computational domain extends from 0 ∼ 350δ ω (0) in the streamwise x-direction, and from −300δ ω (0) ∼ 300δ ω (0) in the crosswise y-direction. In this computational domain, the streamwise range 0 ∼ 200δ ω (0) is the physical computational domain, the downstream flow range 200δ ω (0) ∼ 350δ ω (0) is an additional buffer to avoid the pollution of the flow field caused by the reflected wave generated by the exit boundary at x = 350δ ω (0) [56]. The 3D computational domain will be discussed in detail in Section 2.2.4.
The original grid points in the x and y directions are 576 and 575, totaling about 300,000 grid points. The grid is stretched from the origin along the y-direction and symmetric with respect to the centerline of the computational domain. The grid stretching ratio in the y-direction is 1.01. In the physical computational domain, the number of grid points in the x-direction is 500, grid uniform distribution and grid point spacing ∆x = 0.4δ ω (0). In the buffer zone, the grid points in x-direction grid is 77, along the x-direction grid stretch, stretch ratio of 1.04. This grid is called original grid.
The computational domain and the original grid as well as the freestream parameters in Section 2.1 are the same as those in the studies by Martha et al. [56] and Uzun [57], making it easy to compare the numerical simulation results.

Boundary Conditions and Inlet Forcing
The inlet boundary adopts the hyperbolic tangent mean velocity profile and the velocity forcing is superposed in the cross direction. The mean streamwise velocity at the inlet is specified as: u(y) = (U 1 + U 2 )/2 + (U 1 − U 2 )/2tanh(2y/δ ω (0)) (6) The velocity profile is a good approximation of the downstream flow of the splitter plate after the wake effect disappears. The mean cross-stream velocity at the entrance v(y) = 0 The cross-stream velocity forcing at the inlet is defined as: v(y) = εαU c exp(− y 2 ∆y 0 2 ) (8) where ε is uniform random value between −1 and 1, α = 0.0045, 0.045, 0.45 (different values represent different perturbations), ∆y 0 = 0.165δ ω (0). The outlet boundary condition is adopted for pressure outlet boundary, and the pressure is constant at 101,325 Pa. The top and bottom boundaries are modeled as slip walls.

Mathematical Model and Numerical Method
Ansys Fluent 19.0 is utilized for the computation. The main governing equations of LES are the unsteady incompressible Navier-Stokes equations after filtering: where σ ij is the viscosity stress, and τ ij is the sub-grid scale (SGS) stress which is defined by The SGS stress is modeled as where ν t is the SGS eddy-viscosity and is determined from the dynamic Smagorinsky model. The second order precision bounded central difference (BCD) scheme is used for spatial discretization, and the pressure-implicit with splitting of operators (PISO) scheme is used to coupling pressure and velocity. The PRESTO scheme is employed for pressure interpolation in the momentum equation, and the 'Green-Gauss cell-based' method is used for gradient calculation [58]. In order to maintain high temporal accuracy, a secondorder-accurate implicit time-stepping scheme is used for time advancement with 20 inner iterations. The time step is 1 × 10 −7 s, and the CFL number is 0.25.
The initial conditions of LES are obtained by using the following methods. First, the RANS method based on the standard k-ε turbulence model is used for steady simulation, and the results of convergence are calculated. Then, the spectrum synthesizer method [59] is used to overlay the velocity fluctuation of RANS results on the velocity field, and finally the synthetic flow field is taken as the initial condition to start LES.
A total of 42,000 time-steps are calculated in LES, which corresponds to 12 flowthrough cycles (FTCs). The FTC is determined by the residence time of the fluid particle moving at convective velocity in the computational domain (including buffer zone). One FTC is about 3500 time-steps. The turbulence statistics involved in the result analysis are carried out for the last eight FTCs.
The whole computation is carried out on the Tianhe-2 supercomputer system of Guangzhou National Supercomputer Center.

Verification and Validation
In order to verify and validate the large eddy simulation method for the turbulent mixing layer, the 2D original case computation, the time step study and the grid resolution investigation are carried out in this paper. Then the results are compared with those of experiments and computations in literatures and those of our 3D computation.
The parameters of the original case are described in Sections 2.2.1 and 2.2.2, where α in the inlet disturbance is taken as 0.0045. Figure 1 shows the comparison of the vorticity thickness growth curve of the turbulent mixing layer with the relevant computational results in reference [56] (Figure 1 includes two kinds of inlet cases without forcing and with forcing, single precision and double precision, sp represents single precision, dp represents double precision). Table 1 shows the comparison of Reynolds stress and turbulent mixing layer growth rate with the experimental and computational results in references [56,57,[59][60][61][62]. It can be seen that these results match well.
forcing, single precision and double precision, sp represents single precision, dp represents double precision). Table 1 shows the comparison of Reynolds stress and turbulent mixing layer growth rate with the experimental and computational results in references [56,57,[59][60][61][62]. It can be seen that these results match well.  In order to study the suitability of the time step, a new time step computational case is carried out. The new time step size of 5 × 10 −8 s is half of the original time step size of 1 × 10 −7 s. The computational settings except for the time step are consistent with the original calculation. The computational results are shown in Table 2. It can be seen that the results of the new time step are very close to the original one, indicating that the time step of 1 × 10 −7 s is suitable. Since the cross-wise size of the computational domain of the original case is too large, it can be optimized. The cross-wise range of the optimized computational domain is −25 ω δ (0)~25 ω δ (0) . Firstly, the corresponding optimized grid is generated according to the grid point distribution method of the original grid. Then, the x-direction is refined by 2 times, y-direction is refined by 2 times, x, y-direction is refined by 1.414 times, x and y-direction is refined by 2 times at the same time to generate four sets of refined grids. The computational settings of the new grid are consistent with those of the original one. The computational results are shown in Table 3. It can be seen that both x and y direction mesh have a great influence on the results. From 1.414-fold to 2-fold, the results converge gradually. Especially, the Reynolds stress in the downstream direction is very close to the previous DNS results.  In order to study the suitability of the time step, a new time step computational case is carried out. The new time step size of 5 × 10 −8 s is half of the original time step size of 1 × 10 −7 s. The computational settings except for the time step are consistent with the original calculation. The computational results are shown in Table 2. It can be seen that the results of the new time step are very close to the original one, indicating that the time step of 1 × 10 −7 s is suitable. Since the cross-wise size of the computational domain of the original case is too large, it can be optimized. The cross-wise range of the optimized computational domain is −25δ ω (0)~25δ ω (0). Firstly, the corresponding optimized grid is generated according to the grid point distribution method of the original grid. Then, the x-direction is refined by 2 times, y-direction is refined by 2 times, x, y-direction is refined by 1.414 times, x and y-direction is refined by 2 times at the same time to generate four sets of refined grids. The computational settings of the new grid are consistent with those of the original one. The computational results are shown in Table 3. It can be seen that both x and y direction mesh have a great influence on the results. From 1.414-fold to 2-fold, the results converge gradually. Especially, the Reynolds stress in the downstream direction is very close to the previous DNS results.
In order to further verify the above 2D computations, 3D studies are also conducted in this paper. The 3D computational domain extends from 0 to 350δ ω (0) in the streamwise x-direction, from −25δ ω (0) to 25δ ω (0) in the crosswise y-direction and from 0 to 10δ ω (0) in the spanwise z-direction. The grid points in the x, y and z directions are 1151, 145 and 51, respectively. The total grid points are of about 8 million (it takes about 50,000 CPU hours to finish the computation of 12 flow-through cycles). In this 3D grid, the minimum grid spacings in the x, y and z directions at the mixing layer centerline are 0.2δ ω (0), 0.165δ ω (0) and 0.2δ ω (0), respectively. The 3D grid also has a buffer zone, which is the same as that in the 2D grid. The spanwise boundaries in the 3D grid are treated as periodic. The rest of the boundary conditions and the numerical schemes are consistent with the 2D cases. To generate three-dimensional disturbance, in addition to adding random disturbance in the y-direction, we also set the velocity fluctuation algorithm of the spectral synthesizer [58] at the inlet. Table 4 demonstrates the 3D and 2D computational results. In this study, the 3D results have been found to be close to the 2D ones, which is consistent with the observations in Martha et al. [56].

Dynamic Mode Decomposition Method
Here, the standard DMD algorithm is briefly introduced. The n matrix snapshots {x 1 , x 2 , x 3 , · · · , x n } obtained by experiment or numerical simulation can be written into a snapshot sequence matrix X and Y. the time interval between any two snapshots is t.
It is assumed that the flow field xi+1 can be expressed by the linear mapping of flow field where A is the system matrix of high dimensional flow field. If the dynamic system itself is nonlinear, then the process is a linear estimation process. According to the assumed linear mapping relation, matrix A can reflect the dynamic characteristics of the system. Due to the high dimension of A, it is necessary to calculate a from the data sequence by order reduction. Therefore For matrix X, a matrix A can be provided to replace the high dimensional matrix A, and the two matrices are similar. In order to find the orthogonal subspace of similarity transformation, the singular value decomposition of X is used to obtain the following: The matrix Σ is a diagonal matrix, and the diagonal elements contain r singular values. In the process of singular value decomposition, only r main singular values can be retained and the remaining small singular values can be truncated, so as to reduce the numerical noise. The unitary matrices U and V obtained by SVD satisfy U H U = I and V H V = I. The calculation process of matrix A can be regarded as the minimization problem of Frobenius norm Then we can approximate A by Since the matrix A is a similar transformation of A, the matrix A contains the main eigenvalues of A. The jth eigenvalue is λ j and the eigenvector is w j . Then the jth DMD mode is defined as The growth rate g j and frequency ω j corresponding to the jth mode are defined as follows: g j = Re lg(λ j )/∆t (20) ω j = Im lg(λ j )/∆t (21) When judging the stability, if the growth rate is positive, the corresponding mode is unstable; if the growth rate is negative, the corresponding mode is stable; if the growth rate is zero, the corresponding mode is periodic; if the eigenvalue falls within the unit circle, it means the stable mode, and vice versa.
The dynamic modes of the flow field can be extracted by the above-mentioned DMD method, and the evolution process of the flow field can be further estimated according to the reduced order matrix A. By the singular value decomposition (15), the high-dimensional system x i can be mapped to the subspace z i The governing equation of the reduced order system is obtained as follows: Let W be a matrix whose column vector is the eigenvector w j , and let N be a diagonal matrix A containing singular values, then the feature decomposition can be expressed as follows: Therefore, the snapshot at any time can be estimated as Each column of Φ is defined as a DMD mode, According to Formula (15), there is The modal amplitude α is defined as where α i is the amplitude of the ith mode, which represents the contribution of the mode to the initial snapshot x 1 . For the standard DMD method, the DMD modes are sorted according to the amplitude. If Equations (26) and (27) are brought into Equation (25), the flow field can be predicted at any time

Characteristics of the Turbulent Mixing Layer
After verifying and confirming the large eddy simulation (LES) method of the turbulent mixing layer, this paper selects the grid which is refined twice in x and y directions to study the basic characteristics and disturbance effects of turbulent mixing layers. The results show that the grid number in x direction is 1151, the mesh number in buffer is 151, the stretch ratio is 1.04, the mesh number in y direction is 145, and the stretch ratio is 1.1. The time step is 1 × 10 −7 s. Double precision is used in the calculation. In the study of basic characteristics, α in the inlet disturbance is obtained as 0.0045 (basic case); in the study of disturbance influence, α in inlet disturbance is obtained as 0, 0.0045, 0.045 and 0.45y. Figure 2a shows the vorticity thickness growth rate of the basic case. It can be seen from the figure that the growth process of the turbulent mixing layer mainly presents two Appl. Sci. 2021, 11, 12127 9 of 26 different stages. In the first stage, the flow is in laminar flow state, and the growth rate is relatively small, and the growth curve is approximately linear. The second stage flow is in turbulent state after transition, and the growth rate of turbulent state is significantly higher than that of the first stage, and the growth curve is also approximately linear.

Basic Characteristics of Turbulent Mixing Layer
Appl. Sci. 2021, 11, x FOR PEER REVIEW 11 of 28 is a big difference between the peak shear Reynolds stress at ω δ = 190 (0) x and ω δ = 160 (0) x , and there is no self-similar state in the mixing layer; however, it can be seen from Figure 6b-d that when there is disturbance, the Reynolds shear stress profile downstream of the flow field is in close proximity, and the greater the disturbance is, the more forward the flow direction of self-similarity appears.    The velocity center line of the turbulent mixing layer y c (x) is obtained according to the following two equations: where U is the average streamwise velocity. Figure 2b shows the velocity center line of the turbulent mixing layer, from which it can be seen that it has a significant bias to the low-speed side. Corresponding to the vorticity thickness growth curve in Figure 2a, there are two main sections in the velocity center line of Figure 2b, and the turning positions of the curves in the two graphs are the same, which indicates that the laminar/turbulent flow pattern also affects the degree of the velocity center line of the turbulent mixing layer to the low-speed side.
The Reynolds shear stress is dimensionless as follows: The normal position is dimensionless as follows: According to these two dimensionless quantities, the self-similarity of turbulent mixing layer can be observed. Figure 2c shows theσ xy − ε curve. According to theσ xy − ε curves of the five flow directions, the turbulent mixing layer is self-similar at about threequarters of the flow direction in the physical computational domain. The Reynolds normal stress curves in x and y directions have similar phenomena, which are not provided here due to the length of paper. Figure 3 shows the vorticity evolution contours in the physical domain of the last FTC (41,300~42,000 steps, with 100 steps interval). It can be seen that the turbulent mixing layer begins to lose stability at about one third of the flow direction in the physical domain, and then generates vortex rolling, falling, pairing and merging.

Influence of Inlet Disturbance on Turbulent Mixing Layer
As mentioned above, in the study of disturbance influence, α in the inlet disturbance is taken as 0, 0.0045, 0.045 and 0.45 respectively. If α = 0.0045 corresponds to the reference disturbance, then α = 0 corresponds to no disturbance, α = 0.045 corresponds to 10 times disturbance, and α = 0.45 to 100 times disturbance. Figure 4 shows the transient instantaneous vorticity contours of the last time of the 12 FTCs under different disturbances. It can be seen that the inlet disturbance significantly advances the flow direction of the mixing layer instability, and the influence on the instability increases with the increase of the inlet disturbance intensity. ω δ = 160 (0) x , and there is no self-similar state in the mixing layer; however, it can be seen from Figure 6b-d that when there is disturbance, the Reynolds shear stress profile downstream of the flow field is in close proximity, and the greater the disturbance is, the more forward the flow direction of self-similarity appears.

Decomposition of Flow Field
In DMD analysis of flow field, 11~12 FTCs (35,000-42,000 steps) are sampled ev 20 steps. A snapshot sample of 351 snapshots will therefore be obtained.

DMD Analysis of Basic Case
In this section, DMD analysis is carried out for x-direction velocity U, y-direction locity V and vorticity ω.
Firstly, the x-direction velocity U is analyzed and as shown in Figure 7, the singu value arrangement, growth rate, frequency diagram, eigenvalue distribution diagram a part of the DMD modal contours (showing only the physical computational dom (0~0.0254 m)) is obtained. Figure 7a shows the order of singular values. The DMD modes are arranged in order of singular values. As can be seen from Figure 7a, the first-order singular value more than one order of magnitude higher than other singular values. In Table 5, the s ond column gives the singular values of the first 10 modes. The first-order singular va is utterly dominant, implying that after flow field decomposition, the first-order mode the most important mode. From Figure 7a, it can be seen that the first mode is consist with the average flow field.  Figure 5b shows the variation curves of the velocity center line of the turbulent mixing layer under different disturbances. It can be seen that the greater the inlet disturbance, the greater the degree of the velocity center line deflecting to the low speed side. Under various disturbances, the turning position of velocity center line curve is always the same as that of vorticity thickness growth curve. Figure 6 shows the Reynolds shear stress profile of turbulent mixing layer under different disturbances. It can be seen from Figure 6a that when there is no disturbance, there is a big difference between the peak shear Reynolds stress at x = 190δ ω (0) and x = 160δ ω (0), and there is no self-similar state in the mixing layer; however, it can be seen from Figure 6b-d that when there is disturbance, the Reynolds shear stress profile downstream of the flow field is in close proximity, and the greater the disturbance is, the more forward the flow direction of self-similarity appears.

Decomposition of Flow Field
In DMD analysis of flow field, 11~12 FTCs (35,000-42,000 steps) are sampled every 20 steps. A snapshot sample of 351 snapshots will therefore be obtained.

DMD Analysis of Basic Case
In this section, DMD analysis is carried out for x-direction velocity U, y-direction velocity V and vorticity ω.
Firstly, the x-direction velocity U is analyzed and as shown in Figure 7, the singular

Decomposition of Flow Field
In DMD analysis of flow field, 11~12 FTCs (35,000-42,000 steps) are sampled every 20 steps. A snapshot sample of 351 snapshots will therefore be obtained.

DMD Analysis of Basic Case
In this section, DMD analysis is carried out for x-direction velocity U, y-direction velocity V and vorticity ω.
Firstly, the x-direction velocity U is analyzed and as shown in Figure 7, the singular value arrangement, growth rate, frequency diagram, eigenvalue distribution diagram and part of the DMD modal contours (showing only the physical computational domain (0~0.0254 m)) is obtained. Figure 7a shows the order of singular values. The DMD modes are arranged in the order of singular values. As can be seen from Figure 7a, the first-order singular value is more than one order of magnitude higher than other singular values. In Table 5, the second column gives the singular values of the first 10 modes. The first-order singular value is utterly dominant, implying that after flow field decomposition, the first-order mode is the most important mode. From Figure 7a, it can be seen that the first mode is consistent with the average flow field. Figure 7b shows the relationship between the modal growth rate and the frequency. The modal growth rate is symmetrical with respect to the frequency. The growth rate of the mode is equal to 0, which indicates that the mode does not shift. The mode growth rate is less than 0, which indicates that the mode is decaying and stable. When the mode growth rate is greater than 0, the mode is unstable. Figure 7b shows that most of the modal growth rates are below 0, which indicates that the whole system is stable. Figure 7c shows the distribution of eigenvalues, where the horizontal axis is the real part of the modal eigenvalues and the vertical axis corresponds to the imaginary part. Most of the eigenvalues fall on and within the unit circle, and a few points fall outside the unit circle. Among them, mode 1 is located on the unit circle and the growth rate is 0, so mode 1 is static mode; the mode corresponding to the eigenvalue falling in the unit circle is stable mode; the mode corresponding to the eigenvalue falling outside the unit circle is unstable mode. Figure 7d shows the relationship between modal amplitude and frequency. Modal amplitude is symmetrical with respect to frequency. The amplitude of the first mode is the largest, which indicates that the contribution of the first mode is the largest. Figure 7e shows the contours of the first 10 modes, and Table 5 shows the singular values, growth rates and frequencies of the first 10 modes. It can be seen that except for mode 1, the other modes are paired. The singularities of the paired modes are close, the growth rates are equal, and the frequencies are opposite to each other. This is because the DMD modes are conjugate modes, which are essentially the same mode. The results show that the first 10 modal structures are strip shape in the flow direction, and the strip structure mainly exists in the middle and rear of the flow field, which is close to the vortex position of the mixing layer flow field, indicating that the structure is closely related to the formation of vortex. Figure 8 shows the DMD analysis results of velocity V in the y-direction. As can be seen from Figure 8a, there is no dominant singular value of velocity V. The singular values, growth rates and frequencies of the first 10 DMD modes are given in Table 6. From the growth rates and frequencies, it can be seen that the first 10 DMD modes are five pairs of conjugate modes, that is, the growth rates are equal and the frequencies are opposite to each other. The growth rates of the first six modes are all greater than 0, and they are all unstable modes. Figure 8b shows the relationship between the mode growth rate and frequency. The growth rate of the first mode is 1756, which is greater than 0, and the frequency is not 0. It indicates that it is unstable and not the primary mode of flow field. None of the modes meets the frequency and growth rate of 0 at the same time. Figure 8c shows the eigenvalues basically fall near the unit circle, but the first mode does not fall at (1,0) and none of the modes fall at (1,0). Figure 8d shows that the amplitude of the first mode is much smaller than that of the maximum mode.   It can be seen from Figure 8e that the first 10 modes are identical, indicating that the corresponding modes are conjugate modes, and the flow field structure presents an alternating circular structure.
The results of DMD analysis of vorticity ω are shown in Figure 9. The DMD result of vorticity ω is similar to that of velocity U in x-direction. It is visible from Figure 9a that the dominant singular value of vorticity ω exists, and the first singular value is greater than the other singular values by more than one order of magnitude. It can be seen from Table  7 that except for the first mode, the others are conjugate modes. The growth rates of the first 10 modes are less than 0, indicating that they are stable modes. In Figure 9b, the first mode's growth rate and frequency are both equal to zero. It shows that it is stable and the primary mode of flow field. Figure 9c shows that the eigenvalue basically falls near the unit circle, and the first mode falls at point (1, 0). Figure 9d shows the amplitude of the first mode is the largest. Figure 9e shows the first 10 modes, which is similar to the first 10 modes of velocity U. A slender and large-scale structure is shown by low-frequency modes 2-9.

DMD Analysis of Inlet Disturbance Characteristics
DMD analysis is carried out for different inlet disturbance conditions and the influence of disturbance amplitude is studied. The previous basic case is random disturbance condition α = 0.0045 (1 forcing). Here, two groups of random disturbance condition α = 0.045 (10 forcing), 0.45 (100 forcing) and one group of no forcing condition α = 0 are added. DMD analysis was carried out on the x-direction velocity U of the three groups of working conditions, and 351 snapshots were also sampled.

DMD Analysis of Inlet Disturbance Characteristics
DMD analysis is carried out for different inlet disturbance conditions and the influence of disturbance amplitude is studied. The previous basic case is random disturbance condition α = 0.0045 (1 forcing). Here, two groups of random disturbance condition α =   Figure 10a shows the order 2~350 singular value arrangement under four disturbance conditions. The first-order singular values of no forcing, forcing, 10 forcing and 100 forcing are 1,040,287.5, 1,045,597.5, 1,046,331.7 and 1,046,807.3, respectively. It can be seen that the larger the disturbance is, the larger the first-order singular value is. However, the difference between the three perturbed singular values is smaller than that between the perturbed and undisturbed singular values. As can be seen from Figure 10a, the first-order singular value is more than one order of magnitude larger than other singular values. The singular value of DMD with forcing is obviously larger than that without forcing. There is no significant difference in the singular values of the three disturbed conditions. It demonstrates that the composition of the flow field of the three kinds of conditions of disturbance is similar.
l. Sci. 2021, 11, x FOR PEER REVIEW 19 of the larger the disturbance is, the larger the first-order singular value is. However, the d ference between the three perturbed singular values is smaller than that between the p turbed and undisturbed singular values. As can be seen from Figure 10a, the first-ord singular value is more than one order of magnitude larger than other singular values. T singular value of DMD with forcing is obviously larger than that without forcing. Th is no significant difference in the singular values of the three disturbed conditions demonstrates that the composition of the flow field of the three kinds of conditions disturbance is similar. Figure 10b shows the relationship between the growth rate of the four disturban and the frequency. Since the growth rate is symmetrical about the frequency, Figure 1 only shows the right half for easy observation. It can be seen that the growth rates of four disturbances are mainly concentrated in the region of −5000~0, indicating that m of the modes are stable, and there is a point at the origin under each condition, and corresponding dynamic mode at this point corresponds to the mean flow structure. general, the larger the disturbance is, the farther away the points are from the origin shows that the nonlinear intensity of the disturbed flow field is directly related to magnitude of the disturbance.

Reduced Order Reconstruction and Prediction of Basic Case
In this section, the flow field under basic case is reduced, and the first 150 orders truncated (the sum of the first 150 singular values accounts for more than 90% of the to singular values). The reduced order reconstruction and prediction effect of DMD on turbulent mixing layer flow field with strong nonlinearity are studied. The x-direct velocity U, y-direction velocity V and vorticity ω are reconstructed and predicted, resp tively. Figure 11a shows the comparison between the reconstructed contours and the inst taneous contours. The physical computational domain (0~0.0254 m) is compared. Acco ing to the vorticity contours, the flow field is divided into three regions: the region bef  Figure 10b shows the relationship between the growth rate of the four disturbances and the frequency. Since the growth rate is symmetrical about the frequency, Figure 10b only shows the right half for easy observation. It can be seen that the growth rates of the four disturbances are mainly concentrated in the region of −5000~0, indicating that most of the modes are stable, and there is a point at the origin under each condition, and the corresponding dynamic mode at this point corresponds to the mean flow structure. In general, the larger the disturbance is, the farther away the points are from the origin. It shows that the nonlinear intensity of the disturbed flow field is directly related to the magnitude of the disturbance.

Reduced Order Reconstruction and Prediction of Basic Case
In this section, the flow field under basic case is reduced, and the first 150 orders are truncated (the sum of the first 150 singular values accounts for more than 90% of the total singular values). The reduced order reconstruction and prediction effect of DMD on the turbulent mixing layer flow field with strong nonlinearity are studied. The x-direction velocity U, y-direction velocity V and vorticity ω are reconstructed and predicted, respectively. Figure 11a shows the comparison between the reconstructed contours and the instantaneous contours. The physical computational domain (0~0.0254 m) is compared. According to the vorticity contours, the flow field is divided into three regions: the region before instability, the rolling up region and the vortex interaction region. In the upstream of the flow field, there is no obvious fluid exchange and vortex structure. The vorticity concentrates on the shear line, and its width increases gradually along the flow direction. We call this region the pre-instability region. In the middle of the flow field (i.e., downstream of the region before instability), there is obvious material exchange between the upper and lower flow fields. An obvious single vortex (a series of vortex structures) can be observed from the vorticity contours. The distance between these vortex structures is approximately proportional to their own scale, and gradually increases along the flow direction. We call this region the roll up region. In the downstream of the flow field (that is, the downstream of the roll up region), the interaction between the vortices of the flow field occurs. We call this region the vortex interaction region. proportional to their own scale, and gradually increases along the flow direction. We call this region the roll up region. In the downstream of the flow field (that is, the downstream of the roll up region), the interaction between the vortices of the flow field occurs. We call this region the vortex interaction region. It can be seen that in the basic case, whether it is velocity U, V or vorticity ω, the reconstruction effect of the flow field is poor, and only the region before instability (about 0~0.006 m region) can be reproduced. For the roll up region (about 0.006-0.015 m region) and vortex interaction region (about 0.015~0.0254 m region), the reconstruction effect of DMD is not good because of the strong nonlinearity of the region.
Six points (0.005,0), (0.008,0), (0.011,0), (0.014,0), (0.017,0) and (0.02,0) on the middle line of the physical computational domain are selected to further analyze the reconstruction and prediction of DMD. Figure 11b is the reduced order reconstruction and prediction relative error diagrams of perturbed velocity U. Among them, 0.0035 s~0.0042 s is the time range within the sample, and 0.0042 s-0.0049 s is the predicted time range.
The relative error δ is defined as where a is the reconstructed and predicted value of DMD and A is the actual value.
As can be seen more clearly from Figure 11b, the error of reconstruction segment is obviously less than that of prediction segment, most of the errors of reconstruction segment are less than 10%, a considerable part of the errors of prediction segment are more than 20%, and individual errors are more than 40%. The effect of velocity V and vorticity ω is similar, as such there is no additional illustration.  Figure 11b is the reduced order reconstruction and prediction relative error diagrams of perturbed velocity U. Among them, 0.0035 s~0.0042 s is the time range within the sample, and 0.0042 s-0.0049 s is the predicted time range.
The relative error δ is defined as where a is the reconstructed and predicted value of DMD and A is the actual value.
As can be seen more clearly from Figure 11b, the error of reconstruction segment is obviously less than that of prediction segment, most of the errors of reconstruction segment are less than 10%, a considerable part of the errors of prediction segment are more than 20%, and individual errors are more than 40%. The effect of velocity V and vorticity ω is similar, as such there is no additional illustration.

Influence of Dominant Mode Selection Method on Reduced Order Reconstruction and Prediction
After the flow field is decomposed, the selection of the dominant mode is the key step of flow field reduction reconstruction and prediction. At present, there is no specific method to select the dominant mode of the flow field after DMD decomposition. For different flow characteristics, the most suitable mode selection method may be different. In this section, four modal selection methods are compared: (1) according to the size of singular value; (2) according to the size of modal amplitude; (3) according to the size of modal energy [63]; (4) according to the size of the integral modal amplitude of sampling time [64], to investigate the mode selection method of the turbulent mixing layer suitable for this study.
Here, four methods are briefly introduced. The first and second methods are easier to understand. The expression of the third method is as follows: where ϕ i is the energy contained in the ith mode and ϕ i is the flow velocity of each node in the ith mode.
In the fourth method, the integral modal amplitude is defined as [64]: DMD analysis is carried out for the stochastically perturbed velocity U. Figure 12a shows the relationship between modal energy and frequency, from which it can be seen that high modal energy (greater than 0.7) is mainly concentrated near low frequency and high frequency, and the modal energy of middle frequency is almost about 0.7. The first three modes with the largest energy are shown in Figure 12a, which is named mode 1, mode A and mode B. The modal energy of mode 1 is the largest, and the frequency is 0 Hz. Mode 1 is the first-order mode arranged according to the eigenvalue, which also shows that the first-order mode is the most important mode. Mode A is the second-highest energy mode with a frequency of 1,567,540 Hz. Mode B is the third-highest energy mode with a frequency of 1,570,796 Hz. Modes A and B are high-frequency modes in the flow field. Figure 12b shows the relationship between the integral mode amplitude and frequency, and marks the three modes with the largest amplitude of integral mode which are also mode 1, mode A and mode B. Different from modal energy, modal B is larger than modal A. This shows that different methods extract the dominant mode differently. Compared with Figure 7d, we can see that the morphology has changed significantly. The integral amplitude of the first-order mode is still the largest, which shows once again that the first-order mode is the most important. Figure 13 shows the contours of mode A and mode B. compared with modes 1~10 in Figure 7, it can be seen that the high-frequency modes present row upon row of small-scale structures. The low-frequency modes present large-scale structures. modal A. This shows that different methods extract the dominant mode differently. Compared with Figure 7d, we can see that the morphology has changed significantly. The integral amplitude of the first-order mode is still the largest, which shows once again that the first-order mode is the most important. Figure 13 shows the contours of mode A and mode B. compared with modes 1~10 in Figure 7, it can be seen that the high-frequency modes present row upon row of smallscale structures. The low-frequency modes present large-scale structures.  In order to ensure a certain accuracy and reduce the ord methods are used to select 150 modes to reconstruct and predict The comparison results are shown in Figures 14 and 15. It can be the reconstructed contours by methods 2 and 3 are close to th and can reconstruct a smaller scale structure. The reconstructed and 4 are not ideal, only some large structures can be reconstruc ence is large, which is very different from the instantaneous flo the comparison between the predicted contours and the instanta seen that the predicted results of all methods have some deviatio of method 1 and method 3 are almost the same, both of which ov The result of method 4 is still not optimal. In contrast, the predi are closer to the transient flow field, which is better than the oth fore, method 2 is selected as the modal selection method in this In order to ensure a certain accuracy and reduce the order of the flow field, four methods are used to select 150 modes to reconstruct and predict the velocity U flow field. The comparison results are shown in Figures 14 and 15. It can be seen from Figure 14 that the reconstructed contours by methods 2 and 3 are close to the instantaneous contours and can reconstruct a smaller scale structure. The reconstructed contours by methods 1 and 4 are not ideal, only some large structures can be reconstructed, and the shape difference is large, which is very different from the instantaneous flow field. Figure 15 shows the comparison between the predicted contours and the instantaneous contours. It can be seen that the predicted results of all methods have some deviation. The prediction results of method 1 and method 3 are almost the same, both of which overestimate the flow field. The result of method 4 is still not optimal. In contrast, the prediction results of method 2 are closer to the transient flow field, which is better than the other three methods. Therefore, method 2 is selected as the modal selection method in this study. ence is large, which is very different from the instantaneous flow field. Figure 15 shows the comparison between the predicted contours and the instantaneous contours. It can be seen that the predicted results of all methods have some deviation. The prediction results of method 1 and method 3 are almost the same, both of which overestimate the flow field. The result of method 4 is still not optimal. In contrast, the prediction results of method 2 are closer to the transient flow field, which is better than the other three methods. Therefore, method 2 is selected as the modal selection method in this study.  In this section, method 2 is used to select 150 modes for reconstruction and prediction to analyze the influence of disturbance existence (the difference between the undisturbed condition and the disturbed condition) and disturbance form (the difference between random disturbance and periodic disturbance).

Influence of Inlet Disturbance on Reduced Order Reconstruction and Prediction
In this section, method 2 is used to select 150 modes for reconstruction and prediction to analyze the influence of disturbance existence (the difference between the undisturbed condition and the disturbed condition) and disturbance form (the difference between random disturbance and periodic disturbance). Figure 16a shows the comparison between the undisturbed reconstructed contours (right side) and the instantaneous contours (left side). Figure 17a shows the comparison between the reconstructed contours (right side) and the instantaneous contours (left side). From top to bottom are U, V, and ω. Under the undisturbed condition, the flow field in the computational domain belongs to the region before instability, and the reconstructed contours are almost the same as the instantaneous contours. Compared with Figure 17a, it can be seen that the reconstruction effect of undisturbed DMD is better than that of perturbed DMD.
Appl. Sci. 2021, 11, x FOR PEER REVIEW 24 of 28 Figure 18a shows the comparison between the reconstructed contours of 150 modes selected by method 2 and the instantaneous contours. It is obvious that the reconstruction effect is ideal, and the reconstructed contours are almost the same as the instantaneous contours. The shape of the whole flow field, the size and the position of the vortex are almost the same. From the comparison of the instantaneous contours in Figures 17a and  18a, it can be seen that the flow field under the condition of random disturbance develops faster than that under the condition of periodic disturbance. The results show that the flow field under the random disturbance condition consists of the pre-instability region, the rolling up region, the vortex interaction region, and the rolling up the position of the vortex is about 0.007 m in the x-direction. However, there is no obvious vortex interaction zone in the flow field under the condition of periodic disturbance, and the vortex curling position is about 0.017 m in the x-direction.    Figure 18a shows the comparison between the reconstructed contours of 150 modes selected by method 2 and the instantaneous contours. It is obvious that the reconstruction effect is ideal, and the reconstructed contours are almost the same as the instantaneous contours. The shape of the whole flow field, the size and the position of the vortex are almost the same. From the comparison of the instantaneous contours in Figures 17a and  18a, it can be seen that the flow field under the condition of random disturbance develops faster than that under the condition of periodic disturbance. The results show that the flow field under the random disturbance condition consists of the pre-instability region, the rolling up region, the vortex interaction region, and the rolling up the position of the vortex is about 0.007 m in the x-direction. However, there is no obvious vortex interaction zone in the flow field under the condition of periodic disturbance, and the vortex curling position is about 0.017 m in the x-direction.  As before, the same six points for reconstruction and prediction were selected. The relative error is shown in Figures 16b and 17b. The fluctuation of each undisturbed point value with time is much smaller. In particular, the actual velocity of the points (0.005,0) and (0.008,0) near the entrance is almost unchanged. The relative error is within 1%. In general, the reconstruction and prediction of undisturbed DMD are better than that of perturbed DMD.
The influence of disturbance form is studied below. The disturbance conditions in the previous sections are all random disturbances. We added a group of periodic disturbance conditions for comparative study. The computational settings are consistent with that in Section 2. The sampling method is also consistent with random disturbance, and 351 snapshots are obtained. The periodic disturbance formula is as follows: v(y) = αU c exp(−y 2 /∆y 0 2 )sin(6 .28 × 14, 900 × t ) (37) where α, ∆y 0 consistent with Section 2.2.2, t is the time variable. Figure 18a shows the comparison between the reconstructed contours of 150 modes selected by method 2 and the instantaneous contours. It is obvious that the reconstruction effect is ideal, and the reconstructed contours are almost the same as the instantaneous contours. The shape of the whole flow field, the size and the position of the vortex are almost the same. From the comparison of the instantaneous contours in Figures 17a and 18a, it can be seen that the flow field under the condition of random disturbance develops faster than that under the condition of periodic disturbance. The results show that the flow field under the random disturbance condition consists of the pre-instability region, the rolling up region, the vortex interaction region, and the rolling up the position of the vortex is about 0.007 m in the x-direction. However, there is no obvious vortex interaction zone in the flow field under the condition of periodic disturbance, and the vortex curling position is about 0.017 m in the x-direction. Similarly, the six selected points are reconstructed and predicted, and the relative error is shown in Figure 18b. Except for the point (0.02, 0), the error of other points is less than 5%. Compared with the results of random disturbance in Figure 17, it can be seen that DMD has poor reconstruction and prediction effect for highly nonlinear and aperiodic conditions, but good reconstruction and prediction effect for periodic disturbance conditions.

Conclusions
In the present paper, large eddy simulation (LES) and dynamic mode decomposition (DMD) were used to study turbulent mixing layers. The results allow for the following conclusions to be drawn: (1) The LES method used in this paper can accurately calculate the vorticity thickness growth rate of the turbulent mixing layer, and the peak values of each component of Reynolds stress at several streamwise locations are close to the results of DNS in literature, indicating that the present LES provides accurate big data for DMD.
(2) For the velocity in x-direction, the first-order mode is the dominant mode of the flow field, which represents the average flow field; the low-frequency mode presents a strip-shaped large-scale structure, and the high-frequency mode presents a small-scale structure, which indicates that the small-scale vortex formation in the flow field is mainly related to the high-frequency mode. For the velocity in y-direction, the first-order mode Similarly, the six selected points are reconstructed and predicted, and the relative error is shown in Figure 18b. Except for the point (0.02, 0), the error of other points is less than 5%. Compared with the results of random disturbance in Figure 17, it can be seen that DMD has poor reconstruction and prediction effect for highly nonlinear and aperiodic conditions, but good reconstruction and prediction effect for periodic disturbance conditions.

Conclusions
In the present paper, large eddy simulation (LES) and dynamic mode decomposition (DMD) were used to study turbulent mixing layers. The results allow for the following conclusions to be drawn: (1) The LES method used in this paper can accurately calculate the vorticity thickness growth rate of the turbulent mixing layer, and the peak values of each component of Reynolds stress at several streamwise locations are close to the results of DNS in literature, indicating that the present LES provides accurate big data for DMD.
(2) For the velocity in x-direction, the first-order mode is the dominant mode of the flow field, which represents the average flow field; the low-frequency mode presents a strip-shaped large-scale structure, and the high-frequency mode presents a small-scale structure, which indicates that the small-scale vortex formation in the flow field is mainly related to the high-frequency mode. For the velocity in y-direction, the first-order mode is unstable, which cannot represent the average flow field; there is no mode in which frequency and growth rate are both zero, which indicates that the flow field of velocity in y-direction is more complex than that of velocity in x-direction.
(3) The dominant-mode selection criterion based on mode amplitude is more suitable for turbulent mixing layer flow compared with the other three criteria based on singular value, modal energy and integral modal amplitude.
(4) The standard DMD method reconstructs and predicts the periodically-perturbed mixing layer well. In contrast, for the mixing layer with random disturbance, the standard DMD method could only accurately reconstruct and predict the region before instability happens but is not qualified in the regions after that as nonlinearity is stronger, which implies that improved dynamic mode decomposition methods need to be utilized or developed for the future dynamic mode decomposition of turbulent mixing layers.
Author Contributions: Investigation, computation, visualization, formal analysis, and writingoriginal draft, Y.C.; conceptualization, methodology, investigation, formal analysis, writing-original draft, writing-review, and supervision, Q.C. All authors have read and agreed to the published version of the manuscript.