Nonlinear Unsteady Aerodynamics Reduced Order Model of Airfoils Based on Algorithm Fusion and Multifidelity Framework

A reduced order modeling method based on algorithm fusion and multifidelity framework for nonlinear unsteady aerodynamics is proposed to obtain a low-cost and high-precision unsteady aerodynamic model. This method integrates the traditional algorithm, intelligent algorithm, and multifidelity data fusion algorithm. In this method, the traditional algorithm is based on separated flow theory, the intelligent algorithm refers to the nonlinear autoregressive (NARX) method, and the multifidelity data fusion algorithm uses different fidelity data for aerodynamic modeling, which can shorten the time cost of data acquisition. In the process of modeling, firstly, a multifidelity model with NARX description provides a general intelligent algorithm framework for unsteady aerodynamics. Then, based on the separated flow theory, the correction equation from low-fidelity model to highfidelity result is constructed, and the cuckoo algorithm based on chaos optimization is used to identify the parameters. In order to verify the effectiveness of the method, an unsteady aerodynamic model of NACA0012 airfoil is established. Three kinds of data with low, medium, and high fidelity are used for modeling. The low-fidelity and medium-fidelity data is obtained from the CFD-Euler solver and CFD-RANS solver, respectively, while the high-fidelity data comes from the experimental results. Then, the model is established, and its prediction of unsteady aerodynamic coefficients is in good agreement with the CFD results and the experimental data. After that, the model is applied to a two-dimensional aeroelastic system, and the bifurcation and limit cycle response analysis are compared with the experimental results, which further shows that the model can accurately capture the main flow characteristics in the flow range of low speed and high angle of attack. In addition, the convergence of the model is studied; the accuracy and generalization ability as well as applicability scope of the model are compared with other aerodynamic models and finally discussed.


Introduction
Most of the problems in aircraft aerodynamic design are closely related to unsteady aerodynamic forces. Aerodynamic modeling is a major foundation for flight overall parameter design, flight trajectory design, and flight maneuverability and stability analysis. The accuracy of the aerodynamic model directly affects the performance of the control system and the reliability of flight simulation [1]. There are many factors that can affect the aerodynamic force such as flight altitude, airspeed, movement form, and rudder interference, which leads to a strong nonlinear aerodynamic system. However, the traditional aerodynamic modeling method is based on the assumption of quasisteady and linearization theory, which cannot meet the requirements for advanced civil and military aircraft design. So, in recent years, unsteady and nonlinear aerodynamic modeling methods have been paid more and more attention. There are two major problems in the unsteady aerodynamic modeling, one is the calculation of dynamic derivative, and the other is the aerodynamic modeling at a high angle of attack. As the former directly determines flight quality analysis and flight control system design, the latter is one of the important ways to evaluate the performance of aircraft at a high angle of attack. Therefore, it is of great value to study the flight dynamic derivative identification methods as well as nonlinear and unsteady aerodynamic modeling at high angles of attack.
The unsteady aerodynamic modeling method appeared at the beginning of the 20th century, and a variety of modeling methods had been widely developed since the 1980s. These methods can be divided into two categories: traditional mathematical expression method and intelligent learning methods. The traditional mathematical expression methods are a kind of method which obtain mathematical relationships between aerodynamic force, flight states, control inputs, and other parameters according to the physical characteristics and statistical laws of aerodynamic force. The traditional aerodynamic model is based on a dynamic stability derivative at the beginning, which is used to predict the change of flight load after the flight parameters change. These kinds of models regard each dynamic stability derivative as a constant value independent of time and other variables, so these models are linear aerodynamic models and are only applicable to the linear flight states with a small angle of attack [2]. In order to expand the application range of these models, the nonlinear algebraic model [3], Fourier series model [4][5][6], integral equation model [7], state space model [8,9], etc., were developed. However, generally speaking, the traditional mathematical expression method of aerodynamic modeling is more suitable for aerodynamic modeling with limited flight state range changes, such as the linear/steady aerodynamic modeling under certain altitude, speed, angle of attack, and other conditions. The relationship between aerodynamic force and flight height, velocity, angle of attack, sideslip angle, angular velocity, rudder angle, and other parameters is highly nonlinear for the aircraft with large changes in airspace, airspeed, flight attitude, and control surface deflection range. It is difficult for traditional mathematical modeling methods to give an accurate description of all the relationships with acceptable errors.
With the rapid development of computer technology and the integration of interdisciplinary ideas, many other mathematical methods have also been introduced into the field of unsteady aerodynamic modeling, such as the autoregressive model [10][11][12], Volterra series model [13][14][15], fuzzy logic method [16], neural network method [17][18][19][20][21], and support vector machine method [22][23][24], collectively referred to as the intelligent algorithm. The unsteady aerodynamic model based on these algorithms does not pay attention to the physical principles; that is to say, the relationship between the aerodynamic force and the input state quantity is regarded as a "black box." The aerodynamic model is established by a large number of high-quality sample data, so this kind of model is also called the data-driven model.
Although these methods can give accurate prediction within the range of sample flight states, the ability of extrapolation is usually very poor. At present, the research of aerodynamic modeling based on the intelligent algorithm is mostly based on the results of the wind tunnel dynamic test or CFD calculation. However, the model based on wind tunnel and flight test is difficult to reflect the aerodynamic characteristics of all the flight states or motion forms we are interested in, and the high-fidelity CFD calculation may take a lot of time to obtain enough data, so it is difficult to directly use these aerodynamic models in flight simulation. In addition, when the parameters in the aerodynamic system increase, the number of sample flight states and the calculation requirement for modeling also increase, which makes many parameters cannot be identified. Compared with the traditional methods, these intelligent algorithm models are based on "sample data" and do not care about the flow characteristics, so they do not have a physical background and cannot provide more information for the subsequent application. Therefore, it is necessary to improve the intelligent algorithm to make it suitable for practical application.
To sum up, the aerodynamic model based on the traditional mathematical method needs full understanding of the flow characteristics. The stability of the model can meet the requirements, but the modeling process is very complex. Although the aerodynamic model based on the intelligent algorithm does not rely on physical equations, the generalization ability of the intelligent model is often unsatisfactory. Both of the two methods have their own shortcomings, which makes them difficult to be applied to engineering practice.
To solve this problem, this paper proposes a reduced order model (ROM) for unsteady aerodynamics, which combines the intelligent algorithm with a traditional mathematical model and can consider the aerodynamic nonlinearity at a high angle of attack. Specifically, the intelligent algorithm based on NARX is used to build the structure of nonlinear unsteady aerodynamic ROM. Then, the correction functions based on the aerodynamic physical characteristic assumption of the modified Leishman-Beddoes dynamic stall (LB for short) model [25] are added into the intelligent model, and the parameters of the correction functions are identified by aerodynamic data and intelligent algorithm. Finally, the nonlinear aerodynamic ROM is obtained.
For most of the existing data-driven aerodynamic ROM, the dataset comes from the same data source: numerical simulation, experiment, or flight testing. However, aerodynamic data often comes from different sources in practice, and the cost and accuracy are also different. In order to handle the balance between model accuracy and data generation cost, data-driven modeling driven by data fusion is receiving more and more attention [26]. Here, data fusion refers to the process of combining data and information from multiple sources in order to refine, estimate, and get a better understanding of the data [27], and data-driven modeling based on data fusion mainly refers to the multifidelity modeling method based on aerodynamic data with different fidelity. For aerodynamic data, generally, high-fidelity aerodynamic data are obtained through the flight test, wind tunnel test, or direct numerical simulation, while low-fidelity data usually comes from numerical simulation and simplified, such as rough discretization, relaxation convergence tolerance, low precision, and omission of physics [28]. For data with different fidelity, the cost of data acquisition becomes higher with the increase of data fidelity, and due to the limitation of cost, the amount of available data decreases with the increase of fidelity. In addition, it is difficult for high-fidelity data to include all interesting information, which makes the aerodynamic modeling based on these data expensive and unsatisfactory. For example, most experimental data come from harmonic motion rather than random motion. That is 2 International Journal of Aerospace Engineering because random motion is vulnerable to external noise interference, and the data quality of harmonic motion is better. However, the dynamic information contained in harmonic motion is limited. Based on the above problem, the multifidelity method has been developed and widely used, which provides a general framework for fusing data with different fidelity in the modeling process. The basic idea of the commonly used multifidelity framework is to use low-fidelity data or models (physically driven or data-driven) to provide the dynamic trend, while high-fidelity data is used to correct the model to better reproduce the high-fidelity results. In conclusion, the multifidelity model can obtain high-fidelity results with less cost. Multifidelity methods have been used in many fields, such as mechanism analysis, optimization design [29][30][31], statistical inference, and uncertainty quantification [32]. In the field of fluid dynamics, some scholars have applied the multifidelity method to flapping wing dynamics analysis [33,34], aerodynamic optimization [35][36][37][38], flight simulation [28,39,40], hypersonic aerodynamic load prediction [41], low-fidelity turbulence model correction [42,43], and uncertainty quantification of fluid dynamics system [44][45][46][47]. It is worth noting that most of these studies focus on the evaluation of steady-state aerodynamics to achieve rapid simulation and optimization, while there are few studies on unsteady multifidelity aerodynamic modeling. Because the generation of unsteady aerodynamic data is much more expensive than steady aerodynamic data, it is particularly important to explore the effectiveness of the unsteady aerodynamic multifidelity method. Ghoreyshi et al. [48] firstly showed that using both Euler and RANS data for modeling can improve the prediction results of aerodynamic ROM, and the results are better than that of the model based on a single data source. High-fidelity, lowfidelity, and multifidelity ROM are built by the NARX framework. The multifidelity ROM is constructed by adding low-fidelity time delay prediction as additional inputs of high-fidelity ROM. After that, Kou and Zhang [49] took the cokriging concept as the underlying theory and systematically introduced the problem of unsteady aerodynamic data fusion modeling. They built the multifidelity aerodynamic ROM with improved stability of NACA0012 airfoil under transonic conditions by using the NARX framework. The correction term of the model is constructed by the multicore neural network. Euler and unsteady RANS solvers are used to provide low-fidelity and high-fidelity data, respectively. Results show that only three groups of high-fidelity training cases of harmonic motion are needed to establish a high-fidelity model and obtain high-fidelity results.
In this paper, in order to improve the modeling efficiency of the ROM, the multifidelity data fusion modeling method is also introduced, using a small amount of highfidelity data and a large amount of low-fidelity data to build the high-fidelity model. Using data fusion modeling is aimed at improving the efficiency of data acquisition as the cost of high-fidelity data is usually pretty high, which can also improve the modeling efficiency at the same time. In this paper, three kinds of data with different fidelity are used for fusion modeling. The low-fidelity data is obtained by the CFD solver based on Euler. It should be noted that these low-fidelity data can also reflect the real situation under certain restrictions. That is to say, although it is low-fidelity data, it still meets the accuracy requirements. The medium-fidelity data is obtained by the CFD solver based on RANS. The high-fidelity data comes from the existing experimental data. For the modeling process, firstly, the low-fidelity aerodynamic ROM is established by using lowfidelity data. Then, the nonlinear correction function is constructed based on the flow separation principle [25] and is added to the low-fidelity ROM. After that, the parameters of the nonlinear correction function are preliminarily identified by using the medium-fidelity data and the intelligent parameter identification method. Finally, the high-fidelity data from experiments is used to optimize the correction parameters, and the final high-fidelity nonlinear aerodynamic ROM is obtained. In order to verify the proposed high-fidelity nonlinear aerodynamic ROM, the lift and moment coefficients of NACA0012 airfoil in pitching and plunging motion are predicted and compared with the experimental data and CFD-RANS results. Then, the highfidelity ROM is applied to the two-dimensional aeroelastic system. After that, the convergence of the model is analyzed, and its performance is compared with the traditional modified LB model and the intelligent aerodynamic model based on the RBF neural network.

Multifidelity Framework for Nonlinear Unsteady Aerodynamic ROM Materials and Methods
In this part, the multifidelity framework for nonlinear unsteady aerodynamic ROM based on traditional and intelligent algorithm fusion is introduced. It can be divided into the following parts. First of all, it is necessary to select a suitable multifidelity modeling method to develop the relationship between the low-fidelity model and the high-fidelity model. Secondly, the low-fidelity model is established. Thirdly, correction functions are constructed for the difference between the high-fidelity model and the low-fidelity model, and the parameters are identified. Finally, the nonlinear unsteady multifidelity aerodynamic model suitable for a certain speed range and high angle of attack can be obtained.

Multifidelity Modeling
Method. The multifidelity model, that is, the models with different fidelity and computational efficiency, is combined in a certain way, so that the new model can obtain the accuracy of the high-fidelity model with less computational cost. The low-fidelity model can be divided into three categories, including the simplified model, data fitting model, and projection-based model [32]. Here, we use the data fitting model which is also a kind of datadriven model. To combine the low-fidelity model and the high-fidelity model, the method based on correction is adopted; that is, the low-fidelity model is modified by the data generated by the high-fidelity model. Currently, the main correction methods can be briefly divided into three categories [50]: 3 International Journal of Aerospace Engineering (1) Additive and multiplicative corrections: (2) Comprehensive corrections: (3) Space mapping (input correction): where y HF ðxÞ is the high-fidelity model and y LF ðxÞ is the low-fidelity model. ρðxÞ represents the multiplicative corrector while zðxÞ represents the additive corrector. x is the input vector. In this paper, comprehensive corrections are used. However, there are some differences with the standard form. We divide the input factors into two categories according to their effects on the nonlinear response of the system: x = ½ x 1 ; x 2 . We think x 1 represents the input factor which can cause strong nonlinearity of the system response, x 2 is also a part of input factors, but the nonlinear response of the system caused by x 2 is very small. According to the above assumption, we divide the high-fidelity model and the lowfidelity model into two parts, respectively: where y HF ðxÞ and y LF ðxÞ are the total response of the system calculated by the high-fidelity model and the low-fidelity model, respectively. y 1 HF ðx 1 Þ and y 1 LF ðx 1 Þ are the response of the system caused by x 1 , and y 2 HF ðx 2 Þ and y 2 LF ðx 2 Þ are the response of the system caused by x 2 . Superscripts 1 and 2 are identifiers without mathematical meaning which represent the response part caused by x 1 and x 2 , respectively.
Because the previous assumption that the nonlinear response of the system caused by x 2 is very small, we think that the main difference between the high-fidelity model and the low-fidelity model is the difference between y 1 HF ðx 1 Þ and y 1 LF ðx 1 Þ, so only this part needs to be modified. And there is little difference between y 2 LF ðx 2 Þ and y 2 HF ðx 2 Þ, so we regard y 2 LF ðx 2 Þ = y 2 HF ðx 2 Þ. Then, the relationship between the high-fidelity model and the low-fidelity model can be written as the following equation: where the multiplicative corrector ρðx 1 Þ and the additive corrector zðx 1 Þ are both time-varying functions, which can be identified from high-fidelity data.  [51,52]. For MIMO discrete-time systems, the NARX model can be written as follows:

Reduced Order
where k is the time scale, yðkÞ is the system output at time k, and uðkÞ is the system input at time k. NARX is aimed at establishing the nonlinear relationship between the outputs and the inputs at the current and several previous time steps.
In the current study, the nonlinear aerodynamic ROM for two-dimensional airfoil with pitching and plunging freedom at low speed and high angle of attack is modeled; the twodimensional aeroelastic system is shown in Figure 1. In this case, yðkÞ = ½C l ðkÞ, C m ðkÞ T . C l ðkÞ and C m ðkÞ represent the lift coefficient and pitching moment coefficient at time k, respectively. uðkÞ = ½θðkÞ, _ θðkÞ, € θðkÞ, hðkÞ, _ hðkÞ, € hðkÞ T , θ represents the pitching angle of airfoils, while h represents the plunging displacement. Let αðkÞ = ½θðkÞ, hðkÞ T , uðkÞ can be also written as uðkÞ = ½αðkÞ, _ αðkÞ, € αðkÞ T .

Model
Structure. It is considered that the nonlinearity of the aerodynamic model mainly comes from the nonlinear relationship between the input and output of the static part, while the dynamic part can be represented by a linear system. Therefore, the input variables in formula (7) are divided into two groups, and the input variables related to the static part are separated: It should be noted that when M in x 1 ðkÞ is sufficiently larger, the time delayed in the firstand second-order derivatives in x 2 ðkÞ may not be needed. However, in the modeling process, with the increase of the value of M, the model structure will become more complicated, and the parameters that need to be identified will increase too, which will influence the modeling efficiency and model stability. Therefore, generally, M of the model is not necessarily a large value. At this time, the first-order and second-order derivatives are still needed as supplementary parameters to show the dynamic characteristics of the system. 4 International Journal of Aerospace Engineering According to the multifidelity modeling method, lowfidelity and high-fidelity models can be written as follows: where ρðw, x 1 ðkÞÞ and zðw, x 1 ðkÞÞ are the multiplicative and additive correction functions of the high-fidelity model, respectively, and w is the set of correction parameters to be identified. ϕ LF ðx 1 ðkÞÞ and ϕ LF ðx 2 ðkÞÞ are low-fidelity linear models, and ϕ HF ðx 1 ðkÞÞ is high-fidelity nonlinear models.
In the current study, we can further expand them as the following functions: where j are the parameters to be identified corresponding to different input variables. N and M represent the order of the model, reflecting the unsteady effect of the system. Formula (9) can be written as It should be noted that in order to improve the accuracy and efficiency of identification, the parameters of the low-fidelity model and the high-fidelity model are identified, respectively. The parameters of the low-fidelity model are identified firstly. The parameter set W l of the low-fidelity ROM can be written as follows: And the output of the low-fidelity ROM can be written as where C LF l ðx 1 Þ and C LF l ðx 2 Þ represent the lift coefficient generated by input parameters x 1 and x 2 , respectively, while C LF m ðx 1 Þ and C LF m ðx 2 Þ represent the pitching moment coefficient generated by input parameters x 1 and x 2 , respectively.

Nonlinear Correction Model Based on Flow Separation.
For the multifidelity model, some scholars use a neural network to model the nonlinear correction term [49]. However, data-driven models generally have the problems of easy overfitting, poor generalization ability, and poor robustness, and a large number of high-precision calculation data are required to ensure the accuracy and generalization ability of the model, which reduces the modeling efficiency, so it is difficult to be applied to engineering practice. The traditional mathematical expression modeling method can avoid the above problems, because it is based on the physical characteristics of the flow field. Therefore, according to the physical characteristics of the flow field, the nonlinear aerodynamic correction terms ρ ðw, x 1 ðkÞÞ and zðw, x 1 ðkÞÞ at low speed and high angle of attack are modeled.
For the case of a high angle of attack, its nonlinearity is mainly caused by the separation of airflow. The influence of flow separation on airfoil aerodynamics can be divided into two aspects. The first is that the flow separation can directly change the aerodynamics, and the second is that flow separation causes the generation of vortex which leads to additional aerodynamics. ρðw, x 1 ðkÞÞ and zðw, x 1 ðkÞÞ correspond to the two aspects, respectively. In addition, when the pitching motion with a high angle of attack is carried out at a low Mach number, the separation and reattachment between the flow and wing surface will be delayed, which will bring additional aerodynamic force called "overshoot" [25,53]. Here, the influence of "overshoot" is represented by z lma ðw, x 1 ðkÞÞ. To sum up, for the high-fidelity ROM model at a higher Mach number (generally Ma > 0:3 and Ma < 0:8), it can be written as formula (11). For the highfidelity ROM model at a lower Mach number (generally Ma < 0:3), it can be written as Figure 1: Two-dimensional aeroelastic system.

International Journal of Aerospace Engineering
ρðw, x 1 ðkÞÞ, zðw, x 1 ðkÞÞ, and z lma ðw, x 1 ðkÞÞ can be further expanded according to the outputC l ðkÞ, C m ðkÞ: (1) Modeling for ρðw, x 1 ðkÞÞ. Some scholars have found that the aerodynamic correction term directly caused by flow separation can be written as a function of the separation point [25]. The form of the function is shown in formula (16). Here, we directly use the form of this function, but the parameters in the function are not fixed and need to be identified according to the actual situation: where P 1~P4 are the parameters to be identified related to the airfoil section. f ″ ðkÞ is the separation point at time k, which can be written as the following iterative form: where V is the velocity of the incoming flow, dt is the time step, and c is the chord length of the airfoil. T 1 is the time constant of the dynamic separation point and is one of the time parameters to be identified. f ′ ðkÞ can be written in the following forms: where P 5~P7 are the parameters to be identified according to different airfoils, and P 6 corresponds to the critical separation angle of attack of airfoils. F 1 , F 2 , and F 3 are the parameters to be identified related to the separation point. α f ðkÞ is the effective angle of attack at time k.
(2) Modeling for zðw, x 1 ðkÞÞ. The correction term zðw, x 1 ðk ÞÞ of the additional lift and moment generated by the separation vortex can also be written as a function of the separation point, which is also affected by the time parameter T 2 . The additional lift z cl ðw, x 1 ðkÞÞ generated by the separated vortex can be calculated by the following equation: where τ v is the time parameter of separation vortex failure. T 3 is the time constant of the vortex moving on the whole airfoil.
The correction term z cm ðw, x 1 ðkÞÞ for the pitching moment coefficient generated by vortex can be given by the following formula: where P 8 is the parameter to be identified.
(3) Modeling for z lma ðw, x 1 ðkÞÞ. For the special "overshoot" lift at a low Mach number, the lift correction factor z lma cl ðw , x 1 ðkÞÞ of this part is also the function of the separation point, which can be written as International Journal of Aerospace Engineering where P 9 is the parameter to be identified related to the airfoil, α min is the initial angle at which the airflow readheres to the airfoil when the angle of attack decreases, V x and V xr are the shape functions of the lift coefficient generated by the "overshoot." The expression is as follows [54]: The additional moment coefficient z lma cm ðw, x 1 ðkÞÞ on the wing surface caused by the overshoot effect can be calculated by the following formula: where P 10 and T 4 are the parameters to be identified. To sum up, the nonlinear correction parameter set w can be written as w = ½P 1 , P 2 , ⋯, P 10 , F 1 , F 2 , F 3 , T 1 , T 2 , ⋯, T 4 T , and the parameters to be identified can be divided into three parts: airfoil parameter P, separation point equation parameter F, and separation-related time coefficient T. Let Then, the parameter set to be identified for high-fidelity ROM is W n .
In conclusion, the construction of the multi-fidelity model can be divided into two steps. Firstly, identify W l by low-fidelity data and then obtain the parameter set W n of the correction equation based on high-fidelity data. In addition, it should be also noted that when identifying W l and W n , airfoil chord length c and flow velocity V are used for dimensionless processing. Therefore, chord length c and flow velocity V are also used as additional input parameters of ROM, which broadens the application scope of the model. The parameter identification algorithm can be used to identify the above parameters. The next section will introduce the parameter identification algorithm in detail.

Parameter Identification Methodology.
In this paper, a parameter optimization algorithm based on cuckoo search (CS) is used for parameter identification. CS, proposed by Yang and Deb in 2009 [55], is a typical swarm intelligence optimization algorithm with iterative search characteristics. The CS algorithm originated from egg parasitism of cuckoo species, and it simplified the biological behavior into three idealized hypotheses: (I) each cuckoo lays one egg at a time and places it in randomly selected host nests, (II) high-quality eggs and the best nest will be passed on to the next generation, and (III) the number of available host nests is fixed, and the probability that host birds may find foreign eggs is P ∈ ð0, 1Þ.
Each egg is represented as a solution, and the last hypothesis indicates that some host nests will be replaced by new random nests. For the maximization problem, the quality or fitness of the solution can be simply proportional to the value of the objective function. Other forms of fitness can be defined in a way similar to fitness function in a genetic algorithm. In terms of the search algorithm, the CS algorithm does not search data through a simple isotropic random walk but updates data through the combination of local random walk and global random walk. Controlled by switch parameters p a , local random walk can be written as follows: where x t j and x t k are two solutions generated by random walk, H is the heaviside function, ε is the random number with uniform distribution, and s is the step size. α is the step size scaling factor. The global random walk adopts the Levy flight mode, which can be written as follows: where Lðs, λÞ can be written as However, the above algorithm uses a random number to initialize the position of the bird's nest, which is blindness. The position of the random solution may be far away from the correct solution. Meanwhile, the data search and update mechanism of Levy flight reduces the convergence speed in the iterative process, and it may not be able to find the solution with the required accuracy in a certain period of time. In order to solve this problem, chaos optimization is introduced. The randomness and ergodicity of chaos can avoid the problem of stopping searching when the local minimum is found, which can overcome the shortcomings of traditional optimization algorithms. In this paper, the logistic map, which is commonly used in chaos optimization, is introduced into the solution updating process of the CS algorithm: where u ∈ ð2, 4 is the chaos control variable. Before chaos optimization, the variables to be optimized need to be mapped to ½0, 1.
The cuckoo search algorithm with chaos optimization can be written in the following steps: (1) Define the number of initial nests N and the initial solution x 0i ði = 1, 2, ⋯, NÞ (4) Project x best1 to ½0, 1 and obtain the chaotic sequence y chaom by Logistic mapping. Then inversely map the chaotic sequence to the original solution space to obtain a new solution x chaom and calculate the fitness of the new solution. Repeat the above process several times and find the optimal solution in x chaom ðm = 1, 2, ⋯, MÞ as x best2 . By comparing x best1 with x best2 , the final best solution of this iteration x best is obtained (5) Discard a certain number of nests by probability P and replace them with new nests. Retain the optimal solution x best after the replacement (6) Judge whether the optimal solution x best meets the accuracy requirements. If it meets the requirements, stop the iteration process and output x best . Otherwise, return to step 2 The whole process can also be illustrated by the flow chart shown in Figure 2.
It should be noted that for parameters with physical meaning, the optimization range of these parameters will be given according to the conclusions of some literatures such as Ref. [25,54] to meet the physical meaning.

Model Training and Validation
The NACA0012 airfoil is taken as an example for model training and validation. Firstly, different fidelity data sources are introduced; then, the model training process is represented, and finally, the model is verified in aerodynamic prediction and aeroelastic analysis.
3.1. High-Fidelity, Medium-Fidelity, and Low-Fidelity Data. Theoretically, establishing a multifidelity model only requires two kinds of datasets with low and high fidelity and two steps. For the first step, the parameters of the lowfidelity model are identified through the low-fidelity data. For the second step, the parameters of the high-fidelity model are obtained according to the high-fidelity data. However, for the multifidelity aerodynamic model in this paper, the high-fidelity data obtained by experiments is limited, and the airfoil motion corresponding to the experimental data is usually harmonic motion, which makes the dynamic   International Journal of Aerospace Engineering information of the system in high-fidelity data insufficient.
In addition, there are many parameters to be identified in the high-fidelity model; directly using experimental data to identify the nonlinear parameters may cause overfitting. The generalization ability of the model will be very poor. So, it is not suitable to use limited experimental data alone for nonlinear modeling. Therefore, it is necessary to obtain more high-fidelity aerodynamic data with more dynamic information of the system for modeling. The CFD-RANS model is a good choice to obtain more aerodynamic data with relatively high fidelity. On the one hand, the CFD results based on the RANS solver can reflect the nonlinear aerodynamic characteristics of a high angle of attack to a certain extent. On the other hand, the airfoil motion form can be given randomly according to the actual requirements, which is not limited by the experimental conditions. As a result, CFD data based on the RANS model is chosen as the third kind of data with relatively high fidelity (here defined as medium fidelity) for multifidelity aerodynamic model modeling, which can be used for the preliminary identification of nonlinear high-fidelity model parameters.
For the multifidelity aerodynamic model established in this paper, the modeling process is divided into three stages. Firstly, the low-fidelity model is identified through the lowfidelity data; then, the high-fidelity model parameters are identified for the first time through the medium-fidelity data, and finally, based on the results of the first identification, the high-fidelity model parameters are modified for the second time through the high-fidelity data from experiments. As just described, in this study, three kinds of data with different fidelity are used in modeling, which are called highfidelity data, medium-fidelity data, and low-fidelity data, respectively. Here, low-fidelity data is obtained by CFD-Euler solvers. Medium-fidelity data is obtained by the CFD-RANS solver. High-fidelity data is obtained by experiments. It should be noted that the final model will be modified by high-fidelity experimental data, so the CFD-RANS model only needs to reflect the trend of aerodynamic force at a high angle of attack. The result of the RANS solver does not need to be in perfect agreement with the experimental data, which further saves the time to adjust the parameters   9 International Journal of Aerospace Engineering and mesh quality of the RANS solver, as well as improve the modeling efficiency. Therefore, the aerodynamic force calculated by RANS is taken as the relatively high-fidelity data and here is defined as medium-fidelity data.
In this study, the ALE (Arbitrary Lagrange Euler) finite volume method is used as the governing equation. In this way, the motion of the material structure boundary can be effectively tracked in the calculation process, and the grid position can be appropriately adjusted according to the defined parameters in the calculation process without serious distortion. The double time step implicit time discretization method is used for time discretization, which can

10
International Journal of Aerospace Engineering improve computational efficiency. For the turbulence model, SST k-ω is selected, which can be applied to the calculation of flow separation and reverse pressure gradient. The scheme of the computational grid deformation is based on diffusion, and the aerodynamic grid for CFD-Euler and CFD-RANS is shown in Figure 3. For the mesh for the RANS solver, the grid height of the first boundary layer is 10 -3 mm, and the growth rate is 1.2.

Model
Training. The training data of the unsteady aerodynamic model are usually obtained by two methods. One is to design a training signal with specified excitation (such as random signal, 3211 signal, and chrip signal), which needs to contain the dynamic information of the main motion state in the actual flight situation. The other is to select some representative states in the sampling space, where each state corresponds to a specific motion type, such as harmonic motion with constant frequency and amplitude [49]. Both of the two methods can be easily applied to numerical simulation, but the first method cannot be applied to the experiment due to the noise. For the nonlinear unsteady aerodynamic ROM established in this paper, on the one hand, it needs a large amount of data with a wide range of characteristics to meet the generalization ability requirement; on the other hand, it needs experimental data for secondary correction, so both of the two methods are used to obtain the training data.
As described in Section 3.1, the model training process can be divided into three stages. The first stage is the lowfidelity linear ROM training. The Euler solver is used to obtain the aerodynamic training data. It should be noted that for the Euler solver which ignores the viscosity, the motion form corresponding to the training signal should 13 International Journal of Aerospace Engineering be limited in a certain range. Otherwise, due to the influence of viscous force, the obtained aerodynamic coefficients cannot meet the actual situation and cannot be used in aerodynamic modeling. The RANS solver does not have this problem because it considers the influence of viscosity. However, whether the training data is obtained by the Euler solver or RANS solver, it should meet the following requirement to achieve a better training result; that is, the training signal needs to contain a wider range of angle of attack (initial angle of attack is also included) and reduced frequency (the reduced frequency is defined as k = ðω × cÞ/2V, where ω is the angular frequency of pitching and plunging motion). Therefore, several groups of chirp signals with variable amplitude as well as initial value are selected as training signals, one of which is shown in Figure 4. The frequency range of the signal is 0-5 Hz, and the amplitude range is set according to the actual situation; for the Euler solver, the amplitude of the effective angle of attack (defined as α e = θ + arctan ð _ h /VÞ) is limited to ½−5°, 5°. Each group of training signals also corresponds to multiple working conditions with 14 International Journal of Aerospace Engineering different incoming velocities, which can further widen the reduced frequency range. The aerodynamic training data corresponding to the signal in Figure 4 is shown in Figure 5. The second stage is the first time training of the highfidelity nonlinear ROM, which is also the initial identification of nonlinear correction function parameters. In this stage, the RANS solver is used to obtain the training data. The training signals are the same as described above which are also the chirp signals with variable amplitude, but the effective angle of attack changes from -30°to 30°. CFD data based on the RANS model is used for the preliminary identification of nonlinear correction parameters, and one or more groups of nonlinear correction parameters with stronger generalization ability can be obtained. A group of training data is shown in Figure 6. It can be seen that under the same inflow velocity, the nonlinearity of aerodynamic force is stronger at the motion with high amplitude and low frequency. In addition, the aerodynamic nonlinearity becomes more significant with the increase of the velocity, which may be the result of the decrease of the reduced frequency.
The third stage is the second time training of the highfidelity nonlinear ROM, which is also a second time of optimization for nonlinear correction function. Based on the nonlinear parameters identified in the second stage, the     [25], and Leishman's [57] paper is used as the training data, and the cuckoo search algorithm with chaos optimization is used for the optimization. Finally, the high-fidelity nonlinear ROM with high accuracy and strong generalization ability is obtained. For low-fidelity and medium-fidelity training datasets, six velocity points are selected as the training velocity, which are 10 m/s, 20 m/s, 40 m/s, 60 m/s, 80 m/s, and 100 m/s, respectively. There are one Euler training dataset and two RANS training datasets at each training speed. The training signal of each dataset is the chirp signal with variable amplitude, and the maximum amplitude of the signal is 5°, 20°, and 30°, respectively. The unit time step is 0.001 s, and the number of time points of each training set is 5000. The range of the Reynolds number of the training datasets is 4 * 10 5 5 * 10 6 . For high-fidelity experimental training datasets, three groups of experimental data are used for modeling, and the parameters contained in the experimental data are listed in Table 1, corresponding to Case 3, Case 5, and Case 6. The reason for selecting Case 3, Case 5, and Case 6 as the training dataset is that the high-fidelity model needs to reflect the strong nonlinear characteristics of aerodynamic force. Therefore, three sets of data with stronger nonlinearity are selected as the training datasets. The dataset of Case 4 with relatively weak nonlinearity is used as the test dataset to test the generalization ability of the model.
Through the above training progress, the parameters of the low-fidelity model and high-fidelity model can be identified. For the low-fidelity model, the delay orders M and N must be determined first. Here, M = N = 2 is determined by a trial-and-error method. After that, the parameters in Equations (12) and (26) can be identified by the cuckoo search algorithm with chaos optimization. In total, there are 39 parameters needed to be identified by multifidelity aerodynamic data. After the identification of the parameters, the low-fidelity and high-fidelity models are established.

Model Validation
3.3.1. Validation in Aerodynamic Prediction. By estimating the aerodynamic loads of NACA0012 airfoil under different pitching and plunging motions at low speed, the effectiveness of the unsteady nonlinear multifidelity modeling method based on data fusion is verified.
Firstly, several simple harmonic pitching motions are selected as test cases. The information and the working parameters of the test cases are listed in Table 1. Among them, case3-6 have experimental data and case1-2 are the test cases without experimental data. Use RMSE (root mean squared error) to measure the performance of different aerodynamic models. The calculation formula of RMSE is as follows: where i is the time scale and N means the number of time points of a group of test data. y H and y L represent the aerodynamic coefficients obtained by high-fidelity models/data and low-fidelity models/data, respectively. Table 2 lists the RMSE of test case1-2 compared with the CFD-RANS result. Table 3 lists the RMSE of test case3-6, which is obtained by comparing with the experimental data. Figures 7 and 8 show the comparison of aerodynamic coefficients predicted by different aerodynamic models.
It can be seen from Figure 7 that for linear test Case 1 and Case 2 which pitching angle changes from [-10°, 10°] and reduced frequency is 0.46 and 0.23, respectively, the shape of the hysteresis loop is almost elliptical, and all models can effectively capture the unsteady characteristics of CL (lift coefficient) and CM (moment coefficient), which verifies the low-fidelity ROM and the high-fidelity ROM in the linear range. In addition, it can be seen from Table 2 and Figure 7 that the results calculated by the low-fidelity model and CFD-Euler model are almost the same. And Figures 7(b) and 7(d) show that the curve of the lowfidelity model is closer to the result of CFD-Euler, while the curve of the high-fidelity model corrected by CFD-RANS data is closer to the result of CFD-RANS, which also proves the effectiveness of the nonlinear correction model to a certain extent.
For Cases 3-6, the aerodynamic nonlinearity is very strong, and the low-fidelity model cannot meet the accuracy requirements. CFD-RANS can reflect the trend of nonlinear aerodynamics to a certain extent, but they are not accurate enough to fit the details. The proposed high-fidelity ROM improves this defect, as its parameters are optimized based on experimental data. Results show that the prediction results of the proposed high-fidelity ROM have the best consistency with the experimental data, which can be seen from the fitting degree of the hysteresis loop in Figure 8 as well as the RMSE data in Table 1. In addition, we also found from Table 1 and Figure 8 that with the increase of the initial angle of attack, the aerodynamic nonlinear characteristics will become more significant, and the consistency between CFD results and experimental data will become worse. This

16
International Journal of Aerospace Engineering may imply that when the initial angle of attack or total motion amplitude increases to a certain value, the computational accuracy of CFD will no longer meet the requirements. The proposed high-fidelity ROM based on multifidelity data can avoid the above problem to a certain degree, because it can modify itself based on the experimental data, which makes it more extensible for different conditions. In order to prove the generalization ability of the model, more random test signals are needed. We choose a group of fwgn signals with random parameters as the test signal. The motion of the test case is a combination of pitching motion and plunging motion. The test signal is shown in Figure 9. Figure 10 shows the comparison among the output results from different models and CFD solvers. It can be seen from Figure 10 that the established high-fidelity ROM well reflects the nonlinear aerodynamic characteristics, which proves its good generalization and extrapolation ability.
As shown in the test case, the proposed aerodynamic modeling method only needs a group of low-fidelity and medium-fidelity CFD data and a few high-fidelity data with a specified harmonic motion to establish high-fidelity aerodynamic ROM with satisfactory generalization ability.

Verification in Two-Dimensional Aeroelastic Response.
In order to further verify the validity of the model, the nonlinear ROM is applied to two-dimensional aeroelastic analysis, and the results are compared with the experimental results in Ref. [58].
The two-dimensional aeroelastic system is already shown in Figure 1, which has two degrees of freedom: pitching motion and plunging motion. The important experiment parameters in Ref. [58] are listed in Table 3.
The dynamic equation of the system can be written as follows: where θ is the pitch angle and h is the amplitude of plunge motion. K h and K θ are spring stiffness of plunging and pitching freedom, respectively. C h and C θ represent structural damping. m is the mass per unit length. S θ represents the mass static moment, and I θ is the moment of inertia at quarter chord length. Q h and Q θ are external aerodynamic forces and aerodynamic moments, which can be calculated by the proposed nonlinear aerodynamic ROM with the following equation: where C L and C M are aerodynamic coefficient and aerodynamic moment coefficient, respectively, which are calculated by nonlinear aerodynamic ROM. c is chord length and l is  The Runge-Kutta method is used to solve Equation (32) in the time domain, and the bifurcation and limit cycle analysis can be carried out, as shown in Figures 11 and 12. The x -axis parameter α in Figure 12 can be expressed by pitching angle θ and plunging velocity _ h : α = θ + arctan ð _ h/VÞ. It can be seen from Figures 11 and 12 that the results of two-dimensional aeroelastic bifurcation and limit cycle analysis based on nonlinear aerodynamic ROM proposed in this paper are in good agreement with the experimental results, which further verifies the effectiveness of the aerodynamic reduced order model established in this paper.

Model Analysis and Comparison
In the previous section, the proposed high-fidelity aerodynamic ROM based on the multifidelity framework is verified by the accurate prediction of the nonlinear unsteady aerodynamic coefficients as well as the simulation of twodimensional limit cycle response induced by nonlinear aerodynamics. However, the performance of the model still needs to be investigated. This section studies the performance of the model from three aspects, respectively. Firstly, the convergence of model error with the velocity interval of training samples is studied. Then, the calculation accuracy of the proposed model is compared with the modified dynamic stall model (also called the modified LB model proposed in Ref. [25]) and nonlinear aerodynamic model based on RBF neural network. After that, the application scope and generalization ability of the proposed ROM model is explored and compared with other aerodynamic models.
4.1. Convergence Analysis. As described in Section 3.2, the modeling process of the nonlinear ROM can be divided into three stages. The first stage is to identify the low-fidelity model, and the second and the third stages are to identify the high-fidelity correction equation. However, due to the limited training data in the third stage, only the convergence of the first two stages is studied. The study of the convergence is divided into two parts also, according to different training stages.
The study is conducted as follows. Firstly, we need to define the velocity range of the training model. In this study, the velocity range is 10-60 m/s. Then, the training cases should be grouped according to different velocity intervals. Here, four training groups with different velocity intervals are set up; the corresponding speed intervals are 50 m/s, 25 m/s, 10 m/s, and 5 m/s. The velocity included in each training group is shown in Table 4. After that, the training data of each group are calculated, and the different ROM models are obtained after training. In order to eliminate the influence of signal forms, the training signals for each velocity point are the same, as shown in Figure 4. Finally, a set of test cases are selected, and the aerodynamic coefficients are predicted by different ROM models. The signal    Figure 13, and the velocities of the test cases are listed in Table 5.
As can be seen from Figure 14, for the low-fidelity ROM, the RMSE of different test cases all decreases with the decreasing of velocity interval. However, it should be noted that reducing the velocity interval means more training data need to be obtained, which will increase the time cost for modeling. So, it is necessary to balance the accuracy and efficiency. For the case of this paper, the change of the RMSE is very limited when the interval of training velocity points is less than 10 m/s. Therefore, it is the best to select training group 3 as the training case for low-fidelity ROM.
For the high-fidelity ROM, as shown in Figure 15, RMSE of different test cases change differently with the decreasing of training velocity interval. For test Cases 3~5 where inlet velocity is less than 30 m/s, RMSE increases first and then reduces, while for other test cases with higher inlet velocity, the opposite is true. This is because the influence of gas

20
International Journal of Aerospace Engineering viscosity on aerodynamic force is different from the change in flow velocity. When the flow velocity is pretty low, the viscous effect of the gas will be stronger and the influence on the aerodynamic force will be greater. Therefore, for the same pitching motion at 10 m/s and 60 m/s, there are some differences in aerodynamic characteristics. If the training data contains more cases at a higher speed, the aerodynamic characteristics reflected by the training data are more consistent with the aerodynamic situation at a higher speed. Therefore, for test Cases 1 and 2, the result of training group 2 is the best, as it contains only one case whose inlet velocity is less than 30 m/s, while the other two cases are both at a higher speed. After reducing the velocity interval and increasing lower incoming flow speed as training points, the model is adjusted to fit the aerodynamic characteristics at a lower velocity, so that the fitting degree at 30-60 m/s is reduced. However, compared with the improvement of model generalization ability, the loss of accuracy is acceptable. In addition, it can be seen from Figure 15(a) that for the lift coefficient, RMSE increases slightly when the velocity interval decreases to 5 m/s. This may be due to the overfitting of the model, which should be avoided as far as possible. Therefore, for the high-fidelity ROM, the best training velocity interval is 10 m/s, corresponding to training group 3.     Figure 16 and Table 6. It can be seen from Figure 16 and Table 6 that the aerodynamic ROM based on the RBF neural network has the best fitting degree for the test case. This is because the RBF neural network has the ability to fit any nonlinear function, and the test cases are also in the training set, so the results of the RBF model are the best. However, the generalization ability of the RBF model needs to be further verified because

23
International Journal of Aerospace Engineering of the limited training data, which will be explained in the next subsection.
The fitting degree of the proposed algorithm of test cases is basically equivalent to that of the traditional modified LB algorithm, which can be seen from Figure 16 and Table 6. However, this traditional modified LB model is based on a large amount of experimental data, and parameter determination is based on subjective experience, which costs a lot of time and brings a lot of uncertainty for modeling. Moreover, after changing the experimental parameters and airfoil conditions, it is necessary to redesign the experiment and obtain the data, so the modeling process of the LB model will cost a pretty long time. The aerodynamic ROM proposed in this paper only needs a group of CFD data with low-fidelity or medium-fidelity and few experimental data for modeling, which improves the modeling efficiency.

Generalization Ability and Application Scope Analysis and Comparison.
In this section, the generalization ability and the application scope of the proposed high-fidelity aerodynamic ROM are explored and compared with the RBF model. Firstly, the training data are obtained in a certain speed range (10-60 m/s), and the unsteady aerodynamic ROM is established by the RBF neural network algorithm and the algorithm proposed in this paper, respectively. The signal of training data is the same as that in Figure 4. Five velocities (10 m/s, 20 m/s, 30 m/s, 40 m/s, 50 m/s, and 60 m/s) are selected as training points. Then, the two kinds of ROM are applied to the test cases beyond the training speed range, and the prediction results are compared with the CFD-RANS calculation results. There are two main reasons for taking CFD-RANS calculation results as the reference. Firstly, the experimental data at a higher speed is insufficient. CFD-RANS data is a good supplement as the relatively high-fidelity data, and it has been proved in the article that the CFD-RANS model can reflect the trend of nonlinear aerodynamic force to a certain extent. And the generalization ability can be also seen as the ability to follow the response of the aerodynamic trend under different conditions. Therefore, we think it is reasonable to take the CFD-RANS data as the reference, because CFD-RANS data can reflect the trend of aerodynamic force. Another reason for taking CFD data as a reference is that CFD data is the easiest to obtain. In order to more clearly reflect the variation of aerodynamic nonlinear characteristics with different incoming flow velocities, the chirp signal with variable amplitude is also used as the test case. Figures 17 and 18 show the comparison of the aerodynamic coefficients predicted by the CFD-RANS solver and the two kinds of ROM under different inlet velocities of 100-250 m/s. It can be seen that the aerodynamic coefficients predicted by high-fidelity ROM proposed in this paper are in good agreement with the CFD-RANS results, which shows that the model has a strong generalization ability in the speed range of 100-250 m/s. However, the fitting degree of the RBF model to CFD data decreases rapidly with the increase of airspeed, which indicates that the generalization ability of the RBF model is really poor, and the RBF model needs more data in modeling. Therefore, the modeling effi-ciency of the proposed ROM is higher than that of the aerodynamic model based on the RBF neural network.
It is worth noting in Figures 17 and 18 that with the increase of the incoming flow velocity, the fitting degree of ROM proposed in this paper will become lower, which means that the ROM proposed in this paper may no longer be applicable when the incoming flow velocity is higher than 250 m/s. This is because when the velocity of the incoming flow continues to grow, shock waves will be generated locally on the wing surface, which will affect the aerodynamic characteristics. The proposed aerodynamic models cannot reflect the effect of shock wave, which means it is no longer applicable. Moreover, when the inflow velocity is less than 10 m/s, the adaptability of the ROM model also becomes worse due to the enhanced viscous effect of the gas. In general, the prediction accuracy of the lift coefficient is higher than that of the moment coefficient, because the nonlinearity of the moment coefficient is stronger. Based on the above analysis, we can draw the conclusion that the ROM based on algorithm fusion and multifidelity framework established in this paper is suitable for the airspeed condition at 0.03~0.8 Ma.

Conclusions
A multifidelity modeling framework for low-speed nonlinear unsteady aerodynamic ROM is presented in this paper, which combines traditional and intelligent aerodynamic modeling methods. The intelligent algorithm based on NARX is used to build the nonlinear aerodynamic model framework, and the traditional LB separation equation based on the flow characteristics is used to build the modified equation. The framework integrates data with different fidelity and extends the multifidelity method of nonlinear system identification to simulate the unsteady effects of nonlinear aerodynamics. The low-fidelity data and the mediumfidelity data are obtained by the CFD-Euler solver and CFD-RANS solver, respectively, and the high-fidelity data is obtained from the experimental results. The main conclusions drawn from the results can be summarized as follows: (1) The proposed high-fidelity ROM performs well in aerodynamic coefficient prediction and twodimensional nonlinear aeroelastic analysis of NACA0012 airfoil at low speed and high angle of attack, which indicates that the proposed aerodynamic modeling method only needs a group of low-fidelity or medium-fidelity CFD data and few high-fidelity data with a specified harmonic motion to establish the high-fidelity ROM with satisfactory generalization ability. The results are in good agreement with the experimental data (2) Compared with the traditional modified LB aerodynamic model, the proposed high-fidelity ROM is more flexible and more adaptable (3) Compared with the aerodynamic model based on the RBF neural network, the proposed high-fidelity 24 International Journal of Aerospace Engineering ROM has stronger generalization ability, higher stability, and modeling efficiency (4) Data fusion can further reduce the modeling cost and improve modeling efficiency. Thus, the efficiency of calculation and simulation is greatly improved compared with high-precision CFD and experiments (5) Future work includes improving the modeling and identification algorithm to better represent the correction term and involving it in three-dimensional configurations

Data Availability
The numerical simulation result data used to support the findings of this study are available from the corresponding author upon request.

Disclosure
This paper is a part of a research project supported by the Department of Aeroelasticity, School of Aeronautic Science and Engineering, Beihang University.