Nonlinear Time-Domain System Identification of Aircraft Stability and Control Derivatives

Nonlinear Time-Domain System Identification of Aircraft

Stability and Control Derivatives

∗

Nhan Nguyen

NASA Ames Research Center, Moffett Field, CA 94035

†

Juntao Xiong

KBR Wyle, Moffett Field, CA 94035

This paper presents a nonlinear time-domain system identification method to estimate the stability derivatives for NASA-Boeing Transonic Truss-Braced Wing configuration.

Nonlinear aerodynamics in large-amplitude forced oscillations can cause significant fre- quency interactions with the linear response at the input frequency. This is called a non- linear spillover effect. In our previous work, nonlinear frequency-domain methods have been developed for estimation of stability and control derivatives. In this work, we are proposing to develop a nonlinear time-domain system identification method using a recur- sive least-squares technique. The proposed method is applied to compute the linear and nonlinear stability derivatives of the Transonic Truss-Braced Wing.

I. Introduction

Stability and control are an important aircraft design requirement. Vehicle stability requires that an aircraft design be statically and dynamically stable and controllable in all three axes of the motion in roll, pitch, and yaw. Vehicle control is accomplished by aerodynamic control surfaces or other control effectors.

Aircraft are designed to have sufficient stability and control margin in roll, pitch, and yaw axes. Static stability of an aircraft requires the information on the stability derivatives C , C , and C . Similarly, l m n β α β the roll, pitch, and yaw damping derivatives C , C , and C contribute to the dynamic stability of an l m n p q r aircraft. In extreme maneuvers or turbulence encounters, the excursion in the angle of attack well outside the nominal range could adversely affect stability and control characteristics. Flight control is designed to provide stability augmentation and pilot command tracking in nominal flight operations. When the aircraft stability is degraded in high-alpha deep stall conditions, flight control may not be able to provide sufficient stability augmentation to maintain good handling qualities. This could result in flight safety issues.

Numerical prediction of dynamic stability and control derivatives is a current research topic. Recent work 1–3 in this area can be found in many references. Many of these technique do not adequately capture the various terms of the dynamic stability derivatives. More importantly, the effects of the frequency response on the dynamic stability derivatives is not well addressed especially for aircraft operating in transonic flight regime. Murman uses reduced-frequency approach to address the frequency response, but the approach lacks sufficient details to enable a systematic approach to determining the dynamic stability derivatives. A 4–7 frequency-domain dynamic stability estimation has recently been developed and applied to the Mach 0.8 Transonic Truss-Braced Wing (TTBW) aircraft, shown in Figure 1. Using unsteady RANS CFD simulation results, the time histories of the aerodynamic coefficients are transformed into the frequency domain by a Fourier series analysis as a function of the reduced frequency. A frequency-domain transfer function is designed to capture both the steady state and dynamic derivatives. The transfer function is formulated based on the Theodorsen’s theory. A frequency-domain regression analysis is then performed to estimate the transfer function. The proposed method can estimate those dynamic stability derivatives that are usually not captured in many analyses but do exist in the Theodorsen’s theory such as C due to the apparent m ˙ q mass effect.

∗ Technical Group Lead and Senior Research Scientist, nhan.t.nguyen@nasa.gov † Research Engineer, juntao.xiong@nasa.gov 1 of 16 American Institute of Aeronautics and Astronautics Figure 1. Boeing Mach 0.8 Transonic Truss-Braced Wing Aircraft Concept This study continues the investigation of the stability and control estimation using high-fidelity CFD. In this study, we develop a time-domain estimation of the stability derivatives for the TTBW taking into account significant nonlinearity in the aerodynamic response of an aircraft to large-amplitude forced oscillations. In our prior work, we initially developed a nonlinear frequency-domain system identification method for stability th 5 control derivatives using a 6 -degree polynomial. In that study, the prescribed polynomial is proved to be sufficient for control derivative estimation. In another study, we extended the nonlinear frequency-domain system identification method to handle arbitrary order polynomial models in order to accurately capture the nonlinear aerodynamics. In that study, the input is restricted to a sinusoidal input.

System identification is a process of estimation of the system dynamic response to a prescribed input.

Sinusoidal inputs are a preferred method as it provides the greatest excitation energy at a given input frequency. Swept sine method is used to speed up the system identification process. Other methods involve using multi-sine inputs. In our previous study, we explore the use of a square wave input which in theory can excite all the odd harmonics of the input frequency. While the computational process can be reduced, the excitation energy is spread over all the frequencies in an inversely proportional relationship with the frequency, thereby resulting in the lack of sufficient excitation energy at high frequencies.

In this study, we propose a nonlinear time-domain system identification method that can be broadly applicable to arbitrary polynomial degrees and arbitrary inputs such as the square-wave input which causes difficulty in the frequency-domain system identification method.

II. Nonlinear Frequency-Domain System Identification of Stability Derivatives

Consider the unsteady lift of an aircraft under a harmonic motion about the pitch axis described by θ = θ sin ωt = θ sin kτ (1) 0 0 where θ ( τ ) is the pitch angle as a function of the non-dimensional time or distance traveled in semi-chord 2 U t ω ¯ c τ = and k = is the reduced frequency.

¯ c 2 U In our previous work, the frequency response of the unsteady lift in nonlinear flow is expressed as the n -th power function of the pitch angle N N X X n n n ˜ C = ( Q + iS ) θ = ( Q + iS ) θ sin kτ (2) L n n n n n =1 ... n =1 ...

2 of 16 American Institute of Aeronautics and Astronautics The complex number representation captures the magnitude and phase of the frequency response of the n ˜ partial derivative of C with respect to θ ( τ ) . That is, L ˜ ∂ C L ˜ = C = Q + iS (3) L n n n θ ∂θ n Since the pitch motion is sinusoidal, the n -th power of the sine function contributes to the frequency response at an n -th order harmonic as well as all the lower odd or even harmonics. This implies that the odd power of the sine function contributes to the frequency response at the fundamental frequency of the pitch motion. If one were to extract the fundamental frequency response of the unsteady lift and use it to estimate the stability derivatives, it would have resulted in errors. The fundamental frequency response not only contains the first-order sinusoidal linear response but the contributions from all the other odd-power sinusoidal nonlinear responses. This is called a “spillover” effect of nonlinear response which contaminates the linear response.

This frequency interaction can be better understood by a further analysis of the frequency response of the power of the sine function.

Using the Euler’s formula, the sine function can be cast as ikτ − ikτ e − e sin kτ = (4) 2 i Consider the case when n is even. Then, n − 1 n 2 m − X 2 1 n ! ( − 1) n !

n sin kτ = + cos ( n − 2 m ) kτ (5)
n 2 n − 1 n 2 2 m ! ( n − m )!

!

m =0 This implies that the frequency response of an even-power of the sine function is a cosine series at even harmonics from 2 to n . The even harmonics do not interact with the linear frequency response at the fundamental frequency.

Consider the case when n is odd. Then, n − 1 n − 1 m − X 2 ( − 1) n !

n sin kτ = sin ( n − 2 m ) kτ (6) n − 1 2 m ! ( n − m )!

m =0 The frequency response of an odd-power of the sine function is a sine series at odd harmonics from 1 to n . The odd power of the sine function contributes to the fundamental frequency response. Therefore, if the fundamental frequency response is used for the linear stability derivative estimation, the frequency response would include the nonlinear contributions which would result in errors.

With this in mind, we proceed with the unsteady lift analysis. The unsteady lift is expressed as n − 1 n − 1 N N 2 2 X X X X n n ˜ C = ( Q + iS ) θ a sin ( n − 2 m ) kτ + ( Q + iS ) θ b cos ( n − 2 m ) kτ (7) L n n mn n n mn 0 0 n =1 , 3 ,... m =0 n =2 , 4 ,... m =0 where n − 1 m − ( − 1) n !

a = (8) mn n − 1 2 m ! ( n − m )!

n m − ( − 1) n !

b = (9) mn n − 1 2 m ! ( n − m )!

Note that the bias term of the even-power of the sine function is removed since we are only interested in the alternating signal.

3 of 16 American Institute of Aeronautics and Astronautics   s ¯ c iω ¯ c d We define the non-dimensional Laplace variable ¯ s = = = ik = L . Therefore, we obtain 2 U 2 U dτ 1 d i = . Hence, k dτ n − 1 n − 1 N N 2 2 X X X X n n ˜ C = Q θ a sin ( n − 2 m ) kτ + S θ a ( n − 2 m ) cos ( n − 2 m ) kτ L n mn n mn 0 0 n =1 , 3 ,... m =0 n =1 , 3 ,... m =0 n n − 1 − 1 N N 2 2 X X X X n n + Q θ b cos ( n − 2 m ) kτ − S θ b ( n − 2 m ) sin ( n − 2 m ) kτ (10) n mn n mn 0 0 n =2 , 4 ,... m =0 n =2 , 4 ,... m =0 Using the orthogonality property, the Fourier sine and cosine series coefficients for the odd p -th harmonic are evaluated as n − 1 2 Npπ 2 Npπ Z Z N X X k k n ˜ C sin pkτ dτ = Q θ a sin ( n − 2 m ) kτ sin pkτ dτ L n mn 0 0 n =1 , 3 ,... m =0 N X N π p n = Q θ a (11) n pn k n =1 , 3 ,...

n − 1 2 Npπ 2 Npπ Z Z N 2 k k X X n ˜ C cos pkτ dτ = S θ a ( n − 2 m ) cos ( n − 2 m ) kτ cos pkτ dτ L n 0 mn 0 0 n =1 , 3 ,... m =0 N X N π p n = S θ a p (12) n pn k n =1 , 3 ,...

n − p for p = n − 2 m or m = where N is the number of periods and p p − 1 − ( − 1) n !

a = , n, p odd (13) pn n − p n + p n − 1 2 ! !

2 2 For the even p -th harmonic, the Fourier sine and cosine series coefficients are evaluated as n 2 Npπ 2 Npπ − 1 Z Z N X X k k n ˜ C sin pkτ dτ = − S θ b ( n − 2 m ) sin ( n − 2 m ) kτ sin pkτ dτ L n mn 0 0 n =2 , 4 ,... m =0 N X N π p n = − S θ b p (14) n pn k n =2 , 4 ,...

n 2 Npπ 2 Npπ − 1 Z Z X k k n ˜ C cos pkτ dτ = Q θ b cos ( n − 2 m ) kτ cos pkτ dτ L n mn 0 0 m =0 N X N π p n = Q θ b (15) n pn k n =2 , 4 ,...

n − p for n − 2 m = p or m = where p − ( − 1) n !

b = , n, p even (16) pn n − p n + p n − 1 ! !

2 2 4 of 16 American Institute of Aeronautics and Astronautics We express in Fourier sine and cosine series coefficients in a matrix form as ( ) 2 Npπ   Z n k N πθ p ˜ { a } Q = C sin pkτ dτ , n, p odd (17) pn n L k ( ) 2 Npπ   Z n k N πθ p ˜ { a p } S = C cos pkτ dτ , n, p odd (18) pn n L k ( ) 2 Npπ   Z n k N πθ p ˜ { b } Q = C cos pkτ dτ , n, p even (19) pn n L k ( ) 2 Npπ   Z n k N πθ p ˜ − { b p } S = C sin pkτ dτ , n, p even (20) pn n L k The frequency response of the power of the sine function can be obtained by solving for Q and S by n n the matrix inversion of Eqs. (17)-(20). Consider the case when p = 1 which corresponds to the fundamental harmonic, Q and S are the frequency response of the linear stability derivative since 1 1 ˜ C = Q + iS (21) L 1 1 θ The matrix solution gives 2 Npπ 2 Npπ Z Z N X k k k k ˜ ˜ Q = C n sin nkτ dτ = C (sin kτ + 3 sin 3 kτ + 5 sin 5 kτ + · · · ) dτ 1 L L N πθ N πθ p 0 p 0 0 0 n =1 , 3 ,...

(22) 2 npπ 2 Npπ Z Z N X k k k k ˜ ˜ S = C cos nkτ dτ = C (cos kτ + cos 3 kτ + cos 5 kτ + · · · ) dτ (23) 1 L L N πθ N πθ p 0 0 p 0 0 n =1 , 3 ,...

Consider the linear contribution of unsteady lift which is given by   ˜ ˜ C = C θ = H (¯ s ) ( c + c ¯ s ) + d ¯ s + d ¯ s θ (24) L L 0 1 1 2 θ The aerodynamic transfer function H (¯ s ) has the form P N d i a ¯ s i i =1 H (¯ s ) = 1 + (25) P N − 1 d N i d ¯ s + b ¯ s i i =0 As a note for incompressible flow, the aerodynamic transfer function is known as the Theodorsen’s function. Therefore, for any general linear flow, the aerodynamic transfer function can be represented by H (¯ s ) . A system identification can be performed to identify H (¯ s ) and the coefficients c , c , d , and d 0 1 1 2 ˜ ˙ ˙ ¨ which are the partial derivatives of C with respect to θ and θ due to the circulatory lift and θ and θ due to L 4, 7 the so-called non-circulatory lift, respectively. This has been developed in previous studies.

If the flow is quasi-steady with k ≈ 0 , then we write   2 2 dθ dθ d θ dθ d θ  * ˜  C = H (0) c θ + c + d + d = C θ + C + C (26) L 0 , 1 1 , 1 1 , 1 2 , 1 L L L 1 α q ˙ q  2 2 dτ dτ dτ dτ dτ whereupon we let θ = α , hence C = C . Otherwise, we need to estimate H (¯ s ) through a system L L θ α identification procedure . Thus, we see that c = C , c + d = C , and d = C .

0 , 1 L 1 , 1 1 , 1 L 2 , 1 L α q ˙ q Using the non-dimensional Laplace transform, the linear contribution to the unsteady lift coefficient is expressed as
˜ C = C + iC k − C k θ (27) L L L L 1 α q ˙ q The linear contribution to the unsteady lift coefficient expressed in terms of Q and S is 1 1 ˜ C = ( Q + iS ) θ (28) L 1 1 5 of 16 American Institute of Aeronautics and Astronautics Therefore, the linear stability derivatives can be estimated from Q and S by a linear regression which 1 1 minimizes the following cost function over a set of the reduced frequencies k ≈ 0 , i = 1 , . . . , , M > 1 : i M n o X     2 2 J = C − C k − Q ( k ) + C k − S ( k ) (29) L L i 1 i L i 1 i α ˙ q q i =1 The regression yields the equation       P P M M M 0 − k C Q ( k ) L 1 i i α i =1 i =1 P P       M M = (30)  0 k 0   C   k S ( k )  L i 1 i i q i =1 i =1 P P P M M M 2 4 2 − k 0 k C − k Q ( k ) L 1 i i i ˙ q i i =1 i =1 i =1 Therefore, the stability derivatives are obtained as P P P P M M M M 4 2 2 k Q ( k ) − k k Q ( k ) 1 i 1 i i i i i =1 i =1 i =1 i =1 C = (31) L   α P P M M 4 2 M k − k i =1 i i =1 i P M k S ( k ) i 1 i i =1 C = (32) P L q M k i =1 i P P P M M M 2 2 k Q ( k ) − M k Q ( k ) 1 i 1 i i i i =1 i =1 i =1 C = (33) L   ˙ q P P M M 4 2 M k − k i =1 i i =1 i Note that the condition of invertibility of matrix in Eq. (30) requires at least two orthogonal sine signals with distinct frequencies. Otherwise, the stability derivatives C and C are non-unique due to the fact L L α ˙ q d θ that θ and are simply the same signal but with different scaling factors. Thus, for a single frequency dτ input, the stability derivative C cannot be estimated and the stability derivative C in fact includes the L L ˙ q α contribution of the stability derivative C which is generally small and therefore can be neglected.

L ˙ q The unsteady lift contribution by nonlinear aerodynamics is modeled as n ˜ C = ( c + c ¯ s ) θ (34) L 0 ,n 1 ,n n for n > 1 .

The unsteady lift is expressed as n n dθ dθ n n ˜ C = c θ + c = C θ + C (35) L 0 ,n 1 ,n L n L n n α dθ dτ dτ dτ n q ¯ c dθ n − 1 dθ n − 1 Note that = nθ = nθ . Thus, C is the nonlinear stability derivative with respect to L n dτ dτ 2 U dθ dτ n − 1 nθ q .

From the frequency response, the unsteady lift is obtained as n S dθ n n n ˜ C = ( Q + iS ) θ = Q θ + (36) L n n n n k dτ Therefore, the nonlinear stability derivatives can be estimated by a linear regression which minimizes the following cost function:   M h i X 2 J = [ C − Q ( k )] + k C − S ( k ) (37) L n n i i L n n i α dθ dτ i =1 The nonlinear stability derivatives are obtained as P M Q ( k ) n i i =1 C = (38) L n α M 6 of 16 American Institute of Aeronautics and Astronautics P M k S ( k ) i n i i =1 C = (39) L n P dθ M dτ k i =1 i If H (¯ s ) is taken to be unity which implies a quasi-steady aerodynamic assumption, we express the frequency-domain representation of the unsteady lift in the time domain as   N n 2 X dθ d θ n ˜ C = C θ + C + C (40) L L n L n L α dθ ˙ q dτ dτ dτ n =1 dθ n d θ The terms accounts for the expression iS θ and the term is due to the apparent-mass contribution.

n 2 dt dt n n Given an input θ ( t ) , we can compute the various input terms θ , q , and ˙ q . The unsteady lift can also expressed as ! !

N N X X q ¯ c ˙ q ¯ c n − 1 n − 1 ˜ C = C + C θ θ + C + C nθ + C (41) n n L L L L L L α α q dθ ˙ q dτ 2 U 4 U n =2 n =2 In this form, we see that the stability derivatives with respect to the angle of attack α = θ and the pitch rate q are nonlinear functions of the angle of attack.

We cast the unsteady lift in a compact form ⊤ ˜ C = Θ Φ ( t ) (42) L h i ⊤ C C · · · C C C where Θ = L L L L L is a parameter vector of unknown coefficients which α q N N ˙ q α dθ dτ h i ⊤ N 2 dθ N dθ d θ represent the stability derivatives and Φ = is called a regressor vector.

θ · · · θ dτ dτ dτ The estimation of the parameter vector Θ can be obtained by minimizing the cost function N t X min J = ∥ ϵ ( t ) ∥ (43) i ˆ n × m Θ ∈ R i =1 where ⊤ ˜ ˆ ϵ ( t ) = C ( t ) − Θ Φ ( t ) (44) i L i i The sum in the cost function implies that the minimization is performed over the entire time series data all at once, hence the term batch least-squares. The regression is performed after all the system identification data have been acquired.

ˆ ˜ When J is minimized, the approximation error is also minimized. Then, C approximates C in a L L least-squares sense. Thus, the parameter estimation problem is posed as an optimization or minimization problem.

The necessary condition is given by N N t t h i X X ∂J ∂ ϵ ( t ) i ⊤ ⊤ ˆ ˜ = ∇ J = ϵ ( t ) = Φ ( t ) Φ ( t ) Θ − C ( t ) = 0 (45) ˆ i i i L i Θ ˆ ˆ ⊤ ∂ Θ ∂ Θ i =0 i =0 ˆ where ∇ J is called the gradient of J with respect to Θ .

ˆ Θ We obtain the least-squares regression equation N N t t X X ⊤ ˆ ˜ Φ ( t ) Φ ( t ) Θ = Φ C ( t ) (46) i i L i i =0 i =0 P P N N t ⊤ t ˜ Let A = Φ ( t ) Φ ( t ) and b = Φ ( t ) C ( t ) . The least-squares linear regression equation is i i i L i i =0 i =0 ˆ A Θ = b (47) which yields − 1 ˆ Θ = A b (48) 7 of 16 American Institute of Aeronautics and Astronautics assuming A is a non-singular matrix if there are sufficient and unique data. This is also called the persistent excitation condition.

Instead of the summing cost function, the cost function could be defined at only time t i ⊤ J = ϵ ( t ) ϵ ( t ) (49) i i The minimization of this cost function leads to the second-order gradient method
− 1 ∗ 2 ˆ ˆ Θ ( t ) = Θ ( t ) − ∇ J ( t ) ∇ J ( t ) (50) i i ˆ i ˆ i Θ Θ where h i ⊤ ˆ ˜ ∇ J ( t ) = Φ ( t ) Φ ( t ) Θ − C ( t ) (51) ˆ i i i L i Θ The estimation is recursive in that the estimation takes place as the new data becomes available. The second-order gradient method can also be expressed as
− 1 ∗ 2 ⊤ ˆ ˆ Θ ( t ) = Θ ( t ) − ∇ J ( t ) Φ ( t ) ϵ ( t ) (52) i i ˆ i i i Θ The Hessian matrix ∇ J is required to be invertible. It is noted that the inverse of the Hessian matrix ˆ Θ does not always exist and is generally numerically intensive. So, a first-order approximation can be made by 2 ∗ ˆ ˆ recognizing that ∇ J ≈ η ≥ 0 if Θ ( t ) is in the neighborhood of the minimum Θ ( t ) , where η is a small ˆ i i Θ valued positive definite matrix. This leads to the steepest descent method given by   ∗ ˆ ˆ ˆ Θ ( t ) = Θ ( k ) − η ∇ J Θ ( k ) (53) i ˆ Θ ( k ) A more general formulation of the unsteady lift should take into account the aerodynamic transfer function H (¯ s ) . The unsteady lift when accounting for the aerodynamic transfer function H (¯ s ) is represented by a more complicated expression in the time domain given by N d X w ¯ s i H (¯ s ) = 1 + (54) ¯ s − p i i =1 where p are the poles of the unsteady aerodynamic transfer function which are required to be stable such i that p < 0 .

i The unsteady lift is then expressed as     N n 2 X dθ dθ d θ n ˜ C = H (¯ s ) C θ + C + C θ + C + C L L L L n L n L α q α dθ ˙ q dτ dτ dτ dτ n =2   N N d n 2 2 X X dθ d θ a ¯ sθ + b ¯ s θ i i n = C θ + C + C + (55) L n L n L α dθ ˙ q dτ dτ ¯ s − p dτ i n =1 i =1 We define the unsteady aerodynamic state variables ¯ sθ q = (56) i ¯ s − p i ¯ s θ r = (57) i ¯ s − p i Consider a general multi-sine input N f X θ ( τ ) = θ sin k τ (58) 0 ,j j j =1 In the time domain, the state variables q and r are described by i i N f X dq dθ i = p q + = p q + θ k cos k τ (59) i i i i 0 ,j j j dτ dτ j =1 8 of 16 American Institute of Aeronautics and Astronautics N f X dr d θ i = p r + = p r − θ k sin k τ (60) i i i i 0 ,j j j dτ dτ j =1 The steady-state solutions of q and r are obtained as i i N f X θ k 0 j q = ( k sin k τ − p cos k τ ) (61) i j j i j 2 2 k + p j i j =1 N f X θ k j r = ( p sin k τ + k cos k τ ) (62) i i j j j 2 2 k + p j i j =1 The unsteady lift is then given by   N N d n 2 X X dθ d θ n ˜ C = C θ + C + C + ( a q + b r ) (63) n n L L L L i i i i α dθ ˙ q dτ dτ dτ n =1 i =1 h i ⊤ C C · · · C C C a b · · · a b which can be cast in a form of Eq. (42) where Θ = L L L L L 1 1 N N α q N N ˙ q d d α dθ dτ h i ⊤ N 2 dθ N dθ d θ and Φ = .

θ · · · θ q r · · · q r 2 1 1 N N d d dτ dτ dτ The minimization problem is posed as N X min min J = ∥ ϵ ( t ) ∥ (64) i N − d ˆ n × m 2 p ∈ R Θ ∈ R i =1 h i ⊤ The poles p = < 0 are identified by the outer loop minimization with p > − max ( k , . . . , k ) .

p · · · p 1 N i 1 M d

III. Nonlinear Stability and Control Derivative Estimation for Transonic

Truss-Braced Wing

A forced oscillation numerical experiment is conducted for the TTBW to support a low speed wind tunnel test at NASA LaRC 12-Ft Low-Speed Tunnel. The numerical experiment is performed using FUN3D for a 4 % scale wind tunnel model. The CFD mesh has 120 million nodes. A series of forced oscillations in roll, ◦ ◦ pitch, and yaw is performed for an amplitude of 10 over a range of the mean angle of attack from − 10 to ◦ ◦ ◦ 60 at three frequencies 0 . 25 Hz, 0 . 5 Hz, and 1 Hz, and two sideslip angles 0 and 15 . The corresponding reduced frequencies are k = 0 . 005 , 0 . 01 , and 0 . 02 , respectively. Since the reduced frequencies are small, the quasi-steady aerodynamic assumption is reasonable. The test condition is at a dynamic pressure of 4 psf corresponding to Mach 0 . 052 at sea level. The Turkel low-speed preconditioning scheme is used in the simulations.

A. Pitch Stability Derivative Estimation Figure 2 is a visualization of the streamlines over the TTBW in pitch oscillation at zero angle of attack.

A Fourier series decomposition is performed on all the time history data to determine the frequency con- tents. The data generally exhibit high frequency contents. The lift and pitching moment data contain predominantly odd harmonics, whereas the drag data contains even harmonics. In our previous work, the frequency-domain system identification of the linear and nonlinear stability derivatives is performed using th 10 an 8 -power sine wave to approximate the unsteady aerodynamic force and moment responses.

9 of 16 American Institute of Aeronautics and Astronautics ◦ Figure 2. Pitch Oscillation of Transonic Truss Braced-Wing at ¯ α = 0 and 0 . 25 Hz ˜ Figure 3 shows the plots of the unsteady lift coefficient C ( t ) versus the pitch angle θ ( t ) at the mean L zero angle of attack and the frequency of 0 . 25 Hz computed by FUN3D and by both the frequency-domain and time-domain system identification methods. The agreement between the CFD simulation data and the reconstructed data from both the system identification methods is excellent. The lift coefficient is nearly ◦ ◦ linear in the pitch angle or the angle of attack range from − 4 to 4 .

◦ ˜ Figure 3. C ( t ) vs. θ ( t ) at ¯ α = 0 and 0 . 25 Hz L ˜ Figures 4 shows the plots of the unsteady drag coefficient C ( t ) versus the pitch angle θ ( t ) . Excellent D agreement is noted between the FUN3D data and the reconstructed data from both the system identification methods.

10 of 16 American Institute of Aeronautics and Astronautics ◦ ˜ Figure 4. C ( t ) vs. θ ( t ) at ¯ α = 0 and 0 . 25 Hz D ˜ Figures 5 shows the plots of the unsteady pitching moment coefficient C ( t ) versus the pitch angle m θ ( t ) . Excellent agreement is noted between the FUN3D data and the reconstructed data from the system identification.

˜ ◦ Figure 5. C ( t ) vs. θ ( t ) at ¯ α = 0 and 0 . 25 Hz m Table 1 lists the linear and nonlinear pitch axis stability derivatives at 0 . 25 Hz estimated by the frequency- domain method. For comparison, Table 2 lists the corresponding stability derivatives estimated by the time-domain method. The agreement is excellent between the two methods. Comparing the two methods, the time-domain method is much simpler to implement than the frequency-domain method. The nonlinear spillover effect in the frequency-domain method is completely eliminated in the time-domain method. Note that the stability derivatives with respect to the pitch acceleration ˙ q are not estimated because there is only one frequency in the prescribed pitch angle.

It should be noted that if the linear stability derivatives were to be evaluated without considering the nonlinear spillover effect, the linear stability derivatives for the lift coefficient would have been evaluated 11 of 16 American Institute of Aeronautics and Astronautics n 1 2 3 4 5 6 7 8 C 7 . 474 − 0 . 6771 39 . 13 − 159 . 6 − 5611 1 . 283 e 4 7 . 442 e 4 − 2 . 768 e 5 n L α C 10 . 74 44 . 03 − 2252 − 5150 2 . 176 e 5 7 . 479 e 4 − 3 . 598 e 6 2 . 683 e 6 L n dθ dτ C 0 . 008150 2 . 651 5 . 741 − 248 . 3 − 928 . 9 2 . 146 e 4 2 . 308 e 4 − 3 . 829 e 5 D n α C 0 . 2082 5 . 825 − 38 . 76 1617 8821 − 1 . 564 e 5 − 2 . 734 e 5 3 . 064 e 6 D n dθ dτ C − 2 . 875 2 . 732 − 20 . 51 305 . 9 2745 − 2 . 205 e 4 − 4 . 082 e 4 4 . 825 e 5 m n α C − 123 . 5 − 43 . 37 625 . 3 6547 − 6 . 669 e 4 − 4 . 975 e 5 1 . 593 e 6 8 . 667 e 6 n m dθ dτ Table 1. Pitch Axis Stability Derivatives at 0 . 25 Hz Estimated by Frequency-Domain Method n 1 2 3 4 5 6 7 8 C 7 . 474 − 0 . 6324 39 . 13 − 164 . 9 − 5611 1 . 306 e 4 7 . 442 e 4 − 2 . 801 e 5 L n α C 10 . 63 44 . 03 − 2240 − 5150 2 . 170 e 5 7 . 479 e 4 − 3 . 589 e 6 2 . 683 e 6 L n dθ dτ C 0 . 008150 2 . 645 5 . 741 − 247 . 7 − 928 . 9 2 . 143 e 4 2 . 308 e 4 − 3 . 825 e 5 D n α C 0 . 2211 5 . 825 − 40 . 17 1617 8888 − 1 . 564 e 5 − 2 . 745 e 5 3 . 064 e 6 D n dθ dτ C − 2 . 875 2 . 724 − 20 . 51 306 . 9 2745 − 2 . 210 e 4 − 4 . 082 e 4 4 . 832 e 5 m n α C − 123 . 4 − 43 . 37 622 . 2 6547 − 6 . 655 e 4 − 4 . 975 e 5 1 . 591 e 6 8 . 667 e 6 n m dθ dτ Table 2. Pitch Axis Stability Derivatives at 0 . 25 Hz Estimated by Time-Domain Method using the Fourier series coefficients at the fundamental frequency given by 2 Npπ Z k k ˜ Q = C sin kτ dτ (65) 1 L N πθ p 0 0 2 Npπ Z k k ˜ S = C cos kτ dτ (66) 1 L N πθ p 0 In so doing, we would have obtained C = 6 . 264 , C = 29 . 87 , C = − 0 . 04264 , C = 0 . 2117 , C = L L D D m α q α q α − 2 . 382 , and C = − 123 . 2 as opposed to those shown in Table 1. Without considering the nonlinear m q ˜ ˜ spillover effect, it is tantamount to fitting a linear regression through the data sets of C ( t ) , C ( t ) , and L D ˜ C ( t ) . This would have resulted in errors in the linear stability derivative estimation.

m B. Elevator Control Derivative Estimation In our previous work, the elevator control surface oscillation is simulated in FUN3D for the Mach 0 . 8 condition at the trim angle of attack for the full-scale configuration. The control surface oscillation is ◦ ◦ prescribed by two truncated square waves of 1 and 20 amplitudes, each with 13 harmonics or M = 25 in the reduced frequency range from k = 0 . 02 to k = M k = 0 . 5 . The truncated square wave is a 0 max 0 multi-sine signal described by M X δ ( t ) = δ sin mk τ (67) e 0 0 m m =1 , 3 , ··· odd The nonlinear spillover effect can be isolated by the frequency-domain method for sinusoidal inputs. For complicated multi-sine inputs, the nonlinear response creates multiple overlapping harmonics that make it difficult to isolate the linear response in the frequency-domain method. The time-domain method resolves this difficulty.

◦ Figure 6(a) shows the pressure contour plots for Mach 0 . 8 and the elevator deflection at 20 . Shock- induced boundary layer separation at the hinge line is observed. This creates a highly nonlinear aerodynamic behavior in the forces and moments generated by the elevator.

12 of 16 American Institute of Aeronautics and Astronautics ◦ Figure 6. Instantaneous Pressure Contour for Mach 0 . 8 and Elevator Deflection at 20 th To estimate the control derivatives of the elevator, we employ the time-domain method with an 8 -degree nd polynomial and a 2 -order unsteady aerodynamic model having two poles. The plots of the unsteady lift ˜ ˜ ˜ coefficient C ( t ) , unsteady drag coefficient C ( t ) , and unsteady pitching moment coefficient C ( t ) versus L D m the elevator deflection computed by FUN3D and estimated by the time-domain system identification method are shown in Figures 7, 8, and 9, respectively. The agreement between the FUN3D data and the reconstructed data from the time-domain method is excellent. Table 3 lists the linear and nonlinear control derivatives of the elevator.

n 1 2 3 4 5 6 7 8 C 0 . 7306 − 0 . 03978 0 . 7373 3 . 685 − 11 . 34 − 66 . 15 18 . 37 3 . 561 L n δ e C 2 . 156 − 0 . 04702 1 . 207 0 . 3403 14 . 15 0 . 8486 − 83 . 88 − 24 . 27 n L dδ e dτ C 0 . 004987 0 . 4385 − 0 . 1458 0 . 4423 2 . 568 − 34 . 59 − 12 . 59 215 . 7 D n δ e C − 0 . 1006 − 0 . 02055 − 0 . 03471 1 . 581 0 . 5406 − 1 . 642 − 2 . 215 − 16 . 51 D n dδ e dτ C − 5 . 431 1 . 096 − 10 . 87 − 27 . 75 209 . 8 454 . 3 − 864 . 8 − 2410 n m δ e C 3 0 . 3612 − 9 . 962 − 0 . 8388 − 113 . 7 5 . 529 669 . 1 88 . 76 m n dδ e dτ d δ e a a b b p p 2 1 2 1 2 1 2 dτ C 0 . 1334 − 0 . 06934 − 0 . 375 − 1 . 875 − 0 . 4375 − 0 . 009575 − 0 . 2939 L C 0 . 01348 − 0 . 03003 0 . 04517 0 . 3613 − 0 . 2598 − 0 . 1092 − 0 . 1155 D C − 1 . 129 1 . 875 − 2 7 . 5 − 9 − 0 . 1595 − 0 . 3492 m Table 3. Elevator Control Derivatives Estimated by Time-Domain Method 13 of 16 American Institute of Aeronautics and Astronautics ˜ Figure 7. C ( t ) vs. δ ( t ) at Mach 0.8 L e ˜ Figure 8. C ( t ) vs. δ ( t ) at Mach 0.8 D e 14 of 16 American Institute of Aeronautics and Astronautics ˜ Figure 9. C ( t ) vs. δ ( t ) at Mach 0.8 m e ˜ ∂ C L Figure 10 shows the linear control derivative as a function of the reduced frequency k with and ∂δ e without the effect of unsteady aerodynamics. The magnitude of the frequency response is the square root of the real part squared and the imaginary part squared. The phase shift is the arc tangent of the ratio of the imaginary part to the real part. For a physically real flow, the phase shift should be negative which corresponds to a positive transport time delay. Without considering the unsteady aerodynamics, the lift curve slope indicates that the lift force would respond ahead of the elevator input which is non-physical. For k ≈ 0 where the quasi-steady aerodynamic assumption is valid, the differences in the linear control derivative ˜ ∂ C L is small, but as the frequency increases the differences become significant.

∂δ e ˜ ∂ C L Figure 10. vs. k at Mach 0.8 ∂δ e 15 of 16 American Institute of Aeronautics and Astronautics Conclusions This paper presents a time-domain system identification method for estimation of the linear and nonlinear stability and control derivatives for general aircraft applications. In our previous work on the frequency- domain system identification, a nonlinear spillover effect is identified whereby the nonlinearity produces contributions to the aerodynamic response at the input frequency. A procedure has been developed to remove the nonlinear spillover effect for sinusoidal inputs. For general inputs which are not sinusoids, the nonlinear spillover effect can cause difficulty to the frequency-domain method. The time-domain method resolves this difficulty effectively and is well suited for general applications. The time-domain method is applied to estimate the stability and control derivatives of the Transonic Truss Braced Wing aircraft configuration. Comparison to the frequency-domain method shows that the time-domain method produces virtually identical results but is simpler to implement for stability and control derivative estimation.

Acknowledgment

The authors wish to acknowledge NASA Advanced Air Transport Technology project for the funding support of this work. The authors also acknowledge Boeing Research and Technology and in particular Christopher Droney, Neal Harrison, Michael Beyar, Eric Dickey, and Anthony Sclafani, along with the NASA technical POC, Gregory Gatlin, for their research conducted under the NASA BAART contracts NNL10AA05B and NNL16AA04B. The research published in this paper is made possible by the technical data and wind tunnel test data furnished under these BAART contracts.

References

Wang, F., and Chen, L., “Numerical Prediction of Stability Derivatives for Complex Configurations,” Procedia Engineering 99 (2015) 1561-1575.

Ghoreyshi M., Bergeron K., Lofthouse, A., and Cummings R., “CFD Calculation of Stability and Control Derivatives For Ram-Air Parachutes,” AIAA Applied Aerodynamics Conference, AIAA-2016-1536, 2016.

Murmann, S., “A Reduced-Frequency Approach for Calculating Dynamic Derivatives,” AIAA Aerospace Sciences Meeting, AIAA-2005-0840, 2005.

Nguyen, N. and Xiong, J., “CFD-Based Frequency Domain Method for Dynamic Stability Derivative Estimation with Application to Transonic Truss-Braced Wing,” AIAA Aviation Forum, Applied Aerodynamics, AIAA-2022-3596, June 2022.

Nguyen, N. and Xiong, J., “Nonlinear Dynamic Control Derivative Analysis for Aircraft Application to Transonic Truss- Braced Wing,” AIAA Applied Aerodynamics Conference, AIAA-2024-0251, January 2024.

Nguyen, N. and Xiong, J., “Frequency Domain Method for Dynamic Control Derivative Estimation with Application to Transonic Truss-Braced Wing,” AIAA Applied Aerodynamics Conference, AIAA-2023-3945, June 2023.

Nguyen, N. and Xiong, J., “High-Fidelity Flight Dynamic Analysis of Transonic Truss-Braced Wing,” AIAA SciTech Forum, Applied Aerodynamics, AIAA-2023-0998, January 2023.

Harrison, N. A., Hoffman, K., Lazzara, D. S., Reichenbach, E. Y., Sclafani, A. J., and Droney, C. K., “Subsonic Ultra Green Aircraft Research: Phase IV Final Report - Volume I Mach 0.80 Transonic Truss-Braced Wing High-Speed Design Report,” NASA/CR-20220016017, October 2023.

Theodorsen, T., “General Theory of Aerodynamic Instability and the mechanism of Flutter”, NACA Report No. 496, 1949.

Nguyen, N., Xiong, J., and Webb, B., “Nonlinear System Identification of Aircraft Stability Derivatives in Numerical Forced Oscillation Experiment,” AIAA Applied Aerodynamics Conference, AIAA-2025-2184, January 2025.

Nguyen, N. and Xiong, J., “Computationally Efficient Frequency Domain Method for Dynamic Control Derivative Esti- mation with Application to Transonic Truss-Braced Wing,” AIAA Applied Aerodynamics Conference, AIAA-2024-4348, July 2024.

16 of 16 American Institute of Aeronautics and Astronautics