# An Integrated Test-Based Nonlinear Aeroelastic Analysis Method to Freeplay-Induced Limit-Cycle Oscillations

Authors: Yukai Sun; Zhigang Wu; Chao Yang; Xiangyu Wang; Chuansheng Wang
Journal: International Journal of Aerospace Engineering
 (2022)

---

## Abstract

The nonlinear limit-cycle oscillations induced by aileron rotational freeplay are investigated by numerical simulation and flutter wind tunnel tests. An integrated nonlinear aeroelastic analysis method and procedure based on the nonlinear system identification technology is proposed. Based on the characteristics of the structure, a low-speed nonlinear flutter wing-aileron model that considers both the linear flutter test and the nonlinear limit-cycle oscillation test with aileron rotational freeplay is designed. The wind tunnel tests are conducted with the model, and the stable limit-cycle oscillation is observed by accelerometers under the linear flutter boundary with a 0° angle of attack. Higher-order harmonic vibration is also observed with a low proportion compared to the fundamental oscillation, and the limit-cycle oscillation grows gradually with the increase of airspeed until flutter occurs. The comparison of numerical and experimental investigation shows high consistency.


---

## Body

## 1. Introduction

Freeplay that often occurs naturally in gears, bearings, control linkages, and actuators as a result of many factors during the aircraft’s life cycle is a source of aeroelastic instability and can cause limit-cycle oscillations (LCOs) [1]. Even though it is not as catastrophic as flutter, LCOs may lead to fatigue problems or structural failure, which would reduce the operating time and capacity. The motivation to research nonlinear aeroelastic phenomena, especially LCOs, on control-surface freeplay is quite significant known in plenty of literature from basic 3-degree-of-freedom (3-DOF) airfoil [2] to a real complex aircraft [3]. Several reviews, such as by Lee et al. [4], Dowell et al. [5, 6], Xiang et al. [7], and Yang et al. [8], draw the developments and advances in the aeroelastic analysis considering such structural nonlinearity.Investigations operated on an airfoil section would provide fundamental insights of LCOs induced by freeplay and had been widely studied. Tang and Dowell [9] explored a 3-DOF airfoil with freeplay in pitch and summarized that the limit-cycle motion was affected by airspeed, initial angle of attack, magnitude of gap, and initial conditions. Li et al. [2] investigated the LCOs of a 3-DOF airfoil based on the state-dependent Riccati equation method and designed an effective control law for suppressing the oscillation. Vasconcellos et al. [10] analyzed a 3-DOF airfoil with control-surface freeplay by a continuous approximation based on hyperbolic tangent representation. And an aeroelastic test was carried out on a foam-aluminum-fiberglass rigid wing in an open-circuit blower-like wind tunnel, as well as nonlinear time series analysis served to investigate the behavior and bifurcations of freeplay. A series of studies on flutter suppression technique was also evaluated by wind tunnel tests on a 3-DOF airfoil model [11–13].However, not only more complex phenomena may arise when simulating or testing a real-world aeroelastic system, the analyzing methods of which differ from ones of an airfoil, but also more essential roles the numerical computational burden plays under most circumstances. A numerical investigation on a three-dimensional multiple-actuated wing with freeplay in both leading edge and trailing edge control surfaces was performed by Huang et al. [14]. The nonlinear mathematical model established by the fictitious mass approach and two reduced-order optimal controllers were designed and proved to be effective. The simulation showed that the freeplay in control surfaces led to a lower flutter boundary. Recently, a considerable amount of literature has discussed around the theme of nonlinear aeroelastic wind tunnel tests. Bae et al. [15] performed linear and nonlinear aeroelastic analysis of a fighter-type wing with a control surface by both frequency-domain and time-domain methods. A change from a bending-flapping mode to a torsion-flapping mode was observed during the nonlinear LCO analysis, with different amplitudes of hinge angle, yielding the coexistence of different stable LCOs at the same airspeed. The nonlinear response also contains chaotic motion which cannot be observed directly by frequency-domain analysis. This research provides a general procedure of nonlinear aeroelastic analysis. Tang and Dowell [16, 17] designed an all-movable horizontal tail with a freeplay gap in the actuating mechanism and investigated the flutter response, nonlinear LCO, and gust response. The computational model was conducted and validated by experiments, whose aerodynamic force was calculated with a three-dimensional vortex lattice approach. The computation reached an excellent quantitative agreement to the experiments reasonably. Fichera and Ricci [18] investigated a T-tail with freeplay in the control chain by two numerical methods and wind tunnel tests. The comparison between the numerical and experimental results shows good agreement. Xie et al. [19] analyzed LCO problem of a fluttering plate in both two and three dimensions. Yang et al. [20] performed subsonic wind tunnel tests of folding fins with freeplay nonlinearity, which concluded that in a proper magnitude of gap, a larger freeplay would induce a higher divergent speed with a zero angle of attack. The numerical analysis was carried out by time-marching approach. Recently, Shi et al. and Liu et al. [21, 22] considered the control problem of a system with multisource complex disturbances, which may cause plenty of kinds of nonlinear phenomena, and the nonlinear behavior is analyzed and characterized [23–25]. Cooper et al. [26] have developed a pragmatic approach to integrate test-based system identification and finite element (FE) modeling into a nonlinear structure that is mainly concerned with structural dynamics.Up to now, system identification and prediction of LCO remain challenging to engineers dealing with complicated aeroelastic structures in the real world. In most literature, firstly, the assumption is that the prior knowledge of the nonlinear parameters of nonlinearity was taken for granted, and far too little attention has been paid to the modeling of a nonlinear aeroelastic system based on the experimental nonlinear system identification. Secondly, there are not many representative experimental wing-aileron model with aileron rotational freeplay that can provide stable LCO phenomena induced by freeplay in wind tunnel tests.The contributions of the present work are listed as follows: (1) An integrated nonlinear aeroelastic analysis methodology is proposed combining the substructure synthesis method and the nonlinear system identification technology, which is suitable for freeplay during engineering applications. (2) We also provide an integrated procedure to solve the nonlinear aeroelastic problems from model construction and nonlinear system identification to numerical calculation of structures including nonlinearity. (3) A low-speed nonlinear flutter experimental model that takes into account both the linear flutter test and the nonlinear LCO test with freeplay was designed. Stable LCO induced by freeplay can be obtained during the wind tunnel test, as well as the critical boundary of the linear flutter. The influence of aileron rotational freeplay on LCO was numerically and experimentally researched by numerical simulation and wind tunnel tests. The results show that the method in this paper has high accuracy.This paper has been organized in the following way. Section2 presents the methodology containing the nonlinear model construction and the numerical simulation methods. Section 3 describes the construction of an experimental nonlinear wind tunnel test model. Section 4 describes the detailed procedure on the ground vibration test and nonlinear system identification. Section 5 discusses the results of linear and nonlinear numerical analysis with a comparison to the linear flutter test and the nonlinear LCO test. Section 6 concludes the study.


## 2. Methodology of Integrated Test-Based Nonlinear Aeroelastic Analysis

A proper nonlinear mathematical model is necessary to predict the behavior of the nonlinear LCO induced by freeplay. The methodology of integrated test-based nonlinear aeroelastic analysis is depicted in Figure1. Firstly, the linear FE model is constructed using commercial software MSC.Nastran [27], and then is validated and updated based on the experimental results of the Ground Vibration Test (GVT). Then, the nonlinear system identification is performed to identify the nonlinear types and parameters. As for the nonlinear model, the extended free-interface component mode synthesis (CMS) method considered in this paper is an effective tool to establish the nonlinear governing equation. Finally, the time-marching approach is used to perform the numerical investigation on the nonlinear flutter and LCO induced by freeplay. In the present work, the rectangular wing-aileron test model was divided into two components, a main wing and an aileron connected by a nonlinear rotational hinge. Then, the nonlinear mathematical model was implemented by the combination of linear substructures and nonlinearity. In this section, some of the critical theories are briefly introduced.Figure 1
Procedure of integrated nonlinear aeroelastic analysis.


### 2.1. Extended Free-Interface Component Mode Synthesis Method

The CMS method involves dividing the structure into several substructures, obtaining reduced-order models of components, and then assembling a reduced-order model of entire structure [28]. The CMS method was first introduced by Hurty [29, 30] and extended by Craig and Bampton [31] as the fixed-interface CMS method with widely engineering applications. Then, an alternative type of coupling procedure called free-interface CMS method was proposed by Hou [32] and developed by Rubin [33] and Craig and Chang [34]. Both fixed-interface and free-interface CMS methods were extended for the system with structural nonlinearities by Kim et al. [35] and Yang et al. [20], respectively.The structural dynamic equation of an undamped multidegree of freedom substructure in physical coordinates constructed by the finite-element method is [28]
(1)Mu¨+Ku=f=BTfJ,where M and K are the mass matrix and stiffness matrix, respectively. f represents the external force. Depending on the connecting condition of nodes in FE model, the nodes can be divided into two different types, namely, inner node, represented by subscript I, and boundary node, represented by subscript J. Thus the external force also contains the inner node force fI and boundary node force fJ, while in the procedure of the CMS method, fI=0 is often considered. B is the projection matrix from all nodes to boundary nodes.Take the influence of high-order deleted free-interface normal modes into consideration by introducing the residual flexibility attachment modes [28] Ψd. The displacement response u can be represented by the kept free-interface normal modes Φk and Ψd as
(2)u=Φkqk+ΨdfJ,where qk is the generalized coordinate of Φk. The residual flexibility attachment modes are written as
(3)Ψd=ΓTXK¯−1XTΓBTfJ−ΦkΛk−1ΦkTBT,where Γ=I−MΦRΦRT is the projection matrix to balance the effect of inertia force and ΦR is the rigid-body model matrix. X is a constraint matrix, by which the rigid motions are restricted and therefore K¯=XTKX is nonsingular.Then, the main wing and aileron can be assembled with the compatibility conditions on the interface concerning displacements and boundary force as(4)uJr,W=uJr,A,uJe,W=uJe,A+δ,fJ,W=fJ,A,where the subscripts W and A imply the component of the main wing and aileron, respectively. And subscripts Jr and Je represent the rigid connection and the elastically connected situation on the interface. The boundary force can be expanded as fJ=fJrTfJeTT.The displacement responses of the main wing and aileron can be written as Equation (2).
(5)uW=Φk,Wqk,W+Ψd,WfJ,W,uA=Φk,Aqk,A+Ψd,AfJ,A.Substituting Equation (4) to Equation (5) and the entire assemble generalized coordinates can be transformed into
(6)qk,WfJ,Wqk,AfJ,A=Tqk,Wqk,Aδ,where T is the transform matrix.The transformation of Equation (6) is applied to the integrated dynamic equation of the rectangular wing with aileron deriving the reduced-order equation as
(7)Mqq¨+Kqq=TTfJ,WA,where
(8)Mq=TTM~W00M~AT,(9)Kq=TTK~W00K~AT,(10)q=qk,Wqk,Aδ,(11)fJ,WA=Φk,WTBWTfJ,WΨd,WTBWTfJ,WΦk,ATBATfJ,AΨk,ATBATfJ,A,where M~ and K~ are generalized mass matrix and stiffness matrix, respectively. fJ,WA is the generalized force. However, it has been proved that the generalized force corresponding to the coordinates qk,W and qk,A are zero. The generalized force related to δ is the elastic restoring force to the relative displacement δ as
(12)TTfJ,WA=fJδ=00−fJe,W.In general, the integrated dynamic equation of the entire structure can be written as(13)Mqq¨+Kqq+fJδ=0.For the linear system, the generalized force satisfiesfJe,W=fJe,A=keδ, where ke is the linear stiffness. However, for nonlinear connections, the restoring force is a nonlinear function of the relative displacement δ, i.e., fJe,W=knonδ·δ. Substituting the restoring force expression to Equation (13) and therefore the entire structure is reduced under generalized coordinates q.


### 2.2. Nonlinear System Identification

The accurate representation of nonlinear features in the FE model would be of extreme benefit during the process of nonlinear simulation. However, principally for complicated aerospace structures, the structural nonlinearities may not be measured directly by static methods and characterized accurately. Hence, nonlinear system identification is important in many structural dynamic applications [26]. The overall objective of the nonlinear system identification is to distinguish and describe the nonlinearities mathematically from experimental data.With the wide approval from researchers and structural engineers, Kerschen et al. [36] summarized that the nonlinear system identification involves three steps, namely, detection, characterization, and parameter estimation. The detection and characterization stages were not relevant in this paper based on the fact that nonlinear rocker arm is responsible for the freeplay nonlinearity.In the parameter estimation stage, time-domain nonlinear subspace identification (TNSI) was performed to the parameter estimation of freeplay. Subspace-based methods are commonly recognized as powerful identification tools for linear and time-invariant systems because they are noniterative and naturally applicable to multi-input-multi-output (MIMO) systems [37]. The subspace-based algorithms have been applied to a wide variety of systems and firstly extended to nonlinear systems by Lacy and Bernstein [38]. Recently, the TNSI implementation was accomplished by Marchesiello and Garibaldi [39] using the idea of nonlinear identification through feedback of the outputs (NIFO) [40], similar to the conditioned reverse path (CRP) [41] methods. TNSI methods also used to identify the clearance nonlinearities [42], yielding superior accuracy.As illustrated in Figure2, TNSI methods rely on the feedback interpretation of nonlinear vibrating systems governed by Equation (13) with external force, which recast as follows:
(14)Mqq¨+Kqq+fJδ=ft,where fJδ representing the nonlinear restoring force is modeled using a summation of the form in the nonlinear identification scheme [40]
(15)fJ=∑j=1pμjLnjgjt.Figure 2
Feedback interpretation of nonlinear vibrating systems.The nonlinear term is expressed as the sum ofp components. Lnj indicates the location of the nonlinear element and μj is an unknown nonlinear coefficient corresponding with the jth nonlinear functional form gjt. By moving the nonlinear term of Equation (14) to the right-hand side,
(16)Mqq¨t+Kqqt=ft−∑j=1pμjLnjgjt.Then, the state-space formulation of the equation of motion corresponds to an input vector asu=ftT−g1t⋯−gptT.
(17)q̇q¨=Acqq̇+Bcu,y=Cqq̇+Du,where
(18)Ac=0I−Mq−1Kq0,Bc=00⋯0Mq−1Mq−1μ1Ln1⋯Mq−1μpLnp,C=I0,D=00⋯0.The nonlinear identification procedure is based on the computation of the system “extended” frequency response function (FRF) matrixHEω.
(19)HEω=D+CiωI−Ac−1Bc.The matrixHEω can be recast as
(20)HEω=0+I0P11P12P21P22,(21)00⋯0Mq−1Mq−1μ1Ln1⋯Mq−1μpLnp=P12Mq−1Iμ1Ln1⋯μpLnp,where
(22)N=iωI−Ac=N11N12N21N22,P=iωI−Ac−1=P11P12P21P22.Since the block matrixP12 is
(23)P12=1iωiωI+Mq−1Kq1iω−1=Mq−1Kq−ω2I−1=Kq−ω2Mq−1−1Mq,and the underlying linear system FRF matrix is
(24)Hω=Kq−ω2Mq−1,then Equation (21) finally becomes
(25)HEω=HHμ1Ln1⋯HμpLnp,and especially when ω=0,
(26)HE0=D−CAc−1Bc=Kq−1Kμ1Ln1⋯KμpLnp,from where the unknown coefficient can be estimated by utilizing Equation (26).


### 2.3. Nonlinear Time-Marching Approach

The aeroelastic system described using the extended free-interface CMS method is(27)Mqq¨+Kqq+fJδ=12ρV2Qq,where Q is the generalized aerodynamic force influence coefficient matrix computed in the reduced-frequency k-domain by double-lattice method (DLM), ρ is the air density, and V is the airspeed. If the whole system can be regarded as a linear one, i.e., fJδ=0, then the standard p-k and v-g methods can be performed to acquire the flutter characteristics.The structural freeplay, however, may defeat the basic assumptions of the traditional linear flutter analysis methods. The freeplay restoring force is depicted in Figure3, where δ0 is the rotational gap of wing-aileron hinge and k0 is the outer stiffness of freeplay.Figure 3
Restoring force of freeplay nonlinearity.The time-marching approach is performed to consider the nonlinearity, which, however, involved more considerable burden on computational capability. In Equation (27), the generalized aerodynamic forces can be fitted into a state-space system by using minimum-state rational aerodynamic approximations [43] in the following form:
(28)Q~p=A0+A1p+A2p2+DIp−R−1Ep,where p is the nondimensional complex Laplace variable and A0, A1, A2, D, and E are real number aerodynamic coefficient matrices, respectively. R is the aerodynamic root matrix, the order of which is equal to the number of aerodynamic state variables na. Then, in time-domain, the aeroelastic equation can be written as
(29)Mqq¨+Kqq+fJδ=12ρV2A0q+12ρVLA1q̇+12ρL2A2q¨+12ρV2DsI−RVL−1Eq̇.Introducing an augmenting aerodynamic state vector of dimensionna defined by its Laplace transform
(30)y=sI−RVL−1Esq,and then, Equation (30) can be recast into
(31)ẏ=RVLy+Eq̇.Combining Equation (29) and Equation (31), finally, the state-space aeroelastic equation is established
(32)ddtqq̇y=0I0KSCSDS0ERVLqq̇y−0MS−1fJδ0,where the coefficient matrices are listed as follows:
(33)MS=Mq−12ρL2A2,CS=MS−112ρVLA1,KS=MS−112ρV2A0−Kq,DS=12ρV2MS−1D.For a given airspeed and initial condition, then the determined state-space aeroelastic system described by Equation (32) can be directly simulated by direct time integration method, such as Runge-Kutta scheme or Newmark family, by which the time-domain responses of the generalized coordinates q can be obtained. The responses under physical coordinates then can be transformed from q by the usage of Φk and Ψd in Equation (5) and the transform matrix T in Equation (6). The amplitude and portraits of limit-cycle motion can also be depicted by the simulation, and fast Fourier transform (FFT) analysis of the responses can illustrate the frequency-domain characteristics of LCO. In this paper, the 4-5 order Runge-Kutta simulation method ode45 in MATLAB is used in numerical investigation.


## 2.1. Extended Free-Interface Component Mode Synthesis Method

The CMS method involves dividing the structure into several substructures, obtaining reduced-order models of components, and then assembling a reduced-order model of entire structure [28]. The CMS method was first introduced by Hurty [29, 30] and extended by Craig and Bampton [31] as the fixed-interface CMS method with widely engineering applications. Then, an alternative type of coupling procedure called free-interface CMS method was proposed by Hou [32] and developed by Rubin [33] and Craig and Chang [34]. Both fixed-interface and free-interface CMS methods were extended for the system with structural nonlinearities by Kim et al. [35] and Yang et al. [20], respectively.The structural dynamic equation of an undamped multidegree of freedom substructure in physical coordinates constructed by the finite-element method is [28]
(1)Mu¨+Ku=f=BTfJ,where M and K are the mass matrix and stiffness matrix, respectively. f represents the external force. Depending on the connecting condition of nodes in FE model, the nodes can be divided into two different types, namely, inner node, represented by subscript I, and boundary node, represented by subscript J. Thus the external force also contains the inner node force fI and boundary node force fJ, while in the procedure of the CMS method, fI=0 is often considered. B is the projection matrix from all nodes to boundary nodes.Take the influence of high-order deleted free-interface normal modes into consideration by introducing the residual flexibility attachment modes [28] Ψd. The displacement response u can be represented by the kept free-interface normal modes Φk and Ψd as
(2)u=Φkqk+ΨdfJ,where qk is the generalized coordinate of Φk. The residual flexibility attachment modes are written as
(3)Ψd=ΓTXK¯−1XTΓBTfJ−ΦkΛk−1ΦkTBT,where Γ=I−MΦRΦRT is the projection matrix to balance the effect of inertia force and ΦR is the rigid-body model matrix. X is a constraint matrix, by which the rigid motions are restricted and therefore K¯=XTKX is nonsingular.Then, the main wing and aileron can be assembled with the compatibility conditions on the interface concerning displacements and boundary force as(4)uJr,W=uJr,A,uJe,W=uJe,A+δ,fJ,W=fJ,A,where the subscripts W and A imply the component of the main wing and aileron, respectively. And subscripts Jr and Je represent the rigid connection and the elastically connected situation on the interface. The boundary force can be expanded as fJ=fJrTfJeTT.The displacement responses of the main wing and aileron can be written as Equation (2).
(5)uW=Φk,Wqk,W+Ψd,WfJ,W,uA=Φk,Aqk,A+Ψd,AfJ,A.Substituting Equation (4) to Equation (5) and the entire assemble generalized coordinates can be transformed into
(6)qk,WfJ,Wqk,AfJ,A=Tqk,Wqk,Aδ,where T is the transform matrix.The transformation of Equation (6) is applied to the integrated dynamic equation of the rectangular wing with aileron deriving the reduced-order equation as
(7)Mqq¨+Kqq=TTfJ,WA,where
(8)Mq=TTM~W00M~AT,(9)Kq=TTK~W00K~AT,(10)q=qk,Wqk,Aδ,(11)fJ,WA=Φk,WTBWTfJ,WΨd,WTBWTfJ,WΦk,ATBATfJ,AΨk,ATBATfJ,A,where M~ and K~ are generalized mass matrix and stiffness matrix, respectively. fJ,WA is the generalized force. However, it has been proved that the generalized force corresponding to the coordinates qk,W and qk,A are zero. The generalized force related to δ is the elastic restoring force to the relative displacement δ as
(12)TTfJ,WA=fJδ=00−fJe,W.In general, the integrated dynamic equation of the entire structure can be written as(13)Mqq¨+Kqq+fJδ=0.For the linear system, the generalized force satisfiesfJe,W=fJe,A=keδ, where ke is the linear stiffness. However, for nonlinear connections, the restoring force is a nonlinear function of the relative displacement δ, i.e., fJe,W=knonδ·δ. Substituting the restoring force expression to Equation (13) and therefore the entire structure is reduced under generalized coordinates q.


## 2.2. Nonlinear System Identification

The accurate representation of nonlinear features in the FE model would be of extreme benefit during the process of nonlinear simulation. However, principally for complicated aerospace structures, the structural nonlinearities may not be measured directly by static methods and characterized accurately. Hence, nonlinear system identification is important in many structural dynamic applications [26]. The overall objective of the nonlinear system identification is to distinguish and describe the nonlinearities mathematically from experimental data.With the wide approval from researchers and structural engineers, Kerschen et al. [36] summarized that the nonlinear system identification involves three steps, namely, detection, characterization, and parameter estimation. The detection and characterization stages were not relevant in this paper based on the fact that nonlinear rocker arm is responsible for the freeplay nonlinearity.In the parameter estimation stage, time-domain nonlinear subspace identification (TNSI) was performed to the parameter estimation of freeplay. Subspace-based methods are commonly recognized as powerful identification tools for linear and time-invariant systems because they are noniterative and naturally applicable to multi-input-multi-output (MIMO) systems [37]. The subspace-based algorithms have been applied to a wide variety of systems and firstly extended to nonlinear systems by Lacy and Bernstein [38]. Recently, the TNSI implementation was accomplished by Marchesiello and Garibaldi [39] using the idea of nonlinear identification through feedback of the outputs (NIFO) [40], similar to the conditioned reverse path (CRP) [41] methods. TNSI methods also used to identify the clearance nonlinearities [42], yielding superior accuracy.As illustrated in Figure2, TNSI methods rely on the feedback interpretation of nonlinear vibrating systems governed by Equation (13) with external force, which recast as follows:
(14)Mqq¨+Kqq+fJδ=ft,where fJδ representing the nonlinear restoring force is modeled using a summation of the form in the nonlinear identification scheme [40]
(15)fJ=∑j=1pμjLnjgjt.Figure 2
Feedback interpretation of nonlinear vibrating systems.The nonlinear term is expressed as the sum ofp components. Lnj indicates the location of the nonlinear element and μj is an unknown nonlinear coefficient corresponding with the jth nonlinear functional form gjt. By moving the nonlinear term of Equation (14) to the right-hand side,
(16)Mqq¨t+Kqqt=ft−∑j=1pμjLnjgjt.Then, the state-space formulation of the equation of motion corresponds to an input vector asu=ftT−g1t⋯−gptT.
(17)q̇q¨=Acqq̇+Bcu,y=Cqq̇+Du,where
(18)Ac=0I−Mq−1Kq0,Bc=00⋯0Mq−1Mq−1μ1Ln1⋯Mq−1μpLnp,C=I0,D=00⋯0.The nonlinear identification procedure is based on the computation of the system “extended” frequency response function (FRF) matrixHEω.
(19)HEω=D+CiωI−Ac−1Bc.The matrixHEω can be recast as
(20)HEω=0+I0P11P12P21P22,(21)00⋯0Mq−1Mq−1μ1Ln1⋯Mq−1μpLnp=P12Mq−1Iμ1Ln1⋯μpLnp,where
(22)N=iωI−Ac=N11N12N21N22,P=iωI−Ac−1=P11P12P21P22.Since the block matrixP12 is
(23)P12=1iωiωI+Mq−1Kq1iω−1=Mq−1Kq−ω2I−1=Kq−ω2Mq−1−1Mq,and the underlying linear system FRF matrix is
(24)Hω=Kq−ω2Mq−1,then Equation (21) finally becomes
(25)HEω=HHμ1Ln1⋯HμpLnp,and especially when ω=0,
(26)HE0=D−CAc−1Bc=Kq−1Kμ1Ln1⋯KμpLnp,from where the unknown coefficient can be estimated by utilizing Equation (26).


## 2.3. Nonlinear Time-Marching Approach

The aeroelastic system described using the extended free-interface CMS method is(27)Mqq¨+Kqq+fJδ=12ρV2Qq,where Q is the generalized aerodynamic force influence coefficient matrix computed in the reduced-frequency k-domain by double-lattice method (DLM), ρ is the air density, and V is the airspeed. If the whole system can be regarded as a linear one, i.e., fJδ=0, then the standard p-k and v-g methods can be performed to acquire the flutter characteristics.The structural freeplay, however, may defeat the basic assumptions of the traditional linear flutter analysis methods. The freeplay restoring force is depicted in Figure3, where δ0 is the rotational gap of wing-aileron hinge and k0 is the outer stiffness of freeplay.Figure 3
Restoring force of freeplay nonlinearity.The time-marching approach is performed to consider the nonlinearity, which, however, involved more considerable burden on computational capability. In Equation (27), the generalized aerodynamic forces can be fitted into a state-space system by using minimum-state rational aerodynamic approximations [43] in the following form:
(28)Q~p=A0+A1p+A2p2+DIp−R−1Ep,where p is the nondimensional complex Laplace variable and A0, A1, A2, D, and E are real number aerodynamic coefficient matrices, respectively. R is the aerodynamic root matrix, the order of which is equal to the number of aerodynamic state variables na. Then, in time-domain, the aeroelastic equation can be written as
(29)Mqq¨+Kqq+fJδ=12ρV2A0q+12ρVLA1q̇+12ρL2A2q¨+12ρV2DsI−RVL−1Eq̇.Introducing an augmenting aerodynamic state vector of dimensionna defined by its Laplace transform
(30)y=sI−RVL−1Esq,and then, Equation (30) can be recast into
(31)ẏ=RVLy+Eq̇.Combining Equation (29) and Equation (31), finally, the state-space aeroelastic equation is established
(32)ddtqq̇y=0I0KSCSDS0ERVLqq̇y−0MS−1fJδ0,where the coefficient matrices are listed as follows:
(33)MS=Mq−12ρL2A2,CS=MS−112ρVLA1,KS=MS−112ρV2A0−Kq,DS=12ρV2MS−1D.For a given airspeed and initial condition, then the determined state-space aeroelastic system described by Equation (32) can be directly simulated by direct time integration method, such as Runge-Kutta scheme or Newmark family, by which the time-domain responses of the generalized coordinates q can be obtained. The responses under physical coordinates then can be transformed from q by the usage of Φk and Ψd in Equation (5) and the transform matrix T in Equation (6). The amplitude and portraits of limit-cycle motion can also be depicted by the simulation, and fast Fourier transform (FFT) analysis of the responses can illustrate the frequency-domain characteristics of LCO. In this paper, the 4-5 order Runge-Kutta simulation method ode45 in MATLAB is used in numerical investigation.


## 3. Nonlinear Wing-Aileron Test Benchmark

A low-speed nonlinear flutter experimental model was designed to investigate the self-excited aeroelastic LCO induced by control-surface freeplay. As shown in Figure4(a), this experimental model includes two main components: a main wing and an aileron. The chord of the main wing is Cm=400 mm and the span is bm=1400 mm. A single spar with seven ribs composed the internal metal structure of the main wing, as shown in Figure 4(b). An actuator and an angle sensor are located at the tip two of the main wing beam. The whole metal structure is made of aluminum.Figure 4
Freeplay nonlinear wind tunnel flutter test benchmark.
(a)
The main wing of the nonlinear wind tunnel flutter test benchmark(b)
Inner structures of the nonlinear wind tunnel flutter test benchmarkSince in most instances, it is very desirable for the model stiffness properties to be supplied by an internal structure, and the external shell must be capable of supporting the air loads while distorted without modifying seriously the stiffness properties of the model. Thus, the external shape of the rectangular wing is divided into 7 pieces of 200 mm wide wing boxes, which are enforced by 3D printing wing shells thickening 2 mm. Each wing box consists of a pair of corresponding robs with NACA 0015 airfoil. The wing boxes contributed no markedly to the bending and torsional stiffness since they were not attached to the main wing beam directly. There are 2 different types of wing sections based on the initialization of aileron as in Figure4(a). The aileron is mounted in the place of two wing boxes on the wingtip with the chord Ca=114 mm and the span ba=400 mm. The aileron cover is made of balsa wood. The total mass of the experimental model is 5.828 kg.Linear and nonlinear hinge between the main wing and aileron is implemented by linear and freeplay rotational arms as in Figures5(a) and 5(b), respectively. The freeplay arms consist of actuator flange, bearing, and rocker arm. The bearing guarantees the free rotation of the rocker arm. The gap between groove and block after assembling all three parts above provides the freeplay at aileron rotational hinge. The flange is mounted on the actuator by bolts, then. The aileron rotational stiffness is simulated by a square-section elastic beam made of aluminum. The base of the elastic beam is connected with the aileron beam. The actual gap and rotational stiffness of freeplay are discussed later.Figure 5
The design of rotational arms of the actuator.
(a)
Linear(b)
NonlinearThe FE model was conducted mathematically as presented in Figure6, which is generated using commercial MSC. Nastran software. The main wing beam and robs were conducted using beam finite elements, and the 3D printed shells of wing boxes were retrieved into several concentrated masses attached to the rob for a better computational speed. The aileron was also divided into concentrated masses mounted on the aileron beam. To simulate the rotational motion of aileron, three connections were conducted with the rotational direction of axis Y free, while the DOFs of other directions were connected rigidly. One torsional spring was established at the middle connection.Figure 6
The FE model of the nonlinear wind tunnel flutter test benchmark.To construct the nonlinear computational model using the extended CMS method, the rectangular wing was divided into the main wing and aileron. By the procedure of the extended CMS method, the linear FE model of the main wing was reduced to 6 DOFs under the description of kept modal coordinates, and the mathematical model of the aileron was described using 6 rigid modes. However, the nonlinear rotational DOF was kept while assembling the substructures, the rotational stiffness of which was set to 2.90 N·m/rad. In total, the nonlinear rectangular wing flutter test model was represented by a reduced-order model (ROM) containing 13 DOFs. The comparison of the original FE model and the ROM model is shown in Table1.Table 1
Modal analysis of ROM and the original FE model.
Mode shapeROMNastranError (%)Frequency (Hz)Frequency (Hz)1st bending mode1.891.890.00In-plane mode9.199.180.11Aileron rotating mode11.6911.700.092nd bending mode12.6012.600.001st torsion mode19.0619.060.003rd bending mode28.7028.700.00


## 4. Ground Vibration Test and Nonlinear System Identification

### 4.1. Validation of the Main Wing

As illustrated in Figure7, the first experimental test was conducted on the main wing based on burst random excitation. The main wing was fixed to the ground by a support base. The test model was constructed by 26 measure points (Point.2-27) and a fixed point (Point.1) located at the root of the main wing. PCB-333B30 accelerometers were used to measure the response of the main wing at these 26 measure locations. Point.8 on the second bottom wing box was excited by an excitation shaker.Figure 7
The test model of the main wing (Point.1 is the fixed point and Point.2 to 27 are measured points, while the location of excitation is Point.8).The GVTs were performed using the Spectral Test module in LMS Test.Lab. Based on the design and preliminary analysis in Section3, it is appropriate to select the excitation frequency range of interest of 0-64 Hz, which satisfies the engineering requirements. The main wing was excited at several excitation levels. The purpose of the test was to observe if there was any shift of FRFs. A selection of the FRFs and coherence functions of interests obtained from different levels of excitation test are presented in Figure 8, showing hardly any shifts in the natural frequencies or distortion in the response curve. Therefore, the main wing could be regarded as a linear entity within the exciting range.Figure 8
The sum FRF measured at different level excitation of the main wing.Then, the initial FE model of the main wing was correlated and updated using the GVT test results by using FEMTools [44] developed by Dynamic Design Solutions (DDS). The updated natural frequencies are compared in Table 2, and the first 4 modes obtained from the GVTs are shown in Figure 9. The cross-modal assurance criterion (MAC) matrix in Figure 10 shows an excellent correlation between these four modes.Table 2
Main wing natural frequencies and mode shapes.
Updated FE modelExperimental modelError (%)Mode shapeFrequency (Hz)Mode shapeFrequency (Hz)1st bending mode1.921st bending mode1.930.522nd bending mode11.962nd bending mode11.672.481st torsion mode17.421st torsion mode16.893.143rd bending mode29.033rd bending mode28.183.02Figure 9
First four mode shapes of GVTs (blue dot: updated FE model; red dot: test model).
(a)
1st bending mode(b)
2nd bending mode(c)
1st torsion mode(d)
3rd bending modeFigure 10
Cross MAC of numerical-experimental model.After the validation of the main wing, the aileron of the flutter test model had also been validated. However, the 1st elastic mode of aileron is located far away outside of the frequency range of interest with a natural frequency of 84.2 Hz. Hence, the aileron can be regarded as a rigid component that able to meet the satisfaction of the requirement of accuracy.


### 4.2. Nonlinear System Identification

The nonlinear system identification experimental setup is presented in Figure11(a). To eliminate the influence of main wing motion and increase the accuracy of parameter estimation, the aileron was mounted on a rigid beam with supports fixed to the ground. Displacement, velocity, and acceleration were measured at the trailing edge of aileron with a sampling frequency as 1 kHz. Figure 11(b) illustrates the details of the mechanical linkage between the aileron and actuator on the main wing.Figure 11
Experimental setup of nonlinear system identification of freeplay.
(a)
Experimental setup(b)
Details of the mechanical leakageThe aileron was excited by a shaker at the location of measured point using random excitation since all the frequencies of interests must be excited sufficiently. For the parameter estimation of freeplay using TNSI method, eight nonlinear terms were chosen defined as(34)gjt=0,forδt≤δj,gjt=μjδ−δj·sgnδ,forδt>δj,where δj is a user-defined gap term and μj is the strength of the jth nonlinearity.In order to estimate the correct rotational gap at the wing-aileron hinge. Firstly, the positive displacement range0,maxδt was divided into eight equally spaced intervals δj−1,δj with δ0=0. The TNSI method detected the most important contributions at the first estimation and then narrows the displacement range and repeated until the results meet the satisfaction. Figure 12 shows the FRF of the underlying linear system obtained by TNSI method. On the comparison, FRF was also estimated by utilizing the traditional H2 estimation method. The contributions of each nonlinear term are presented in Figure 13(a), and the reconstruction of nonlinear restoring force is shown in Figure 13(b). The freeplay estimation by TNSI method corresponds to the zone near the correct rotational gap of aileron, i.e., the gap between the groove and block in Figure 5. The rotational gap was also measured using a static method by caliper and represented between two red dashed lines in Figure 13(b). The estimation by TNSI method achieved a significant agreement to the true freeplay, and the parameters estimated are presented in Table 3.Figure 12
FRF of aileron rotational motion (dashed line: traditionalH2 estimation; solid line: TNSI method).Figure 13
Parameter estimation by TNSI method.
(a)
Force contributions of nonlinear terms(b)
Nonlinear restoring forceTable 3
Nonlinear parameters of FE model and experimental model.
ParametersTrueEstimationError (%)Freeplay gap (°)1.701.731.76Linear-part stiffness (N·m/rad)2.902.910.34After the validations of the main wing and aileron, as well as the nonlinear system identification of freeplay, the construction of the nonlinear governing equation is finally in reach. By setting the rotational gap to zero using linear arm as in Figure5(a), the nonlinear system retrieved to a linear one. The result of the modal analysis shows a great agreement between the ROM and the test model as summarized in Table 4. The first four out-plane modes are shown in Figure 14.Table 4
Comparison of natural frequency between the ROM and GVTs.
Mode shapeGVTs (Hz)ROM (Hz)Error (%)1st bending mode1.931.892.072nd bending mode11.9511.692.18Aileron rotating mode12.5812.600.161st torsion mode19.5019.062.26Figure 14
First four mode shapes of the linear FE model.
(a)
1st bending mode(b)
2nd bending mode(c)
Aileron rotating mode(d)
1st torsion mode


## 4.1. Validation of the Main Wing

As illustrated in Figure7, the first experimental test was conducted on the main wing based on burst random excitation. The main wing was fixed to the ground by a support base. The test model was constructed by 26 measure points (Point.2-27) and a fixed point (Point.1) located at the root of the main wing. PCB-333B30 accelerometers were used to measure the response of the main wing at these 26 measure locations. Point.8 on the second bottom wing box was excited by an excitation shaker.Figure 7
The test model of the main wing (Point.1 is the fixed point and Point.2 to 27 are measured points, while the location of excitation is Point.8).The GVTs were performed using the Spectral Test module in LMS Test.Lab. Based on the design and preliminary analysis in Section3, it is appropriate to select the excitation frequency range of interest of 0-64 Hz, which satisfies the engineering requirements. The main wing was excited at several excitation levels. The purpose of the test was to observe if there was any shift of FRFs. A selection of the FRFs and coherence functions of interests obtained from different levels of excitation test are presented in Figure 8, showing hardly any shifts in the natural frequencies or distortion in the response curve. Therefore, the main wing could be regarded as a linear entity within the exciting range.Figure 8
The sum FRF measured at different level excitation of the main wing.Then, the initial FE model of the main wing was correlated and updated using the GVT test results by using FEMTools [44] developed by Dynamic Design Solutions (DDS). The updated natural frequencies are compared in Table 2, and the first 4 modes obtained from the GVTs are shown in Figure 9. The cross-modal assurance criterion (MAC) matrix in Figure 10 shows an excellent correlation between these four modes.Table 2
Main wing natural frequencies and mode shapes.
Updated FE modelExperimental modelError (%)Mode shapeFrequency (Hz)Mode shapeFrequency (Hz)1st bending mode1.921st bending mode1.930.522nd bending mode11.962nd bending mode11.672.481st torsion mode17.421st torsion mode16.893.143rd bending mode29.033rd bending mode28.183.02Figure 9
First four mode shapes of GVTs (blue dot: updated FE model; red dot: test model).
(a)
1st bending mode(b)
2nd bending mode(c)
1st torsion mode(d)
3rd bending modeFigure 10
Cross MAC of numerical-experimental model.After the validation of the main wing, the aileron of the flutter test model had also been validated. However, the 1st elastic mode of aileron is located far away outside of the frequency range of interest with a natural frequency of 84.2 Hz. Hence, the aileron can be regarded as a rigid component that able to meet the satisfaction of the requirement of accuracy.


## 4.2. Nonlinear System Identification

The nonlinear system identification experimental setup is presented in Figure11(a). To eliminate the influence of main wing motion and increase the accuracy of parameter estimation, the aileron was mounted on a rigid beam with supports fixed to the ground. Displacement, velocity, and acceleration were measured at the trailing edge of aileron with a sampling frequency as 1 kHz. Figure 11(b) illustrates the details of the mechanical linkage between the aileron and actuator on the main wing.Figure 11
Experimental setup of nonlinear system identification of freeplay.
(a)
Experimental setup(b)
Details of the mechanical leakageThe aileron was excited by a shaker at the location of measured point using random excitation since all the frequencies of interests must be excited sufficiently. For the parameter estimation of freeplay using TNSI method, eight nonlinear terms were chosen defined as(34)gjt=0,forδt≤δj,gjt=μjδ−δj·sgnδ,forδt>δj,where δj is a user-defined gap term and μj is the strength of the jth nonlinearity.In order to estimate the correct rotational gap at the wing-aileron hinge. Firstly, the positive displacement range0,maxδt was divided into eight equally spaced intervals δj−1,δj with δ0=0. The TNSI method detected the most important contributions at the first estimation and then narrows the displacement range and repeated until the results meet the satisfaction. Figure 12 shows the FRF of the underlying linear system obtained by TNSI method. On the comparison, FRF was also estimated by utilizing the traditional H2 estimation method. The contributions of each nonlinear term are presented in Figure 13(a), and the reconstruction of nonlinear restoring force is shown in Figure 13(b). The freeplay estimation by TNSI method corresponds to the zone near the correct rotational gap of aileron, i.e., the gap between the groove and block in Figure 5. The rotational gap was also measured using a static method by caliper and represented between two red dashed lines in Figure 13(b). The estimation by TNSI method achieved a significant agreement to the true freeplay, and the parameters estimated are presented in Table 3.Figure 12
FRF of aileron rotational motion (dashed line: traditionalH2 estimation; solid line: TNSI method).Figure 13
Parameter estimation by TNSI method.
(a)
Force contributions of nonlinear terms(b)
Nonlinear restoring forceTable 3
Nonlinear parameters of FE model and experimental model.
ParametersTrueEstimationError (%)Freeplay gap (°)1.701.731.76Linear-part stiffness (N·m/rad)2.902.910.34After the validations of the main wing and aileron, as well as the nonlinear system identification of freeplay, the construction of the nonlinear governing equation is finally in reach. By setting the rotational gap to zero using linear arm as in Figure5(a), the nonlinear system retrieved to a linear one. The result of the modal analysis shows a great agreement between the ROM and the test model as summarized in Table 4. The first four out-plane modes are shown in Figure 14.Table 4
Comparison of natural frequency between the ROM and GVTs.
Mode shapeGVTs (Hz)ROM (Hz)Error (%)1st bending mode1.931.892.072nd bending mode11.9511.692.18Aileron rotating mode12.5812.600.161st torsion mode19.5019.062.26Figure 14
First four mode shapes of the linear FE model.
(a)
1st bending mode(b)
2nd bending mode(c)
Aileron rotating mode(d)
1st torsion mode


## 5. Wind Tunnel Tests and Validation

Linear and nonlinear flutter experiments were achieved after the model mounted vertically in the wind tunnel with a zero angle of attack, which is presented in Figure15. The airspeed increased gradually until nearing the critical flutter velocity. At each airspeed, signals of all the transducers were measured when reaching stability. However, unexpected phenomena might occur during the increase of airspeed when the oscillation was significantly visible. Vibration signals were also measured and saved during the velocity change of incoming flow. All measured acceleration signals were filtered by a low-pass filter to decrease the influence of noise.Figure 15
Experimental setup of linear and nonlinear wind tunnel flutter tests.
(a)
Layout(b)
ModelA photograph of accelerometer arrangement is shown in Figure16. The rotational angle of the aileron was measured by a rotational displacement sensor, as shown in Figure 4. Three PCB accelerometers were mounted to the aileron and to the tip span of main wing as declared in Table 5. Z-axis is perpendicular to the plane of rectangular wing. These measured signals are directly recorded on a computer with data acquisition and analysis software National Instrument LabVIEW.Figure 16
Locations of accelerometers in the wind tunnel flutter test.Table 5
Definition and location of measure sensors during the wind tunnel flutter tests.
Measured quantityDirectionSensor locationAccelerationZ-axisThe trailing edge of aileronAccelerationZ-axisRoot of aileronAccelerationZ-axisTip span of the main wingAngular displacement—The rotational axis of aileron


### 5.1. Linear Flutter Analysis and Test

Firstly, the linear flutter analysis was performed byg-method [45] using aeroelastic analysis software ZONA ZAERO [46] before the wind tunnel flutter tests. The unsteady aerodynamic force was calculated by the subsonic DLM [47] to obtain Q in Equation (27). The aerodynamic computational grid was divided into two components as shown in Figure 17(a). The aileron rotational stiffness is linear with a quantity of 2.91 N·m/rad as mentioned beforehand. The numerical V-g and V-f flutter diagrams are shown in Figures 17(b) and 17(c). The aileron rotating mode crosses the zero damping line at airspeed 47.6 m/s with an oscillation frequency of 13.4 Hz. The flutter mode was the interaction of aileron rotating mode and the 1st torsion mode.Figure 17
The results of linear flutter analysis.
(a)
Aerodynamic computational grid(b)
V-g plot(c)
V-f plotThe linear flutter test was carried out to investigate the characteristics of the rectangular wing. In the test process, the airspeed increased gradually up to 45.5 m/s starting from 30 m/s. The vibrating responses enhanced consistently with the increasing incoming flow. Approximately, the experimental flutter airspeed was around 45.5 m/s. Fourier transform was performed and the main oscillation frequency is 13.4 Hz, as shown in Figure18. To avoid model suffering catastrophic damage, however, the experimental critical flutter velocity is modestly lower than the linear computational values.Figure 18
FFT of the acceleration responses at the trailing edge of the aileron with the airspeed of 45.5 m/s.


### 5.2. Nonlinear LCO Analysis and Test

By substituting the nonlinear parameters estimated by the TNSI method in the process of the extended CMS method, the time-marching approach described in the previous section was performed. The simulations of the nonlinear ROM were performed by the Runge-Kutta scheme with a variable time step of which the upper limit was 1 ms. The simulated time interval should be sufficiently long enough until the response reaching stable equilibrium. A proper initial condition of the nonlinear system was essential in the numerical simulation by the time-marching approach. However, the experimental model did not require any external excitation since any small perturbations on model during wind tunnel testing would evoke LCO [18].The simulation airspeed range is between 10.0 m/s and 50.0 m/s. The LCO begins at 31.6 m/s as shown in Figures19(a) and 19(b). By using the Fourier transform, the fundamental oscillating frequency was 15.5 Hz. At 40.0 m/s, the time history of aileron rotation angle is depicted in Figure 19(c). The oscillation amplitude gradually grew as the airspeed increased. From the Fourier analysis in Figure 19(d), the fundamental oscillating frequency was 15.3 Hz, with the second primary oscillating frequency of 45.6 Hz. The higher-order harmonic component was caused by the existence of freeplay in the aileron rotational hinge. Continue to increase the airspeed to 47.0 m/s, under which the LCO approached the critical boundary as depicted in Figure 19(e). The fundamental oscillating frequency was 14.2 Hz, as presented in Figure 19(f) by using a Fourier transform. The system is diverging at 47.5 m/s with a divergence frequency of 13.9 Hz as shown in Figures 19(g) and 19(h).Figure 19
The angular displacement response of aileron rotational DOF (the freeplay interval in time history plot is depicted by two red dashed lines).
(a)
Time history at 31.6 m/s(b)
Portraits at 31.6 m/s(c)
Time history at 40.0 m/s(d)
Portraits at 40.0 m/s(e)
Time history at 47.0 m/s(f)
Portraits at 47.0 m/s(g)
Time history at 47.5 m/s(h)
Portraits at 47.5 m/sIn total, the portraits of aileron rotational DOF at different airspeeds are shown in Figure20 at 32.0 m/s, 40.0 m/s, 44.0 m/s, 46.0 m/s, and 47.0 m/s, respectively. With the airspeed increasing to 35.0 m/s, the portraits expand and distort apparently because of the higher-order oscillation. However, the shape of the portraits does not change abruptly, which indicates that the form of LCO does not change with airspeed.Figure 20
Portraits of aileron rotational DOF at different airspeeds.The nonlinear wind tunnel flutter tests were performed to investigate the characteristics of LCO induced by aileron rotational freeplay and to validate the integrated test-based nonlinear aeroelastic analysis methodology. The acceleration responses at the trailing edge of aileron were used to characterize the LCO portraits, and first-order integration was performed to acquire the velocity. The nonlinear test setup was similar to the linear one as presented in Figure15 with 0° angle of attack. The airspeed increased gradually and the LCO was observed during the wind tunnel test. Unfortunately, the square-section elastic beam, which offers the aileron rotational stiffness, suffered plastic deformation because of the long-time vibration of LCO, and the flutter test ended at 40 m/s to avoid further model damage, which was far from the linear flutter boundary.The comparison between the results of simulation and experiments was performed on accelerometer 1 at the trailing edge of aileron and accelerometer 3 at the tip span of the main wing. The aerodynamic noise and disturbances increased with the increasing air speed during the wind tunnel tests, which were reflected on the sensors. However, the effect of the noise was not completely considered during the simulation. Moreover, the abrupt increase in amplitude can be regarded as the starting point of the LCO, based on the freeplay-induced LCO characteristics [8]. The results of nonlinear LCO test revealed one particular oscillating pattern of LCO happened and grew gradually with the increase of airspeed as shown in Figure 21. The system enters nonlinear LCO at 33.5 m/s with the acceleration amplitude escalating sharply, on the contrast to the linear flutter test.Figure 21
Results of the nonlinear LCO test compared to the numerical analysis.
(a)
Accelerometer 1(b)
Accelerometer 3Figure22 gives the results of the nonlinear LCO test at 36.0 m/s and 40.0 m/s. The experimental response in Figures 22(a) and 22(b) and the portraits in Figures 22(c) and 22(d) indicate a stable and significant nonlinear LCO phenomena. However, the shape of portraits is slightly shifted because of the higher-order oscillations and the complicated environment noise during the wind tunnel test. The time-frequency analysis in Figures 22(e) and 22(f) and FFT analysis in Figures 22(g) and 22(h) are utilized to the experimental response, informing the fundamental oscillation frequency of 15.1 Hz, and other two peaks in frequency-domain are at 30.0 Hz and 45.2 Hz, which near the integer multiples of fundamental frequency. Even though these higher-order harmonic components depicted by the distortion of LCO portraits existed, the magnitude of which was considerably small compared to the fundamental frequency.Figure 22
Nonlinear LCO test results at the trailing edge of aileron (accelerometer 1).
(a)
Time history at 36.0 m/s(b)
Time history at 40.0 m/s(c)
Portrait at 36.0 m/s(d)
Portrait at 40.0 m/s(e)
Time-frequency analysis at 36.0 m/s(f)
Time-frequency analysis at 40.0 m/s(g)
FFT at 36.0 m/s(h)
FFT at 40.0 m/sClearly, the integrated test-based method described in this paper could predict the appearance of LCO and track the development of limit-cycle motion. The error between the analysis and test is under 5% as in Table6, confirming that the contribution of present work paves an alternative way to analyze nonlinear aeroelastic phenomena.Table 6
Comparison on amplitude of the numerical analysis and wind tunnel tests.
ItemNumerical analysisWind tunnel testsError (%)Linear flutter boundarySpeed47.5 m/s45.5 m/s4.41Frequency13.9 Hz13.4 Hz3.71Nonlinear LCO starting boundarySpeed31.6 m/s33.0 m/s4.24Frequency15.5 Hz15.14 Hz2.38


## 5.1. Linear Flutter Analysis and Test

Firstly, the linear flutter analysis was performed byg-method [45] using aeroelastic analysis software ZONA ZAERO [46] before the wind tunnel flutter tests. The unsteady aerodynamic force was calculated by the subsonic DLM [47] to obtain Q in Equation (27). The aerodynamic computational grid was divided into two components as shown in Figure 17(a). The aileron rotational stiffness is linear with a quantity of 2.91 N·m/rad as mentioned beforehand. The numerical V-g and V-f flutter diagrams are shown in Figures 17(b) and 17(c). The aileron rotating mode crosses the zero damping line at airspeed 47.6 m/s with an oscillation frequency of 13.4 Hz. The flutter mode was the interaction of aileron rotating mode and the 1st torsion mode.Figure 17
The results of linear flutter analysis.
(a)
Aerodynamic computational grid(b)
V-g plot(c)
V-f plotThe linear flutter test was carried out to investigate the characteristics of the rectangular wing. In the test process, the airspeed increased gradually up to 45.5 m/s starting from 30 m/s. The vibrating responses enhanced consistently with the increasing incoming flow. Approximately, the experimental flutter airspeed was around 45.5 m/s. Fourier transform was performed and the main oscillation frequency is 13.4 Hz, as shown in Figure18. To avoid model suffering catastrophic damage, however, the experimental critical flutter velocity is modestly lower than the linear computational values.Figure 18
FFT of the acceleration responses at the trailing edge of the aileron with the airspeed of 45.5 m/s.


## 5.2. Nonlinear LCO Analysis and Test

By substituting the nonlinear parameters estimated by the TNSI method in the process of the extended CMS method, the time-marching approach described in the previous section was performed. The simulations of the nonlinear ROM were performed by the Runge-Kutta scheme with a variable time step of which the upper limit was 1 ms. The simulated time interval should be sufficiently long enough until the response reaching stable equilibrium. A proper initial condition of the nonlinear system was essential in the numerical simulation by the time-marching approach. However, the experimental model did not require any external excitation since any small perturbations on model during wind tunnel testing would evoke LCO [18].The simulation airspeed range is between 10.0 m/s and 50.0 m/s. The LCO begins at 31.6 m/s as shown in Figures19(a) and 19(b). By using the Fourier transform, the fundamental oscillating frequency was 15.5 Hz. At 40.0 m/s, the time history of aileron rotation angle is depicted in Figure 19(c). The oscillation amplitude gradually grew as the airspeed increased. From the Fourier analysis in Figure 19(d), the fundamental oscillating frequency was 15.3 Hz, with the second primary oscillating frequency of 45.6 Hz. The higher-order harmonic component was caused by the existence of freeplay in the aileron rotational hinge. Continue to increase the airspeed to 47.0 m/s, under which the LCO approached the critical boundary as depicted in Figure 19(e). The fundamental oscillating frequency was 14.2 Hz, as presented in Figure 19(f) by using a Fourier transform. The system is diverging at 47.5 m/s with a divergence frequency of 13.9 Hz as shown in Figures 19(g) and 19(h).Figure 19
The angular displacement response of aileron rotational DOF (the freeplay interval in time history plot is depicted by two red dashed lines).
(a)
Time history at 31.6 m/s(b)
Portraits at 31.6 m/s(c)
Time history at 40.0 m/s(d)
Portraits at 40.0 m/s(e)
Time history at 47.0 m/s(f)
Portraits at 47.0 m/s(g)
Time history at 47.5 m/s(h)
Portraits at 47.5 m/sIn total, the portraits of aileron rotational DOF at different airspeeds are shown in Figure20 at 32.0 m/s, 40.0 m/s, 44.0 m/s, 46.0 m/s, and 47.0 m/s, respectively. With the airspeed increasing to 35.0 m/s, the portraits expand and distort apparently because of the higher-order oscillation. However, the shape of the portraits does not change abruptly, which indicates that the form of LCO does not change with airspeed.Figure 20
Portraits of aileron rotational DOF at different airspeeds.The nonlinear wind tunnel flutter tests were performed to investigate the characteristics of LCO induced by aileron rotational freeplay and to validate the integrated test-based nonlinear aeroelastic analysis methodology. The acceleration responses at the trailing edge of aileron were used to characterize the LCO portraits, and first-order integration was performed to acquire the velocity. The nonlinear test setup was similar to the linear one as presented in Figure15 with 0° angle of attack. The airspeed increased gradually and the LCO was observed during the wind tunnel test. Unfortunately, the square-section elastic beam, which offers the aileron rotational stiffness, suffered plastic deformation because of the long-time vibration of LCO, and the flutter test ended at 40 m/s to avoid further model damage, which was far from the linear flutter boundary.The comparison between the results of simulation and experiments was performed on accelerometer 1 at the trailing edge of aileron and accelerometer 3 at the tip span of the main wing. The aerodynamic noise and disturbances increased with the increasing air speed during the wind tunnel tests, which were reflected on the sensors. However, the effect of the noise was not completely considered during the simulation. Moreover, the abrupt increase in amplitude can be regarded as the starting point of the LCO, based on the freeplay-induced LCO characteristics [8]. The results of nonlinear LCO test revealed one particular oscillating pattern of LCO happened and grew gradually with the increase of airspeed as shown in Figure 21. The system enters nonlinear LCO at 33.5 m/s with the acceleration amplitude escalating sharply, on the contrast to the linear flutter test.Figure 21
Results of the nonlinear LCO test compared to the numerical analysis.
(a)
Accelerometer 1(b)
Accelerometer 3Figure22 gives the results of the nonlinear LCO test at 36.0 m/s and 40.0 m/s. The experimental response in Figures 22(a) and 22(b) and the portraits in Figures 22(c) and 22(d) indicate a stable and significant nonlinear LCO phenomena. However, the shape of portraits is slightly shifted because of the higher-order oscillations and the complicated environment noise during the wind tunnel test. The time-frequency analysis in Figures 22(e) and 22(f) and FFT analysis in Figures 22(g) and 22(h) are utilized to the experimental response, informing the fundamental oscillation frequency of 15.1 Hz, and other two peaks in frequency-domain are at 30.0 Hz and 45.2 Hz, which near the integer multiples of fundamental frequency. Even though these higher-order harmonic components depicted by the distortion of LCO portraits existed, the magnitude of which was considerably small compared to the fundamental frequency.Figure 22
Nonlinear LCO test results at the trailing edge of aileron (accelerometer 1).
(a)
Time history at 36.0 m/s(b)
Time history at 40.0 m/s(c)
Portrait at 36.0 m/s(d)
Portrait at 40.0 m/s(e)
Time-frequency analysis at 36.0 m/s(f)
Time-frequency analysis at 40.0 m/s(g)
FFT at 36.0 m/s(h)
FFT at 40.0 m/sClearly, the integrated test-based method described in this paper could predict the appearance of LCO and track the development of limit-cycle motion. The error between the analysis and test is under 5% as in Table6, confirming that the contribution of present work paves an alternative way to analyze nonlinear aeroelastic phenomena.Table 6
Comparison on amplitude of the numerical analysis and wind tunnel tests.
ItemNumerical analysisWind tunnel testsError (%)Linear flutter boundarySpeed47.5 m/s45.5 m/s4.41Frequency13.9 Hz13.4 Hz3.71Nonlinear LCO starting boundarySpeed31.6 m/s33.0 m/s4.24Frequency15.5 Hz15.14 Hz2.38


## 6. Conclusion

An integrated test-based nonlinear aeroelastic analysis method for freeplay-induced LCO is demonstrated in a real-world case in the present work. The nonlinear model construction of this approach combines linear modeling techniques and nonlinear system identification. Then, both the numerical and experimental investigations on a wing-aileron system with aileron rotational freeplay were performed. Compared to the experimental tests, the prediction of nonlinear LCO induced by the aileron rotational freeplay is feasible to an acceptable degree of engineering accuracy.A low-speed nonlinear flutter experimental model was designed, constituted by metal structures and 3D printing wing shells. We dexterously designed the freeplay rotational arm that could introduce a small gap at the aileron rotational freeplay, while the main wing and the aileron remain to be linear. The wind tunnel tests indicated that the experimental model is capable to obtain the nonlinear limit-cycle oscillations induced by freeplay.A limitation of the proposed method is that it does not apply to all types of structural nonlinearities in an aeroelastic system. Only concentrated parametric nonlinearities, such as cubic stiffness and freeplay, can use this method. For some kind of high aspect ratio wing with geometric parametric nonlinearities, another alternative approach will be needed. In the future, we plan to (1) focus on controlling and suppressing the LCO to increase the nonlinear flutter speed as we mounted an actuator inside the model, (2) design a experimental wing with two control surfaces, considering freeplay in either or both hinges of the control surface, and (3) design the friction generating device in the wing-aileron hinge and study the effect of friction on the LCO.



---

*Source: 5195030-2022-09-26.xml*