A NON-PARAMETRIC TECHNIQUE FOR AERODYNAMIC STATE IDENTIFICATION FOR FLUTTER PREDICTION

The objective of the present work is to study forms of identifying the aerodynamic states, for state space aeroelastic stability analyses in the frequency domain, using highfidelity CFD codes in an efficient fashion. In particular, there is interest in the use of nonparametric identification techniques in order to obtain the necessary aerodynamic transfer functions from a single unsteady CFD computation. The CFD tool used is based on the 2-D Euler equations, which are discretized using a finite volume approach for unstructured grids. A centered scheme, with added artificial dissipation, is used for spatial discretization and explicit Runge-Kutta methods are employed for time marching. The work implements a single combined input method in order to simultaneously excite the aerodynamic responses in all the system natural modes and, at the same time, allow the output to be split into the contribution of each individual mode to the aerodynamic transfer function. The proposed approach is applied to a NACA 0012 airfoil-based typical section model in transonic flight. Results are compared to a previously implemented identification capability in which the aerodynamic transfer functions are computed with an unsteady simulation for each modal movement.


INTRODUCTION
It is well-known that, in recent years, computational fluid dynamics (CFD) techniques have played an increasingly important role in many applications which require understanding the aerodynamics of a given configuration.This statement also holds for aeroelastic analyses.However, the use of high-fidelity CFD solvers, based on the Euler or the Navier-Stokes equations, usually leads to fairly long computational times, especially for the unsteady calculations that would be required in aeroelasticity.Hence, the actual issue addressed in the present work is the matter of reducing the computational cost of such CFD calculations such that the overall cost of an aeroelastic study would be acceptable in industrial conditions.The present effort builds upon previous work by Marques and Azevedo [1,2], which used unsteady CFD solutions to impulsive and indicial perturbations in each structural mode of the configuration in order to build the required aerodynamic transfer functions and, then, identify aerodynamic states that would allow aeroelastic stability analyses in the frequency domain.The drawback with such approach is that the construction of an aeroelastic stability root locus, for a given flight condition, costs one steady CFD run, in order to obtain the initial steady state solution at the flight condition of interest, plus as many unsteady CFD runs as the number of mode shapes considered in the analysis.For a real engineering problem, in which typically some 30 or more modes are considered, this cost is not yet acceptable.Actually, it would probably be more cost effective, for such a case, to construct the stability root locus for that flight condition by straightforward simultaneous integration in time of the aerodynamic and structural dynamic equations, for a sequence of increasing dynamic pressures.
Therefore, the objective of the present work is to study forms of identifying the aerodynamic states, for state space aeroelastic stability analyses in the frequency domain, using high-fidelity CFD codes in an efficient fashion.In particular, there is interest in the use of non-parametric identification techniques in order to obtain the necessary aerodynamic transfer functions from a single unsteady CFD computation.Hence, the construction of an aeroelastic stability root locus could be accomplished with only two CFD runs, one steady calculation, to set up the mean steady flow condition, and one unsteady simulation, which would yield the aerodynamic response to an appropriate input.The approach adopted here follows the work of Raveh [3,4], Kim et al. [5] and Silva [6].The paper discusses the details of the methodology adopted in as much detail as the available space allows, and it applies the proposed approach to a NACA 0012 airfoil-based typical section model in transonic flight.
The CFD tool used here is based on the 2-D Euler equations.These equations are discretized using a finite volume approach for unstructured grids.A centered scheme, with added artificial dissipation, is used for spatial discretization and explicit Runge-Kutta methods are employed for time marching.The paper concentrates initially in obtaining impulsive and indicial aerodynamic responses for the typical section model.In this initial stage, discrete step and "step like" motions motions are applied individually to each modal coordinate of the aeroelastic system, and the corresponding unsteady aerodynamic responses are calculated individually.Afterwards, the work implements a single combined input method in order to simultaneously excite the aerodynamic responses in all the natural modes of the configuration and, at the same time, allow the output to be split into the contribution of each individual mode to the aerodynamic transfer function.The technique for identifying these transfer functions is based on the calculation of the power spectral densities of the inputs and outputs of the system.As will be further discussed along the text, if the system inputs are uncorrelated, the transfer functions can be estimated by a division, in the frequency domain, of the cross power spectral density of the corresponding input and output, and the power spectral density of the input.Results are compared to the previous capability [1,2] in which the aerodynamic transfer functions are computed with an unsteady simulation for each modal movement.

FLOW SOLVER
The CFD tool applied in this work is based on the 2-D Euler equations, which represent two-dimensional, compressible, rotational, inviscid and nonlinear flows.Therefore, it is completely capable of capturing the shock waves present in transonic flows.Due to the use of unstructured meshes and the adoption of the finite volume approach, these equations are writ-ten in Cartesian form.Furthermore, as usual in CFD applications, flux vectors are employed and the equations are nondimensionalized.Hence, they can be written as In Eq. 1, Ω represents the volume of the control volume or, more precisely, its area in the two-dimensional case.S is its surface, or its side edges in 2-D.Q is the vector of conserved properties of the flow, given by E and F are the inviscid flux vectors in the x and y directions, respectively, defined as The nomenclature adopted here is the usual one in CFD: ρ is the density, u and v are the Cartesian velocity components and e is the total energy per unity of volume.The pressure, p, is given by the perfect gas equation, written as Once again, as usual, γ represents the ratio of specific heats.The contravariant velocity components, U and V , are determined by where x t and y t are the Cartesian components of the mesh velocity in the unsteady case.The aerodynamic equations are discretized using a cell centered, finite volume scheme.The code was developed to be used with unstructured meshes composed of triangles.Spatial discretization is equivalent to a centered scheme [7,8].As usual with centered schemes, artificial dissipation terms must be added to the discretized equations in order to maintain numerical stability, especially near shock waves.In the present case, the artificial dissipation operator is constructed as a blend of second and fourth difference terms [9] with coefficients that are proportional to local pressure gradients.This dissipation operator uses a scaling of both undivided Laplacian and biharmonic operators which is similar to the one presented by Mavriplis [8].Such scaling is constructed such that both operators are multiplied by a term proportional to the integral of the maximum eigenvalue of the Euler equations in the direction normal to the cell faces (edges).
After the complete spatial discretization and the inclusion of the artificial dissipation terms, the Euler equations can be written for the i-th volume as where C and D represent the convective and the artificial dissipation operators, respectively.Equation 6 is advanced in time using a second-order accurate, 5-stage, explicit, hybrid Runge-Kutta time stepping scheme [7,8].The convective operator is evaluated at every stage of the integration process, but the artificial dissipation operator is only evaluated at the two initial stages [9].Moreover, convergence acceleration techniques, such as local time stepping and implicit residual smoothing [9,10], are employed in steady calculations in order to guarantee an acceptable computational efficiency in steady-state mode.Further details on the approach here used for the numerical discretization of the Euler equations and on the construction of the present code can be seen in Ref. [1].

STATE SPACE FORMULATION
A state space representation of the typical section system can be written as where {x ( t)} is the system state vector, which can be written according to Here, {η ( t)} represents the two degrees of freedom of the typical section system, i.e., the plunge and pitch modes, respectively.The [ M ] and matrices are defined in Ref. [2].The forcing term is where { Qa ( t) } is the appropriately normalized vector of the generalized modal aerodynamic forces [2].
If one is interested in representing the aerodynamic forces in the frequency domain, the system can be conveniently studied in the Laplace domain.Applying the Laplace transform to Eq. 7, one obtains where s = s ω r (13) is the dimensionless Laplace transform complex variable and ω r is a reference circular frequency.The idea behind this procedure is to evaluate the aerodynamic influence coefficient matrix over a reduced frequency range of interest and, by making use of the analytical continuation principle [11], to extend such result to the entire complex plane.The data, that would come from the CFD solver, however, even after an appropriate identification of the aerodynamic transfer functions in the frequency domain, would consist of sets of numerical values of the aerodynamic coefficients as a function of (reduced) frequency.This format of the aerodynamic data is not convenient for the solution of Eq. 12.One approach for dealing with the problem consists in approximating these data using interpolating polynomials, which then would lead to the creation of (new) aerodynamic states.The literature reports on a number of different interpolating polynomials that could be used in the present case, and the interested reader is referred to Ref. [2] for a detailed discussion of this issue.In any event, regardless of which specific polynomial is used for the representation of the aerodynamic response in the Laplace domain, the resulting state space representation of the typical section aeroelastic system can be written as Here, {χ ( t)} is the new state vector that results from the addition of the aerodynamic state variables, and [D] is the system dynamic matrix.This matrix is defined in terms of the system mass and stiffness matrices, and of other matrices which appear from the interpolating polynomial representation of the aerodynamic forces.Further details can, again, be seen in Ref. [2].Finally, the aeroelastic stability analysis of the system can, then, be reduced to the classical eigenvalue problem, for each value of some characteristic parameter, as where [I] is the identity matrix, N = 2n a + 4, and n a is the number of added aerodynamic state variables.Typically, the characteristic parameter is a reduced speed (U * ) or a reduced dynamic pressure.

IMPULSIVE RESPONSES AND DISCRETE SOLUTIONS
Linear time-invariant continuous-time systems present the very interesting property of possessing elementary solutions which are representative of the entire frequency content of the system response.The most elementary solution is the impulsive one [12].In a frequency domain point of view, the Fourier transform of an impulse function assumes a unit value over the entire domain.In other words, an impulse input has the capability of uniformly exciting a system in the entire frequency domain.Hence, the system response to an arbitrary input can be obtained by the convolution integral of the impulsive response with the given input.Similarly, the step, or indicial, response is another elementary solution.Just as reasoned before, the convolution integral of the indicial response, but at this time with the input time derivative, is equal to the system response to that particular excitation.Therefore, the determination of the impulsive or indicial response of an aerodynamic system is sufficient to know all the unsteady behavior of that system.Such approach is very interesting for engineering situations in which the aerodynamic behavior is coupled with other dynamic systems, such as in aeroelasticity and aeroservoelasticity.By doing so in the aeroelastic problem, for example, the governing equations of the structural dynamic problem may be solved separately from the aerodynamic equations.In this case, the aerodynamic influence is represented in terms of the aerodynamic forces, which, in turn, are known for an arbitrary deformation of the structure [13,14].
On the other hand, it has been very common in the literature to look at CFD solvers as mere approximations to continuous-time systems.This traditional form of reasoning has led many authors [14,15,16,17] to justify the use of smooth pulse excitations, since the theoretical continuous-time impulse and indicial excitations would not be numerically feasible.However, Silva [18] demonstrates that, once the governing equations have been discretized, the resulting numerical scheme is actually a discrete-time system, which has its own properties and peculiarities.Hence, Ref. [18] suggests the use of the equivalent discrete-time excitations, i.e., unit sample and discrete step [19].The unit sample sequence is given by and the discrete step sequence is or, equivalently, The unit sample and discrete step sequences formally hold properties very similar to the ones attributed to the impulse and step functions, respectively.Namely, if a linear timeinvariant discrete-time system is subjected to a unit sample excitation, then the corresponding response will contain all the information about the system, and the response to every other input is given by the convolution sum [19], where in[n] is a generic input sequence, out[n] is the respective response, and h[n] is the system response to a unit sample input.It is very important to emphasize that, although the convolution sum given in Eq. 19 resembles the convolution integral [12], it is not an approximation to such integral, but a formally defined discrete operation.Moreover, the unit sample rigorously excites uniformly the complete frequency domain.The discrete step response also characterizes a discrete-time system since one can reproduce the unit sample response from it.As demonstrated in Ref. [20], the backward difference of the discrete step response yields the unit sample one, i.e., where S[n] designates the system response to a discrete-step input.Thus, theoretically, it would be more convenient, and even computationally cheaper, because the transient solution should die out more rapidly, to acquire CFD results submitting the body to either unit sample or discrete step type perturbations than to smooth inputs.In addition, Ref. [3] states that simulations carried out with the discrete step excitation tend to be less influenced by the simulation parameters than those using the unit sample.Another possibility, also described in Ref. [1], consists in employing exponentiallyshaped signals which attempt to approximate smooth impulsive functions.As the exponentiallyshaped input is not a unit sample excitation, the real unit sample response, in this case, is evaluated using a well-known property of the convolution theorem [19], where h[n] represents the time response to a unit sample excitation, and g [n] is the response to the sampled exponentially-shaped excitation, f p [n].The sequences in capital letters are the discrete Fourier transforms of the corresponding sequences in lower case letters.Therefore, the frequency domain unit sample responses can be obtained by dividing the discrete Fourier transform (DFT) of the responses to the sampled exponentially-shaped pulse by the DFT of the input sequence.Although the input is not the exact unit sample excitation, it is capable of exciting the reduced frequencies of interest in aeroelastic studies.
The corresponding frequency domain points resulting from the all previously described discrete procedures can be written as [14] Here, f designates the frequency, a ∞ is the undisturbed flow sound speed, c is the chord length, ∆t is the time step of the aerodynamic solution, and N is the maximum number of discrete iterations in which this solution is known.Moreover, the bar designates a dimensionless time step, defined as Equation 22, rewritten in terms of the reduced frequency based on the half-chord length (b), stands as where ω is the circular frequency, U ∞ is the undisturbed flow speed, and M ∞ is the undisturbed flow Mach number.

SYSTEM IDENTIFICATION
As discussed, indicial and impulsive responses are typically used for the identification of the dynamic characteristics of many systems.On the other hand, the use of a high-fidelity CFD solver is justified for the transonic flight regime due to its inherent nonlinear behavior.However, as shown in Ref. [1], there is an amplitude range in which the aerodynamic response due to the unsteady motion can be considered linear with respect to the mode shapes, even in the transonic regime.Therefore, assuming that the motion amplitude is sufficiently small, one can formulate the problem as a linear dynamic system in the frequency domain and, as such, write the input-output relation as where G ij is the transfer function of the input j to the output i.Alternatively, G ij could be seen as the generalized response in i-th mode due to an impulsive input in the j-th mode.
The traditional approach to identifying this system is to provide an excitation signal in one of the inputs and measure the outputs.Each transfer function is, then, obtained by the division, in the frequency domain, of the output by the input.This, however, is a single-inputmultiple-output (SIMO) approach to the multiple-input-multiple-output (MIMO) problem and requires, in the applications of interest here, a number of unsteady simulations equal to the number of modes considered in the system.If, on the other hand, one multiplies the convolution sum equation for a given output by one of the inputs, considers the expected value and uses the Wiener-Khintchin relations [21], one can rewrite any line of Eq. 27 as where P y k u i is the cross power spectral density of the output y k and the input u i and P u j u i is the cross power spectral density of the inputs u j and u i .The aforementioned steps would not be of much help unless the input signals are uncorrelated.If such, then, any cross power spectral density of two different inputs would be zero.Hence, assuming that uncorrelated input signals are employed, the system can be excited simultaneously in all of the modes considered in the problem and the transfer functions can be estimated by One possible set of uncorrelated input signals can be obtained with the use of Walsh functions [6].These are a set of block step functions which form an orthogonal basis of the square-integrable functions.As noted by Silva [6], this family of functions has a similarity to step inputs and, therefore, embodies the impulsive nature of step inputs with regard to frequency bandwidth.Each function takes positive and negative unitary values and each step block period is an integer multiple of a 2 n division of the function length.Each one of these functions may be seen as a line, or a column, of a Hadamard [22] matrix of order 2 n .In order to construct these matrices, the authors have employed the method proposed by Sylvester [23], in which a Hadamard matrix H 2 n is defined recursively by where H 1 is the unitary identity matrix.Each element of a row is, then, assigned as the magnitude of the corresponding step block in a Walsh function.For obvious reasons, one must construct a Hadamard matrix of order greater or equal to the number of inputs desired for the system.Furthermore, one should be careful when choosing the length of the Walsh function so that each block period is capable of exciting the relevant frequencies in the system.

MESH GENERATION AND MOVEMENT
The meshes used in the present work are generated with the commercial grid generator ICEM CFD c ⃝ .Figure 1 shows two views of the mesh around a NACA 0012 airfoil which is employed to obtain the results discussed here.Although other meshes have been used in the  present work, the one shown in Fig. 1, which has about 16000 triangular control volumes and just over 290 points along the airfoil, provides adequate resolution of the phenomena of interest.Moreover, the experience developed in the work reported in Ref. [1] was also helpful in defining the computational meshes for the present case.Unsteady calculations involve body motion and, therefore, the computational mesh should be somehow adjusted to take this motion into account.Two approaches for moving the computational grid mesh are available in the present CFD code.The first one follows the ideas presented in Ref. [24], according to which the sides of the triangles that constitute the mesh are modeled as springs with stiffness constants inversely proportional to their lengths.The other way of accounting for the motion of the body is to rigidly move the mesh together with the airfoil.In the present work, only results using the second approach are presented, because a standard typical section model is of interest.Therefore, it is simpler and less expensive to move the mesh rigidly with the airfoil.For the present application, as demonstrated in Ref. [1], the results obtained with either approach are equivalent.It is important, however, to notice that, in this approach, the far-field boundary is no longer kept fixed, but it moves together with the rest of the mesh.Hence, the far-field boundary conditions need to be adjusted in order to account for this motion [1].Moreover, regardless of the approach used to move the grid, mesh velocity components can be calculated as The new cell areas are determined employing the geometric conservation law concept [25], where the total area conservation is imposed in a similar manner to the other flow conserved variables.

RESULTS
The test case considered throughout this paper involves a NACA 0012 airfoil at M ∞ = 0.8 and α 0 = 0.In order to perform unsteady calculations, an initial steady-state solution must be obtained prior to initiating the unsteady CFD runs.Such solution is obtained by the same CFD code, running in steady mode, i.e., with a variable time step method, which is converged to machine zero.Since the numerical method implemented in the code provides steady-state solutions which are independent of the time step, such converged solution can, now, be used as the initial condition for all the unsteady calculations in the present work.
Pressure coefficient results for the initial steady-state solution for the present case are shown in Fig. 2. As one can see, the present results are capable of representing the physical phenomena expected for the flows of interest, including the strong shock waves that occur over the airfoil.Figure 2   Comparisons with the numerical results available in Ref. [26] (Lit.Num.) and experimental ones found in Ref. [27] (Lit.Exp.) are also shown in Fig. 2 (b).The current results present a small oscillation in the pressure coefficient distribution in the region prior to the shock wave and in the vicinity of the trailing edge.These discrepancies, however, do not affect the overall adequateness of the solution and this is rather confirmed by the good agreement between the current solution and that of Ref. [26].Moreover, for a NACA 0012 airfoil with α 0 = 0, one would expect zero lift and moment, but numerical errors lead to small residual values.Furthermore, although one cannot expect Euler results to exactly match experimental measurements, Fig. 2 (b) shows that, except for the aforementioned regions, the numerical pressure coefficient distribution is very close to the experimental one.Therefore, one may consider the current computational results to be in good agreement with the experimental and numerical data in the literature.
From the converged stationary solution, the mesh is rigidly moved in a prescribed pattern and a total of 100, 000 time steps of unsteady flow are computed with a constant ∆ t of 0.003 dimensionless time units.Discrete step inputs and Walsh functions are considered in order to prescribe the airfoil movement.The step inputs are used in a mode-by-mode fashion whereas the Walsh functions are used to simultaneously excite all the system modes.Furthermore, regardless of the motion being prescribed, the maximum amplitudes considered here are 0.000001 c for the plunging mode and 0.0001 deg.for the pitching degree of freedom.The reason for choosing such low amplitudes, as discussed earlier, is to remain within the "quasilinear" region around the nonlinear steady solution and to allow the CFD code to accurately propagate the disturbances from the discrete motion [3].Further discussion of the effects of the amplitudes of modal motions, in the context of the present CFD solver, are presented in Ref. [1].
Moreover, the Walsh input considered here was modified in order to present a blank region at the end of the unsteady run as illustrated in Fig. 3.This zero-motion region attempts to guarantee that the computation starts and ends with a solution without any perturbations from the prescribed motion.The size of this blank region also influences the size of the window function to be considered and, therefore, the resolution in the frequency domain.
Many configurations have been studied by the authors with sets of different sizes of zero-motion regions and locations, some of which are presented in Ref. [28].The authors have chosen to present the results from the set illustrated in Fig. 3 in the present paper because such results have a better resolution in the frequency domain than those presented in Ref. [28].Either way, the simulations with the set of Walsh functions require only one unsteady CFD run, whereas those with a discrete step require two, i.e., one for each modal movement.Figure 4 shows the time history of the solution in terms of the aerodynamic coefficients C ℓ and C m for a discrete step (DS) input in plunge.C ℓ and C m denote the lift and moment coefficients, respectively, and the latter refers to the quarter-chord point and it is positive in the nose-up direction.Although the total integration time is 300 dimensionless time units, only the initial part of the solution is shown in order to emphasize the most relevant transient response.A similar simulation is performed for a pitching motion and the corresponding time responses are given in Fig. 5.Although not specifically shown in the figures, these unsteady    results were compared to similar calculations presented in Ref. [1] and they are in very good agreement.
A better visualization of the perturbations created by the different motions, for the present unsteady simulations, can be inferred from a plot of the residue time histories.It is true that the concept of residue does not really apply to unsteady calculations, but nevertheless it is helpful in the present case.Figure 6 presents the time history of the L ∞ norm of the density residues for the unsteady calculations using the set of Walsh function inputs here considered.When the body is moved according to the prescribed pattern shown in Fig. 3, a  seen in the residue curve.As the aerodynamic perturbations from the motion die out, there is a marked decrease in the residue of the calculation, as indicated in Fig. 6.This figure also shows that there is a significant decrease of the perturbations prior to the next motion being prescribed in the set of Walsh functions used.The fact that the residues quickly decay after a given perturbation is an indication that the amplitudes of motion here used are adequate and the CFD code is capable of handling the perturbation velocities so generated.
Reference [1] presents frequency domain responses for the transonic configuration in question when the airfoil is submitted to both plunging and pitching motions.These responses are obtained numerically with a solver very similar to the one used in the present work.The present authors attempt to reproduce these solutions using discrete step and the set of Walsh function inputs in order to validate the present capability of obtaining the aerodynamic transfer functions.Hence, initially, results for the impulsive response in the frequency domain, for both C ℓ and C m , obtained using the concept embodied in Eq. 29 are compared to the transfer functions obtained using the approach described in Ref. [1].The latter approach corresponds to the application of Eq.21 to the fast Fourier transform (FFT) of the input and output signals, for instance, to an impulsive excitation of one of the modes of the problem.Figure 7 compares, then, the C ℓ and C m transfer functions in the frequency domain obtained by the two different approaches, using a discrete step (DS) input in plunge.The aerodynamic coefficients are shown in terms of real and imaginary parts of the frequency domain transfer functions.A similar comparison can be seen in Fig. 8 for the a DS input in pitch.It is clear from these figures that the transfer functions obtained from the two approaches are nearly identical.The plots only show the lower portion of the frequency range, because this is the range of reduced frequencies of interest for flutter analyses.Although not shown in the figures, similar good agreement also holds for the entire frequency range.It should be observed that the maximum frequency which can be computed with this approach is associated to the sampling frequency of the CFD code, i.e., it is equal to the frequency defined by the unsteady calculation time step, as indicated in Eq. 22.However, only half of this spectrum has physical relevance because the time history results are not periodic.
The validation of the procedure implemented, that uses the Walsh functions, can be achieved by comparing the transfer functions obtained with a single unsteady CFD run with the ones previously calculated using the mode by mode approach.Figure 9 shows the comparison of the G C ℓ ,h and G Cm,h transfer functions obtained using the set of Walsh functions and the corresponding transfer functions calculated by Eq. 29 from a discrete step in plunge.The latter results are the same already presented in Fig. 7.A similar comparison for the G C ℓ ,α and G Cm,α transfer functions is shown in Fig. 10, where the reference result comes from a discrete step in pitch as presented in Fig. 8.As before, only the lower portion of the frequency range is shown in Figs. 9 and 10.Except for some isolated discrepancies for very low values of reduced frequency, particularly for the transfer functions due to pitching motion, the agreement of all curves is excellent.It should be noticed that only an extreme enlargement of the plots near the origin would allow a visualization of such issues.These discrepancies occur due to a reduction in frequency resolution caused by the windowing function that must be applied when calculating the power spectral densities.The window length is defined by the size of an information block, according to the definition of the set of Walsh functions being used, and this length is constant for a given power spectral density calculation.In other words, the information block is defined by the length of the input string where no variations occur.Hence, for the calculations using the WF input set, the window length is N/3 (rounded down), where N is the total number of samples in the time history of the aerodynamic coefficients.
The last stage in preparing the data generated by the CFD calculations in order to be able to perform aeroelastic stability analyses in the frequency domain consists in approximating the string of transfer function values by an interpolating polynomial.In the present work, the interpolating polynomials proposed by Eversman and Tewari [29] without any provision for the treatment of repeated, or very close, poles are used.These are referred to as the first form of the Eversman and Tewari polynomials in Ref. [2].Here, all calculations considered these interpolating polynomials with the addition of 6 poles.Moreover, as described in Ref. [2], the interpolating polynomials are obtained by a least-squares fitting of the transfer functions generated from the CFD results.In the present case, this fitting is performed for the reduced frequency range of 0 ≤ κ ≤ 3.
Figure 11 presents the comparison for the results for the C ℓ transfer functions.The transfer functions are shown for the usual 0 ≤ κ ≤ 3 reduced frequency range, in order to maintain consistency with previous figures, and in terms of an enlargement of the very low frequency region, for κ ≤ 0.3.Moreover, for all figures, the real and imaginary parts of the approximated transfer functions given by the interpolating polynomials, obtained from the transfer functions generated using the WF inputs, are compared to the transfer functions obtained directly from the CFD solution using the DS approach.It should be observed that these approximated transfer functions are the data in which a flutter stability analysis is based on.A similar comparison for the C m transfer functions is shown in Fig. 12.
As one can see in Figs.11 and 12, in general, there is good agreement between the sets of results from the approximating polynomials and the original CFD data.Discrepancies occur in the extremely low portion of the frequency range analyzed and in regions in which the there is a rapid variation in the CFD data.Such behavior is consistent with the results reported in Ref. [1] and, therefore, it has no relation with the procedure for generating the transfer functions.It is simply a statement that the polynomial is an approximation of the actual transfer function and, hence, it has difficulties when there are rapid changes in the frequency domain response.

CONCLUDING REMARKS
The present paper discusses different forms of using CFD techniques for generating the aerodynamic operator for aeroelastic stability analysis.The major motivation for the work is the development of the capability of generating unsteady aerodynamic transfer functions in an efficient way, such that CFD results could be readily incorporated into aeroelastic analysis   procedures.In other words, the objective is that the aerodynamic information for aeroelastic stability analyses should cost, at most, two CFD runs, i.e., one steady and one unsteady solution.Moreover, a discussion is also presented regarding different forms of exciting the dynamic system in order to achieve such CFD computational cost reduction.It should be also clear that the major interest of the approaches here discussed lies in the calculation of transonic flutter boundaries, in which the nonlinearities in the representation of the aerodynamic operator are relevant for the accurate prediction of the phenomena of interest.Furthermore, since a linearization with respect to the modal displacements is implicit in the present approach, the calculation procedure is adequate for flutter analyses, but it may not be suitable for limit-cycle oscillation calculations if the structural displacements become larger.
Results are presented for a NACA 0012 airfoil-based typical section model, for freestream conditions M ∞ = 0.8 and α 0 = 0. Steady aerodynamic results are compared to the literature data and they are shown to be in good agreement.From this stationary condition, the system is perturbed in a prescribed pattern in order to obtain unsteady CFD results.Initially, the non-parametric system identification technique, which employs the calculation of power spectral densities, is validated by comparing the aerodynamic transfer functions obtained with such technique with those calculated from impulsive or indicial excitations in each mode sep-   arately.Afterwards, the use of a simultaneous excitation signal is employed and it is shown that the system identification routine is capable of splitting the output into the contribution of each individual mode to the corresponding aerodynamic transfer functions.The values in the frequency domain originated by the simultaneous excitation signal match those obtained in a mode by mode fashion.Finally, it is shown that the discrete strings of transfer function values can be conveniently approximated using selected interpolating polynomials.These interpolating polynomials contain, therefore, the aerodynamic states, which can be used for direct aeroelastic stability analyses in the frequency domain.Since the approximating polynomials here obtained are in good agreement with similar data in the literature, this means that the procedure implemented accomplishes the desired goal of obtaining the aerodynamic operators for aeroelastic analyses with a single unsteady CFD calculation.
Complete mesh, with 8323 nodes and 16232 volumes.

Figure 1 .
Figure 1.Mesh around NACA 0012 profile with 292 wall points and 50 chord-lengths of distance between the body and the far-field.
(a) presents pressure coefficient contours for the flowfield near the Pressure coefficient distribution over airfoil.

Figure 3 .
Figure 3. Normalized motions prescribed by the Walsh function (WF) input.
C m time history.

Figure 4 .
Figure 4. Initial part of the time response to a DS in plunge.
C m time history.

Figure 5 .
Figure 5.Initial part of the time response to a DS in pitch.

Figure 6 .
Figure 6.Maximum density residue throughout the unsteady calculations.
C ℓ frequency domain response.
C m frequency domain response.

Figure 7 .
Figure 7. Low reduced frequency impulsive response of a NACA 0012 airfoil at M ∞ = 0.8 and α 0 = 0 to a DS input in plunge.Comparison between current method and Ref. [1].
C m frequency domain response.

Figure 8 .
Figure 8.Low reduced frequency impulsive response of a NACA 0012 airfoil at M ∞ = 0.8 and α 0 = 0 to a DS input in pitch.Comparison between current method and Ref. [1].
C ℓ frequency domain response.
C m frequency domain response.

Figure 9 .
Figure 9. Low reduced frequency G C ℓ ,h and G Cm,h transfer functions for a NACA 0012 airfoil at M ∞ = 0.8 and α 0 = 0 obtained from a DS input in plunge and from the Walsh functions.
C ℓ frequency domain response.
C m frequency domain response.

Figure 10 .
Figure 10.Low reduced frequency G C ℓ ,α and G Cm,α transfer functions for a NACA 0012 airfoil at M ∞ = 0.8 and α 0 = 0 obtained from a DS input in pitch and from the Walsh functions.
solution) Im (CFD solution) Re (Interpolated from WF) Im (Interpolated from WF) (a) G C ℓ ,h , normal region of study.solution) Im (CFD solution) Re (Interpolated from WF) Im (Interpolated from WF) (b) G C ℓ ,h , detail for very low reduced frequency.

Figure 11 .
Figure 11.G C ℓ ,h and G C ℓ ,α transfer functions at low reduced frequencies.
solution) Im (CFD solution) Re (Interpolated from WF) Im (Interpolated from WF) (a) G Cm,h , normal region of study.solution) Im (CFD solution) Re (Interpolated from WF) Im (Interpolated from WF) (b) G Cm,h , detail for very low reduced frequency.
solution) Im (CFD solution) Re (Interpolated from WF) Im (Interpolated from WF) (c) G Cm,α , normal region of study.solution) Im (CFD solution) Re (Interpolated from WF) Im (Interpolated from WF) (d) G Cm,α , detail for very low reduced frequency.

Figure 12 .
Figure 12.G Cm,h and G Cm,α transfer functions at low reduced frequencies.