Active Flutter Suppression of Stochastic Airfoil with Uncertain Pitch Stiffness

Abstract

Flutter suppression of a typical airfoil with uncertain pitch stiffness by feedback control of its trailing-edge surface is analyzed. The pitch stiffness coefficient is modeled by a quadratic function of a bounded random variable with the bounded probability density function, thus the airfoil system becomes a stochastic structure depending on its random parameter. It is found that just like its deterministic counterpart the flutter suppression of stochastic airfoil system can also be realized by feedback control. In analysis the feedback controlled stochastic airfoil system is first reduced into its equivalent deterministic system by Gegenbauer polynomial approximation. The critical flutter speed of the stochastic airfoil system with control can be determined by Hopf-bifurcation point of its equivalent system. Then the problem of suppressing limit cycle oscillation of stochastic airfoil system by feedback control can be studied numerically via the equivalent deterministic system. Numerical results show that the critical flutter speeds of a stochastic airfoil system can be reasonably lifted by properly control its trailing-edge surface. And this kind of feedback control is robust.

Keywords

Stochastic airfoil, Flutter suppression, Feedback control, Robust control, Hopf-bifurcation, Limit cycle oscillation, λ-PDF, Gegenbauer polynomial approximation

Nomenclature

ξη: Bounded Random Variable; p_λ(ξ): Probability Density Function; λ: Parameter; τ(λ): Gamma Function; $G_n^{\lambda }(\xi )$: Gegenbauer Polynomials; (λ)_k: Pochhammer Symbol; _h: Plunge Displacement; θ: Pitch Angle; t: Time; m: Mass about the Elastic Axis per Unit Span; I_θ: Moment of Inertia about the Elastic Axis per Unit Span; x_θ: Non-Dimensional Distance Between Center of Gravity and the Elastic Axis; c_h,c_θ: Viscous Damping Coefficients; k_h: Linear Plunge Stiffness; k_θ(θ): Nonlinear Pitch Stiffness; b: Semi-Chord Length of the Airfoil; L: Aerodynamic Lift; M: Aerodynamic Moment; ρ: Air Density; U: Speed; $c_{l_{\theta }}$: Lift Coefficient per Unit Angle of Attack; $c_{m_{\theta }}$: Lift Moment per Unit Angle of Attack; $c_{l_{\phi }}$: Lift Coefficient per Unit Angle of Attack of the Control Surface; $c_{m_{\phi }}$: Lift Moment per Unit Angle of Attack of the Control Surface; ab: Non-dimensional Distance Between the Elastic Axis and the Middle Point of Chord; ϕ: Angle of Attack; δ₁, δ₂: Control Parameters; k_θ k'_θ: Random Parameters.

Introduction

Aircraft wing flutter is a kind of self-excited oscillation resulting from interactions of inertia, elastic, and aerodynamic forces on the structure [1]. Wing flutter may deteriorate structural integrity and even lead to failure of wings. Thus, the prediction of flutter onset and the active suppression of flutter are of the utmost importance in aircraft design.

In the past few decades as a main topic of nonlinear aero-elasticity the studies on mechanism of flutter and various strategies for suppressing flutter have been paid on a great amount of attention. An active flutter suppression strategy for a nonlinear airfoil was examined by Block and Strganac [2], which effectively extends the stable flight region. Based on a back-stepping design technique, Xing and Singh suggested an adaptive controller for the control of an aero-elastic system by output feedback [3]. An analytical and experimental study of a control surface with distributed piezoelectric actuators for active flutter suppression of a wing model was reported by Chen, et al. [4], whose test results show that compared with the open loop system, the close loop system increases its critical flutter speed by 12%. Analytical and experimental studies on active flutter suppression of a wing section by both the leading and the trailing surfaces were reported by Platanitis and Strganac [5], where uncertainty in the nonlinear pitch stiffness was examined. A semi-active flutter suppression strategy for wing section through trailing edge surface control was suggested by Sun, Wen and Hu [6] by using an actuator equipped with Magneto-Rheological (MR) damper. Their simulation results show that the on-off control of MR damper can increase the critical flutter speed by 17% or so.

However, most of the above-mentioned works are limited to deterministic airfoil model only. In fact, some parameters of an airfoil system may be uncertain, for example the pitch stiffness. For the sake of simplification and emphasis on the main effect of the deterministic structure, these uncertainties used to be neglected in analysis. Yet it is more reasonable to model an uncertain physical parameter as a random variable with a given probability density function (λ-PDF for short). Thus, the system considered becomes a stochastic system, whose characteristics depend upon the random variable. And the dynamical behavior in a nonlinear stochastic system is much more complicated than in the corresponding deterministic system. One has to deal with the problem of stochastic bifurcation and even stochastic chaos.

A current mathematical tool to solve the dynamical problem of a stochastic system is the orthogonal polynomial approximation, which is forwarded by Spanos and Ghanem [7], and further developed later on by Li [8]. In view of stochastic structures with bounded random parameters Fang, et al. first used Chebyshev polynomial approximation to solve the evolutionary random response problem of a stochastic system with random variables of an arch-like PDF [9]. Later on Gegenbauer polynomial approximation was suggested by Fang, et al. to solve the dynamical problem of stochastic systems with bounded random variables of λ-PDF [10-15], including the problem of stochastic bifurcation and stochastic chaos and its control or synchronization.

In this paper Gegenbauer polynomial approximation is used to analyze the active flutter suppression of an airfoil system with bounded random parameter of λ-PDF. The stochastic airfoil system is first transformed into its equivalent deterministic system. And the critical flutter speed of the stochastic airfoil system with control can be determined by the Hopf-bifurcation point of its equivalent system. Then the active suppression of limit cycle oscillation of the stochastic airfoil system can be studied through numerical simulations of the equivalent system. Simulation results show that the critical flutter speed of a stochastic airfoil system can be reasonably lifted by proper control of its trailing-edge surface.

λ-PDF and Gegenbauer Polynomial Functions

The probability density functions of a broad type of bounded random variable ξ can be modeled by the so-called λ-PDF, which can be expressed as follows

$$p_{\lambda }(\xi )=\left\{\begin{array}{cc}{q}_{\lambda }{(1-{\xi }^2)}^{\lambda -\frac{1}{2}}& \left|\xi \right|\le 1\\ 0& \left|\xi \right|>1\end{array}\right\}\text{ (1)}$$

where λ ≥ 0 is a fixed parameter, and ρ_λ is a normalizing constant expressed by

$${\text{ρ}}_{\text{λ}}\text{ = }\frac{\text{Γ(λ+1)}}{\text{Γ(}\frac{\text{1}}{\text{2}}\text{)Γ(λ+}\frac{\text{1}}{\text{2}}\text{)}}\text{}\text{}\text{}\text{}\text{}\text{}\text{}(2)$$

where τ(λ) is a Gamma function. The graphs of λ-PDF for different values of λ can be shown as Figure 1.

One can see that when λ = 0, λ-PDF is a pitfall-like distribution

$$p_0(\xi )=\frac{1}{\pi \sqrt{1-{\xi }^2}}\xi \in [-1,1](3)$$

When λ = 0.5, λ-PDF is a uniform distribution

$$p_{1/2}(\xi )\text{ = }1/2\xi \in [-1,1](4)$$

And when λ = 1, λ-PDF is an arch-like distribution

$$p_1(\xi )\text{ = }\frac{2}{\pi }\sqrt{1-{\xi }^2}\xi \in [-1,1](5)$$

One can also see that any bounded random parameter symmetrically defined on [-1, 1], with a mono-peak (or mono-valley) might be modeled approximately by ξ with one of the λ-PDF's most close to it. In addition even some un-symmetrically bounded random parameters can also be approximately modelled by certain polynomial functions of ξ [15].

It is well-known that in orthogonal polynomial approximation the choice of orthogonal basis depends on the probability density function of the random parameter considered. For example, Hermite polynomials for Gaussian random parameter, Legendre polynomials for bounded uniformly distributed random parameter, and so on. As for random parameters with λ-PDF, Gegenbauer polynomials are the unique right choice for orthogonal basis. The first few Gegenbauer polynomials are as follows.

$$\begin{array}{l}{G}_0^{\lambda }\left(\xi \right)\text{ = 1}\\ G_1^{\lambda }\left(\xi \right)\text{ = 2}\lambda \xi \\ G_2^{\lambda }\left(\xi \right)\text{ = 2}\lambda \left(1\text{ + }\lambda \right){\xi }^2\text{ - }\lambda \\ G_3^{\lambda }\left(\xi \right)\text{ = }\frac{4}{3}\lambda \left(1\text{ + }\lambda \right)\left(\text{2 + }\lambda \right){\xi }^3\text{ - 2}\lambda \left(1\text{ + }\lambda \right)\xi \\ G_3^{\lambda }\left(\xi \right)\text{ = }\frac{2}{3}\lambda \left(1\text{ + }\lambda \right)\left(\text{2 + }\lambda \right)\left(\text{3 + }\lambda \right){\xi }^4\text{ - 2}\lambda \left(1\text{ + }\lambda \right)\left(\text{2 + }\lambda \right){\xi }^2\text{ + }\frac{1}{2}\lambda \left(1\text{ + }\lambda \right)(6)\end{array}$$

The orthogonal relationships of Gegenbauer polynomials may be expressed as

$$\underset{-1}{\overset{1}{{\displaystyle \int }}}{p}_{\lambda }(\xi )G_l^{\lambda }(\xi )G_n^{\lambda }(\xi )d\xi =\left\{\begin{array}{c}{b}_n^{\lambda }\\ 0\end{array}\begin{array}{c}l=n\\ l\ne n\end{array}\right\}\text{ (7)}$$

Where

$$b_n^{\lambda }\text{ = }\frac{\Gamma (1/2)}{2^{2\lambda -1}n!}\cdot \frac{\lambda \Gamma (2\lambda +n)}{(\lambda +n)\Gamma (\lambda )\Gamma (\lambda +1/2)}(8)$$

The first and the second order recurrent formulas for Gegenbauer polynomials are expressed as [15]

$$\xi G_n^{\lambda }(\xi )\text{ = }{\alpha }_n^{\lambda }{G}_{n-1}^{\lambda }(\xi )+{\beta }_n^{\lambda }{G}_{n+1}^{\lambda }(\xi )(9)$$

and

$${\xi }^2{G}_n^{\lambda }(\xi )\text{ = }{\gamma }_n^{\lambda }{G}_{n-2}^{\lambda }(\xi )+{\rho }_n^{\lambda }{G}_n^{\lambda }(\xi )+{\sigma }_n^{\lambda }{G}_{n+2}^{\lambda }(\xi )(10)$$

where

$$\begin{array}{l}{\alpha }_n^{\lambda }=\frac{2\lambda +n-1}{2(\lambda +n)}, {\beta }_n^{\lambda }=\frac{n+1}{2(\lambda +n)}\\ {\gamma }_n^{\lambda }=\frac{(2\lambda +n-1)(2\lambda +n-2)}{4(\lambda +n)(\lambda +n-1)}\text{, }{\sigma }_n^{\lambda }\text{ = }\frac{(n+1)(n+2)}{4(\lambda +n)(\lambda +n+1)}\\ {\rho }_n^{\lambda }\text{ = }\frac{1}{4(\lambda +n)}\left\{\frac{(n+1)(2\lambda +n)}{\lambda +n+1}+\frac{n(2\lambda +n-1)}{\lambda +n-1}\right\}(11)\end{array}$$

Equation (7) implies a kind of weighted averaging property with λ-PDF as a weighting function. Owing to the orthogonal relationships of Gegenbauer polynomials, any f(ξ) ⊂ L², defined on [-1, 1] can be expressed as a series of $G_n^{\lambda }(\xi )$, namely

$$f(\xi )\text{ = }{\displaystyle \sum _{n=0}^N{s}_n}{G}_n^{\lambda }(\xi )(12)$$

where

$$s_n=\frac{2^{2\lambda -1}(\lambda +n)n!{\Gamma }^2(\lambda )}{\pi \Gamma (2\lambda +n)}{\displaystyle {\int }_{-1}^1{(1-{\xi }^2)}^{\lambda -1/2}}f(\xi )G_n^{\lambda }(\xi )d\xi (13)$$

It is worth noting here that Eq. (12) is valid only for taking the sum over infinite number of items. If in practice only a limited number of items are taken in this equation, the result is merely an approximation with a minimal mean square residual.

Dynamic Equation of a Stochastic Airfoil System

A nonlinear stochastic airfoil system permitting plunge h and pitch θ motion about the elastic axis, and equipped with a trailing-edge control surface is shown in Figure 2. The dynamical equation of the system, neglecting control surface dynamics, may be described as

$$\begin{array}{l}m\ddot{h}+m{x}_{\theta }b\ddot{\theta }+c_h\dot{h}+k_{h}h\text{ = }-L\\ m{x}_{\theta }b\ddot{h}+I_{\theta }\ddot{\theta }+c_{\theta }\dot{\theta }+k_{\theta }(\theta )\theta \text{ = }M(14)\end{array}$$

where h and θ stand for system responses; an overhead dot on h or θ denotes d/dt , and t is the time ; m and I_θ are mass and moment of inertia about the elastic axis per unit span, respectively; x_θ is non-dimensional distance between center of gravity and the elastic axis; c_h and c_θ are viscous damping coefficients for plunge and pitch motion respectively; k_h and k_θ(θ) are linear plunge stiffness and nonlinear pitch stiffness respectively; b is the semi-chord length of the airfoil; and L and M are quasi-steady aerodynamic lift and moment, which can be approximately expressed as

$$\begin{array}{l}L\text{ = }\rho U^{2}b{c}_{l_{\theta }}(\theta +\frac{\dot{h}}{U}+(\frac{1}{2}-a)b\frac{\dot{\theta }}{U})+\rho U^{2}b{c}_{l_{\phi }}\phi \\ M\text{ = }\rho U^2{b}^2{c}_{m_{\theta }}(\theta +\frac{\dot{h}}{U}+(\frac{1}{2}-a)b\frac{\dot{\theta }}{U})+\rho U^2{b}^2{c}_{m_{\phi }}\phi (15)\end{array}$$

where ρ and U represent air density and the speed of coming incompressible flow respectively; $c_{l_{\theta }}$, $c_{m_{\theta }}$ are coefficient of lift and moment per unit angle of attack of the airfoil; $c_{l_{\theta }}$, $c_{m_{\theta }}$ are coefficient of lift and moment per unit angle of attack of the control surface; ab is the non-dimensional distance between the elastic axis and the middle point of chord; and ϕ is the angle of attack of the control surface to carry out the feedback control law.

It is well-known that for a conventional airfoil system the state feedback control of its trailing edge surface is an effective way to suppress airfoil flutter. To illustrate this flutter suppression strategy applies as well to a stochastic airfoil system, we simply use the proportional feedback control law, namely,

$$\phi \text{ = }-{\delta }_1\theta +{\delta }_{2}h(16)$$

Where δ₁, δ₂ are independent control parameters.

In general the coefficient of nonlinear pitch stiffness k_θ(θ) is difficult to obtain analytically. One could only find its approximate expression through curve fitting the experimental data. Now we assume k_θ(θ) may be expressed in the form

$$k_{\theta }(\theta )\text{ = }{k}_{\theta }+{k^{\prime }}_{\theta }{\theta }^2(17)$$

Where k_θ and k'_θ are to be determined yet. Since there is always some uncertainty in measuring and manufacturing, we would rather take k_θ and k'_θ as random parameters. Without losing generality, we assume that they can be modeled by a quadratic function of a bounded random variable ξ with a given λ-PDF, namely

$$\begin{array}{l}{k}_{\theta }=K_{\theta }(1+{\sigma }_1\xi +{\sigma }_2{\xi }^2)\\ {k^{\prime }}_{\theta }={K^{\prime }}_{\theta }(1+{{\sigma }^{\prime }}_1\xi +{{\sigma }^{\prime }}_2{\xi }^2)(18)\end{array}$$

Where k_θ, k'_θ, σ₁, σ₂, σ'₁, σ'₂ are properly selected deterministic constants. For different values of these constants, k_θ, k'_θ may have different probability distributions. For example, let η = k_θ/k_θ, σ₁ = 0.2, and σ₂ = 0.1, then the probability density function of η can be shown as Figure 3. One can see that p_λ(η) even can be un-symmetric. Substituting equations (15) - (18) into Eq. (14), we have the following dynamical equation for the stochastic airfoil system

$$\begin{array}{l}{m}_{11}\ddot{h}+m_{12}\ddot{\theta }+c_{11}\dot{h}+c_{12}\dot{\theta }+k_{11}h+k_{12}\theta \text{ = }0\\ m_{21}\ddot{h}+m_{22}\ddot{\theta }+c_{21}\dot{h}+c_{22}\dot{\theta }+k_{21}h+k_{22}\theta +k_{\theta }\theta +{k^{\prime }}_{\theta }{\theta }^3\text{ = }0(19)\end{array}$$

Where $m_{11}\text{ = }m\text{, }{m}_{12}\text{ = }m{x}_{\theta }b,m_{21}\text{ = }m{x}_{\theta }b,m{}_{22}\text{= }{I}_{\theta },c_{11}\text{ = }{c}_h+\rho Ub{c}_{l_{\theta }},c_{12}\text{ = }\rho U{b}^2{c}_{l_{\theta }}(\frac{1}{2}-a)$,

$c_{21}\text{ = }-\rho U{b}^2{c}_{m_{\theta }},c_{22}\text{ = }{c}_{\theta }-\rho U{b}^3{c}_{m_{\theta }}(\frac{1}{2}-a),k_{11}=k_h+{\delta }_2\rho U^{2}b{c}_{l_{\phi }}$,

$k_{12}\text{ = }\rho U^{2}b{c}_{l_{\theta }}-{\delta }_1\rho U^{2}b{c}_{l_{\phi }},k_{21}\text{ = }-{\delta }_2\rho U^2{b}^2{c}_{m_{\phi }},k_{22}\text{ = }{\delta }_1\rho U^2{b}^2{c}_{m_{\phi }}-\rho U^2{b}^2{c}_{m_{\theta }}$

This is what we need the dynamic equation of a 2-DOF nonlinear stochastic airfoil system with feedback control, which can be further processed by the Gegenbauer polynomial approximation.

Gagenbauer Polynomial Approximation for a Stochastic Airfoil System

Since the system (19) itself depending upon a bounded random parameter ξ, the response of the system should be a function of both t and ξ, namely

$$\begin{array}{l}h\text{ = }h(t,\xi )\\ \theta \text{ = }\theta (t,\xi )(20)\end{array}$$

In Gegenbauer polynomial approximation these responses can be approximated by

$$\begin{array}{l}h(t,\xi )\text{ = }{\displaystyle \sum _{l=0}^N{h}_l}(t)G_l^{\lambda }(\xi )\\ \theta (t,\xi )\text{ = }{\displaystyle \sum _{l=0}^N{\theta }_l}(t)G_l^{\lambda }(\xi )(21)\end{array}$$

where $G_l^{\lambda }(\xi )$ represents the l^th order Gegenbauer polynomial, and N stands for the highest order we have taken. When N → ∞, ${\displaystyle {\sum }_{l=0}^N{h}_l}(t)G_l^{\lambda }(\xi )$ and ${\displaystyle {\sum }_{l=0}^N{\theta }_l}(t)G_l^{\lambda }(\xi )$ are exact solutions of Eq. (19), otherwise they are just approximate solutions. In practice we can choose a proper N to meet the required accuracy.

Substituting Eq. (21) into Eq. (19), we have

$$\begin{array}{l}{\displaystyle \sum _{l=0}^N[m_{11}{\ddot{h}}_l(t)+m_{12}\ddot{\theta }{}_l(t)+c_{11}{\dot{h}}_l(t)+c_{12}{\dot{\theta }}_l(t)+k_{11}{h}_l(t)+k_{12}{\theta }_l(t)]G_l^{\lambda }}(\xi )\text{ = }0\\ {\displaystyle \sum _{l=0}^N[m_{21}{\ddot{h}}_l(t)+m_{22}{\ddot{\theta }}_l(t)+c_{21}{\dot{h}}_l(t)+c_{22}{\dot{\theta }}_l(t)+k_{21}{h}_l(t)+k_{22}}{\theta }_l(t)]G_l^{\lambda }(\xi )+\\ +k_{\theta }{\displaystyle \sum _{l=0}^N{\theta }_l(t)G_l^{\lambda }}(\xi )+{k^{\prime }}_{\theta }{[}^{{\displaystyle \sum _{l=0}^N{\theta }_l(t)G_l^{\lambda }}}\text{ = }0(22)\end{array}$$

By using the recurrent formulas of Gegenbauer polynomials, Eq. (9) - Eq. (10), we can eliminate the random variables explicitly appearing in the second equation of Eq. (22). Thus, the linear term $k_{\theta }{\displaystyle {\sum }_{l=0}^N{\theta }_l}(t)G_l^{\lambda }(\xi )$ can be rewritten as

$$\begin{array}{l}{k}_{\theta }{\displaystyle \sum _{l=0}^N{\theta }_l(t)G_l^{\lambda }}(\xi )\text{ = }{K}_{\theta }(1+{\sigma }_1\xi +{\sigma }_2{\xi }^2){\displaystyle \sum _{l=0}^N{\theta }_l(t)G_l^{\lambda }}(\xi )\\ \text{ = }{K}_{\theta }{\displaystyle \sum _{l=0}^N\{{\theta }_l}(t)+\sigma {}_1[_{{\theta }_{l+1}}^{(}{}_t{}^{)}{}_{{\alpha }_{l+1}^{\lambda }}{}^{+}{}_{{\theta }_{l-1}}{}^{(}{}_t{}^{)}{}_{{\beta }_{l-1}^{\lambda }}{}^{]}\\ +{\sigma }_2[{\gamma }_{l+2}^{\lambda }{\theta }_{l+2}(t)+{\rho }_l^{\lambda }{\theta }_l(t)+{\upsilon }_{l-2}^{\lambda }{\theta }_{l-2}(t)]\}{G}_l^{\lambda }(\xi )(23)\end{array}$$

While deriving Eq. (23), it is assumed that θ_l(t) = 0，when l > N or l < 0, which is also valid in the following derivation.

The nonlinear term $[{\displaystyle {\sum }_{l=0}^N{\theta }_l(t)G_l^{\lambda }(\xi ){]}^3}$ in Eq. (22) can be expanded as

$${[{\displaystyle \sum _{l=0}^N{x}_l(t)G_l^{\lambda }}(\xi )]}^3\text{ = }{\displaystyle \sum _{i=0}^{3N}{\displaystyle \sum _{k=0}^i{\displaystyle \sum _{j=0}^k{x}_j}}}(t)x_{k-j}(t)x_{i-k}(t)G_j^{\lambda }(\xi )G_{k-j}^{\lambda }(\xi )G_{i-k}^{\lambda }(\xi )(24)$$

By Eq. (6) the product $G_i^{\lambda }(\xi )G_j^{\lambda }(\xi )G_l^{\lambda }(\xi )$ can be rewritten as

$$G_i^{\lambda }(\xi )G_j^{\lambda }(\xi )G_l^{\lambda }(\xi )\text{ = }{\displaystyle \sum _{m=0}^{i+j+l}{c}_{i,j,l,m}^{\lambda }}{(\frac{\xi -1}{2})}^m(25)$$

Where

$$c_{i,j,l,m}^{\lambda }\text{ = }{\displaystyle \sum _{n=0}^m{\displaystyle \sum _{k=0}^n{a}_{i,k}^{\lambda }}}{a}_{j,n-k}^{\lambda }{a}_{l,m-n}^{\lambda }(26)$$

Where $a_{i,k}^{\lambda }$ can be written as $a_{i,k}^{\lambda }\text{ = }\frac{1}{k!(i-k)!}\frac{{(2\lambda )}_i{(2\lambda +i)}_k}{{(\lambda +1/2)}_k}$ in which (λ)_k is the Pochhammer symbol, which stands for ${(\lambda )}_k\text{ = }\lambda (\lambda +1)\cdots (\lambda +k-1)$

It is worth noting that in Eq. (26) we have $a_{n,m}^{\lambda }=0$, when m > n. By Eq. (6) the term [(ξ - 1)/2]^m can be expressed as a linear combination of $G_i^{\lambda }(\xi )$ as

$${(\frac{\xi -1}{2})}^m\text{ = }{\displaystyle \sum _{i=0}^m{b}_{m,i}^{\lambda }}{G}_i^{\lambda }(\xi )(27)$$

where

$$\begin{array}{l}{b}_{i,i}^{\lambda }=\frac{1}{a_{i,i}}\\ b_{i+1,i}^{\lambda }\text{ = }-\frac{a_{i+1,i}}{a_{i+1,i+1}{a}_{i,i}}\\ b_{i+2,i}^{\lambda }\text{ = }-\frac{a_{i+2,i}{a}_{i+1,i+1}-a_{i+1,i}{a}_{i+2,i+1}}{a_{i,i}{a}_{i+1,i+1}{a}_{i+2,i+2}}\\ b_{i+3,i}^{\lambda }\text{ = }-\frac{a_{i+3,i}a{}_{i+2,i+2}a{}_{i+1,i+1}-a_{i+3,i+1}{a}_{i+2,i+2}{a}_{i+1,i}-a_{i+3,i+2}{a}_{i+2,i}{a}_{+1,i+1}+a_{i+3,i+2}{a}_{i+2,i+1}{a}_{i+1,i}}{a_{i,i}{a}_{i+1,i+1}{a}_{i+2,i+2}{a}_{i+3,i+3}}(28)\end{array}$$

By substituting Eq. (27) into Eq. (25), the product $G_i^{\lambda }(\xi )G_j^{\lambda }(\xi )G_l^{\lambda }(\xi )$ can be expressed as linear combinations of individual Gegenbauer polynomials, namely

$$G_i^{\lambda }(\xi )G_j^{\lambda }(\xi )G_l^{\lambda }(\xi )\text{ = }{\displaystyle \sum _{n=0}^{i+j+l}{d}_{i,j,l,n}^{\lambda }{G}_n^{\lambda }(\xi )}(29)$$

where

$$d_{i,j,l,n}^{\lambda }\text{ = }{\displaystyle \sum _{m=n}^{i+j+l}{c}_{i,j,l,m}}{b}_{m,n}^{\lambda }(30)$$

Substituting Eq. (29) into Eq. (24), we have

$$\begin{array}{l}{[}^{{\displaystyle \sum _{l=0}^N{x}_l(t)G_l^{\lambda }}}\text{ = }{\displaystyle \sum _{i=0}^{3N}{\displaystyle \sum _{k=0}^i{\displaystyle \sum _{j=0}^k{x}_j}}}(t)x_{k-j}(t)x_{i-k}(t){\displaystyle \sum _{n=0}^i{d}_{j,k-j,i-k,n}^{\lambda }}{G}_n^{\lambda }(\xi )\\ \text{ = }{\displaystyle \sum _{n=0}^{3N}{g}_n^{\lambda }(t)}{G}_n^{\lambda }(\xi )(31)\end{array}$$

Where

$$g_n^{\lambda }(t)\text{ = }{\displaystyle \sum _{i=n}^{3N}{\displaystyle \sum _{k=0}^i{\displaystyle \sum _{j=0}^k{d}_{j,k-j,i-k,n}^{\lambda }{x}_j}}(t)x_{k-j}(t)x_{i-k}(t)}(32)$$

Similarly, by the recurrent formulas Eq. (9) - Eq. (10), the nonlinear term ${k^{\prime }}_{\theta }{[{\displaystyle {\sum }_{l=0}^N{\theta }_l(t)}{G}_l^{\lambda }(\xi )]}^3$ can be reduced into

$$\begin{array}{l}{k^{\prime }}_{\theta }{[}^{{\displaystyle \sum _{l=0}^N{\theta }_l(t)G_l^{\lambda }}}\text{ = }{K^{\prime }}_{\theta }(1+{{\sigma }^{\prime }}_1\xi +{{\sigma }^{\prime }}_2{\xi }^2){[}^{{\displaystyle \sum _{l=0}^N{\theta }_l(t)G_l^{\lambda }}}\\ \text{ = }{K^{\prime }}_{\theta }(1+{{\sigma }^{\prime }}_1\xi +{{\sigma }^{\prime }}_2{\xi }^2){\displaystyle \sum _{l=0}^{3N}{g}_l^{\lambda }}(t)G_l^{\lambda }(\xi )\\ \text{ = }{K^{\prime }}_{\theta }{\displaystyle \sum _{l=0}^{3N}\{g_l^{\lambda }}(t)+{{\sigma }^{\prime }}_1[g_{l+1}^{\lambda }(t){\alpha }_{l+1}^{\lambda }+g_{l-1}^{\lambda }(t){\beta }_{l-1}^{\lambda }]\\ +{{\sigma }^{\prime }}_2[{\gamma }_{l+2}^{\lambda }{g}_{l+2}^{\lambda }(t)+{\rho }_l^{\lambda }{g}_l^{\lambda }(t)+{\upsilon }_{l-2}^{\lambda }{g}_{l-2}^{\lambda }(t)]\}{G}_l^{\lambda }(\xi )(33)\end{array}$$

Substituting Eq. (23) and Eq. (33) into Eq. (22), we have

$$\begin{array}{l}{\displaystyle \sum _{l=0}^N[m_{11}{\ddot{h}}_l(t)+m_{12}\ddot{\theta }{}_l(t)+c_{11}{\dot{h}}_l(t)+c_{12}{\dot{\theta }}_l(t)+k_{11}{h}_l(t)+k_{12}{\theta }_l(t)]G_l^{\lambda }}(\xi )\text{ = }0\\ {\displaystyle \sum _{l=0}^N[m_{21}{\ddot{h}}_l(t)+m_{22}{\ddot{\theta }}_l(t)+c_{21}{\dot{h}}_l(t)+c_{22}{\dot{\theta }}_l(t)+k_{21}{h}_l(t)+k_{22}}{\theta }_l(t)]G_l^{\lambda }(\xi )+K_{\theta }{\displaystyle \sum _{l=0}^N\{{\theta }_l}(t)+\\ \sigma {}_1[_{{\theta }_{l+1}}^{(}{}_t{}^{)}{}_{{\alpha }_{l+1}^{\lambda }}{}^{+}{}_{{\theta }_{l-1}}{}^{(}{}_t{}^{)}{}_{{\beta }_{l-1}^{\lambda }}{}^{]}+{\sigma }_2[{\gamma }_{l+2}^{\lambda }{\theta }_{l+2}(t)+{\rho }_l^{\lambda }\theta {}_l(_t^{)}+{\upsilon }_{l-2}^{\lambda }\theta {}_{l-2}(_t^{)}]\}{G}_l^{\lambda }(\xi )+{K^{\prime }}_{\theta }{\displaystyle \sum _{l=0}^{3N}\{g_l^{\lambda }}(t)\\ +{{\sigma }^{\prime }}_1[g_{l+1}^{\lambda }(t){\alpha }_{l+1}^{\lambda }+g_{l-1}^{\lambda }(t){\beta }_{l-1}^{\lambda }]+{{\sigma }^{\prime }}_2[{\gamma }_{l+2}^{\lambda }{g}_{l+2}^{\lambda }(t)+{\rho }_l^{\lambda }g{}_l{}^{\lambda }(_t^{)}+{\upsilon }_{l-2}^{\lambda }{g}_{l-2}^{\lambda }(t)]\}{G}_l^{\lambda }(\xi )\text{ = }0(34)\end{array}$$

Multiplying $G_i^{\lambda }(\xi )$, (i = 0,1,….,N) onto both sides of Eq. (34) in sequence, then taking expectations with respect to ξ, owing to the orthogonality relationships of Gegenbauer polynomials, we finally obtain

$$\left\{\begin{array}{l}{m}_{11}{\ddot{h}}_0(t)+m_{12}\ddot{\theta }{}_0(_t^{)}+c_{11}{\dot{h}}_0(t)+c_{12}{\dot{\theta }}_0(t)+k_{11}{h}_0(t)+k_{12}{\theta }_0(t)\text{ = }0\\ m_{21}{\ddot{h}}_0(t)+m_{22}{\ddot{\theta }}_0(t)+c_{21}{\dot{h}}_0(t)+c_{22}{\dot{\theta }}_0(t)+k_{21}{h}_0(t)+k_{22}{\theta }_0(t)\\ +K_{\theta }\{{\theta }_0(t)+\sigma {}_1\theta {}_1(t){\alpha }_1^{\lambda }+{\sigma }_2[{\gamma }_2^{\lambda }\theta {}_2(_t^{)}+{\rho }_0^{\lambda }{\theta }_0(t)]\}+\\ {K^{\prime }}_{\theta }\{g_0^{\lambda }(t)+{{\sigma }^{\prime }}_1{g}_1^{\lambda }(t){\alpha }_1^{\lambda }+{{\sigma }^{\prime }}_2[{\gamma }_2^{\lambda }{g}_2^{\lambda }(t)+{\rho }_0^{\lambda }g{}_0{}^{\lambda }(_t^{)}]\}\text{ = }0\\ m_{11}{\ddot{h}}_1(t)+m_{12}\ddot{\theta }{}_1(_t^{)}+c_{11}{\dot{h}}_1(t)+c_{12}{\dot{\theta }}_1(t)+k_{11}{h}_1(t)+k_{12}{\theta }_1(t)\text{ = }0\\ m_{21}{\ddot{h}}_1(t)+m_{22}{\ddot{\theta }}_1(t)+c_{21}{\dot{h}}_1(t)+c_{22}{\dot{\theta }}_1(t)+k_{21}{h}_1(t)+k_{22}{\theta }_1(t)\\ +K_{\theta }\{{\theta }_1(t)+\sigma {}_1[_{{\theta }_2}^{(}{}_t{}^{)}{}_{{\alpha }_2^{\lambda }}{}^{+}{}_{{\theta }_0}{}^{(}{}_t{}^{)}{}_{{\beta }_0^{\lambda }}{}^{]}+{\sigma }_2[{\gamma }_3^{\lambda }\theta {}_3(_t^{)}+{\rho }_1^{\lambda }{\theta }_1(t)]\}\\ {K^{\prime }}_{\theta }\{g_1^{\lambda }(t)+{{\sigma }^{\prime }}_1[g_2^{\lambda }(t){\alpha }_2^{\lambda }+g_0^{\lambda }(t){\beta }_0^{\lambda }]+{{\sigma }^{\prime }}_2[{\gamma }_3^{\lambda }{g}_3^{\lambda }(t)+{\rho }_1^{\lambda }g{}_1{}^{\lambda }(_t^{)}]\}\text{ = }0\\ \cdots \cdots \\ m_{11}{\ddot{h}}_N(t)+m_{12}\ddot{\theta }{}_N(_t^{)}+c_{11}{\dot{h}}_N(t)+c_{12}{\dot{\theta }}_N(t)+k_{11}{h}_N(t)+k_{12}{\theta }_N(t)\text{ = }0\\ m_{21}{\ddot{h}}_N(t)+m_{22}{\ddot{\theta }}_N(t)+c_{21}{\dot{h}}_N(t)+c_{22}{\dot{\theta }}_N(t)+k_{21}{h}_N(t)+k_{22}{\theta }_N(t)+K_{\theta }\{{\theta }_N(t)+\\ \sigma {}_1[_{{\theta }_{N+1}}^{(}{}_t{}^{)}{}_{{\alpha }_{N+1}^{\lambda }}{}^{+}{}_{{\theta }_{N-1}}{}^{(}{}_t{}^{)}{}_{{\beta }_{N-1}^{\lambda }}{}^{]}+{\sigma }_2[{\gamma }_{N+2}^{\lambda }\theta {}_{N+2}(_t^{)}+{\rho }_N^{\lambda }{\theta }_N(t)+{\upsilon }_{N-2}^{\lambda }{\theta }_{N-2}(t)]\}\\ +{K^{\prime }}_{\theta }\{g_N^{\lambda }(t)+{{\sigma }^{\prime }}_1[g_{N+1}^{\lambda }(t){\alpha }_{N+1}^{\lambda }+g_{N-1}^{\lambda }(t){\beta }_{N-1}^{\lambda }]+{{\sigma }^{\prime }}_2[{\gamma }_{N+2}^{\lambda }{g}_{N+2}^{\lambda }(t)+{\rho }_N^{\lambda }{g}_N^{\lambda }(t)\\ +{\upsilon }_{N-2}^{\lambda }{g}_{N-2}^{\lambda }(t)]\}\text{ = }0\end{array}\right\}\text{ (35) }$$

Now Eq. (35) is a deterministic nonlinear one obtained through some weighted averaging by Gegenbauer polynomial approximation on the stochastic airfoil system with feedback control. We call it an equivalent deterministic system, which plays a significant role in the following analysis. Since Eq. (35) is deterministic, it can be solved by any available effective conventional theory or method, including Runge-Kutta method, MATLAB, MAPLE and so on. Once h_i(t), θ_i(t) （i = 0~N） are obtained through Eq. (35), the approximate stochastic responses of the stochastic airfoil system can be readily obtained through Eq. (21), namely

$$\begin{array}{l}h(t,\xi )\approx {\displaystyle \sum _{l=0}^N{h}_l}(t)G_l^{\lambda }(\xi )\\ \theta (t,\xi )\approx {\displaystyle \sum _{l=0}^N{\theta }_l}(t)G_l^{\lambda }(\xi )(36)\end{array}$$

And the ensemble averaging responses are simply as

$$\begin{array}{l}E[h(t,\xi )]\approx {\displaystyle \sum _{l=0}^N{h}_l}(t)E[G_l^{\lambda }(\xi )]\text{ = }{h}_0(t)\\ E[\theta (t,\xi )]\approx {\displaystyle \sum _{l=0}^N{\theta }_l}(t)E[G_l^{\lambda }(\xi )]\text{ = }{\theta }_0(t)(37)\end{array}$$

For any sample random variable $\xi =\overline{\xi }$, the sample responses of the stochastic system are as

$$\begin{array}{l}h(t,\overline{\xi })\approx {\displaystyle \sum _{l=0}^N{h}_l}(t)G_l^{\lambda }(\overline{\xi })\\ \theta (t,\overline{\xi })\approx {\displaystyle \sum _{l=0}^N{\theta }_l}(t)G_l^{\lambda }(\overline{\xi })(38)\end{array}$$

Critical Flutter Speed of Stochastic Airfoil System

By introducing the following state variables

$$Y(t)\text{ = }{[y_1,y_2,y_3,y_4,\cdots ,y_{2N+1},y_{2N+2},y_{2N+3},y_{2N+4},\cdots ,y_{4N+3},y_{4N+4}]}^T$$

$$\text{ = }{[h_0,{\theta }_0,h_1,{\theta }_1,\cdots ,h_N,{\theta }_N,{\dot{h}}_0,{\dot{\theta }}_0,\dots ,{\dot{h}}_N,{\dot{\theta }}_N]}^T$$

the equivalent deterministic system can be rewritten into a system of first order differential equations

where

$$\begin{array}{l}{\kappa }_1\text{ = }(c_{21}{m}_{12}-c_{11}{m}_{22})/M,{\kappa }_2\text{ = }(c_{22}{m}_{12}-m_{22}{c}_{12})/M,\kappa {}_3\text{= }(k_{21}{m}_{12}-k_{11}{m}_{22})/M,\\ {\kappa }_4\text{ = }(m_{12}{k}_{22}+m_{12}{K}_{\theta }-m_{22}{k}_{12})/M,{\kappa }_5\text{ = }{m}_{12}{K}_{\theta }{\sigma }_2/M,{\kappa }_6\text{ = }{m}_{12}{K}_{\theta }\sigma {}_1/M\\ {\kappa }_7\text{ = }{m}_{12}{K}_{\theta }\sigma {}_1/M,{\kappa }_8\text{ = }{m}_{12}{K}_{\theta }{\sigma }_2/M,{\kappa }_9\text{ = }{m}_{12}{K}_{\theta }{\sigma }_2/M,{\kappa }_{10}\text{ = }{m}_{12}/M\\ \eta {}_1\text{= }-(m_{11}{c}_{21}-m_{21}{c}_{11})/M,{\eta }_2\text{ = }-(m_{11}{c}_{22}-c_{12}{m}_{21})/M,{\eta }_3\text{ = }-(m_{11}{k}_{21}-k_{11}{m}_{21})/M\\ {\eta }_4\text{ = }-(m_{11}{k}_{22}+m_{11}{K}_{\theta }-m_{21}{k}_{12})/M,{\eta }_5\text{ = }-m_{11}{\sigma }_2{K}_{\theta }/M,{\eta }_6\text{ = }-m_{11}{K}_{\theta }\sigma {}_1/M\\ {\eta }_7\text{ = }-m_{11}{K}_{\theta }\sigma {}_1/M,\eta {}_8\text{= }-m_{11}{K}_{\theta }{\sigma }_2/M,{\eta }_9\text{ = }-m_{11}{K}_{\theta }{\sigma }_2/M,{\eta }_{10}\text{ = }-m_{11}/M\end{array}$$

$$\begin{array}{l}{f}_i(y_2,\cdots y_{2N+2})\text{ = }{K}_{\theta }\{g_i^{\lambda }(t)+\sigma {}_1[_{g_{i+1}^{\lambda }}^{(}{}_t{}^{)}{}_{{\alpha }_{i+1}^{\lambda }}{}^{+}{}_{g_{i-1}^{\lambda }}{}^{(}{}_t{}^{)}{}_{{\beta }_{i-1}^{\lambda }}{}^{]}\\ +{\sigma }_2[{\gamma }_{i+2}^{\lambda }{g}_{i+2}^{\lambda }(t)+{\rho }_i^{\lambda }g{}_i{}^{\lambda }(_t^{)}+{\upsilon }_{i-2}^{\lambda }{g}_{i-2}^{\lambda }(t)]\}\end{array}$$

$$M=m_{12}{m}_{21}-m_{11}{m}_{22}$$

Eq. (39) can be written in short as

$$\dot{Y}\text{ = }F(Y,U),Y\in R^{4N+4}(40)$$

The critical flutter speed can be determined by the Hopf-bifurcation point of Eq. (40), which is the transition point of equilibrium solution into self-excited oscillatory solution [16]. The Hopf-bifurcation point of Eq. (40) can be determined as follows.

Let the equilibrium solution of Eq. (40) be Y_E, then the transient solution in its neighborhood may be assumed as

$$Y\text{ = }{Y}_E+\epsilon P{e}^{i\omega t}(41)$$

where ε is an infinitesimal, and P = P₁ + iP₂ is a complex eigenvector with respect to the imaginary eigenvalue iω. Substituting Eq. (41) into Eq. (40), and expanding it into a Taylor series around Y_E, we have

$$i\omega \epsilon P{e}^{i\omega t}\text{ = }F(Y_E,U)+\epsilon F_Y(Y_E,U)P{e}^{i\omega t}+0({\epsilon }^2)(42)$$

At the equilibrium point Y_E, we have

$$F(Y_E,U)\text{ = }0(43)$$

Substituting Eq. (43) into Eq. (42), and neglecting the second and higher order terms, we have

$$i\omega P-F_Y(Y_E,U)P\text{ = }0(44)$$

Then inserting P = P₁ + iP₂ into Eq. (44) and letting the imaginary and the real parts of the equation equal zero, we have

$$F_Y{P}_1+\omega P_2\text{ = }0,F_Y{P}_2-\omega P_1\text{ = }0(45)$$

To guarantee Eq. (45) having unique solutions P₁ and P₂, we need the following additional restriction equations

$$q^T{P}_1\text{ = }0,q^T{P}_2\text{ = }1(46)$$

where q is a constant vector [16]. Combining Eq. (43), Eq. (45) and Eq. (46) into one, we have the following equation

$$\Gamma (Z)=\left[\begin{array}{l}{F}_Y{P}_1+\omega P_2\\ q^T{P}_1\\ F_Y{P}_2-\omega P_1\\ q^T{P}_2-1\\ F(Y_E,U)\end{array}\right]\text{ = 0 }Z\text{ = }\left[\begin{array}{l}P{}_1\\ U\\ P_2\\ \omega \\ Y\end{array}\right](47)$$

Eq. (47) can be solved by the Newtonian iterative method.

Numerical Results and Discussions

The effectiveness of the proposed flutter suppression strategy for a stochastic airfoil system can be shown through illustrative numerical examples. In calculation the deterministic data of airfoil system are taken as follows [17]

a = -0.6, b = 0.135 m, k_h = 284404 N/m, c_h = 27.43 Ns/m, $c_{\theta }=0.036Ns$, $\rho =1.225kg/m^3$, $c_{l_{\theta }}=6.28$, $c_{l_{\phi }}=3.358$, $c_{m_{\theta }}=(0.5+a)c_{l_{\theta }}$, $m=12.387kg$, $I_{\theta }=0.065kg{m}^2$, $x_{\theta }=[0.0873-(b+ab)]/b$