# Control of Limit Cycle Oscillation in a Three Degrees of Freedom Airfoil Section Using Fuzzy Takagi-Sugeno Modeling

Authors: Douglas Domingues Bueno; Luiz Carlos Sandoval Góes; Paulo José Paupitz Gonçalves
Journal: Shock and Vibration
 (2014)

---

## Abstract

This work presents a strategy to control nonlinear responses of aeroelastic systems with control
surface freeplay. The proposed methodology is developed for the three degrees of freedom typical
section airfoil considering aerodynamic forces from Theodorsen’s theory. The mathematical model
is written in the state space representation using rational function approximation to write the
aerodynamic forces in time domain. The control system is designed using the fuzzy Takagi-Sugeno modeling to compute a feedback control gain. It useds Lyapunov’s stability function and linear matrix inequalities (LMIs) to solve a convex optimization problem. Time simulations
with different initial conditions are performed using a modified Runge-Kutta algorithm to compare
the system with and without control forces. It is shown that this approach can compute linear
control gain able to stabilize aeroelastic systems with discontinuous nonlinearities.


---

## Body

## 1. Introduction

The requirement for more accurate tools for predictions of nonlinear effects has motivated many research groups to investigate aeroelastic systems considering nonlinearities. In particular, the problem involving freeplay in control surfaces has called attention of various researchers because it can be a cause of limit cycle oscillation (LCO) leading to serious consequences such as fatigue, pilot handling/ride quality, confined manoeuvrings envelope, weapon aiming of military aircraft, and induced flutter.Another motivation to consider control surface freeplay is that the requirements for aircraft design according to military specification can be quite difficult to achieve in practice, increasing the manufacturing and maintenance costs. Considerable experimental and analytical efforts have been devoted to obtain representative aeroelastic models and develop methodologies to study the freeplay problem.In this context, many works in literature have presented studies to understand and characterize nonlinear aeroelastic behaviour. Conner et al. [1, 2] presented results for a typical airfoil section based on time domain simulations. The authors showed accuracy between numerical and experimental data. Tang and colleagues published theoretical and experimental results considering an aeroelastic apparatus and the high order Harmonic Balance methods [3]. This method was introduced by Kryloff and Bogoliuboff in 1947 and it has been studied by different researchers as shown in [4–8].Kholodar and Dickinson studied the effects of aileron freeplay in different configurations of a real aircraft [9]. Time domain simulations were used to confirm the limit cycles previously predicted using the Harmonic Balance method. Also considering similar approaches, Anderson and Mortara presented results for the F-22 aircraft including control surface freeplay [10]. The authors discussed the limits of freeplay to keep the system stable. Recently, Abdelkefi and colleagues performed numerical and experimental investigations using a pitch and plunge rigid airfoil supported by a torsional spring. The authors studied different mathematical representations of freeplay nonlinearity, such as polynomial expansion and hyperbolic tangent [11, 12]. Also considering a two degrees of freedom airfoil, Guo and Chen employed multivariable and Floquet theories to detect the fold bifurcation and amplitude jump phenomenon in supersonic flow [13].Several control strategies with focus on LCO suppression on aeroelastic systems have been developed in these last years. Kurdila et al. [14] presented an extensive review of nonlinear control methods for high energy LCO. Experimental works on controlling aeroelastic apparatus are presented in [15, 16]. In [17] Li et al. designed a suboptimal controller using the state-dependent Riccati equation considering cubic nonlinearity. Adaptive filters cut-off frequency with feedback gain has also been used to suppress limit cycles and chaotic motions [18]. The authors investigated an augmented controller with time delay parameter to determine regions of instabilities in closed-loop configurations.This paper proposes a control strategy based on fuzzy Takagi-Sugeno (FTS) solved using linear matrix inequalities (LMIs) to control the LCO of a nonlinear aeroelastic system. Techniques based on LMIs have been used to solve linear aeroelastic problems, mainly considering structural uncertainties [19, 20]. However, their application for solving nonlinear aeroelastic problems is rarely found in the literature. The advantage of using LMIs to design a control strategy is based on the robust interior point algorithms that provides a guarantee for finding optimal solutions, if that exists.The aeroelastic system is formulated in state space form and is integrated in time domain using the forth order Runge-Kutta method and Hénon’s technique [21]. Hénon’s technique is used to locate the switching points in the procedure of numerical integration, as discussed herein. Theodorsen aerodynamic forces are transformed to time domain using rational function approximation. However, the proposed approach is valid for other aerodynamic theories. Finally, numerical simulations are performed on the benchmark airfoil problem to demonstrate that LMIs combined with FTS modeling can be used to design controllers for nonlinear aeroelastic problems.


## 2. Aeroelastic System with Freeplay

The typical section airfoil shown in Figure1 with a trailing edge control surface is normally used to represent an aeroelastic system with three degrees of freedom that includes pitch α
(
t
), plunge h
(
t
), and control surface rotation β
(
t
). This modeling was previously proposed in [22]. The structural properties of this system are represented by the springs k
α, k
h, and k
β, structural damping, and inertial properties. Theodorsen theory is used to compute the aerodynamic forces; however, the same control strategy could be adopted for other aerodynamic theories [22, 23].Figure 1
Typical section airfoil.The nonlinear discontinuity (freeplay) is considered to occur in the control surface springk
β. So that the equation of motion for this system can be represented by
(1)
M
u
¨
(
t
)
+
D
u
˙
(
t
)
+
F
(
K
,
u
,
t
)
=
q
Q
u
(
t
)
+
B
o
c
u
c
(
t
)
,
where M and D are, respectively, the structural mass and damping matrices, Q is a matrix of aerodynamic coefficients that depend on the airfoil geometry, q
=
(
1
/
2
)
ρ
V
2 is the dynamic pressure, and ρ is the air density and V is the airspeed. B
o
c is a matrix of input (forces or moments). The vector u
(
t
)
=
{
h
(
t
)
α
(
t
)
β
(
t
)
}
T represents the physical displacements and u
c
(
t
) is the control force. The vector F represents the elastic restoring moment which depends on the control surface restoring moment T
β
(
t
) shown in Figure 2.Figure 2
Freeplay nonlinearity.By considering a freeplay amplitude of2
δ, the elastic restoring moment can be written as
(2)
F
(
K
,
u
,
t
)
=
{
F
h
(
t
)
T
α
(
t
)
T
β
(
t
)
}
T
such that
F
h
(
t
)
=
k
h
h
(
t
)
,
T
α
(
t
)
=
k
α
α
(
t
)
,
T
β
(
t
)
=
0
if
|
β
(
t
)
|
≤
δ
,
T
β
(
t
)
=
k
β
[
β
(
t
)
-
δ
]
if
|
β
(
t
)
|
>
δ
.The expressions in (2) are rearranged in a matrix form and the nonlinear function f
nl is introduced such that the stiffness is defined as a nonlinear structural stiffness matrix K
nl given by
(3)
K
nl
=
[
k
h
0
0
0
k
α
0
0
0
k
β
f
n
l
]
⟹
F
(
K
,
u
,
t
)
=
K
nl
u
(
t
)
,
where
(4)
f
nl
=
0
if
|
β
(
t
)
|
≤
δ
,
f
nl
=
1
-
δ
β
(
t
)
if
|
β
(
t
)
|
>
δ
.


### 2.1. Representation in State Space Form

In order to obtain a nonlinear state space model, a rational function is used to write the aerodynamic forces in time domain. Jones in 1940 solved this problem using a rational function approximation to approximate unsteady aerodynamic loads for the typical section [24].After Jones’ work, different formulations have been used to approximate generalized aerodynamic forces for arbitrary motion, as, for instance, the approach based on Chebyshev’s polynomials introduced by Dinu and colleagues with focus on aeroservoelasticity [25]. Karpel proposes the Minimum-State method with accuracy per model order superior to previous works, but it is more complicated and computing-time consuming because it involves nonlinear problems [26]. Vepa proposes a numerical technique based on Padè approximation and according to the author the main advantage of this method is that it can be generalized to three-dimensional lifting surfaces [27]. Recently, Biskri and colleagues present a very interesting method based on a combination of the Least Squares (or Roger’s approach) and Minimum-State methods, of which the main idea is to minimize the number of lag terms without passing through a long iterative algorithm [28].In this paper Roger’s approach is used. This approximation involves identifying every matrixQ
j shown in (5) using a Least Square algorithm as proposed by Roger [29] and summarized in [30]. The approximation contains a polynomial part representing the forces on the typical section acting directly connected to the displacements u
(
t
) and their first and second derivatives. Also, this equation has a rational part representing the influence of the wake acting on the section with a time delay. Consider
(5)
Q
(
s
)
≈
[
∑
j
=
0
2
‍
Q
j
s
j
(
b
V
)
j
+
∑
j
=
1
n
lag
‍
Q
(
j
+
2
)
(
s
s
+
(
b
/
V
)
β
j
)
]
u
(
s
)
,
where s is the Laplace variable, n
lag is the number of lag terms, β
j is the jth lag parameter (
j
=
1
,
…
,
n
lag
), and b is the aerodynamic semichord.Substituting (5) into (1) and considering (3) it is possible to write the equation of motion for the aeroelastic system in state space format:
(6)
x
˙
(
t
)
=
A
nl
x
(
t
)
+
B
u
c
(
t
)
,
y
(
t
)
=
C
x
(
t
)
,
where x
(
t
)
=
{
u
˙
(
t
)
u
(
t
)
u
a
(
t
)
}
T is the state vector and u
a
(
t
) are states of lags required for the approximation of matrix Q. The matrix of outputs C
=
[
C
v
C
d
C
z
] has dimension 2
m
×
m
(
2
+
n
lag
), where C
v and C
d are, respectively, the velocity and displacement output matrices and the submatrix C
z has only zeros to complete the matrix dimension, and m is the number of degrees of freedom. Matrix A
nl is presented in the following form:
(7)
A
nl
=
[
-
M
a
-
1
D
a
-
M
a
-
1
K
a
(
nl
)
q
M
a
-
1
Q
3
⋯
q
M
a
-
1
Q
(
2
+
n
lag
)
I
0
0
⋯
0
I
0
(
-
V
b
)
β
1
I
0
⋯
⋮
⋮
0
⋱
⋯
I
0
⋮
⋯
(
-
V
b
)
β
n
lag
I
]
,
where
(8)
M
a
=
M
-
q
(
b
V
)
2
Q
2
,
D
a
=
D
-
q
(
b
V
)
Q
1
,
K
a
(
nl
)
=
K
nl
-
q
Q
0
,
B
=
[
M
a
-
1
B
o
c
0
]
T
,
where M
a and D
a are, respectively, the aeroelastic mass and damping matrices, K
a
(
nl
) is the nonlinear aeroelastic stiffness matrix, B is the input matrix, I is an identity matrix, and 0 is a matrix of zeros with appropriated dimension.


### 2.2. Time Integration

To calculate the time response of the system with freeplay a modified 4th order Runge-Kutta algorithm using Hénon method is used. Hénon method is necessary in order to identify changes in the stiffness due to the freeplay. This approach is used to minimize integration errors mainly with respect to the phase shifts.The main idea considered in Hénon’s method is to change the independent variable (time) to the degree of freedom with freeplay (a spacial variable) always that the stiffness changes with the freeplay region. In these cases, the time becomes the dependent variable and the integration step is done in the degree of freedom related to the freeplay. The size of the step is the amount necessary for this degree of freedom to coincide with the transition points shown in Figure2. Once the control surface position is the transition point, the system of equations is rewritten considering the time as the independent variable. See complementary information in [2, 21, 31].


## 2.1. Representation in State Space Form

In order to obtain a nonlinear state space model, a rational function is used to write the aerodynamic forces in time domain. Jones in 1940 solved this problem using a rational function approximation to approximate unsteady aerodynamic loads for the typical section [24].After Jones’ work, different formulations have been used to approximate generalized aerodynamic forces for arbitrary motion, as, for instance, the approach based on Chebyshev’s polynomials introduced by Dinu and colleagues with focus on aeroservoelasticity [25]. Karpel proposes the Minimum-State method with accuracy per model order superior to previous works, but it is more complicated and computing-time consuming because it involves nonlinear problems [26]. Vepa proposes a numerical technique based on Padè approximation and according to the author the main advantage of this method is that it can be generalized to three-dimensional lifting surfaces [27]. Recently, Biskri and colleagues present a very interesting method based on a combination of the Least Squares (or Roger’s approach) and Minimum-State methods, of which the main idea is to minimize the number of lag terms without passing through a long iterative algorithm [28].In this paper Roger’s approach is used. This approximation involves identifying every matrixQ
j shown in (5) using a Least Square algorithm as proposed by Roger [29] and summarized in [30]. The approximation contains a polynomial part representing the forces on the typical section acting directly connected to the displacements u
(
t
) and their first and second derivatives. Also, this equation has a rational part representing the influence of the wake acting on the section with a time delay. Consider
(5)
Q
(
s
)
≈
[
∑
j
=
0
2
‍
Q
j
s
j
(
b
V
)
j
+
∑
j
=
1
n
lag
‍
Q
(
j
+
2
)
(
s
s
+
(
b
/
V
)
β
j
)
]
u
(
s
)
,
where s is the Laplace variable, n
lag is the number of lag terms, β
j is the jth lag parameter (
j
=
1
,
…
,
n
lag
), and b is the aerodynamic semichord.Substituting (5) into (1) and considering (3) it is possible to write the equation of motion for the aeroelastic system in state space format:
(6)
x
˙
(
t
)
=
A
nl
x
(
t
)
+
B
u
c
(
t
)
,
y
(
t
)
=
C
x
(
t
)
,
where x
(
t
)
=
{
u
˙
(
t
)
u
(
t
)
u
a
(
t
)
}
T is the state vector and u
a
(
t
) are states of lags required for the approximation of matrix Q. The matrix of outputs C
=
[
C
v
C
d
C
z
] has dimension 2
m
×
m
(
2
+
n
lag
), where C
v and C
d are, respectively, the velocity and displacement output matrices and the submatrix C
z has only zeros to complete the matrix dimension, and m is the number of degrees of freedom. Matrix A
nl is presented in the following form:
(7)
A
nl
=
[
-
M
a
-
1
D
a
-
M
a
-
1
K
a
(
nl
)
q
M
a
-
1
Q
3
⋯
q
M
a
-
1
Q
(
2
+
n
lag
)
I
0
0
⋯
0
I
0
(
-
V
b
)
β
1
I
0
⋯
⋮
⋮
0
⋱
⋯
I
0
⋮
⋯
(
-
V
b
)
β
n
lag
I
]
,
where
(8)
M
a
=
M
-
q
(
b
V
)
2
Q
2
,
D
a
=
D
-
q
(
b
V
)
Q
1
,
K
a
(
nl
)
=
K
nl
-
q
Q
0
,
B
=
[
M
a
-
1
B
o
c
0
]
T
,
where M
a and D
a are, respectively, the aeroelastic mass and damping matrices, K
a
(
nl
) is the nonlinear aeroelastic stiffness matrix, B is the input matrix, I is an identity matrix, and 0 is a matrix of zeros with appropriated dimension.


## 2.2. Time Integration

To calculate the time response of the system with freeplay a modified 4th order Runge-Kutta algorithm using Hénon method is used. Hénon method is necessary in order to identify changes in the stiffness due to the freeplay. This approach is used to minimize integration errors mainly with respect to the phase shifts.The main idea considered in Hénon’s method is to change the independent variable (time) to the degree of freedom with freeplay (a spacial variable) always that the stiffness changes with the freeplay region. In these cases, the time becomes the dependent variable and the integration step is done in the degree of freedom related to the freeplay. The size of the step is the amount necessary for this degree of freedom to coincide with the transition points shown in Figure2. Once the control surface position is the transition point, the system of equations is rewritten considering the time as the independent variable. See complementary information in [2, 21, 31].


## 3. Development of Fuzzy Takagi-Sugeno Controller

A fuzzy model uses rules, which are linguistic IF-THEN statements involving fuzzy sets, fuzzy logic, and fuzzy inference. These rules play a key role in representing expert control/modeling knowledge and experience and in linking the input variables of fuzzy controllers/models to output variable (or variables). To explain the procedure, consider the open loop nonlinear aeroelastic system described in following form:(9)
x
˙
i
(
t
)
=
∑
j
=
1
N
‍
f
i
j
(
x
(
t
)
)
x
j
(
t
)
,
where u
c
=
0, N
=
n
(
2
+
n
lag
), i
=
1
,
…
,
N, and f
i
j
(
x
(
t
)
)
=
f
k represents the kth nonlinear function, where k
=
1
,
…
,
N
nl
≤
N
2.The nonlinear system described by (9) can be represented by the Takagi-Sugeno model using the following rule:
(10)
Rule
i
:
if
z
1
(
t
)
is
M
i
1
and
…
and
z
p
(
t
)
is
M
i
p
Then
x
˙
(
t
)
=
A
i
x
(
t
)
,
y
(
t
)
=
C
x
(
t
)
,
i
=
1
,
…
,
r
,
where M
i
j is the fuzzy set and r is the number of model rules; z
1
(
t
)
,
…
,
z
p
(
t
) are known premise variables that in general may be functions of the state variables, external disturbances, and/or times. Each linear model represented by A
i is called a subsystem [32].Taniguchi et al. [33] present a simple method to identify the subsystems. The basic idea is to write each nonlinear function f
k as a linear combination of its maximum f
k
max
⁡ and minimum f
k
min
⁡ values both given, respectively, by
(11)
f
k
max
⁡
=
max
⁡
[
f
k
(
x
(
t
)
)
]
,
f
k
min
⁡
=
min
⁡
[
f
k
(
x
(
t
)
)
]
.From these maximum and minimum values, each nonlinear function can be represented by(12)
f
k
(
x
(
t
)
)
=
[
g
k
min
⁡
(
x
(
t
)
)
]
f
k
min
⁡
+
[
g
k
max
⁡
(
x
(
t
)
)
]
f
k
max
⁡
,
where
(13)
0
≤
g
k
min
⁡
(
x
(
t
)
)
,
g
k
max
⁡
(
x
(
t
)
)
≤
1
,
g
k
min
⁡
(
x
(
t
)
)
+
g
k
max
⁡
(
x
(
t
)
)
=
1
,
g
k
min
⁡
(
x
(
t
)
)
=
[
f
k
(
x
(
t
)
)
-
f
k
max
⁡
]
[
f
k
min
⁡
-
f
k
max
⁡
]
then
,
g
k
max
⁡
(
x
(
t
)
)
=
1
-
g
k
min
⁡
(
x
(
t
)
)
.If[
g
k
min
⁡
(
x
(
t
)
)
+
g
k
max
⁡
(
x
(
t
)
)
=
1
], ∀
k
=
1
,
…
,
N
nl, then each nonlinear function can be conveniently written as
(14)
f
k
(
x
(
t
)
)
=
∏
j
=
1
,
j
≠
k
N
nl
‍
{
g
j
min
⁡
(
x
(
t
)
)
+
g
j
max
⁡
(
x
(
t
)
)
}
d
d
d
d
d
d
d
d
×
{
[
g
k
min
⁡
(
x
(
t
)
)
]
f
k
min
⁡
+
[
g
k
max
⁡
(
x
(
t
)
)
]
f
k
max
⁡
}
,
or
(15)
f
k
(
x
(
t
)
)
=
∑
j
=
1
r
1
‍
N
j
(
x
(
t
)
)
f
k
min
⁡
+
∑
j
=
r
1
+
1
(
2
N
nl
)
‍
N
j
(
x
(
t
)
)
f
k
max
⁡
,
where r
1
=
2
(
N
nl
-
1
) and N
j
(
x
(
t
)
) is given by
(16)
N
j
(
x
(
t
)
)
=
∏
k
=
1
N
nl
‍
g
k
(
·
)
(
x
(
t
)
)
,
∑
j
=
1
2
N
nl
‍
N
j
(
x
(
t
)
)
=
1
,
and the superscript (
·
) indicates the combination between maximum and minimum values for each kth function g
k
(
·
)
(
x
(
t
)
).Considering thatA
nl
=
A
nl
(
f
1
,
f
2
,
…
,
f
N
nl
) and substituting each f
k according to (15), the nonlinear equation of motion (17) is rewritten after some rearrangement as
(17)
x
˙
(
t
)
=
(
∑
j
=
1
2
N
nl
‍
N
j
(
x
(
t
)
)
A
j
)
x
(
t
)
+
(
∑
j
=
1
2
N
nl
‍
N
j
(
x
(
t
)
)
B
j
)
u
c
(
t
)
,
where each linear matrix A
j is the original nonlinear matrix at which the N
nl functions f
k are substituted by the combination of maximum and minimum values f
k
max
⁡ and f
k
min
⁡. In this case B
j
=
B
∀
j.Particularly for the typical section including freeplay nonlinearity, there are three nonlinear functions(
N
nl
=
3
) into the matrix A
nl and they are given by
(18)
f
1
=
a
nl
(
1,6
)
=
-
M
a
(
1,3
)
-
1
[
k
β
f
nl
-
ρ
0.5
V
2
Q
0
(
1,3
)
]
,
f
2
=
a
nl
(
2,6
)
=
-
M
a
(
2,3
)
-
1
[
k
β
f
nl
-
ρ
0.5
V
2
Q
0
(
2,3
)
]
,
f
3
=
a
nl
(
3,6
)
=
-
M
a
(
3,3
)
-
1
[
k
β
f
nl
-
ρ
0.5
V
2
Q
0
(
3,3
)
]
,
where a
nl
(
i
,
j
) indicates an element into matrix A
nl and the subscripts (
i
,
j
) indicate the ith row and jth column in each respective matrix. These functions f
k depend on f
nl which is illustrated in Figure 3.Figure 3
The nonlinear functionf
nl represented in a particular range of β
(
t
).


### 3.1. Closed-Loop Takagi-Sugeno Controller

In this section a feedback control strategy is used to control the amplitude of oscillation in the control surface. The feedback force is obtained by applying a feedback gain to the state systems states, such that the control inputs are given byu
c
(
t
)
=
-
G
x
(
t
), where G is the feedback gain matrix. The feedback gain is calculated using linear matrix inequality to solve Lyapunov’s function for stability, in this case
(19)
x
˙
(
t
)
=
[
∑
j
=
1
2
N
nl
‍
N
j
(
x
(
t
)
)
(
A
j
-
B
j
G
)
]
x
(
t
)
,
(20)
V
˙
L
=
(
∑
j
=
1
2
N
nl
‍
N
j
(
x
(
t
)
)
[
A
j
-
B
j
G
]
x
(
t
)
)
T
P
x
(
t
)
+
x
(
t
)
T
P
(
∑
j
=
1
2
N
nl
‍
N
j
(
x
(
t
)
)
[
A
j
-
B
j
G
]
x
(
t
)
)
<
0
,
where V
L
=
x
T
P
x is the Lyapunov function. The stability of this system is assured if there is a positive-definite matrix P such that the inequality of (20) is true. After some rearrangement this equation can be rewritten such as [32]
(21)
x
T
(
t
)
∑
j
=
1
2
N
nl
‍
N
j
(
x
(
t
)
)
s
s
s
s
s
s
s
s
s
×
[
(
A
j
-
B
j
G
)
T
P
+
P
(
A
j
-
B
j
G
)
]
x
(
t
)
<
0
,
and finally, the solution of this inequality is equivalent to the solution of all the following inequalities satisfied simultaneously [34]:
(22)
(
A
j
-
B
j
G
)
T
P
+
P
(
A
j
-
B
j
G
)
<
0
j
=
1
,
…
,
2
N
nl
.
Inequality (22) is not linear: then considering X
=
P
-
1, G
x
=
G
X, it is possible to write the following linear matrix inequality:
(23)
X
A
j
T
-
G
x
T
B
j
T
+
A
j
X
-
B
j
G
x
<
0
j
=
1
,
…
,
2
N
nl
such that
G
=
G
x
X
-
1
.Although the control gainG has been computed using inequality (23), the required control force u
c
(
t
) can exceed desired limits. In this case, the control gain is multiplied by a constant g between the interval [
0
1
], conveniently adjusted during the control design. The control force is rewritten as
(24)
u
¯
c
(
t
)
=
g
u
c
(
t
)
or
u
¯
c
(
t
)
=
-
g
G
x
(
t
)
.


## 3.1. Closed-Loop Takagi-Sugeno Controller

In this section a feedback control strategy is used to control the amplitude of oscillation in the control surface. The feedback force is obtained by applying a feedback gain to the state systems states, such that the control inputs are given byu
c
(
t
)
=
-
G
x
(
t
), where G is the feedback gain matrix. The feedback gain is calculated using linear matrix inequality to solve Lyapunov’s function for stability, in this case
(19)
x
˙
(
t
)
=
[
∑
j
=
1
2
N
nl
‍
N
j
(
x
(
t
)
)
(
A
j
-
B
j
G
)
]
x
(
t
)
,
(20)
V
˙
L
=
(
∑
j
=
1
2
N
nl
‍
N
j
(
x
(
t
)
)
[
A
j
-
B
j
G
]
x
(
t
)
)
T
P
x
(
t
)
+
x
(
t
)
T
P
(
∑
j
=
1
2
N
nl
‍
N
j
(
x
(
t
)
)
[
A
j
-
B
j
G
]
x
(
t
)
)
<
0
,
where V
L
=
x
T
P
x is the Lyapunov function. The stability of this system is assured if there is a positive-definite matrix P such that the inequality of (20) is true. After some rearrangement this equation can be rewritten such as [32]
(21)
x
T
(
t
)
∑
j
=
1
2
N
nl
‍
N
j
(
x
(
t
)
)
s
s
s
s
s
s
s
s
s
×
[
(
A
j
-
B
j
G
)
T
P
+
P
(
A
j
-
B
j
G
)
]
x
(
t
)
<
0
,
and finally, the solution of this inequality is equivalent to the solution of all the following inequalities satisfied simultaneously [34]:
(22)
(
A
j
-
B
j
G
)
T
P
+
P
(
A
j
-
B
j
G
)
<
0
j
=
1
,
…
,
2
N
nl
.
Inequality (22) is not linear: then considering X
=
P
-
1, G
x
=
G
X, it is possible to write the following linear matrix inequality:
(23)
X
A
j
T
-
G
x
T
B
j
T
+
A
j
X
-
B
j
G
x
<
0
j
=
1
,
…
,
2
N
nl
such that
G
=
G
x
X
-
1
.Although the control gainG has been computed using inequality (23), the required control force u
c
(
t
) can exceed desired limits. In this case, the control gain is multiplied by a constant g between the interval [
0
1
], conveniently adjusted during the control design. The control force is rewritten as
(24)
u
¯
c
(
t
)
=
g
u
c
(
t
)
or
u
¯
c
(
t
)
=
-
g
G
x
(
t
)
.


## 4. Numerical Application

To illustrate the method, numerical simulations were performed using the three degrees of freedom airfoil section for which the equations of motion are presented by Theodorsen et al. in [23]. The structural mass, stiffness, and aerodynamic forces matrices can be found in [17, 22, 35].Initially, the linear flutter boundary was found extracting the eigenvalues from the state space dynamic matrix without freeplay. After that preliminary verification, the first order harmonic balance method (HBM) was used to predict the LCO amplitudes. Different researchers have used HBM methods to study limit cycles oscillations in aeroelastic systems, as shown in [3, 6, 36, 37].


### 4.1. Linear Flutter Boundary

Figure1 illustrates the model and its physical and geometric properties are presented in Table 1. Figures 4 and 5 show, respectively, the classical V-f and V-g diagrams for the linear flutter solution. According to these results flutter speed is equal to 12.7 m/s.Table 1
Physical and geometric properties of the 2D airfoil.
Parameter
Value
Semichord—b
0.15 m
Airfoil mass—M
5.0 kg
Air density—ρ
1.225 kg/m3
c
0.6%
a
-
0.4%
x
-
θ
0.2 m/m
r
-
β
(
6.25
×
1
0
-
3
)
-
1
/
2 m/m
x
-
β
0.0125 m/m
r
-
θ
(
0.25
)
-
1
/
2 m/m
Plunge frequency—f
h
3.0 Hz
Pitch frequency—f
θ
4.5 Hz
Control surface deflection—f
β
12.0 Hz
Parameters of lag—β
j
(
j
=
1
,
…
,
4
)
0.2, 1.2, 1.6, 1.8
Reduced frequency—k
[0.1, 2.0],Δ
k
=
0.1Figure 4
Linear flutter analysis:V-f diagram.Figure 5
Linear flutter analysis:V-g diagram.


### 4.2. LCO Preliminary Predictions

After the preliminary linear analysis to identify the flutter boundary, the first order HBM was used to predict the LCO amplitudes. The results shown in Figures6 and 7 were obtained by extracting the eigenvalues from the matrix A
nl defined for different values of equivalent stiffness k
eq, that is, the control surface stiffness assuming values from zero to k
β. The values of k
eq
/
k
β are shown in Figure 8, where the nonlinear function f
nl assumed a unitary value and the control force u
c
(
t
)
=
0. The HBM also provided an estimate for the first harmonic of the system response according to Figure 9.Figure 6
LCO prediction for|
δ
|
=
0
.
4
∘.Figure 7
LCO prediction for|
δ
|
=
0
.
6
∘.Figure 8
Linear equivalent stiffness from the HBM.Figure 9
LCO frequencies obtained from the HBM.


### 4.3. Controller Design

Using the methodology described in Section4, three nonlinear functions were used to describe the aeroelastic system with freeplay (f
k,  k
=
1,2
,
3). These functions were computed assuming the maximum oscillation amplitude equal to five degrees; that is, -
5
°
≤
β
(
t
)
≤
5
°. Figure 10 shows a comparison between their actual values and computed values by (12).Figure 10
Values off
k: actual versus computed using TS approach.This leads to eight(
2
3
) dynamic matrices that are written combining the maximum and minimum values of these functions, such that
(25)
A
1
=
A
(
f
1
max
⁡
,
f
2
min
⁡
,
f
3
min
⁡
)
,
A
2
=
A
(
f
1
max
⁡
,
f
2
max
⁡
,
f
3
min
⁡
)
,
A
3
=
A
(
f
1
min
⁡
,
f
2
max
⁡
,
f
3
min
⁡
)
,
A
4
=
A
(
f
1
min
⁡
,
f
2
min
⁡
,
f
3
min
⁡
)
,
A
5
=
A
(
f
1
max
⁡
,
f
2
min
⁡
,
f
3
max
⁡
)
,
A
6
=
A
(
f
1
max
⁡
,
f
2
max
⁡
,
f
3
max
⁡
)
,
A
7
=
A
(
f
1
min
⁡
,
f
2
max
⁡
,
f
3
max
⁡
)
,
A
8
=
A
(
f
1
min
⁡
,
f
2
min
⁡
,
f
3
max
⁡
)
.The results shown in this section were obtained considering two different freeplay amplitudes and a time step equal to 0.001 seconds. To compute the control gain, a column input matrix was used to represent an actuator applying force on the control surface (B
o
c
=
{
0
0
1
}
T). Also, it was considered that the parameter g is equal to 1.0 for all cases. However, in particular for experimental applications where limitations exist on the actuator force, the parameter g can be conveniently chosen between ]0,1
[.For the cases in which the freeplay amplitude was|
δ
|
=
0
.
4
°, the initial conditions were β
(
0
)
=
0
.
9
° and 2
.
5
°. Similarly, for the cases in which |
δ
|
=
0
.
6
°, the initial conditions were β
(
0
)
=
1
.
3
° and 3
.
5
°. These values were defined based on the results obtained from the HBM, as discussed in Section 4.2.Two different conditions of airspeed were considered to demonstrate the method. Different controllers were designed for each case. In first case was chosen a condition of stable limit cycle for each freeplay amplitude. In the second one, unstable limit cycles were chosen in order to evaluate the effectiveness of the approach. Consider if all neighboring trajectories approach the LCO, it is stable; otherwise, the LCO is unstable. Definitions and a detailed discussion about stable and unstable limit cycles can be found in [38].Case 1 (stable LCO).
In the first caseV
=
12.5 m/s was considered for both freeplay amplitudes. Figures 11 and 12 show that the uncontrolled system response exhibits a limit cycle oscillation with first harmonic around 3.7 Hertz and amplitude approximately equal to 1 degree. It is possible to note in those figures that although the first order HBM cannot predict, the system response exhibits a dominant component around 15 Hertz. In addition, the plunge and pitch degrees of freedom behavior can be observed in the phase plan shown in Figures 13 and 14. In these cases of stable LCOs the designed controllers were able to suppress the limit cycles. The appendix presents details involved to compute the Fourier transform when Hénon’s technique is used for the time integration process.Figure 11
Case1: control surface rotation (|
δ
|
=
0
.
4
° and V
=
12.5 m/s).Figure 12
Case1: control surface rotation (|
δ
|
=
0
.
6
° and V
=
12.5 m/s).Figure 13
Case1: phase plan (|
δ
|
=
0
.
4
° and V
=
12.5 m/s).Figure 14
Case: phase plan (|
δ
|
=
0
.
6
° and V
=
12.5 m/s).Case 2 (unstable LCO).
In this case was consideredV
=
13.1 m/s for both freeplay amplitudes. Figures 15 and 16 show that the uncontrolled system response exhibits an unstable limit cycle. However, the controller gains computed using the proposed approach are able to suppress those nonlinear responses. To improve the comprehension, Figure 17 shows a zoom selection area for the control surface rotation considering the second freeplay amplitude.Figure 15
Case2: control surface rotation (|
δ
|
=
0
.
4
° and V
=
13.1 m/s).Figure 16
Case2: control surface rotation (|
δ
|
=
0
.
6
° and V
=
13.1 m/s).Figure 17
Case2: zoom selection area for the control surface rotation (|
δ
|
=
0
.
6
° and V
=
13.1 m/s).For this second case the plunge and pitch degrees of freedom behavior can be observed in the phase plan shown in Figures18 and 19. Finally, Figure 20 shows the control force for δ
F
=
0
.
6
° and similar results were obtained for the other cases. Note that the parameter g
<
1.0 can be chosen in order to decrease the required control force mainly for the first time steps.Figure 18
Case2: phase plan (|
δ
|
=
0
.
4
° and V
=
13.1 m/s).Figure 19
Case2: phase plan (|
δ
|
=
0
.
6
° and V
=
13.1 m/s).Figure 20
Case2: control force (|
δ
|
=
0
.
6
° and V
=
13.1 m/s).


## 4.1. Linear Flutter Boundary

Figure1 illustrates the model and its physical and geometric properties are presented in Table 1. Figures 4 and 5 show, respectively, the classical V-f and V-g diagrams for the linear flutter solution. According to these results flutter speed is equal to 12.7 m/s.Table 1
Physical and geometric properties of the 2D airfoil.
Parameter
Value
Semichord—b
0.15 m
Airfoil mass—M
5.0 kg
Air density—ρ
1.225 kg/m3
c
0.6%
a
-
0.4%
x
-
θ
0.2 m/m
r
-
β
(
6.25
×
1
0
-
3
)
-
1
/
2 m/m
x
-
β
0.0125 m/m
r
-
θ
(
0.25
)
-
1
/
2 m/m
Plunge frequency—f
h
3.0 Hz
Pitch frequency—f
θ
4.5 Hz
Control surface deflection—f
β
12.0 Hz
Parameters of lag—β
j
(
j
=
1
,
…
,
4
)
0.2, 1.2, 1.6, 1.8
Reduced frequency—k
[0.1, 2.0],Δ
k
=
0.1Figure 4
Linear flutter analysis:V-f diagram.Figure 5
Linear flutter analysis:V-g diagram.


## 4.2. LCO Preliminary Predictions

After the preliminary linear analysis to identify the flutter boundary, the first order HBM was used to predict the LCO amplitudes. The results shown in Figures6 and 7 were obtained by extracting the eigenvalues from the matrix A
nl defined for different values of equivalent stiffness k
eq, that is, the control surface stiffness assuming values from zero to k
β. The values of k
eq
/
k
β are shown in Figure 8, where the nonlinear function f
nl assumed a unitary value and the control force u
c
(
t
)
=
0. The HBM also provided an estimate for the first harmonic of the system response according to Figure 9.Figure 6
LCO prediction for|
δ
|
=
0
.
4
∘.Figure 7
LCO prediction for|
δ
|
=
0
.
6
∘.Figure 8
Linear equivalent stiffness from the HBM.Figure 9
LCO frequencies obtained from the HBM.


## 4.3. Controller Design

Using the methodology described in Section4, three nonlinear functions were used to describe the aeroelastic system with freeplay (f
k,  k
=
1,2
,
3). These functions were computed assuming the maximum oscillation amplitude equal to five degrees; that is, -
5
°
≤
β
(
t
)
≤
5
°. Figure 10 shows a comparison between their actual values and computed values by (12).Figure 10
Values off
k: actual versus computed using TS approach.This leads to eight(
2
3
) dynamic matrices that are written combining the maximum and minimum values of these functions, such that
(25)
A
1
=
A
(
f
1
max
⁡
,
f
2
min
⁡
,
f
3
min
⁡
)
,
A
2
=
A
(
f
1
max
⁡
,
f
2
max
⁡
,
f
3
min
⁡
)
,
A
3
=
A
(
f
1
min
⁡
,
f
2
max
⁡
,
f
3
min
⁡
)
,
A
4
=
A
(
f
1
min
⁡
,
f
2
min
⁡
,
f
3
min
⁡
)
,
A
5
=
A
(
f
1
max
⁡
,
f
2
min
⁡
,
f
3
max
⁡
)
,
A
6
=
A
(
f
1
max
⁡
,
f
2
max
⁡
,
f
3
max
⁡
)
,
A
7
=
A
(
f
1
min
⁡
,
f
2
max
⁡
,
f
3
max
⁡
)
,
A
8
=
A
(
f
1
min
⁡
,
f
2
min
⁡
,
f
3
max
⁡
)
.The results shown in this section were obtained considering two different freeplay amplitudes and a time step equal to 0.001 seconds. To compute the control gain, a column input matrix was used to represent an actuator applying force on the control surface (B
o
c
=
{
0
0
1
}
T). Also, it was considered that the parameter g is equal to 1.0 for all cases. However, in particular for experimental applications where limitations exist on the actuator force, the parameter g can be conveniently chosen between ]0,1
[.For the cases in which the freeplay amplitude was|
δ
|
=
0
.
4
°, the initial conditions were β
(
0
)
=
0
.
9
° and 2
.
5
°. Similarly, for the cases in which |
δ
|
=
0
.
6
°, the initial conditions were β
(
0
)
=
1
.
3
° and 3
.
5
°. These values were defined based on the results obtained from the HBM, as discussed in Section 4.2.Two different conditions of airspeed were considered to demonstrate the method. Different controllers were designed for each case. In first case was chosen a condition of stable limit cycle for each freeplay amplitude. In the second one, unstable limit cycles were chosen in order to evaluate the effectiveness of the approach. Consider if all neighboring trajectories approach the LCO, it is stable; otherwise, the LCO is unstable. Definitions and a detailed discussion about stable and unstable limit cycles can be found in [38].Case 1 (stable LCO).
In the first caseV
=
12.5 m/s was considered for both freeplay amplitudes. Figures 11 and 12 show that the uncontrolled system response exhibits a limit cycle oscillation with first harmonic around 3.7 Hertz and amplitude approximately equal to 1 degree. It is possible to note in those figures that although the first order HBM cannot predict, the system response exhibits a dominant component around 15 Hertz. In addition, the plunge and pitch degrees of freedom behavior can be observed in the phase plan shown in Figures 13 and 14. In these cases of stable LCOs the designed controllers were able to suppress the limit cycles. The appendix presents details involved to compute the Fourier transform when Hénon’s technique is used for the time integration process.Figure 11
Case1: control surface rotation (|
δ
|
=
0
.
4
° and V
=
12.5 m/s).Figure 12
Case1: control surface rotation (|
δ
|
=
0
.
6
° and V
=
12.5 m/s).Figure 13
Case1: phase plan (|
δ
|
=
0
.
4
° and V
=
12.5 m/s).Figure 14
Case: phase plan (|
δ
|
=
0
.
6
° and V
=
12.5 m/s).Case 2 (unstable LCO).
In this case was consideredV
=
13.1 m/s for both freeplay amplitudes. Figures 15 and 16 show that the uncontrolled system response exhibits an unstable limit cycle. However, the controller gains computed using the proposed approach are able to suppress those nonlinear responses. To improve the comprehension, Figure 17 shows a zoom selection area for the control surface rotation considering the second freeplay amplitude.Figure 15
Case2: control surface rotation (|
δ
|
=
0
.
4
° and V
=
13.1 m/s).Figure 16
Case2: control surface rotation (|
δ
|
=
0
.
6
° and V
=
13.1 m/s).Figure 17
Case2: zoom selection area for the control surface rotation (|
δ
|
=
0
.
6
° and V
=
13.1 m/s).For this second case the plunge and pitch degrees of freedom behavior can be observed in the phase plan shown in Figures18 and 19. Finally, Figure 20 shows the control force for δ
F
=
0
.
6
° and similar results were obtained for the other cases. Note that the parameter g
<
1.0 can be chosen in order to decrease the required control force mainly for the first time steps.Figure 18
Case2: phase plan (|
δ
|
=
0
.
4
° and V
=
13.1 m/s).Figure 19
Case2: phase plan (|
δ
|
=
0
.
6
° and V
=
13.1 m/s).Figure 20
Case2: control force (|
δ
|
=
0
.
6
° and V
=
13.1 m/s).


## 5. Conclusions

This paper presents a methodology to control limit cycle oscillations in a typical section airfoil with control surface freeplay. The idea was to use the fuzzy Takagi-Sugeno modeling to describe the nonlinear aeroelastic system through linear matrix inequalities. The closed-loop problem was written using a convex space and a linear control force was computed using convex optimization. A linear flutter analysis was performed to identify the flutter boundary. Also, the first order harmonic balance method was solved to predict the limit cycle oscillation envelope. After these first predictions two freeplay amplitudes were defined and two airspeeds were chosen to demonstrate stable and unstable limit cycles. Finally, numerical simulations presented different nonlinear aeroelastic responses comparing the controlled and uncontrolled system. The results show that fuzzy TS modeling is an efficient tool for solving this kind of problem.



---

*Source: 597827-2014-06-10.xml*