Estimation of Unsteady Aerodynamic Models from Dynamic Wind Tunnel Data

NATO UNCLASSIFIED + SWE

Estimation of Unsteady Aerodynamic Models from Dynamic Wind Tunnel Data

Patrick C. Murphy NASA Langley Research Center, Hampton, VA, 23681-2199 Patrick.C.Murphy@nasa.gov Vladislav Klein National Institute of Aerospace, NASA Langley Research Center, Hampton, VA, 23681-2199 Vladislav.Klein@nasa.gov

ABSTRACT

Demanding aerodynamic modelling requirements for military and civilian aircraft have motivated researchers to improve computational and experimental techniques and to pursue closer collaboration in these areas. Model identification and validation techniques are key components for this research. This paper presents mathematical model structures and identification techniques that have been used successfully to model more general aerodynamic behaviours in single-degree-of-freedom dynamic testing. Model parameters, characterizing aerodynamic properties, are estimated using linear and nonlinear regression methods in both time and frequency domains. Steps in identification including model structure determination, parameter estimation, and model validation, are addressed in this paper with examples using data from one- degree-of-freedom dynamic wind tunnel and water tunnel experiments. These techniques offer a methodology for expanding the utility of computational methods in application to flight dynamics, stability, and control problems. Since flight test is not always an option for early model validation, time history comparisons are commonly made between computational and experimental results and model adequacy is inferred by corroborating results. An extension is offered to this conventional approach where more general model parameter estimates and their standard errors are compared.

NOMENCLATURE

A, B, C = transfer function coefficients J(  ) = sum of squared residuals A , B = Fourier coefficients j j k = reduced frequency, 2 / f V  a , b , c = indicial function parameters = half chord or half span, ft b = wing span, ft m = No. of harmonics in Fourier expansion C = rolling-moment coefficient l N = number of data points C = normal-force coefficient N p, q = roll and pitch rates, rad/sec c = mean aerodynamic chord, ft R = squared multiple correlation coefficient , C C = in-phase components of C s = estimated standard error a a a   t = time, sec , C C = out-of-phase components of C a a a q p V = airspeed, ft/sec D = time domain differential operator y = output variable d , d = constants z = measured variable 0 1 f = frequency, Hz Senior Research Engineer, Dynamic Systems & Control Branch, Mail Stop 308, Associate Fellow.

Professor Emeritus, Dynamic Systems & Control Branch, Mail Stop 308, Associate Fellow.

RTO-MP-AVT-189 PAPER NBR - 1 NATO UNCLASSIFIED + SWE NATO UNCLASSIFIED + SWE Estimation of Unsteady Aerodynamic Models from Dynamic Wind Tunnel Data  = angle of attack, rad or deg  = dummy integration variable  = mean angle of attack, rad or deg 1 V    = non - dimensional time constant,     = sideslip angle, rad or deg b    = state variable  = residuals  = vector of unknown parameters  = angular frequency, rad/sec , 2 f    = variance Subscrip ts Aerodynamic Derivatives A = amplitude C C   l N ( ) , ( ) C C C C       l l N N a = N or l         Superscripts C C   l N ( ) , ( ) C C C C       l l N N ~ = complex variable p p q q V p V q   ~* = complex conjugate

1.0 INTRODUCTION

ince the early days of flight, aerodynamicists have focused on the problem of finding adequate

S

aerodynamic models with good prediction capability. The problem has become more acute over time with modern military aircraft requirements expanding the flight envelope, demanding more capability in nonlinear unsteady flight regimes, and civilian aircraft requirements for safety, demanding more capability in similar adverse aerodynamic conditions. Over time researchers have endeavored along two paths to address this nonlinear unsteady modelling problem: a numerical path using high fidelity Computational Fluid Dynamics (CFD) technology and an experimental path using direct measurements of aircraft responses in flight or scale models in wind tunnels along with System Identification (SID) technology to extract adequate mathematical models from the measured data. Advances in both numerical and experimental technologies have created opportunities for the combination of these technologies to make significant improvements in handling more difficult aircraft modelling problems.

In order to capitalize on these opportunities, model validation and identification should be carefully considered. Although it is not common practice for wind tunnel or CFD results to be presented with associated error bounds for validation comparisons, some efforts, Refs. [1-4], are being made to improve this shortcoming. Model identification, considered in this paper, is a research area that involves a number of disciplines in modelling and measurement science. The subject of aircraft system identification, in particular, ® is documented in Ref. [5] and MATLAB software is included for a variety of system identification methods commonly used in aerospace engineering practice. The identification process includes a number of key steps such as, experiment design, model structure determination, parameter estimation, and model validation. These steps lead to specific demands on test techniques and test facility capabilities. The demand for more general, higher fidelity, models has highlighted some limitations of the conventional approaches. For example, in ground-based dynamic wind-tunnel experiments where single frequency sinusoidal inputs are typically used, wide-band inputs have been shown to be more effective and less costly since fewer runs are required, Ref [6- 7]. Advancement in these areas is hampered by the expense of developing advanced test facilities capable of more general motions. CFD offers an opportunity to ameliorate a number of these limitations.

As an initial step toward improved integration of CFD and SID technologies, this paper presents identification methods used at NASA Langley Research Centre that can be applied when nonlinear unsteady aerodynamic PAPER NBR - 2 RTO-MP-AVT-189 NATO UNCLASSIFIED + SWE NATO UNCLASSIFIED + SWE Estimation of Unsteady Aerodynamic Models from Dynamic Wind Tunnel Data responses are present. The model identification problem for conventional one-degree-of-freedom motion in wind tunnel forced-oscillation experiments is described. Mathematical model structures are presented that include conventional static and rotary dynamic terms but are extended by including indicial functions to represent unsteady responses. Although other model structures have been studied, the authors’ goal is to maintain the conventional stability and control derivative structure to capitalize on the large engineering knowledge base built around that structure. Model parameters describing aerodynamic properties of the system under test can be estimated using a least-squares principle. The resulting estimators form linear or nonlinear regression problems that can be formulated in the time or frequency domain. Five parameter estimation methods are described and demonstrated on experimental data from wind tunnel and water tunnel experiments. These methods are harmonic analysis, nonlinear regression, two-step linear regression, equation error, and output error. Model validation of these techniques is demonstrated with examples using the NASA Generic Transport Model (GTM), Fig. 1a, and two scale models of the F-16XL: a 2.5% scale water tunnel model, Fig. 1b, and an 18% scale wind tunnel model not shown.

2.0 MATHEMATICAL MODEL STRUCTURES

Estimation of aerodynamic parameters from wind tunnel data requires that a mathematical model of the aircraft is postulated. The mathematical model includes both the aircraft equations of motion and the equations for aerodynamic forces and moments, known as the aerodynamic model equations. For most of the nominal operating envelope, where the aerodynamic flows are attached and quasi-steady in behaviour, aerodynamic model equations are usually developed using time-invariant linear terms (stability and control derivatives). In off-nominal cases these terms are generalized by including time-dependent and nonlinear terms. Results from numerous studies, such as Ref. [6], have shown the dependence of aerodynamic parameters on frequency. This dependency contradicts the basic assumption that stability derivatives are constants. Because the effect of frequency is related to unsteady aerodynamics, the aerodynamic model equations will be formulated in terms of indicial functions or “step - response functions” , Refs. [8-9]. As an example, the normal force coefficient is considered as t t c ( ) (0) ( ) ( ) ( ) ( ) C t C C t d C t q d               N N N N q  2 V 0 0 (1) where ( ) C t and ( ) C t are the indicial functions. For further development of Eq. (1) it will be assumed N N  q that a) only the increments to steady conditions are considered b) the effect of ( ) q t on the normal force can be neglected c) the indicial functions can take the form of a simple exponential b t  ( ) (1 ) C t a e c    N  b t  ( ) C ae    N  where the indicial function is separated into a steady and deficiency function, Ref. [10]. Consequently, this RTO-MP-AVT-189 PAPER NBR - 3 NATO UNCLASSIFIED + SWE NATO UNCLASSIFIED + SWE Estimation of Unsteady Aerodynamic Models from Dynamic Wind Tunnel Data model can be simplified as t c ( ) b t    ( ) ( ) ( ) ( ) ( ) ( ) C t C t C q t a e d           N N N q  2 V (2a) or in operator form as a ( ) ( ) ( ) ( ) ( ) C t C t C q t D        N N N q  V D b  (2b) By introducing t ( ) b t    ( ) ( ) t e d       (3) the state-space form of Eq. (2) is ( ) ( ) ( ) t b t t       ( ) ( ) ( ) ( ) ( ) ( ) C t C t C q t a t        N N N q  V (4) Applying the Laplace transform to Eq. (4), the transfer function for the coefficient is obtained as ( ) C s As Bs C   N  ( ) s s b   (5) where s is the Laplace transform parameter and ( ) A C   N q V ( ) ( ) B C a b C      1 N N q  V (6) ( ) C b C   1 N  is the half chord or half span as required. A different model structure suitable for the analysis of oscillatory data, with amplitude  , follows from the steady solution of Eq. (2) as A ( ) sin( ) cos( ) C t C t kC t       N A N A N q  (7) PAPER NBR - 4 RTO-MP-AVT-189 NATO UNCLASSIFIED + SWE NATO UNCLASSIFIED + SWE Estimation of Unsteady Aerodynamic Models from Dynamic Wind Tunnel Data where / k cf V   is the reduced frequency. For the conventional model form with linear aerodynamics, ( ) C t is assumed to be functions of , , , q   and q . The in-phase and out-of-phase components, C and N N  C , can then be expressed in terms of steady flow and acceleration terms N q ( ) C C k C    N N N q   (8a) ( ) C C C    N N N q q  (8b) These acceleration terms are the conventional equivalent to the unsteady terms. In order to explain the variation of the in-phase and out-of-phase components with frequency these components were expressed in terms of deficiency functions, as in Eq. (2). In Refs. [10-11], the model for these components has the form 2 2 k  ( ) C C a    N N 2 2   1 k   (9a)  ( ) C C a    N N 2 2 q q 1 k   (9b) where 2 / V cb   is an aerodynamic time constant. In the aerodynamic model equations presented here there 1 1 are six unknown parameters: C , C , a , b (or  ), C , and C . Techniques for estimating these N N 1 1 N N  q  q unknown parameters are provided in the next section.

3.0 ESTIMATION METHODS

The parameters in the aerodynamic model equations can be estimated in various ways, e.g., U.S. Air Force DATCOM, based on an extensive set of rules for aerodynamic dependencies and aircraft geometry, can produce basic aerodynamic parameters but these are limited to basic geometries and steady-flow dynamics.

Alternatively, analytic methods such as CFD, where general physics-based equations define the mathematical model, can provide general aerodynamic solutions but application to stability and control problems is still an area of research. Although these methods work well to varying degrees, in particular for low angle of attack and low rotational rate cases, the best aerodynamic predictions are still obtained using experimental methods.

Flight experiments provide the most direct measurements of aerodynamic behaviours but this is the most expensive approach and not available in early design phases. The wind tunnel provides the next best approach with current technology; however, this approach has a number of difficulties related to low speeds, similitude scaling, and tunnel effects. In this paper the aerodynamic parameters will be estimated from wind tunnel dynamic tests to demonstrate modeling methodology that can also be applied to appropriate CFD simulations.

Five methods for estimation of parameters in aerodynamic model equations are presented in this paper: Harmonic Analysis, Nonlinear Regression, Two-Step Linear Regression, Equation Error, and Output Error.

They are based on application of a least-squares principle to experimental data assuming that the aerodynamic RTO-MP-AVT-189 PAPER NBR - 5 NATO UNCLASSIFIED + SWE NATO UNCLASSIFIED + SWE Estimation of Unsteady Aerodynamic Models from Dynamic Wind Tunnel Data model structure is known. Wind tunnel dynamic tests usually include forced oscillations about the aircraft body axes at selected angles of attack, sideslip, frequencies, amplitudes, and Reynolds numbers. During these tests the inputs (Euler angles) and the outputs (aerodynamic coefficients of forces and moments) are recorded.

All the estimation methods presented are demonstrated using data from wind tunnel experiments using the NASA GTM and 18% scale F-16XL, and from water tunnel experiments using a 2.5% scale F-16XL model.

3.1 Harmonic Analysis A method of harmonic analysis, Ref. [12], is applied to measured aerodynamic coefficients to allow estimation of the in-phase and out-of-phase components. For the development of this method it is assumed that a periodic function ( ) ( 2 ) y t y t    with the period  2 is represented by a series of discrete values ( ) y i at 2 i  ( ) , 1,2,..., t i i N   N It is assumed that ( ) y i can be approximated by a trigonometric series m m

     i j B i j A A i y sin cos  

 

j j 0 0 0   j j 1 1 (10) 2  where   . This series uses an orthogonal-basis of harmonic sinusoids. It is further assumed that the N measurements of ( ) y i are obtained as ( ) ( ) ( ) z i y i v i   (11) where ( ) z i are the measured values and ( ) v i is white measurement noise with zero mean and constant variance  . The parameters , A and B , are Fourier coefficients that can be estimated from the measurements by A 0 j j minimizing the least squares (LS) criterion N ( ) J z i y i     

    

LS   1 i  (12) where  is the vector of unknown parameters, ( A , A , B , …, A , B ). Then the LS estimates of parameters in 0 1 1 m m (10) are PAPER NBR - 6 RTO-MP-AVT-189 NATO UNCLASSIFIED + SWE NATO UNCLASSIFIED + SWE Estimation of Unsteady Aerodynamic Models from Dynamic Wind Tunnel Data N ˆ

   i z A



N  i 1 N ˆ

   i j i z A cos 

 j 0

N  i 1 N 2 (13) ˆ

   i j i z B sin 

 j 0

N  i 1 The parameter variances of these estimates are 2 2 ˆ s A  

 

N 2 2 2 ˆ ˆ s A s B   

   

j j (14) N for all j. The estimate of variance   is N ˆ s z i y i    

    

  N 1 i  (15) ˆ where   i y follows from Eq. (10) by replacing the parameters by their estimates.

The adequacy of the model given by Eq. (10) can be assessed by the coefficient of determination that indicates how much variation in the data is explained by the model. R is expressed as N ˆ

    z i y i  

 

2 2 1 i  1 , 0 1 R R     N z i z  

   

(16) 1 i  where z is the mean value of z i .

 

For the pitch oscillation case, considering only the first harmonic and a model with linear aerodynamics, the in-phase and out-of-phase components of ( ) C t , can be expressed in terms of coefficients A and B , 1 1 N respectively, as B A 1 1 and C C   N N q  k   A A (17) where / k cf V   is the reduced frequency. Eq. (17) can be easily modified for the roll and yaw oscillation case, Ref. [12].

An example of harmonic analysis performed on pitch oscillatory data is given in Fig. 2 for the NASA Generic RTO-MP-AVT-189 PAPER NBR - 7 NATO UNCLASSIFIED + SWE NATO UNCLASSIFIED + SWE Estimation of Unsteady Aerodynamic Models from Dynamic Wind Tunnel Data Transport Model (GTM). Ordinate values were removed in order to maintain proprietary agreements. In this figure the in-phase and out-of-phase components, and the coefficient of determination are plotted against angle of attack for seven frequencies and one selected amplitude. This experiment was conducted in the NASA Langley 14x22 Wind Tunnel. Model geometry for the GTM is given in Fig. 1a.

3.2 Nonlinear Regression Harmonic analysis results shown in Fig. 2, indicate that the in-phase component, C , does not vary with N  frequency and, although not shown, its values agree with C determined from steady measurements. On the N  other hand, the out-of-phase component, C , shows variation with frequency for f ≤ 0.43 Hz and 10° ≤  ≤ N q 50°. Values of R indicate that a linear aerodynamic model may be adequate for  ≤ 35°. In this case, three unknown parameters, , , and C a  , in Eq. (9b), can be estimated by a nonlinear regression method, 1 N q suggested in Ref. [13], using the frequency dependent measurements.

A general model for this method can be formed as ( ) [ ( ), ] ( ) 1,2,..., z j g x j v j j m     (18) where ( ) x j is a vector of regressors computed from measured data at the jth data point, g is a nonlinear function of , and unknown parameters are given by the vector,  . The least-squares estimator can be ( ) x j obtained by minimizing the sum of squared errors m ( ) { ( ) [ ( ), ]} J z j g x j      NR 1 j  (19) To demonstrate this estimation method the out-of-phase component, C , shown in Fig. 2, will be used as an N q example. In this case the regression equation is  ( ) ( ) C j C a    N N 2 2 q q 1 ( ) k j   (20) To check the regression, parameter estimates can be used to compute an estimated C using Eq. (12).

N q Measured C and estimated C are plotted in Fig. 3. These results show a good fit to the measured data.

N N q q For a demonstration of model prediction capability and model validation, Fig. 4 shows the variation of C N with  at f = 0.86 Hz,  = 10°, and two nominal values of  = (14°, 26°). In both cases the shape of the A 0 predicted data forms a regular ellipse reflecting the underlying linear aerodynamic model structure. The measured C shows some deviation from a regular ellipse, or linear behaviour, in the pre-stall and stall N regions.

PAPER NBR - 8 RTO-MP-AVT-189 NATO UNCLASSIFIED + SWE NATO UNCLASSIFIED + SWE Estimation of Unsteady Aerodynamic Models from Dynamic Wind Tunnel Data 3.3 Two-Step Linear Regression An aerodynamic model of an aircraft performing a one degree-of-freedom oscillatory motion about one of its body axes can be formulated in terms of the in-phase and out-of-phase components, Ref. [12]. Using subscript a to represent the appropriate force or moment, the model can be written for rolling motion as ( )sin sin C C af      1 a a   sin C C af    

 

0 a a p p (21a) for pitching motion as

  C C af   

1 a a   C C af   

 

0 a a q q (21b) and for yawing motion as cos cos C C af     

 

1 a a   cos C C af    

 

0 a a r r (21c) where 2 2 k  f  2 2 1 k    f  (22) 2 2 1 k   With the expression, 2 2 1 k  1   2 2 2 2 1 1 k k     1 1 (23) equations (21a), (21b) or (21c) can be rearranged into a set of equations for m different values of k as , 1,2,..., y j a a x j j m   

   

0 1 (24) where for rolling oscillations RTO-MP-AVT-189 PAPER NBR - 9 NATO UNCLASSIFIED + SWE NATO UNCLASSIFIED + SWE Estimation of Unsteady Aerodynamic Models from Dynamic Wind Tunnel Data , x C y C   a a p  ( ) sin a C a a C      

 

  0 1 a a

p  (25) a    1 1 and similarly for the pitching and yawing oscillations. In the first step a linear regression is used in estimation of parameters and in Eq. (25) from measured in-phase and out-of-phase components at m different a a 0 1 values of k , m > 2.

The second step of regression follows from equations (21a) and (22) replacing  by its estimated value. The resulting equations are , 1,2,..., y j d d x j j m   

   

1 0 1 1 , 1,2,..., y j c d x j j m   

   

2 0 1 2 (26) where for rolling oscillations these terms are

  , j y j C y j C j  

     

1 2 a a p 

        sin , sin x j f j x j f j      

1 1 2 0 sin , sin , C d a c C d       

   

0 1 0 a a (27) p  Similar expressions are obtained for pitching and yawing oscillations. More discussion on the development of regression equations and estimator properties can be found in Ref. [12].

Two-Step Regression was applied to roll oscillatory data from NASA Langley 14x22 Wind Tunnel experiment using an 18% scale F-16XL aircraft. The results of harmonic analysis are plotted in Fig. 5. Both in-phase and out-of-phase components indicate frequency dependence for 25 ° ≤  ≤ 45 °. The coefficient of determination confirms linearity of the aerodynamic model over the set of  considered with the exception of  = 45°. The linear dependence of both components is shown in Fig. 6. The slope of these data is equal to the time constant  . All parameter estimates for  = 34° are presented in Table I. Variation of these parameters 1 0 with  is given in Fig. 7, and a comparison of measured and predicted rolling moment coefficient is shown in Fig. 8.

Table I. Two-Step Regression Parameter estimates for unsteady model at  = 34° and  = 10°.

0 A Step 1 Step 2  a C C lp   l  12.8 0.331 - 0.09 9 - 0.039 ( 0.22 ) ( 0.0041 ) ( 0.0074 ) ( 0.0084 ) 3.4 Equation Error Method The Equation Error (EE) Method is in principle a linear regression. In general, it can be developed either in PAPER NBR - 10 RTO-MP-AVT-189 NATO UNCLASSIFIED + SWE NATO UNCLASSIFIED + SWE Estimation of Unsteady Aerodynamic Models from Dynamic Wind Tunnel Data the time domain or the frequency domain. However, considering the state-space representation in Eq. (4), the time domain approach is not possible because  and  are not measureable variables.

For the frequency domain approach the model is given by Eq. (5) after expressing s as i  . In this case the model has the form A C Bi      ( ) ( ) C     N b i   (28a) ( ) ( ) ( ), 1,2,..., z j C j v j j m    N (28b) where ( ) C  and ( )   are the Fourier transforms of ( ) C t and ( ) t  , ( ) v j is the measurement noise, m is the N N number of frequencies at which the transformed input output data are known, and  is the angular frequency.

In the Equation Error formulation unknown parameters are estimated from minimization of m ( ) ( )( ) ( ) ( ) , 1,2,..., J C j b i A C Bi j j m            1 EE N j j j



1 j  (29) In this formulation ( ) v j are residuals which encompass measurement noise and equation errors. Parameters A , B , and C are related to the aerodynamic coefficients by Eq. (6).

An example of this method, Ref. [6], was applied to measurements obtained from pitch oscillation tests of a 2.5% scale F-16XL model in a water tunnel. Also an alternative to the conventional single-frequency forced- oscillation test input was introduced in Ref. [6] by using Schroeder sweeps. This approach was found to provide substantially more efficient and more effective input for identification of unsteady models and was validated in Ref. [7]. Practical identification using a single-frequency input would require 6 runs, each at a different frequency, to ensure the four parameter model is adequately identified. For this example, a Schroeder sweep is used to allow model identification in one forced-oscillation run at  = 42.5°. Figure 9 shows the  input, ( ) t  , and output, ( ) C t , displacements from their starting values at  = 42.5°. The power spectrum  N of the Schroeder sweep is shown in Fig. 10. Schroeder sweeps provide a flat power spectrum over a specified frequency band with low peak-to-peak amplitudes. Amplitude of the Schroeder sweep is controlled by properly phasing each harmonic. Harmonic analysis provided the in-phase and out-of-phase components shown in Fig. 11. In this case, both components present strong frequency dependence for 30° ≤  ≤ 70°, indicating an unsteady model is needed for that angle of attack region. Final estimates from the EE approach are shown with the Output Error (OE) results in the next section. In this study the EE estimates were used as initial estimates or starting values for the optimization process inherent in the OE method.

3.5 Output Error Method (Frequency Domain) The output error method in the frequency domain can be applied directly using Eq. (28). In this case parameter estimates are obtained by minimization of the mean square output error.

RTO-MP-AVT-189 PAPER NBR - 11 NATO UNCLASSIFIED + SWE NATO UNCLASSIFIED + SWE Estimation of Unsteady Aerodynamic Models from Dynamic Wind Tunnel Data m * ( ) ( ) ( ) J v j v j    OE 1 j  (30) where ( ) ( ) A j C Bi j      ( ) ( ) ( ) v j z j j    ( ) b i j   (31) The OE approach was applied to the same wide-band forced-oscillation data as the EE method discussed previously. Figure 12 shows the parameter estimates and their standard errors over the range of angle of attack where unsteady behaviours were observed. These estimates are consistent with the values obtained from the EE method. Estimate from the two methods in the frequency domain (f.d.) are compared in Table II for  =  42.5°.

Table II. Parameter estimates for unsteady model at  = 42.5°.

Equation Error (f.d.) Output Error (f.d.) Output Error (t.d.)

ˆ ˆ ˆ ˆ ˆ ˆ  ( ) s  ( ) s  ( ) s      A 0.939 0.063 0.895 0.082 B 2.571 0.055 2.521 0.064 n ot applicable C - 0.090 0.051 - 0.101 0.019 b 0.138 0.026 0.144 0.013 C - 0.65 0.39 - 0.70 0.15 - 0.35 0.025 N  C 2. 79 0.19 2.66 0.24 2.60 0.045 Nq a - 3.09 0.39 - 3.09 0.16 - 2.80 0.025  21.53 4.11 20.57 1.90 17.48 0.29  Both the parameters directly estimated for the transfer function given by Eq. (28a) and those indirectly computed using Eq. (6) are shown in Table II. The two methods show very good agreement for the parameter estimates and the output error method obtained slightly better standard errors. Validation of these estimates is demonstrated by application to test data not used for estimation. Fig. 13 shows the output error estimates applied to ramp-and-hold data for a model validation test. The time histories show a very good match.

3.6 Output Error Method (Time Domain) The estimation methods presented so far have been applied to linear aerodynamic models. These methods can be used to estimate all damping and cross derivative terms, however, an OE estimation method is needed when more general nonlinear aerodynamic model structures are required. In Ref. [14], a general mathematical model structure and modelling methodology was presented to address the nonlinear unsteady case. The methodology suggests using the least complex model that provides an adequate representation of the PAPER NBR - 12 RTO-MP-AVT-189 NATO UNCLASSIFIED + SWE NATO UNCLASSIFIED + SWE Estimation of Unsteady Aerodynamic Models from Dynamic Wind Tunnel Data aerodynamic response. With this approach, a progression is made toward more complex models and correspondingly more complex experiment design. The simplifying assumptions associated with the least complex model are removed in progression, only as required to achieve model adequacy. In general, this approach should lead to the best predictive model.

An example from Ref. [14] using measurements from a single-axis forced-oscillation experiment in roll with a 2.5% scale F-16XL model is presented here to demonstrate the method. An additional example using a transport configuration during single-axis yaw oscillations can be found in Ref. [13]. For the F-16XL example, measurements were obtained at six non-dimensional frequencies, k = [0.066, 0.095, 0.131, 0.160, 0.197, 0.262], for angle of attack of 37.5°, during large amplitude (  =30°) roll oscillation tests. As shown in A previous examples harmonic analysis provides key diagnostic information on the adequacy of a linear aerodynamic model. In this example some insight is also provided on the degree of nonlinear behaviour. In Fig. 14 the top plot shows the harmonic analysis for this case over a large range of angle of attack. Frequency dependent behaviour is observed for 25 ° ≤  ≤ 40°, indicating an unsteady model is needed for that angle of attack region. The middle plot presents R for a first harmonic model (the linear case). This plot indicates that for the lower frequencies the linear model does not completely explain the total variation present in the data.

2 st nd rd 2 The lower plot presents R for 1 , 2 , and 3 order harmonic models. A dramatic improvement in R occurs rd 2 st only when a 3 order harmonic model is used. For the k = 0.066 case, R changes from 0.78, with the 1 rd harmonic model, to 0.97 when a 3 order harmonic model is used. One contribution to the nonlinear behaviour is from the static terms in the model. Fig. 15 shows the fairly severe nonlinear character of the roll moment over angle of attack and sideslip. The surface is relatively smooth and linear at low and very high angles of attack; however in the mid angle-of- attack range, 20° ≤  ≤ 50°, fairly severe static nonlinearities occur.

As discussed in Ref. [14], for a single-degree-of-freedom roll test it is reasonable to assume that C can be l expressed as ( , ) C p  and that a slightly more general form than Eq. (4) can be expressed as l ( ) ( ) ( ) ( ) t b t a         l ( ) ( ; ) ( ; ) ( ) ( ) C t C C p t t         a a a p V (32) A broad class of nonlinear unsteady responses can be modelled with this basic model structure where the static term, steady-flow damping term, and indicial term parameters can vary nonlinearly with  . In this example, it is assumed that each of the four unknown parameters ( ( ; ), ( ; ), ( ), ( ) C C a b       ) can be l lp represented by polynomials in  . In general, experiments can be designed to identify models for the static, damping, and indicial terms separately, Ref. [14], although few facilities are capable of executing the motions required. For the OE formulation in the time domain, unknown parameters are estimated through minimization of N ( ) ( ) J v i    OE 1 i  l ( ) ( ) ( ; ( )) ( ; ( )) ( ) ( ) v i C i C i C i p i i          a a a where p (33) V and  i  is computed from the state equation in Eq. (32). In this example only conventional single-frequency, RTO-MP-AVT-189 PAPER NBR - 13 NATO UNCLASSIFIED + SWE NATO UNCLASSIFIED + SWE Estimation of Unsteady Aerodynamic Models from Dynamic Wind Tunnel Data single-axis, forced-oscillation test data is used. Figure 16 shows sample time histories of the large amplitude oscillatory data at  = 37.5° and  = 0°. The figure shows four cycles of oscillation for four frequencies, k = 0 0 [0.066, 0.095, 0.131, 0.160]. Each frequency of the input and output measurements are stacked to ensure all frequency content is included in the analysis. Integration of Eq. (30) must be performed respecting the different conditions at each frequency. The first cycle of predicted response is eliminated to remove any initial transients and ensure that steady harmonic motion is used in calculating residuals. The amplitude,  , is the A same for each frequency; however, the resulting  and  oscillations change during roll oscillations according to the kinematic relationships: 1  sin sin( )sin( )cos( ) cos( )sin( )        

 

0 0 0 sin( )sin( )     1  tan cos( ) tan( )        cos( )cos( )   (34) 0 0   where  and  are the nominal pitch and yaw offsets, respectively. The distortion of the rolling moment  0 response from a steady harmonic wave, in particular for the low frequency cases, indicates nonlinear behaviour.

Under the constraint of planar harmonic data, an iterative 2-stage process shown in Fig. 17 can be used to identify the complete model. Before starting the process, initial parameter values for a linear model are estimated using one of the methods described previously. Small amplitude motions tend to be well represented by linear models. Polynomials of the static values should be well modelled since these values can be directly measured in the wind tunnel; however for this example, the static terms were estimated from the dynamic measurements.

To start the 2-stage process, first model structure determination of the unsteady term polynomials, a(  ) and b(  ) , is performed using stepwise regression (SR) described in Ref. [5]. SR is an extension of linear regression that includes identification of statistically significant mathematical model structures. To setup the regression equation, a variable ( ) y t is formed by subtracting the static term and approximated steady-flow damping term from the force or moment measurement, C (t) , as a ( ) ( ) ( ; ) ( ; ) ( ) y t C t C C p t        a a a p V (35) Then the regression equation, showing the discrete measured values at time t(i) , is given as ( ) [ ( ( )) ( ) ( ( )) ( )] ( ), 1,2,..., y i b i y i a i i i i N          1 E E E E E  (36) where index E indicates the measured values and  (i) is an equation error. Derivatives of measurement, y ,  E are formed using a locally smoothed numerical derivative, Ref. [5]. Parameter estimates are then updated with output error estimation. State and measurement equations for OE estimation are PAPER NBR - 14 RTO-MP-AVT-189 NATO UNCLASSIFIED + SWE NATO UNCLASSIFIED + SWE Estimation of Unsteady Aerodynamic Models from Dynamic Wind Tunnel Data ( ) ( ) ( ) ( ) ; ( 0) (0) t b t a t             1 E ( ) ( ) ( ; ( )) ( ; ( )) ( ) ( ), 1,2,..., y t C i C i C i p i i i N           E a a E a E E E p V (37) A second stage identification is done since ( ; ) C   is only known approximately and the estimates of a(  ) a p and b(  ) are based on this a priori information. To update the damping term, a second variable is formed ( ) ( ) ( ; ) ( ) z t C t C t       a a ( ; ) ( ) C p t    a p V (38) If the measured values and new estimates for the unsteady term are substituted into Eq. (38), a new regression equation is obtained as ( ) ( ; ( )) ( ) ( ), 1,2,..., z i C i p i i i N       E a E E z p V (39) The model structure and parameters for the steady-flow damping term can be Table III . Nonlinear model for C at  = 37.5°,  = 0°.

0 0 l estimated using Eq. (39). This Parameter Nonlinear Model completes the first iteration of the two   stage identification. In the next ( ; ) C   0.007 +0.317       l   iteration, a(  ) and b(  ) can be     2 4 estimated again using the new values ( ; ) C    0.786  1.791  +97.616  lp and model structure for the steady-flow (0.0045) (0.34) (4.48)    damping term. This process is ( ) a  0.527     continued until parameter estimates   have converged. 2 ( ) b  0.624 -8.109  (0.011) (0.14) Since harmonic analysis inferred that at ( 0)    3.63 least a cubic nonlinearity is present and (0.027) the static component is a relatively R 0.97 large component in this case, a cubic representation of the static data was th initially implemented. It was found, however, that a 5 order polynomial was required for an adequate static model. Table III shows the final estimated model from the 2-stage process applied at  = 37.5° and the corresponding overall R . The parameter covariance matrix for this case contained some large pair-wise correlations indicating some identifiability issues. Figure 18 shows the measured and predicted rolling moments for the four frequencies considered. These graphs all show distortion from the regular elliptic shape associated with linear response. For each frequency, the nonlinear model prediction matches the response well. At the lowest frequency, k = 0.066, the distortion or nonlinearity is greatest, but the model match is still very good. Although the fit is very good the model in Table III is not completely satisfactory due to some inconsistency with small-amplitude linear analysis. Ideally the nonlinear model should represent the linear case as a subset then, for example, at  = 0, the nonlinear model estimate, ( ; 0) C    , should be close to the lp RTO-MP-AVT-189 PAPER NBR - 15 NATO UNCLASSIFIED + SWE NATO UNCLASSIFIED + SWE Estimation of Unsteady Aerodynamic Models from Dynamic Wind Tunnel Data corresponding linear model analysis of the small amplitude data, a value of -0.13. The differences may reflect a number of possible sources of error besides the model structure for the damping and unsteady terms, such as data accuracy and information content of the data.

4.0 CONCLUDING REMARKS

This paper presented several methods for estimating mathematical models useful for stability and control analysis that can be applied to both CFD simulations and wind tunnel measurements. These methods can accommodate general aerodynamic behaviors and characterize model parameter uncertainty. Uncertainty information significantly improves model validation tests by providing more information than conventional time history comparisons. Conventional stability and control derivatives are retained in the recommended model structures, even for the nonlinear case, to take advantage of the experience and significant knowledge base built up through years of aerospace engineering practice.

Collaboration between computational and experimental researchers offers a significant opportunity to advance the technology in both disciplines. A number of modelling research areas present opportunity for collaboration that can be further developed. Dynamic test techniques in wind tunnels are limited to current facility capabilities and most facilities only allow planar harmonic motion. Using CFD simulations offers an opportunity to further the design of both dynamic inputs for experiments and desirable facility capabilities.

Efficacy of simple improvements such as using wide-band inputs, oscillatory coning, and plunging dynamics can be relatively quickly and inexpensively assessed with CFD experiments. Some flight conditions are not easily reproduced in a tunnel and some conditions are either too transient or too dangerous for flight test such as with the study of transport loss of control problems. Collaboration between CFD and SID can facilitate both discovery and analysis of these difficult-to-model flight conditions. This approach can lower concerns for pilot and vehicle safety and potentially reduce costs.

REFERENCES

[1] Morrison, Joseph H., and Hemsch, Michael J ., “Statistical Analysis of CFD Solutions from the Third AIAA Drag Prediction Workshop,” AIAA 2007 -254, January 2007.

[2] Hemsch, Michael J., “Statistical Analysis of Computational Fluid Dynamics Solutions from the Drag Prediction Workshop,” Journal of Aircraft , Vol. 41, No. 1, January-February, 2004.

[3] Vicroy, D. and Loeser, T., “SACCON Dynamic Wind Tunnel Tests at DNW - NWB and 14´x22´ NASA LaRC”, AIAA-2010-4394, June 2010.

[4] Murphy, Patrick C, Klein, Vladislav, Frink, Ne al T., and Vicroy, Daniel D, “System Identification Applied to Dynamic CFD Simulation and Wind Tunnel Data”, AIAA 2011 -xxxx, August 2011.

st [5] Klein, Vladislav and Morelli, Eugene, “Aircraft System Identification: Theory and Practice,” 1 edition, AIAA Inc., Reston, Virginia, 2006.

[6] Murphy, P.C. and Klein, V., “Estimation of Aircraft Unsteady Aerodynamic Parameters from Dynamic Wind Tunnel Testing,” AIAA 2001 -4016, August 2001.

[7] Murphy, P.C. and Klein, V., “Validation of Methodology for Estimating Aircraft Unst eady Aerodynamic PAPER NBR - 16 RTO-MP-AVT-189 NATO UNCLASSIFIED + SWE NATO UNCLASSIFIED + SWE Estimation of Unsteady Aerodynamic Models from Dynamic Wind Tunnel Data Parameters from Dynamic Wind Tunnel Tests,” AIAA 2003 -5397, August 2003.

[8] Jones, Robert T. and Fehlner, Leo F., “Transient Effects of the Wing Wake on the Horizontal Tail,” NACA TN-771, 1940.

[9] Tobak, Murray, “On the Use of Indicial Function Concept in the Analysis of Unsteady Motions of Wings and Wing-Tail Combinations, NACA Rep. 1188, 1954.

[10] Klein, Vladislav and Norderer, Keith D., “Modeling of Aircraft Unsteady Aerodynamic Characteristics.

Part I – Postulated Models,” NASA TM 109120, May, 19 94.

[11] Klein, Vladislav, “Modeling of Longitudinal Unsteady Aerodynamics of a Wing - Tail Combination,” NASA CR-1999-209547, September, 1999.

[12] Klein, Vladislav, Murphy, P.C., and Szyba, Nathan M., “Analysis of Wind Tunnel Lateral Oscillatory Data of the F-16XL Aircraft,” NASA/TM -2004-213246, 2004.

[13] Murphy, P.C. and Klein, V., “Transport Aircraft System Identification Using Roll and Yaw Oscillatory Wind Tunnel Data,” AIAA 2010 -6122, August 2010.

[14] Murphy, P.C. and Klein, V., “Progressive Aerodynamic Model Identific ation From Dynamic Water Tunnel Test of the F- 16XL Aircraft,” AIAA 2004 -5277, August 2004.

RTO-MP-AVT-189 PAPER NBR - 17 NATO UNCLASSIFIED + SWE NATO UNCLASSIFIED + SWE Estimation of Unsteady Aerodynamic Models from Dynamic Wind Tunnel Data

FIGURES

S = 5.90 ft , c = 0.915 ft, b = 6.85 ft Figure 1b. Three - view drawing of 2.5% F - 16XL water - Figure 1a. Model geometry for NASA tunnel model .

GTM, experimental sub - scale aircraft.

f (Hz) C N q measured 0.05 computed 0.08 -10 0 10 20 30 40 50 C N q 0.12 0.23 C N  0.43 0.54 -10 0 10 20 30 40 50 0.86 0.95 R 0.9 0.85 0 0.005 0.01 0.015 0.02 0.025 0.03 -10 0 10 20 30 40 50 k  , deg Figure 3. Results of nonlinear regression for measured Figure 2. Harmonic analysis for normal - force coefficient, out - of - phase components of NASA GTM,   °, NASA GTM configuration,   °.     °.

 PAPER NBR - 18 RTO-MP-AVT-189 NATO UNCLASSIFIED + SWE NATO UNCLASSIFIED + SWE Estimation of Unsteady Aerodynamic Models from Dynamic Wind Tunnel Data C l p -1 f (Hz) 0.3 -2 0 20 40 60 80 C (t) 0.5 N 0.2 0.7 f = 0.86 Hz C 0.78 l  0.9 C static N -0.2 0.95 0 20 40 60 80 C measured N 1.1 1 2 C predicted N R -10 0 10 20 30 40 50 0.5  , deg 0 20 40 60 80  , deg   Figure 4. Comparison of measured and predicted Figure 5. Har monic analysis for roll - moment coefficient, normal - force coefficient, NASA GTM, f = 0.86 Hz, 18% F - 16XL configuration,     °.

    °.

 0.2 -0.2 measured C 0 l p estimated -0.4 -0.2 C l p -0.6 -0.4 0 10 20 30 40 50 60 70 80 -0.8 -1 C l  -0.1 -1.2 -0.2 -1.4 0.6 -0.2 -0.18 -0.16 -0.14 -0.12 -0.1 -0.3 0 10 20 30 40 50 60 70 80 0.4 C l   , deg   0.2 Figure 6. Measured vs Estimated roll - Figure 7 a . Two - Step Regression estimates, 18% F - 16XL ,   °.

moment c oeff icients, 18% F - 16XL ,   °.   -0.2 0 10 20 30 40 50 60 70 80 0.6 0.4  a 1 0.2 -0.2 -5 0 10 20 30 40 50 60 70 80 0 10 20 30 40 50 60 70 80  , deg  , deg   , deg    Figure 7 c . Two - Step Regression estimates, 18% F - 16XL , Figure 7 b . Two - Step Regression estimates, 18% F - 16XL ,   °.

  °.

  -5 0 10 20 30 40 50 60 70 80  , deg RTO-MP-AVT-189 PAPER NBR - 19 NATO UNCLASSIFIED + SWE NATO UNCLASSIFIED + SWE Estimation of Unsteady Aerodynamic Models from Dynamic Wind Tunnel Data 0.05 0.05 C (t) C (t) l l 0 0 f = 0.3 Hz,  = 20 f = 1.1 Hz,  = 20 0 0 -0.05 -0.05 -10 -5 0 5 10 -10 -5 0 5 10  , deg  , deg measured 0.05 0.05 predicted C (t) C (t) l l 0 0 f = 0.3 Hz,  = 30 f = 1.1 Hz,  = 30 0 0 -0.05 -0.05 -10 -5 0 5 10 -10 -5 0 5 10  , deg  , deg 0.05 0.05 C (t) C (t) l l 0 0 f = 0.3 Hz,  = 40 f = 1.1 Hz,  = 40 0 0 -0.05 -0.05 -10 -5 0 5 10 -10 -5 0 5 10  , deg  , deg Figure 8 . Two - Step Regr ession model of roll - moment predicted and mea sured data , 18% F - 16XL ,   °.

 10 40 ( ) t  -5 0 50 100 150 200 250 300 350 0.4 -20 0.2 PSD Magnitude (dB) ( ) C t -40 N -0.2 -60 0 0.2 0.4 0.6 0.8 1 0 50 100 150 200 250 300 350 time (sec) frequency, Hz Figure 9 . Measurements of  perturbations from Figur e 10 . Harmonic content of  wide - band  = 42.5° and normal force during wide - band input, 2.5% F - 16XL in water tunnel,   °.

 experiment with 2.5% F - 16XL in water tunnel .

PAPER NBR - 20 RTO-MP-AVT-189 NATO UNCLASSIFIED + SWE NATO UNCLASSIFIED + SWE Estimation of Unsteady Aerodynamic Models from Dynamic Wind Tunnel Data 4 5 C N  C (  ) Nq k=0.05 0 10 20 30 40 50 60 70 k=0.1 k=0.15 k=0.2 C N q k=0.25 -1 -2 -3 30 40 50 60 70 0 10 20 30 40 50 60 70 (deg)  Figure 11 . In - phase and out - of - phase coefficients from wide - band experiments, 2.5% F - 16XL,   °.

 C (  ) N  a -1 -1 -2 -2 -3 -3 30 40 50 60 70  (deg) -4 Figure 12 a . Normal Force model parameters and 2 -  bounds , using OE , 2.5% F - 16XL,   °.

 -5 30 40 50 60 70  (t) (deg) q(t) (deg/sec) -2 0 10 20 30 40 50 0.1  (t)  1 0.05 0 10 20 30 40 50 predicted 1.8 C (t) measured N 1.6 -10 0 10 20 30 40 50 30 40 50 60 70 time (sec)  (deg) Figure 12 b . Normal f orce model parameters Figure 13. Model validation of C response to N and 2 -  bounds , using OE , 2.5% F - 16XL, ramp (  =40° - 45 °), non - dim pitch rate of 0.03,   °. 2.5% F - 16XL.

 RTO-MP-AVT-189 PAPER NBR - 21 NATO UNCLASSIFIED + SWE NATO UNCLASSIFIED + SWE Estimation of Unsteady Aerodynamic Models from Dynamic Wind Tunnel Data 0.06 k=0.066 0.04 k=0.095 k=0.131 0.02 k=0.16 C l k=0.197 k p A C k=0.262  0 l -1 C / l B -0.02 -2  =30  ,  =0  A 0 -0.04 -3 0 10 20 30 40 50 60 70 80 90 -0.06 1 -20  , deg k=0.066 -10 k=0.095 0.8 k=0.131 k=0.16 R 40 k=0.197 10 0.6 ) 20   deg l  , deg k=0.262 C   deg  , deg ( 20 0 R 0.4 0.2 Figure 15. Surface plot of roll moment static  =30  ,  =0  A measurements , 2.5% F - 16X L . 0 0 10 20 30 40 50 60 70 80 90  , deg Order=1 Order=2 0.8 Order=3 R ) 0.6 l   t  C ( 0.4 (t), deg R   =30  ,  =0  , k=0.066 A -10 0.2 -20 0 10 20 30 40 50 60 70 80 90  , deg   deg)  0 p (t) Figure 14. Harmonic A nalysi s of 2.5% F - 16XL p(t), deg/s -10 during large amplitude roll oscillat i ons ,  = 30° .

A -20 0.1 0.05 ( ) C t l (t) 0 l C -0.05 -0.1 0 50 100 150 200 250 300 350 400 450 time, s t ime (sec) Figure 16. Time histories of the large amplitude oscillatory data at  = 37.5°,  = 0°, and four f requencies, k = [0.066, 0 0 0.095, 0.131, 0.160] , 2.5% F - 16XL .

PAPER NBR - 22 RTO-MP-AVT-189 NATO UNCLASSIFIED + SWE NATO UNCLASSIFIED + SWE Estimation of Unsteady Aerodynamic Models from Dynamic Wind Tunnel Data two stage optimization two stage optimization yes yes

  y ˆ ˆ   y ˆ ˆ

   b    b    b    b E z E z i no i no E E SR ML Converged? SR ML Converged? OE eqn A eqn A SR SR    a ˆ    a ˆ    a ˆ    a ˆ ˆ ˆ C C a a p p

  y   y

E E 1  i 1  i eqn B eqn B   eqn A: ; , p z t C t    

     

0 a E   p V     eqn B: ; , ; , p y t C t C C t         

         

0 0 E   a a a V p   Figure 1 7. Block diagram of model identification procedure using stepwise regression (SR) and an output error (OE ) estimation .

0.1 0.1 0.05 0.05 C C l l 0 0 -0.05 -0.05 k =0.066 k = 0.095 -0.1 -0.1 -20 -10 0 10 20 -20 -10 0 10 20  , deg  , deg 0.1 0.1 meas.

meas.

comp.

pred.

0.05 0.05 C C l l 0 0 -0.05 -0.05 k = 0.160 k = 0.131 -0.1 -0.1 -20 -10 0 10 20 -20 -10 0 10 20  , deg  , deg Figure 18. Comparison of measured and computed rolling - moment coefficient from time domain OE method at  = 37.5°.

RTO-MP-AVT-189 PAPER NBR - 23 NATO UNCLASSIFIED + SWE