Whirl Flutter Suppression of Tiltrotor Aircraft Using Actively Controlled Aileron

Abstract

Whirl flutter of a tiltrotor aircraft is a complex aeroelastic phenomenon and it can result in catastrophic consequences. The deflection of an aileron mounted on a wing has the potential to solve this fatal problem. Whirl flutter suppression using an actively controlled aileron is studied in this study. Firstly, a semi-span aeroelastic model is established for the whirl flutter problem using the Hamilton principle. This model is composed of three parts: a rigid rotor, a rigid nacelle and a flexible wing, and the effect of the aileron deflection on the aeroelastic responses is also taken into consideration through a quasi-steady aerodynamic model. In addition, the accuracy of this aeroelastic model is validated with the results of two different wind-tunnel tests. Then, an LQR controller is developed to control the dynamic deflection of the aileron, and a full-dimensional state observer is built to estimate the state of the time-invariant system of a tiltrotor aircraft. Finally, simulations are carried out using the aeroelastic model and the LQR controller at different flight conditions to study the influence of the aileron deflection on whirl flutter. The simulation results demonstrate that the flutter boundary speed can be improved by 18.1% with the active deflection of the aileron, compared with the uncontrolled condition.

1. Introduction

Tiltrotor aircraft combine the advantages of helicopters and fixed-wing propeller-driven aircraft with hovering and high-speed forward flight. This type of aircraft, however, suffers from dynamic instability problems resulting from aeroelastic coupling between the rotor and flexible wing. In high-speed flight conditions, the in-plane aerodynamics resulting from blade flap movement couple with the flapping and torsional deformation of the flexible wing. This coupling causes the instability of whirl flutter and limit the high-speed advantage of tiltrotor aircraft.

In general, there are two different ways to suppress the whirl flutter of tiltrotor aircraft, namely active control and passive control. Typical passive control methods, such as blade geometry optimization and aeroelastic tailoring [1,2,3,4,5,6], enable favorable coupling between the rotor and the flexible wing to improve the stability boundary of whirl flutter. In recent years, new concepts of passive control, including extensions and winglets, have received increasing attention [7,8,9], and their principle is to generate a restoring force which is opposite to the twist motion of the flexible wing, using these additional aerodynamic surfaces, thus suppressing whirl flutter.

Active control has shown its peculiar advantages in vibration alleviation and aeroelastic stability augmentation, and its application in the whirl flutter suppression of tiltrotor aircraft dates back to the 1970s. Depending on the technical implementation, active control methods for whirl flutter suppression can be generally divided into two categories, including active control through the swashplate of the rotor and active control through ailerons mounted on a flexible wing. In 1973, Boeing added feedback control to the Boeing Model 222 rotor [10], and the modal damping of the wing vertical bending mode, which was critical for the whirl flutter phenomenon, was increased through the active control of the collective and cyclic pitch angles of the rotors. In 1986, Nasu et al. at the NASA Ames Research Center used APC (Active Pitch Control) to control the longitudinal cyclic control of the rotor and selected the chord bending velocity as the feedback signal [11]. It was shown that this method was capable of improving the boundary velocity of whirl flutter. In 2004, Mueller et al. of Eurocopter conducted an aeroelastic study by combining multibody dynamics and a GPC (Generalized Predictive Control) algorithm, in which the wing bending deformation was used as the feedback signal to control the input of the swashplate, and the simulation results demonstrated that the GPC control was able to extend the whirl flutter boundary by 21% [12]. In 2017, Matthew et al. carried out a simulation study of the influence of active wing tips on whirl flutter stability [13]. The deflection of the wing tips, which was controlled by a PID controller, was capable of increasing the boundary velocity from 160 kts to over 200 kts, demonstrating that a movable aerodynamic surface had the potential to improve the boundary velocity of whirl flutter. Paik et al. examined the effectiveness of two different methods of active stability augmentation, including wing-flaperon and swashplate actuation, and feedback control based on the wing states was used instead of full-state control for simplicity [14]. A new wind-tunnel test model named Tiltrotor Aeroelastic Stability Testbed (TRAST) was developed by the US Army and NASA, and an active stability augmentation test was conducted with this testbed [15,16]. Ivanco et al. performed a simulation of active stability augmentation and vibration reduction on a tiltrotor model; in this study, an aeroelastic model was established for the TRAST model using the Rotorcraft Comprehensive Analysis System (RCAS), and a GPC algorithm was used for the blade pitch [17].

The aforementioned studies have proven the effectiveness of active stability augmentation methods. In order to further improve the performance of active whirl flutter suppression, an active control method based on an LQR algorithm and full-order state observer was built and tested through simulation. A comprehensive aeroelastic model composed of a rigid rotor, a flexible wing and a rigid nacelle was established. Aerodynamics generated by the aileron were calculated using a quasi-steady model. The experiment results of two different wind-tunnel tests were used to validate the accuracy of this aeroelastic model. An active controller was designed based on an LQR (Linear Quadratic Regulator) algorithm to control the deflection movement of the aileron. A full-order state observer was built to obtain the unmeasurable variables of the flexible wing, which were used as the feedback signals in the active control. The simulations were conducted at different flight conditions to suppress whirl flutter using the established aeroelastic model combined with the LQR controller. Based on the simulation results, a semi-span model equipped with an actively controlled aileron will be designed and tested in a wind tunnel for further validation of the active stability augmentation method.

2. Aeroelastic Model

2.1. Definition of Coordinate Systems

According to numerical research results, the symmetric modal coupling between the rotor and flexible wing of a tiltrotor aircraft is more prone to whirl flutter than the antisymmetric modal coupling [18], making a semi-span model sufficient for whirl flutter simulation. In addition, this type of model can dramatically reduce the computational cost. Therefore, a semi-span aeroelastic model was built in this study.

A series of coordinate systems were established to describe the motion and deformation of the semi-span model. The flexible wing and the coordinate systems defined on the rotor are shown in Figure 1.

As shown in Figure 1, the origin of the inertial coordinate system ${\left\{\begin{array}{ccc}{I}_i& J_i& K_i\end{array}\right\}}^T$ is defined at the root of the wing elastic axis. The origin of the undeformed coordinate system ${\left\{\begin{array}{ccc}{I}_d& J_d& K_d\end{array}\right\}}^T$ is located at the wing elastic axis. Three twist angles are generated between the deformed coordinate system ${\left\{\begin{array}{ccc}{I}_w& J_w& K_w\end{array}\right\}}^T$ and the undeformed coordinate system ${\left\{\begin{array}{ccc}{I}_d& J_d& K_d\end{array}\right\}}^T$ of the flexible wing due to its elastic deformation. The nacelle coordinate system ${\left\{\begin{array}{ccc}{I}_p& J_p& K_p\end{array}\right\}}^T$ is used to describe the movement of the rigid nacelle, with the origin defined on the tilt hinge, and ${\left\{\begin{array}{ccc}{x}_{pw}& y_{pw}& z_{pw}\end{array}\right\}}^T$ are the biases of the origin with respect to the wing elastic axis. ${\alpha }_p$ is the tilt angle of the nacelle, with the angle corresponding to the helicopter mode defined as 0°. ${\left\{\begin{array}{ccc}{I}_h& J_h& K_h\end{array}\right\}}^T$ is the rotor hub coordinate system, with its origin defined at the center of the hub, and ${\left\{\begin{array}{ccc}{x}_h& y_h& z_h\end{array}\right\}}^T$ are the biases of the origin with respect to the tilt hinge. The rotor rotating coordinate system ${\left\{\begin{array}{ccc}{I}_r& J_r& K_r\end{array}\right\}}^T$ is generated by rotating the rotor hub coordinate system an azimuth of $\Psi \mathrm{around}\mathrm{the}{K}_h\mathrm{axis}$. The blade rigid flap coordinate system ${\left\{\begin{array}{ccc}{I}_g& J_g& K_g\end{array}\right\}}^T$ is established by rotating the ${\left\{\begin{array}{ccc}{I}_r& J_r& K_r\end{array}\right\}}^T$ an angle of ${\beta }_G^i+{\beta }_p$ around the $I_r$ axis, with ${\beta }_G^i$ being the flap angle of the $i\mathrm{th}$ blade, and ${\beta }_p$ being the precone angle of the rotor hub. The blade pitch coordinate system ${\left\{\begin{array}{ccc}{I}_b& J_b& K_b\end{array}\right\}}^T$, which is used to describe the rigid flapping and pitching motion of the rotor blades, is attached to the pitching hinge.

2.2. Hamilton Principle

The aeroelastic dynamic equations of the semi-span model are derived using the Hamilton principle [19,20,21,22] as expressed in Equation (1):

$$\mathsf{\delta }\mathsf{\Pi }=\underset{t_1}{\overset{t_2}{{\displaystyle \int }}}\left(\delta U-\delta T-\delta W\right)dt=0$$

(1)

In Equation (1), $\delta U$ is the virtual elastic potential energy of the model, including the elastic wing and the rotor hub, given by:

$$\delta U=\mathsf{\delta }{U}_{wing}+\delta U_{hub}$$

(2)

$\delta T$ is the virtual kinetic energy, which combines the virtual kinetic energy of the elastic wing, the nacelle and the rotor blades, and $\delta T$ is calculated by:

$$\delta T=\mathsf{\delta }{T}_{wing}+\delta T_{nacelle}+{\displaystyle \sum }_{i=1}^{N_b}\delta T_{b_i}$$

(3)

$\delta W$ is the virtual work of the aerodynamics. In this study, the aerodynamics of the wing and the rotor blades are taken into consideration. $\delta W$ is estimated as:

$$\delta W=\mathsf{\delta }{W}_{wing}+{\displaystyle \sum }_{i=1}^{N_b}\delta W_{b_i}$$

(4)

2.3. Flexible Wing Model

The flexible wing is one of the most fundamental parts of the semi-span model, and its elastic deformation is taken into consideration, while the movement of the fuselage is ignored. Therefore, the wing can be modeled as a cantilever beam as shown in Figure 2.

The motion modes of the flexible wing include the flap bending mode $q_1$, chordwise bending mode $q_2$ and torsional mode $p$, and mode shape functions are used to describe the elastic deformation of the wing of these three modes. Since the whirl flutter boundary of the tiltrotor aircraft is primarily determined by the lower modes of the flexible wing, only the lowest order of the motion modes is retained, so that the elastic deformation of the wing can be described as:

$$z_w\left(y_w,\mathrm{t}\right)=q_1\left(t\right){\eta }_1\left(y_w\right)$$

(5)

$$x_w\left(y_w,\mathrm{t}\right)=q_2\left(t\right){\eta }_2\left(y_w\right)$$

(6)

$${\theta }_w\left(y_w,\text{t}\right)=p\left(t\right)\epsilon \left(y_w\right)$$

(7)

where $z_w$, $x_w$ and ${\theta }_w$ are the flap bending displacement, the chordwise bending displacement and the twist angle of the wing, respectively.

After defining the deformation of the flexible wing, its elastic potential energy can be calculated by:

$$U_{wing}=\frac{1}{2}\underset{0}{\overset{y_{TW}}{{\displaystyle \int }}}\left[E{I}_X{\left(\frac{{\partial }^2{z}_w}{{\partial }^2{y}_w}\right)}^2+E{I}_Z{\left(\frac{{\partial }^2{x}_w}{{\partial }^2{y}_w}\right)}^2+GJ{\left(\frac{{\partial }^2{\theta }_w}{{\partial }^2{y}_w}\right)}^2\right]d{y}_w$$

(8)

Then, the stiffness matrix related to the elastic potential energy can be determined using:

$$K_{w_{ij}}=\frac{{\partial }^2\delta U_{wing}}{\partial \delta q_i\delta q_j}$$

(9)

The kinetic energy of the wing is given as:

$$T_{wing}=\frac{1}{2}\underset{0}{\overset{y_{TW}}{{\displaystyle \int }}}\left[m_z{\left(\frac{\partial z_w}{\partial t}\right)}^2+m_x{\left(\frac{\partial x_w}{\partial t}\right)}^2+I_{\theta w}{\left(\frac{\partial {\theta }_w}{\partial t}\right)}^2\right]d{y}_w$$

(10)

The virtual kinetic energy of the wing can be described as generalized forces as:

$$\delta T_{wing}={\displaystyle \sum }_1^i\frac{\partial \left(\delta T_{wing}\right)}{\partial q_i}\delta q_i$$

(11)

The mass matrix related to the kinetic energy of the wing can be evaluated by:

$$M_{w_{ij}}=\frac{{\partial }^2\left(-\delta T_{wing}\right)}{\partial \delta q_i\delta {\ddot{q}}_j}$$

(12)

The work of the aerodynamics acting on the wing is expressed as:

$$W_{wing}=\underset{0}{\overset{y_{TW}}{{\displaystyle \int }}}\left[F_Z·z_w+F_X·x_w+\left(M_P·{\theta }_w-F_Z·\left({\overline{X}}_A·c_w·{\theta }_w\right)\right)\right]d{y}_w$$

(13)

In Equation (13), ${\overline{X}}_A$ is the non-dimensioned bias of the aerodynamic center with respect to the elastic axis. The virtual work of the wing related to the aerodynamics can be described as generalized forces as:

$$\delta W_{wing}={\displaystyle \sum }_1^i\frac{\partial \left(\delta W_{wing}\right)}{\partial q_i}\delta q_i$$

(14)

The damping matrix and stiffness matrix related to the aerodynamic work are as follows:

$$C_{w_{ij}}=\frac{{\partial }^2\left(-\delta W_{wing}\right)}{\partial \delta q_i\delta {\dot{q}}_j}$$

(15)

$$K_{w_{ij}}=\frac{{\partial }^2\left(-\delta W_{wing}\right)}{\partial \delta q_i\delta q_j}$$

(16)

The structural damping of the flexible wing can be expressed as:

$$C=2\zeta \sqrt{M_w{K}_w}$$

(17)

where $C$ is the modal damping, $M_w$ is the modal mass, $K_w$ is the modal stiffness and $\zeta$ is the modal damping ratio.

2.4. Nacelle Model

The nacelle is mainly composed of several large components, such as the engine and rotor shaft. The elastic deformation of these components is not strongly related to whirl flutter; hence, the nacelle can be described by a rigid model. Any point on the nacelle can be expressed as:

$$R_{pi}=\left\{\begin{array}{c}{u}_{tip}\\ v_{tip}\\ w_{tip}\end{array}\right\}+{\left(T_{dw}{T}_{wi}\right)}^T\left\{\begin{array}{c}{x}_{pw}\\ y_{pw}\\ z_{pw}\end{array}\right\}+{\left(T_p{T}_{pw}{T}_{wd}{T}_{di}\right)}^T\left\{\begin{array}{c}{x}_p\\ y_p\\ z_p\end{array}\right\}$$

(18)

where ${\left\{\begin{array}{ccc}{x}_{pw}& y_{pw}& z_{pw}\end{array}\right\}}^T$ is the bias of the nacelle connection point with respect to the elastic axis of the wing, ${\left\{\begin{array}{ccc}{x}_p& y_p& z_p\end{array}\right\}}^T$ is the coordinate of any point located in the nacelle coordinate system and ${\left\{\begin{array}{ccc}{u}_{tip}& v_{tip}& w_{tip}\end{array}\right\}}^T$ is the displacement vector of the wing tip. The velocity of any point related to the nacelle can be calculated by:

$$V_{pi}=\frac{d{R}_{pi}}{dt}$$

(19)

Then, yields the virtual kinetic energy of the nacelle:

$$\delta T_{nacelle}={{\displaystyle \int \int \int }}^{}{\rho }_p\delta V_{pi}^T{V}_{pi}dV$$

(20)

The degrees of freedom (dof) of the nacelle are identical to those of the wing, and $q={\left\{\begin{array}{ccc}{q}_1& q_2& p\end{array}\right\}}^T$. Hence, the linearized mass matrix of the nacelle can be expressed as:

$$M_{P_{ij}}=\frac{{\partial }^2\left(-\delta T_{nacelle}\right)}{\partial \delta q_i\delta {\ddot{q}}_j}$$

(21)

2.5. Rotor Model

In this aeroelastic model, the flap bending stiffness of the rotor hub with a universal joint is taken into consideration, while each rotor blade is treated as a rigid body. Then, the virtual elastic potential energy of the ith rotor blade resulting from the flap motion can be expressed as:

$$\delta U_i=\delta {\beta }_G^i\left(t\right)k_{\beta G}{\beta }_G^i\left(t\right)$$

(22)

where $k_{\beta G}$ is the flap bending stiffness of the rotor hub. The flapping angle of the ith rotor blade ${\beta }_G^i\left(t\right)$, which is related to the rigid flapping dof of the rotor hub, can be expressed in terms of the flapping dof of the rotor hub and the azimuth ${\Psi }_i\left(t\right)$ of the ith rotor blade:

$${\beta }_G^i\left(t\right)=\left[cos{\Psi }_i\left(t\right)-sin{\Psi }_i\left(t\right)\right]\left[\begin{array}{c}{\beta }_{GC}\left(t\right)\\ {\beta }_{GS}\left(t\right)\end{array}\right]$$

(23)

Adding the virtual elastic potential energy of each blade together yields the virtual elastic potential energy of the rotor:

$$\delta U_{hub}={\displaystyle \sum }_{i=1}^{N_b}\delta U_i=\left[\begin{array}{cc}\delta {\beta }_{GC}\left(t\right)& \delta {\beta }_{GS}\left(t\right)\end{array}\right]\left[\begin{array}{cc}\frac{N_b{k}_{\beta G}}{2}& 0\\ 0& \frac{N_b{k}_{\beta G}}{2}\end{array}\right]\left[\begin{array}{c}{\beta }_{GC}\left(t\right)\\ {\beta }_{GS}\left(t\right)\end{array}\right]$$

(24)

The dof of the rigid rotor blade, which combines the convected motion resulting from the wing elastic deformation and the blade’s rotating dof ${\beta }_G^i\left(t\right)$, is expressed as $q={\left\{\begin{array}{cccc}{q}_1& q_2& p& {\beta }_G^i\end{array}\right\}}^T$. The position vector of any point in the blade coordinate system can be expressed in the inertial coordinate system as:

$$\begin{gathered} \begin{array}{c}{R}_{Ri}=\left\{\begin{array}{c}{u}_{tip}\\ v_{tip}\\ w_{tip}\end{array}\right\}+{\left(T_{dw}{T}_{wi}\right)}^T\left\{\begin{array}{c}{x}_{pw}\\ y_{pw}\\ z_{pw}\end{array}\right\}+{\left(T_x{T}_p{T}_{wi}^T{T}_{dw}{T}_{wi}\right)}^T\left\{\begin{array}{c}{x}_h\\ y_h\\ z_h\end{array}\right\}\\ \end{array} \\ +{\left(T_{bg}{T}_{gr}{T}_{rp}{T}_{hp}{T}_p{T}_{pw}{T}_{wd}{T}_{di}\right)}^T\left\{\begin{array}{c}{x}_b\\ y_b\\ z_b\end{array}\right\} \end{gathered}$$

(25)

The velocity of this point is calculated by:

$$V_{Ri}=\frac{d{R}_{Ri}}{dt}$$

(26)

Then, the virtual kinetic energy of each blade is given by:

$$\delta T_{blade}={{\displaystyle \int \int \int }}^{}{\rho }_b\delta V_{Ri}^T{V}_{Ri}dV$$

(27)

Using the zero-time end condition and exchanging the order of integration, we can obtain:

$$\delta T_{blade}={{\displaystyle \int \int \int }}^{}{\rho }_b\delta R_{Ri}^T{\ddot{R}}_{Ri}dV$$

(28)

The virtual kinetic energy of each blade can be described as generalized forces as:

$$\delta T_{blade}={\displaystyle \sum }_1^i\frac{\partial \left(\delta T_{blade}\right)}{\partial q_i}\delta q_i$$

(29)

The linearized mass matrix of the rotating blade can be represented as:

$$M_{b_{ij}}=\frac{{\partial }^2\left(-\delta T_{blade}\right)}{\partial \delta q_i\delta {\ddot{q}}_j}$$

(30)

Then, the stiffness matrix and damping matrix related to the blade virtual kinetic energy are given by:

$$K_{b_{ij}}=\frac{{\partial }^2\left(-\delta T_{blade}\right)}{\partial \delta q_i\delta q_j}$$

(31)

$$C_{b_{ij}}=\frac{{\partial }^2\left(-\delta T_{blade}\right)}{\partial \delta q_i\delta {\dot{q}}_j}$$

(32)

The transformation from the rotating frame to the non-rotating frame is expressed by:

$$\left\{\begin{array}{c}{q}_1\left(t\right)\\ q_2\left(t\right)\\ p\left(t\right)\\ {\beta }_G^i\left(t\right)\end{array}\right\}=\left[\begin{array}{ccccc}1& 0& 0& 0& 0\\ 0& 1& 0& 0& 0\\ 0& 0& 1& 0& 0\\ 0& 0& 0& cos{\Psi }_i\left(t\right)& -sin{\Psi }_i\left(t\right)\end{array}\right]\left\{\begin{array}{c}{q}_1\left(t\right)\\ q_2\left(t\right)\\ p\left(t\right)\\ {\beta }_{GC}\left(t\right)\\ {\beta }_{GS}\left(t\right)\end{array}\right\}$$

(33)

The sectional aerodynamic loads acting on the rotor blade are given by:

$$F_L=\frac{1}{2}\rho \left(V_T^2+V_P^2\right)c{c}_l$$

(34)

$$F_D=\frac{1}{2}\rho \left(V_T^2+V_P^2\right)c{c}_d$$

(35)

Decomposing these loads into components perpendicular and parallel to the rotor disk, we will obtain:

$$\begin{array}{c}{F}_x=0\\ F_y=-F_{L}sin\varphi -F_{D}cos\varphi \\ F_z=F_{L}cos\varphi -F_{D}sin\varphi \end{array}$$

(36)

The sectional aerodynamic loads of the rotor blade are denoted as:

$$F_{cb}={\left\{\begin{array}{ccc}{F}_x& F_y& F_z\end{array}\right\}}^T$$

(37)

They can be expressed in the inertial coordinate system as:

$$F_{ci}={\left(T_{gr}{T}_{rp}{T}_{hp}{T}_p{T}_{pw}{T}_{wd}{T}_{di}\right)}^T·F_{cb}$$

(38)

Then, the virtual aerodynamic work of each blade can be calculated by:

$$-\delta W_{blade}=-{{\displaystyle \int }}_0^R{\left(T_{gr}{T}_{rp}{T}_{hp}{T}_p{T}_{pw}{T}_{wd}{T}_{di}\right)}^T·F_{cb}·\delta R_{Ri}dr$$

(39)

The virtual work of each blade related to the aerodynamics can be described as generalized forces as:

$$\delta W_{blade}={\displaystyle \sum }_1^i\frac{\partial \left(\delta W_{blade}\right)}{\partial q_i}\delta q_i$$

(40)

Finally, the stiffness matrix and damping matrix related to the aerodynamic loads can be determined as:

$$K_{b_{ij}}=\frac{{\partial }^2\left(-\delta W_{blade}\right)}{\partial \delta q_i\delta q_j}$$

(41)

$$C_{b_{ij}}=\frac{{\partial }^2\left(-\delta W_{blade}\right)}{\partial \delta q_i\delta {\dot{q}}_j}$$

(42)

2.6. Aerodynamics of the Aileron

According to quasi-steady aerodynamic theory [23,24], the additional aerodynamic loads generated by the aileron, including lift $\Delta L$ and moment $\Delta M$, can be calculated using the pitch angle ${\delta }_f$, angular velocity ${\dot{\delta }}_f$ and angular acceleration ${\ddot{\delta }}_f$ of the aileron as follows:

$$\begin{array}{c}\Delta L=\frac{1}{2}\rho c\left[c_{l{\delta }_f}{\delta }_f{U}_t^2+c_{l{\dot{\delta }}_f}{\dot{\delta }}_f{U}_t^2+c_{l{\ddot{\delta }}_f}{\ddot{\delta }}_f\right]\\ \Delta M=\frac{1}{2}\rho c\left[c_{m{\delta }_f}{\delta }_f{U}_t^2+c_{m{\dot{\delta }}_f}{\dot{\delta }}_f{U}_t^2+c_{m{\ddot{\delta }}_f}{\ddot{\delta }}_f\right]\end{array}$$

(43)

These loads are denoted as:

$$F^a=\left[D_2\right]\ddot{\delta }+\left[D_1\right]\dot{\delta }+\left[D_0\right]\delta$$

(44)

In this study, the deflection motion of the aileron is selected as the control input. Therefore, the linearized governing equation of the semi-span model is given in the following form:

$$\left[M\right]\ddot{q}+\left[C\right]\dot{q}+\left[K\right]q=\left[D_2\right]\ddot{\delta }+\left[D_1\right]\dot{\delta }+\left[D_0\right]\delta$$

(45)

where $q={\left\{\begin{array}{ccccc}{q}_1& q_2& p& {\beta }_{GC}& {\beta }_{GS}\end{array}\right\}}^T$. $q_1$, $q_2$ and $q_2$ represent the flap bending mode, chordwise bending mode and twisting mode of the flexible wing, respectively. ${\beta }_{GC}$ and ${\beta }_{GS}$ are the flapping angles of the rotor hub described in the non-rotating coordinate system of the rotor hub.

2.7. Model Validation I

Firstly, the wind-tunnel test, which was conducted by Bell Helicopter with a full-scaled semi-span tiltrotor model [25], was chosen to validate the accuracy of the established model in this study. This test was performed at the 40-by-80-foot wind tunnel at the Ames Research Center, and the fundamental parameters of this test model are listed in

Table 1. Fundamental parameters of the Bell test model [25].

Parameters	Value
Number of blades $N_b$	3
Rotor radius $R$	3.81 m
Rotor solidity $\sigma$	0.089
Rotor blade pitch/flap coupling $K_p$	−0.268
Lift curve slope $a$	5.7
Rotor rotational speed in hover ${\Omega }_h$	565 r/min
Rotor rotational speed in cruise ${\Omega }_c$	458 r/min
Blade flap inertia $I_{\beta }$	142 kg/m2
Wing length (semi-span) $L_w$	1.333 R
Wing chord length $c$	0.413 R
Rotor mast height $h$	0.261 R
Pylon center of gravity location $h_p$	0.05 R
Pylon mass	655 kg
Pylon pitch moment of inertia $I_y$	257 kg/m2
Pylon yaw moment of inertia $I_x$	231 kg/m2
Wing sweep angle ${\delta }_w$	−6.5°

.

The modal damping and the modal frequencies of this semi-span model at different flight speeds were calculated by the aeroelastic model built in this study. As shown in Figure 3, the calculated values are in good agreement with the experiment results given in reference, indicating that the established aerodynamic model is of high accuracy.

2.8. Model Validation II

In order to study the aeroelastic characteristics of the tiltrotor aircraft and to validate the accuracy of the numerical methods, we performed a wind-tunnel test with a scaled semi-span model, as shown in Figure 4. The basic parameters of this model are listed in

Table 2. Fundamental parameters of the scaled semi-span model.

Parameters	Value
Number of blades $N_b$	3
Rotor blade chord length $c_b$	31 mm
Rotor rotational speed in cruise ${\Omega }_c$	1200–1800 r/min
Wing length (semi-span) $L_w$	690 mm
Wing chord length $c$	240 mm
Airfoil of wing	NACA0024
Number of wing sections $N_s$	7
Span length of wing section $L_s$	80 mm
Wing dihedral angle ${\delta }_d$	0°
Wing sweep angle ${\delta }_s$	0°
Wing preinstall angle ${\delta }_{pre}$	0°
Rotor mast height $h$	150 mm

, and more detail about the wind-tunnel test can be found in [26].

In the wind-tunnel tests, the nacelle vibrations were measured at different wind speeds using accelerometers, which were fixed on the nacelle. The vibration responses at wind speeds of 15 m/s and 20.5 m/s are shown in Figure 5. It can be seen that, at a wind speed of 20.5 m/s, the nacelle vibration increased dramatically, and whirl flutter phenomenon appeared.

Numerical simulations were carried out with the established aeroelastic model, in which a transient response analysis method was employed. The torsion responses at about 20 m/s are demonstrated in Figure 6. The wing torsion response converges slowly at a wind speed of 20 m/s; however, if the wind speed increases to 21 m/s, the response cannot converge, illustrating that the whirl flutter boundary of this semi-span model is around 21 m/s, which is in good agreement with the test results.

3. Active Controller and State Observer

3.1. Active Controller Based on LQR Algorithm

For the time-invariant system of the semi-span tiltrotor model, its state equation and output equation can be expressed as:

$$\begin{array}{c}\dot{x}\left(t\right)=Ax\left(t\right)+Bu\left(t\right)\\ y\left(t\right)=Cx\left(t\right)+Du\left(t\right)\end{array}$$

(46)

where $x={\left\{\begin{array}{cccccccccc}{\dot{q}}_1& {\dot{q}}_2& \dot{p}& {\dot{\beta }}_{GC}& {\dot{\beta }}_{GS}& q_1& q_2& p& {\beta }_{GC}& {\beta }_{GS}\end{array}\right\}}^T$, $y={\left\{\begin{array}{cccc}{\ddot{q}}_1& {\ddot{q}}_2& \ddot{p}& \delta \end{array}\right\}}^T$.

The objective function is defined as:

$$J=\frac{1}{2}{{\displaystyle \int }}_0^{\infty }\left[x^T\left(t\right)Qx\left(t\right)+u^T\left(t\right)Ru\left(t\right)\right]\mathrm{d}t$$

(47)

In Equation (47), $Q$ and $R$ are the weight matrices. In performing the Laplace transformation on Equation (46), the following is obtained:

$$sX\left(s\right)=AX\left(s\right)+BU\left(s\right)+x_0$$

(48)

Then:

$$X\left(s\right)={\left(sI-A\right)}^{-1}{x}_0+{\left(sI-A\right)}^{-1}BU\left(s\right)$$

(49)

Applying the inverse Laplace transformation to Equation (49) yields:

$$x\left(t\right)=L^{-1}\left[X\left(s\right)\right]=e^{At}{x}_0\left(t\right)+\underset{0}{\overset{t}{{\displaystyle \int }}}{e}^{A\left(t-\tau \right)}Bu\left(\tau \right)d\tau$$

(50)

Given that the feedback control signal satisfies the following relationship:

$$u\left(t\right)=-Kx\left(t\right)$$

(51)

Here, the state feedback control is defined, as represented by $K$ in Figure 7.

Substituting Equation (51) into Equation (46) yields:

$$\dot{x}\left(t\right)=A_{c}x\left(t\right)$$

(52)

where

$$A_c=A-BK$$

(53)

With the state feedback control, the natural vibration characteristics of the semi-span model are determined by matrix $A_c$. The solution of Equation (52) has the following form:

$$x\left(t\right)=e^{A_{c}t}{x}_0$$

(54)

Then, the objective function defined by Equation (47) becomes:

$$J=x_0^T\left[{{\displaystyle \int }}_0^{\infty }{e}^{A_c^{T}t}\left(Q+K^{T}RK\right)e^{A_c^{T}t}dt\right]x_0=x_0^{T}X{x}_0$$

(55)

According to the optimal control theory, the optimal controller can be obtained as:

$$K=R^{-1}{B}^{T}X$$

(56)

$X$ should satisfy the following Riccati equation:

$$A^{T}X+XA-XB{R}^{-1}{B}^{T}X+Q=0$$

(57)

According to Equation (51), the controller output is determined by state variable $x\left(t\right)$; however, it is not practically measurable. Therefore, a state observer is needed which will be described in the following section.

3.2. State Observer

Since it is convenient to measure the acceleration responses of the semi-span model in wind-tunnel tests, the accelerations are selected in this study to reconstruct the state variable $x\left(t\right)$ using a state observer. The dynamic equations of a full-order state observer can be expressed as:

$$\{\begin{array}{c}\dot{\widehat{x}}\left(t\right)=A\widehat{x}\left(t\right)+Bu\left(t\right)-H\left[\widehat{y}\left(t\right)-y\left(t\right)\right]\\ \widehat{y}\left(t\right)=C\widehat{x}\left(t\right)+Du\left(t\right)\end{array}$$

(58)

Then,

$$\begin{array}{r}\hfill \dot{\widehat{x}}\left(t\right)=A\widehat{x}\left(t\right)+Bu\left(t\right)-H\left[\widehat{y}\left(t\right)-y\left(t\right)\right]\\ \hfill =\left(A-HC\right)\widehat{x}\left(t\right)+Bu\left(t\right)+Hy\left(t\right)\end{array}$$

(59)

In Equation (59), $\left(A-HC\right)$ is defined as the system matrix of the state observation, and $H$ is an n-by-q matrix.

The estimation error of the state observer can be given by combining Equation (46) and Equation (59):

$$x\left(t\right)-\widehat{x}\left(t\right)=e^{\left(A-HC\right)\left(t-t_0\right)}\left[x\left(t_0\right)-\widehat{x}\left(t_0\right)\right]$$

(60)

The convergence of the estimation error relies on the configuration of matrix $A-HC$. Provided that the real parts of its eigenvalues are negative, the estimation error will converge rapidly, and the convergence rate depends on the pole configuration of the observer. The observer designed in this study is shown in Figure 8.

4. Simulation Results

4.1. Model Configuration

A simulation model was established by adding an actively controlled aileron to the Bell semi-span model as shown in Figure 9. The distance between the pitching axis of the aileron and the aerodynamic center of the airfoil is $l_2=0.5c$, while the distance between the wing root and the inner end of the aileron is $l_3=0.5R$. The maximum deflection amplitude of the aileron is limited to ±6°. Other characteristics of the simulation model, such as the stiffness and inertia configuration, are identical to the Bell model.

4.2. Whirl Flutter Suppression

A root locus was applied to evaluate the stability of the semi-span model. As shown in Figure 10a, if the LQR controller is turned off, the root locus, which corresponds to the wing flapping mode, crosses the imaginary axis and moves into the right half of the s-plane at a forward flight speed of 828 km/h, meaning whirl flutter occurs.

Figure 10b demonstrates the root locus of the semi-span model with the LQR controller turned on. In this condition, the root locus remains in the left half of the s-plane, illustrating that the system is stable and the actively controlled aileron combined with the LQR controller is effective in improving the whirl flutter stability of the tiltrotor aircraft.

Transient response analyses were performed at forward flight speeds of 920 km/h, 970 km/h and 990 km/h. At the 10th second, an acceleration disturbance was applied at the end of the wing to stimulate its whirl flutter, and the deflection angle of the aileron is shown in Figure 11. It could be seen that the aileron responded rapidly to suppress whirl flutter, demonstrating the effectiveness of the aileron as well as the established active controller.

5. Conclusions

In this study, a semi-span aeroelastic model was established for a tiltrotor aircraft to evaluate the influence of an aileron on whirl flutter phenomenon. In this model, the elastic deformation of the wing and rotor hub was taken into consideration, while the rotor blades and nacelle were treated as a rigid body. An active controller combined with a full-order state observer was built to control the deflection of the aileron. Numerical simulations were carried out at different forward flight speeds. The principal conclusions are as follows:

Firstly, the established aeroelastic model is capable of predicting the whirl flutter boundary of a tiltrotor aircraft. Secondly, a wing-mounted aileron, which was controlled by a controller based on an LQR algorithm and full-order state observer, is effective in increasing the whirl flutter stability. For the simulation model in this study, its boundary speed of whirl flutter was increased from 828 km/h to above 978 km/h.

Future works will focus on experimental research on whirl flutter suppression using an actively controlled aileron. A scaled model equipped with an aileron driven by piezoelectric actuators was proposed and wind-tunnel tests were scheduled in the following years.