Disturbance Observer Based Backstepping Terminal Sliding Mode Aeroelastic Control of Airfoils

1. Introduction

Nonlinear aeroelastic phenomena, such as limit cycle oscillation [1] and free-play flutter [2], not only threaten the aircraft structure, but can even result in loss of control. This challenge attracts many researchers to develop aeroelastic control for flutter suppression and avoidance of LCO [3].Usually passive flutter control has been widely used for wing flutter suppression in practical engineering, such as mass balance, local stiffness enhancement, etc., however, this method requires extra mass and makes the airplane heavier [4]. Therefore, active flutter control has been proposed and the objectives are to suppress flutter and to make the airplane flight exceed its flutter speed.

Active flutter control is a hot topic to suppress wing flutter by using modern control theories. Many researchers proposed useful schemes for active flutter suppression. Ouyang et al. studied the influence of flap camber on the wing flutter bounds, and proposed a Linear Quadratic Guassian (LQG) based control to actively suppress the aeroelastic response by the flap deflection, and the flutter speed and frequency of the wing will be changed by the orders of the reduced order model of the morphing flap camber [5]. Na et. al. proposed a LQG control method with sliding mode control (SMC) observer to active control of composite thin-walled beam structures in compressible flow, the dynamic aeroelastic responses of the beam structures subject to sharp-edged gust, blast loads and sonic-boom are investigated by this method [6]. Compared with linear aeroelasticity, nonlinear aeroelastic issue is more attractive, which includes freeplay, time delay and cubic nonlinearities in aeroelastic systems. Neil et al. studied structure nonlinear effect, analyzed and experimented by a rigid wing supported by a nonlinear spring [7]. Lee et al. presented the effect of a cubic structure restoring force on the flutter of a two dimensional airfoil placed in an incompressible flow, and the results showed the dependence of the divergence flutter boundary on initial conditions for a soft spring, while the independence of the divergence flutter boundary on initial conditions for the hard spring [8]. To stabilize a nonlinear aeroelastic system, Singh et al. proposed an adaptive backstepping output feedback control scheme, and trajectory control of the pitch angle is achieved by this method [9]. To deal with model uncertainties of aeroelastic systems, Lin et al. suggested a decoupled fuzzy sliding mode control (DFSMC), and the DFSMC system was designed to control the plunge and pitch motions simultaneously [10].Chen et. al. proposed a terminal SMC (TSMC)for aeroelastic systems, the LCO suffered in the aeroelastic system is removed by the TSMC [11]. Mukhopadhyay introduced historical perspective on analysis and control of aeroelasticity over the century, and showed some experiments of aeroelastic control for different aircrafts [12]. Bouma et al. investigated structural and aerodynamic nonlinearities effects for three degrees of freedom (DOF) aeroelastic systems, and stall effects were introduced using quasi-steady and unsteady aerodynamic representations [13]. Vishal et al. studied the synchronization characteristics of a pitch–plunge aeroelastic system, possessing discontinuous nonlinearities in both structural and aerodynamic fronts. Meanwhile, a multi-physics fluid–structure interaction problem possessing coupled non-smooth nonlinearities can be alleviated by the use of synchronization theory [14]. Livne gave a survey on aircraft active flutter suppression (AFS) through state of the art, and AFS wind tunnel test and flight-test programs were illustrated [15].

Flutter suppression theory has also been developed by the researchers. Considering delayed control problem, a single time delay was determined according to both the critical stability condition and the stability margins, Zou et al. proposed a delayed proportional–differential (PD) controller for robust suppression the body-freedom flutter (BFF) of a flying-wing drone [16].Trapani and Guo proposed a rudderless aeroelastic fin for a flexible blended-wing-body UAV to improve divergence and flutter margin [17]. Although these methods are useful for the linear time-invariant systems, they cannot guarantee the closed-loop performances of the nonlinear and time-varying aeroelastic system. For the structural nonlinearities and unsteady flow, Zhang and Behal proposed a continuous robust controller to realize aeroelastic vibration control of a two dimensional airfoil [18]. Lee and Singh proposed an L₁ adaptive control for the aeroelastic system with unsteady aerodynamics and gust load [19].For the nonlinear aeroelastic system, Platanitis and Strganac proposed a feedback linearization combined with model reference, and it can adaptively suppress LCO of a typical wing section with LE and TE control surfaces [20]. Prime set up a linear parameter varying (LPV) model for the three DOF Nonlinear Aeroelastic Test Apparatus (NATA), and developed a H₂ robust control by linear matrix inequalities [21]. Li et al. proposed a state-dependent Riccati equation based method to control an airfoil section with control surface nonlinearities, but it depended on the accurate model [22]. To reduce dependence on the model, Tang et al. proposed a neural network based active flutter suppression under control input constraints [23], but neural networks need extensive calculations. Recently Gao introduced sliding mode control into an aeroelastic system and presented a finite-time fault-tolerant control for wing flutter suppression [24].Chen et al. proposed a free-model deep reinforcement learning algorithm for the active flutter suppression through aeroelsatic wind tunnel testing, thus the jet flow intensity can be intelligently selected according to the real-time state of the flexible wing [25].

Since the aircrafts fly in a varying environment, the influence of model uncertainties and external disturbances on aeroelastic systems have also been considered by the researchers. Pettit briefly described aerodynamic and structure uncertainties and applications of uncertainty with design optimization [26]. For the nonlinear perturbations, Hao et al. analyzed stochastic nonlinear flutter behaviors of 3-DOF wing model by the stochastic p-bifurcation analysis method [27].Lee and Singh proposed a robust higher-order super-twisting control to deal with aeroelastic system with unsteady aerodynamics [28]. To estimate unknown parameters, Prabhu and Srinivas proposed an intelligent observer and then robust control of a three degrees of freedom aeroelastic model [29], but the neural network based observer has a heavy computational load. To improve the convergence speed of the aeroelastic system, Chen et al. proposed a terminal sliding mode control to reach the target within a finite time frame [30], and Yuan et al. proposed a sliding mode observer control for a two dimensional aeroelastic system with gusts [31]. Xu et al. also applied Backstepping Sliding Mode Control (BSMC) design to suppress the LCO with uncertainties and external disturbances [32]. All these methods can suppress or attenuate external disturbances and achieve the desired robust performances, however an alternative approach is available, the so-called disturbance observer based control (DOBC). The DOBC lumps all internal uncertainties and external disturbances acting on the aeroelastic system into a single term and then identifies them by the DO. The DO can provide fast, excellent control performances and smooth control actions without using large feedback gains [33]. Liu and Whidborne have proposed a nonlinear disturbance observer to deal with bounded unknown disturbances [34], but this observer does not cope well with variable disturbances, and the output tracking responses by this DO-BSMC cannot be guaranteed to converge within finite time.

Motivated by Xu’s [32] and Chen’s work [33], this paper considers trajectory tracking control for the two-dimensional airfoil-flap aeroelastic system. The major contributions of this work are as follows:

Considering quasi-steady and unsteady aerodynamic loads, nonlinear dynamic models are set up for an aeroelastic two-dimensional airfoil with both trailing-edge and leading-edge control surfaces.

To suppress flutter, a backstepping terminal sliding-mode control is proposed for the two-dimensional airfoil-flap aeroelastic system. The benefits of this approach have TSMC merits of high robustness, fast transient response and finite time convergence, as well as backstepping control in terms of globally asymptotic stability.

To overcome the limitation of the BTSMC being dependent on prior knowledge, a nonlinear disturbance observer is designed for coping with model uncertainties and variable external bounded gust disturbances, which is different from Backstepping SMC for flutter suppression [32].

To adjust the speed where the aeroelastic response converges to zero, a learning rate is incorporated into the BTSMC design.

Simulation results show the proposed DO-BTSMC method is effective and has advantages in LCO suppression and flutter control for a two-dimensional wing with LE and TE control surfaces.

This paper is organized as follows. Section 2 establishes the nonlinear aeroelastic dynamics model of the two-dimensional airfoil-flap system. Section 3 proposes a DO-BTSMC design, and stability is analyzed for the closed-loop control system. Simulations and analysis are in Section 4. Some conclusions are given in Section 5.

2. Aeroelastic System Dynamic Model

The studied model is a two-DOF pitch-plunge wing section with both LE and TE control surfaces and is shown in Figure 1, both control surfaces are used as control inputs.

The physical parameters of the wing are listed as follows, a is nondimensional distance from midchord to elastic axis position, b and c denote semichord and full chord length of wing section respectively.

The aeroelastic governing models are established as follows [20]

$$\left[\begin{array}{cc}{m}_T& m_w{x}_{\alpha }b\\ m_w{x}_{\alpha }b& I_{\alpha }\end{array}\right]\left[\begin{array}{c}\ddot{h}\\ \ddot{\alpha }\end{array}\right]+\left[\begin{array}{cc}{c}_h& 0\\ 0& c_{\alpha }\end{array}\right]\left[\begin{array}{c}\dot{h}\\ \dot{\alpha }\end{array}\right]+\left[\begin{array}{cc}{k}_h& 0\\ 0& k_{\alpha }(\alpha )\end{array}\right]\left[\begin{array}{c}h\\ \alpha \end{array}\right]=\left[\begin{array}{c}-L\\ M\end{array}\right],$$

(1)

where

$$L=\rho V^{2}b{s}_p{C}_{L\alpha }\left(\alpha +{\displaystyle \frac{\dot{h}}{V}}+\left({\displaystyle \frac{1}{2}}-a\right)b{\displaystyle \frac{\dot{\alpha }}{V}}\right)+\rho V^{2}b{s}_p{C}_{L\beta }\beta +\rho V^{2}b{s}_p{C}_{L\gamma }\gamma ,$$

(2)

$$M=\rho V^2{b}^2{s}_p{C}_{m\alpha ,eff}\left(\alpha +{\displaystyle \frac{\dot{h}}{V}}+\left({\displaystyle \frac{1}{2}}-a\right)b{\displaystyle \frac{\dot{\alpha }}{V}}\right)+\rho V^2{b}^2{s}_p{C}_{m\beta ,eff}\beta +\rho V^2{b}^2{s}_p{C}_{m\gamma ,eff}\gamma ,$$

(3)

C_{mα, eff} , C_{mβ, eff} , and C_{mγ, eff} are the effective dynamic and control moment derivatives due to angle of attack and trailing- and leading-edge control surface deflection, respectively, about the elastic axis and are defined as follows

$$C_{m\alpha ,eff}=\left(0.5+a\right)C_{L\alpha }+2C_{m\alpha },C_{m\beta ,eff}=\left(0.5+a\right)C_{L\beta }+2C_{m\beta },$$

$$C_{m\gamma ,eff}=\left(0.5+a\right)C_{L\gamma }+2C_{m\gamma }$$

(4)

m_T denote total mass of pitch–plunge system, c_α and c_h denote damping of pitch and plunge respectively, k_α and k_h denote pitch and plunge stiffness respectively, s_p is wing section span, I_α denotes total pitch moment of inertia about elastic axis,ρ denotes air density, h is plunge displacement, V is freestream velocity, α is angle of attack (AOA),β and γ trailing- and leading-edge control surface deflection respectively, C_L, C_Lα, C_Lβ, C_Lγ denote wing section lift coefficients and C_m, C_mα, C_mβ, C_mγ denote wing section pitch moment coefficients at quarter-chord.

Denote the state variable vector is

$$x_1={\left[\begin{array}{cc}h& \alpha \end{array}\right]}^T,x_2={\dot{x}}_1={\left[\begin{array}{cc}\dot{h}& \dot{\alpha }\end{array}\right]}^T,x={\left[\begin{array}{cc}{x}_1^T& x_2^T\end{array}\right]}^T,$$

(5)

$${\dot{x}}_2(t)=f_2(x)+g_2(x)U(t),$$

(6)

$$f_2(x)=\left[\begin{array}{cc}-k_1& -k_2{V}^2+p(x_2)\\ -k_3& -k_4{V}^2+q(x_2)\end{array}\right]x_1+\left[\begin{array}{cc}-c_1& -c_2\\ -c_3& -c_4\end{array}\right]x_2,U=\left[\begin{array}{c}\beta \\ \gamma \end{array}\right],g_2(x)=\left[\begin{array}{cc}{g}_{2,11}& g_{2,12}\\ g_{2,21}& g_{2,22}\end{array}\right],$$

$$\Xi =m_T{I}_a-{(m_W{x}_{\alpha }b)}^2,k_1=I_a{k}_h/\Xi ,k_3=-m_W{x}_{\alpha }b{k}_h/\Xi ,$$

$$k_2=\left(I_a\rho s_{p}b{C}_{L\alpha }+m_W{x}_{\alpha }\rho s_p{b}^3{C}_{m\alpha }\right)/\Xi ,k_4=-\left(m_W{x}_{\alpha }\rho s_p{b}^2{C}_{L\alpha }+m_T\rho s_p{b}^2{C}_{m\alpha ,eff}\right)/\Xi ,$$

$$c_1=\left(I_{\alpha }\left(c_h+\rho V{s}_{p}b{C}_{L\alpha }\right)+m_W{x}_{\alpha }\rho V{s}_p{b}^3{C}_{m\alpha ,eff}\right)/\Xi ,$$

$$c_2=\left(I_{\alpha }\rho V{s}_p{b}^2{C}_{L\alpha }\left(0.5-a\right)-m_W{x}_{\alpha }b{c}_{\alpha }\right.\left.+m_W{x}_{\alpha }\rho V{s}_p{b}^4{C}_{m\alpha ,eff}\left(0.5-a\right)\right)/\Xi ,$$

$$c_3=\left(-m_W{x}_{\alpha }b\left(c_h+\rho Vb{s}_p{C}_{L\alpha }\right)-m_T\rho V{s}_p{b}^2{C}_{m\alpha ,eff}\right)/\Xi ,$$

$$c_4=\left[m_T\left(c_{\alpha }-\rho V{s}_p{b}^3{C}_{m\alpha ,eff}\left(0.5-a\right)\right)\right.\left.-m_W{x}_{\alpha }\rho V{s}_p{b}^3{C}_{L\alpha }\left(0.5-a\right)\right]/\Xi ,$$

$$g_{2,11}=\left(-I_{\alpha }\rho b{s}_p{C}_{L\beta }-m_W{x}_{\alpha }\rho s_p{b}^3{C}_{m\beta ,eff}\right)V^2/\Xi ,$$

$$g_{2,12}=\left(-I_{\alpha }\rho b{s}_p{C}_{L\gamma }-m_W{x}_{\alpha }\rho s_p{b}^3{C}_{m\gamma ,eff}\right)V^2/\Xi ,$$

$$g_{2,21}=\left(m_W{x}_{\alpha }\rho s_p{b}^2{C}_{L\beta }+m_T\rho s_p{b}^2{C}_{m\beta ,eff}\right)V^2/\Xi ,$$

$$g_{2,22}=\left(m_W{x}_{\alpha }\rho s_p{b}^2{C}_{L\gamma }+m_T\rho s_p{b}^2{C}_{m\gamma ,eff}\right)V^2/\Xi ,$$

(7)

and uncertain functions of plunge and AOA motions are

$$q(x_2)=m_T{k}_{\alpha }/\Xi ,p(x_2)=-m_W{x}_{a}b{k}_a/\Xi .$$

(8)

The servo motor dynamics for the TE and LE can be modeled as a second-order system of the form

$$\left\{\begin{array}{l}{I}_{\beta }\ddot{\beta }+c_{\beta s}\dot{\beta }+k_{\beta s}\beta =k_{\beta s}{\beta }_c\\ I_{\gamma }\ddot{\gamma }+c_{\gamma s}\dot{\gamma }+k_{\gamma s}\gamma =k_{\gamma s}{\gamma }_c\end{array}\right.$$

(9)

where$c_{\beta s}$and$c_{\gamma s}$are the damping of the TE and LE respectively, $k_{\beta s}$ and$k_{\gamma s}$ are stiffness of the TE and LE respectively, and${\beta }_c$and${\gamma }_c$denote the desired positions of the servo motors for the TE and LE respectively.$m_{\alpha }$and $I_{\alpha 0}$denote the rotational wingbody of mass and inertia respectively, $I_{\alpha }=I_{\alpha 0}+m_{\alpha }{x}_{\alpha }^2{b}^2$,$I_{\beta }=I_{\beta 0}+m_{\beta }{r}_{\beta }^2$ and $I_{\gamma }=I_{\gamma 0}+m_{\gamma }{r}_{\gamma }^2$are the moments of inertia of the TE and LE about its pivot respectively , $m_{\beta }$,$m_{\gamma }$,$I_{\beta 0}$ and $I_{\gamma 0}$are the rotational TE and LE section of mass and inertia respectively. Let$\overline{m}=m_h+m_{\alpha }+m_{\beta }+m_{\gamma }$,$m_h$denotes a translational body of mass.

When the aerodynamic performance of the airfoil is unsteady, the equations of motion for the entire system can be set up by Lagrangian energy theory,

$${\displaystyle \frac{d}{dt}}\left[\begin{array}{c}\partial \mathcal{L}/\partial h\\ \partial \mathcal{L}/\partial \alpha \\ \partial \mathcal{L}/\partial \beta \\ \partial \mathcal{L}/\partial \gamma \end{array}\right]-\left[\begin{array}{c}\partial \mathcal{L}/\partial h\\ \partial \mathcal{L}/\partial \alpha \\ \partial \mathcal{L}/\partial \beta \\ \partial \mathcal{L}/\partial \gamma \end{array}\right]=-\left[\begin{array}{c}0\\ 0\\ k_{\beta s}\beta \\ k_{\gamma s}\gamma \end{array}\right]-\left[\begin{array}{c}{c}_h\dot{h}\\ c_{\alpha }\dot{\alpha }\\ c_{\beta s}\dot{\beta }\\ c_{\gamma s}\dot{\gamma }\end{array}\right]+\left[\begin{array}{c}-L\\ M_{\alpha }\\ M_{\beta }\\ M_{\gamma }\end{array}\right]+\left[\begin{array}{c}0\\ 0\\ k_{\beta s}{\beta }_c\\ k_{\gamma s}{\gamma }_c\end{array}\right]$$

(10)

where $\mathcal{L}$denotes Lagrangian function. The expanded form of (10) is

$$\left[\begin{array}{cccc}\overline{m}& {\xi }_1& m_{\beta }{r}_{\beta }& m_{\gamma }{r}_{\gamma }\\ *& {\xi }_2& I_{\beta }+x_{\beta }{m}_{\beta }{r}_{\beta }& I_{\gamma }+x_{\gamma }{m}_{\gamma }{r}_{\gamma }\\ *& *& I_{\beta }& 0\\ *& *& 0& I_{\gamma }\end{array}\right]\left[\begin{array}{c}\ddot{h}\\ \ddot{\alpha }\\ \ddot{\beta }\\ \ddot{\gamma }\end{array}\right]+\left[\begin{array}{c}{k}_{h}h+c_h\dot{h}+{\xi }_3\\ k_{\alpha }\alpha +c_{\alpha }\dot{\alpha }+{\xi }_4\\ k_{\beta s}\beta +c_{\beta \mathrm{s}}\dot{\beta }+{\alpha }^2{m}_{\beta }{r}_{\beta }{x}_{\beta }\sin \beta \\ k_{\gamma s}\gamma +c_{\gamma \mathrm{s}}\dot{\gamma }+{\alpha }^2{m}_{\gamma }{r}_{\gamma }{x}_{\gamma }\sin \gamma \end{array}\right]=\left[\begin{array}{c}-L\\ M_{\alpha }\\ M_{\beta }\\ M_{\gamma }\end{array}\right]+\left[\begin{array}{c}0\\ 0\\ k_{\beta s}{\beta }_c\\ k_{\gamma s}{\gamma }_c\end{array}\right]$$

(11)

wherea * denotes symmetric matrix element,

$$\begin{gathered} {\xi }_1=(m_{\alpha }{r}_{\alpha }+m_{\beta }{r}_{\beta }+m_{\gamma }{r}_{\gamma })\cos \alpha +m_{\beta }{r}_{\beta }\cos (\alpha +\beta )+m_{\gamma }{r}_{\gamma }\cos (\alpha +\gamma ), \\ {\xi }_2=I_{\alpha }+I_{\beta }+I_{\gamma }+m_{\beta }{x}_{\beta }^2+m_{\gamma }{x}_{\gamma }^2+2x_{\beta }{m}_{\beta }{r}_{\beta }\cos \beta +2x_{\gamma }{m}_{\gamma }{r}_{\gamma }\cos \gamma , \\ {\xi }_3=-{(\alpha +\beta )}^2{m}_{\beta }{r}_{\beta }\sin (\alpha +\beta )-{(\alpha +\gamma )}^2{m}_{\gamma }{r}_{\gamma }\sin (\alpha +\gamma )-{\alpha }^2\sin \alpha (m_{\alpha }{r}_{\alpha }+m_{\beta }{x}_{\beta }+m_{\gamma }{x}_{\gamma }), \\ {\xi }_4=-\beta (2\alpha +\beta )m_{\beta }{r}_{\beta }{x}_{\beta }\sin (\beta )-\gamma (2\alpha +\gamma )m_{\gamma }{r}_{\gamma }{x}_{\gamma }\sin (\gamma ) \end{gathered}$$

(12)

The aerodynamic force and moment in (8) can be calculated by Theodorsen method [35,36]

$$L=\pi \rho b^2{s}_p\left(\ddot{h}+V\dot{\alpha }-ba\ddot{\alpha }-(V/\pi )T_4\dot{\beta }-(b/\pi )T_1\ddot{\beta }+(V/\pi )T_4\dot{\gamma }-(b/\pi )T_1\ddot{\gamma }\right)+2\pi \rho Vb{s}_{p}C(k)\left(V\alpha +\dot{h}+b(0.5-a)\dot{\alpha }+(V/\pi )T_{10}\beta +(b/2\pi )T_{11}\dot{\beta }+(V/\pi )T_{10}\gamma +(b/2\pi )T_{11}\dot{\gamma }\right),$$

(13)

$$M_{\alpha }=\pi \rho b^2{s}_p\left(ba\ddot{h}-Vb(0.5-a)\dot{\alpha }-b^2(1/8+a)\ddot{\alpha }-(V^2/\pi )T_{15}(\beta \right.+\gamma )\left.-(Vb/\pi ){\mathrm{T}}_{16}(\dot{\beta }+\dot{\gamma })-2(b^2/\pi )T_{13}(\ddot{\beta }+\ddot{\gamma })\right)+2\pi \rho V{b}^2{s}_p(a+0.5)C(k)\left(V\alpha +\dot{h}+b(0.5-a)\dot{\alpha }+(V/\pi )T_{10}(\beta +\gamma )+(b/2\pi )T_{11}(\dot{\beta }+\dot{\gamma })\right),$$

(14)

$$M_{\beta }=-\rho b^2{s}_p\left(\left(-2T_9-T_1+T_4(a-0.5)\right)Vb\dot{\alpha }+2T_{13}{b}^2\ddot{\alpha }+(1/\pi )V^2(\beta +\gamma )(T_5-T_4{T}_{10})\right.\left.-(1/2\pi )Vb(\dot{\beta }+\dot{\gamma })T_4{T}_{11}-(1/\pi )T_{3}b(\ddot{\beta }+\ddot{\gamma })-T_{1}b\ddot{h}\right)-\rho V{b}^2{s}_p{T}_{12}C(k)\left(V\alpha +\dot{h}+b(0.5-a)\dot{\alpha }+(1/\pi )T_{10}V(\beta +\gamma )+b(1/2\pi )T_{11}(\dot{\beta }+\dot{\gamma })\right)$$

(15)

where T_i (i =1,2,…,14) are the Theodorsen constants, see Appendix B [37]. $M_{\gamma }$is derived from those by Throdorsen [35,36,37].

The aerodynamic force and moments in (13)-(15) are dependent on Theodorsen function C(k), where k is the non-dimensional reduced frequency of harmonic oscillation, C(k) can be approximated by Jones method,

$$C(s)=1-{\displaystyle \frac{0.165s}{s+0.0455V/b}}-{\displaystyle \frac{0.335s}{s+0.3V/b}}={\displaystyle \frac{1}{2}}+{\displaystyle \frac{c_{1}V/b}{s+c_{3}V/b}}+{\displaystyle \frac{c_{2}V/b}{s+c_{4}V/b}},$$

(16)

where c₁=0.0075, c₂=0.10055, c₃=0.0455, c₄=0.3. C(s) can be presented as a state-space model as follows

$$\left\{\begin{array}{l}\left[\begin{array}{c}{\dot{x}}_{a1}\\ {\dot{x}}_{a2}\end{array}\right]={\displaystyle \frac{V}{b}}\left[\begin{array}{cc}-c_3& 0\\ 0& -c_4\end{array}\right]\left[\begin{array}{c}{x}_{a1}\\ x_{a2}\end{array}\right]+\left[\begin{array}{c}1\\ 1\end{array}\right]r\\ y_a={\displaystyle \frac{V}{b}}\left(\begin{array}{cc}{c}_1& c_2\end{array}\right)\left[\begin{array}{c}{x}_{a1}\\ x_{a2}\end{array}\right]+0.5r\end{array}\right.,$$

(17)

where

$$r=V\alpha +\dot{h}+b(0.5-a)\dot{\alpha }+(V/\pi )T_{10}\beta +b(1/2\pi )\dot{\beta },$$

(18)

The aerodynamic force and moments can be described as the sum of non-circulatory and circulatory loads,

$$F_{ac}=F_{nc}+F_c,$$

(19)

where F_nc and F_c are non-circulatory load and circulatory load respectively. Denote${\overline{x}}_1={\left[\begin{array}{cccc}h& \alpha & \beta & \gamma \end{array}\right]}^T$, x_a = [x_a1, x_a2]^T, then we get

$$F_{nc}=-M_{nc}{\ddot{\overline{x}}}_1-V{C}_{nc}{\dot{\overline{x}}}_1-V^2{K}_{nc}{\overline{x}}_1,$$

(20)

$$F_c=-\rho Vb{s}_p{R}_c\left({\displaystyle \frac{V}{b}}\left(\begin{array}{cc}{c}_1& c_2\end{array}\right)x_a+0.5\left(S_1{\dot{\overline{x}}}_1+V{S}_2{\overline{x}}_1\right)\right),$$

(21)

where x_α is nondimensional distance from elastic axis to centre of mass, and meeting x_α= r_cg /b, r_cg is distance from elastic axis to centre of mass.

$$\begin{gathered} M_{nc}=\rho b^2{s}_p\left[\begin{array}{cccc}\pi & -\pi ba& -b{T}_1& -b{T}_1\\ -\pi b\overline{a}& \pi b^2(1/8+a^2)& 2b^2{T}_{13}& -b^2{T}_{13d}\\ -b{T}_1& 2b^2{T}_{13}& -T_3{b}^2/\pi & 0\\ -b{T}_1& 2b^2{T}_{13d}& 0& -T_3{b}^2/\pi \end{array}\right], \\ {K}_{nc}=\rho b^2{s}_p\left[\begin{array}{cccc}0& 0& 0& 0\\ 0& 0& T_{15}& T_{15}\\ 0& 0& T_{18}/\pi & 0\\ 0& 0& 0& T_{18}/\pi \end{array}\right], \\ {R}_c={\left[\begin{array}{cccc}2\pi & -2\pi b(0.5+a)& T_{12}& T_{12}\end{array}\right]}^T, \\ {S}_1=\left[\begin{array}{cccc}1& b(0.5-a)& b{T}_{11}/(2\pi )& b{T}_{11}/(2\pi )\end{array}\right], \\ {S}_2=\left[\begin{array}{cccc}0& 1& {\mathrm{T}}_{10}/\pi & {\mathrm{T}}_{10}/\pi \end{array}\right], \end{gathered}$$

(22)

$T_{13}$and$T_{13d}$are the Theodorsen constants for the TE and LE respectively. Substitute (17) , (13) - (15) into (11), then yields

$$\left(M_s+M_{nc}\right){\ddot{\overline{x}}}_1+\left(C_s+V{C}_{nc}+V{C}_c\right){\dot{\overline{x}}}_1+\left(K_s+V^2{K}_{nc}+V^2{K}_c\right){\overline{x}}_1=-V^2{E}_c{x}_a+G_s{u}_c$$

(23)

where(3)

$$C_{nc}=\rho b^2{s}_p\left[\begin{array}{cccc}0& \pi & -T_4& T_4\\ 0& \pi b(0.5-\overline{a})& b{T}_{16}& b{T}_{16}\\ 0& b{T}_{17}& b{T}_{19}/\pi & 0\\ 0& b{T}_{17}& 0& b{T}_{19}/\pi \end{array}\right],C_s=\left[\begin{array}{cccc}{c}_h& 0& 0& 0\\ 0& c_{\alpha }& 0& 0\\ 0& 0& c_{\beta \mathrm{s}}& 0\\ 0& 0& 0& c_{\gamma \mathrm{s}}\end{array}\right],K_s=\left[\begin{array}{cccc}{k}_h& 0& 0& 0\\ 0& k_{\alpha }& 0& 0\\ 0& 0& k_{\beta s}& 0\\ 0& 0& 0& k_{\gamma s}\end{array}\right],G_s=\left[\begin{array}{cc}0& 0\\ 0& 0\\ k_{\beta s}& 0\\ 0& k_{\gamma s}\end{array}\right],M_s=\left[\begin{array}{cccc}\overline{m}& {\xi }_1& m_{\beta }{r}_{\beta }& m_{\gamma }{r}_{\gamma }\\ *& {\xi }_2& I_{\beta }+x_{\beta }{m}_{\beta }{r}_{\beta }& I_{\gamma }+x_{\gamma }{m}_{\gamma }{r}_{\gamma }\\ *& *& I_{\beta }& 0\\ *& *& 0& I_{\gamma }\end{array}\right],C_c={\displaystyle \frac{1}{2}}\rho b{s}_p{R}_c{S}_1,K_c={\displaystyle \frac{1}{2}}\rho b{s}_p{R}_c{S}_2,E_c=\rho s_p{R}_c\left[\begin{array}{cc}{c}_1& c_2\end{array}\right],u_c=[{\beta }_c,{\gamma }_c{]}^{\mathrm{T}}$$

(24)

Denote the full state vector$\overline{x}={\left[\begin{array}{ccc}{{\overline{x}}_1}^T& {{\dot{\overline{x}}}_1}^T& x_a^T\end{array}\right]}^T$, the overall state-space model is

$$\dot{\overline{x}}=A\overline{x}+Bu$$

(25)

where the A and B matrix are as follows

$$A=\left[\begin{array}{ccc}0& I& 0\\ -M^{-1}K(V)& -M^{-1}C(V)& -V^2{M}^{-1}{E}_c\\ V{K}_a& C_a& V{A}_a\end{array}\right],B=\left[\begin{array}{c}0\\ M^{-1}{G}_s\\ 0\end{array}\right],$$

(26)

where

$$M=M_{nc}+M_s,K(V)=K_s+V^2{K}_{nc}+V^2{K}_c,C(V)=C_s+V{C}_{nc}+V{C}_c,K_a=\left[\begin{array}{c}1\\ 1\end{array}\right]S_2,C_a=\left[\begin{array}{c}1\\ 1\end{array}\right]S_1,A_a={\displaystyle \frac{1}{b}}\left[\begin{array}{cc}-c_3& 0\\ 0& -c_4\end{array}\right].$$

(27)

Since there exists wind disturbance for the aeroelastic system, then the external disturbance D is considered and added into the system of (6) and (25) as

$${\dot{x}}_2(t)=f_2(x)+g_2(x)U(t)+D(t),$$

(28a)

$${\ddot{\overline{x}}}_1(t)=-M^{-1}\left(K(V){\overline{x}}_1+C(V){\dot{\overline{x}}}_1+V^2{E}_c{x}_a\right)+M^{-1}{G}_{s}u+\overline{D}(t),$$

(28b)

where$D\in R^2$and $\overline{D}(t)\in R^4$denote external disturbances and bound by$‖D‖<D_U$, D_U is a real value meeting D_U≥0. Denote${\overline{f}}_2({\overline{x}}_1)=-M^{-1}\left(K(V){\overline{x}}_1+C(V){\dot{\overline{x}}}_1+V^2{E}_c{x}_a\right)$ ,

${\overline{g}}_2({\overline{x}}_1)=M^{-1}{G}_s$. The aerodynamic loads due to the wind disturbances are defined as [38]

$$L_g=\rho V^{2}b{s}_p{C}_{L\alpha }{\omega }_G(\tau )/V=\rho Vb{s}_p{C}_{L\alpha }{\omega }_G(\tau ),$$

(29)

$$M_g=(0.5+a)b{L}_g,$$

(30)

where${\omega }_G(\tau )$denotes the disturbance velocity, τ is a dimensionless time variable defined as $\tau =Vt/b$, so $D={\left[\begin{array}{cc}-L_g/m_T& M_g/I_{\alpha }\end{array}\right]}^T$.Suppose that the above disturbance can be generated by a linear exogenous system [33],

$$\left\{\begin{array}{l}{\dot{\varsigma }}_d(t)=W_d{\varsigma }_d(t)\\ D(t)=V_d{\varsigma }_d(t)\end{array}\right.,$$

(31)

where${\varsigma }_d\in R^p$,W_d, V_d are coefficient matrices with corresponding dimensions.

Airfoil aeroelastic vibration Control Problem: The control objective is to suppress the aeroelastic vibrations, so the desired smooth trajectory x_d can be chosen to be zero for all time. Alternatively the desired trajectory x_d along with the actual pitching and plunging variables is driven toward the origin.

3. DO-BTSMC Aeroelastic Control Design

This section gives an overview of the DO-BTSMC control for the two dimensional airfoil.

3.1. Nonlinear Disturbance Observer Design

Considering unknown disturbances in the aeroelastic system (28), a nonlinear disturbance observer (NDO) is proposed to improve the tracking precision [33],

$$\left\{\begin{array}{l}{\dot{z}}_d(t)=(W_d-l(x)V_d)z_d(t)+W_d{p}_d(x)-l(x)\left(V_d{p}_d(x)+f_2(x)+g_{2}U(t)\right)\\ {\widehat{\varsigma }}_d(t)=z_d(t)+p_d(x)\\ \widehat{D}(t)=V_d{\widehat{\varsigma }}_d(t)\end{array}\right.,$$

(32)

where$\widehat{D}$is the estimated disturbance, l (x) is the vector of observer gains to be tuned for performance, z_d is the internal state vector, and p_d(x) is a nonlinear function to be designed. The NDO gain l(x) is determined by

$$l(x)={\displaystyle \frac{\partial p_d(x)}{\partial x}}$$

(33)

Integrating l(x₂) with respect to the aeroelastic system state x yields,

$$p_d(x)=L_{d}Li{e}_f^{r-1}\hslash (x)$$

(34)

where L_d is a gain vector, Lie_f denotes a Lie derivative of function f(x), and r is the relative degree. Select$\hslash (x_2)=x_2$. The nonlinear disturbance estimation error is$d=D-\widehat{D}$.

It has been shown that the estimation$\widehat{D}$of the NDO approaches the disturbance D(t) exponentially if the observer gain l(x) is chosen such that (32) is global exponentially stable for all $x\in R^n$ [39]. The estimation of harmonic disturbance, $\widehat{D}$, yielded by the DO of (32) approaches the disturbance, D, and is global exponentially stable if the following inequalities are satisfied [33]

$$V={\displaystyle \frac{1}{2{\gamma }_d}}{d}^{T}d>0,\dot{V}<0$$

(35)

Hence the disturbance compensator can be designed as

$$u_d(t)=-g_2^{-1}\widehat{D}(t)$$

(36)

By using the nonlinear disturbance observer, $\widehat{D}$will approach D, so the disturbance can be attenuated through the disturbance compensator (36) and the system (28) can be transformed as follows

$${\dot{x}}_2(t)=f_2(x_1,x_2)+g_{2}U(t)+D=f_2(x_1,x_2)+g_2\left(U_{BTSMC}(t)+u_d(t)\right)+D(t)=f_2(x)+g_2(x)U_{BTSMC}(t)+d$$

(37)

The observable disturbances usually can be measured and estimated by the observer, such as 1-cosine winds. It has been shown that the estimation$\widehat{D}(t)$ of the NDO approaches the disturbance D(t) exponentially if the error dynamics for d(t) is global exponentially stable [33]. By the similar method, the disturbances $\overline{D}(t)$for unsteady loads can also be estimated.

3.2. DO-BTSMC Based Aeroelastic Control Design

To drive the aeroelastic output x to converge to the desired value vector x_d, the aeroelastic controller is designed as follows.

The tracking error vector of the plunge displacement and pitch angle is defined as

$$z_1=x_1-x_{1d}$$

(38)

and its derivative is

$${\dot{z}}_1={\dot{x}}_1-{\dot{x}}_{1d}=x_2-{\dot{x}}_{1d}$$

(39)

Define the stabilizing function${\varphi }_1(x)$, as the virtual speed vector of${\dot{x}}_{1d}$

$${\varphi }_1(x_1)={\dot{x}}_{1d}-K_1{z}_1$$

(40)

where${\dot{x}}_{1d}$represents the reference or virtual speed vector, K₁>0 is often chosen as a diagonal matrix to simplify the design, i.e., K₁ =diag(k₁₁, k₁₂), k_1i (i = 1, 2) is constant value.

Construct Lyapunov function V₁ as follows

$$V_1={\displaystyle \frac{1}{2}}{z}_1^T{z}_1$$

(41)

Define the speedtracking error vector for the plunge and pitch dynamics as

$$z_2={\dot{x}}_1-{\varphi }_1(x_1)=x_2-{\varphi }_1(x_1)$$

(42)

and its derivative is

$${\dot{z}}_2={\ddot{x}}_1-{\ddot{x}}_{1d}+K_1{\dot{z}}_1$$

(43)

The detailed derivation is given in the Appendix A.

Remark 1 If there are uncertainties for the system (28), then f₂ and g₂ will include uncertain terms, the uncertain terms can be regarded as disturbances and added into D, and then the observer (32) can be used with nominal f₂ and g_2.

Therefore, the structure of the DO-BTSMC controller for aeroelastic control is designed as Figure 2, which includes aeroelastic dynamics system of the two dimensions airfoil, an adaptive BTSMC used for aeroelastic flutter suppression, a nonlinear DO used to estimate observable disturbances, and an adaptive disturbance estimator used to observe the unobservable disturbances. A learning rate is applied to speed up the convergence of the proposed sliding mode control.

4. Discussion

The studied aeroelastic system is the NATA model and the physical parameters of the aeroelastic system are listed in Table 1 [20,21,40].

4.1. Flutter Suppression

In this section, flutter suppression of the airfoil is studied. For the two dimensional airfoil, the polynomial model of the stiffness is found by static measurements on the nonlinear pitch cam, k_α(α) =12.77+53.47α+1003α²Nm/rad [20]. When the airspeed is increased up to 11.4m/s, the LCO phenomenon appears. The output responses are shown in Figure 3. To suppress this LCO and flutter vibrations, different scenarios are studied as follows. First the proposed BTSMC is used to suppress this flutter. The initial conditions of the two-dimensional aeroelastic airfoil are h(0)= −0.01m, α(0)=0.1rad,$\dot{h}(0)=-0.1$(m/s), $\dot{\alpha }(0)=-2$(rad/s), the nominal airspeed V =11.4m/s. The parameters of the BTSMC are designed after several design iterations as follows,K₁=diag(0.2,0.1), h_s=diag(0.5,0.5), ς=diag(30,45),φ_s=0.4,β_s=1,p=5, q=3, β=diag(1,1), α₀=10, β₁=3, λ=1.8*diag(1,1). By using the proposed BTSMC design, the simulation results refer to Figure 4.

It can be seen from Figure 4 that the dynamic responses of the proposed nonlinear aeroelastic system flutter under the initial input in the first 2 sec, and the resulting dynamic responses under the BTSMC are stable after 2sec. The FFT results of plunge displacements and AOA with and without BTSMC control are shown in Figure 3 and Figure 4, the vibration frequency is 2.7Hz (see Figure 3), and the amplitude of h_FFT changes from 1.5×10^-3 to 1.25×10^-4, while the amplitude of α_FFT changes from 0.1 to 7.4×10^-3, which shows the flutter is suppressed by the proposed BTSMC.

For high stiffness nonlinearity in pitch, k_α(α) can be presented as the following polynomial type [21],

The flutter speed is related to its stiffness, the flutter equation for the p-k method can be written as follows by Eq. (23) [41],

$$\left[\left({\displaystyle \frac{V}{L_{ref}}}\right)\mathrm{M}{p}^2+\mathrm{K}-{\displaystyle \frac{1}{2}}\rho V^2\left({\mathrm{Q}}^R(M,k)+{\displaystyle \frac{p}{k}}{\mathrm{Q}}^I(M,k)\right)\right]\cdot \left\{\mathrm{q}\right\}=0$$

(44)

where M and K are the generalized mass and stiffness matrices, respectively; Q^R and Q^I are the real and imaginary parts of the generalized aerodynamic force matrix that depends on Mach number M and reduced frequency k. V is the velocity of undisturbed flow, L_ref is the reference length. q denotes eigenvector of flutter equation, p = g_s+iω_s is non-dimensional Laplace complex parameter and is also refered to as an eigenvalue of Eq.(44). The imaginary part of p must meet the following consistency condition, $k=\left({\displaystyle \frac{L_{ref}}{V}}\right)\cdot \mathrm{Im}(p)$, i.e.,$k=\left({\displaystyle \frac{L_{ref}}{V}}\right)\cdot \omega$, and $\omega =\sqrt{{\displaystyle \frac{\mathrm{K}}{\mathrm{M}}}}$. When the damping curve {(V_s,g_s)} for the vibration mode intersects with g_s =0, the associated speed is flutter speed. Therefore, for a given k, when the stiffness K is increased, the frequence is increased, and the flutter speed Vs is increased. Here the high stiffness is applied, the flutter speed is increased from V = 11.4m/s to V=13.2m/s. By using this stiffness the simulation results refer to Figure 5 and Figure 6,

It can be seen from Figure 5 and Figure 6 that the flutter speed is increased from V=11.4 m/s to V=13.2m/s by the high stiffness nonlinearity, and the output responses of plunge and pitch motions are quickly convergent by the proposed BTSMC control.

Considering high angle of attack (AOA) in the post-flutter behavior, a command of 10deg square wave is tracked and the output responses are shown as Figure 7,

It can be seen from Figure 7 that 10deg AOA command has been accurately tracked by the proposed BTSMC, but it is easy to make the control surfaces saturation and results in small oscillations of the output responses.

4.2. Flutter Suppression under Unknown Wind Disturbances

Case I: for quasi-steady aerodynamic loads

There are two kinds of gusts to be studied here, that is, a triangular-type gust and a sine gust as follows [32,38]

$$\begin{gathered} w_G(\tau )=2w_0{\displaystyle \frac{\tau }{{\tau }_G}}\left(H(\tau )-H(\tau -{\displaystyle \frac{{\tau }_G}{2}})\right)+2w_0\left({\displaystyle \frac{\tau }{{\tau }_G}}-1\right)\left(H(\tau -{\tau }_G)-H(\tau -{\displaystyle \frac{{\tau }_G}{2}})\right) \\ \left(\mathrm{triangular\; gust}\right) \end{gathered}$$

(44)

where H(⸱) denotes the Heaviside step function to describe the typical velocity distribution time-history corresponding to the gust loads, this triangular gust lasts 0.5s from t =0 to t = 0.5s, i.e. t_G = 0.05s, w₀ =0.3m/s, ${\tau }_G=V{t}_G/b$.

$$w_G(\tau )=H(\tau )w_0\sin \left(0.5\tau \right)\left(\mathrm{sine\; gust}\right)$$

(45)

where w₀ =0.3m/s, and this sine gust lasts 120 s from t =0 to t = 120s.

Note that the triangular gust is very similar to the traditional 1-cosine gust-type function, both of which can be classified as ephemeral disturbances. Since gusts will result in the lift force and pitch moment change, which then drive the vertical speed and pitch rate change, for simplification, we regard them as the external disturbances by

$$D(t)={\Gamma }_w{w}_G(t)$$

(46)

where${\Gamma }_w$is the effectiveness matrix of gust, which can be approximate by Eqs.(29) and (30),

$${\Gamma }_w={\left[\begin{array}{cc}-\rho {(V+D_w)}^{2}b{s}_p{C}_{L\alpha }/(V{m}_T)& \rho {(V+D_w)}^2{b}^2{s}_p{C}_{m\alpha }/(V{I}_{\alpha })\end{array}\right]}^T$$

(47)

The disturbance parameters of (31) are

$$W_d=\left[\begin{array}{cc}0& 0.5\\ -0.5& 0\end{array}\right],V_d=\left[\begin{array}{cc}2& 0\end{array}\right]$$

(48)

The associated DO gains of (33) are designed as follows

$$l(h,\dot{h})={\left[\begin{array}{cc}20& 6\end{array}\right]}^{T}l(\alpha ,\dot{\alpha })={\left[\begin{array}{cc}5& 0.4\end{array}\right]}^T,\left(\mathrm{for\; sine\; gust}\right)$$

(49)

$$l(h,\dot{h})={\left[\begin{array}{cc}2.5& 0.17\end{array}\right]}^{T}l(\alpha ,\dot{\alpha })={\left[\begin{array}{cc}5& 0.4\end{array}\right]}^T,\left(\mathrm{for\; triangular\; gust}\right)$$

(50)

Remark 2 In this paper single (1-cos) gusts (i.e. triangular gust) and continuous sinusoidal (harmonic) gusts are considered. It can be seen from Eqs (49) and (50) that the observer gains of l (h,$\dot{h}$

) are different for different kinds of gusts, which shows the DO is adaptive for the gusts. Different kinds of gusts can be detected by the short-range light detection and ranging (Lidar) wind measurement.

The parameter of the adaptive observer (A19) is set as γ_d = 0.2. To illustrate our method in trajectory tracking control, a BSMC design from [42] is used for comparison. The sliding mode surface is designed as

$$s={\lambda }_1{z}_1+z_2$$

(51)

and the BSMC controller is designed as follows

$$U_{BSMC}={(g_2)}^{-1}\left(-(K_1+{\lambda }_1){\dot{z}}_1-f_2+{\ddot{x}}_{1d}-\widehat{d}-\left(hs+\varsigma \mathrm{sgn}(s)\right)\right)$$

(52)

And a PID controller is applied as follows

$$U_{\mathrm{PID}}={(g_2)}^{-1}\left(-(K_{p,z_1}{z}_1+K_{I,z_1}{\displaystyle {\int }_0^t{z}_{1}dt})\right.\left.-(K_{p,z_2}{z}_2+K_{I,z_2}{\displaystyle {\int }_0^t{z}_{2}dt})-f_2+{\ddot{x}}_{1d}-\widehat{d}\right)$$

(53)

The parameters of the studied controllers are designed as Table 2 after several design iterations.

By using the proposed DO-BTSMC, BSMC and PID design, the simulation results refer Figure 8, Figure 9, Figure 10, Figure 11, Figure 12, Figure 13 and Figure 14 and Table 3 and Table 4.

Example 1 Firstly, we consider different cases of the BTSMC design with and without learning rate α_r under triangular gust of (44), and the DO is applied, the output responses are shown in Figure 8 and Figure 9 and Table 3, where the reference inputs of root mean square error (RMSE) are h_cmd =0, α_cmd=0.

Example 2 Secondly, we consider the case of the BTSMC design with and without DO under sine gust of (45), and there is no learning rate, the output responses are shown in Figure 10 and Figure 11 and Table 3.

Considering triangular gust input, it can be seen from Figure 8 and Figure 9 and Table 3 that the resulting dynamic responses of the proposed nonlinear aeroelastic system are stable, and the rise-time of plunge responses for BTSMC with learning rate and DO is shorter than that of the BTSMC without learning rate but with DO, which shows the tracking speed of the BTSMC with learning rate for the plunge motions is faster than that of the BTSMC without learning rate, but LR can produce a small overshoots, see Figure 8. The root mean square errors of the plunge and pitch motion for the BTSMC with LR and DO are smaller those for the BTSMC without LR. This shows the learning rate of the BTSMC can improve the control performance.

Now considering sine gust input, Figure 10, Figure 11 and Table3 show that the resulting dynamic responses of the proposed nonlinear aeroelastic system are stable, and the steady-state errors of the BTSMC with DO for the plunge and pitch motions are less than those of the BTSMC without DO. The associated control inputs for the LE control surface are smaller than that of the BTSMC without DO. This shows the disturbance observer of (32) and (A19) can effectively estimate the unknown time-varying gust and the disturbance compensator reduces the disturbance effect for the aeroelastic system.

Finally, we consider the different flutter suppression approaches of the BTSMC, BSMC and PID under triangular gust of (44).By using the proposed design with the DO of (32) and adaptive observer (A19), the output responses are shown in Table 4 and Figure 12, Figure 13 and Figure 14. The sign function of the SMC design is replaced by the tangent hyperbolic function with a boundary layer (where φ_s =0.4),

It can be seen from Figure 12. and Table 4 that the resulting dynamic responses of the proposed nonlinear aeroelastic system are convergent to the desired values, and the rise-time of plunge responses for DO-BTSMC is shorter than those of the DO-BSMC and DO-PID control, which shows the BTSMC response is the fastest among them. The output responses of the PID control easily oscillate under triangular gust while the responses of the BTSMC and BSMC design are steady after transition process, which shows the robust performances of the BTSMC and BSMC. The root mean square errors of the plunge and pitch motion for the DO-BTSMC are the smallest, while those for the PID control are the highest, and RMSEs for the DO-BSMC are in the middle of them, which show tracking precision of the DO-BTSMC is the highest, the second one is the DO-BSMC, the last one is the PID control. Figure 10 shows that the triangular gust of (44) can be effectively estimated and the unobservable disturbances including model uncertainties are estimated by the adaptive observer (A19), see Figure 14, which reduce effect of the triangular gust on the aeroelastic system.

Case II: for Unsteady aerodynamic load

As the post-flutter behavior of the airfoil is intrinsically unsteady, so the unsteady aerodynamic loads are calculated by (10)-(27), and the parameters of the studied controllers are designed as Table 5 after several design iterations.

The associated DO gains of (33) are designed as follows

$$\begin{gathered} l(h,\dot{h})={\left[\begin{array}{cc}25& 0.25\end{array}\right]}^{T}l(\alpha ,\dot{\alpha })={\left[\begin{array}{cc}50& 0.4\end{array}\right]}^T,l(\beta ,\dot{\beta })={\left[\begin{array}{cc}0.5& 0.04\end{array}\right]}^T,l(\gamma ,\dot{\gamma })={\left[\begin{array}{cc}0.5& 0.04\end{array}\right]}^T, \\ \left(\mathrm{for\; sine\; gust}\right) \end{gathered}$$

(54)

$$\begin{gathered} l(h,\dot{h})={\left[\begin{array}{cc}25& 0.25\end{array}\right]}^{T}l(\alpha ,\dot{\alpha })={\left[\begin{array}{cc}50& 0.4\end{array}\right]}^T,l(\beta ,\dot{\beta })={\left[\begin{array}{cc}5& 0.4\end{array}\right]}^T,l(\gamma ,\dot{\gamma })={\left[\begin{array}{cc}5& 0.4\end{array}\right]}^T, \\ \left(\mathrm{for\; triangular\; gust}\right) \end{gathered}$$

(55)

Example 3 Firstly, we consider different cases for the BTSMC design with and without learning rate α_r under triangular gust of (44), and the external disturbances by(46), where${\Gamma }_w$is the gust effectiveness matrix, which can be approximate as

$${\Gamma }_w={\left[\begin{array}{cccc}-\rho {(V+D_w)}^{2}b{C}_{L\alpha }/(V{m}_T)& M_D& M_D/2& M_D/2\end{array}\right]}^T$$

(56)

where$M_D=\rho {(V+D_w)}^2{b}^2{s}_p{C}_{m\alpha }/(V{I}_{\alpha })$, the disturbance parameters are the same as Example 1 with w₀ =3m/s, and the DO is applied, the output responses are shown in Figure 15, Figure 16 and Figure 17.

The control inputs refer to Figure 18.

The wind disturbances are estimated by the proposed DO (32) and refer to Figure 19.

Example 4 Secondly, we consider the case of the BTSMC design with and without DO under sine gust of (45), the disturbance parameters are the same as Example 2, and there is no learning rate, the output responses are shown in Figure 20, Figure 21, Figure 22 and Figure 23 and Table 6.

The control inputs refer to Figure 23.

Considering triangular gust input, it can be seen from Figure 15, Figure 16 and Figure 17 and Table 6 that the resulting dynamic responses of the proposed nonlinear aeroelastic system are stable, and the rise-time of plunge responses for BTSMC with LR and DO is shorter than that of the BTSMC without LR but with DO, which shows the tracking speed of the BTSMC with learning rate and DO for the plunge motions are faster than that of the BTSMC without learning rate, but LR can produce a small overshoot, see Figure 15, Figure 16 and Figure 17. The RMSE of the plunge and pitch motion for the BTSMC with learning rate and DO are smaller than those for the BTSMC without LR. It shows the learning rate of the BTSMC can improve the control performances.

Now considering sine gust input, Figure 20, Figure 21, Figure 22 and Figure 23 and Table3 show that the resulting dynamic responses of the proposed nonlinear aeroelastic system are stable, and the steady-state errors of the BTSMC with DO for the plunge and pitch motions are less than those of the BTSMC without DO, and the associated control inputs are smaller than those of the BTSMC without DO. It shows the disturbance observer of (32) and (A19) can effectively estimate the unknown time-varying gust and the disturbance compensator reduces the disturbance effect for the aeroelastic system. Figure 19 shows that the external disturbances can be estimated quickly and accurately.

4.3. Robust Flutter Suppression with Parameter Uncertainties and Disturbances

Case I: for quasi-steady aerodynamic loads

Example 5 Considering parameter uncertainties of the aeroelastic system, we set them as follows according to parameters perturbation when the work conditions change, uncertainties of k_h, k_α, and c_h are ±1%, uncertainty of c_α and C_Lα are ±30% and ±7% respectively. Meanwhile the external disturbance is set as a sine gust of (45). The controller parameters are designed and listed in Table 2, and the sign function of the SMC design is replaced by a continuous approximation of sign function $sgn\left(s\right)\approx s/\sqrt{s^{T}s+{\beta }_s}$with β_s =1. Since the DO-BTSMC can deal with model uncertainties through adaptive disturbance estimation and gain regulation, then we use the proposed DO-BTSMC design, and the simulation results are shown in Table 7 and Figure 24, Figure 25, Figure 26 and Figure 27.

It can be seen from Figure 24 and Table 7 that the resulting dynamic responses of the proposed nonlinear aeroelastic system are convergent, and the rise-time of plunge responses for DO-BTSMC is the shortest, which shows the BTSMC response is the fastest among them, and demonstrates the DO-BTSMC has fast transient response and finite time convergence. But BTSMC will produce larger overshoots than those of BSMC and PID design, this is because more inputs are required for finite time convergence in the transition process for the DO-BTSMC. The RMSEs of the plunge and pitch motion show tracking precision of the DO-BTSMC are the best, and then is for the DO-BSMC, the last one is for the PID control. Meanwhile, the output responses of the PID control are vibrating under the sine gust, which shows that the robust performance of the PID controller is less than that of the BTSMC and BSMC. Figure 25 shows that the sine gust of (45) can be effectively estimated and the unobservable disturbances are also estimated by the adaptive observer of (A19) as Figure 26, which reduce the effect of the sine gust on the aeroelastic system. Figure 27 shows that the control surface deflections of DO-BTSMC are the smallest and fastest, but the control inputs of the DO-PID control are the largest and easily drive the control surface saturation, which also demonstrates the BTSMC robust performances.

Case II: for Unsteady aerodynamic load

Example 6 For unsteady aerodynamic loads with disturbances and model parameter uncertainties, the parameters of disturbances and uncertainties are the same as Example 5, we use the proposed DO-BTSMC design. The simulation results are shown in Table 8 and Figure 28, Figure 29, Figure 30, Figure 31, Figure 32 and Figure 33.

It can be seen from Figure 28–30 and Table 8 that the resulting dynamic responses of the proposed nonlinear aeroelastic system are convergent, and the rise-time of plunge responses for DO-BTSMC is the shortest, which shows the BTSMC response is the fastest among them, and illustrates the DO-BTSMC has fast transient response and finite time convergence. But BTSMC will produce larger overshoots than BSMC and PID design. This is because more inputs are required for finite time convergence in the transition process for the DO-BTSMC. The root mean square errors of the plunge and pitch motion for the DO-BTSMC are 0.0010m and 0.0183rad respectively, and those for the DO-BSMC are 0.0010m and 0.0151rad respectively, while those for the PID control are 0.0017m and 0.0199rad respectively, which show tracking precision of the DO-BSMC is the best, then is for the DO-BTSMC, the last one is for the PID control. Although the tracking precision of the DO-BSMC is better than that for the DO-BTSMC, the rise-time of plunge motion response for DO-BSMC is very slow, which shows its vulnerable damping.

The disturbance can be estimated by the proposed DO (32) and (A19), the outputs refer to Figure 31 and Figure 32.

The control inputs refer to Figure 33.

Figure 31 shows that the sine gust of (45) can be effectively estimated and the unobservable disturbances are also estimated by the adaptive observer of (A19) as Figure32, which reduces the effect of the sine gust on the aeroelastic system. Figure 33 shows that the control surface deflections of DO-BTSMC are the smallest and fastest. And the control inputs of the DO-BSMC control are the largest and easily drive the control surface saturation, which also demonstrates the DO-BSMC has small damping and reduce its tracking performances.

5. Conclusion

A nonlinear backstepping nonsingular terminal SMC method is used to design a controller to suppress limit cycle oscillation (LCO) and flutter. The developed controller stabilizes the plunge and pitch motions of the two dimensional airfoil with LE and TE control surfaces within finite-time convergence. Stability analysis shows that the closed-loop aeroelastic dynamics are globally asymptotically stable. Six examples of the variable bounded disturbances of triangular gust and sine gust and model parameter uncertainties are simulated. Compared with the BSMC and PID controllers, the DO-BTSMC with learning rate achieves better flutter suppression performances in fast response (improved by 18.6%) and root mean square error( reducing 75% for α). Therefore, the effectiveness and availability of the DO- BTSMC design are demonstrated.