Neuro-Fuzzy Network-Based Reduced-Order Modeling of Transonic Aileron Buzz

: In the present work, a reduced-order modeling (ROM) framework based on a recurrent neuro-fuzzy model (NFM) that is serial connected with a multilayer perceptron (MLP) neural network is applied for the computation of transonic aileron buzz. The training data set for the speciﬁed ROM is obtained by performing forced-motion unsteady Reynolds-averaged Navier Stokes (URANS) simulations. Further, a Monte Carlo-based training procedure is applied in order to estimate statistical errors. In order to demonstrate the method’s ﬁdelity, a two-dimensional aeroelastic model based on the NACA65 1 213 airfoil is investigated at different ﬂow conditions, while the aileron deﬂection and the hinge moment are considered in particular. The aileron is integrated in the wing section without a gap and is modeled as rigid. The dynamic equations of the rigid aileron rotation are coupled with the URANS-based ﬂow model. For ROM training purposes, the aileron is excited via a forced motion and the respective aerodynamic and aeroelastic response is computed using a computational ﬂuid dynamics (CFD) solver. A comparison with the high-ﬁdelity reference CFD solutions shows that the essential characteristics of the nonlinear buzz phenomenon are captured by the selected ROM method.


Introduction
Unsteady aerodynamic and aeroelastic phenomena, such as flutter, buffet and buzz, determine the boundaries of the flight envelope of an aircraft. The latter phenomenon, also referred to as aileron buzz, represents an aeroelastic phenomenon occurring in transonic flights. Aileron buzz is denoted as a single degree-of-freedom flutter, characterized by self-excited and sudden vibration of the aileron, yielding large amplitudes and potentially permanent damage of the aircraft structure. The associated flow field typically involves shock waves dynamics and boundary layer separation. For viscous flows, aileron buzz was classified by Lambourne [1] as Type A, B and C, depending on the free stream Mach number and the location of the shock wave. Type A is caused by the interaction of shock-waves with the boundary layer, typically encountered at Mach numbers slightly higher than the critical Mach number. Type B results from the interaction of the shock-waves with the aileron movement at higher Mach numbers, compared to Type A. Since for Type B no significant intervention with the boundary layer is observed, it is commonly referred to as non-classical aileron buzz [2]. Type C buzz is associated with supersonic flow over the entire control surface and shock wave development at the trailing edge of the aileron [1,3].
In order to accurately capture the flow physics of aileron buzz, high-fidelity numerical methods based on unsteady Euler or unsteady Reynolds-averaged Navier Stokes (URANS) simulations must be employed [2,4]. However, due to the necessity to resolve the underlying flow mechanism associated transonic buzz condition. Therefore, a high-fidelity model of the aeroelastic system is defined by coupling the CFD aerodynamic model with the dynamics of the rigid aileron.

Reduced-Order Model Approach
In order to identify the underlying characteristics of the aerodynamic and aeroelastic system, a data set including the characteristic features of the system must be defined for an efficient ROM construction. For this purpose, unsteady CFD simulations are performed for predefined flow conditions, while the selected rigid body degree of freedom (DoF) of the structure, in particular the aileron structure, is excited by means of an user-defined signal, covering the amplitude and frequency range of interest. Due to the external motion, a time-series of forces and moments acting on the investigated structure, results. Combining the training signal, which represents the system input, with the CFD-based output, the merged data set can be employed for ROM training purposes [5].
The applied NFM-MLP-based ROM is based on a time-discrete external dynamics filtering approach. Based on this approach, a nonlinear function F of current system inputs combined with previous, time-delayed system inputs and outputs can be defined for the representation of the underlying dynamic system [5]. For the representation of the unknown function F, the NFM and MLP neural network are arranged in a series connection. An illustration of the connected ROM is given in Figure 1. The training of the nonlinear system identification approach is basically composed of three main steps, which are briefly discussed in the following. In the first stage of the NFM-MLP training, the NFM is independently trained by means of the local linear model tree (LOLIMOT) algorithm [25]. A local linear neuro-fuzzy model including N local linear sub-models can be expressed as follows: withŷ defining the scalar output of the NFM related to the input vector x. The model parameters of the ith LLM are given by the respective weights (w i0 , w i ), the centers c i and basis function widths Σ i .
Further, ψ i is defined as a fuzzy membership function composed of Gaussians based on the Euclidean distance between the input vector and the center c i of the ith neuron. [12] Each LLM is connected to a single membership function. Based on the LOLIMOT algorithm, the structure and the parameters of the NFM are optimized. In particular, the structure optimization is accomplished by defining a sufficient number of LLM and neurons. The LOLIMOT training procedure can be divided in five steps, which are briefly introduced in the following.

1.
Model Initialization: As a first step, a global linear model (ψ 1 = 1, N = 1) is calculated by estimating the NFM weights (w i0 , w i ) by means of a local linear least-squares algorithm [25]. Therefore, the available CFD training data set is applied. 2.
LLM error estimation: Within the next step, the worst performing local linear model is located by means of a locally-defined loss function, which is evaluated for all available models (i = 1, .., N). The LLM yielding the highest prediction error is selected for the subsequent splitting procedure [12]. If only the global model is available, it is automatically chosen for further refinement. 3.
LLM refinement: The LLM which shows the lowest performance is divided into two models using an axis-orthogonal split [25]. For each resulting model, the centers (c i ), widths (σ) and weights (w i0 , w i ) must be re-computed. The centers are defined as the centers of the corresponding jth hyperrectangle, while the widths are determined by defining the input space extension of the LLM scaled by a factor k σ [12]. According to Nelles [25], k σ is chosen as 1/3. The linear weights are determined by the application of the local weighted least-squares optimization.

4.
Error evaluation: In order to evaluate the best splitting configuration, a loss function is applied. Therefore, in contrast to the weight estimation procedure in training step one, the available validation data set is applied. Based on the error evaluation, the partition-setup yielding the lowest error is chosen for the last training step.

5.
Termination: The aforementioned training steps are repeated until the relative change of the error as calculated in the previous step becomes smaller than a user-defined value. As an alternative, the splitting process can be terminated by defining a maximum number of LLMs [12].
Considering time-series prediction tasks, the NFM input vector x is defined based on (time-delayed) inputs (q(k),...,q(k − m)) and outputs y, whereat the latter are used for initial parameter estimation. Therefore, no iterative feedback is performed within the training process of the NFM.
Within the second step of the NFM-MLP training, a NFM simulation is performed based on the available training data set. Therefore, within each considered time step the NFM is applied in order to obtain the next prediction output. By means of the input vector x, the output of the NFMŷ is processed in a recurrent manner. Based on the training and simulation step, two output data sets containing the original training response and a multi-step ahead time-series prediction are available for the following application [12].
The last step of the overall NFM-MLP training process includes the individual training of the MLP neural network. A MLP network with a single hidden layer composed of M neurons can be expressed as follows:ỹ In Equation (2),ỹ is denoted as the output of the MLP network based in the input vector v. Considering the serial connected ROM, the input vector of the MLP network is obtained by combining the (time-delayed) inputs and the (time-delayed) outputs of the NFM. Since the complete output of the NFM is already available, the information of the current considered time step (k)ŷ(k) is summarized in the MLP input vector. Therefore, the corresponding MLP output includes only known information. The respective linear and nonlinear weights (g i and G i ) of the MLP network are initialized using the Nguyen-Widrow method and are optimized by means of a backpropagation approach combined with a Levenberg-Marquardt optimization algorithm. Due to the serial model structure, a dynamically stable response is guaranteed, since the MLP outputỹ is not processed in a recurrent manner. Therefore, the MLP can be considered as a nonlinear correction of the NFM output [12].
With respect to the application case described in Section 3.1, the system input q of the combined NFM-MLP ROM is defined by the aileron deflection angle δ, whereas the aerodynamic force acting on the aileron, represented by the hinge moment coefficient C M H , defines the system outputỹ. Therefore, the system identification approach for the representation of the predicted hinge moment coefficientC M H can be written as follows, with k being defined as the current considered discrete time In Equation (3), m ∈ N and n ∈ N refer to the maximum dynamic delay-orders for the respective input and output quantities. In order to avoid stability issues, the maximum output delay should be chosen larger or equal to the maximum input delay [5].

Structural Model
In order to investigate the performance of the NFM-MLP-based ROM concerning the prediction of aileron buzz aerodynamics, the aforementioned system identification procedure is applied to the NACA65 1 213 airfoil with an integrated aileron. The aeroelastic model in the present work represents the coupling between the structural degree of freedom of the aileron, defined by the aileron deflection angle δ and the aerodynamic response, given by the hinge moment M H . The structural model is represented by a fixed two-dimensional wing section with a reference length of c re f = 1 m. The aileron is hinged by the three-quarter chord location (x H = 75%) and its degree of freedom is modeled as rigid. Therefore, structural elastic and dissipative contributions are neglected in the present study. The motion of the aileron is described by the aileron deflection angle δ(t) about the hinge point, whose dynamics are expressed by the following one-degree of freedom equation: Here, I H denotes the aileron mass moment of inertia about the hinge line and k is defined as a spring constant. In the present study, k is assumed to be zero, since the aileron is only constrained by the hinge point. Further, a downward deflection of the aileron is indicated by a negative δ and vice versa.
For the aeroelastic investigation, the structural model including the aileron deflection is implemented in ANSYS Fluent by means of a local mesh deformation mechanism. In particular, the deformation is accomplished by building up the airfoil shape by overlaying the thickness distribution of the airfoil with a parametrically defined camber line y c according to [26]. The camber line associated to the aileron is defined by a third-order polynomial function, which is parameterized in order to allow for an update of the control surface deflection and the control of its minimum and maximum amplitude. The resulting parameterization is defined as follows: with t and t H defining the current time and the start time of the deflection, respectively. The current position of surface nodes on the aileron is indicated by x, whereas the deflection angle is defined by δ. In Figure 2, the geometry of the NACA65 1 213 airfoil with the implemented positive and negative aileron deflection is visualized. In Figure 3, the flow field including the respective aileron motion is visualized by Mach number contour plots. Therefore, four time steps of a single buzz period T Buzz are selected. In order to verify the implemented buzz modeling approach, a numerical test case has been conducted and compared to experimental results as obtained by Erickson [27]. Therefore, the buzz flow conditions has been defined as Ma ∞ = 0.82, Re = 20.7 × 10 6 and α = −1 • with a minimum and maximum deflection angle of δ min,max = ±9 • . In Figure 4, a comparison between the experimental and numerical results in terms of a variation of the local shock position with varying deflection angle is visualized for a single aileron deflection period. As it can be seen, a sufficient agreement is achieved.

CFD-Solver
In order to account for viscous effects and boundary layer interaction, an URANS approach is selected in the present study. The simulations capturing the buzz phenomenon are performed using ANSYS Fluent. Both the training and validation data time-series of the integral aerodynamic quantities are computed using the same CFD setup. ANSYS Fluent solves the URANS equations in conservation form using a shock-capturing finite-volume method. The temporal integration is performed using a first order implicit dual-time stepping approach, while the spatial discretization is achieved using a Roe flux-difference splitting (FDS) scheme. The reconstruction of the gradients is accomplished by means of a least squares cell-based scheme. In order to account for efficient turbulence modeling, the standard Spalart-Allmaras turbulence model is applied. The mesh deformation for the respective aileron motion of the NACA65 1 213 test case is implemented by means of a dynamic mesh deformation, which is defined by a user defined function.

Training Data Generation
For an efficient training of the ROM, an amplitude-modulated pseudo-random binary signal (APRBS) [12] is selected for the excitation of the aileron structure. The signal used in the present work is shown in Figure 5. For the following application, the minimum and maximum amplitude (δ) of the signal are chosen as δ = ±12 • . According to a study proposed by Steger and Bailey [4], aileron buzz of the NACA65 1 213 airfoil appears between a minimum and maximum aileron deflection amplitude of δ min = −12 • and δ max = +9 • , respectively. Therefore, the selected training signal covers the amplitude range of interest. Further, according to Steger and Bailey [4], aileron buzz of the NACA65 1 213 airfoil is characterized by varying freestream conditions. Therefore, in the present study three different flow conditions are considered, defined by a freestream Mach number of Ma ∞ = [0.8, 0.82, 0.83], a Reynolds number of Re = 20.7 × 10 6 and an angle of attack of α = −1 • . For the first and third condition the initial aileron deflection is constrained to δ start = 0 • , whereas the deflection of the second condition is defined by δ start = 4 • . Further, the corresponding reduced frequencies of the selected buzz conditions are summarized in Table 1: For the training of the ROM the flow condition defined by Ma ∞ = 0.82, Re = 20.7 × 10 6 , α = −1 • and δ start = 0 • is selected, whereas the remaining conditions are used for the subsequent application process.  The hybrid computational grid for the viscid aerodynamic and aeroelastic investigations was built using ICEM CFD. The grid is composed of structured cells around the airfoil surface and in the farfield, while the cells around the aileron are unstructured in order to allow for the implementation of large grid deformations. In order to guarantee a solution that is independent from the spatial resolution, a grid sensitivity study was conducted based on the steady state computation (Ma ∞ = 0.8, Re = 20.7 × 10 6 , α = 0 • , δ start = 0 • ) of C M H (see Figure 6). In this regard, in total four grids with the respective halved and doubled number of elements have been investigated. Since the relative error between the grids containing 12.48 × 10 5 and 24.96 × 10 5 elements is given as 0.65% with respect to C M H , the grid with 12.48 × 10 5 was chosen to be adequate for the following study. Based on a subsequent convergence study, the nondimensional time step for the simulations has been defined as ∆τ = 0.135, corresponding to a physical time step of ∆t = 0.0005 s.

Nonlinear System Identification
Based on the computed input and output CFD data set, the NFM-MLP-based ROM is trained according to the workflow discussed in Section 2. Therefore, several hyperparameters are defined for the NFM-MLP ROM. The number of neurons for the MLP network containing two hidden layers has been chosen as [5,2] for the respective first and second hidden layer. Based on the LOLIMOT algorithm, the splitting procedure of the local linear models (LLMs) is terminated using an error threshold of 1%, resulting in a maximum number of 20 LLMs. The maximum dynamic delay orders (m, n) for the input and output quantities are defined by applying the Lipschitz index as proposed by He and Asada [28]. In the present study, the iteration process for both delays is terminated due to a relative change of the Lipschitz index over one iteration of 1%. Based on this approach, both m and n are defined as 150.
In the present study, the available data set is segmented prior to the training procedure. Therefore, 70% of the data are randomly chosen and used for the estimation of initial parameters, whereas the remaining data is exploited for validation purposes in order to guarantee an accurate model. In order to circumvent statistical errors originating from this random segmentation procedure, a Monte Carlo-based training method is applied [12]. Therefore, the training of the NFM-MLP ROM is performed by means of N MC = 10 Monte Carlo iterations. In the present work, the meanȳ as well as the standard deviation σ are considered in order to analyze the data obtained by the Monte Carlo procedure. In order to guarantee a clear illustration of the results, all following diagrams only include the mean response of all parallel trained NFM-MLP ROMs.

Aerodynamic System Identification
In a first step, the NFM-MLP ROM is applied for an aerodynamic investigation, without considering any structural characteristics (I H = 0). Therefore, the flow condition is defined as the training condition (Ma ∞ = 0.82, Re = 20.7 × 10 6 , α = −1 • and δ start = 0 • ) with a minimum and maximum deflection amplitude of δ min = −12 • and δ max = +9 • , respectively. However, prior to the aerodynamic application, the ROM is evaluated based on the training data itself by performing multi-step ahead predictions. In Figure 7 the resulting hinge moment coefficient trend due to the excitation as computed by the ROM compared to the CFD reference solution, is shown.
For error evaluation purposes, the fit factor Q i [12,29] is introduced in Equation (6). Therefore, Q i is computed by defining the mean ROM response as the model outputỹ. A computed fit factor of 100% indicates an exact agreement of the CFD and ROM results. The resulting fit factor for the multi-step ahead prediction of the hinge moment coefficient is calculated as Q(C M H ) = 92.4%. Subsequent, the ROMs performance is investigated by means of harmonic pitch oscillations of the aileron. Therefore, varying reduced excitation frequencies are considered (k red,Ex = [0.5, 0.6, 0.7, 0.8]), which also include the reduced frequencies of the selected buzz cases (see Table 1). The reference full-order CFD computations are employed using the aforementioned ANSYS Fluent setup. In order to ensure a solution without any transient influence, several excitation periods have been computed with both the CFD solver and the ROMs for each reduced frequency. Analog to the previous study, the Monte Carlo procedure is applied in order to obtain the results presented below. Therefore, as mentioned in the previous section, all following diagrams show the mean response of the NFM-MLP ROMs. In Figure 8, a comparison of the ROM and the reference CFD results are demonstrated for each considered reduced frequency. In addition, frequency domain responses using a Fast-Fourier-Transformation (FFT) are visualized in Figure 9 for the selected excitation frequencies. As shown, the trained model is clearly able to reproduce the resulting hinge moment coefficient of the test cases. Considering the fit factors summarized in Table 2, the high prediction capability is emphasized.

Aeroelastic System Identification
Subsequent to the aerodynamic investigation, the ROM is employed for a coupled aeroelastic investigation. Therefore, according to [30], the aileron moment of inertia is defined as I H = 0.554. In the first step, the system identification is applied in order to represent a stable aeroelastic response. Therefore, according to Fusi [31], the selected flow condition is defined by a free stream Mach number of Ma ∞ = 0.8, a Reynolds number of Re = 10 6 and an angle of attack of α = 0 • . The recorded output of the CFD simulation and the computation of the ROMs are visualized in Figure 10. The corresponding fit factor is calculated as Q(C M H ) = 94.78%, which emphasizes a good correlation of both simulation results. In a second step, the prediction capability of the ROM is employed at the selected buzz conditions as described in Section 4.1. For the investigation, the CFD-based aeroelastic response and the ROM output are compared in terms of aileron deflection angle and hinge moment coefficient. The first considered condition is defined as the training condition, however the initial aileron deflection angle is set to δ start = 4 • . A comparison of the CFD and ROM results is visualized in Figure 11. As it can be seen, a slight discrepancy is indicated for the first deflection periods, however, after the first time steps an accurate agreement is visible. Considering the corresponding fit factor (see Table 3), a high prediction capability is emphasized.
Subsequent, the trained ROM is applied towards the remaining considered buzz cases in order to investigate the performance regarding varying Mach numbers which are not included in the training range. In Figure 12, the CFD and ROM result for Ma ∞ = 0.83, Re = 20.7 × 10 6 , α = −1 • and δ start = 0 • are visualized. Compared to the test case included in the training range, a larger disagreement between the CFD and ROM result is visible for the initial deflection periods. However, with proceeding simulation time, a high agreement is clearly indicated, which is also emphasized by the corresponding fit factor. For the third test case condition, also a slight discrepancy is indicated for the first deflection periods, while the compliance of the CFD and ROM results improves after the initial deflection periods (see Figure 13). However, in comparison to the other test cases, a higher disagreement is indicated by the corresponding fit factors. This issue might result from the selected Mach number, which deviates more from the training condition (Ma ∞ = 0.82) than the test case with Ma ∞ = 0.83. However, based on the results it is clearly visible that the NFM-MLP ROM is able to reproduce buzz conditions defined by varying freestream parameters, which are not included in the training data, with sufficient accuracy.

Computational Effort
All computations presented in this study have been performed on a workstation using a single Intel Xeon 2.2 GHz processor. The CFD computations have been conducted using sixteen cores, whereas the ROM results have been obtained using a single core. The computation of the APRBS forced-motion CFD simulation for the generation of the training data set required approximately 1150 CPU hours. In contrast to that, the training of the NFM-MLP ROM including the Monte-Carlo-based training procedure was performed within two CPU hours. Thus, the overall computational effort for obtaining the ROM training results is summed up to approximately 1152 CPU hours. Consequently, only the unsteady CFD simulations of the training signal have a noteworthy contribution to the overall ROM training process.
Due to the dependency of the oscillation period on k red,Ex , a different number of computed time steps result for each of the harmonic motions with varying reduced frequencies (k red,Ex = [0.5, 0.6, 0.7, 0.8]). Therefore, an averaged computation time of 48 CPU hours with ANSYS Fluent is assumed for each simulation. Considering the number of applied cores, a total computational time of 3072 CPU hours result. In addition, the computation of the coupled aeroelastic simulations for each of the considered buzz test cases was carried out within approximately 36 CPU hours, resulting in a total computational effort of 1728 CPU hours.
In contrast, the computation of the NFM-MLP ROM required on average two hours regarding each harmonic aileron motion, resulting in a total computation time of 8 CPU hours. For the aeroelastic simulations approximately 2.5 CPU hours were needed, yielding a total computation time of 10 CPU hours. Therefore, in the present study a reduction of computational time for the aerodynamic and aeroelastic investigations by a factor of 240 and 57.6 was achieved, considering the number of test cases and the APRBS training effort. Regarding the application of the ROM itself without taking the training effort into account, a reduction of computational time by two orders of magnitude is achieved.

Conclusions
In the present study, an identification-based reduced-order modeling approach based on a recurrent neuro-fuzzy model (NFM) that is serially connected with a multilayer perceptron (MLP) neural network has been applied to predict aerodynamic and aeroelastic characteristics associated to non-classical aileron buzz. Therefore, suitable training and validation data sets with respect to various flow conditions were provided based on the NACA65 1 213 airfoil test case. The structural model is represented by a fixed two-dimensional airfoil geometry with an integrated, rigid aileron. A one-degree of freedom aeroelastic model is implemented in the solver, which allows for an appropriate aeroelastic coupling for the representation of non-classical aileron buzz.
The trained ROM was applied for an aerodynamic and subsequent aeroelastic investigation. In particular, the aim of this study was the efficient prediction of the aileron deflection angle and the corresponding response, defined by the hinge moment coefficient. Based on the results it was indicated that the ROM is able to reproduce the motion-induced characteristics with sufficient accuracy. Further, buzz test conditions outside the training range have been predicted by the ROM with a high degree of accuracy. Besides, a reduction in computational time by a factor of 240 and 57.6 for the aerodynamic and aeroelastic investigations, respectively, was achieved. Regarding the ROM application itself without considering the training effort, computational time was reduced by a at least two orders of magnitude for the aerodynamic and aeroelastic simulations.