Airfoil buffet aerodynamics at plunge and pitch excitation based on long short-term memory neural network prediction

In the present study, a nonlinear system identification approach based on a long short-term memory (LSTM) neural network is applied for the prediction of transonic buffet aerodynamics. The identification approach is applied as a reduced-order modeling (ROM) technique for an efficient computation of time-varying integral quantities such as aerodynamic force and moment coefficients. Therefore, the nonlinear identification procedure as well as the generalization of the ROM are presented. The training data set for the LSTM–ROM is provided by performing forced-motion unsteady Reynolds-averaged Navier–Stokes simulations. Subsequent to the training process, the ROM is applied for the computation of the aerodynamic integral quantities associated with transonic buffet. The performance of the trained ROM is demonstrated by computing the aerodynamic loads of the NACA0012 airfoil investigated at transonic freestream conditions. In contrast to previous studies considering only a pitching excitation, both the pitch and plunge degrees of freedom of the airfoil are individually and simultaneously excited by means of an user-defined training signal. Therefore, strong nonlinear effects are considered for the training of the ROM. By comparing the results with a full-order computational fluid dynamics solution, a good prediction capability of the presented ROM method is indicated. However, compared to the results of previous studies including only the airfoil pitching excitation, a slightly reduced prediction performance is shown.


Introduction and motivation
Unsteady dynamic aeroelastic and aerodynamic phenomena, such as flutter and buffet, are of paramount importance concerning safety and efficiency requirements of passenger aircraft.Transonic buffet, also referred to as shock buffet, is characterized by shock-boundary layer interaction, resulting in self-sustained cycles of shock movement and partial flow separation.However, even in the absence of any structural excitation, these phenomena can occur.A comprehensive overview of computational and experimental studies of transonic buffet is given by Giannelis et al. [1].Several studies in recent years revealed that transonic shock buffet oscillations are characterized by a characteristic frequency, which is a dependent on airfoil shape, Mach number, and angle of attack [2,3].
Due to the self-sustained flow characteristics, the resulting aerodynamic forces and moments are characterized by unsteady variations.In an industrial context, unsteady Reynolds-averaged Navier-Stokes (URANS) simulations must be performed to capture the strongly nonlinear flow physics associated with shock formation and shock movement [2].In addition, the applied computational method must account for viscous effects to represent boundary layer interaction and flow separation.However, to resolve these relevant mechanisms, the effort of computational methods is still time-and cost-consuming, resulting in a high demand for fast and accurate surrogate methods to decrease computational time and costs.Therefore, a strong focus has been set on the further development of system identification and reduced-order model (ROM) methods in the past years.Considering aerospace applications in particular, a short overview of linear and nonlinear ROM approaches is given in the following: several studies concerning aerodynamic ROM approaches based on linear system identification have been conducted.Examples are based on the eigensystem realization algorithm (ERA) [4] applied by Silva and Bartels [5], the auto-regressive with moving average (ARMA) approach applied by Raveh [6], or the auto-regressive with exogeneous input (ARX) architecture used by Zhang and Ye [7].However, nonlinear system identification methods also capturing separated flow conditions and large-amplitude motions are still under development.Reliable approaches for the computation of, e.g., flutter and limit-cycle oscillations (LCO) are based on radial basis function (RBF) neural networks [8] or multilayer perceptron (MLP) neural networks [9,10].Moreover, ROMs based on Kriging interpolation [11] and Wiener type models [12] produce accurate predictions.Furthermore, Zhang et al. [13] applied a multi-kernel RBF neural network for modeling unsteady aerodynamics considering different flow conditions.In addition, ROMs based on fuzzy logic [14][15][16] yield accurate and reliable results for nonlinear system identification purposes.
However, in contrast to the application of ROMs for flutter and LCOs, challenges arise when considering buffet aerodynamics.Due to the self-induced unsteady flow condition even in the absence of any structural excitation, a transient output based on steady system inputs must be covered by the applied ROM method.Therefore, a feedback connection between the system input and output must be implemented [17].
So far, only a limited amount of studies dealing with ROMs in the context of transonic buffet is available.Bourguet et al. [18] performed buffet investigations using a proper orthogonal decomposition (POD) approach.Kou et al. [19] applied dynamic mode decomposition (DMD) to identify the dominant flow modes associated with transonic buffet.Furthermore, Winter and Breitsamter developed a nonlinear ROM approach based on a neuro-fuzzy model (NFM) that is serially connected with an MLP neural network [20].The application of this approach yields accurate results for the prediction of integral aerodynamic quantities associated with transonic buffet [17].
In the present paper, a nonlinear ROM approach based on a recurrent long short-term memory (LSTM) neural network is applied for an aerodynamic investigation of the transonic buffet phenomenon.The LSTM network is trained and utilized for the computation of integral aerodynamic quantities, in particular the lift and pitching moment coefficient, induced by transonic buffet.Therefore, a reference numerical data set is computed by means of URANS simulations.As a test case, the NACA00012 airfoil is investigated at transonic freestream conditions.Although the following training and application procedure has already been applied by the authors toward the NACA0012 airfoil undergoing user-defined pitch motion [21], the present study includes the investigation of strong nonlinear effects due to a userdefined individual and simultaneous excitation of the pitch and plunge DoFs of the NACA0012 airfoil.Compared to previous studies, which only considers a pitching excitation, a slightly reduced prediction performance of approximately 10% is indicated for all considered test cases.In addition, the training performance of the pitch-plunge-trained ROM is reduced in terms of training error and number of training iterations.However, since strong nonlinear effects are considered due to a simultaneous DoF excitation, the prediction performance of the trained LSTM is still sufficient.

Computational methods
In the following section, the applied ROM approach as well the computational methods for the generation of the training and validation data are briefly described.Therefore, a short introduction of basic recurrent neural networks (RNNs) is given, followed by a more detailed description of the applied LSTM neural network and the corresponding training procedure.The last subsection covers a short introduction of the applied CFD solver.

Recurrent neural network
RNNs are defined as a class of neural networks, characterized by internal self-connections, which allow for sequential information processing.Since all elements of the input sequence are processed in the same manner, the network is able to develop a memory based on outputs of current and previous time-steps.This characteristic allows the network to store information as memory data [22].
A common RNN architecture includes an input and output layer, which are characterized by feed-forward connections.Between these two layers, one or more hidden layers are included, which process data in a recurrent way.Within each layer, information is processed by means of activation functions, including weights and biases [23].
The training of the RNN is commonly performed by modifying its parameters using a gradient descent optimization [23].Therefore, a loss function quantifying the accuracy of the neural network is minimized to achieve a high network performance.Considering a prediction of time-series data, the loss function evaluates and minimizes the difference between the output as obtained by the network and a reference full-order solution.By applying gradient descent optimization, a loss function is computed at each time-step, followed by an update of the corresponding weights.This weight modification is accomplished by computing the gradient of the loss function by means of the backpropagation algorithm.Since the algorithm is applied at each time-step, it is also referred to as backpropagation through time (BPTT) [24].However, if the number of hidden layers is increased, the vanishing gradient problem can arise during the training of an RNN.Due to this issue, information from previous time-steps is not taken into account in an efficient way, possibly resulting in a loss of information provided by the data at the current considered time-step [25].

Long short-term memory neural network
To avoid the aforementioned vanishing gradient problem, a recurrent LSTM neural network architecture is applied in the present study to represent the nonlinear system dynamics associated with the transonic buffet phenomenon.The LSTM architecture was first published by Hochreiter and Schmidhuber [26] and represents a special type of RNNs, designed for time-series prediction based on both short-and long-term dependencies.In Fig. 1, the node architecture of an LSTM network, which is commonly referred to as a cell, is illustrated.
Within the LSTM network architecture, a system output y t is computed based on the input of the current time-step x t , the output from a previous time-step y t−1 , as well as memory data c t−1 from several previous time-steps.Since the filter- ing of memory data is accomplished by linear interactions, the LSTM network is able to take information from several previous time-steps into consideration for output prediction.This characteristic eliminates the aforementioned vanishing gradient issue [23].In contrast to a classical RNN architecture, an LSTM cell processes information by means of three different gates, namely the forget gate f t , the input gate i t , and the output gate o t (see Fig. 1).Within the gate structure, data processing is accomplished by either a sigmoid ( ) or tan h activation function and pointwise multiplication.
Within the forget gate f t , the output from the previous time step y t−1 as well as the input of the current time-step x t are collected (see Eq. 1) In Eq. 1, W f is defined as the corresponding gate weight matrix, while b f represents a bias vector.By means of an activation function, the data are either further processed or discarded from the previous cell state.In the present study, the sigmoid activation function is chosen (see Eq. 2) Within the input gate i t , input from the present cell state is processed by the sigmoid function (see Eq. 3) Subsequently, a new vector c t is created by means of a tanh activation function, including information that is potentially chosen to be added to the present cell state.Subsequently, the present cell state is updated by multiplying the selected data to the state vector c t (see Eq. 4) At the output gate o t , the output y t of the cell is defined by adding the cell data filtered by the activation function (see Eq. 5) to the aforementioned state vector (1) (3) Finally, the transferred memory data of the single cell is passed to consequent cells or the output layer of the LSTM network (see Eq. 6) Based on a numerical reference data set including the input and output system information, the parameters are updated in the training process using a backpropagation through time (BPTT) [24] approach.For minimizing the training error, which is defined as the mean-square error (MSE) in the present study, the adaptive moment estimation (ADAM) [27] algorithm is applied.By applying ADAM, the overall training data set is split into several units, also referred to as batches, which are processed to the LSTM network in a sequential manner.Furthermore, the algorithm includes a learning rate , which must be defined prior to the training process.However, in comparison to different gradient descent methods such as the Stochastic Gradient Descent (SGD), ADAM computes different learning rates for each weight of the network by means of gradient moment estimation.

Unsteady aerodynamic reduced-order modeling
In the present study, the LSTM neural network is applied for the identification process of the input and output system dynamics associated with the transonic buffet phenomenon.
Considering k as the current discrete time-step, the system identification can be defined as follows: The training and validation of the LSTM-ROM are accomplished based on a numerical data set, representing the characteristic system features.Therefore, CFD simulations are performed for a certain flow condition, while the structural DoFs are excited by means of an user-defined training signal.For a sufficient application, the signal must be defined to cover the amplitude and frequency range, which should be predicted by the ROM.Due to the numerical computation based on external motions, varying time-series of aerodynamic forces and moments result.By combining the excitation signal, which represents the system input, and the resulting CFD calculation output, the merged data set can be employed for the training and validation of the ROM.

CFD solver
The URANS simulations for the computation of the buffet phenomenon are conducted using the TAU code developed (5) by the German Aerospace Center (DLR) [28].The computation of the training, validation, and test data set are accomplished using the same settings.TAU solves the URANS equations in conservation form using a shock-capturing finite volume scheme.For the spatial discretization, a central scheme with matrix dissipation is applied.The discretization of the convective fluxes is accomplished using a secondorder central scheme, while the gradients are reconstructed using a Green-Gauss scheme.The temporal integration is performed using a backward Euler implicit scheme, while the embedded pseudo-time solution is computed by means of a lower upper symmetric Gauss-Seidel (LU-SGS) algorithm.To accelerate the computations, a multi-grid approach is applied.For turbulence modeling, the Spalart-Allmaras turbulence model with Edwards modification is applied.The mesh deformation for the respective pitch and plunge motion of the NACA0012 airfoil is implemented using a TAU-Python interface.

Test case: NACA0012 airfoil
To investigate the prediction capability of the LSTM-ROM, the training procedure introduced in Sect.2.3 is applied for the computation of the lift ( C L ) and pitching moment coef- ficient ( C My ) of the NACA0012 airfoil.Based on numerical studies performed by Raveh [2], the transonic buffet condition of the examined airfoil is defined by a free stream Mach number of Ma ∞ = 0.72 , a Reynolds number of Re = 10 7 , and an angle of attack of = 6 • .In Fig. 2, the distribution of the pressure coefficient at buffet condition, divided by four time-steps within a single buffet period T Buffet , is shown.Even without any structural excitation, the cyclic shock motion and the change in separation intensity are clearly visible.Based on the analysis of the resulting lift coefficient response at buffet condition, a reduced buffet frequency of k red = 2 f ⋅ c ref ∕U ∞ = 0.43 has been identified.For the pre- sent investigation, the chord length of the NACA0012 airfoil has been defined as c ref = 1 m, while the pitching axis is located at 25% of the chord length of the airfoil.
The employed, block-structured reference CFD grid is discretized by 28 ⋅ 10 4 elements using ICEM CFD.To ensure the independence of the solution from the spatial discretization, a grid-sensitivity study was conducted.Therefore, several grid refinement steps have been performed.In Fig. 3, the time-series of the lift coefficient is presented for the applied grid levels.Since the relative error between the grids discretized by 56 ⋅ 10 4 and 28 ⋅ 10 4 elements is computed as 0.13% and 0.12% with respect to the minimum and maximum amplitude of C L , the grid with 28 ⋅ 10 4 cells is cho- sen for the following study.To define the time-step for the simulations, a time resolution study has been additionally performed.Based on the study, the time-steps are chosen as = 0.11, corresponding to a physical time-step of t = 5 ⋅ 10 −3 s, since the relative error compared to a smaller time- step of t = 1 ⋅ 10 −3 s is given as 0.07%.However, increasing the time-step to t = 1 ⋅ 10 −2 s results in a relative error of almost 5%.

Generation of training data
For the training of the LSTM-ROM, the pitch and plunge DoFs of the NACA0012 airfoil are excited by means of an amplitude-modulated pseudo-random binary signal (APRBS).This signal type is characterized by a high information content per signal length, since a large range of various amplitudes and frequencies can be covered [20,29].Therefore, the application of an APRBS enables a reduction in simulation time and costs compared to other excitation signals.
For the following study, three different APRBS have been randomly defined, including 15,000 time-steps each.For the application, the signals have been scaled to limit the pitch amplitude to = ±1 • .The plunging motion is also limited to h = ±1 ⋅ c ref .
In addition, the scaled signals are smoothed to mitigate large gradients in the original signal.In Fig. 4, Analogous to the remaining two signals, the first and last parts of the APRBS are constructed for a separate excitation of the pitch (gray line) and plunge (blue dotted line) DoF, respectively.In between these two time instances, a combined excitation of both DoFs is accomplished to capture the nonlinear interactions caused by a combined input.The separate excitation for the pitch and plunge DoF is defined by 6000 time-steps each, whereas the combined excitation includes 4500 time-steps.

Nonlinear system identification
With respect to the application case in Sect.3, the input of the LSTM-ROM is defined by the pitching angle and the plunging degree of freedom h, whereas the lift and pitching moment coefficient denote the system output.Therefore, the ROM represents a multiple-input/multiple-output system, given by the following input ( x t ) and output vector ( y t ): Based on the numerical input and output data set, the LSTM-ROM is trained according to the aforementioned training procedure (see Sect. 2.3).Prior to the ROM training, the overall numerical data are divided into a training and validation data set.Therefore, two of the three available APRBS are applied for the training, whereas the remaining one is employed for ROM validation during training.
To define the hyperparameters for the training process, an extensive parameter study is conducted.Based on the study, the following hyperparameters of the LSTM network are defined: the number of hidden layers is varied between one and three; however, two hidden layers are chosen for the recent study.A single-layer LSTM indicates a lower performance quality, whereas the implementation (8) of three hidden layers intensively increases the training time.For each of the two hidden layers, the following number of neurons is considered: n neurons = [50, 100, 200].Since the implementation of 100 and 200 neurons also results in an increasing training time without any considerable performance improvement, the number of neurons per layer is set to 50.As already stated in Sect.2.2, the ROM is trained using gradient descent with an initial learning rate of = 0.001.An increase of the learning rate to = 0.01 results in an acceleration of convergence; however, the training performance decreases.Reducing the learning rate to = 0.0001 leads to a deceleration of convergence and therefore an considerable increase in training time.Since two of the three available APRBS are chosen for LSTM training, the training batch size is defined as two.The number of time-steps is equal to the number of time-steps covered by the training signals.Furthermore, instead of using the entire number of training samples in each batch, a sequence length of 100 time-steps is selected.Since both the excitation of the pitch and plunge DoF are considered, the number of input features is defined as two.Furthermore, as already stated in Sect.2.2, a tan h is selected as the state activation function, while the gate activation function is chosen as .The convergence trends of the LSTM-ROM during the training process are visualized in Fig. 5 Prior to the application of the LSTM network for unknown data prediction, the LSTM network is applied to the training data itself by performing multi-step ahead predictions.In Fig. 6, the results of the lift coefficient due to the excitation of the APRBS, as shown in Fig. 4, are given.The results as computed by the LSTM-ROM are compared to the reference CFD solution.
For error quantification, the fit factor Q i [30], as defined in Eq. ( 9), is computed.Therefore, the model output ŷ is represented by the ROM response.A fit factor of 100% defines an exact agreement between the ROM and the CFD results [20] Fig. 4 Time-series of the scaled and smoothed training signal (APRBS) for a separate and combined excitation of the pitch and plunge degree of freedom Based on the results of the reference CFD simulation and the ROM prediction, fit factors of Q i (C L ) = 86.36%and Q i (C My ) = 85.12% are computed for the lift and pitching moment coefficient, respectively.As it can be seen, the ROM prediction of the lift and pitching moment coefficient yields a good agreement with the CFD solution, which is a minimum requirement for further application tests.

Application of the LSTM-ROM
. pitching amplitude is defined as = ± 1 • , whereas the plunge motion is limited to h = ± 0.1 ⋅ c ref .
The refer- ence full-order CFD simulations are computed based on the aforementioned TAU setup.To compare the ROM and CFD solutions without the influence of any initial transient errors, in total, 14 excitation periods have been computed for each DoF with both the CFD solver and the LSTM-ROM.Furthermore, the combined excitation period is represented by seven cycles.To provide a clear visualization of the results, time-domain as well as frequency-domain responses are presented in the following.The frequency-domain responses are visualized using a Fast-Fourier Transformation (FFT), plotted over the reduced frequency ( k red ).
Evaluating the time-and frequency-domain response of the lift coefficient (see Figs. (8,9), a nonlinear and frequency-dependent interaction between the airfoil motion and the buffet phenomenon is visible.The same interaction is also visible in the resulting pitching moment coefficient (see Fig. 10).This interaction has already been studied by Raveh [2] and Winter and Breitsamter [17] for the considered NACA0012 airfoil test case.For a reduced frequency of k red,Ex = 0.2, which is smaller than the corresponding shock buffet reduced frequency ( k red,Buffet = 0.43), an inter- action of the buffet frequency and the airfoil motion frequency is indicated.The influence of the shock buffet can be clearly distinguished from the effects resulting from the combined airfoil excitation.Furthermore, the buffet frequency appears to be less distinctive than the excitation frequency of k red,Ex = 0.2.In addition, also other frequency peaks at around k red = 0.4 and k red = 0.6 occur, which might be a result of the high nonlinear interaction.
In contrast, considering an excitation frequency of k red,Ex = 0.4, which is very close to the buffet frequency, the influence of buffet decreases to a large extent and the system is predominantly governed by the applied airfoil motion [17].This characteristic is known as the "lock-in" effect [2].Considering the remaining application cases, also a decrease in the buffet frequency is indicated.For both cases, most of the excitation response is represented by the airfoil excitation motion.However, at an airfoil excitation of k red,Ex = 0.8, the peak of the buffet frequency ( k red,Ex = 0.43) is larger compared to the excitations of k red,Ex = 0.4 and k red,Ex = 0.6.Although the underlying flow physics of the buffet phenomenon are nonlinear to a high extent, a comparison of the reference CFD results with the results of the LSTM-ROM yields a good prediction capability of the selected test cases.Some less dominant frequency peaks are not fully captured by the LSTM network; however, the most dominant peaks are represented correctly in frequency and amplitude.Considering the computed fit factors, which are summarized in Table 1, a sufficient performance quality is further emphasized.In addition, the LSTM-ROM is able to capture the lock-in effect originating from the separate and combined excitation of the pitch and plunge DoF.

Computational effort
In the present work, the computations for the CFD reference solutions have been performed on the SuperMUC-NG of the Leibniz Supercomputing Centre, using 6 nodes with 48 cores each, resulting in 288 cores.In contrast, the ROM computations have been conducted on a state-of-the-art workstation using a single Intel Xeon 2.2 GHz processor.
The computation of the APRBS forced-motion CFD simulation required approximately 72 h on the SuperMUC-NG.Therefore, considering the number of applied cores, an overall computational time of 20,736 CPU hours results.In contrast, the LSTM-ROM training process was conducted within approximately one CPU hour on the workstation.Thus, the overall computational effort for the ROM training is summed up to 20,737 CPU hours, while the computations of the APRBS response have the highest share on the overall ROM training process.
Due to the dependency of the oscillation period on k red,Ex , a different number of computed time-steps is set for each of the harmonic motions with varying frequencies.Therefore, an averaged computation time of 48 h with the TAU solver is assumed for each simulation, resulting in a total computational time of 55,296 ( 4 * 48 * 288 ) CPU hours.In contrast, the computation of the LSTM-ROM required on average 2 h for each harmonic motion excitation, resulting in a total computational time of 8 CPU hours.Therefore, considering the number of selected test cases and the APRBS training computation, a reduction of computational time by a factor of 2.66 was achieved.However, without considering the computational time necessary for the computation of the training data, for every additional test simulation, a reduction by four orders of magnitude is possible.

Conclusion
In this paper, an unsteady aerodynamic reduced-order model (ROM) based on a long short-term memory (LSTM) neural network was applied to predict integral aerodynamic forces and moments associated with transonic buffet.Therefore, a training data set computed by means of unsteady Reynoldsaveraged Navier-Stokes simulations has been provided.To demonstrate the prediction capability of the ROM, the NACA0012 airfoil was investigated at a transonic buffet condition, undergoing individual and simultaneous pitch and plunge motions.Based on the results, a precise performance quality for predicting both the self-sustained unsteadiness and the characteristics resulting from the external excitation was emphasized.Furthermore, the overall computational effort was reduced by four orders of magnitude compared to a full-order computational fluid-dynamics solution for every additional simulation, neglecting the production time of the training data.However, compared to previous studies considering only an individual pitching excitation, a decrease  in prediction performance of approximately 10% was indicated for all considered test cases.In addition, the training performance of the LSTM-ROM was reduced.However, since the trained ROM has been applied for the prediction of strong nonlinear effects due to a simultaneous DoF excitation, the resulting performance quality is still sufficient.
For future studies, it is intended to apply the LSTM-ROM to three-dimensional aircraft configurations.Furthermore, coupling the ROM with a structural model will allow for the computation of aeroelastic buffeting.

Fig. 1
Fig. 1 Architecture of a LSTM cell

For
the final model generalization test, sine signals representing different reduced frequencies k red,Ex = [0.2,0.4,0.6,0.8] are defined for the excitation of both DoFs.The test signal representing a reduced frequency of k red,Ex = 0.2 is exem- plary shown in Fig. 7. Analogous to the training signal, the (9)