Robust aeroelastic instability suppression of an advanced wing with model uncertainty in subsonic compressible flow field

Aerospace Science and Technology 25 (2013) 242–252

Contents lists available at SciVerse ScienceDirect

Aerospace Science and Technology www.elsevier.com/locate/aescte

Robust aeroelastic instability suppression of an advanced wing with model uncertainty in subsonic compressible flow field Ji-Seok Song a , Jeonghwan Choo a , Seog-Ju Cha a , Sungsoo Na a,∗ , Zhanming Qin b a b

Korea University, Seoul 136-701, Republic of Korea Xi’an Jiaotong University, Xi’an, 710049, People’s Republic of China

a r t i c l e

i n f o

a b s t r a c t

Article history: Received 21 August 2010 Received in revised form 28 September 2011 Accepted 15 January 2012 Available online 20 January 2012 Keywords: Aeroelastic instability Model uncertainty Subsonic compressible flow Sliding mode control

This paper addresses the robust aeroelastic instability suppression of an advanced aircraft wing modeled as a thin-walled beam structure with fiber-reinforced composite material featuring a circumferentially asymmetric stiffness (CAS) configuration. The unsteady aerodynamic loads in subsonic compressible flow are derived through the indicial function approach. Aeroelastic instability suppression is achieved via sliding mode control (SMC) based on sliding mode observer (SMO). To demonstrate the robustness of the SMC based on SMO, linear quadratic Gaussian (LQG) methodology is compared with respect to the model uncertainty. To this end, the obtained numerical simulation results emphasize the efficiency of the sliding mode control based on sliding mode observer to control the unstable aeroelastic phenomenon in conjunction with the model uncertainty. © 2012 Elsevier Masson SAS. All rights reserved.

1. Introduction The dynamical aeroelastic instability phenomenon is caused by the interactions among aerodynamic loads, static aeroelasticity, structural elasticity and dynamics. As a self-excited dynamic instability, it may decrease performance of the flight vehicles or even lead to fatigue failure of the structure. To overcome these problems, many research works on the flutter suppression have been conducted in various ways [18,21,32]. Among these investigations, advanced thin-walled beam structures, which are made of anisotropic composite material, are used in the design of new generation flight vehicle’s wing structure because of their high structural efficiency. In the previous research by Umansky [36], the theory of thin-walled beam structure has dealt with a new part of strength of materials and airplanes. The directional characteristic in the fibrous constituent materials and static aeroelastic tailoring is studied by Krone [10]. And Rehfield and Atilgan [31] developed the elastic coupling mechanism in composite thin-walled beams as rotor blades. Furthermore, the vibration of the thin-walled beam structures was studied by Librescu and Song [16] and Song and Librescu [35] and non-classical effects such as transverse shear strain effect, warping constraint, three-dimensional strain effect and anisotropy of the thin-walled beam structure have been addressed by Librescu et al. [13,15]. Qin et al. [30] developed a unified aeroelastic model

*

Corresponding author. E-mail address: [email protected] (S. Na).

1270-9638/$ – see front matter doi:10.1016/j.ast.2012.01.012

© 2012 Elsevier Masson

SAS. All rights reserved.

and presented subcritical aeroelastic response in the compressible subsonic flight speed range. The related research works regarding active flutter suppression strategy of the aircraft wing have been developed to overcome the deficiencies of the passive flutter suppression methodology [26]. In particular, piezoelectric actuators are applied to improve the aerodynamic and aeroelastic performance. Heeg [8] and McGowan [7] analytically and experimentally examined the capabilities of piezoelectric plate actuators for flutter suppression performance. Librescu et al. [14] addressed piezoelectric actuators spread over the entire wing span which efficiently suppressed the occurrence of the dynamic aeroelastic instability by LQG regulator based on SMO. Also, thermally induced dynamic response and active control methodology of the rotating composite thin-walled beam structure are investigated by Na et al. [25]. For the two-dimensional airfoil, the linear quadratic Gaussian (LQG) regulator based on sliding mode observer (SMO) technique efficiently alleviates the unstable aeroelastic response in incompressible and compressible flow fields [22–24]. Degaki and Suzuki [3] adopted sliding mode control (SMC) for suppressing incompressible flow flutter instability in which aerodynamic uncertainties and structural nonlinearity are incorporated. SMC based on SMO are used to expand the flutter boundary in supersonic flow with existence of uncertainty by Lee et al. [11]. The control of three-dimensional flight vehicle’s wing is of practical importance, however, research works on the control of the dynamic aeroelastic instability of flight vehicles in subsonic compressible flows are scarce. In this connection, we present an advanced aircraft wing which is modeled as anisotropic thin-walled beam structure in

J.-S. Song et al. / Aerospace Science and Technology 25 (2013) 242–252

243

Nomenclature

(x, y , z) Global Cartesian coordinate system (n, y , s) Local coordinate system u 0 , v 0 , w 0 Three translations in x, z, y directions, respectively θx , φ, θz Rotations of the cross-section with respect to x, z, y axes, respectively F w (s), n · a(s) Primary warping and secondary warping functions ε0y y , εny y On-contour and off-contour axial normal strains, respectively γsy0 Contour-wise shear strain γ yz , γ yx Transverse shear measures of the cross-section ψ(s) Torsional function θ Ply angle D / Dt Substantial derivative with respect to dimensional time Flight Mach number, U ∞ /a∞ , a∞ is the speed of Mf sound U ∞ , U n Stream-wise and chord-wise free stream speeds, respectively ai j One-dimensional stiffness coefficients N sy , N y y , N ny Stress resultants

Mx, M y

One-dimensional stress couples

ˆx M

Piezoelectrically induced bending moment at the wing tip L ae , L g , L b Unsteady aerodynamic loads, lift due to gust and lift due to blast load, respectively T ae , T g Twist moments about the reference axis induced by the unsteady aerodynamic and gust load ˆ 0 , φˆ w Non-dimensional rigid body translations in z direction and rotation about the reference axis ta , S a Piezoactuator thickness and height, respectively K, L Linear quadratic Gaussian (LQG) regulator gain matrix and Kalman filter gain matrix H l , H nl Linear and nonlinear sliding mode observer gain matrix, respectively ζ (τ ), ω(τ ) Plant disturbance and sensor noise Plant disturbance intensity and sensor noise intensity Z, Ω matrix, respectively Φ, η Nonlinearity scale factor of the SMO and SMC, respectively Tr , Tl Orthogonal and non-singular transformation matrix, respectively

Fig. 1. Configuration of a three-dimensional advanced aircraft wing structure.

subsonic compressible flow and propose robust state estimation methodology via comparison between conventional Kalman filter (KF) and SMO with respect to parametric uncertainty in system matrix. Furthermore, the control performance to suppress or postpone the occurrence of the dynamic aeroelastic instability by using robust estimation control methodology is investigated when aeroelastic system is subjected to various type of blast pulses. 2. Aeroelastic system modeling

coordinate and the (n, s) denote the thickness and tangential-wise coordinate of the cross-section contour line. In order to obtain the mathematical representation of the thinwalled beam structure, the three-dimensional displacement quantities are defined as follows [29]:

u (x, y , z, t ) = u 0 ( y , t ) + z · φ( y , t ) w (x, y , z, t ) = w 0 ( y , t ) − x · φ( y , t )



v (x, y , z, t ) = v 0 ( y , t ) + x(s) − n ·

dz ds

(1b)

 · θz ( y , t )

  dx · θx ( y , t ) + z(s) + n · ds   − F w (s) + n · a(s) · φ  ( y , t )

2.1. Structural modeling We consider an advanced aircraft wing structure consisting of a single-cell, closed cross-section and fiber-reinforced composite thin-walled beam. Fig. 1 presents the geometric configuration of the three-dimensional aircraft wing structure. Two coordinate systems are defined; one is the global Cartesian system (x, y , z) and the other is the local coordinate system (n, y , s) where the (x, z) denote the cross-section of the beam while y-axis is the span-wise

(1a)

(1c)

where

θx ( y , t ) = γ yz ( y , t ) − w 0 ( y , t ) 

θz ( y , t ) = γxy ( y , t ) − u 0 ( y , t )

(2a) (2b)

244

J.-S. Song et al. / Aerospace Science and Technology 25 (2013) 242–252

Fig. 3. Cross-section specification of the wing structure.

Fig. 2. Circumferentially asymmetric stiffness (CAS) configuration.

  dz dx a(s) = − z · +x· ds

(2c)

ds

In Eqs. (1a)–(1c), the kinetic variables, which constitute the basic unknowns of the problem, u 0 ( y , t ), w 0 ( y , t ), v 0 ( y , t ), θx ( y , t ), θz ( y , t ) and φ( y , t ) represent three translations in the x, z, y directions and three rotations of the plane of the cross-section about the axes x, z and the twist about the y-axis. Herein, γ yz ( y , t ) and γxy ( y , t ) denote the transverse shear measures of the crosssection. The strain quantities that contribute to the potential energy can be summarized as:



ˆ0 ˙ˆ 0 − φˆ + 2 ∂ w w a (ˆx, η, τ ) = U n · 2w · tan Λe A R ∂η   1 ∂ φˆ − xˆ · φ˙ˆ + · tan Λe A R ∂η ˆ − xˆ · φ˙ˆ } = U · {w n

aT

0 y y (s,

y, t ) + n · ε

n y y (s,

y, t )

(3a)

Tangential shear strain:

γsy (s, y , t ) = γsy0 (s, y , t ) + ψ(s) · φ  ( y , t ) Transverse shear strain: dz dx γny (s, y , t ) = −γxy + γ yz ds ds  dz    dx   = − u 0 + θz · + w 0 + θx · ds ds

(3b)

(3c)

where ε 0y y (s, y , t ) and εny y (s, y , t ) denote on-contour and offcontour axis normal strains, respectively. A circumferentially asymmetric stiffness (CAS) configuration, shown in Fig. 2, is selected as to induce bending–torsion elastic coupling for thin-walled beam structure. It implies the ply angle distribution [17,25]

θ(x) = −θ(−x)



ˆ0 2 ∂w A R ∂η

 · tan Λe

  1 ∂ φˆ · tan Λe φ˙ˆ a P (η, τ ) ≡ φ˙ˆ + A R ∂η

(6a) (6b)

The indicial functions are determined numerically for certain Mach numbers. These quantities, such as (Φc )0 (τ ); (ΦcM )0 (τ ) and (Φcq )0 (τ ); (ΦcMq )0 (τ ), represent the lift and moment indicial functions due to the unit step change of the vertical velocity and pitching rate at the leading edge, respectively. The aerodynamic lift and moment about the mid-chord (reference axis, see Fig. 3) associated with the vertical velocity are expressed are:

L T (η, τ ) = −2π · ρ∞ · U n2 · b2

·



+



˙

ˆ aT (η, 0) + φˆ a P (η, 0) · (Φc )0 (τ ) w

τ

ˆ aT (η, σ ) + φ˙ˆ a P (η, σ )) ∂( w · (Φc )0 (τ − σ ) dσ ∂σ



0

(7a) T yT (η, τ ) = −4π · ρ∞ ·

·



2.2. Subsonic aerodynamic loads in compressible flow The indicial function approach gives a general and convenient way to describe the compressible unsteady aerodynamics within the linear aerodynamic theory and arbitrary small motion. These indicial functions can be approximated via computational fluid dynamics or experimental results. Therefore, in this paper, the strip theory aerodynamics and the two-dimensional indicial functions

(5)

where the functions of dimensionless span-wise coordinate (η ) and time (t) are defined as [30]:

(4)

Within the CAS configuration, the vibrations of the thin-walled beam structure completely split into two groups, which consist of vertical bending–twist–vertical transverse shear motions and extension–lateral bending–lateral transverse shear motions. The vertical bending (w 0 : plunging) and twist (φ : pitching) motions are mainly concerned in the present study.



aP

˙ˆ 0 − φˆ + ˆ aT (η, τ ) ≡ U n · 2 w w

Span-wise strain:

ε y y (n, s, y , t ) = ε

developed by Leishman are adopted to derive the aeroelastic behavior of the three-dimensional aircraft wing [2,12,19,20]. The downwash velocity in dimensionless form can be expressed as [30]:

·b

˙

2



ˆ aT (η, 0) + φˆ a P (η, 0) · (ΦcM )0 (τ ) w

τ +

U n2

ˆ aT (η, σ ) + φ˙ˆ a P (η, σ )) ∂( w · (ΦcM )0 (τ − σ ) dσ ∂σ



0

+ b · L T (η, τ )

(7b)

In addition, the aerodynamic lift and moment about the midchord associated with the pitching rate can be described as:

J.-S. Song et al. / Aerospace Science and Technology 25 (2013) 242–252

L q (η, τ ) = 4π · ρ∞ ·

  · b · φ˙ˆ a P (η, 0) · (Φcq )0 (τ )

U n2

∂(φ˙ˆ a P (η, σ )) + · (Φcq )0 (τ − σ ) dσ ∂σ 0

˙  2 2 T yq (η, τ ) = 8π · ρ∞ · U n · b · φˆ a P (η, 0) · (ΦcMq )0 (τ )

τ

τ +

∂(φ˙ˆ a P (η, σ )) · (ΦcMq )0 (τ − σ ) dσ ∂σ

0

(8a)



T ae (η, τ ) = T yT (η, τ ) + T yq (η, τ )

(9b)

Herein, the analytical indicial functions that are approximated in quasi-polynomial forms are expressed as:



A ci · exp −βic · τ

t2 (δ T − δ V + δ W nc ) dt = 0 t1

with δ u 0 = δ w 0 = δ v 0 = δθx = δθz = δφ = 0 at t = t 1 , t 2

Kinetic energy: T=



A cM · exp −βicM · τ i

cq

(Φcq )c /2 (τ ) = A 0 +



cq

cq

A i · exp −βi · τ

·



i =1

(ΦcMq )c /2 (τ ) =

cMq A0

+

l

cMq Ai

· exp



cMq −βi

·τ



(10b)

V =

In the aeroelastic design of an advanced aircraft wing structure, the effect of aeroelastic response due to sharp-edged gust and time-dependent external pressure pulses should be considered in a preliminary stage of design. The indicial functions and the Duhamel’s convolution integral are used to describe the sharp-edged gust that is traveling at the same speed. Herein the discrete type gust model is adopted as below [19]:

(11)

where the peak gust velocity, V G , represents the gust intensity, which is identically applied along the whole wing span and H (τ ) is the Heaviside step function. The aerodynamic lift and moment about the reference axis due to the sharp-edged gust are:



L g (τ ) = C L φ · b · U n · w G (0) · ψ K (τ )

+ 0

T g (τ ) =

1 2

∂ w G (τ0 ) · ψ K (τ − τ0 ) dτ0 ∂ τ0 

∂u ∂t



C L φ · b2 · U n · w G (0) · ψ K (τ )

ρ (k)

(k)

2

 +

L  ml 

1 2

(12a)

k=1h

c

0

2.3. Sharp-edged gust and blast loads

w G (τ ) = H (τ ) · V G

k=1h

c

∂w ∂t

2

 +

∂v ∂t

2  dn · ds · dy

(15)

Strain energy:

i =1

τ

2



(10a)

i =1 l

L  ml 

1

0



(14)

where all three quantities are defined as:



i =1 l

(13)

(8b)

(9a)

(ΦcM )c /2 (τ ) = A cM 0 +

τ >0

The extended Hamilton’s principle can be used to obtain the system equations of motion and boundary conditions.

L ae (η, τ ) = L T (η, τ ) + L q (η, τ )

(Φc )c /2 (τ ) = A c0 +

(12b)

2.4. Governing aeroelastic systems

where the subscript q means the pitching rate. The total unsteady aerodynamic lift (positive for upwards) and the moment (positive for nose-up) in subsonic compressible flow about the mid-chord are summarized as follows [28]:

l



where the Kussner’s function, ψ K (τ ) can be developed by following exponential approximation [2]:

ψ K (τ ) = 1 − 0.5e −0.13τ − 0.5e −τ ,

0

+ b · L q (η, τ )

∂ w G (τ0 ) · ψ K (τ − τ0 ) dτ0 ∂ τ0

+

τ

245

[σ y y · ε y y + σsy · γsy

(k)

+ σny · γny ]h(k) dn · ds · dy

(16a)

Virtual work of non-conservative forces:

L δ W nc =





p z ( y , t ) · δ w 0 ( y , t ) + m y ( y , t ) · δφ( y , t ) dy

0

ˆ x δθx ( L , t ) +M

(16b)

in which the total lift density (positive upward), p z = L ae + L g and the total twist moment density (positive for nose-up), m y = T ae + T g . Herein L ae , T ae are defined earlier in Eqs. (9a), (9b), while L g are aerodynamic lift due to gust defined in Eq. (12a); T g is an aerodynamic moment due to gust defined in Eq. (12b). Because of the features for the CAS configuration (see Section 2.1), the total equations of motion and the boundary conditions are completely decoupled and for present bending–torsion coupled problem, the corresponding equilibrium equations are represented as:

δ w 0: δφ :

¨0=0 Q z + L ae + L g − b1 · w

M y − B w + T ae + T g − (b4 + b5 ) · φ¨ + (b10 + b18 ) · φ¨ 

=0 δθx :

(17a)

M x − Q z − (b4 + b14 ) · θ¨x = 0 

(17b) (17c)

246

J.-S. Song et al. / Aerospace Science and Technology 25 (2013) 242–252

Fig. 4. Distribution of piezoactuators.

Boundary conditions: At y = 0;

w 0 = 0, φ = 0, φ  = 0, θ x = 0

At y = L ;

ˆx Q z = 0, B w = 0, M x = M − B w + M y + (b10 + b18 ) · φ¨  = 0

(18)

In addition, the force–displacement relations are described as below:

⎧ M ⎪ ⎪ x ⎪ ⎨Q z

⎪ Bw ⎪ ⎪ ⎩

My

⎫ ⎪ ⎪ ⎪ ⎬

⎡

a33

⎢0 ⎢ =⎢ ⎪ ⎣0 ⎪ ⎪ ⎭

a37

0

0

⎫ ⎤⎧  θx ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎬ ⎨  0 ⎥ ⎥ ( w 0 + θx ) ⎥ ⎪ φ  0 ⎦⎪ ⎪ ⎪ ⎪ ⎪ ⎭ ⎩  a77 φ a37

a55

a56

a56

a66

0

0

(19)





¨ 0 = 0 (20a) a55 · w 0 + θx + a56 · φ  + L ae + L g − b1 · w 



δφ :

  a37 · θx + a77 · φ  − a56 · w  0 + θx − a66 · φ

δθx :

+ T ae + T g − (b4 + b5 ) · φ¨ + (b10 + b18 ) · φ¨  = 0   a33 · θx + a37 · φ  − a55 · w 0 + θx − a56 · φ  − (b4 + b14 ) · θ¨x = 0

(20b)

(20c)

Boundary conditions: At y = 0;

w 0 = 0, φ = 0, φ  = 0, θ x = 0

At y = L ;

δ w 0:

  −a56 · w 0 + θx − a66 · φ  + a37 · θx

δφ  :

+ a77 · φ  = −(b10 + b18 ) · φ¨    −a56 · w 0 + θx − a66 · φ  = 0

δθx :

ˆx a33 · θx + a37 · φ  = M

M x ( y, t ) =

z · N yy + L yy ·

dx

(21)



ds

c

  Q z ( y, t ) =

N sy · c



B w ( y, t ) = −

 M y ( y, t ) =



dz ds

+ Nny ·

ds

(23a)

 ds



c

(22)

(23b)

where x˜ is the estimated state of x, a 2m(1 + 2) × 1 vector, with  and m being three aerodynamic lag terms and five structural modes taken in the present study, respectively, and K and L denote the optimal state-feedback control gain matrix and Kalman filter gain matrix, respectively. Herein the plant disturbance ζ (τ ) and sensor noise ω(τ ) are Gaussian white noise, uncorrelated zeromean random vectors satisfy the following correlation matrices [4].

 E

ds

F w (s) · N y y + a(s) · L y y ds

N sy · ψ(s) ds c

dx



u (τ ) = − K · x˜ (τ )

where b1 , b4 , b5 , b10 , b14 , b15 , b18 denote the one-dimensional inertia coefficients and M x , Q z , B w , M y are the one-dimensional stress resultants and stress couple measures that are defined as:

 



x˙˜ (τ ) = A · x˜ (τ ) + B · u (τ ) + L · y˜ (τ ) − C · x(τ )

with

a55 · w 0 + θx + a56 · φ  = 0

δφ :

Within the legitimacy that all state vectors are not measurable accurately, the state estimator, or observer, should be constructed in order to estimate unmeasurable state variables based on the output measurements and control variables. The feedback control scheme based on the estimated states is described in Fig. 5, which represents the linear quadratic Gaussian (LQG) regulator with a Kalman filter. In this sense, the state-space representation of the closed-loop regulator is expressed as:

y˜ (τ ) = C · x˜ (τ )





3. Robust control methodology 3.1. Linear quadratic Gaussian regulator

The equations of motion and boundary conditions in terms of the six basic unknowns are:

δ w0:

ˆ x denotes the piezoelectrically induced bendFurthermore, M ing moment that appears solely in the boundary conditions at the beam tip. As a result, the control is achieved via the piezoelectrically induced boundary bending moment (see Refs. [16,17,35]). In this context, piezoactuators featuring in-plane isotropic properties, spread over the entire span of the beam and bonded symmetrically on the outer and inner faces of the beam, but activated out-ofphase, generate a bending moment at the beam tip in response to one of the mechanical quantities characterizing the system response (see Fig. 4).

ζ (τ1 )

ω(τ1 )



   Z 0 ζ (τ2 ) ω (τ2 ) = δ(τ1 − τ2 ) 0 Ω



T

T

(24)

Within Eq. (24), the plant disturbance intensity matrix Z and the sensor noise intensity matrix Ω are positive definite. Also, δ denotes the Dirac’s delta function and E {·} represents the expected value, which has zero if its inner two quantities are uncorrelated. The optimal control gain matrix K is obtained by solving the algebraic Riccati equation (ARE) developed to minimize the performance index. In the same manner, the Kalman filter gain matrix L is obtained from the filter algebraic Riccati equation (FARE). These two equations are represented respectively, as follows:

J.-S. Song et al. / Aerospace Science and Technology 25 (2013) 242–252

247

Fig. 5. Block diagram of the LQG regulator system.

Fig. 6. Block diagram of the sliding mode system.

A T · P + P · A + M − P · B · N −1 · B T · P = 0 T

T

A· R + R · A +Z− R ·C ·Ω

−1

·C ·R =0

(25a) (25b)

where P and R are the positive definite solution of the above equations, respectively. Also, the disturbance and noise intensity matrices are dual of the performance index matrices M and N in order. Herein the matrices M and N minimize the performance index expressed as below:

J LQG = E

∞ 

x

T

u

T

    M 0 x 0

N

u

dτ

(26)

The optimal control gain K and Kalman filter gain L matrices are obtained from

K = − N −1 · B T · P

(27a)

−1

(27b)

L = −R · C · Ω





x˙˜ (τ ) = A · x˜ (τ ) + B · u (τ ) + H l · y˜ (τ ) − C · x(τ ) + H nl · σ (28)

0

T

their neglect may cause both control and observation spillover. Spillover effects are undesirable and may cause system instability [1] and reduction in performance [9]. It is desirable to reduce the observation spillover as to remove the potential to generate any instability. To this end, a sliding mode observer is introduced so that it reduces the effect of observation spillover from the unmeasured states. In this sense, a sliding mode observer is designed to provide estimation of the states. Considering the controlled modes of dynamic system, assuming those are observable, and x˜ denoting as an estimate of x [27], the sliding mode observer in the given aeroelastic system has the form (see Fig. 6) [6]:

3.2. State estimation via sliding mode observer In the present study the controller and the observer are designed based on the mathematical model that incorporates both the measured and the estimated states only. Since unmeasured states are not considered in the controller and observer design,

where H l and H nl denote the linear and nonlinear observer gain matrix, respectively, while the discontinuous switching vector σ is defined by:



σ=

−Φ ·

F ∗ ·C ·e  F ∗ ·C ·e 

0

if e = 0 otherwise

(29)

In Eq. (29), Φ is a positive scalar function as a nonlinearity scale factor of the observer, e is error between real state and estimated state, and the symmetric positive definite matrix F ∗ is a solution of the follow Lyapunov equation:



F ∗ · Ap + Ap

T

· F ∗ = − Qˆ

(30)

248

J.-S. Song et al. / Aerospace Science and Technology 25 (2013) 242–252

where Qˆ is a positive definite matrix which is design parameter, taken as an identity matrix in this study. In addition, a stable design matrix A p comes from the linearly transformed dynamical systems shown in Ref. [6]. A sliding motion, which converges toward the equilibrium points in the error space, takes place on the hyper-plane given by:

Se = e ∈ Rn: C · e = 0

T By considering Q¯ 21 = Q¯ 12 , performance index can be expressed ∗ in terms of y -coordinate system.

J SMC =

1 2

∞ 

 y 1∗ T · Q¯ ∗ · y ∗1 + υ T · Q¯ 22 · υ dτ

(39a)

τf

with

!

(31)

in which n means dimension of present aeroelastic system matrix. The corresponding error dynamics of the system is described as:

e˙ (τ ) = [ A − H l · C ] · e (τ ) + H nl · σ − B e · Q g (τ )

(32)

where Q g denotes the generalized gust vector with dimension m × 1, B e is a (2 + ) · m × m matrix, with  and m being three aerodynamic lag terms and five structural modes taken in the present study, respectively. 3.3. Robust sliding mode control

−1 ¯ Q¯ ∗ = Q¯ 11 − Q¯ 12 · Q¯ 22 · Q 21

(39b)

−1 ¯ υ = y ∗2 + Q¯ 22 · Q 21 · y ∗1

(39c)

The positive definite solution P ∗ is obtained from the algebraic Riccati equation

¯ + A¯ T · P ∗ − P ∗ · A 12 · Q¯ −1 · A T · P ∗ + Q¯ ∗ = 0 P∗ · A 12 22

(40)

¯ = A 11 − A 12 · Q˜ −1 · Q¯ 21 , which is called the modified where A 22 constraint equation, is developed by eliminating y ∗2 (τ ) among Eqs. (35a) and (39c). Then the optimal υ minimizing the performance index, Eq. (39a), is given by:

The design process of this control method is composed of two parts: (1) the design of the hyper-plane S, so that the reducedorder sliding motion appears when the states are confined to the hyper-plane which is selected by the designer; (2) the design of the control law which enforces the states to remain on the hyperplane in finite time, τ f . Consider the given aeroelastic system as follows:

The feedback matrix can be obtained by substituting Eq. (41) into Eq. (39c) as follows:

x˙ (τ ) = A · x(τ ) + B · u (τ ) + B e · Q g (τ )

To design the control law, a second transformation is defined by [5]:

(33)

In order to design the proper hyper-plane S, state equation (28) is transformed into a regular form by [5,6]:

y ∗ (τ ) = T r · x(τ )

(34)

where T r is an orthogonal matrix computed by Q R decomposition and it confirms T r · B = [0 B 2 ] T . The transformed state equations are expressed as:

y˙ ∗1 (τ ) = A 11 · y ∗1 (τ ) + A 12 · y ∗2 (τ )

(35a)

y˙ ∗2 (τ ) = A 21 · y ∗1 (τ ) + A 22 · y ∗2 (τ ) + B 2 · u (τ )

(35b)

In the same manner, the switching function can be partitioned as S · T r−1 = [ S 1 S 2 ] and it can be re-written as S 1 · y ∗1 (τ ) + S 2 · y ∗2 (τ ) = 0 in y ∗-coordinate. By letting S 2 = I κ , the hyper-plane matrix S will be determined by:



S= M

Iκ



· Tr

(36)

where M is the feedback matrix in Eq. (35b) and the subscription κ represents the dimension of the control vectors, i.e. u (τ ) ∈ R q . This procedure has the advantage of decreasing the possibility of numerical errors by reducing the amount of calculation in designing the hyper-plane. To find the feedback matrix M, the optimal sliding mode method, which minimizes the quadratic performance index, is applied. The quadratic performance index is expressed as follows [33,34]:

J SMC =

1 2

∞

x T (τ ) · Q¯ · x(τ ) dτ

(37)

τf

where Q¯ is a transformed symmetric definite matrix, defined as:



T r · Q¯ · T rT =

Q¯ 11

Q¯ 12

Q¯ 21

Q¯ 22



(38)

−1 T υ = − Q¯ 22 · A 12 · P ∗ · y ∗1

(41)

 −1  T M = Q¯ 22 · A 12 · P ∗ + Q¯ 21

(42)

z∗ (τ ) = T l · y ∗ (τ )

(43)

where T l is the orthogonal and non-singular matrix and the associated transformed state equations are expressed as:

ˆ 11 · z∗ (τ ) + A 12 · z∗ (τ ) z˙ 1∗ (τ ) = A 1 2

(44a)

ˆ 21 · z (τ ) + Aˆ 22 · z (τ ) + B 2 · u (τ ) z2 (τ ) = A 1 2

(44b)

˙∗

∗

∗

ˆ 11 = A 11 − A 12 · M, Aˆ 21 = M · Aˆ 11 + A 21 − A 22 · M, and where A ˆ 22 = M · A 12 + A 22 . A To this end, the linear and nonlinear parts of the control law in the z∗ -coordinate are, respectively

 





1 ∗ ∗ ˆ ˆ · z2∗ ul z∗ = − B − 2 · A 21 · z1 + A 22 − Φ

F¯ · z2∗

 

1 unl z∗ = −η · B − 2 ·

 F¯ · z2∗ 

,

!

(45a)

z2 = 0

(45b)

where Φ is any n × n matrix with left-hand half-plane eigenvalues, η is a positive scalar parameter to be selected and F¯ denote the positive definite unique solution of Lyapunov equation

F¯ · Φ ∗ + Φ ∗ T · F¯ = − I q

(46)

The sliding mode control law is defined as the summation of the linear and nonlinear part of the control law in original coordinate system.

u (x) = ul (x) · x(τ ) + η ·

W · x(τ ) Y · x(τ )

(47)

where



1 ˆ ul (x) = − B − 2 · A 21

−1



W = −B2 · 0

F¯





ˆ 22 − Φ ∗ · T l · T r A · Tl · Tr ,



Y= 0

F¯



(48a)

· Tl · Tr

(48b)

J.-S. Song et al. / Aerospace Science and Technology 25 (2013) 242–252

249

Table 1 Material properties of the composite material (graphite/epoxy). Parameter

Value

E 11 E 22 = E 33 G 12 G 13 = G 23

206.8 × 109 [N/m2 ] 5.17 × 109 [N/m2 ] 3.10 × 109 [N/m2 ] 2.55 × 109 [N/m2 ] 0.25 1.528 × 108 [kg/m3 ]

μ12 = μ13 = μ23 ρ

Table 2 Geometric specifications of the wing structure. Parameter

Value

Length (L) Width (2b) Depth (2d) Aspect ratio Ply angle Sweep angle Wall thickness (h) Number of layers

4.5407 [m] 0.757 [m] 0.0997 [m] 12 105◦ (default), 115◦ 0◦ (default), −5◦ , 30◦ 0.0203 [m] 6

4. Simulation results and discussions 4.1. Dynamic aeroelastic response The dynamical aeroelastic responses of advanced aircraft wing structure are revealed in this section. The wing structure is excited by the initial disturbances which are given such that only the first state has a value of 0.01 and the other states are all zeros. The response of the vertical bending and the twist motions are mainly concerned in the present study. Furthermore, the flutter speed has been verified herein via the solution of both complex eigenvalue problem and from the subcritical aeroelastic response analysis. From the result (not shown here), the aeroelastic response amplitude increased markedly when the flight speed reached M f = 0.7. At this moment, the wing structure becomes aeroelastically unstable that is referred to as flutter phenomenon, whereas the aeroelastic response converges to the equilibrium state at the other flight speeds [30]. The material properties and geometric specifications of the wing structure are listed in Tables 1 and 2, respectively.

Fig. 7. Non-dimensional plunging displacement of the wing structure with Kalman filter (a) and sliding mode observer (b) due to initial disturbance (θ = 105◦ , M f = 0.7).

4.2. Performance of robust estimation In Figs. 7 and 8, the performance of the Kalman filter (KF) and sliding mode observer (SMO) is compared, in which the available sensor output is the first mode measurement only due to possible failures of some sensors. The estimation performance of the sliding mode observer is superior to Kalman filter, for both plunging displacement and pitching displacement. However, it should be noted that initial larger fluctuation shown in Fig. 8(b), which is called chattering, is observed for the case of SMO compared with Kalman filter case since chattering is caused by the discontinuous switching vector in Eq. (29). In this regard, as for chattering alleviation it may be implemented using sigmoid-like function. Figs. 9 and 10 represent the estimation of non-dimensional plunging and pitching displacements, respectively. Through a comparison of plunging and pitching displacements in Figs. 7 and 8 with nominal case, it is found that the parametric uncertainty among the system matrix of the original plant affects the estimation performance. Herein the parametric uncertainty is produced by 10% reduction of the mass-density in the estimator, while the plant system has the nominal mass-density. From Figs. 9 and 10, sliding mode observer generates stable plunging and pitching deflection estimates in spite of the parametric uncertainty.

Fig. 8. Non-dimensional pitching displacement of the wing structure with Kalman filter (a) and sliding mode observer (b) due to initial disturbance (θ = 105◦ , M f = 0.7).

In Figs. 11(a) and 12(a), the previous results are shown again in terms of the dimensionless plunging and pitching displacement errors in the nominal system. It is shown that the error produced by sliding mode observer rapidly converges to zero, while Kalman filter fails. Furthermore, the similar conclusion is reached in Figs. 11(b) and 12(b) in a parametric uncertainty with 10% mass-density reduction.

250

J.-S. Song et al. / Aerospace Science and Technology 25 (2013) 242–252

Fig. 9. Non-dimensional plunging displacement of the wing structure with Kalman filter (a) and sliding mode observer (b) due to initial disturbance considering a 10% mass-density reduction of the observer system (θ = 105◦ , M f = 0.7).

Fig. 10. Non-dimensional pitching displacement of the wing structure with Kalman filter (a) and sliding mode observer (b) due to initial disturbance considering a 10% mass-density reduction of the observer system (θ = 105◦ , M f = 0.7).

4.3. Robust vibration control based on observer From the results shown in previous section, the sliding mode observer generates stable plunging and pitching estimates. To demonstrate its outstanding performance, when it is connected with the controller, various simulations are performed in the subcritical flight speeds.

Fig. 11. Non-dimensional plunging error between actual and estimated displacements with nominal (a) and 10% mass-density reduction (b) systems (θ = 105◦ , M f = 0.7).

Fig. 12. Non-dimensional pitching error between actual and estimated displacements with nominal (a) and 10% mass-density reduction (b) systems (θ = 105◦ , M f = 0.7).

Fig. 13 depicts the parametric uncertainty effect on the dynamic response control based on observer. When there exists a model uncertainty with respect to 10% mass-density reduction, either LQG or SMC based on KF cannot afford to control the vertical bending vibration since the KF does not produce correct estimates as shown as in Figs. 9(a) and 10(a). The simulation results in Fig. 14, in contrast to Fig. 13, reveal that sliding mode control based on sliding mode observer is ro-

J.-S. Song et al. / Aerospace Science and Technology 25 (2013) 242–252

Fig. 13. Uncontrolled and controlled response using LQG and SMC controllers based on Kalman filter with nominal (a) and 10% mass-density reduction (b) systems subjected to sharp-edged gust (θ = 105◦ , M f = 0.7).

Fig. 14. Uncontrolled and controlled response using LQG and SMC controllers based on SMO with nominal (a) and 10% mass-density reduction (b) systems subjected to sharp-edged gust (θ = 105◦ , M f = 0.7).

bust to reduce the vertical bending vibration of the wing structure, although transient response and steady-state error are little influenced by parametric uncertainty. Fig. 15 represents the open-loop and closed-loop non-dimensional plunging and pitching responses of the swept wing structure. It is found that for swept-back (positive sweep angle) wing structure, there is the beneficial influence associated with aeroe-

251

Fig. 15. Implication of the sweep angle on the uncontrolled and controlled nondimensional plunging (a) and pitching (b) displacements subjected to sharp-edged gust (θ = 105◦ , M f = 0.6).

Fig. 16. Implication of the ply angle on the uncontrolled and controlled nondimensional plunging (a) and pitching (b) displacements subjected to sharp-edged gust (Λ = 0◦ , M f = 0.6).

lastic response amplitude compared to swept-forward wing structure. In addition, Fig. 16 displays the plunging and pitching response of a wing structure featuring different ply angles. It reveals that the characteristics of fiber-reinforced composite materials used here play an influence on the reduction of response amplitudes related to the vertical bending and twisting displacement. Actually, the case θ = 105◦ generates the lower oscillation amplitude compared with θ = 115◦ .

252

J.-S. Song et al. / Aerospace Science and Technology 25 (2013) 242–252

5. Conclusion In this paper, a comprehensive aeroelastic model of the threedimensional advanced aircraft wing system operating in subsonic compressible flow field and its robust vibration control based on SMO are presented. The excellent performance of sliding mode methodology, which in contrast to linear quadratic methodology, yields robustness to parametric uncertainty, time-dependent external pressure pulses and restricted sensor measurement due to existent sensor failure or observer spill over effect. Furthermore, its application has potential for extension to the compressible transonic and supersonic fields and nonlinear aeroelastic problem. Acknowledgements This work was supported by Basic Science Research Program through the National Research Foundation of Korea (20110001145). References [1] M.J. Balas, Feedback control of flexible systems, IEEE Transactions on Automatic Control AC-23 (1978) 637–679. [2] R.L. Bisplinghoff, H. Ashley, R.L. Halfman, Aeroelasticity, Dover Publications, NY, 1996. [3] T. Degaki, S. Suzuki, Sliding mode control application for two-dimensional active flutter suppression, Transactions of the Japan Society for Aeronautical and Space Sciences 43 (142) (2001) 174–181. [4] P. Dorato, C. Abdallah, V. Cerone, Linear–Quadratic Control: An Introduction, Prentice Hall, NJ, 1995. [5] C.M. Dorling, A.S.I. Zinober, Two approaches to hyperplane design in multivariable variable structure control systems, International Journal of Control 44 (1) (1986) 65–82. [6] C. Edwards, S.K. Spurgeon, Sliding Mode Control: Theory and Applications, Taylor and Francis, PA, 1998. [7] V. Giurgiutiu, Review of smart-materials actuation solutions for aeroelastic and vibration control, Journal of Intelligent Material Systems and Structures 11 (2000) 525–544. [8] J. Heeg, Analytical and experimental investigation of flutter suppression by piezoelectric actuation, NASA Technical Paper 3241, Langley Research Center, Hampton, VA, 1993. [9] D.J. Inman, Active modal control for smart structures, Philosophical Transactions: Mathematical, Physical and Engineering Sciences 359 (1778) (2001) 205– 219. [10] N.J. Krone Jr., Divergence elimination with advanced composites, AIAA Paper 75-1009, 1975. [11] B.H. Lee, J.H. Choo, S.S. Na, P. Marzocca, L. Librescu, Sliding mode robust control of supersonic three degrees-of-freedom airfoils, International Journal of Control, Automation and Systems 8 (2) (2010) 279–288. [12] J. Leishman, Indicial lift approximations for two-dimensional subsonic flow as obtained from oscillatory measurements, Journal of Aircraft 30 (3) (1993) 340– 351. [13] L. Librescu, A.A. Khdeir, Aeroelastic divergence of swept-forward composite wings including warping restraint effect, AIAA Journal 26 (11) (1988) 1373– 1377.

[14] L. Librescu, S.S. Na, Z. Qin, B.H. Lee, Active aeroelastic control of aircraft composite wings impacted by explosive blasts, Journal of Sound and Vibration 318 (1–2) (2008) 74–92. [15] L. Librescu, O. Song, On the static aeroelastic tailoring of composite aircraft swept wings modeled as thin-walled beams structures, in: Use of Composite in Rotor Craft and Smart Structures (Special Issue), Composite Engineering 2 (5–7) (1992) 497–512. [16] L. Librescu, O. Song, Adaptive response control of cantilevered thin-walled beam structure carrying heavy concentrated masses, Journal of Intelligent Material Systems and Structures 5 (1) (1994) 42–48. [17] L. Librescu, O. Song, Thin-Walled Composite Beams – Theory and Application, Springer, NY, October 2005. [18] E. Livne, Future of airplane aeroelasticity, Journal of Aircraft 40 (6) (2003) 1066–1092. [19] P. Marzocca, L. Librescu, G. Chiocchia, Aeroelasticity of two-dimensional lifting surfaces via indicial function approach, Aeronautical Journal 106 (1057) (2002) 147–153. [20] P. Marzocca, L. Librescu, D.H. Kim, I. Lee, S. Schober, Development of an indicial function for the two-dimensional incompressible/compressible aerodynamic load modeling, Proceedings of the Institution of Mechanical Engineering, Part G – Journal of Aerospace Engineering 221 (G3) (2007) 453–463. [21] V. Mukhopadhyay, Historical perspective on analysis and control of aeroelastic responses, Journal of Guidance, Control, and Dynamics 26 (5) (2003) 673–684. [22] S.S. Na, L. Librescu, M.H. Kim, I.J. Jeong, P. Marzocca, Aeroelastic response of flapped-wing systems using robust estimation control methodology, Journal of Guidance, Control, and Dynamics 29 (1) (2006). [23] S.S. Na, L. Librescu, M.H. Kim, I.J. Jeong, P. Marzocca, Robust aeroelastic control of flapped wing system using a sliding mode observer, Aerospace Science and Technology 10 (2) (2006) 120–126. [24] S.S. Na, L. Librescu, P. Marzocca, G.C. Yoon, C. Rubillo, K. Bong, Robust aeroelastic control of two-dimensional supersonic flapped wing systems, Acta Mechanica 192 (1–4) (2007) 37–47. [25] S.S. Na, G.C. Yoon, B.H. Lee, P. Marzocca, L. Librescu, Dynamic response control of rotating composite booms under solar radiation, Journal of Thermal Stresses 32 (1) (2009) 21–40. [26] J. Njuguna, Flutter prediction, suppression and control in aircraft composite wings as a design prerequisite: A survey, Structural Control and Health Monitoring 14 (2007) 715–758. [27] A. Preumont, Vibration Control of Active Structures, Kluwer Academic Publishers, 2002. [28] Z. Qin, Vibration and aeroelasticity of advanced aircraft wings modeled as thin-walled beams-dynamics, stability and control, Ph.D. thesis, VPI&SU, USA, 2001. [29] Z. Qin, L. Librescu, Aeroelastic instability of aircraft wings modeled as anisotropic composite thin-walled beams in incompressible flow, Journal of Fluids and Structures 18 (1) (2003) 43–61. [30] Z. Qin, P. Marzocca, L. Librescu, Aeroelastic instability and response of advanced aircraft wings at subsonic flight speeds, Aerospace Science and Technology 6 (3) (2002) 195–208. [31] L.W. Rehfield, A.R. Atilgan, Analysis, design and elastic tailoring of composite rotor blades, NASA Technical Report, NASA-CR-181234, 1987. [32] M.A. Shubov, Flutter phenomenon in aeroelasticity and its mathematical analysis, Journal of Aerospace Engineering 19 (1) (2006) 1–12. [33] A. Sinha, Linear Systems: Optimal and Robust Control, CRC Press, 2007. [34] A. Sinha, D.W. Miller, Optimal sliding-mode control of a flexible spacecraft under stochastic disturbances, Journal of Guidance, Control, and Dynamics 18 (3) (1995) 486–492. [35] O. Song, L. Librescu, Bending vibration of cantilevered thin-walled beams subjected to time-dependent external excitation, The Journal of the Acoustical Society of America 98 (1) (1995) 313–319. [36] A.A. Umansky, Torsion and Bending of Thin-Walled Aircraft Structures, Oboronghiz, Moscow, 1939.