Unsteady aerodynamic identification based on recurrent neural networks

Bo Zhang1 , Jinglong Han2 , Tuanyuan Zhang3 , Ruiqun Ma4

1, 2, 4State Key Laboratory of Mechanics and Control of Mechanical Structures, Nanjing University of Aeronautics and Astronautics, Nanjing, 210016, China

3China Academy of Aerospace Aerodynamics, Beijing, 100074, China

2Corresponding author

Journal of Vibroengineering, Vol. 23, Issue 2, 2021, p. 449-458. https://doi.org/10.21595/jve.2020.21612
Received 23 July 2020; received in revised form 12 October 2020; accepted 12 November 2020; published 15 December 2020

Copyright © 2020 Bo Zhang, et al. This is an open access article distributed under the Creative Commons Attribution License, which permits unrestricted use, distribution, and reproduction in any medium, provided the original work is properly cited.
Creative Commons License
Table of Contents Download PDF Acknowledgements References
Cite this article
Views 81
Reads 33
Downloads 378
WoS Core Citations 0
CrossRef Citations 0
Abstract.

The dynamic stall at high angle of attack is an important aerodynamic problem faced by aircraft, and it has always been a hotspot of aerodynamic research. The traditional reduced order model (ROM) methods needs to establish an accurate model, and has a high demand for experience. In this paper, a novel nonlinear aerodynamic identification method based on recurrent neural networks (RNNs) is proposed. The computational fluid dynamics (CFD) method is used to calculate the unsteady aerodynamic parameters of the NACA0012 airfoil. A group of sinusoidal chirp signals with variable amplitude and frequency are adopted as the excitation signals, and the obtained data are used to train the recurrent neural networks, and the ROM of the nonlinear aerodynamic model of high angle of attack dynamic stall is obtained. Finally, the aerodynamic parameters of a group of composite sinusoidal motion signals different from the training signals are predicted by the trained neural networks model and compared with the CFD results.

Keywords: unsteady aerodynamic identification, recurrent neural networks, computational fluid dynamics, high angle of attack, reduced order model.

1. Introduction

Dynamic stall refers to the unsteady separation of the flow field around an oscillating or other unsteady airfoils when it exceeds the critical angle of attack, which leads to stall. For example, a dynamic stall occurs in the flow of a fighter’s post-stall maneuver, helicopter blade rotation, aeroengine’s turbine and fan blade flutter. Its main feature is the existence of complex unsteady separation vortex and large-scale structure vortex in the flow field, the vortex turbulence separates from the leading edge or trailing edge area of the airfoil, and the airfoil keeps stalling until the air flow re-adheres to the airfoil. Compared with the static stall, the dynamic stall is more complex. The lift, moment and drag characteristics of the airfoil vary with the different motion forms of the airfoil, and there will be obvious nonlinear hysteresis phenomenon at the maximum angle of attack when the dynamic stall occurs. Once dynamic stall occurs, the consequences are more serious and the duration is longer [1]. The aerodynamic coefficient may greatly deviate from the static value. So the simple static stall research can’t meet the needs of dynamic stall research.

At present, the research of dynamic stall mainly adopts numerical simulation or wind tunnel test methods. Wind tunnel test can measure the aerodynamic force and pressure of the reduced ratio airfoil in the wind tunnel, which can intuitively observe the hysteresis effect of aerodynamic force in the process of dynamic stall of airfoil. However, the wind tunnel test method is difficult to capture the generation, separation and reattachment of vortices in the oscillating process of airfoil, its preliminary preparation and specific wind tunnel test process are very complex, expensive and time-consuming, which can only be used for principle research of specific working conditions. The numerical simulation method can overcome the shortcomings and difficulties of the wind tunnel test [2]. It can simulate the dynamic stall of airfoil under various conditions and capture the whole falling off process of the detached vortices. This method has become an important way to study the unsteady dynamic stall characteristics.

Numerical methods can be divided into semi empirical model based on test data and computational fluid dynamics (CFD) method. Researchers have developed a series of semi empirical dynamic stall theoretical models, such as Gangwani model [3], ONERA model [4], Leishman-Beddoes (L-B) model [5]. The Gangwani model is the early semi empirical model. The ONERA model and Leishman Beddoes (L-B) model are internationally recognized as two-dimensional aerodynamic force and widely used. Although the above two methods have achieved the engineering application accuracy and greatly reduced the calculation time, the semi empirical model is not a strict prediction method, which is only applicable to specific airfoil, and there is a big deviation between the calculated value and the test value of the peak value in the whole stall process and the aerodynamic force during the reattachment after the flow separation. In addition, the semi empirical model can only calculate the aerodynamic force effectively, but it can't simulate the details of the vortex flow in the complex separation and reattachment phenomenon in the unsteady flow field of the airfoil. It’s difficult to carry out further research for semi empirical model due to the above reasons. The computational fluid dynamics method can capture the detail of the formation, development, shedding and reattachment process of the vortex in the dynamic stall process, which is an important means to study the dynamic stall problem. Dubuc [6] used the implicit dual-time step method to solve the Euler equation to study the unsteady motion of a NACA0012 airfoil with small amplitude transonic oscillation and accurately simulate the hysteresis effect of the airfoil in the shallow stall state. Later, a grid deformation algorithm based on transfinite interpolation method was developed and verified by the forced oscillation numerical simulation of flap validity of the method [7].

The CFD method has a high calculation accuracy, but it needs a lot of time and computing resources. Therefore, the reduced order model (ROM) methods for complex nonlinear systems have been widely concerned, such as proper orthogonal decomposition (POD) [8], Volterra series [9], dynamic mode decomposition (DMD) [10, 11], etc. However, the shortcomings of these methods in nonlinear and multi-scale simulation limit their application. Lucia [12, 13] developed a kind of unsteady aerodynamic reduced order model based on CFD technology, which was successfully used to study the flutter and gust response dynamics of aircraft at transonic speed or high angle of attack. However, most of these ROM methods are dynamic linear models, which could take into account the static nonlinear characteristics of the flow at transonic speed or high angle of attack, but cannot be used for the aerodynamic modeling research with large amplitude shock motion, unsteady separation and other dynamic nonlinear characteristics.

The traditional ROM methods has a high demand for experience, which is difficult to give a suitable explicit expression due to the complexity of the nonlinear flow field. The neural networks can deal with high-dimensional problems and does not need to give the explicit mathematical expression between the input and output of the identification system, which is an ideal tool for nonlinear aerodynamic prediction. Marques [14] identified the unsteady aerodynamic forces of transonic airfoils using an artificial neural networks with multi-layer functions. Zhang [15, 16] proposed a nonlinear unsteady aerodynamic layered reduced order model based on radial basis function networks (RBFNN), and used the model to predict the unsteady aerodynamic force, limit cycle oscillation and flutter characteristics of AGARD wing 115.6 and NACA 64A010 airfoils. Winter [17] proposed a nonlinear aerodynamic reduction scheme based on the series connection of a recurrent local linear neural fuzzy model and multilayer perceptron neural networks. Wang [18] uses a long-short term memory (LSTM) neural networks to construct a set of hypersurfaces representing the reduced fluid dynamic system. Afshar [19] proposes a flow field prediction model based on convolutional neural networks (CNNs), which is successful to predict the velocity and pressure field in geometry given the pixelated shape of the object and unseen flow conditions. Mannarino [20] develops a reduced order method based on continuous time recurrent neural networks (CTRNN) identification. Their researches show that deep neural networks identification is efficient for dealing with fluid dynamic problems with strong nonlinearity.

In this paper, a novel reduced order model method based on recurrent neural networks is presented to predict the nonlinear unsteady aerodynamic forces at high angle of attack. The unsteady aerodynamic force of NACA0012 airfoil at deep stall is calculated by the CFD method. The results are compared with the experimental data, and the effectiveness of the CFD model is verified. Two kinds of unsteady aerodynamic forces are calculated by CFD method: a variable frequency and amplitude motion is used as the identification training input signals; a composite sine signal is used to verify the effectiveness of the identification system. The trained recurrent neural networks is used for unsteady aerodynamic force identification, and the results predicted by neural networks are compared with the results of CFD.

2. Recurrent neural networks

2.1. Networks structure

Artificial neural networks can be divided into forward neural networks (FNNs) and recurrent neural networks (RNNs) according to whether there is feedback nodes. The architecture of basic forward neural networks is shown in Fig. 1.

Fig. 1. Architecture of forward neural networks

 Architecture of forward neural networks

There is not delay layer and feedback nodes in the FNNs, which cannot deal with the time series problems and dynamic problems. Its mapping relationship between input and output can be expressed as:

(1)
y = f x = f 2 W 2 f 1 W 1 x + b 1 + b 2 ,

where y is the output of the system; f represents the transfer function of the system; x is the input of the system; f1 and f2 denote the transfer functions of the neuron, respectively; W1 and W2 represent the weight of the neuron, respectively; and b1 and b2 are the biase of the neuron, respectively.

The structure of recurrent neural networks is shown in the Fig. 2.

Its mapping relationship between input and output can be expressed as:

(2)
y t = f y t - 1 , y t - 2 , , y t - p , x t - 1 , x t - 2 , , x t - q ,

where yt denotes the output of the system at t moment; and xt represents the input of the system at t moment.

2.2. Optimization algorithm

The commonly used optimization algorithms of neural networks include: gradient descent method; Newton method; Levenberg Marquardt (L-M) method; conjugate gradient method. Gradient descent method decreases rapidly in the first few steps, but with the approach to the optimal value, the gradient tends to zero, which makes the objective function decline slowly; Newton method can produce an ideal search direction near the optimal value; which could converge fast; L-M method combines the advantages of Newton method and gradient descent method, which could get the best solution quickly; conjugate gradient method needs less storage, but its convergence speed is slower than L-M method and Quasi Newton method. This paper adopts the L-M method as the optimization algorithm.

Fig. 2. Architecture of recurrent neural networks

 Architecture of recurrent neural networks

The L-M method is an iterative technique, which can find the minimum value of multivariate function expressed as the sum-of-squares of nonlinear real valued functions [21]:

(3)
E x = i = 1 m F i 2 x ,

where Ex is the minimum value of multivariate function, i is the training pattern index and Fix is the error for pattern i, which can be defined by the following formula:

(4)
F x = y x , t 1 - φ t 1 y x , t 2 - φ t 2 y x , t m - φ t m ,

where yx,ti and φti are scalar functions.

The weights of RNNs are updated by the following solutions of linear equations:

(5)
d k = - J x k T J x k + λ k I - 1 J x k T F x k ,

where Jxk is the Jacobian matric of Fxk; λk is a scalar regulation which controls the modulus sum direction of dk and I is the identity matrix. When λk is zero, the direction dk is the same as the Newton method. As λk tends to infinity, dk tends to the steepest descent direction and the mode tends to zero, which means that when λk is large enough, the term Fxk+dk<Fxk holds. Therefore, the term λk can be controlled to ensure the descent even if the second-order term which will limit the efficiency of the Newton method is encountered. Hence, the search direction of L-M method combines the advantages of Newton direction and steepest descent direction.

3. CFD simulation of aerodynamic at high angle of attack

3.1. Grid model

The spring grid method and transfinite interpolation method, which are widely used at present, can effectively avoid the appearance of negative volume and grid distortion, and have certain accuracy for small amplitude oscillation motion. However, it is difficult to ensure the grid orthogonality and smoothness when the airfoil with large deformation and large amplitude oscillation motion, which will result in grid distortion and negative volume and cannot be calculated. The nested grid allows the grid in different regions to perform independent parallel computing without changing the grid shape, and coupling by interpolation. Therefore, nested grid technology is often used to study the dynamic stall of high angle of attack. This work adopt the nested grid.

In this paper, the ICEM is used to generate the computational grid. Firstly, the body fitted grid around the airfoil is generated, which is a C-H type foreground structure grid with the boundary condition set to overset. Then, the fixed H type background structure grid is generated. The two types of grid are merged by interpolation. The distance between the boundary of the background grid and the airfoil is 15 times of the chord length. The thickness of the first layer grid on the airfoil surface is 0.001 mm, and the growth rate of grid thickness is 1.1. The total number of elements is 365482. Fig. 3 shows the structure of the nested grid and the boundary condition setting. From the area around the airfoil surface, it can be seen that the grid has better performance of orthogonality and smoothness, which provide a good basis for the Reynolds averaged Navier-Stokes (N-S) equation to capture the viscous and high angle of attack flow separation. The user-defined functions (UDF) is used to define the region motion of the foreground grid in the nested grid, which is used to control the motion of the airfoil.

Fig. 3. Nested mesh and local enlargement diagram of NACA0012 airfoil

 Nested mesh and local enlargement diagram of NACA0012 airfoil

3.2. Numerical simulation of unsteady dynamic stall of airfoil

In this section, the CFD model of NACA0012 airfoil is established, nested grid is employed to divide flow field, Transition SST kω model is chosen as the viscous model. According to the experiments in reference [22], the airfoil makes large amplitude oscillation around 1/4 chord length point. The move equation is shown as bellow:

(6)
α = α 0 + α 1 sin k τ ,

where α0 is the average angle of attack; α1 is the amplitude angle of attack; k is the reduction frequency; and τ is dimensionless time (reduction time). k and τ can be defined as follows [23]:

(7)
k = π f c u ,
(8)
τ = 2 u t c ,

where u is the flow velocity in the far field; t is the time; f is the oscillation frequency; and c is the reference chord length of the airfoil. The Reynolds number based on the chord length is defined by the following formula [23]:

(9)
R e = ρ u c μ ,

where ρ is the density of inflow air; and μ is dynamic viscosity.

3.3. Unsteady aerodynamic verification

The above method is used to simulate the unsteady aerodynamic force of NACA0012 low speed large amplitude oscillation. The calculation conditions are shown in Table 1. The test values of this state are given in reference [20]. Other calculation conditions are: solver type is density-based; pressure far-field boundary condition is adopted; the boundary of the overlapping region is set to overset; and the wall condition is slip free boundary. Firstly, the steady flow field of airfoil at the average angle of attack is calculated; then, set the calculated value of the steady flow field as the initial value of the unsteady calculation to ensure the accuracy of the transient calculation results. The calculation results are shown in Fig. 4.

Fig. 4. CFD results compared with experiments: a) lift coefficient, b) moment coefficient

 CFD results compared with experiments: a) lift coefficient, b) moment coefficient

a)

 CFD results compared with experiments: a) lift coefficient, b) moment coefficient

b)

Fig. 4 shows the variation law of lift coefficient and moment coefficient in the process of airfoil motion. It can be seen from the figure that both lift coefficient and moment coefficient are basically the same as the experiment values [20]. At the same time, the results of CFD can reflect the secondary stall which is not reflected in the test results. The results indicate that the calculation grid and numerical method adopted in this paper can prediction the nonlinear unsteady aerodynamic of airfoil at large angle of attack.

Table 1. Nonlinear aeroelastic system parameters

Parameters
Definition
Value
u
Far field velocity
14 m/s
c
Reference chord
1
k
Reduction frequency
0.1
α 0
Average angle of attack
10°
α 1
Amplitude angle of attack
15°
R e
Reynolds number
1.35×105
ρ
Density of air
1.226 kg/m3
T
Temperature
288.15 K
t
Time step
0.001 s

4. Results and discussion

4.1. Training data preparation

A group of sinusoidal chirp signals with variable amplitude and frequency are used as input signals for neural networks training. The maximum amplitude of pitch angle is 15°, and the frequency range is 0-3 Hz. The variation law of pitch angle is as follows:

(10)
α 1 = 0.0222 t 2 - 0.667 t + 15 sin 0.314 t 2 .

This signal is used to establish the motion input of the airfoil, and the unsteady aerodynamic force of the airfoil under this motion is calculated by the CFD method as the training parameter of the neural networks. The trained RNNs is used to predict the lift and moment coefficient of the airfoil. The identification results of neural networks are shown in Fig. 5 and Fig. 6.

Fig. 5. Lift coefficient identification results of chirp input

 Lift coefficient identification results of chirp input

The dotted line shows the response of computational fluid dynamics simulation, and the solid line shows the results of recurrent neural networks prediction. The two results are agreement very good, and it is hard to distinguish one from the other, which means that the trained networks could be used to predict the force of the unsteady aerodynamic.

4.2. Model verification

In order to verify the effectiveness of the identification model, a new set of composite sinusoidal motion signals, which are different from the training signals, are used as the verification input signals:

(11)
α 2 = 5 sin 2 π t + 3 sin 3 π t + 2 sin 5 π t .

Fig. 6. Moment coefficient identification results of chirp input

 Moment coefficient identification results of chirp input

Its maximum amplitude is 10°. Then, the trained neural networks model is used to predict the aerodynamic force of the airfoil in this motion. The predicted results are compared with the CFD results. The results are shown in Fig. 7 and Fig. 8.

Fig. 7. Lift coefficient prediction results of verification signal

 Lift coefficient prediction  results of verification signal

Fig. 8. Moment coefficient prediction results of verification signal

 Moment coefficient prediction  results of verification signal

The dotted line represents the CFD simulation results, and the solid line represents the RNNs prediction results. It can be seen from the figures that the two results are consistent, which indicates that the neural networks identification method can approach the nonlinear aerodynamic system very well.

The mean squared error (MSE) of the prediction data is used to evaluate the performance of the RNNs, which is defined as:

(12)
M S E = i = 1 N y ^ i - y i 2 N ,

where y^i is the prediction results of the RNNs; yi represent the CFD results; and N is the number of calculated dates. The MSE of the lift coefficient prediction is 4.0217×10-7, and the MSE of the lift coefficient prediction is 7.5282×10-10. The results show that RNNs can accurately predict the lift coefficient and moment coefficient, but the prediction of moment coefficient is more accurate.

5. Conclusions

In this work, the aerodynamic parameters of NACA0012 airfoil with sinusoidal motion around its 1/4 chord length are calculated by nested grid method, and compared with the test results, which verify the effectiveness of CFD method. A recurrent neural networks with time-delay element is used to approach the nonlinear aerodynamic system of large amplitude oscillation airfoil. A group of aerodynamic force of chirp motion with variable amplitude and frequency are obtained from the CFD method. The parameters of motion and aerodynamic force are used as input and output signals of the neural networks training. The trained networks is used to predict the aerodynamic force of composite sinusoidal motion. The results show that the neural networks ROM can accurately identify the nonlinear unsteady aerodynamic force of airfoil with large amplitude motion.

Acknowledgements

This work was supported by National Natural Science Foundation of China (Grant No. 11472133).

References

  1. Ericsson L. E., Reding J. P. Fluid mechanics of dynamic stall part I. Unsteady flow concepts. Journal of Fluids and Structures, Vol. 2, Issue 1, 1988, p. 1-33. [Publisher]
  2. Lu K., Xie Y. H., Zhang D. Numerical study of large amplitude, nonsinusoidal motion and camber effects on pitching airfoil propulsion. Journal of Fluids and Structures, Vol. 36, 2013, p. 184-194. [Publisher]
  3. Gangwani S. T. Prediction of dynamic stall and unsteady airloads for rotor blades. Journal of the American Helicopter Society, Vol. 27, Issue 4, 1982, p. 57-64. [Publisher]
  4. Peters D. A. Toward a unified lift model for use in rotor blade stability analyses. Journal of the American Helicopter Society, Vol. 30, Issue 3, 1985, p. 32-42. [Publisher]
  5. Leishman J. G., Beddoes T. S. A semi-empirical model for dynamic stall. Journal of the American Helicopter Society, Vol. 34, Issue 3, 1989, p. 3-17. [Search CrossRef]
  6. Dubuc L., Cantariti F., Woodgate M., Gribben B., Badcock K. J., Richards B. E. Solution of the unsteady euler equations using an implicit dual-time method. AIAA Journal, Vol. 36, Issue 8, 1998, p. 1417-1424. [Publisher]
  7. Dubuc L., Cantariti F., Woodgate M., Gribben B., Badcock K. J., Richards B. E. A Grid deformation technique for unsteady flow computations. International Journal for Numerical Methods in Fluids, Vol. 32, Issue 3, 2000, p. 285-311. [Publisher]
  8. Hall K. C., Thomas J. P., Dowell E. H. Proper orthogonal decomposition technique for transonic unsteady aerodynamic flows. AIAA Journal, Vol. 38, Issue 10, 2000, p. 1853-1862. [Publisher]
  9. Raveh D. E. Reduced-order models for nonlinear unsteady aerodynamics. AIAA Journal, Vol. 39, Issue 8, 2001, p. 1417-1429. [Publisher]
  10. Schmid P. J. Dynamic mode decomposition of numerical and experimental data. Journal of Fluid Mechanics, Vol. 656, Issue 10, 2008, p. 5-28. [Search CrossRef]
  11. Schmid P. J., Li L., Juniper M. P., Pust O. Applications of the dynamic mode decomposition. Theoretical and Computational Fluid Dynamics, Vol. 25, 2011, p. 249-259. [Publisher]
  12. Dowell E. H., Hall K. C. Modeling of fluid-structure interaction. Annual Review of Fluid Mechanics, Vol. 33, 2001, p. 445-490. [Publisher]
  13. Lucia D. J., Beran P. S., Silva W. A. Reduced-order modeling: new approaches for computational physics. Progress in Aerospace Sciences, Vol. 40, 2004, p. 51-117. [Publisher]
  14. Marques F. D., Anderson J. Identification and prediction of unsteady transonic aerodynamic loads by multi-layer functionals. Journal of Fluids and Structures, Vol. 15, Issue 1, 2001, p. 83-106. [Publisher]
  15. Zhang W., Wang B., Ye Z., Quan J. Efficient method for limit cycle flutter analysis based on nonlinear aerodynamic reduced-order models. AIAA Journal, Vol. 50, Issue 5, 2012, p. 1019-1028. [Publisher]
  16. Kou J., Zhang W. Layered reduced-order models for nonlinear aerodynamics and aeroelasticity. Journal of Fluids and Structures, Vol. 68, 2017, p. 174-193. [Publisher]
  17. Winter M., Breitsamter C. Nonlinear identification via connected neural networks for unsteady aerodynamic analysis. Aerospace Science and Technology, Vol. 77, 2018, p. 802-818. [Publisher]
  18. Wang Z., Xiao D., Fang F., Govindan R., Pain C. C., Guo Y. Model identification of reduced order fluid dynamics systems using deep learning. International Journal for Numerical Methods in Fluids, Vol. 86, 2018, p. 255-268. [Publisher]
  19. Bhatnagar S., Afshar Y., Pan S., Duraisamy K., Kaushik S. Prediction of aerodynamic flow fields using convolutional neural networks. Computational Mechanics, Vol. 64, 2019, p. 525-545. [Publisher]
  20. Mannarino A., Mantegazza P. Nonlinear aeroelastic reduced order modeling by recurrent neural networks. Journal of Fluids and Structures, Vol. 48, 2014, p. 103-121. [Publisher]
  21. Fu X., Li S., Fairbank M., Wunsch D. C., Alonso E. Training recurrent neural networks with the levenberg-marquardt algorithm for optimal control of a grid-connected converter. IEEE Transactions on Neural Networks and Learning Systems, Vol. 26, Issue 9, 2015, p. 1900-1912. [Publisher]
  22. Lee T., Gerontakos P. Investigation of flow over an oscillating airfoil. Journal of Fluid Mechanics, Vol. 512, 2004, p. 313-341. [Publisher]
  23. Roe P. L. Approximate Riemann solvers, parameter vectors and difference schemes. Journal of Computational Physics, Vol. 43, 1981, p. 357-372. [Publisher]