Abstract
Volterra series is one of the powerful system identification methods for representing the nonlinear dynamic system behavior. The methods of step response and impulse response are commonly applied to a discrete aerodynamic Computational Fluid Dynamic (CFD) to identify the first and secondorder Volterra kernels. A critical problem, however, is the difficulty of identifying the secondorder Volterra kernels correctly in CFDbased method. In this paper the secondorder Volterra kernel function is expanded in terms of Chebyshev functions to reduce the size of the problem and the accuracy of the identification is also improved based on a thirdorder reduced model of Volterra series.
1. Introduction
Nowadays, the most powerful and sophisticated tools for predicting nonlinear unsteady aerodynamic characteristics are being developed in the field of computational fluid dynamics (CFD). The application of CFD codes involves, in general, the application of the discretized NavierStokes (NS) equations. For the case where flowinduced nonlinearities at transonic flight conditions are significant in the aerodynamic response, time integration of the nonlinear equations is necessary. Using CFD codes the computational costs become prohibitively expensive, both in terms of CPU costs and turnaround time. This can be on the order of hours or days, depending on the user’s computer, the complexity of the structure model and the condition of flow field. This disadvantage really restricts its applications combining with other disciplines; especially in the design of aircraft control systems and aeroelastic systems.
Attempts to address this problem include the development of transonic indicial responses and impulse responses. Another approach to reduce the computational cost of CFD codes is to linearize the response about a nonlinear steadystate condition, obtain a linear statespace representation of the system of interest. If a linearized aerodynamic model from a CFD model is preferred, the method of the step or pulse input can be applied. In the identification fields, other methods or manners can also be employed. For examples: wavelet approximation and scale transformation, POD method, Volterra variational approach, neural network, incremental harmonic balance methods for airfoil flutter, etc [15].
The discrete impulse input was first suggested by Silva [6] in his dissertation to identify the Volterra kernels of nonlinear aerodynamic systems using the CFL3D code. Three conditions should be satisfied that 1) the kernels, input function, and subsequently, the output function are realvalued functions; 2) the system is causal; 3) the system is time invariant. Continuing on the work of Silva, Raveh generated a reducedorder model (ROM) by using step responses as the ROM’s kernels, for evaluation of generalized aerodynamic force (GAF) matrices. Raveh’s work [7] showed a better accuracy as compared to the Volterra ROM based on impulse response. CFD analysis shows that a step response, in contrast with an impulse method, is less prone to numerical problem and the identification of steptype kernels is less sensitive to such error. Raveh’s work also showed that it was difficult to identify the secondorder Volterra kernels correctly for a particular aerodynamic problem.
There are mainly two reasons listed below:
a) The numerical limitations restrict the impulse/step maximum input amplitude and the minimum time step which can be accurately handled by the CFD numerical scheme. Whichever a pulseresponse or a stepresponse method is used, the data needed for the identification of secondorder kernel function is huge.
b) Based on the secondorder ROM the secondorder kernel function is difficult to be identified correctly.
In this paper, we first set up a thirdorder Volterra ROM to improve the accuracy of identification and then try to identify the secondorder kernel function by using an alternative method, called Chebyshev expansion method. This identification technique does not strongly rely on the entire input and output data produced by CFD simulation. That means a colossal waste of time can be saved.
In Section 2, we briefly introduce the conceptions of Chebyshev polynomials, Volterra series and discrete impulse input. Next, we introduce first and secondorder ROM respectively. In Section 3, we put forward to the thirdorder ROM. Finally, we propose the identification method of Chebyshev polynomials. A CFD numerical scheme is employed to compute the aerodynamic forces of NACA0012 airfoil as the data for identification in transonic flow.
2. Preliminaries
2.1. Chebyshev polynomials
Chebyshev polynomial method is not new; however fewer authors investigated and attempted to apply it into identification in aerodynamics [8]. This method can be classified as parameter identification. Using it needs two steps to fulfill the identification of secondorder kernels of a nonlinear system. The first step is called secondkernels decomposition. The secondorder kernels can be represented in two orthogonal Chebyshev polynomials’ product and then sum up.
Chebyshev polynomials can be expressed in the form below:
All the Chebyshev polynomials form a complete orthogonal set on the interval with respect to the weighting function. It can be show below:
The secondorder kernels can be expressed in term of the orthogonal basis set as:
where $a$ are coefficients need to be identified.
2.2. Volterra series
In 1887, Volterra first studied some nonlinear problems by employing series. The Volterra theory was developed in 1930. The concepts of the identification of linear or nonlinear systems were definitely proposed by Wiener in 1958. Consequently, Volterra series has been applied into various fields, such as system identification, aerodynamic systems, signal and image processing system, adaptive control systems, defect diagnosis, noise suppressing, predictive modeling, etc.
Volterra series can be applied into modeling a timevariant or timeinvariant system in time domain or frequent domain. In this paper, we focus on the timedomain Volterra theory because CFD analyses are typically performed in the time domain. We only study on timeinvariant nonlinear aerodynamic systems [913].
If a timeinvariant nonlinear system can be modeled by Volterra series, it has the continuous time form below:
where $y$ is the response of the nonlinear system to the input $u$, an arbitrary input; ${h}_{0}$ is a steady value which can be computed by executing CFD code; ${h}_{n}$ is the nth kernel of the nonlinear system.
Using the impulse response method to kernels’ identification, the socalled unit impulse input in a continuous or discrete system is defined by:
2.3. Firstorder ROM and secondorder ROM
Based on the Volterra theory the firstorder ROM [1416] in a continuous or discrete system has the form below:
where $h$ is called the linear kernel function.
By employing the impulse input $u={\xi}_{0}\delta $, the linear kernel can be identified in a continuous or discrete form below:
Similarly, the secondorder ROM can be modeled in a continuous form below:
where ${h}_{1}$, ${h}_{2}$ are called first and secondorder kernels respectively.
In a continuous or discrete form, by employing impulse input $u={\xi}_{0}\delta $, one can obtain:
Similarly, employing impulse input $u={\xi}_{1}\delta $, one can obtain:
For simplicity, in the following sections, only discrete expression will be shown.
Let:
Then, the firstkernels in a discrete form is represented by:
In order to identify the secondorder kernels of a discrete nonlinear system, the socalled double impulse function is defined as:
where ${\xi}_{0}$ is a constant, an arbitrary real number; $T$ represents the numbers of time interval between the two pulses in a discrete system.
Feeding Eq. (13) into the counterpart of Eq. (8), namely equation in a discrete form, one can obtain:
$+{\xi}_{0}^{2}\cdot {h}_{2}\left[nT+1,nT+1\right]+2{\xi}_{0}^{2}\cdot {h}_{2}\left[n,nT+1\right],$
where ${y}_{T}\left[n\right]$ is the response to the input ${u}_{T}\left[n\right]$, $\stackrel{}{{y}_{T}}\left[n\right]={y}_{T}\left[n\right]{h}_{0}$.
Under the timeinvariant hypothesis, by using Eqs. (9) and (11), Eq. (14) can be represented as:
Hence, the secondorder kernel's function has the form below:
In a summary, in the firstorder ROM of a nonlinear system of interest, there is only a linear kernel function which needs to be identified. In the secondorder ROM, the firstorder kernels’ function is quite different from the linear kernels' function in the firstorder ROM. One can see it from Eqs. (7) and (12). That implies the kernels identified from the higherorder ROM is more near to the true kernels of nonlinear systems.
3. Thirdorder ROM
As mentioned in the section of introduction, the secondorder kernels sometimes can’t be identified correctly. In order to improve it, the higherorder ROM is proposed here. In a discrete system, it has the form below:
$+\sum _{i=1}^{n}\sum _{j=1}^{n}\sum _{k=1}^{n}{h}_{3}\left[ni+1,nj+1,nk+1\right]\cdot u\left[i\right]\cdot u\left[j\right]\cdot u\left[k\right].$
Similarly to the process of identification of secondorder ROM, firstly let:
Next, feed them into Eq. (17). It follows:
where $\left\right{u}_{i}\left\right={2}^{(i1)}{\xi}_{0}$, $i=\text{1,}\text{}\text{2,}\text{}\text{3}$.
The firstorder kernels and the components of the secondorder kernels obtained by solving Eq. (19) have the forms below:
Feeding ${u}_{T}\left[n\right]={\xi}_{0}\cdot \delta \left[n\right]+{\xi}_{0}\cdot \delta [nT+1]$ into Eq. (17), it follows:
$+{\xi}_{0}^{2}\cdot {h}_{2}\left[nT+1,nT+1\right]+2\cdot {\xi}_{0}^{2}\cdot {h}_{2}\left[n,nT+1\right],$
where those thirdorder kernels’ terms are omitted since they are very smaller.
One can obtain the secondorder kernels in the form as following:
${\xi}_{0}^{2}\cdot {h}_{2}[nT+1,nT+1])/(2\cdot {\xi}_{0}^{2}),$
where $\mathrm{}{h}_{1}\left[n\right]$, ${h}_{2}[n,n]$ are defined by Eqs. (20) and (21) respectively; ${h}_{2}[nT+1,nT+1]$ can be evacuated by a time shift of ${h}_{2}[n,n]$ according to the timeinvariant hypothesis.
4. Example – nonlinear circuit
In the study of aerodynamics, the common methods include impulseresponse method and stepresponse method. All the data for identification come from executing CFD numerical scheme. Before executing CFD code, the amplitude of an input must be chosen carefully whatever using an impulseor stepresponse method since it may arise in different nonlinearities.
In order to identify the secondorder kernels correctly, the thirdorder ROM is introduced here rather than the loworder ROMs studied by many other authors. One of them, Raveh [2], ever claimed in his paper that the problem to correctly identify the secondorder kernels of nonlinear systems was quite difficult to a particular application.
For a simple nonlinear circuit system [1, 18], the governing equation is given below:
where $v$ is an input of voltage; $i$ is the output of the current around the circuit.
This nonlinear system can be remodeled by the secondorder ROM as seen in Eq. (8). Using Eqs. (12) and (16), the first and secondorder kernels can be evacuated in an impulseresponse method. They are presented in Figs. 1 and 2. As can be seen, the secondorder kernel function is quite small in magnitude as compared to the firstorder kernel function and goes to zero very quickly as well as the firstorder kernel. Figs. 1 and 2 indicate that nonlinear effects are quite small.
Fig. 1Firstkernel function of Riccati equation
Fig. 2Secondkernel curve of Riccati equation
Using Chebyshev polynomials method for identification of the second kernel mentioned as in Eq. (3), the first thing is that we have to determine how many basis functions are properly chosen? Here, by testing after several times we find out that the number of basis functions is between 10 and 12. Hence, there are fiftyfive coefficients need to be identified if we choose the number 10. Based on the orthogonal properties of Chebyshev polynomials, see in Eq. (2), all the coefficients can be evacuated. Substituting them into Eq. (3), one can obtain the secondorder kernel function which is presented in Fig. 3.
Let $s=$–1, the secondorder kernels identified respectively by using an impulseresponse method and Chebyshev polynomials method are presented in Fig. 4.
This example illustrates Chebyshev polynomials method can be used for deification of the secondorder kernel function fully accurate in contrast to that identified by using an impulseresponse method.
Fig. 3Secondkernel by Chebyshev polynomial
Fig. 4Comparison of secondkernels respectively by Chebyshev method and impulse method
5. Example identification in transonic flows
During the last several decades there has been considerable interest in the accurate prediction of steady and unsteady transonic flows. Since transonic flow is inherently nonlinear and involves the formation of shock waves, accurate simulation of transonic flow requires the use of a nonlinear field equation, and for unsteady phenomena the solution needs to be time accurate. This paper employs a CFD (Computational Fluid Dynamics) code which is based on the Euler equation.
The twodimensional Euler equations in Cartesian coordinates ($x$, $y$) can be written in nondimensional conservative from as following [17]:
where $W=\left(\begin{array}{c}\rho \\ \rho u\\ \rho v\\ \rho E\end{array}\right)$ is the vector; $F=\left(\begin{array}{c}\rho U\\ \rho uU+p\\ \rho vU\\ U\left(\rho E+p\right)+{x}_{t}p\end{array}\right)$ and $G=\left(\begin{array}{c}\rho V\\ \rho uV\\ \rho vV+p\\ V\left(\rho E+p\right)+{y}_{t}p\end{array}\right)$ are the Cartesian components of the convective fluxes.
In the above $\rho $ is the fluid density, $u$ and $v$ are the $x$ and $y$ components of the fluid velocity, $E$ the total energy and $p$ is the pressure which is obtained from:
where $\gamma $ is the ratio of specific heats of the fluid.
The contravariant velocities $U$ and $V$ are defined by:
where ${x}_{t}$, ${y}_{t}$ are the grid speeds in the $x$ and $y$ directions respectively.
The Eq. (25) are discretized using a cellcentered finite volume method that transforms the partial differential equations into a set of ordinary differential equations that can be written as:
where ${R}_{i,j}\left(W\right)$ is the flux residual for the cell ($i$, $j$) that contains all the terms arising from the spatial discretization.
For the solution of the unsteady Euler equations, an implicit dualtime method can be employed. The initial conditions are chosen as the steady flow before the object moves. The boundary condition at the body is subject to the kinematic condition requiring the body to be impenetrable to the flow. Namely, the flow remains tangent to the body surface. The farfield boundary conditions are treated by the characteristic boundary method based on the Riemann invariants.
Fig. 5 is an Otype grid with 181×48 grid points for NACA0012. The steady inviscid comprssible flows about this airfoil is computed at Mach number ${M}_{\infty}=\text{0.8}$ and angle of attack $\alpha =\text{0}$ deg. The model has a NACA0012 airfoil secontion and a rectangular planform with a span of 0.8128 m and a chord of 0.2032 m. The rigidbody plunge mode consists only of vertical translation and the rigidbody pitch mode consists onlyof rotation about the midchord [19].
In aerodynamics, using Chebyshev polynomial method for the identification of kernels of the aerodynamic forces is seldom seen. Much more of its application can be found in theoretical research. The popular identification methods to a nonlinear system, say pulse and steptype methods, highly rely on the entire data produced by executing CFD code again and again. As mentioned in the section of introduction, the cost of timeconsumption is prohibitively expensive.
If the secondorder kernels can not be identified correctly, the result of the aerodynamic forces of predicting can not be accurate and even badly wrong. Therefore, the thirdorder ROM is first suggested to be used here. The process of identification is presented in the section of thirdorder ROM above. The mainly advantages here are a) the first and secondorder kernels are steadily to be evacuated by Eqs. (20) and (23); b) Since the thirdorder kernels are not necessary to be identified, this model is appropriate to be applied; c) The highorder ROM helps to improve the identification accuracy of secondorder kernels of nonlinear systems of interest.
Fig. 5Naca0012 airfoil and aerodynamic grid
In order to investigate the kernels' accuracy from the two models of ROMS of the secondorder and the thirdorder, the first and secondorder kernels have been computed out respectively from Eqs. (20) and (23) in the thirdorder ROM and from Eqs. (12) and (16) in the secondorder ROM, here ${\xi}_{0}=\text{0.05}$ is chosen and impulseresponse method is employed.
The secondorder kernels can also be computed out from Eq. (3) by Chebyshev Polynomials method based on the coefficients identification which needs use the least squares estimation to a number of data produced by executing CFD code. For simplicity, we call the kernels identification of 2ROM Silva Method and call the Chebyshev Polynomials identification of 3ROM our New Method.
Fig. 6 shows the difference of the two firstkernels identified by using New Method and Silva Method respectively. Fig. 7 shows the prediction of aerodynamic forces by using CFD directly computational method and using convolution techniques to the same sinusoid input with the firstkernels of 2ROM and 3ROM respectively. It illustrates that the higher of the ROM is used the more accuracy the kernels are.
Fig. 6First kernels obtained in Silva Method and New Method
Fig. 7Lift predicted by using CFD method and convolution techniques respectively with firstkernels of 2ROM and 3ROM to the same sinusoid input
Note that the accuracy of the prediction of aerodynamic forces using CFD directly computational method is dependent on the choice of the time step. Fig. 8 shows the difference of the accuracy when executing CFD code with the step time 0.001 and 0.0001 respectively. The input is given by $u\left(t\right)=0.1\mathrm{s}\mathrm{i}\mathrm{n}\left(108.36t\right)$. However, it does not a matter to consider the selection of time step size in the view of identification which is merely concerned about data of inputs and outputs.
Fig. 8CFD Results based on time step 0.001, 0.0001
Using CFD method, impulseresponse method (Silva method based on 2ROM) and Chebyshev polynomials method (our new method based on 3ROM), the response curves is given respectively in Fig. 9.
Fig. 9Response curves computed by CFD, Chebyshev method and Silva method respectively
6. Conclusions
In order to improve the accuracy of identification, the thirdorder ROM is first proposed to model the nonlinear and unsteady aerodynamic forces. The second example shows the kernels of thirdorder ROM are much more accurate than those of a secondorder ROM.
In order to overcome the prohibitively expensive timeconsumptions by executing CFD code to get the fully data for identification in an impulse or stepresponse method, Chebyshev polynomials method is suggested and investigated in details. This method of identification is classified as a parameter method. Since only parts of the data are used for identification, the number of executing CFD code can be reduced at least one order of magnitude. That is extremely timesaving.
References

Jeffrey P. T., Dowell E. H., et al. Threedimensional transonic aeroelasticity using proper orthogonal decomposition based on reduced order models. American Institute of Aeronautics and Astronautics, Paper No. 20011526, 2001.

Jeffrey P. T., Dowell E. H., et al. A harmonic balance approach for modeling threedimensional nonlinear unsteady aerodynamics and aeroelasticity. International Mechanical Engineering Congress and Exposition, Paper No. 200232532, 2002, p. 13231334.

Jawad K., Wu Z. G., Yang C. Volterra kernel identification of MIMO aeroelastic system through multiresolution and multiwavelets. Computational Mechanics, Vol. 49, Issue 4, 2012, p. 431458.

Bommanahal M., et al. Nonlinear unsteady aerodynamic modeling by Volterra variational approach. American Institute of Aeronautics and Astronautics, Paper No. 20124654, 2012.

Antonios K. A., Achilleas D. Z. Wavelet Neural Networks: A Practical Guide. Neural Networks, Vol. 42, 2013, p. 127.

Silva W. A. Discretetime linear and nonlinear aerodynamic impulse responses for efficient CFD analyses. Ph. D. Thesis., College of William & Mary Williamsburg, 1997.

Raveh D. E., Mavris D. N. Reducedorder models based on CFD impulse and step responses. American Institute of Aeronautics and Astronautics, Paper No. 20011527, 2001.

Wang Y. H., Han J. L. Approach to identification of a secondorder Volterra kernel of nonlinear systems by Chebyshev polynomials method. Research Journal of Applied Sciences Engineering and Technology, Vol. 5, Issue 20, 2013, p. 49504955.

Lind R., Prazenica R. J., et al. Identifying parameterdependent Volterra kernels to predict aeroelastic instabilities. American Institute of Aeronautics and Astronautics Journal, Vol. 43, Issue 3, 2005, p. 24962502.

Marzocca P., Librescu L., Silva W. A. Volterra series approach for nonlinear aeroelastic response of 2D lifting surfaces. American Institute of Aeronautics and Astronautics, Paper No. 20011459, 2001.

Mirri D., Iuculano G., et al. Nonlinear dynamic system modeling based on modified Volterra series approaches. Measurement, Vol. 23, 2003, p. 921.

Lind R. D., Prazenica R. J., et al. Estimating nonlinearity using Volterra kernels in feedback with linear models. American Institute of Aeronautics and Astronautics, Paper No. 20031406, 2003.

Silva W. A. Identification of nonlinear aeroelastic systems based on the Volterra theory: progress and opportunities. Nonlinear Dynamics, Vol. 39, 2005, p. 2562.

Raveh D. E., Reduced order models for nonlinear unsteady aerodynamics. American Institute of Aeronautics and Astronautics, Paper No. 20004785, 2000.

Muteanu S., et al. Volterra kernel reducedorder model approach for nonlinear aeroelastic analyses. American Institute of Aeronautics and Astronautics, Paper No. 20051854, 2005.

Silva W. A. Recent enhancements to the development of CFDbased aeroelastic reducedorder models. American Institute of Aeronautics and Astronautics, Paper No. 20072051, 2007.

Dubuc L., Cantariti F., Woodgate M., Gribben B., Badcock K. J., Richards B. E. Solution of the unsteady Euler equations using an implicit dualtime method. American Institute of Aeronautics and Astronautics Journal, Vol. 36, Issue 8, 1998, p. 14171424.

Wang Y. H., Han J. L., Zhang T. T. Computation of Volterra kernels’ identification to Riccati nonlinear equation. International Conference on Computer Science and Information Technology, Vol. 6, 2010, p. 68.

Rivera J. A., Dansberry B. E., Bennett R. M., et al. NACA 0012 benchmark model experimental flutter results with unsteady pressure distributions. American Institute of Aeronautics and Astronautics, Paper No. 19922396Cp, 1992.
About this article
This paper is supported by Specialized Research Fund for the Doctoral Program of Higher Education of China (No. 20113218110002) and a project funded by the Priority Academic Program Development of Jiangsu Higher Education Institutions of China.