An effort saving method to establish global aerodynamic model using CFD

Abstract

Purpose

The typical approach of modeling the aerodynamics of an aircraft is to develop a complete database through testing or computational fluid 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 preconfigured 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 flying 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 flight 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-fidelity models during primary design, allowing for more precise flight dynamics analysis.

Citation

Xie, J., Huang, J., Song, L., Fu, J. and Lu, X. (2022), "An effort saving method to establish global aerodynamic model using CFD", Aircraft Engineering and Aerospace Technology, Vol. 94 No. 11, pp. 1-19. https://doi.org/10.1108/AEAT-10-2021-0299

1. 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 data-driven 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 data-driven models that generalize beyond the training data (Brunton et al., 2020). 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 (Brunton et al., 2020).

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.

2. Methods involved in modeling

2.1 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.

2.2 Aerodynamic coefficients involved

The aerodynamic force coefficients and coordinate system are defined identically to their classic definitions in wind tunnel tests:

(1) [CDCYCL]=1qS[−cos⁡ α cos⁡βsin⁡βsin⁡ α cos⁡β−cos⁡ α sin⁡βcos⁡β−sin⁡ α sin⁡βsin⁡α0−cos⁡ α][XYZ]

The aerodynamic moment coefficients and coordinate system are defined identically to their classic definitions in flight dynamics:

(2) [ClCmCn]=1qS[1/b1/c1/b][LMN]

where the forces X, Y and Z and the moments L, M, and N are in the body frame.

2.3 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 (Brunton et al., 2020), 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 (Brunton et al., 2020).

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

(3) {Cs(x)=Cs(−x)Ca(x)=−Ca(−x)

where C_s = C_D, C_L, C_M, C_a = C_Y, C_l, C_N and x = β, δ_a, δ_r. 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 higher-order 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.

2.4 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:

(4) y=xTW+εr

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, σ²) according to the central limit theorem. The logarithmic likelihood function is:

(5) ln ⁡L=−n2ln⁡(2π)−n2ln⁡σ2−εr′εr2σ2

where ε_r is the n × 1 vector of ε_r. lnL was maximized when W was as follows:

(6) W^=(XTX)−1XTy

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^TX is invertible. If a sample point was added to the regressor matrix, Xn+1=[Xnxn+1T], and the value vector, yn+1=[ynyn+1], then the estimation of W became:

(7) W^n+1=(XnTXn+xn+1xn+1T)−1(XnTyn+xn+1Tyn+1)

We had already known that:

(8) (A+uvT)−1=A−1−A−1u(1+vTA−1u)−1vTA−1

Therefore, we had:

(9) W^n+1={(XnTXn)−1−(XnTXn)−1xn+1[1+xn+1T(XnTXn)−1xn+1]−1xn+1T(XnTXn)−1}(XnTyn+xn+1Tyn+1)

Define the matrix G_n = (X^TX)⁻¹ and the real number k=1+xn+1TGnxn+1:

(10) W^n+1=W^n+Gnxn+1k−1(yn+1−xn+1TW^n)

This was the parameter estimation of the recursive maximum likelihood method.

2.5 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, 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:

(11) p0=1, pk+1=xk+1−∑i=1k(xk+1,pi)(pi,pi)pi

where the x_k was the kth variable. The orthogonal vector was defined as:

(12) Pi=[pi,0,pi,1,⋯,pi,m−1]T

where p_i,j represented the jth orthogonal function with the explainable variables at sample point i.

Then, the matrix X^TX became:

(13) XTX=[∑i=1n(Pi,0,Pi,0)∑i=1n(Pi,0,Pi,1)⋯∑i=1n(Pi,0,Pi,m−1)∑i=1n(Pi,1,Pi,0)∑i=1n(Pi,1,Pi,1)⋯∑i=1n(Pi,1,Pi,m−1)⋮⋮⋱⋮∑i=1n(Pi,m−1,Pi,0)∑i=1n(Pi,m−1,Pi,1)⋯∑i=1n(Pi,m−1,Pi,m−1)]

Because P was an orthogonal function system:

(14) XTX=diag{∑i=1n(Pi,0,Pi,0),∑i=1n(Pi,1,Pi,1),⋯,∑i=1n(Pi,m−1,Pi,m−1)}

Therefore, we had:

(15) Gn=diag{1∑i=1n(Pi,0,Pi,0),1∑i=1n(Pi,1,Pi,1),⋯,1∑i=1n(Pi,m−1,Pi,m−1)}

The parameter estimation was:

(16) W^n+1=W^n+GnPn+1k−1(yn+1−Pn+1TW^n)

2.6 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:

(17) y^−y|=xT(W^−W)+ε|=xT[(XTX)−1XTy−W]+ε|=xT[(XTX)−1XT(XW+ε)−W]+ε|=xTGnXTε+ε

where ε and ε are remainders consisting of regression residuals and measurement noises because of the inaccuracy of CFD, which also obeys normal distributions N (0, σ²) according to the central limit theorem, and thus E(y^−y)=0; i.e. the prediction of y was unbiased.

The variance of prediction error for a series of linear unrelated points was:

(18) D(y^−y)=D(xTGnXTε+ε)

Note that the first term on the right side of the above equation has a form resembling Σa_iδ_i; additionally, we know that ε and ε are independent but obey the same distribution, and thus:

(19) D(y^−y)=xTGnXTσ2+σ2

If a sample point was added, matrix G became:

(20) Gn+1=diag{1∑i=1n(Pi,0,Pi,0)+(Pn+1,0,Pn+1,0),⋯,1∑i=1n(Pi,m−1,Pi,m−1)+(Pn+1,m−1,Pn+1,m−1)}

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:

(21) xTGn+1x<xTGnx

That is:

(22) D(y^n+1−y)<D(y^n−y)

This equation means that the variance of the prediction error decreases as the sample size increases. In addition, it is easy to see that ‖G‖², the 2-norm of G approaches zero as the sample size n increases. Thus, the prediction variance error converges to σ². Because the regression residual and CFD inaccuracy are unrelated, σ² equals to the sum of their variances.

2.7 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⁴) 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.

3. Case study: a flying wing

3.1 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 the aircraft was 0.82 m and the airspeed was 70 m/s; thus, the Reynolds number was approximately 3.9 × 10⁶. 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.

3.2 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_∞ = 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.

3.3 Determination of the aircraft model structure

The symmetric terms (AOA, δ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 (δa, δ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.

3.4 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 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 + 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 σ_d.

3. CFD results that were 3σ_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 σ_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.

3.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.

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 3σ.

Figure 7 depicted the aerodynamics at α = 6°, δa = δe = 0° on AOS vs δ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_lδ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.

3.5.2 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.

3.5.2.1 Steady-state flight.

The importance of steady-state flight cannot be overstated. It stands for trim capability in the realm of aircraft design, which is important to aircraft excellence. It serves as an initial condition for the simulation of all other situations in the realm of flight simulation. 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 β = 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:

(23) TIC=1nt⋅n0∑i=1nt[(ydata−y)TW(ydata−y)1nt⋅n0∑i=1nt(yTWy)+1nt⋅n0∑i=1nt(ydataTWydata)

where n_t is the number of time-domain data points, n₀ 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

(24) y=[αvqθβprϕ]T

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.

3.5.2.3 Chirp response.

Chirp signals are another sort of stimulation that are commonly used. The chirp stimulation used here lasted 10 s and ranged from 0.1 to 10 Hz. The open-loop responses were presented in Figures 16 and 17, with TICs of 0.0974 and 0.0776 for each case, respectively.

3.5.2.4 Coordinate turning.

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_nδ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. 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.

4. 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.

Figures

Figure 1

Symmetry of the combination of the AOS and the aileron deflection

Figure 2

2-D example of recursive sampling using a Voronoi diagram

Figure 3

Control surface layout of the aircraft in the case study

Figure 4

Surface mesh of the aircraft

Figure 5

Flow chart of the automatic computing procedure

Figure 6

Prediction error of the model

Figure 7

Predicted values and prediction errors

Figure 8

6-DOF flight simulation framework

Figure 9

Control system for trimming

Figure 10

Trimming control and AOA with AOS = 0°

Figure 11

Trimming control and AOA with AOS = 5°

Figure 12

Response of the elevator doublet

Figure 13

Time history of the longitudinal coefficients

Figure 14

Response of aileron doublet

Figure 15

Time history of lateral coefficients

Figure 16

Response of elevator chirp

Figure 17

Response of aileron chirp

Figure 18

Control system of coordinate turning

Figure 19

Control required for making coordinate turnings

Figure 20

Flight path comparison of turning with a radius of 1,200 m

Table 1

Computed states

AOA	AOS	Aileron	Elevator	Rudder
−1	−8	−30	−20	−40
1	−4	−25	−10	−30
2	0	−20	0	−20
3	4	−10	10	−10
4	8	0	20	0
5	–	10	25	10
7	–	20	30	20
10	–	25	–	30
15	–	30	–	40

Table 2

Verification of grid independence

Mesh density	Number of elements	Cd (at α = 5°, β = 0°)	CL (at α = 5°, β = 0°)	CM (at α = 5°, β = 0°)
Coarse	6 × 106	0.0139	0.2365	−0.2621
Medium	1.1 × 107	0.0142	0.2377	−0.2615
Fine	2.3 × 107	0.0142	0.2377	−0.2615

Table 3

Regressors of the model

Orders	Regressors
	Symmetric					Antisymmetric
0	1	–	–	–	–	–	–	–	–	–
1	α	δe	α2	β2	δa2	β	δa	δr	–	–
2	δe2	δr2	αδe	βδa	βδr	αβ	αδa	αδr	βδe	δaδe
	δaδr	–	–	–	–	δeδr	–	–	–	–
3	α3	δe3	α2δe	β2δe	δa2δe	β3	δa3	δr3	α2β	α2δa
	αβ2	αδa2	αδe2	αδr2	δeδr2	α2δr	β2δa	β2δr	δa2δr	δe2δr
	αβδa	αβδr	αδaδr	βδaδe	βδaδr	βδa2	βδe2	βδr2	δaδe2	δaδr2
	βδeδr	δaδeδr	–	–	–	αβδe	αδaδe	αδeδr	βδaδr	–
4	α4	δa4	δe4	δr4	αδe3	αβ3	αδa3	αδr3	βδe3	δaδe3
	βδa3	βδr3	δaδr3	α3δe	β3δa	δeδr3	α3β	α3δa	α3δr	β3δe
	β3δr	δa3δr	α2β2	α2δa2	α2δe2	δa3δe	δe3δr	αβδa2	αβδe2	αβδr2
	α2δr2	β2δa2	β2δe2	β2δr2	δa2δe2	αδaδe2	αδaδr2	βδaδr2	βδeδr2	δaδeδr2
	δa2δr2	δe2δr2	αδeδr2	βδaδe2	βδaδr2	αβ2δa	αβ2δr	αδa2δr	αδe2δr	βδa2δe
	αβ2δe	αδa2δe	βδa2δr	βδe2δr	δaδe2δr	βδa2δr	α2βδe	α2δaδe	α2δeδr	β2δaδe
	α2βδa	α2βδr	α2δaδr	β2δaδr	αβδaδe	β2δaδr	β2δeδr	δa2δeδr	αβδaδr	βδaδeδr
	αβδaδr	αβδeδr	αδaδeδr	βδaδeδr	–	δa5	δr5	βδa4	βδe4	βδr4
5	αδa4	αδe4	αδr4	δeδr4	α4δe	δaδe4	δaδr4	α4β	α4δa	α4δr
	δa4δe	α2δe3	β2δe3	δa2δe3	α3β2	δa4δr	δe4δr	α2β3	α2δa3	α2δr3
	α3δa2	α3δe2	α3δr2	δe3δr2	αβδa3	β2δa3	β2δr3	δa2δr3	δe2δr3	β3δa2
	αβδr3	αδaδr3	βδaδe3	βδaδr3	βδeδr3	β3δe2	β3δr2	δa3δe2	δa3δr2	αβδe3
	δaδeδr3	αβ3δa	αβ3δr	αδa3δr	βδa3δe	αδaδe3	αδeδr3	βδaδr3	αβ3δe	αδa3δe
	βδa3δr	βδe3δr	δaδe3δr	α3βδa	α3βδr	αδe3δr	βδa3δr	α3βδe	α3δaδe	α3δeδr
	α3δaδr	β3δaδe	β3δaδr	β3δeδr	δa3δeδr	β3δaδr	α2β2δa	α2β2δr	α2δa2δr	α2δe2δr
	α2β2δe	α2δa2δe	β2δa2δe2	β2δa2δr	α2δeδr2	β2δa2δr	β2δe2δr	δa2δe2δr	α2βδa2	α2βδe2
	β2δaδr2	β2δeδr2	δa2δeδr2	αβ2δa2	αβ2δe2	α2βδr2	α2δaδe2	α2δaδr2	β2δaδe2	β2δaδr2
	αβ2δr2	αδa2δe2	αδa2δr2	αδe2δr2	βδa2δr2	βδa2δe2	βδa2δr2	βδe2δr2	δaδe2δr2	α2βδaδr
	α2βδaδe	α2βδaδr	α2βδeδr	α2δaδeδr	β2δaδeδr	β2δaδeδr	αβ2δaδe	αβ2δaδr	αβ2δeδr	αδa2δeδr
	αβ2δaδr	βδa2δeδr	αβδa2δr	αβδe2δr	αδaδe2δr	βδa2δeδr	αβδa2δe	αβδa2δr	βδaδe2δr	αβδaδr2
	βδaδe2δr	αβδaδe2	αβδaδr2	βδaδeδr2	αβδaδeδr	αβδeδr2	αδaδeδr2	βδaδeδr2	αβδaδeδr