A transition prediction method for flow over airfoils based on high-order dynamic mode decomposition

Mengmeng WU, Zhonghu HAN,*, Hn NIE, Wenping SONG,Soledd Le CLAINCHE, Esten FERRER

a School of Aeronautics, Northwestern Polytechnical University, Xi’an 710072, China

b ETSIAE-UPM, School of Aeronautics, Universidad Polite′cnica de Madrid, Pza Cardenal Cisneros 3, Spain

KEYWORDS

Abstract This article presents a novel approach for predicting transition locations over airfoils,which are used to activate turbulence model in a Reynolds-averaged Navier-Stokes flow solver.This approach combines Dynamic Mode Decomposition (DMD) with eN criterion. The core idea is to use a spatial DMD analysis to extract the modes of unstable perturbations from a steady flowfield and substitute the local Linear Stability Theory (LST) analysis to quantify the spatial growth of Tollmien-Schlichting (TS) waves. Transition is assumed to take place at the stream-wise location where the most amplified mode’s N-factor reaches a prescribed threshold and a turbulence model is activated thereafter. To improve robustness, the high-order version of DMD technique (known as HODMD) is employed. A theoretical derivation is conducted to interpret how a spatial highorder DMD analysis can extract the growth rate of the unsteady perturbations. The new method is validated by transition predictions of flows over a low-speed Natural-Laminar-Flow (NLF) airfoil NLF0416 at various angles of attack and a transonic NLF airfoil NPU-LSC-72613. The transition locations predicted by our HODMD/eN method agree well with experimental data and compare favorably to those obtained by some existing methods . It is shown that the proposed method is able to predict transition locations for flows over different types of airfoils and offers the potential for application to 3D wings as well as more complex configurations.?2019 Chinese Society of Aeronautics and Astronautics.Production and hosting by Elsevier Ltd.This is an open access article under the CC BY-NC-ND license(http://creativecommons.org/licenses/by-nc-nd/4.0/).

1. Introduction

Transition from laminar to turbulent flow over airfoils presents great interest for aircraft manufacturers.The loss of performance and drag penalties associated with the prediction of early transition, has pushed academics and engineers to conceive transition prediction methods that simulate the flow transition at an affordable computational cost.

Industrial simulations of flow over airfoils and wings do not typically consider transition explicitly by means of Direct Numerical Simulation (DNS) or Large Eddy Simulation(LES)computations.1-3Instead,boundary layers are split into laminar and turbulent regions with a prescribed transition location and a Reynolds-Averaged Navier-Stokes (RANS)simulation is adopted thereafter. After the prescribed transition point, the flow is considered fully turbulent and a turbulence model is activated to account for turbulent effects, e.g.increased skin friction. Before the transition point, the flow is laminar and no turbulence modelling is necessary.

To find transition point (or transition line in 3D computations), industrial design relies on well-established techniques such as Linear-Stability Theory (LST) analysis.4The LST relies on the parallel flow (slowly varying flow in the streamwise direction) assumption to derive the Orr-Sommerfeld and Squire equations that govern the spatial growth of Tollmien-Schlichting (TS) waves. When the amplitudes of 2D TS waves grow to a prescribed value eNcrat a point, where Ncrdenotes the critical N-factor, flow transition is considered to take place. After this point, a turbulence model should be switched on. The value of Ncr, having been calibrated empirically during the years,5varies in every application and depends on many factors such as free-stream turbulence or surface roughness. Typically, the commonly accepted value is in the range from Ncr=6 to 12. This technique that combines LST analysis and eNcriterion is denoted by LST/eN. The LST/eNmethod has been proved successful in predicting transition locations for 2D and 3D configurations and is now a well-established and cost-effective method for design of airfoil/wing/wing-body configurations.6-16A detailed historical review for LST/eNmethod has been compiled by Van Ingen.5Despite its success,LST/eNmethod has the limitations inherited from the hypotheses underlying its derivation.Namely, the method relies on the parallel flow assumptions(i.e. relies on linear-stability analysis) and was conceived for TS-driven transition: exponential spatial growth of TS waves and subsequent nonlinear interactions that lead to turbulent flow. In theory, the original LST/eNmethod cannot be used to predict transition upon non-parallel flows,or where physical mechanisms,such as the cross-flow or transient-growth effects,govern the transition. Nonetheless, modifications have been proposed to account for these complex physics.5,17

As a summary of literature survey, we found that despite the success of the existing transition prediction methods,there is still a strong need for the community to develop improved or new methods that could overcome their drawbacks towards more effective transition prediction in industrial application context, which motivates the research of this article.

The objective of this article is to propose an alternative and novel approach, which is called ‘‘HODMD/eN” method. This method is similar to an LST/eNmethod in which a critical N-factor Ncris required to determine the transition point,but it does not rely on local linear-stability analysis. The core idea of the new method is to use a spatial DMD analysis to extract the modes of perturbations from a steady flowfield and substitute the local LST analysis to quantify the spatial growth of TS waves. The implementation of the proposed method may be summarized in three steps. First, a transition-point-fixed computation of flow over an airfoil is performed to obtain a steady laminar flowfield.Second,a spatial High-Order Dynamic Mode Decomposition (HODMD)28is performed on the laminar region of the flowfield to obtain the eigenvalue and eigenmodes of perturbations that indicate the wavenumber and growth information of TS waves. Third,once the spatial growths are extracted through the DMD analysis, an eNcriterion is employed to determine the transition location.

This article continues in Section 2 for presenting the methodology paying particular attention to how to predict the transition with an HODMD/eNmethod and how it is coupled with a RANS flow solver. Besides, the theory and algorithm of DMD and HODMD are briefly introduced. A theoretical derivation is also performed to explain why the DMD analysis of a steady laminar flowfield can provide the growth rate of unsteady perturbations. In Section 3, the capability of HODMD for extracting complex modes is validated through an analytical test case. Then the proposed method is validated and demonstrated by the examples of transition predictions of flows over a low-speed NLF airfoil and a transonic NLF airfoil. The last section gives some conclusions from the current study and presents the outlook for future study.

2. Methodology

2.1. Framework of flow simulation with transition prediction

In this article, the flow simulation considering laminar-toturbulent transition is performed, coupling a RANS solver with a transition prediction technique, as shown in Fig.1.First, RANS simulation with prescribed transition locations is performed to obtain a steady flowfield, from which the modes of unstable perturbations are extracted. Theoretically,a full laminar flow should be simulated. However, this may lead to the difficulty associated with the convergence to a steady state. The researchers in this area have reached a conclusion that prescribed transition with sufficient laminar flow region is necessary. Second, the transition prediction method returns a transition location to the RANS solver. Third, a new flowfield is obtained with the new transition location.The above procedure is repeated until the whole process converges.

Fig.1 Diagram of flow simulation coupling transition prediction technique with a RANS solver.

Fig.2 Diagram of stability analysis and transition prediction by DMD/eN or HODMD/eN method.

In the following sections, we will briefly introduce the flow solver, transition prediction methods (LST and the proposed method), the theory of DMD and HODMD, as well as the interpretation of the proposed method.

2.2. Flow solver

An in-house CFD solver, called ‘‘PMNS2D”,35,36is used in this work. The RANS equations are discretized spatially by the Jameson-Schmidt-Turkel scheme. The LU-SGS scheme is implemented for time integration and a full-approximationstorage multi-grid technique is used to accelerate the convergence to a steady state. The turbulent viscosity is calculated through the Spalart-Allmaras one-equation turbulence model.To improve the efficiency and robustness of a compressible CFD code for low-speed flows or flows with low-speed flow regions, a preconditioning method is adopted within the framework of a Full Approximation Storage (FAS) multigrid method.35

2.3. Transition prediction methods

2.3.1. LST/eNmethod

Let us first review the basic features of an LST/eNmethod, to provide a context for our DMD/eNand HODMD/eNmethods,which will be detailed in the next section.

The LST/eNmethod4-17is based on a LST analysis, and calculates the growth of perturbation amplitudes starting from the stagnation point of boundary layer. In this method, the Orr-sommerfeld and Squire (OS) equations are solved to obtain a neutral curve as well as the evolution of unstable waves of various frequencies. The N-factor relates to the spatial growth rate of the most unstable perturbation as where αidenotes the spatial growth rate of each ith sinusoidal wave, ω is the temporal frequency, and s is the airfoil arclength position.

Fig.3 Schematic diagram of segments for stability analysis and N-factor calculations.

Fig.4 Schematic diagram of snapshot grid for stability analysis via DMD.

Typically,the OS system is solved in time and hence dependent on spatial growth with temporal frequency ω. However,for a TS-driven transition, only the spatial growth is of interest.Transition is assumed to take place where the most amplified mode’s N-factor at a given stream-wise location reaches a prescribed threshold (i.e. the critical N-factor or Ncr). The application of an eNcriterion involves three steps: first to calculate the laminar velocity profiles (the base flow) at different stream-wise locations for a given geometry of interest, second to calculate the neutral curve as well as the local growth rate of the most unstable waves and third to calculate the transition onset location based on the local growth rates integrated along the streamline direction, following Eq. (1).

2.3.2. DMD/eNmethod and HODMD/eNmethod

In our method, DMD or HODMD analysis (see next section for details)is used to extract the modes of small perturbations from the steady flowfield inside the boundary layer, and then the spatial growth and wavenumber of these unstable perturbations are obtained.In this method,DMD or HODMD actually serves as an alternative stability analysis to an LST. The proposed method can be summarized in the following steps(see Fig.2).

First, numerical simulation of flow over an airfoil is performed until it converges to a steady state, typically when the relative averaged-density residual of the compressible Navier-Stokes system is below 10-6.

Second, the boundary layer of the airfoil is divided into L segments along the arc length from the stagnation point to the trailing edge point: [0, 1/L], [0, 2/L], [0, 3/L], ..., [0,(L-1)/L],[0,1]as shown in Fig.3.A‘‘snapshot grid”is generated to extract spatial snapshots for stability analysis. The snapshot grid lines start from equally spaced points over the airfoil surface, following the normal direction of the airfoil and ending at the boundary layer edge, as sketched in Fig.4.It should be mentioned that the CFD grid for flow simulation and the snapshot grid for DMD analysis are generated separately.The relationship between them is that the flow solution is interpolated from the CFD grid to the snapshot grid. As shown in Fig.5, the flow property φ on a ‘‘snapshot grid”point is obtained by linear interpolations of flow solution φ1-φ4along two orthogonal directions.

To ensure that the accuracy of the interpolated data is accurate enough, a convergence study of the ‘‘snapshot grid” size should be conducted using a fine CFD mesh. As can be seen in Fig.6, the interpolated velocity contours on snapshot grids agree perfectly with the original velocity contours on CFD mesh,which proves that our interpolation method is effective.

Third,the stability analysis based on DMD or HODMD is performed for each of the segments,written as[0,1/L],[0,2/L],[0, 3/L], ..., to obtain the local growth/decay rates of spatial perturbations (positive value for growth and negative value for decay). We use all ‘‘snapshot sections” contained in one region to compute the growth rate. The most unstable eigenvalue in the eigenvalue map is taken as an effective eigenvalue.

Fourth, we compute the N-factor associated with the spatial growth in one segment, previously extracted by the DMD/HODMD technique, using the definition:

Fig.5 Schematic diagram of interpolating flow properties from CFD mesh to snapshot grids.

Fig.6 Comparison of original and interpolated velocity contours in a transonic flow (the local velocity u along x direction is nondimensionalized by the farfield pressure p∞and density ρ∞).

Fig.7 Schematic diagram of integrating amplification factor N of perturbations within different segments.

where Δs=s-s0; λris real part of the eigenvalue obtained from the DMD analysis. Note that the integral calculus is based on the assumption that the spatial growth is a constant in the interval Δs.

Finally, the amplification factor N is computed for each segment and a transition location is obtained following the same logic of an eNcriterion, which is depicted in Fig.7, in the figure,N1-N4refers to the amplification factor of each segment, Ncritis the critical factor and Stris the predicted transition location. After the predicted transition location, the flow is considered turbulent and the production term in the turbulence model is switched on.Note that the critical factor Ncritis empirically calibrated and its commonly accepted value is in the range from 6 to 12.

2.4. Dynamic mode decomposition

The DMD technique37-39using flow snapshots can find its theoretical background in Koopman analysis39of nonlinear dynamical systems.This method can extract dynamic information from a flowfield to extract the resulting flow structures by ranking its frequency or wavenumber content.A DMD can be applied in time or space (i.e. frequency or wavenumber), providing information about the temporal or spatial growth of perturbations, respectively. Retaining the spatial framework enables the analysis of the growth of perturbations with different spatial wavenumbers.Many applications and extensions of this method can be found. Compared with other transition prediction methods (e.g. LST), the DMD/eNmethod only requires the snapshots of a flowfield obtained by solving the steady laminar Navier-Stokes equations, without modelling of perturbations or simplification such as the parallel flow or linearization assumption.

where vi(i=1,2,...,N)represents the flowfield snapshot at the ith station.We also assume that the sampling space Δs between the ordered sequences of data is a constant, which is required by a DMD theory. The spatial evolution of flow snapshots is sketched in Fig.8.

We assume that the flow snapshot viand the subsequent flow snapshot vi+1can be linked with a matrix A of a linear mapping (also known as Koopman assumption39), which is of the form

The sequence of the flow snapshots can be formulated as a Krylov sequence

Fig.8 Schematic diagram of flow snapshot evolution in space.

Eq. (6) may be written in a matrix form as

Eq.(8)shows that FDMDis the projection of A onto the subspace of modes of proper orthogonal decomposition(POD)in U. The eigenvalues and eigenvector of FDMDcan be expressed as FDMDyi=μiyi, where yidenotes the ith eigenvector and μidenotes the ith complex eigenvalue of FDMD. The eigenvalues of FDMDapproximate the eigenvalues of the matrix A with the largest growth. Finally, the dynamic modes φican be retrieved for the eigenvectors of FDMD,

which represents the coherent modal structures of the flowfield(e.g. the TS waves). The logarithmic mapping of the eigenvalues provides the spatial growth and wavenumber of modes:

where λrand λidenote the growth and wavenumber,respectively.

2.5. Higher-order dynamic mode decomposition

The Higher-Order Dynamic Mode Decomposition(HODMD)technique28is an extension of a DMD technique. It provides similar information (DMD modes, frequencies and growth rates), but has showed enhanced robustness in extracting dynamic information from the flowfields with noise.40An HODMD/eNmethod also requires snapshots of the flowfield obtained by solving the laminar NS equations to achieve transition predictions, without modelling of perturbations or simplification such as the parallel flow or linearization assumption.

The HODMD technique, also known as DMD-d,28relies on the high order Koopman assumption, in which the flow snapshots vi+dcan be represented by previous ‘‘d” flow snapshots as

where ‘‘d” is a tunable constant that gives the name to the order of an HODMD scheme: DMD-d. The parameter d can be compared with the number of sliding windows used in the well-known technique Power Spectral Density (PSD), which is a more accurate approximation of the classical Fast Fourier Transform (FFT). Thus, similarly to a PSD, which provides better performance than an FFT (more accurate and cleaner results), the performance of the HODMD is also better than the classical DMD. It is evident that when d=1, the method is the same as the classical DMD. It is possible to write the high order Koopman assumption in a matrix form as

Overall,DMD modes containing spatial growths and wave numbers are extracted using the HODMD as it is more robust.For simplicity,the interpretation below considers d=1.But it is readily extended for values of d larger than 1.

2.6. HODMD/eN method for transition prediction

In this subsection, we will make a theoretical derivation to show how a DMD analysis of the steady laminar flowfield can provide the information of spatial growth rate for unsteady small perturbations.

In Section 2.5,the matrix of flow dynamic system is written as Eq. (12).

Considering d=1 for simplicity in the dynamic system, it is assumed that there is a small perturbation δiat the ith station vi,which evolves along the flow direction.At the i+1 station,the perturbation becomes δi+1.If the perturbation is small enough, the flowfield where the small perturbations are added still satisfies the previous dynamic system, namely

Combining Eq. (12) and Eq. (15), we can get

It can be seen that the perturbations between adjacent stations are linearly linked by the matrix A, which represents the dynamics of the base flow. The eigenvalues μ of the matrix A can be expressed as

where i is the complex number unit. The eigenvalue μ and the corresponding eigenvector φ of the matrix A reflect the spatial evolution of the perturbations.

Referring to the LST method, we assume that the expression of a small perturbation is constructed as

where φkis the kth eigenvector of the matrix A,corresponding to a small perturbation of the time frequency ωk,s denotes the stream-wise coordinate along the airfoil surface, and t is the term of time.Substituting Eq.(18)into Eq.(16)and combining Eq. (17) yield

which represents the perturbation at the (i+1)th station.

According to Eq.(18)and Eq.(19),the growth of perturbations between adjacent stations can be written as

Hence, in a certain flow area within the airfoil boundary layer, the growth rate of the perturbation is

From Eq. (21), one can conclude that the real part of the logarithmically mapped eigenvalues obtained by a higherorder DMD analysis is exactly the growth rate of the perturbation at a certain frequency in each of the segments. This is a very critical conclusion, which exposes the theoretical background for the proposed HODMD/eNmethod.With the above derivation,one can explain why a DMD analysis of the steady laminar flowfield can be used to quantify the spatial growth of unsteady small perturbation.

3. Validation and examples

3.1. Validation of HODMD

We use the following test function to validate our HODMD code.

which contains 4 growing spatial waves with fixed wavenumbers. In this test case, the spatial complexity is 1(single point)and the spectral complexity is 4 (4 different modes). Fig.9 shows the spatial evolution (theoretical solution) of this test function. Note that 200 snapshots are used which are located at an equal space of 0.0025 along the x direction for the HODMD analysis.

Fig.9 Diagram of theoretical solution of test function.

From Table 1, we observe that the spatial HODMD technique provides exact results for the growth and wavenumber of 4 spatially growing perturbations. Through comparison of DMD and HODMD methods for the reconstruction of the function (see Fig.10) as well as the relative errors between the reconstructed snapshots and the theoretical solution (see Fig.11), it is shown that the HODMD can accurately reconstruct the true function but the DMD fails in this case. When the spatial complexity is smaller than the spectral complexity,the HODMD technique can capture the complete spatial dimension of the system(based on the Takens’delayed embedding theorem41), which explains the good performance of this higher-order assumption against its low-order counterpart(see more details in Ref. 28). For example, in the analysis of noisy data(i.e.obtained in experimental measurements),some of the DMD modes are employed to represent the noisy spectrum,and consequently altering the real spatial complexity of the solution.42A similar behavior is found in the analysis of transitory solution,in which some of the modes represent the transient dynamics,43and in flows conformed with a large number of small amplitudes-small scale spatial structures(as in transition to turbulent flow).40

An HODMD/eNmethod is based on the evolution of perturbations with small amplitudes, along the flow direction.The HODMD proves to be useful in that, on one hand, the presence of small perturbations alters the spatial complexity of the flow, and on the other hand, it is necessary to use a robust and accurate tool to capture small amplitude modes.These two aspects make the HODMD a suitable and unique technique for this type of analysis.

3.2. Transition prediction of flow over a low-speed NLF airfoil

In this case, we compute the transition location at different angles of attacks for a subsonic airfoil. The flowfield is obtained by solving RANS equations with prescribed transitions over the upper and lower surface to ensure sufficient laminar region.The Mach number is 0.1 and the Reynolds number is 4×106(based on the airfoil chord). The CFD grid for numerical simulation is sketched in Fig.12. The result at the angle of attack α=1° is described in detail. The flowfield is shown in Fig.13.

A convergence study of the grid size for CFD simulation is conducted. The predicted transition locations, the computed lift and drag coefficients and the growth and wavenumber information for the most unstable eigenvalue (with the largest positive growth rate) extracted by an HODMD (d=6) for segment [0, 0.2] are shown in Table 2.

It can be seen that force coefficients, predicted transition locations and the extracted wavenumber and growth rate ofthe most unstable perturbation converge with the increase of grid size. Therefore, the third grid at the size of 641×241 is selected for subsequent computations to ensure the accuracy of flow solution and DMD analysis. The distance of the first layer of grid normal to the airfoil surface is 1×10-6c, where c denotes the chord length, ensuring y+(1)＜1. Typical grids have 513 points along the airfoil and 241 points in the direction normal to the airfoil surface.

Table 1 Comparison of eigenvalues calculated by HODMD with theoretical values.

Fig.10 Comparison of reconstructed snapshots using HODMD and DMD.

Fig.11 Comparison of reconstruction relative error for u(x)using HODMD and DMD.

Fig.12 NLF0416 airfoil: C-type CFD grid.

Fig.13 NLF0416 airfoil: pressure contour and streamline computed by free-transition simulation (the local pressure p is non-dimensionalized by the farfield pressure p∞).

The more accurate the Navier-Stokes solution is, the less numerical noise there will be in DMD analysis. To check if the base flow is sufficiently converged, a convergence study on the flow solution is also conducted, as shown in Table 3.The resulting N factors are shown in Fig.14. It can be seen that for a threshold of averaged relative density residuals(ratio of the latest and initial density residual)of 10-6,the force coefficients, growth rates and predicted transition location change very little.In what follows,a residual tolerance of 10-6will be retained.

A high-resolution grid and a sufficiently converged steady flowfield ensure the accurate capture of TS waves. An estimation of the TS wavelength can be precomputed using the Blasius solution for flows over a laminar flat plate with the relationship of the minimum wavelength of small perturbations λTSto the boundary layer thickness δ: λTS/δ ≈6, see chapter XVI of the monograph by Schlichting and Gersten.44The boundary layer thickness can be estimated by δ=5 (νx/U)1/2≈0.0028c (at x=c=1, where c denotes the chord length, for a velocity U=Ma×340 m/s=34 m/s and Re=4.0×106, which lead to a kinematic viscosity ν=1×10-5m/s). Therefore, we approximately have a TS wavelength λTS≈0.017c. The use of 256 grid points and 200 snapshots along upper or lower surface of the airfoil has a spacing of 0.005c, which is fine enough to capture the minimum TS waves. Additionally, the Blasius solution provides a preliminary wavenumber λi_estimate=2π/λTS≈380 of the perturbations,which helps to locate the important physical eigenvalues when using a DMD method.

Table 2 NLF0416 airfoil: mesh convergence for force coefficients and transition prediction of upper surface.

Table 3 NLF0416 airfoil:influence of flow solution residual on force coefficients and predicted transition locations on upper surface of airfoil.

Fig.14 NLF0416: influence of solution residuals of RANS solver on N-factor obtained by HODMD.

According to a convergence study on the number of segments as shown in Fig.15, it is feasible to divide the upper and lower surface of the airfoil into 10 or more segments.The L segments (L=10) are divided over the airfoil from the stagnation point to the trailing edge point, as [0, 1/L],[0, 2/L], [0, 3/L], ..., [0, (L-1)/L], [0, 1]. The velocity on all segments [0-i/L] (i=1-10) will be inputted for an HODMD analysis to extract the eigenvalues.

The steady flow information is interpolated from the CFD grid to the equally-spaced snapshot grids for HODMD analysis. We choose the velocity (u, v) as the input information for the HODMD analysis. Note that the generation of snapshot grids has already been detailed in Section 2.3.2.To clarify this procedure,we depict the eigenvalues of dynamic modes as well as the growth rate and wavenumber information of perturbations for segment [0, 0.2] in Figs. 16 and 17, respectively. We choose the most unstable eigenvalues for this segment and retain the wavenumber λi=392.7 and the growth rate λr=13.1. Note that the value for the retained wavenumber is close to the estimated value λi_estimate≈380 by the Blasius solution, which convinces us that the results are reliable.

As it has been mentioned previously, we can perform HODMD analysis for each of the segments (i.e. [0, 0.1], [0,0.2], [0, 0.3], ...) and obtain the physical eigenvalues. The Nfactor can then be calculated according to the formulation explained above.

Fig.15 NLF0416:influence of number of segments on predicted transition location by HODMD.

Fig.16 NLF0416 airfoil: diagram of eigenvalues of dynamic modes for segment [0,0.2].

Fig.17 NLF0416 airfoil: growth rate and wavenumber information of perturbations for segment [0,0.2].

Tables 4 and 5 show the perturbation’s amplification factor N at different segments for the upper and lower surfaces of the NLF0416 airfoil, respectively. In the tables, s/c denotes the stream-wise coordinate on the airfoil surface and x/c denotes the stream-wise coordinate along the chord. The N factor represents the largest growth of perturbations in the specific segment. We choose Ncr=6 as the transition threshold.Table 5 shows that the N factor is beyond 6 in the segment[0,0.5],but not in the segment[0,0.6],which means that transition takes place between [0.5, 0.6]. There is no need to perform an HODMD analysis on the region after the transition(x/c ＞0.6) because the flow is modelled turbulent. A more precise location between[0.5,0.6]can be found by subdividing into smaller intervals, e.g. [0, 0.52] and [0, 0.54]. However, we find that using linear interpolation from the N factors obtained in the segments [0, 0.5] and [0, 0.6] provides a reasonable approximation.33,34The growth of perturbations obtained by the HODMD/eNand LST/eNmethods are shown in Figs. 18 and 19. We can find that the envelope of the N factor curve by the HODMD/eNis similar to that by the LST/eNmethod.

Table 6 compares the predicted transition location with the experiments. The computational cost of different methods is also compared. Note that the predicted transition locations agree well with the experiments, which proves the reliability of our method.Additionally,the computational cost is almost identical for both the methods.For completeness,the pressure coefficient Cpand skin friction Cfdistribution resulting from the HODMD/eNmethod are presented in Figs. 20 and 21,respectively, which agree well with the results obtained by the LST/eNmethod.

Finally, the predicted transition location Xtr_upand Xtr_lowat different angles of attack are displayed and compared with the experimental results, the published data45and the results obtained by the LST/eNmethod in Figs.22 and 23.The figures show that the transition locations at different angles of attack agree well with the other results, which suggests that the HODMD/eNmethod is a robust and accurate tool to predict the transition dominated by TS waves.

3.3. Transition prediction of flow over a transonic NLF airfoil

Simulation of flow over an NPU-LSC-72613 airfoil10is performed to demonstrate the proposed HODMD/eNmethod for transonic flows.In this case,the steady flowfield is obtained by solving the RANS equations with transition location fixed at 70% chord length on the upper and lower surfaces of the airfoil, at a free-stream Mach number of 0.72, Reynolds number of 2×107(based on the chord of airfoil) and angle of attack 0.5772°.The SA turbulence model is used after the transition points. A convergence study of grid size for the CFD simulation is also presented.The CFD grid for numerical simulation is illustrated in Fig.24.The distance of the first layer of grid normal to the airfoil surface is 1×10-6c,where c denotes the chord length, ensuring y+(1)＜1. Typical grids have 513 points along the airfoil and 241 points in the direction normal to the airfoil surface. The flowfield is illustrated in Fig.25,where a shock wave can be overserved on the suction side.

Table 4 NLF0416 airfoil: N factor calculated by HODMD/eN method on upper surface (Ma=0.1, Re=4.0×106, α=1).

Table 5 NLF0416 airfoil: N factor calculated by HODMD/eN method on lower surface (Ma=0.1, Re=4.0×106, α=1°).

Fig.18 NLF0416 airfoil at α=1: comparison of N factors by HODMD/eN and LST/eN methods on upper surface.

Fig.19 NLF0416 airfoil at α=1: comparison of N factors by HODMD/eN and LST/eN methods on lower surface.

Fig.20 NLF0416 airfoil at α=1: comparison of computed pressure distributions with experimental data.45

Fig.21 NLF0416 airfoil at α=1 comparison of computed skin friction distributions obtained by different methods.

Table 6 NLF0416 airfoil: comparison of predicted transition locations and time costs for upper and lower surface (Ma=0.1,Re=4.0×106, α=1°).

Fig.22 NLF0416 airfoil: comparison of predicted transition location at different angles of attack on airfoil upper surface.

Fig.23 NLF0416 airfoil: comparison of predicted transition location at different angles of attack on airfoil lower surface.

Fig.24 NPU-LSC-72613 airfoil: C-type CFD grid.

Fig.25 NPU-LSC-72613 airfoil: pressure contours and streamline computed for free-transition using HODMD/eN method.

The results by using the HODMD/eNmethod agree well with the results from the other methods, which validates the proposed HODMD/eNmethod. Transition on the upper surface of the airfoil is attributed to the presence of shock wave.A minor difference is observed for the‘‘FLUENT”prediction,which is attributed to the prediction of shock wave location.The shock-wave-induced transition is typical for laminar airfoils and provides an interesting test case for our method.Note that the lower surface is dominated by TS instability and in this case all methods show comparable transition locations.For completeness,Fig.26 depicts the predicted pressure distributions and transition locations for all methods. This case illustrates the capability of our HODMD/eNmethod to compute transition prediction of compressible flows over transonic airfoils.

4. Conclusions

A novel transition prediction approach called HODMD/eNwas proposed in this article. The new method uses a spatial high order DMD (HODMD) analysis to extract the modes of unsteady small perturbations from a steady laminar flowfield and substitute the LST analysis in an LST/eNmethod to quantify the growth of Tollmien-Schlichting (TS) waves.Once the growth rate and wavenumber are obtained,an eNcriterion is employed to determine the transition point(s), which is used to activate the turbulence model in a RANS simulation.The HODMD code was validated by an analytical test case and the proposed HODMD/eNwere validated and demonstrated by transition predictions of flows over a low-speed NLF0416 airfoil and a transonic NPU-LSC-72613 airfoil.Some conclusions follow:

(1). The proposed HODMD/eNmethod has advantages as it directly extracts the growth and wavenumber information of small perturbations from the simulated flow,instead of modelling the perturbations upon a base flow(as in an LST/eNmethod). In addition, The HODMD/eNmethod does not explicitly rely on the parallel flow assumptions, which can help to generalize the domainof applicability, for example when in studying highly curved surfaces,surface imperfections or secondary flow instabilities.

Table 7 NPU-LSC-72613 airfoil: comparison of transition location and lift/drag coefficient among different solvers (Ma=0.72,Re=2.0×107, α=0.5772°).

Fig.26 NPU-LSC-72613: comparison of pressure distributions and predicted transition locations obtained by different flow solvers and transition prediction methods.

(2). Theoretical derivation proves how a spatial HODMD analysis can extract the eigenvalues and eigenmodes of unsteady small perturbations from a steady laminar flowfield.It is shown that the real part of the logarithmically mapped eigenvalue is exactly the growth rate of the small perturbation at a certain frequency in each of the segments along the boundary of an airfoil, and the real part denotes the wavenumber.

(3). An HODMD is more robust than the classic DMD for the problems whose spatial complexity is smaller than the spectral complexity and more robust to capture physical information embedded in noise.

The test cases treated in this work may be grouped in transitional scenarios dominated by the spatial growth of TS waves. However, since the proposed method relies on direct simulation of laminar flow regions, it is expected that the method can be extended to compute more complex transition scenarios, such as dominated by non-linear interactions, e.g.cross-flow transition over a 3D swept wing. We have done some preliminary study,47which showed interesting results.In addition, the proposed method can be extended to transition prediction for supersonic or hypersonic flows. These will be studied in future work.

Acknowledgements

This work was supported by the National Natural Science Foundation of China (No. 11772261), the Aeronautical Science Foundation of China (No. 2016ZA53011), the ATCFD Project (No. 2015-F-016) and the 111 Project of China (No. B17037). The author would like to thank Prof.Eusebio Valero, Ms. Liu Fangliang, Ms. Han Li and Ms.Wang Shaonan for their guidance and technical support.