Yuting DAI, Haoran RONG, You WU, Chao YANG, Yuntao XU
School of Aeronautic Science and Engineering, Beihang University, Beijing 100083, China
KEYWORDS Deep learning;Dynamic stall;Limit-cycle oscillation;Reduced order model;Stall flutter
Abstract The modeling of dynamic stall aerodynamics is essential to stall flutter, due to the flow separation in a large-amplitude pitching oscillation process.A newly neural network based Reduced Order Model (ROM) framework for predicting the aerodynamic forces of an airfoil undergoing large-amplitude pitching oscillation at various velocities is presented in this work. First, the dynamic stall aerodynamics is calculated by solving RANS equations and the transitional SST-γ model.Afterwards,the stall flutter bifurcation behavior is calculated by the above CFD solver coupled with structural dynamic equation. The critical flutter speed and limit-cycle oscillation amplitudes are consistent with those obtained by experiments. A newly multi-layer Gated Recurrent Unit (GRU) neural network based ROM is constructed to accelerate the calculation of aerodynamic forces. The training and validation process are carried out upon the unsteady aerodynamic data obtained by the proposed CFD method.The well-trained ROM is then coupled with the structure equation at a specific velocity,the Limit-Cycle Oscillation(LCO)of stall flutter under this flow condition is predicted precisely and more quickly.In order to predict both the critical flutter velocity and LCO amplitudes after bifurcation at different velocities,a new ROM with GRU neural network considering the variation of flow velocities is developed. The stall flutter results predicted by ROM agree well with the CFD ones at different velocities.Finally,a brief sensitivity analysis of two structural parameters of ROM is carried out. It infers the potential of the presented modeling method to depict the nonlinearity of dynamic stall and stall flutter phenomenon.
Stall flutter is a self-excited nonlinear aeroelastic phenomenon at large angle of attack.Unlike the classical flutter,stall flutter occurs with flow separation and reattachment, which causes strong nonlinearity on aerodynamic loads. Stall flutter problems exist in many fields. For example, the blades on helicopters are sometimes in harmonic large-amplitude motion to balance the aerodynamic loads, leading to the state of dynamic stall or even stall flutter. High-aspect-ratio wings are usually designed in a high-altitude long-endurance unmanned aerial vehicle to enhance the lift-to-drag ratio.This kind of aircraft may fly slowly with flexible structure and may vibrate with large deformation when encountering gust. Stall flutter is also a focus for a high-altitude long-endurance unmanned aerial vehicle. Stall flutter also exists in machines with rotating blades such as turbomachinery blades.
Due to the strong nonlinear aerodynamic feature, it is difficult and complicated to predict the stall flutter phenomenon exactly and efficiently. Researchers have implemented many experiments on stall flutter and developed numerous stall flutter simulation theories.Dimitriadis and Jing1conducted experiments on NACA0012 wing to study the stall flutter behavior,indicating that stall flutter develops with a subcritical Hopf bifurcation. Bhat and Govardhan2experimentally studied the stall flutter characteristics by sinusoidally oscillating the NACA0012 airfoil at low Reynolds number,and the stall flutter boundaries are distinguished with energy transfer theory.Fagley et al.3applied the cyber-physical systems to stall flutter experiments and explored the areas of the stall flutter parameter space. Culler and Farnsworth4studied the phenomenon of higher frequency with flow separation information by experiment method. Refs. 5-8 implemented a series of experiments to systematically study the LCO characteristics with low Reynolds number effect.
The experimental method often costs much and needs long preparation periods, and aerodynamic calculation methods were developed by researchers. Tang and Dowell9,10used the ONERA aerodynamic model to calculate the stall flutter of a high-aspect-ratio wing with a slender body, and the results were compared with those of experiments. Laxman and Venkatesan11modified the ONERA model,which was coupled with the structure model to analyze the pitch and plunge response and explain the existence of the helicopters vibratory frequencies below the blades rotor frequency. Wang et al.12built an aeroelastic model for two-dimensional airfoil with Leishman-Beddos model based aerodynamic calculation method and the uncertainty of aerodynamic and structure parameters were quantified. The CFD-based method became popular in recent years because of its accuracy. Yabili et al.13developed a Fluid-Structure Interaction (FSI) solver in Open-FOAM to simulate the dynamic stall and stall flutter characteristics of rigid airfoil, which were correlated with the experiments by Dimitriadis. A large amount of researchers14-21performed CFD calculations with different turbulence models on pitching airfoil at various Reynolds numbers, pitch amplitudes and reduced frequencies, in order to simulate the dynamic stall accurately. To sum up, dynamic stall and stall flutter characteristics can be simulated by CFD-based method accurately.
However, the CFD method may not be suitable for many engineering applications because it is time-consuming and sensitive to the parameters such as grids, time step sizes, turbulence models, numerical scheme, etc. In recent years, rapid progress has been achieved in ROM methods. The ROMbased unsteady aerodynamic model is a low-dimensional representation of the full-order aerodynamic system, such as the aerodynamic results obtained by CFD method or experiment.The ROM methods can predict the aerodynamic characteristics accurately and reduce the order of computation costs.Kou and Zhang22summarized the progress in data-driven ROM in recent years. Current ROM algorithms for modeling the unsteady aerodynamics can be classified as feature extraction and system identification.The main idea of feature extraction method is to obtain the main feature of the flow,known as coherent structures or modes, to understand the complex flow physics gradually. According to the difference of decomposition ways of flow modes, there are two methods of feature extraction, named as Proper Orthogonal Decomposition(POD)23and Dynamic Mode Decomposition (DMD),24,25where POD method extracts flow modes by energy content and DMD method obtains flow modes by single frequency and growth rate. System identification method is to treat the aerodynamic system as a black-box without any physical insight,or a grey-box with partial prior knowledge.The inputs of this model are the flow state parameters and structure motion, while the outputs are the aerodynamic loads on the body.The progress of system identification based ROM is usually completed as follows.First,the training data from the fullorder model is obtained.Then,the parameters of the model are optimized to minimize the loss function on the training data.Finally,the unsteady aerodynamic forces for arbitrary motions can be predicted. The unsteady aerodynamic characteristics vary with the degree of unsteadiness and nonlinearity. Specially, it can be distinguished as steady-linear, unsteadylinear and unsteady-nonlinear according to flow conditions and body motions.26The unsteady-linear models, including Autoregressive Moving Average (ARMA) method,27,28are able to simulate the unsteady-linear or weak-nonlinear aerodynamic characteristics, classical flutter response or transonic Limit-Cycle Oscillations (LCO) characteristics for instance.The numerous unsteady-nonlinear models, such as Volterra series,29,30Kriging,31support vector machine.32and Artificial Neural Network (ANN), are able to simulate the nonlinear aerodynamic characteristics, dynamic stall and stall flutter for instance.
Among all the unsteady nonlinear models, ANN-based ROM became prevalent on account of its ability of approximating the strong nonlinear relationship between inputs and outputs.33ANN simulates the way in which the brain processes information. The simple units, named as neurons, are assigned together with different weights to eventually form the entire ANN. An activation function is added between two layers of neurons to improve the nonlinear expression ability of the model. Numerous ANN-based ROMs are used to identify and predict the transonic unsteady aerodynamic characteristics. The categories of these ANN-based ROMs consist of Full Connected Neural Network (FCNN) by Marques and Anderson,34fuzzy neural network with output feedback by Winter and Breitsamter,35Radial Basis Function(RBF)neural networks by Wang36and Kou37et al.,and Long Short-Term Memory(LSTM)neural network by Kou et al.38,39and Halder et al.40.Then,by coupling the structure equations,the ROMs can be used to predict the transonic flutter characteristics of wing section with the variation of structure parameters.The studies on the nonlinear ROM for the dynamic stall and stall flutter are much less. Spentzos et al.41presented FCNN with output feedback based ROM to simulate the three-dimensional dynamic stall of helicopter rotor blades.Qiu and Wang42also developed FCNN with output feedback to predict the stall flutter of a two-dimensional airfoil, and investigated the relationship between the stall flutter characteristics and the structure natural frequency and flow vortex shedding frequency. Zhang et al.43proposed the Recurrent Neural
Fig. 1 Structure of GRU network.
Fig. 3 Aeroelastic model with freedom of pitch.
Fig.2 Structure of GRU recurrent neural network in this work.the structure equation.Moreover,a new ROM considering the variation of flow velocities is developed to predict the bifurcation point and LCO characteristics at different velocities. The bifurcation and LCO results predicted by ROM are compared with the CFD results.Finally,a brief sensitivity analysis of two structural parameters of ROM is carried out.
2.1.Reduced order modeling based on multi-layer GRU network
Recurrent neural network is suitable for time-series predicting problem,in which the sequential information of the time-series data can be preserved and transmitted in the sequence direction.In order to overcome the drawbacks of traditional RNNs such as the vanishing gradient, LSTM network was presented by Hochreiter and Schmidhuber.45The LSTM cell contains three gates, including input gate, output gate and forget gate.Cho et al.46presented GRU by simplifying the structure of LSTM. The structure of GRU cell is shown in Fig. 1.
Note that the symbol ‘‘×” represents Hadamard product.GRU cells employ-two gates, namely reset gate r and update gate z,to control the preservation of the information from previous cells.The equations to compute these states and gates are given as follows:Network (RNN) based ROM to identify and to predict the aerodynamic forces at given motions of dynamic stall.In other terms,Wang et al.44developed the RNN-based unsteady aerodynamic ROM and investigated the effect of introducing crossvalidation and bootstrapping in the training process. To sum up, the researches on the ANN-based ROM in terms of dynamic stall and stall flutter are quite absent, and the angle of attack of the existing researches is not high enough(usually no more than 20°), which does not capture the strong nonlinear effect induced by deep stall. However, the amplitude of stall flutter LCOs in experiments can reach 40°or even higher.8Thus, it is necessary to build the ROM that has the ability to capture the strong nonlinear relationship between aerodynamic forces and structure motions.
In this paper, the stall flutter bifurcation behavior of a pitching NACA0012 airfoil at low Reynolds number is calculated by CFD-based FSI method,and the results are validated with experiments. Then the multi-layer GRU neural network based ROM is presented to simulate the aerodynamics characteristics of an airfoil undergoing large-amplitude harmonic pitch motions(up to 40°),and the stall flutter LCO at a specific flow velocity is predicted by the proposed ROM coupled with
where W and U stand for weight matrices,the subscripts r and z represent the reset gates and update gates respectively, n is the current time step, x denotes the input vector of this cell,φ denotes the activation function that maps the variables into range between 0 and 1,I represents the unit matrix,and h and hˆare the cell hidden state and the candidate state,respectively.The reset gates r computed by Eq. (1) reveal the relevance of the previous state h(n-1)and the current candidate state hˆ(n),meaning that the current candidate state hˆ(n) is independent of the previous state h(n-1)if the components of the reset gate r(n)are close to zeros, as shown in Eq. (3). The update gates z computed by Eq. (2) decide whether to forget the previous memory. If the components of the update gates z(n)are close to zeros, the previous memory, represented by the previous state h(n-1), is preserved to current hidden state h(n).
Table 1 Structure parameters of NACA0012 wing section.
Fig. 4 Computational mesh around a NACA0012 airfoil.
The inputs of neural network and the outputs of each hidden layer are usually far from zeros if not treated, leading to the trouble of vanishing gradient. To address such issue, Ioffe and Szegedy47presented the method of Batch Normalization(BN), accelerating the training process of neural network and preventing it from vanishing gradient. The procedure of BN is shown as follows:
where xirepresents the component of input vector x,and μ and σ represent the average and variance of input vector x,respectively. ε is a small positive factor, preventing the denominator from being close to zeros. λ and β, trained by the Backward Propagation(BP)algorithm,are the factors of linear transformation that converts the input vector ×into the form with an average of β and a variance of λ.
Fig. 5 Computational mesh for validation case.
Fig. 6 Comparison of lift coefficient in current work with experimental and previous CFD results.
The structure of neural network in this work is shown in Fig. 2, consisting of 2 GRU hidden layers and BN layers,and the unit number of each GRU layer is 25. The last cell state is considered in the output of last GRU layer.The Adam optimizer is utilized to minimize the loss function in the training process for BP algorithm.In addition,in order to make full use of the training data set, the input matrix is transformed into a ‘‘ladder” form, shown as Eq. (6). ndand N denote the number of GRU cells, or time delays number, and training sample size,respectively.x(n)and y(n)represent the input vector and output vector at time step n, respectively.Fig.7 Comparison of bifurcation diagram of current work with experiments.
The typical two-dimensional airfoil concerning the freedom of pitch motion is shown in Fig.3.The aeroelastic equations can be written as.
Kα, Cαand ζ are the stiffness, damping and damping ratio of the pitch motion,respectively.Iαis the mass moment of inertia about stiffness center.M and Cmare the aerodynamic moment and moment coefficient about stiffness center, respectively. ρ and U∞are the density and freestream velocity of flow,respectively. c and l are the chord length and span of the wing section, respectively. Convert Eq. (7) into a state space form as Eq. (8), and it can be easily solved by a 4th-order Runge-Kutta method, given an initial condition.
The structure parameters of the aeroelastic model in this paper are given from the stall flutter experiments by Goyaniuk et al.,8as shown in Fig.3 and Table 1.The results of the experiments showed that the freedom of plunge can be negligible in a stall flutter. Thus, the NACA0012 wing section only has the freedom of pitch motion α. In this case, the natural frequency of the pitch motion is 2.76 Hz.
Fig. 8 Typical aeroelastic response at different flow velocities obtained by FSI method.
To obtain the nonlinear aerodynamic characteristic of the airfoil undergoing deep dynamic stall, the two-dimensional,incompressible, Reynolds-Averaged Navier-Stokes (RANS)equations are solved:where U-, P and μeffare the mean velocity, static pressure and effective viscosity of flow respectively, δijis the Kronecker delta,and t and ^xlrepresent the time and the coordinate component respectively. The Reynolds number in this paper is between 6.4 × 104and 9.6 × 104. The intermittency transition model (transitional SST-γ model48) is used to deal with the transition effect. The γ transition model is simplified based on transition SST model(γ-Reθtransition model),which solves only one transport equation for the turbulence intermittency γ.It avoids the need for the second Reθequation. The time step size of 1 × 10-3s is chosen. The mesh used in this paper is illustrated in Fig. 4. The total number of computational cells is around 3 × 104. The thickness of first grid layer is 1.5 × 10-5m and the maximum y + is 1. The sliding mesh technique is used, in which the internal zone moves rigidly in a pitch motion. In this situation, the computational mesh remains a good mesh quality.
Fig. 9 Pitch angle and moment coefficient response at flow velocity U∞= 8 m/s.
Fig. 10 Pressure contours and stream lines around the airfoil at flow velocity U∞= 8 m/s.
Fig. 11 Two kinds of training signals of pitch motion.
Fig. 12 Validation data set.
This study aims to validate the numerical solutions of dynamic stall. To verify the CFD approach, the CFD results of forced pitching motion are compared with the results of experiment49and other numerical simulation results from previous work.16,21The pitch motion is given by.
Fig. 13 Training process of ROM1 (record every 50 epochs).
while ω and k = ωc/(2U∞) represent the angular frequency and reduced frequency respectively.
The inflow velocity is 14 m/s and the Reynolds number is 1.4×105.The computational mesh for the forced motion case is recreated as shown in Fig.5 due to the change of elastic axis and flow Reynolds number. The thickness of first grid layer is 7.5 × 10-6m and the maximum y + of this mesh is 1. The total number of computational cells is around 3.7 × 104. The relative radiuses of the internal and external zones are the same as the mesh in Fig. 4. Other CFD solution setup is consistent with that in Section 2.3.The lift coefficient comparison results are shown in Fig.6.The tendency of aerodynamic loads in current work is generally consistent with that in previous studies,indicating that the CFD solution setup is reasonable to acquire the nonlinear unsteady aerodynamic characteristics of dynamic stall.
The stall flutter experiments by Goyaniuk et al.8were conducted for flow velocity U∞between 5.4 m/s and 12 m/s and different kinds of initial perturbation were applied to the airfoil.It is worth noting that in this experiment, under the same flow velocity condition, small initial perturbation may lead to small amplitude (about 5°) oscillations due to laminar separation, while large initial perturbation may lead to largeamplitude LCO due to stall flutter. In this paper, FSI method is implemented for flow velocity U∞between 6 m/s and 10 m/s,corresponding to the range of Reynolds number between 6.4 × 104and 9.6 × 104. The large initial perturbation of 3 rad/s angular velocity is applied to all flutter computation case to induce stall flutter. Fig. 7 shows the pitch amplitude with different flow velocities obtained by FSI method compared with experimental results. Fig. 8 shows the typical time-domain responses of aeroelastic system for different flow velocities. In current FSI results, the aeroelastic response converges until the flow velocity reaches 6.5 m/s. Then small amplitude oscillation about 4° occurs below the velocity of 7 m/s. Large-amplitude LCO occurs when the flow velocity comes up to 7.5 m/s,indicating that the bifurcation velocity of FSI method is between 7 m/s and 7.5 m/s, and the bifurcation velocity in the experiment is around 7.3 m/s. Then the amplitude of LCO gradually expands as the flow velocity increases.The tendency of aeroelastic phenomena of current work is generally consistent with experiments.
Fig. 14 Identification results.
The LCO and moment coefficient response at flow velocity of 8 m/s are shown in Fig. 9 and the pressure contours and stream lines around the airfoil at six time instances are illustrated in Fig. 10 to reveal the stall flutter process. At point A and point B, in the upstroke process, a strong leadingedge vortex appears, causing the upward moment. In the downstroke process, the vortex breaks down and moves to the trailing edge at point C, resulting in the downward moment.Then the vortex is weakened at point D,which causes a drastic reduction of the downward moment. At point E, a trailing-edge vortex appears, leading to the appearance of the downward moment again.Finally at point F,upward moment arises because of the departure of the vortex.The aerodynamic patterns of the leading-edge vortex and trailing-edge vortex provide energy to the airfoil, forming the condition of stall flutter.
In this section, the neural network based ROM, named as ROM1, is established to predict the stall flutter response at the flow velocity of 8 m/s. It is worth noting that the design of the training signal can be directed by the characteristic of stall flutter response.First,the structural response of stall flutter is usually large-amplitude sinusoidal oscillations.Then,the frequency of stall flutter response is generally around the structure eigen frequency. Thus, two kinds of training signals are designated for consideration. The frequency range is between 2.5 Hz and 3.5 Hz, corresponding to the reduced frequency range from 0.153 to 0.214.The sampling time of the first signal(signal A) is 20 s. It is in the form of a swept-frequency as shown in Eq. (11) and Fig. 11(a).
Fig. 15 Dynamic stall prediction with different amplitudes and reduced frequencies.
Table 2 Identification and prediction errors of two kinds of training signals.
Fig. 16 Generalization ability evaluation of amplitudes and reduced frequencies.
Fig. 17 Flutter prediction by ROM1 at U∞= 8 m/s.
Fig. 18 Moment coefficient of same sinusoidal motion at different flow velocities computed by CFD.
Fig. 19 Training signal of ROM2.
Fig. 20 Training process of ROM2 (record every 50 epochs).
The sampling time of the second signal(signal B)is 10 s.It is in the form of a sinusoidal superposition as shown in Eq.(12) and Fig. 11(b).
where A is the random large amplitude and the frequency is evenly distributed within the above range.
The moment coefficient result of given training signals is acquired by CFD method, which are treated as the training data sets.The time step size is set to 5×10-3s so that the size of training data sets N is 4000, where the training data set of the signal B is replicated to ensure the same sample size. Then another signal is used as the validation data set in Fig. 12,which is intended to observe whether the neural network is overfitting in the training process and is not involved in the update of weights. The structure of the neural network in this work is shown in Section 2.1. The unit number of each GRU layer,or the time delays number ndis 25.The Adam optimizer is utilized to minimize the loss function in the training process for BP algorithm. According to Eq. (6), the input vector and the output vector are x(n)= [α(n)] and y(n)= [Cm(n)], respectively,developing a Single Input Single Output(SISO)system.The Mean Square Error (MSE) is used as the loss function.The training processes of two training signals and the identification results are shown in Fig. 13 and Fig. 14. At 8000 epochs, the MSE of both signals almost stop decreasing and the training processes are concluded.
Then sequences of dynamic stall results with the same amplitude and different reduced frequencies are used to test the generalization ability on reduced frequency,and the results with the same reduced frequency and different amplitudes are used to test the generalization ability on amplitude. In these samples, four cases with extreme amplitudes and reduced frequencies are used to test the generalization ability of the ROM trained by both signals. Every sinusoidal motion case is computed for 2 s, where the last 1 s CFD result is treated as the test data set. The moment coefficient prediction results compared with CFD results for motions of different amplitudes and reduced frequencies are shown in Fig.15.Identification error and prediction error e of ROM are calculated as.
where yCFDand yROMrefer to the output matrix acquired by CFD method and ROM, respectively. The error results of identification and prediction are summarized in Table 2. It is obvious that although strong nonlinear relations exist between pitch angle and moment coefficient when pitch amplitude is large, GRU-based RNN is able to identify and predict the nonlinear aerodynamic force approximately with good accuracy.
The generalization ability evaluation of amplitudes and reduced frequencies are shown in Fig. 16. Apparently, the error curve of signal A is much smoother, indicating that the generalization ability of ROM trained by signal A is stronger,since the amplitude and frequency characteristics in signal A are continuously varying.Meanwhile,in signal B,small amplitude messages are insufficient and the characteristics are discrete. However, when the amplitudes and reduced frequencies of oscillation are in the training set, the prediction errors of ROM trained by signal B are much lower. That is because the full characteristics of amplitudes and frequencies are contained in each segment of the training sample of signal B, and partial frequency and amplitude information are included in each sample segment of signal A. It can be expected that the prediction effect of LCO for this specific frequency and amplitude range is better with the ROM trained with signal B. Moreover, the CFD simulation cost of signal B is only half that of signal A.Therefore,the signal B is finally selected as the training signal for the ROM1, and the subsequent study is carried out.
Fig. 21 Identification result of ROM2.
Fig. 22 Dynamic stall prediction with same amplitude and frequency at different velocities.
Table 3 Identification and prediction errors of two kinds of ROMs.
Then ROM1 trained by signal B is used to predict the flutter characteristic at flow velocity of 8 m/s. The initial condition, initial angular velocity of 3 rad/s, is the same as that in FSI case. The aeroelastic Eq. (8) is solved by Runge-Kutta method and the moment coefficient is predicted by ROM1.The LCO responses predicted by ROM1 are compared with FSI results in Fig.17.The relative error of the LCO amplitude predicted by ROM1 is within 5 %. It shows consistency in terms of the amplitude and frequency of LCO in general,indicating the ability of the multi-layer GRU neural network based ROM to predict the strong nonlinear aerodynamic force and aeroelastic response.
Fig. 23 Stall flutter prediction at various flow velocities.
Fig. 24 Comparison of bifurcation diagram of two ROMs with FSI.
In this section,a new ROM is established to predict the bifurcation point and LCO amplitudes at different velocities. The variation in flow velocity or Reynolds number can cause great changes in aerodynamic properties under the same pitch motion. Fig. 18 shows the results of the moment coefficient of the same sinusoidal motion at different flow velocities calculated by CFD. Apart from some minor variations, the main influence on the aerodynamic force lies in the phase difference with the pitch motion, which results in the variation of the aerodynamic work done per cycle at different velocities. Thus,the ROM considering the input of velocity, named as ROM2,is established in this section.
According to Eq.(6),the input vector and the output vector of ROM2 are x(n)= [α(n), U∞] and y(n)= [Cm(n)] respectively,developing a Multiple Input Single Output (MISO) system.Other neural network structure and training setup are the same as ROM1.The pitch training signal function of ROM2 is also in the form of Eq.(12),the signal B,and flow velocity vector of U∞= [6, 7, 8, 9](m/s) is included, as shown in Fig. 19. The moment coefficient results under each velocity condition by CFD method are collected as the training data set, so that the size of training data set N is 8000.It is shown in Fig.7 that small amplitude oscillations occur at low flow velocities.However,there is little small amplitude message in the training data set. But this is no matter because the bifurcation point and LCO amplitude for different velocities, rather than the smallamplitude oscillation characteristics at low velocities, are concerned in this section. The training process is displayed in Fig. 20 and the identification result of ROM2 is shown in Fig. 21. The moment coefficient results of dynamic stall for various velocities with the same sinusoidal motion (amplitude 0.698 rad, frequency 3.3 Hz) by CFD method are collected as the test data set.To further illustrate the necessity of considering the velocity input in ROM2, the prediction results of ROM2 and ROM1 which is established in Section 3.3 and ignores the velocity input, are compared in Fig. 22. The identification error and prediction error of both ROM are summarized in Table 3. It can be indicated that ROM2 is capable of predicting the nonlinear aerodynamic forces in the range of velocities of the training set, but ROM1 does not work.
Fig. 25 Training processes of ROMs with different structural parameters (record every 20 epochs).
Table 4 Identification and prediction errors of ROMs with different structural parameters.
Then ROM1 and ROM2 are used to predict the aeroelastic responses at different velocities as shown in Fig.23.The bifurcation diagram predicted by two ROMs are compared with the FSI results in Fig. 24. As mentioned previously, ROM2 is not enough to capture the small amplitude oscillation characteristics, but the bifurcation velocity and stable large-amplitude LCO characteristics after the bifurcation point are precisely predicted. ROM1 predicts the large-amplitude LCO well, but it is not able to well simulate the aeroelastic response around the bifurcation point, although the training velocity of ROM1(flow velocity of 8 m/s)is close to the bifurcation point(between 7 m/s and 7.5 m/s).This result emphasizes once again that the aerodynamic characteristics are sensitive to the flow velocity or Reynolds number. It is essential to consider the velocity or Reynolds number as one of the input if the bifurcation velocity is unknown.
In this section, a brief sensitivity analysis of two structural parameters,namely the number of time delays ndand the number of GRU layers, is carried out. The structure of the neural network established in previous sections, containing 25 orders of time delay and 2 GRU layers, is considered as baseline.Then ROMs with time delays of 10 and 40 are built to study the effect of the number of time delays nd. ROMs with GRU layers of 1 and 3 are trained to research the influence of the number of GRU layers, and the same increase or decrease is done for the BN layers. The sinusoidal superposition signal in Section 3.3, as shown in Fig. 11 (b) and Eq. (12), and the corresponding aerodynamic moment by CFD in Fig. 14 (b)are selected as the training data set. The CFD results of dynamic stall at different amplitudes and reduced frequencies in Fig. 15 are chosen as the test data set in this section. For each ROM structural parameter, 5 models are trained to reduce the error caused by training fluctuations, and the ROM output is the average of the predictions of 5 models.The typical training processes of ROMs with the above structural parameters are displayed in Fig. 25. The identification and prediction errors of each ROM are shown in Table 4.Fig.26 and Fig.27 exhibit the generalization ability evaluation of each ROM.
Fig. 26 Generalization ability evaluation of ROMs with 2 GRU layers and different time delays number nd.
Fig. 27 Generalization ability evaluation of ROMs with nd = 25 and different GRU layers.
Compared with the baseline ROM, when the number of time delays is reduced to 10 or the GRU layers are decreased,the MSE decreases more slowly in the training process as can be seen in Fig. 25, resulting in the increase of identification errors, and the prediction errors of dynamic stall at different reduced frequencies also increase. In other words, reducing the complexity of the neural network may lead to underfitting.When the number of GRU layers increases to 3,the identification and prediction performance does not differ much from the baseline ROM. However, when the number of time delays is increased to 40,the prediction error of dynamic stall increases slightly compared to the baseline ROM, but the difference in identification performance is not significant, meaning a slight overfitting. In addition, the training processes cost more if the complexity of the network increases.Meanwhile,it is interesting to note that when the number of time delays or GRU layers decrease, the prediction errors of small amplitude dynamic stall reduce. Reducing complexity of neural network may allow for greater generalization ability to linear part.
Herein, a CFD-based FSI solver is developed to calculate the stall flutter behavior of a NACA0012 airfoil.A ROM based on a multi-layer GRU network is introduced to identify and predict the aerodynamic force of an airfoil undergoing dynamic stall. CFD solver is adopted to generate training data set and test the generalization ability of the proposed ROM.The well-trained model is then coupled with the structural equation to calculate the stall flutter response, and the results are compared with the stall flutter response obtained by the proposed FSI solver for validation. The following conclusions are drawn:
(1) Through the comparison with the experiment, it is seen that the CFD solver using the transitional model could accurately calculate the aerodynamic characteristic of an airfoil undergoing dynamic stall at a low Reynolds number. After coupled with the structural equation,the CFD-based FSI solver could calculate the critical flutter velocity and LCO amplitude after bifurcation.This ensures that the training data generated by the CFD solver and the validation data generated by the FSI solver are more convincible.
(2) The design of training signal has a significant impact on the prediction accuracy of the ROM.A properly designed training signal can improve the ability of the ROM to predict stall flutter in a specific amplitude and frequency range while minimizing the computational cost.
(3) The order of time delay and the number of layers significantly influence the performance of the proposed GRU network.A less layer or time delay may lead to underfitting that is characterized by a high training and predicting error.However,using more layers or time delay may cause overfitting problem.This issue results in the reduction of the ROM’s prediction accuracy.
Declaration of Competing Interest
The authors declare that they have no known competing financial interests or personal relationships that could have appeared to influence the work reported in this paper.
Acknowledgement
This study is supported by the National Natural Science Foundation of China (No. 11672018).
