# Designing a Control System for an Airplane Wing Flutter Employing Gas Actuators

Authors: A. R. Ansari; A. R. B. Novinzadeh
Journal: International Journal of Aerospace Engineering
 (2017)

---

## Abstract

The wing flutter is a dynamic instability of a flight vehicle associated with the interaction of aerodynamic, elastic, and inertial forces (aeroelastics phenomena). In this study, just the primary control is investigated. Also, in order to control the two-dimensional wing flutter, the force jet and pulse width pulse frequency (PWPF) are suggested. The PWPF modulator has the advantage of almost linear operation, low jet gas consumption, flexibility in addressing various needs, and good accuracy in presence of oscillations. This scheme makes use of quasi-steady dynamic premises and incompressible flow, as well as the thin airfoil theory. It should be noted that, to justify the application of the aerodynamic theory, we have speculated that the thruster jet ejected through a nozzle with a diameter smaller than several millimeters has a supersonic regime (with Mach number of the order ofM≈3.5). Consequently, the interference of the thruster jet in the boundary layer, flow, and circulation around the airfoil which are characterized by low speed would be negligible. The operation of the jet as a thruster is handled by the PWPF modulator, and the process output is fed back to the system via a PD controller. In order to control the wing flutter oscillation, the location of installing the actuator on the airfoil is investigated.


---

## Body

## 1. Introduction

Flutter is a result of interaction of aerodynamic, inertial, and elastic forces as is shown in Figure1; discussion and investigation of the phenomenon are related to the field of aeroelasticity [1, 2]. Such an event can occur in somewhat flexible structures, such as airplane wings, telegraph wires, bridges, and buildings [3]. It has not been long since controlling the wing flutter phenomenon has been investigated by aeroelastics researchers. Each has proposed a method to damp or evade flutter.Figure 1
The summarized studies on an aeroelasticity triangle of collar [1].Zhang et al. [4, 5] suggested single- and multiple-input adaptive control solutions for unsteady airfoil while they considered mass unbalanced trailing- and leading-edge flaps in their model; that is, their model captures the dynamics of the attached control surfaces.Utilizing piezoelectric materials and control surfaces is of the suggested methods. Conventional methods for linear and nonlinear control of wing flutter have deployed control surfaces and smart structures such as piezoelectric (PZT) and electrorheological fluids (ERFs) materials. For example, using two control surfaces and feedback linearization techniques, Ko et al. [6] proved global stability of an adaptive control law for an airfoil with nonlinear damping and stiffness operating in unsteady flow.In the innovative scheme proposed in this paper, two-dimensional control of wing flutter is reached via the momentum resulting from the gas exhaled from the jet installed on the wing. The jet fluid used by the actuator can be supplied via the inlet vents in front of the airplane or through the airplanes compressor engine. In this study, ideal air is considered as the fluid operating the jet.In practice, the PWPF modulator decides on the on- and off-time of the jet (thruster) by maintaining the merits of a Bang-Bang controller. Using two thrusters provides the opportunity of changing the direction of force and momentum (acted) on the airfoil. This has been made possible via a command block. Equations solutions and simulations are carried out in MATLAB software. The model and equations governing this system are presented in the following sections.


## 2. Aeroelastic Model and Equations Governing the Two-Dimensional Wing Movement

In aeroelastics, the pitch and plunge oscillations of an airfoil are modeled via linear and torsion spring as illustrated in Figure1. As mentioned before, due to using the thin airfoil theory, the angle of attack (α) is small. Another assumption employed in this scheme is quasi-steady aerodynamics. Dynamic equations of the system are obtained on the basis of Lagrange equations. The dynamic system has two degrees of freedom (i.e., the pitch and the plunge motions). The free stream velocity is V and the positive direction is downward. Accordingly, the upward lift force will be negative. All torque components are considered around the elasticity pivot. The airfoil chord is “c,” and “b” is the mid-chord.In Figure2, Xf is the position of the elasticity pivot, Xc is the position of mass center, and Kh and Kα show the plunge and pitch spring stiffness coefficients, respectively. There e or a is dimensionless distance from the aerodynamic center to the pivot. α is the pitch angle.Figure 2
Two-dimensional model of the airplane wing.For two-dimensional, incompressible flow aerodynamic center is at the airfoil quarter-chord; for supersonic flow it moves back to the half-chord [8].I
α indicates the inertial moment around the elasticity pivot, and the moments are around the elasticity axis. In order to derive equations of aeroelastic airfoil, two sets of equations are coupled, namely, the equations of structural model of airfoil and the equations of the aerodynamic model of the airfoil.In this study, a controlling model in order to indicate the effect of the actuators and monitoring elements is added to the total model. First, the structural equations are derived. The vertical displacement of the structure according to Figure2 is described as in (1)z=h+axf,where xf indicates the place of the elastic center from the leading edge:(2)T=12mh˙2+sh˙α˙+12Iαα˙2.The overall kinetic energy of the structure is [2](3)s=mc2-xf,where(4)Iα=13mc2-3cxf+3xf2.The potential energy for the above scheme is the sum of all energy saved in the spring:(5)V=12khh2+12kαα2.The sum of kinetic and potential energies is the total energy of the system. Here, in order to find the motion equations, the total energy is placed in Lagrange equation [9]:(6)ddt∂T-V∂q˙+∂T-V∂q=Q,q=h,αT.The result of Lagrange equation is the following matrix equation:(7)Mq¨+kq=Q,where “Q” is the resultant vector of external forces (i.e., thruster and aerodynamic forces). The resultant equations of aerodynamic model of airfoil depend upon flow regimes. Here, since the fluid is incompressible, the Mach number is not included. The instantaneous aerodynamic forces on the airfoil depend upon both position and strength of the wake vertices (unsteady aerodynamics), which means the instantaneous aerodynamic forces rest not solely on the current motion of the airfoil, but also upon its past motion from the beginning of the movement. In quasi-steady aerodynamic forces, the simplest, possible modeling includes ignoring the effects of the wake. From the thin airfoil theory for the lift coefficient, the following equation is obtained [2]:(8)cl=2πA0+A12,(9)dzdx=xf-xα˙U.From the thin airfoil theory for moment and lift coefficients and (9), the following equations are obtained [2]:(10)clt=2πα+h˙U+34c-xfα˙U,cmxft=-cπ8Uα˙t+xfc-14clt.The moments are referenced around the elasticity pivot. Hence, the contribution of the lift to the distance from the leading edge (Xf) is added to moment around the flexural axis. As well as lift and moment forces, the force of the virtual mass due to the airfoil acceleration is also exerted to the airfoil. This force affects both the lift and the moment forces.In case of diverging the wing oscillations at a certain speed, the activation of the actuator (thruster), force(Ftjet), and moment (Mtjet) generated by the thruster in two directions (Figure 3) leaves a damping impact on the force and moment. The total motion equations are presented in (11) and (12):(11)Nt=Lt+Ftjet=ρπb2h¨-xf-c2α¨+ρπb2Uα˙+ρπU2cα+h˙U+34c-xfα˙U+Ftjet.Figure 3
Two-dimensional model of an airplane wing, equipped with an actuator capable of on-off operation (positive direction is downward).For the overall moment,(12)Mt=Mxft+Mtjet=xf-c2ρπb2h¨-xf-c2α¨-ρπb28α¨-34c-xfρπb2Uα˙+ρU2ec2πα+h˙U+34c-xfα˙U-116ρUc3πα˙+Mtjet.In (11) and (12) equations, L(t) and Mxf(t) are, respectively, the lift and overall moment, created by aerodynamic forces and added mass, referenced around the elasticity pivot.By substituting (2), (5), (11), and (12) in (7) and after simplification, the aeroelastic equations are derived:(13)mssIαh¨α¨+ρπb21c2-xfc2-xfc2-xf2+b28h¨α¨+ρUcπ134c-xf+c4-ecc2-xf2+34c-xfc4h˙α˙+kh00kαhα+ρU2cπ000-echα=FjetMjet.Equations in (13) are the second-order linear differential equations. In these equations, there are the stiffness, damping, and mass matrices both structurally and aerodynamically. As denoted, the right side of the equations in (13) is nonzero. These equations are heterogeneous.


## 3. Gas Actuator

Controlling the wing oscillations is employed via a jet, as an actuator, to control flutter and damp oscillations. The momentum resulting from this jet affects wing oscillations amplitude. In the simulations, activating the thruster in both up and down directions is programmed. The flow rate and pressure of cold-gas system are provided by the vents in the front of the airplane or through the airplane's engine compressor. The actuator is provoked as the relative speed of the motion reaches the flutter speed. The influence of the thruster on the wing surface contains force and moment. In case the jet resides precisely at the elasticity pivot, it will not have a moment around the pivot. The distance from the jet to the leading edge is (xj), Figure 4, which is a percentage of the chord (C).Figure 4
Three-dimensional perspective of the wing and the placement and operation of the thruster in up direction (symmetrically employed in down direction in appropriate intervals).The force acting upon the wing by the thruster is presented in (14), which is a function of speed, mass flow rate, and the difference between ambience and discharge pressures, as well as jet cross section. The amount of the thrust resulting from the jet equals [10] (14)Fjet=Vem˙e+AePe-Pa.Here,Ve indicates emitted fluid speed, m˙e is the mass flow rate, Pe is the emitted pressure, Pa denotes the ambient pressure, and Ae is the cross section of the jet nuzzle. The acting thruster fluid is the ideal air.The kind of jet fluid in order to produce thrust is of importance. However, the second important characteristic parameter of the thrust source is the Specific Impulse (Isp), indicating the amount of fuel mass which can be converted to thrust stroke. According to definition, the Specific Impulse is(15)Isp=Fjetm˙g,Fjet=Isp·m˙g.The moment produced by the thruster is expressed by (16). All moments are referenced around the elasticity pivot:(16)Mjet=xjet×Fjet.


## 4. Modulator

PWPF modulator is a pulse to pulse frequency-type thruster control system. This type of modulator consists of a two-variable Schmitt-trigger circuit, a first-order filter, and a feedback loop. In such a modulator, the filter precedes Schmitt-trigger, not in the control feedback (Figure5). This modulator produces a sequence of pulses for the thruster through frequency- and width-modulated pulses [6]. The modulator has its merits over other thruster control systems. It also enjoys all the advantages of the Bang-Bang control system. Moreover, it is one of the popular systems of position control in the space.Figure 5
PWPF modulator circuit.In addition, the mentioned modulator is chosen as a candidate for being used in satellite system of European Student Moon Orbiter (ESMO) [11]. Static analysis indicates that the modulator, for a fixed input, has a quasi-linear (almost linear) performance. In order to analyze the stability, a dynamic analysis for a nonlinear input is required, and one can make use of the approach of fractional describing function [12]. In the present paper, the jet performance is controlled via the modulator; that is, the modulator, in appropriate intervals, provides us with the situation to monitor the jet performance (on and off). This modulator runs like digital converters.The jet (or, if needed, “jets”) will operate in case the flight speed reaches the critical speed. In this speed, known as “flutter speed,” the wing oscillations divergence is controlled via the thruster. By the operation of the thruster which is modeled as two-directional one (up and down), the wing oscillation amplitude is controlled in two directions.According to Figure5, the four parameters of this modulator are simply recognized, namely, km, Uon, Uoff, and Tm. km, indicating the gain of the filter, and Tm is the time constant of the filter. Uon and Uoff depict the turn-on and turn-off values of the Schmitt-trigger, respectively. In addition, h is the hysteresis band, which is derived by Uon-Uoff=h.The Schmitt-trigger, as illustrated in Figure5, is an on-off relay composed of a dead-zone and a hysteresis band. When the input to the Schmitt-trigger “f(t)” is greater than Uon, the thruster output is Um (i.e., the thruster switches on). If its input “f(t)” is smaller than Uoff, the output is zero, and when input is in order between Ton and Toff, the thruster sustains its on-state and keeps on with its operation. The system error is the difference between the Schmitt-trigger output Um and the system reference r(t), Figure 6 [13].Figure 6
PWPF modulator response tot fixed input C=0.5.


### 4.1. On- and Off-Time of the Modulator

For a fixed reference signal in Figure5, the error e(t) and the filter output f(t) are correlated as(17)et=rt-ut,(18)ft=kmet1-exp⁡-tTm+f0exp⁡-tTm,where f(0) is the initial condition of the filter. When the reference signal is a step function with amplitude r(t)=C, the modulator output is also a step function u(t)=U (Figure 6). Within the on-time period of the thruster, the filter output falls asymptotically. At the turn-on instant, the filter output reaches Uon; hence, the initial condition is f(0)=Uon.AfterTon passes, the filter output f(t) will stop at Uoff. Substituting these values in (18), one can find the thruster on-time:(19)Ton=-Tmln⁡1-hUon-kmC-U.The modulator on-time is referred to as pulse width (PW) [14]. With analogous procedure, the modulator off-time is obtained as(20)Toff=-Tmln⁡1-hkmC-Uoff.The modulator output frequency can also be expressed as(21)f=1Ton+Toffand the modulator duty cycle as(22)DC=TonTon+Toff.The duty cycle, also referred to as an adjustment factor, is a measure to determine how the modulator responds to an input; for example, to what extent the output tracks the reference input signal [13].


## 4.1. On- and Off-Time of the Modulator

For a fixed reference signal in Figure5, the error e(t) and the filter output f(t) are correlated as(17)et=rt-ut,(18)ft=kmet1-exp⁡-tTm+f0exp⁡-tTm,where f(0) is the initial condition of the filter. When the reference signal is a step function with amplitude r(t)=C, the modulator output is also a step function u(t)=U (Figure 6). Within the on-time period of the thruster, the filter output falls asymptotically. At the turn-on instant, the filter output reaches Uon; hence, the initial condition is f(0)=Uon.AfterTon passes, the filter output f(t) will stop at Uoff. Substituting these values in (18), one can find the thruster on-time:(19)Ton=-Tmln⁡1-hUon-kmC-U.The modulator on-time is referred to as pulse width (PW) [14]. With analogous procedure, the modulator off-time is obtained as(20)Toff=-Tmln⁡1-hkmC-Uoff.The modulator output frequency can also be expressed as(21)f=1Ton+Toffand the modulator duty cycle as(22)DC=TonTon+Toff.The duty cycle, also referred to as an adjustment factor, is a measure to determine how the modulator responds to an input; for example, to what extent the output tracks the reference input signal [13].


## 5. Results and Comparison

For a two-dimensional wing with the following features (Table1), simulations are conducted in two ways: (1) without control at flutter speed; (2) with the thruster enabled by the modulator at the same speed. The simulated block diagram of this process is illustrated in Figure 7. The simulations are carried out by MATLAB software.Table 1
Parameters for a thin symmetrical airfoil.
ρ
1.2254 kgs/m3
I
α
6.1272 kgs⋅m2
k
α
2563.83 N-m/Rad
k
h
1962.433 N/m
m
19.624 kgs
C
=
2
b
1.829 m
x
f
0
.
4
∗
c m
a
−0.2Figure 7
MATLAB block diagram of a two-dimensional wing controller with the modulator.One of the important properties of the modulator is minimizing the acting thruster fluid and its appropriate performance in the case of oscillations. In this paper, merely the control aspect is discussed.In the first case, the thruster is disabled/off and there is no control over the wing. After the program is performed, the system condition is investigated. As mentioned earlier, the quasi-steady equations are applied and the results obtained from the same equations. However, the eigenvalues of the matrix describe the behavior of the system. Natural frequencies and damping coefficients can be calculated by them.Figure8 depicts the damping variations versus speed. At the moment when the damping coefficient gets positive, the flutter speed is measured, which is approximately UFlutter=29.97m/s. The pitch and plunge oscillation amplitudes are displayed in Figure 9. Divergence in both cases occurs after the flutter speed.Figure 8
Variations in two damping coefficients and pitch and plunge frequencies versus speed derived from quasi-steady aerodynamic equations. Flutter speed is aboutUF=29.97 m/s.Figure 9
Open-loop responses (limit cycle) in pitch and plunge oscillation amplitude versus time for two-dimensional wing, in degree and meter, respectively.The program data is offered in Table1. The initial conditions for the wing are expressed in (23)x0=α,h,α˙,h˙T=0.07,0.02,0,0T.In the second case, the system is controlled in order to prevent the structure against damage. The process data is again derived from Table1. Also, modulator parameters are km=16, Tm=0.15, Uon=0.45, h=0.2, and Um=1 in this study. In addition, (PD) controller gains are kp=5000 and kd=3000.According to Figure7, the jet and control blocks are enabled in this state. The time range of jet operation and fluid consumption rate are related to the modulator parameters.Employing the modulator data and applying the force and moment to the right side of aeroelasticity equations (13), together with the help of the PD controller, the divergence of oscillation amplitude is controlled.From (14), the amount of employed force from the thruster in this scheme is Fjet=37.6 N and is placed at the distance of xjet=0.98∗C from the chord. Due to the assumption of the two-dimensional wing (airfoil), the thruster can alter its location just on the x-axis.Figure10 shows controlling the divergence of the oscillation amplitude. In this picture, oscillations in both pitch and plunge are reduced and controlled. This procedure has been employed in a similar paper [7] in which the flutter occurrence of a two-dimensional wing has been investigated via applying a control surface.Figure 10
Closed-loop responses (controlled flutter) in pitch and plunge oscillation amplitude in two- dimensional wing, versus time (in degree and meter, resp.).The results are thoroughly satisfactory in order to control the flutter. As it can be noticed in Figure11, the two-dimensional wing at the flutter speed of UFlutter=15 m/s displays periodic pitch and plunge oscillations amplitude.Figure 11
Results from [7] for a small two-dimensional wing with cycle-limited asymptotic oscillations.Through activating the thruster at this speed, the oscillations amplitude tends to zero att=5 s, Figure 12. The demanded thrust is Fjet=12.45 N and is employed at Xjet=0.10∗C from the leading edge. Also, modulator parameters are km=20, Tm=0.2, Uon=0.4, h=0.2, and Um=1 in this study. In addition, (PD) controller gains are kp=800 and kd=350.Figure 12
Oscillation control of a two-dimensional wing via PWPF modulator for a short time, fromt = 0 s and t = 10 s (according to Table 2).Table 2
Parameters for a thin symmetrical airfoil from [7].
ρ
1.2254 kgs/m3
I
α
0.1065 kgs⋅m2
k
α
Nonlinear N-m/Rad
k
h
2844.4 N/m
m
12.287 kgs
C
=
2
b
0.27 m
x
f
0
.
4
∗
c m
a
-
0
.
6The characteristics of this wing are stated in Table2. The spring stiffness coefficient is decided on by [15](24)kα=2.82∗1-22.1∗α+1315.5∗α2-8580∗α3+17289.7∗α4.The initial conditions for this wing are expressed in(25)x0=α,h,α˙,h˙T=0.2,0.02,0,0T.



---

*Source: 4209619-2017-10-01.xml*