# Efficient Training Data Generation for Reduced-Order Modeling in a Transonic Flight Regime

Authors: Haojie Liu; Yonghui Zhao
Journal: International Journal of Aerospace Engineering
 (2018)

---

## Abstract

In this study, a time-dependent surrogate approach is presented to generate the training data for identifying the reduced-order model of an unsteady aerodynamic system with the variation of mean angle of attack and Mach number in a transonic flight regime. For such a purpose, a finite set of flight samples are selected to cover the flight range of concern at first. Subsequently, the unsteady aerodynamic outputs of the system under given inputs of filtered white Gaussian noise at these flight samples are simulated via CFD technique which solves Euler equations. The unsteady aerodynamic outputs, which are viewed as a time-dependent function of flight parameters, can be approximated via the Kriging technique at each time step. By this way, the training data for any combination of flight parameters in the range of concern can be obtained without performing any further CFD simulations. To illustrate the accuracy and validity of the training data generated via the proposed approach, the constructed data are used to identify the reduced-order aerodynamic models of a NACA 64A010 airfoil via a robust subspace identification algorithm. The unsteady aerodynamics and aeroelastic responses under various flight conditions in a transonic flight regime are computed. The results agree well with those obtained by using the training data of CFD technique.


---

## Body

## 1. Introduction

The techniques of computational fluid dynamics (CFD) have been widely used to simulate both linear and nonlinear flow fields for various flight vehicles. However, it is still time-consuming for any high-fidelity CFD techniques to simulate the unsteady aerodynamic loads due to the broad variation of parameters, such as different combinations of mean angle of attack and Mach number. For example, the linear doublet-lattice method has been frequently used for the aeroelastic analysis of aircraft although CFD techniques offer more accurate numerical simulations. In a transonic flight regime, however, the double-lattice method does not work properly because of the aerodynamic nonlinearities coming from shock waves and flow separations [1]. Thus, various CFD-based reduced-order models (ROMs) have been developed [2] so as to provide an effective way to simulate unsteady aerodynamic loads with a high level of accuracy.In general, the unsteady aerodynamic ROMs can be classified into frequency-domain type and time-domain type. The frequency-domain ROMs mainly contain the proper orthogonal decomposition (POD) approach [3], balanced POD [4], and the harmonic balance approach [5]. On the other hand, the time-domain ROMs include the Eigensystem realization algorithm (ERA) [6], the Volterra theory [7], the auto regressive-moving-average (ARMA) model [8], surrogate models via artificial neural networks [9–11] and Kriging technique [12], and the Wiener-type cascade model [13].Although a significant progress has been made for the unsteady aerodynamic ROMs, almost all aforementioned ROMs are only valid for a set of fixed flight parameters. It is still time-consuming to generate the unsteady aerodynamic ROMs for a range of flight parameters since the same procedure of computation has to be repeated under each flight condition in the range of concern. Recently, several reduced-order modeling approaches have been proposed to predict unsteady aerodynamic loads for a range of flight parameters. For example, a surrogate-based recurrence framework ROM was developed to model the unsteady aerodynamics on a rotating airfoil [14]. Yet, it is time-consuming to generate the training data for each combination of parameters in the parameter space. A ROM adaptation approach based on the interpolation in a tangent space to a Grassmann manifold was proposed to predict the aeroelastic characteristics for an F-16 configuration [15]. To obtain new basis vectors for a new flight condition via interpolation, the size of the basis vectors of different local ROMs should be the same so that the accuracy of the POD approach may be decreased. A ROM approach by combining linear convolution with a nonlinear correction factor was proposed to model the aerodynamic characteristics for multiple Mach regimes [16]. However, the approach gives a distinct phase shift at higher oscillation frequencies. A time-dependent surrogate model to fit the relationship between flight parameters and step functions was proposed to model the unsteady and nonlinear aerodynamic loads [17]. Only one kind of input, the pitching motion of the aircraft in their study, can be taken into account at each time. A Kriging surrogate model was proposed to model the unsteady aerodynamic forces with respect to a range of Mach numbers [18]. However, the computational cost for Kriging model increased significantly in order to take the Mach number into consideration. Recently, a reduced-order modeling approach based on recurrent local linear neuro-fuzzy models was developed to model the generalized aerodynamic forces over a range of Mach numbers [19]. To guarantee the accuracy of the local linear neuro-fuzzy model, the flight parameters for training data should be selected carefully in the flight range of concern. Very recently, a nonlinear interpolation method based on a set of local linear state-space models was proposed to model the generalized aerodynamic forces for an elastic wing with control surfaces [20]. Nevertheless, enough local linear state-space models were required to cover the flight range of concern and only Mach number variations were taken into consideration.Almost all aforementioned unsteady aerodynamic ROMs, which take flight parameter variations into account, require enough training data to capture the dynamic characteristics of the unsteady aerodynamic systems. However, it is time-consuming to generate the training data via any direct CFD simulations when both Mach number and mean angle of attack are taken into consideration. The motivation of this study is to generate the training data efficiently for a range of flight parameters including mean angle of attack and Mach number in a transonic flight regime. For such a purpose, a time-dependent surrogate approach is proposed. Once the Kriging models are constructed via the proposed approach, the training data under an arbitrary flight condition in the range of concern can be obtained without performing any further CFD simulations. The remainder of the paper is organized as follows. In Section2, the theoretical background of the training data generation approach and the robust subspace algorithm is presented. In Section 3, the unsteady aerodynamics and aeroelastic responses of a NACA 64A010 airfoil under different flight conditions in a transonic flight regime are investigated to validate the proposed approach. In Section 4, some conclusions are given.


## 2. Theoretical Background

In this section, the theoretical background of the training data generation approach and the robust subspace algorithm is presented. The flowchart of the reduced-order modeling approach is shown in Figure1. The first step is to obtain the unsteady aerodynamic outputs of the system via direct CFD simulations at a finite set of flight samples selected to cover the flight range of concern. Then, a time-dependent surrogate method based on the Kriging technique is used to approximate the relationship between the unsteady aerodynamic outputs of the system and the flight parameters. Afterwards, with the training data obtained, the robust subspace algorithm is implemented to identify the discrete-time state-space model under an arbitrary flight condition.Figure 1
The flowchart of the reduced-order modeling approach.


### 2.1. Selection of Excitation Signal

As well known, the dynamic characteristics of an unsteady aerodynamic system of concern should be appropriately excited via the input signal in order to construct the training data for successfully identifying the corresponding time-domain ROMs. For this purpose, the filtered white Gaussian noise (FWGN) is used as the input signal to make CFD simulations. Under such a type of excitation, only the frequency response spectrum of concern can be excited and there is no need to repeat computations for each combination of frequency and amplitude.Once the type of input signal is determined, it is necessary to select the sampling rate and the sampling length for the input signal. The sampling ratefS determined by the time step Δt of the CFD solver is given as follows:
(1)fS=1Δt.According to the Nyquist sampling theorem, the Nyquist frequencyfN is fS/2 which must exceed the maximal frequency of concern. Furthermore, the frequency resolution of the input signal with the sampling length N, given by fN/N, should be adequate for the representation of the minimal frequency of concern. Hence, the time step Δt and the sampling length N for the input signal need to be adjusted carefully.Although only one unsteady aerodynamic simulation under FWGN excitation is required under each flight condition, it is still time-consuming to generate the training data to capture the dynamic characteristics of an unsteady aerodynamic system in the mean angle of attack and Mach number space.


### 2.2. Time-Dependent Surrogate Approach for Training Data Generation

To generate the training data for a range of flight parameters efficiently, a time-dependent surrogate method based on the Kriging technique is proposed in this subsection. Similar to the previous study [17], the unsteady aerodynamic outputs of the system under the FWGN excitation and a finite set of flight conditions are considered as a set of time-correlated spatial processes. Under such an assumption, the surrogate techniques can be used to model the time-dependent outputs of the aerodynamic system at each time step as a function of flight parameters. Here, the Kriging technique is implemented since it has good accuracy and robustness with small data sets [21].To approximate the input-output relationship between the flight parameters and the unsteady aerodynamic outputs of the system,M flight samples need to be selected to cover the range of concern at first. The design of experimental studies [22] can be used to select these samples from the mean angle of attack α0 and Mach number M∞ space. It is worthy to note that the selection of the input parameters strongly depends on the flight regime of concern and the aircraft configuration. For example, more flight parameter samples need to be placed at the region where there exist strong aerodynamic nonlinearities in the flow field. Subsequently, the input vector X for the Kriging technique can be defined as
(2)X=M∞,1α0,1M∞,2α0,2⋮⋮M∞,Mα0,M.Under each flight condition defined byM∞,i and α0,i, an unsteady aerodynamic computation needs to be performed via direct CFD simulation under the given FWGN excitation. The corresponding outputs of the aerodynamic system of the M flight samples can be written as
(3)Yo=y1,1y2,1⋯yN,1y1,2y2,2⋯yN,2⋮⋮⋱⋮y1,My2,M⋯yN,M,where yi,j is the aerodynamic output vector and N is the number of time steps of the FWGN excitation. With the input vector X and the output vector Yo known, the Kriging technique can be used to approximate the relationship between the flight parameter and the aerodynamic outputs of the system at each time step.What follows is a brief description of the Kriging technique, where the Kriging model approximates the target function at an untried sitex∗ as
(4)Φx∗=fx∗β+zx∗,Here,fx∗ and β are the regression model and the regression coefficients, respectively. In (4), the first term is the mean value, which can be thought as a globally valid trend function. The second term zx∗ given in (4) is a Gaussian distributed error term with zero mean and variance σ2. To construct the Kriging model, the values of the regression coefficients β can be approximated by solving a least-square regression problem as follows:
(5)β^=FTR−1F−1FTR−1h,where R is the spatial correlation matrix. h and F are the M×1 output vector and the M×M mapping matrix at the sample inputs, respectively. h and F are defined as
(6)h=y1y2⋮yM,F=f1x1f2x1⋯fr+1x1f1x2f2x2⋯fr+1x2⋮⋮⋱⋮f1xMf2xM⋯fr+1xM.With the vectorr=Rx∗,x1Rx∗,x2⋯Rx∗,xMT, the prediction at an unsampled location x∗ can be obtained as
(7)Φ^x∗=fx∗β^+rTR−1h−Fβ^.It has to be noted that one mapping functionΦ^ needs to be constructed for each component of the output vector at each time step. Once the Kriging models are constructed, the unsteady aerodynamic outputs of the system under an arbitrary flight condition in the range of concern can be obtained without performing any further CFD simulations. For example, the corresponding aerodynamic outputs of the system under the same FWGN excitation for the mean angle of attack α0∗ and Mach number M∞∗ can be written as
(8)y1∗y2∗⋯yN∗=Φ^1M∞∗,α0∗Φ^2M∞∗,α0∗⋯Φ^NM∞∗,α0∗.Hence, the training data under an arbitrary flight condition in the range of concern can be obtained efficiently via the proposed approach.


### 2.3. Robust Subspace Algorithm

With the training data obtained via the time-dependent surrogate approach, a time-domain ROM approach can be used to construct the unsteady aerodynamic ROM. In this study, the robust subspace algorithm [23] is implemented to identify the discrete-time state-space models. This algorithm is always numerically stable and convergent and has no identification problems, such as lack of convergence, slow convergence, or numerical instability.A brief description of the algorithm is given as follows. The discrete-time state-space model assumes the following model structure:(9)xak+1=Axak+Buk,yk=Cxak+Duk,where xa, u, and y are the vectors of state, input, and output variables, respectively. Given the training data, four constant matrices A, B, C, and D are estimated such that the unsteady aerodynamic outputs of the system under the FWGN excitation can be reproduced.At first, the input and output Hankel block matrices based on the training data are defined as(10)UpUf=defu0u1⋯uj−1u1u2⋯uj⋯⋯⋯⋯ui−1ui⋯ui+j−2uiui+1⋯ui+j−1ui+1ui+2⋯ui+j⋯⋯⋯⋯u2i−1u2i⋯u2i+j−2,YpYf=defy0y1⋯yj−1y1y2⋯yj⋯⋯⋯⋯yi−1yi⋯yi+j−2yiyi+1⋯yi+j−1yi+1yi+2⋯yi+j⋯⋯⋯⋯y2i−1y2i⋯y2i+j−2,where the number j of columns is given as
(11)j=N−2i+1.In (11), N is the length of training data and i is the number of block rows in the Hankel block matrices. Then, based on the RQ decomposition of the Hankel block matrix, the weighted oblique projection οiΠUf⊥ can be computed. The singular value decomposition of the weighted oblique projection is given as follows:
(12)οiΠUf⊥=USVT.The order of state-space model can be determined by inspecting the singular values inS. Afterwards, U1 and S1 can be obtained according to the order and the matrices. The extended observability matrix Γi can be computed as
(13)Γi=U1S11/2.Here,Γi−1 denotes the matrix Γi without the last p rows. To this end, the system matrices A, B, C, and D can be obtained with the matrices Γi, Γi−1, and the decomposition matrix of the Hankel block matrix. For further details, one can refer [23].In summary, the efficiency of the algorithm enables one to identify the unsteady aerodynamic ROM quickly under an arbitrary flight condition.


## 2.1. Selection of Excitation Signal

As well known, the dynamic characteristics of an unsteady aerodynamic system of concern should be appropriately excited via the input signal in order to construct the training data for successfully identifying the corresponding time-domain ROMs. For this purpose, the filtered white Gaussian noise (FWGN) is used as the input signal to make CFD simulations. Under such a type of excitation, only the frequency response spectrum of concern can be excited and there is no need to repeat computations for each combination of frequency and amplitude.Once the type of input signal is determined, it is necessary to select the sampling rate and the sampling length for the input signal. The sampling ratefS determined by the time step Δt of the CFD solver is given as follows:
(1)fS=1Δt.According to the Nyquist sampling theorem, the Nyquist frequencyfN is fS/2 which must exceed the maximal frequency of concern. Furthermore, the frequency resolution of the input signal with the sampling length N, given by fN/N, should be adequate for the representation of the minimal frequency of concern. Hence, the time step Δt and the sampling length N for the input signal need to be adjusted carefully.Although only one unsteady aerodynamic simulation under FWGN excitation is required under each flight condition, it is still time-consuming to generate the training data to capture the dynamic characteristics of an unsteady aerodynamic system in the mean angle of attack and Mach number space.


## 2.2. Time-Dependent Surrogate Approach for Training Data Generation

To generate the training data for a range of flight parameters efficiently, a time-dependent surrogate method based on the Kriging technique is proposed in this subsection. Similar to the previous study [17], the unsteady aerodynamic outputs of the system under the FWGN excitation and a finite set of flight conditions are considered as a set of time-correlated spatial processes. Under such an assumption, the surrogate techniques can be used to model the time-dependent outputs of the aerodynamic system at each time step as a function of flight parameters. Here, the Kriging technique is implemented since it has good accuracy and robustness with small data sets [21].To approximate the input-output relationship between the flight parameters and the unsteady aerodynamic outputs of the system,M flight samples need to be selected to cover the range of concern at first. The design of experimental studies [22] can be used to select these samples from the mean angle of attack α0 and Mach number M∞ space. It is worthy to note that the selection of the input parameters strongly depends on the flight regime of concern and the aircraft configuration. For example, more flight parameter samples need to be placed at the region where there exist strong aerodynamic nonlinearities in the flow field. Subsequently, the input vector X for the Kriging technique can be defined as
(2)X=M∞,1α0,1M∞,2α0,2⋮⋮M∞,Mα0,M.Under each flight condition defined byM∞,i and α0,i, an unsteady aerodynamic computation needs to be performed via direct CFD simulation under the given FWGN excitation. The corresponding outputs of the aerodynamic system of the M flight samples can be written as
(3)Yo=y1,1y2,1⋯yN,1y1,2y2,2⋯yN,2⋮⋮⋱⋮y1,My2,M⋯yN,M,where yi,j is the aerodynamic output vector and N is the number of time steps of the FWGN excitation. With the input vector X and the output vector Yo known, the Kriging technique can be used to approximate the relationship between the flight parameter and the aerodynamic outputs of the system at each time step.What follows is a brief description of the Kriging technique, where the Kriging model approximates the target function at an untried sitex∗ as
(4)Φx∗=fx∗β+zx∗,Here,fx∗ and β are the regression model and the regression coefficients, respectively. In (4), the first term is the mean value, which can be thought as a globally valid trend function. The second term zx∗ given in (4) is a Gaussian distributed error term with zero mean and variance σ2. To construct the Kriging model, the values of the regression coefficients β can be approximated by solving a least-square regression problem as follows:
(5)β^=FTR−1F−1FTR−1h,where R is the spatial correlation matrix. h and F are the M×1 output vector and the M×M mapping matrix at the sample inputs, respectively. h and F are defined as
(6)h=y1y2⋮yM,F=f1x1f2x1⋯fr+1x1f1x2f2x2⋯fr+1x2⋮⋮⋱⋮f1xMf2xM⋯fr+1xM.With the vectorr=Rx∗,x1Rx∗,x2⋯Rx∗,xMT, the prediction at an unsampled location x∗ can be obtained as
(7)Φ^x∗=fx∗β^+rTR−1h−Fβ^.It has to be noted that one mapping functionΦ^ needs to be constructed for each component of the output vector at each time step. Once the Kriging models are constructed, the unsteady aerodynamic outputs of the system under an arbitrary flight condition in the range of concern can be obtained without performing any further CFD simulations. For example, the corresponding aerodynamic outputs of the system under the same FWGN excitation for the mean angle of attack α0∗ and Mach number M∞∗ can be written as
(8)y1∗y2∗⋯yN∗=Φ^1M∞∗,α0∗Φ^2M∞∗,α0∗⋯Φ^NM∞∗,α0∗.Hence, the training data under an arbitrary flight condition in the range of concern can be obtained efficiently via the proposed approach.


## 2.3. Robust Subspace Algorithm

With the training data obtained via the time-dependent surrogate approach, a time-domain ROM approach can be used to construct the unsteady aerodynamic ROM. In this study, the robust subspace algorithm [23] is implemented to identify the discrete-time state-space models. This algorithm is always numerically stable and convergent and has no identification problems, such as lack of convergence, slow convergence, or numerical instability.A brief description of the algorithm is given as follows. The discrete-time state-space model assumes the following model structure:(9)xak+1=Axak+Buk,yk=Cxak+Duk,where xa, u, and y are the vectors of state, input, and output variables, respectively. Given the training data, four constant matrices A, B, C, and D are estimated such that the unsteady aerodynamic outputs of the system under the FWGN excitation can be reproduced.At first, the input and output Hankel block matrices based on the training data are defined as(10)UpUf=defu0u1⋯uj−1u1u2⋯uj⋯⋯⋯⋯ui−1ui⋯ui+j−2uiui+1⋯ui+j−1ui+1ui+2⋯ui+j⋯⋯⋯⋯u2i−1u2i⋯u2i+j−2,YpYf=defy0y1⋯yj−1y1y2⋯yj⋯⋯⋯⋯yi−1yi⋯yi+j−2yiyi+1⋯yi+j−1yi+1yi+2⋯yi+j⋯⋯⋯⋯y2i−1y2i⋯y2i+j−2,where the number j of columns is given as
(11)j=N−2i+1.In (11), N is the length of training data and i is the number of block rows in the Hankel block matrices. Then, based on the RQ decomposition of the Hankel block matrix, the weighted oblique projection οiΠUf⊥ can be computed. The singular value decomposition of the weighted oblique projection is given as follows:
(12)οiΠUf⊥=USVT.The order of state-space model can be determined by inspecting the singular values inS. Afterwards, U1 and S1 can be obtained according to the order and the matrices. The extended observability matrix Γi can be computed as
(13)Γi=U1S11/2.Here,Γi−1 denotes the matrix Γi without the last p rows. To this end, the system matrices A, B, C, and D can be obtained with the matrices Γi, Γi−1, and the decomposition matrix of the Hankel block matrix. For further details, one can refer [23].In summary, the efficiency of the algorithm enables one to identify the unsteady aerodynamic ROM quickly under an arbitrary flight condition.


## 3. Case Studies

This section presents the study on the unsteady aerodynamic problem and aeroelastic problem of a NACA 64A010 airfoil, respectively, for the mean angle of attackα0 ranging from 0 to 5 deg and the Mach number M∞ ranging from 0.75 to 0.85 so as to validate the proposed approach.The airfoil of concern has two degrees of freedom, as shown in Figure2, where x¯α is the dimensionless distance between the center of gravity and the stiffness center, r¯α is the gyration radius of the airfoil around the stiffness center, ωα and ωh are the uncoupled natural frequencies of the pitch and plunge, and b, μ, V∞, and V∞∗ are the half-chord length, mass ratio, free-stream velocity, and dimensionless velocity, respectively. With the dimensionless time τ=ωα⋅t, the dynamic equation of the airfoil can be established in the following state-space form.
(14)ẋ=As⋅x+Bs⋅fax,t.Figure 2
An airfoil of two degrees of freedom.The vectors and matrices in (14) can be written as
(15)As=02×2I2×2−M−1K02×2,Bs=02×2M−1,fa=V∞∗2π−Cl2Cm,ξ=hbα,x=ξξ̇,where the mass and stiffness matrices are defined as
(16)M=1x¯αx¯αr¯α2,K=ωhωα200r¯α2.The structural parameters of the airfoil are chosen to bex¯α=1.8, r¯α2=3.48, a¯=−2.0, ωh=100rad/sec, ωα=100rad/sec, and μ=60, respectively.


### 3.1. Training Data for the Unsteady Aerodynamic System

The CFD technique which solves the Euler equations [13, 18] is used to compute the unsteady aerodynamics over the airfoil in transonic regime. The spatial discretization of the CFD solver is based on the cell-centered finite-volume approach. The widely used central scheme with artificial dissipation is employed for convective flux calculation. The dual time-stepping approach is employed for unsteady simulations. After performing mesh and time step sensitivity analyses, the unstructured mesh with 4098 cells and dimensionless fluid time step 0.2 are used in the present study. Figure 3 presents the CFD grid consisting of 4098 triangle elements and 2145 nodes. The high-fidelity CFD solver needs to be validated at first. In order to validate the Euler solver used in the present study, the test case CT6 of [24] is studied. For such case, the free stream Mach number is 0.796 and the 64A010 airfoil pitches around its quarter-chord point and the pitching motion are governed as
(17)αt=α0+αAsinωt,where α0=0 deg and αA=1.01 deg are the mean angle of attack and the amplitude of pitching oscillation, respectively. The relationship between the angular frequency ω and the reduced frequency kr=0.202 is described as
(18)kr=ωc2V∞,where c and V∞ are the chord length and the free-stream velocity, respectively. As shown in Figure 4, the CFD lift coefficients agree well with the experimental data. However, there are apparent differences for the moment coefficients between the numerical and experimental results. Such discrepancy may be due to the high sensitivity of the moment coefficients to the location of shock wave.Figure 3
The CFD grid for the NACA64A010 airfoil.Figure 4
Lift and moment coefficient loops for a NACA 64A010 airfoil at Mach number 0.796.In the numerical simulations, the dimensionless bandwidth of the FWGN excitation signal is chosen to be 0.01~0.4. Figure5 presents the time histories of the two FWGN excitation signals selected to drive the pitch and plunge motions of the airfoil, respectively. The length of the input signal in CFD simulations is taken as 2000. Subsequently, 20 flight samples are selected for the mean angle of attack α0 in [0, 5] deg and for the Mach number M∞ in [0.75, 0.85], respectively, as shown in Figure 6. For each flight parameter sample, an unsteady aerodynamic simulation is performed via the CFD solver.Figure 5
The FWGN excitation signals to drive pitch and plunge motions.Figure 6
The flight parameter samples to cover the range of concern.According to the time-dependent surrogate approach, 4000 Kriging models are required since the length of the FWGN excitation is 2000 for two unsteady aerodynamic outputs of the system, that is, the lift coefficientCl and the moment coefficient Cm. To demonstrate the accuracy of the training data, Cl and Cm are predicted via the Kriging models under two flight conditions in a transonic flight regime, one with M∞=0.76 and α0=1 deg and the other with M∞=0.82 and α0=4 deg. As shown in Figures 7 and 8, the results predicted via the Kriging models agree well with those obtained via the direct CFD technique.Figure 7
The training data validation of a NACA 64A010 airfoil atM∞=0.76 and α0=1deg.Figure 8
The training data validation of a NACA 64A010 airfoil atM∞=0.82 and α0=4deg.The CPU time to construct the Kriging models is shown in Table1, where all the computations are performed by using a laptop with two CPUs of 2.4 GHz. Even though the total CPU time cost associated with the construction of the Kriging models is a little bit expensive, the training data under an arbitrary flight condition in the range of concern can be obtained efficiently without performing any further CFD simulation. For example, the computational time of the direct CFD approach which solves Euler equations is about 0.284 h in order to obtain the training data at M∞=0.76 and α0=3deg, while the computational time of the Kriging models for the same case is 0.006 h only.Table 1
The computational cost associated with the time-dependent surrogate approach.
CPU time
Computation of the 20 training cases
4.98 h
Generation of the Kriging models
0.01 h
Total CPU time
4.99 h


### 3.2. Unsteady Aerodynamic Prediction

To demonstrate the application of the training data in unsteady aerodynamic prediction, the constructed data are used to identify the state-space models via the robust subspace algorithm. The unsteady aerodynamic outputs of the system are computed under different excitation signals, and the results are compared with those obtained by using the CFD training data.At first, the FWGN excitation signals, not the same as the training signal during the ROM construction, are used to drive the airfoil atM∞=0.76 and α0=3deg. The pitch and plunge motions of the airfoil are shown in Figure 9. The corresponding unsteady aerodynamic coefficients are shown in Figure 10. For further comparison, the flight condition with M∞=0.82 and α0=5deg is also selected. The excitation signals and the corresponding unsteady aerodynamic results are presented in Figures 11 and 12, respectively. In addition, a harmonic oscillation is chosen as the input signal at M∞=0.84 and α0=0deg. The unsteady aerodynamic results are presented in Figure 12.Figure 9
The pitch and plunge motions under the FWGN excitation for ROM validation atM∞=0.76 and α0=3deg.Figure 10
The ROM validation under the FWGN excitation atM∞=0.76 and α0=3deg (ROM1-training data via CFD; ROM2-training data via the proposed method).Figure 11
The FWGN excitation signals for ROM validation atM∞=0.82 and α0=5deg.Figure 12
The ROM validation under the FWGN excitation atM∞=0.82 and α0=5deg (ROM1-training data via CFD; ROM2-training data via the proposed method).As can be seen from Figures10, 12, and 13, the time histories of Cl and Cm obtained via the ROM2s identified by using the time-dependent surrogate approach agree well with those of ROM1s identified by using the direct CFD simulation data. That is, the proposed approach is able to generate good training data for ROMs in a transonic flight regime.Figure 13
The ROM validation under the harmonic excitation atM∞=0.84 and α0=0deg (ROM1-training data via CFD; ROM2-training data via the proposed method).


### 3.3. Aeroelastic Time-Marching Simulation

To further demonstrate the performance of the proposed approach in the flutter prediction in a transonic flight regime, the ROM2s are coupled with the structural model to make time-marching simulations. At first, the airfoil is forced to make three complete cycles of harmonic motion in pitch at the frequencyωα/2. Then, the aeroelastic system is allowed to evolve by its own self-induced loads.In the present case, the aeroelastic responses of the airfoil predicted via the ROM2s and ROM1s are compared under two flight conditions. As shown in Figure14, the aeroelastic responses of the airfoil with M∞=0.81 and α0=0deg are computed for V∞∗=0.73. Both ROM2 and ROM1 can predict a neutrally stable oscillation. Afterwards, the aeroelastic responses of the airfoil with M∞=0.85 and α0=1deg are evaluated. As shown in Figure 15, the results of ROM2 and ROM1 present a neutrally stable oscillation at V∞∗=0.57.Figure 14
The aeroelastic responses of the airfoil forV∞∗=0.73 at M∞=0.81 and α0=0deg (ROM1-training data via CFD; ROM2-training data via the proposed method).Figure 15
The aeroelastic responses of the airfoil forV∞∗=0.57 at M∞=0.85 and α0=1deg (ROM1-training data via CFD; ROM2-training data via the proposed method).Finally, Figures16 and 17 give the flutter boundaries of the aeroelastic model with the mean angle of attack α0=1deg and α0=3deg, respectively. As shown in these two figures, the computational results based on ROM2s and ROM1s get a good agreement in the flutter boundaries. Nevertheless, the computation of ROM2 is much more efficient than that of ROM1 since no further CFD simulations are required once the Kriging models are constructed.Figure 16
The flutter boundary of the airfoil with the mean angle of attackα0=1deg (ROM1-training data via CFD; ROM2-training data via the proposed method).Figure 17
The flutter boundary of the airfoil with the mean angle of attackα0=3deg (ROM1-training data via CFD; ROM2-training data via the proposed method).As can be seen from the flutter boundaries of the airfoil at different mean angle of attack, the aeroelastic characteristics of the model are very sensitive to the mean angle of attack in a transonic flight regime. Hence, the effects of the variations of the mean angle of attack have to be taken into account for aeroelastic analysis.


## 3.1. Training Data for the Unsteady Aerodynamic System

The CFD technique which solves the Euler equations [13, 18] is used to compute the unsteady aerodynamics over the airfoil in transonic regime. The spatial discretization of the CFD solver is based on the cell-centered finite-volume approach. The widely used central scheme with artificial dissipation is employed for convective flux calculation. The dual time-stepping approach is employed for unsteady simulations. After performing mesh and time step sensitivity analyses, the unstructured mesh with 4098 cells and dimensionless fluid time step 0.2 are used in the present study. Figure 3 presents the CFD grid consisting of 4098 triangle elements and 2145 nodes. The high-fidelity CFD solver needs to be validated at first. In order to validate the Euler solver used in the present study, the test case CT6 of [24] is studied. For such case, the free stream Mach number is 0.796 and the 64A010 airfoil pitches around its quarter-chord point and the pitching motion are governed as
(17)αt=α0+αAsinωt,where α0=0 deg and αA=1.01 deg are the mean angle of attack and the amplitude of pitching oscillation, respectively. The relationship between the angular frequency ω and the reduced frequency kr=0.202 is described as
(18)kr=ωc2V∞,where c and V∞ are the chord length and the free-stream velocity, respectively. As shown in Figure 4, the CFD lift coefficients agree well with the experimental data. However, there are apparent differences for the moment coefficients between the numerical and experimental results. Such discrepancy may be due to the high sensitivity of the moment coefficients to the location of shock wave.Figure 3
The CFD grid for the NACA64A010 airfoil.Figure 4
Lift and moment coefficient loops for a NACA 64A010 airfoil at Mach number 0.796.In the numerical simulations, the dimensionless bandwidth of the FWGN excitation signal is chosen to be 0.01~0.4. Figure5 presents the time histories of the two FWGN excitation signals selected to drive the pitch and plunge motions of the airfoil, respectively. The length of the input signal in CFD simulations is taken as 2000. Subsequently, 20 flight samples are selected for the mean angle of attack α0 in [0, 5] deg and for the Mach number M∞ in [0.75, 0.85], respectively, as shown in Figure 6. For each flight parameter sample, an unsteady aerodynamic simulation is performed via the CFD solver.Figure 5
The FWGN excitation signals to drive pitch and plunge motions.Figure 6
The flight parameter samples to cover the range of concern.According to the time-dependent surrogate approach, 4000 Kriging models are required since the length of the FWGN excitation is 2000 for two unsteady aerodynamic outputs of the system, that is, the lift coefficientCl and the moment coefficient Cm. To demonstrate the accuracy of the training data, Cl and Cm are predicted via the Kriging models under two flight conditions in a transonic flight regime, one with M∞=0.76 and α0=1 deg and the other with M∞=0.82 and α0=4 deg. As shown in Figures 7 and 8, the results predicted via the Kriging models agree well with those obtained via the direct CFD technique.Figure 7
The training data validation of a NACA 64A010 airfoil atM∞=0.76 and α0=1deg.Figure 8
The training data validation of a NACA 64A010 airfoil atM∞=0.82 and α0=4deg.The CPU time to construct the Kriging models is shown in Table1, where all the computations are performed by using a laptop with two CPUs of 2.4 GHz. Even though the total CPU time cost associated with the construction of the Kriging models is a little bit expensive, the training data under an arbitrary flight condition in the range of concern can be obtained efficiently without performing any further CFD simulation. For example, the computational time of the direct CFD approach which solves Euler equations is about 0.284 h in order to obtain the training data at M∞=0.76 and α0=3deg, while the computational time of the Kriging models for the same case is 0.006 h only.Table 1
The computational cost associated with the time-dependent surrogate approach.
CPU time
Computation of the 20 training cases
4.98 h
Generation of the Kriging models
0.01 h
Total CPU time
4.99 h


## 3.2. Unsteady Aerodynamic Prediction

To demonstrate the application of the training data in unsteady aerodynamic prediction, the constructed data are used to identify the state-space models via the robust subspace algorithm. The unsteady aerodynamic outputs of the system are computed under different excitation signals, and the results are compared with those obtained by using the CFD training data.At first, the FWGN excitation signals, not the same as the training signal during the ROM construction, are used to drive the airfoil atM∞=0.76 and α0=3deg. The pitch and plunge motions of the airfoil are shown in Figure 9. The corresponding unsteady aerodynamic coefficients are shown in Figure 10. For further comparison, the flight condition with M∞=0.82 and α0=5deg is also selected. The excitation signals and the corresponding unsteady aerodynamic results are presented in Figures 11 and 12, respectively. In addition, a harmonic oscillation is chosen as the input signal at M∞=0.84 and α0=0deg. The unsteady aerodynamic results are presented in Figure 12.Figure 9
The pitch and plunge motions under the FWGN excitation for ROM validation atM∞=0.76 and α0=3deg.Figure 10
The ROM validation under the FWGN excitation atM∞=0.76 and α0=3deg (ROM1-training data via CFD; ROM2-training data via the proposed method).Figure 11
The FWGN excitation signals for ROM validation atM∞=0.82 and α0=5deg.Figure 12
The ROM validation under the FWGN excitation atM∞=0.82 and α0=5deg (ROM1-training data via CFD; ROM2-training data via the proposed method).As can be seen from Figures10, 12, and 13, the time histories of Cl and Cm obtained via the ROM2s identified by using the time-dependent surrogate approach agree well with those of ROM1s identified by using the direct CFD simulation data. That is, the proposed approach is able to generate good training data for ROMs in a transonic flight regime.Figure 13
The ROM validation under the harmonic excitation atM∞=0.84 and α0=0deg (ROM1-training data via CFD; ROM2-training data via the proposed method).


## 3.3. Aeroelastic Time-Marching Simulation

To further demonstrate the performance of the proposed approach in the flutter prediction in a transonic flight regime, the ROM2s are coupled with the structural model to make time-marching simulations. At first, the airfoil is forced to make three complete cycles of harmonic motion in pitch at the frequencyωα/2. Then, the aeroelastic system is allowed to evolve by its own self-induced loads.In the present case, the aeroelastic responses of the airfoil predicted via the ROM2s and ROM1s are compared under two flight conditions. As shown in Figure14, the aeroelastic responses of the airfoil with M∞=0.81 and α0=0deg are computed for V∞∗=0.73. Both ROM2 and ROM1 can predict a neutrally stable oscillation. Afterwards, the aeroelastic responses of the airfoil with M∞=0.85 and α0=1deg are evaluated. As shown in Figure 15, the results of ROM2 and ROM1 present a neutrally stable oscillation at V∞∗=0.57.Figure 14
The aeroelastic responses of the airfoil forV∞∗=0.73 at M∞=0.81 and α0=0deg (ROM1-training data via CFD; ROM2-training data via the proposed method).Figure 15
The aeroelastic responses of the airfoil forV∞∗=0.57 at M∞=0.85 and α0=1deg (ROM1-training data via CFD; ROM2-training data via the proposed method).Finally, Figures16 and 17 give the flutter boundaries of the aeroelastic model with the mean angle of attack α0=1deg and α0=3deg, respectively. As shown in these two figures, the computational results based on ROM2s and ROM1s get a good agreement in the flutter boundaries. Nevertheless, the computation of ROM2 is much more efficient than that of ROM1 since no further CFD simulations are required once the Kriging models are constructed.Figure 16
The flutter boundary of the airfoil with the mean angle of attackα0=1deg (ROM1-training data via CFD; ROM2-training data via the proposed method).Figure 17
The flutter boundary of the airfoil with the mean angle of attackα0=3deg (ROM1-training data via CFD; ROM2-training data via the proposed method).As can be seen from the flutter boundaries of the airfoil at different mean angle of attack, the aeroelastic characteristics of the model are very sensitive to the mean angle of attack in a transonic flight regime. Hence, the effects of the variations of the mean angle of attack have to be taken into account for aeroelastic analysis.


## 4. Conclusions

In this study, an efficient approach is proposed to generate the training data for constructing the reduced-order model of an unsteady aerodynamic system with the variation of flight parameters in a transonic flight regime. The Kriging technique is used in the approach to approximate the relationship between the unsteady aerodynamic outputs of the system and the flight parameters. No further CFD simulations are required to generate the training data once the Kriging models are constructed. The training data obtained via the proposed approach are used to identify the reduced-order aerodynamic models of a NACA 64A010 airfoil via a robust subspace identification algorithm. The unsteady aerodynamics and aeroelastic responses of the airfoil under various flight conditions in a transonic flight regime are computed. The results agree well with those obtained by using the direct simulations of computational fluid dynamics.When the time-dependent surrogate approach is applied to a three-dimensional aeroelastic model in a transonic flight regime, it is expected to save more computation time compared with two-dimensional cases. Furthermore, the unsteady aerodynamics generated via the approach can be used as the validation data to determine the validation region of an unsteady aerodynamic ROM. Based on the time-dependent surrogate approach, a method by assembling the aerodynamic outputs of a set of precomputed local linear models will be investigated to construct the aerodynamic ROM for a range of flight parameters.



---

*Source: 4083538-2018-09-12.xml*