An effort saving method to establish global aerodynamic model using CFD

Purpose – The typical approach of modeling the aerodynamics of an aircraft is to develop a complete database through testing or computational ﬂ uid dynamics (CFD). The database will be huge if it has a reasonable resolution and requires an unacceptable CFD effort during the conceptional design. Therefore, this paper aims to reduce the computing effort required via establishing a general aerodynamic model that needs minor parameters. Design/methodology/approach – The model structure was a precon ﬁ gured polynomial model, and the parameters were estimated with a recursive method to further reduce the calculation effort. To uniformly disperse the sample points through each step, a unique recursive sampling method based on a Voronoi diagram was presented. In addition, a multivariate orthogonal function approach was used. Findings – A case study of a ﬂ ying wing aircraft demonstrated that generating a model with acceptable precision (0.01 absolute error or 5% relative error) costs only 1/54 of the cost of creating a database. A series of six degrees of freedom ﬂ ight simulations shows that the model ’ s prediction was accurate. Originality/value – This method proposed a new way to simplify the model and recursive sampling. It is a low-cost way of obtaining high-ﬁ delity models during primary design, allowing for more precise ﬂ ight dynamics analysis.


Introduction
Designers typically propose several layouts of aircraft during the conceptual design stage before settling on one as the final plan. The layout selection matrix is a frequent layout comparison tool (Gudmundsson, 2013). This method scores design aspect of each layout, and then uses weighted summation to choose the layout with the highest total score. The matrix also incorporates aircraft dynamics; the designers' extensive experience is used to analyze the benefits and drawbacks of these aspects from design aspects. However, aircraft design is continuously evolving nowadays, and a huge number of aircrafts with novel layouts have been introduced. Even the most expert designers cannot quantitatively determine the dynamic characteristics of a layout without prior knowledge. The inaccuracy in the estimation of performance, stability and control will eventually affect the selection decision. As a result, even in conceptual design, ensuring accuracy as much as feasible is always good.
Using the classical six degrees of freedom (6-DOF) dynamic model of an aircraft plus an aerodynamic module to supply aerodynamic coefficients is a practical way to build a dynamic model in the early design stage. To accurately describe the aerodynamic of an aircraft, we usually build a database. Each explainable variable, such as angle of attack (AOA), angle of sideslip (AOS), deflection of aerodynamic control surfaces, altitude and speed, has one dimension in the database. The size of this database grows exponentially as the number of coupling explainable variables increases. In certain circumstances, it will reach the tens of thousands of entries.
This database can be created in a variety of ways. Wind tunnel or flight experiments can surely produce highly credible data and had been widely used (Morelli and DeLoach, 2003;Ozger, 2007;Morelli, 2012;Grauer and Morelli, 2014;Abramov et al., 2019), but because only one layout will be chosen, conducting wind tunnel or flight testing for all the layouts is too expensive and time consuming for conceptual design. Computational aerodynamics, such as panel technique and the vortex lattice method, is the most often used method. These approaches are quick and, in most situations, provide reasonably accurate results. However, aircraft layouts are continually evolving, with many new layouts having separated flows. Most computational aerodynamics methods are no longer capable of handling these complex flows and therefore are no longer dependable. Computational fluid dynamics (CFD), a more advanced calculation method, can provide a more precise estimate of aerodynamic for complex flows. But, because the computation is relatively slow, using CFD to create aerodynamic databases is also expensive, even if it saves time compared to wind tunnel experiments. If a technique can be developed to speed up the process of obtaining aerodynamic data via CFD, then the aircraft layout selection will have a more accurate reference.
Aerodynamic modeling has been prevalently used in aircraft design, e.g. USAF Datcom is a kind of empirical based model. Nowadays approximate models, surrogation models and reduced order models have been widely used to reduce the workload (Jegarkandi Mohsen et al., 2009;Ghoreyshi et al., 2012). In general, the number of model parameters is one or two orders of magnitude lower than the number of database entries, and in the case of extreme coupling, there are even fewer. Among these models, the ones that do not consider flow physics are referred as datadriven models (Kou and Zhang, 2021). Three methods are commonly used to determine a data-driven model: system identification, data fusion and feature extraction. System identification is one of the cornerstones of control theory and has been widely used in aerospace (Plaetschke and Weiss, 1988;Hamel and Jategaonkar, 1996;Tischler and Remple, 2012;Jategaonkar, 2015;Morelli and Klein, 2016).The definition of data fusion varies in different science communities, Although different scientific communities have varied definitions for data fusion, it is commonly defined as the process of merging data and information from multiple sources. Data fusion aims to balance the cost and accuracy of data with multifidelity to refine or estimate the data (Murphy et al., 2016;Kou and Zhang, 2021;Shi et al., 2021). Feature extraction establishes models based on flow modes, and it is commonly used to solve unsteady flow problems. Proper orthogonal decomposition (Lumley, 1967) and dynamic mode decomposition (Schmid, 2010) are the most commonly used feature extraction modeling methods. The majority of classical system identification approaches, which include the output error method, the filter error method and the equation error method (Jategaonkar, 2015), can likewise be applied to aerodynamic modeling. Some other modeling method based on input and outputs, like Volterra series and Kriging, also have been extensively applied to modeling many nonlinear systems (Han et al., 2010;Cheng et al., 2017).
With the rise of machine learning, other methodologies have emerged as well. Many classical modeling techniques can be classified as machine learning, as they are datadriven models that generalize beyond the training data . There are three types of learning algorithms: supervised, semisupervised and unsupervised learning. Neural networks are probably among the most well-known supervised learning methods. It has universal approximate capability that is able to fit any smooth function exactly, under certain conditions (Kubat, 1999). Long short-term memory, a variant of recurrent neural network, has proved its superiority in time series forecasting problems by overcoming stability problem (Pascanu et al., 2013). Classification methods are perhaps the earliest methods for supervised learning, and the most prevalent methods are support vector machines and random forests. Semisupervised learning includes generative adversarial networks, which had already shown some success in enhancing the resolution of flow fields (Xie et al., 2018), and reinforcement learning, which had been widely used in replicating unstable flow dynamics (Gazzola et al., 2014;Gazzola et al., 2016;Guéniat et al., 2016;Loucks and van Beek, 2017). As for unsupervised learning, it is just another name for feature extraction modeling .
Obtaining information, determining the model and estimating parameters and validating the model are the three steps in a classical system identification workflow. The problem of generating an aerodynamic database could be changed into the task of determining aerodynamic model parameters if a preconfigured model is used to predict the aircraft's aerodynamics. However, several issues must be resolved before the modeling method can be used in the conceptional design. The model must be precise enough to predict complex behavior induced by separated flow but not overly complex to avoid overfitting. The parameter estimation must be recursive to conserve valuable computing resources, which requires a recursive sampling method to evenly disperse the samples at each step. The research proposed a structure preconfigured polynomial model and a unique recursive sampling strategy based on the Voronoi diagram to overcome the challenges mentioned above. This method can greatly reduce calculation time, making it feasible to evaluate the performance, stability and control characteristics of each alternative configuration quantitatively during the conceptual design stage. Furthermore, this technique provides a low-cost approach to estimate new aircraft performance and handling characteristics, making the use of CFD in the development of novel layout aircraft a feasible task. This is especially beneficial for most researchers and small aviation manufacturing companies with limited high-performance computing capabilities, which will surely promote innovation.

Overview
Because the model is deterministic if the flow is steady and all the explainable variables and aerodynamic coefficients (i. e. inputs and outputs) are always available, the requirements of system identification are naturally satisfied (Zadeh, 1962). Sampling, model structure identification, parameter estimation and error analysis are all steps in creating a full aerodynamic model. The structure of the aerodynamic model is not overly sophisticated for most aircraft, and it can be substituted with a preconfigured generic model framework. As a result, the focus of this article is on presenting this preconfigured aerodynamic model framework as well as the associated parameter identification approach, which includes sampling and parameter estimation.
The following requirements must be met by this model framework: The model should have a limited number of parameters. Estimating more model parameters necessitates more input data, which will inevitably increase CFD effort, which is in direct opposition to our goal.
The model must be complex enough to represent complex flows such as stalls and separated flow on certain areas of the aircraft.
It is preferable for the model to have a precise physical meaning to obtain understanding from the parameter estimation results, which will speed up the design process. Meanwhile, the model can be simply created so that simplification can be based on features like symmetry.
The parameter identification approach must also meet the following requirements: The aerodynamic estimating accuracy must be calculable so that we can verify that the result produced is accurate within the range we are interested in. More crucially, once the model gets the accuracy we require, we can terminate the estimation process.
The procedure should be recursive, and the process should be terminated as soon as the required precision is achieved, saving the most time possible.
Because the aerodynamic of the aircraft are unknown, if the sample point distribution was nonuniform, information will be lost in the sparse sample point area, resulting in a decrease in the model's prediction accuracy. Furthermore, if the estimating approach is recursive, the sampling must be "recursively uniform." If a sample set is almost uniform before adding the next sampling point, it will remain about uniform after adding the following sampling point. To achieve these criteria, a whole new sampling procedure is required.

Aerodynamic coefficients involved
The aerodynamic force coefficients and coordinate system are defined identically to their classic definitions in wind tunnel tests: Àcos a cosb sinb sin a cosb Àcos a sinb cosb Àsin a sinb sina 0 Àcos a (1) The aerodynamic moment coefficients and coordinate system are defined identically to their classic definitions in flight dynamics: where the forces X, Y and Z and the moments L, M, and N are in the body frame.

Establishment and simplification of the general model
As previously mentioned, the model itself must be simple enough to reduce calculation time, which is our goal. At the same time, the model should provide us with some insights, allowing us to simplify the model depending on certain features. As a result, despite its popularity in recent years, machine learning is not exactly appropriate here. First is that most universal artificial intelligence models require a large amount of data to train , especially in deep learning (Taira et al., 2020), resulting in an unacceptably high CFD effort. For example, have the C L of an airfoil accurately predicted takes approximately 400 results to train a convolutional neural network (Zhang et al., 2018). Second argument is that for researchers, an artificial intelligence model is more of a "black box." It can produce good results, but we do not know much about its internals, which makes it difficult to simplify. The last disadvantage is that many machine learning architectures still do not readily incorporate physical constraints in the form of symmetries . Physics based models, for example, surrogate models, may has advantages in many scenes. However, with the rapid development of new layout aircraft, many layouts may have different flow physics in some special working conditions. It is challenging to establish physical-based models compatible with all possible working conditions. It is preferable to use a polynomial model. As a data-driven model, the polynomial model does not contain lots of physical information, which improves the generality of the model. Its form is simple, making it easy to design to represent aircraft aerodynamic properties. For example, an antisymmetric moment coefficient such as C l or C n can be simply represented by an odd function. As a result, a polynomial model was adopted.
Each explanatory variable's order should be sufficient to replicate the aerodynamic yet keeping the number of regressors as low as possible. Fewer regressors reduce workload while also preventing overfitting.
The model could be simplified because most aircraft are laterally symmetrical. C D , C L and C M are symmetric about lateral explainable variables, whereas C Y , C l and C N are antisymmetric, i.e.: Moreover, the aerodynamics of lateral symmetrical aircraft is symmetrical about the combination of two lateral explainable variables, as shown in Figure 1. The flight state is symmetrical about the product of the AOS and the aileron deflection. Because in the polynomial model the symmetry is equivalent to even function, we can further determine the symmetry of the higherorder terms from the sum of the order of the AOS and the aileron deflection. If the sum is odd, this term will be antisymmetric. In symmetrical coefficients such as C D , C L and C M , the antisymmetrical terms are set to zero and vice versa. Theoretically, this reduces the number of parameters by half.

Recursive maximum likelihood estimation method
Because it is difficult to predict how many samples will be needed before the estimation, the sample size often exceeds the necessary if they were prepared ahead of time. The burden can be further decreased if the regression approach uses a recursive strategy. The form of a classic linear regression model is: where y is the aerodynamic coefficient, x is the m Â 1 column vector of regressors, W is the column vector of model parameters and « r is the regression residual with a normal distribution N (0, s 2 ) according to the central limit theorem. The logarithmic likelihood function is: where e r is the n Â 1 vector of « r . lnL was maximized when W was as follows: where y is the n Â 1 value vector of aerodynamic coefficients for n sample points, X is the n Â m matrix of regressors and we assume that X T X is invertible. If a sample point was added to the regressor matrix, X n 1 1 ¼ X n x T n 1 1 ! , and the value vector, , then the estimation of W became: We had already known that: Therefore, we had: Define the matrix G n = (X T X) À1 and the real number This was the parameter estimation of the recursive maximum likelihood method.

Orthogonalization procedure of estimation
The matrix is typically badly conditioned when we use explanatory factors as regressors directly. Therefore, an orthogonalization strategy is usually used to solve the problem (Morelli, 1993(Morelli, , 1995. Orthogonalization also made proving the recursive procedure's convergence and unbiasedness in estimations easier. The regressors' column vector was made up of orthogonal functions produced by the Gram-Schmidt process: where the x k was the kth variable. The orthogonal vector was defined as: where p i,j represented the jth orthogonal function with the explainable variables at sample point i. Then, the matrix X T X became: Because P was an orthogonal function system: Figure 3 Control surface layout of the aircraft in the case study Therefore, we had: The parameter estimation was: c W n 1 1 ¼ c W n 1 G n P n 1 1 k À1 y n 1 1 À P T n 1 1 c W n (16)

Proof of unbiasedness and convergence
Now, it is easy to prove thatŷ is an unbiased prediction of y and that the prediction converges to y as the number of sample points increases. The estimate error was: Figure 5 Flow chart of the automatic computing procedure where e and « are remainders consisting of regression residuals and measurement noises because of the inaccuracy of CFD, which also obeys normal distributions N (0, s 2 ) according to the central limit theorem, and thus E b y À y ð Þ¼ 0; i.e. the prediction of y was unbiased.
The variance of prediction error for a series of linear unrelated points was: Note that the first term on the right side of the above equation has a form resembling Ra i d i ; additionally, we know that e and « are independent but obey the same distribution, and thus: If a sample point was added, matrix G became: Because P i,j is orthogonal, (P i,j , P i,j ) > 0, and the matrix G is positive definite, we know that for arbitrary x: That is: This equation means that the variance of the prediction error decreases as the sample size increases. In addition, it is easy to see that kGk 2 , the 2-norm of G approaches zero as the sample size n increases. Thus, the prediction variance error converges to s 2 . Because the regression residual and CFD inaccuracy are unrelated, s 2 equals to the sum of their variances.

Recursive sampling method using Voronoi diagram
If we do not know enough about the aircraft's aerodynamics, the ideal sampling technique is to evenly disperse the sample points to get a more accurate model. Although orthogonal experimental design techniques are extensively used, the majority of them adopt batch processing. Because recursive regression necessitates uniformly distributed samples when each new sample is generated, a new technique must be devised. The center of the largest empty sphere (LES) is considered as the best sampling location in this scenario. The LES problem is defined as follows: find the largest radius hypersphere (within a specific border) that does not overlap with any points in a specific set. Based on LES, the sampling process was completed into two parts. The Latin hypercube technique was used first to obtain the minimum number of samples required for the modeling procedure. The next stage was to solve the LES problem recursively to find the new sample point.
The enumeration approach solves the LES problem, but it has an O(n 4 ) time complexity, which means it will take a long time if a lot of samples are required. The Voronoi diagram has been shown to be a superior tool.
A Voronoi diagram partitions a space into convex hulls (Voronoi cells) and satisfies the following: For a given set of points, each Voronoi cell contains exactly one point belonging to the set, and any point in the cell is closer to this point than any other point from the set.
Points on the hyperplanes between the Voronoi cells are equidistant between the points in each cell.
The Voronoi points, edges and faces are defined as the vertices, edges and faces of the Voronoi cells, respectively. The center of the LES must be a Voronoi point. Solving the LES problem with a Voronoi diagram has a time complexity of only O(nlogn) (Toussaint, 1983).
The additional point would tend to lie on the border because there were no points outside the boundary. To prevent this problem, the maximum convex hull was extracted and mirrored with all the borders. Figure 2 was an example of a two-dimensional sampling, and higher-dimensional sampling is similar. To save processing time, any points in the hyperrectangle with edges double the border edge length were removed after the mirror procedure.

Problem definition
A flying wing working at sea level with a span of b = 1.372 m and an aspect ratio = 1.673 was studied. The aircraft has three types of control surfaces: elevator, aileron and split drag rudder (SDR), as shown in Figure 3. The mean aerodynamic chord of Figure 11 Trimming control and AOA with AOS = 5°G the aircraft was 0.82 m and the airspeed was 70 m/s; thus, the Reynolds number was approximately 3.9 Â 10 6 . The classic wind tunnel axis system was used to deal with the aerodynamic forces and the body axis system was used to deal with the aerodynamic moments.
Assuming that the aircraft is flying at sea level, there are five explanatory variables involved: the AOA, the AOS and three control surface deflections. The ranges of interest of the variables in this work are shown in Table 1. Because we also assume the control surfaces influence each other, i.e. fully coupled, the aerodynamic database has 25,515 records.

Computational fluid dynamics method and validation
The CFD process adopted the finite volume method and the hexahedral grid. The surface mesh of aircraft is shown in Figure 4. The velocity inlet (U 1 = 70 m/s) was set on the front wall of the box calculation domain, whereas the outflow was set on the back wall. All other walls were set to symmetry to prevent boundary layers from forming. All the boundaries were 30b distant from the aircraft. The steady uncompressible reynolds-averaged Navier-Stokes equations and the Spalart-Allmaras turbulence model were used. The result of grid independence test is shown in Table 2. The mesh with approximately 11 million cells produced a satisfactory result.

Determination of the aircraft model structure
The symmetric terms (AOA, d e) are expanded to at least fourth orders to simulate stalling. The AOS terms were expanded to third orders to simulate the slightly nonlinear behavior, whereas the other terms (d a, d r) were expanded to the fifth order to correctly predict some complex phenomena such as nonlinear behavior of SDR or aileron stalling. The explanatory variables were multiplied up to the fifth order in all possible combinations. As mentioned in Section 2.3, the aileron deflection, the SDR deflection and the AOS have the combined symmetry with each other and can be simplified. The final regressors are shown in Table 3.

Automatic computing
The automatic computing progress was built up with sampling procedure and estimating procedure.
In the sampling procedure, a matrix of states was generated. The sampling procedure was carried out within the bounds defined in Table 1. To meet the regression procedure's minimum criterion, a set of 134 points was first constructed using the Latin hyper cube approach. By iteratively solving the Figure 12 Response of the elevator doublet LES issue with a Voronoi diagram, the following 466 points were created. The distances between the Voronoi points and the second-nearest sample points were compared if there were duplicate solutions. This technique would be repeated until no duplicate solutions remained. A total of 600 points were prepared in the end.
The estimating procedure load the state matrix and calculate the result via CFD approach and estimate the model parameters. Yet, a remaining problem was that the values computed via CFD were not always trustworthy. A screening procedure must be carried out to deal with this problem. So, the complete procedure was like the following: 1 Load the first (134 1 30) sample points and calculate corresponding CFD results. The first 134th samples compose the data set, whereas the last 30th samples make up the verification set. 2 Use the data set to estimate model parameters, then determine the model's predict error and standard variation s d . 3 CFD results that were 3s d off from the model prediction were considered outliers and were eliminated from the data set. 4 Repeat steps 2-4 until there are no more outliers. 5 Calculate the verification set's predict error and standard error.
6 The calculating process is finished if the standard variation of the verification set s v is smaller than the absolute error threshold (t ea ) or relative error threshold (t er ). Otherwise, add the first member of the verify set to the data set and load a new sample into the verification set, repeat steps 2-7.
According to the Drag Prediction Workshop through years (Hemsch, 2004;Derlaga and Morrison, 2017), the accuracy of CFD prediction on drag is constantly improving. Drag prediction accuracy is less than 35 counts even at the early results. Our goal is to use CFD in the conceptual design stage. Aerodynamics may change as modifications are introduced in later design stages, such as adding protrusions onto the aircraft surface. Hence, the accuracy requirement of 35 counts was sufficient for this case. The absolute error threshold we selected for all the aerodynamic coefficients was t ea = 0.0035, whereas the relative error was t er = 5%. The complete workflow is shown in Figure 5.

Results
The automatic computing process stops after calculating 470 points, indicating that the model reduced the workload to approximately 1/55 compared to establishing a full database. Figure 13 Time history of the longitudinal coefficients 3.5.1 Verification of aerodynamic prediction A collection of 100 sample points was produced with CFD and compared to the model prediction to verify the model's prediction accuracy. All samples used in this analysis are linearly unrelated to the samples used to calculate the parameters. Figure 6 depicts the predicted error. The horizontal dashed line represented t er (5%) and three times the value of t ea (0.0105). The majority of the absolute errors were less than 0.0105 and those that were not were all less than the 5% relative error level. Outliers were defined as values greater than 3s . Figure 7 depicted the aerodynamics at a = 6°, d a = d e = 0°o n AOS vs d r plane. The model's prediction, displayed as the surface in the figure, fits the calculation result well for most of the coefficients. The C l , which represented the sympathetic derivative, or C ld r , on the other hand, had a low degree of precision. One cause for this was that the coefficient is small and the data for estimation itself had a considerable inaccuracy. This could be demonstrated by comparing the C d and C l values. The increase in C d was not the primary effect of rudder, but it had a larger influence, therefore the prediction was more accurate. The other reason was the estimating process calculates all the parameters of one coefficient at the same time. Although the total error was acceptable, the estimation did not consider the error of a single derivative. As a result, the estimation of sympathetic derivatives was rather imprecise. Nonetheless, the predict error of C l was less than 33%, which is acceptable.

Verification of flight simulation
However, it is widely known that the consistency of the simulation response cannot be guaranteed solely by the closeness of the aerodynamic derivatives. We used the aerodynamic model and the database to run a series of 6-DOF simulations on the same aircraft to further confirm the model's validity. The simulation was created using Simulink and includes four main modules, similar to those found in traditional simulations: control system, aerodynamics, 6-DOF dynamics and atmosphere, as shown in Figure 8. The rigid body dynamics module adopted quaternion, and the numerical integration method was ode45. The controller of each case, including all of the feedback gains, is completely the same. The results of a series of simulations, including steady-state trimming, doublet and chirp responses and coordinated turning, demonstrated that this model and the database would provide identical outcomes. Calculating the elevator deflection angle, throttle position and trim AOA of an airplane at various speeds is one technique to evaluate its trim ability. As shown in Figure 9, a basic proportion integration differentiation (PID) controller was constructed to gather trim information at various speeds.
The trimming result without sideslip is shown in Figure 10. When the speed was low, the model produces more elevator deflection and vice versa. That means the predicted trim capability was lower than the database's capability; yet, throughout the speed range, the elevator deflection was always less than 0.01 radians. The difference in throttle position was small at low speeds, but it grew with speed and stops growing at roughly 3% when the speed hits 70 m/s. The trim AOA difference has been kept small until the speed exceeds 80 m/s, reaching a maximum of 0.12°at 90 m/s. The model performed even better in the trimming case with a sideslip b = 5°, as shown in Figure 11. The model's accuracy was completely demonstrated by the fact that all five variances were very small throughout the speed range.
3.5.2.2 Doublet response Doublets are a type of stimulation widely used in automatic control field. Theoretically, doublets include all the frequency components. Figures 12 and 13 demonstrated the open-loop response and the time history of longitude coefficients. The time history of coefficients showed just a minor change and AOA's response was nearly identical. The differences in speed and height response seemed to be significant, mainly because of the phase difference of Phugoid model at the end of stimulation. Despite this, the differences were only 2 m or 0.1 m/s, which was a little variance.
The open-loop response and the time history of lateral coefficients corresponding to aileron doublet are shown in Figures 14 and 15. There was only a minor difference in the responses. The time histories of C l were essentially identical, although C c and C n differed little. Because aileron primarily creates rolling moment while side force and yawing moment caused by aileron were relatively small, a small absolute error causes a large relative error.
The Theil inequality coefficient (TIC) is a popular model prediction accuracy criterion. TIC is defined as follows: Figure 15 Time history of lateral coefficients where n t is the number of time-domain data points, n 0 is the number of measured output channels, y data is the system output, y was the model output and W is the diagonal matrix of weight factors that is used to balance outputs with various dimensions (Tischler and Remple, 2012). A TIC of less than 0.3 is considered acceptable (Jategaonkar et al., 2004). The measured output was defined as follows: because we assumed the aerodynamics model was fully coupled and the flight simulation used a nonlinear 6-DOF dynamics model Except for weight factor of the airspeed was 0.0573 as suggested, all nonzero elements in the matrix W were one. The elevator doublet reaction had a TIC of 0.184 whereas the aileron doublet response had a TIC of 0.0548, both of which met the requirement. Coordinate turning uses all the control channels of the aircraft, so it is an ideal maneuver to verify if the model precisely predicts the combined control response. The turning was carried out at 100 m height and a consist speed of 70 m/s. The control was performed by a PID controller and consisted of a trimmer, a bank angle retainer and a sideslip eliminator. A switch was used to start the turning, as shown in Figure 18. Figure 19 depicted the controls required to perform a series of coordinate turnings of varying radius. The difference in elevator was just approximately 0.01 rad, whereas the difference in throttle was less than 0.1%. Because the deflection was small, the aileron differential appeared to be a little larger. The difference in rudder was around 17%. One reason was that the sympathetic derivative, or C nd a , was small and was difficult to be solved precisely. As a result, the samples used for parameter estimation were highly noisy, resulting in a substantial inaccuracy in model parameters. As a result, the minor aileron variance resulted in a larger variation in C n .

Figure 16 Response of elevator chirp
Because of the inefficiency of SDR, additional rudder deflection was necessary to eliminate the aileron error, resulting in a relatively large rudder difference. The flight paths were almost coincidence with only a little difference in height, as shown in Figure 20. The overall predict ability of the model is satisfactory.

Conclusions
This work proposed a complete method for constructing an aerodynamic model of an aircraft using the CFD approach. The method aimed on saving CFD effort and used a recursive sampling and parameter estimation method to solve a structure-preconfigured polynomial model. To support the recursive regression approach, a novel sampling process was devised, which included recursive sampling by solving the LES problem via the Voronoi diagram. The feasibility of the recursive estimation method was proven by proving the prediction error of the model decreases with an increasing number of sample points. The model can be simplified because of the lateral symmetry of the aircraft.
To validate the modeling methods, a case study of a flying wing was conducted. Each coefficient in the example contains five explanatory variables and 133 parameters. Compared with creating a database, the modeling method took only 470 samples, or 1/54 of the calculating effort. For all the aerodynamic coefficients, the prediction accuracy reached an absolute error of 0.0105 (or relative error of 5%). The model correctly predicted the dynamic and control of the aircraft in the 6-DOF simulations, which included steady-state flight, response to the doublets and coordinate turning. This model has a higher prediction accuracy for the dominant influence of the control surface than it does for the nondominant effect.