Unsteady aerodynamic identification based on 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. 1Architecture 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:

where $y$ is the output of the system; $f$ represents the transfer function of the system; $x$ is the input of the system; $f_1$ and $f_2$ denote the transfer functions of the neuron, respectively; $W_1$ and $W_2$ represent the weight of the neuron, respectively; and $b_1$ and $b_2$ 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:

where $y\left(t\right)$ denotes the output of the system at $t$ moment; and $x\left(t\right)$ 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. 2Architecture 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]:

where $E\left(x\right)$ is the minimum value of multivariate function, $i$ is the training pattern index and $F_i\left(x\right)$ is the error for pattern $i$, which can be defined by the following formula:

where $y\left(x,t_i\right)$ and $\phi \left(t_i\right)$ are scalar functions.

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

where $J\left(x_k\right)$ is the Jacobian matric of $F\left(x_k\right)$; ${\lambda }_k$ is a scalar regulation which controls the modulus sum direction of $d_k$ and $I$ is the identity matrix. When ${\lambda }_k$ is zero, the direction $d_k$ is the same as the Newton method. As ${\lambda }_k$ tends to infinity, $d_k$ tends to the steepest descent direction and the mode tends to zero, which means that when ${\lambda }_k$ is large enough, the term $F\left(x_k+d_k\right)<F\left(x_k\right)$ holds. Therefore, the term ${\lambda }_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.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. 3Nested 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–\omega$ 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:

where ${\alpha }_0$ is the average angle of attack; ${\alpha }_1$ is the amplitude angle of attack; $k$ is the reduction frequency; and $\tau$ is dimensionless time (reduction time). $k$ and $\tau$ can be defined as follows [23]:

where $u_{\infty }$ 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]:

where $\rho$ is the density of inflow air; and $\mu$ 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. 4CFD results compared with experiments: a) lift coefficient, b) moment coefficient

a)

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.

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:

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. 5Lift 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:

Fig. 6Moment 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. 7Lift coefficient prediction results of verification signal

Fig. 8Moment 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:

where $\widehat{y}\left(i\right)$ is the prediction results of the RNNs; $y\left(i\right)$ 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.