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.
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 Figure 1 with a trailing edge control surface is normally used to represent an aeroelastic system with three degrees of freedom that includes pitch , plunge , and control surface rotation . This modeling was previously proposed in [22]. The structural properties of this system are represented by the springs , , and , 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].
The nonlinear discontinuity (freeplay) is considered to occur in the control surface spring . So that the equation of motion for this system can be represented by where and are, respectively, the structural mass and damping matrices, is a matrix of aerodynamic coefficients that depend on the airfoil geometry, is the dynamic pressure, and is the air density and is the airspeed. is a matrix of input (forces or moments). The vector represents the physical displacements and is the control force. The vector represents the elastic restoring moment which depends on the control surface restoring moment shown in Figure 2.
By considering a freeplay amplitude of , the elastic restoring moment can be written as
The expressions in (2) are rearranged in a matrix form and the nonlinear function is introduced such that the stiffness is defined as a nonlinear structural stiffness matrix given by where
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 matrix 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 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 where is the Laplace variable, is the number of lag terms, is the th lag parameter , and 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: where is the state vector and are states of lags required for the approximation of matrix . The matrix of outputs has dimension , where and are, respectively, the velocity and displacement output matrices and the submatrix has only zeros to complete the matrix dimension, and is the number of degrees of freedom. Matrix is presented in the following form: where where and are, respectively, the aeroelastic mass and damping matrices, is the nonlinear aeroelastic stiffness matrix, is the input matrix, is an identity matrix, and 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 Figure 2. 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: where , , , and represents the th nonlinear function, where .
The nonlinear system described by (9) can be represented by the Takagi-Sugeno model using the following rule: where is the fuzzy set and is the number of model rules; are known premise variables that in general may be functions of the state variables, external disturbances, and/or times. Each linear model represented by 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 as a linear combination of its maximum and minimum values both given, respectively, by
From these maximum and minimum values, each nonlinear function can be represented by where
If , , then each nonlinear function can be conveniently written as or where and is given by and the superscript indicates the combination between maximum and minimum values for each th function .
Considering that and substituting each according to (15), the nonlinear equation of motion (17) is rewritten after some rearrangement as where each linear matrix is the original nonlinear matrix at which the functions are substituted by the combination of maximum and minimum values and . In this case .
Particularly for the typical section including freeplay nonlinearity, there are three nonlinear functions into the matrix and they are given by where indicates an element into matrix and the subscripts indicate the th row and th column in each respective matrix. These functions depend on which is illustrated in Figure 3.
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 by , where is the feedback gain matrix. The feedback gain is calculated using linear matrix inequality to solve Lyapunov’s function for stability, in this case where is the Lyapunov function. The stability of this system is assured if there is a positive-definite matrix such that the inequality of (20) is true. After some rearrangement this equation can be rewritten such as [32] and finally, the solution of this inequality is equivalent to the solution of all the following inequalities satisfied simultaneously [34]: Inequality (22) is not linear: then considering , , it is possible to write the following linear matrix inequality:
Although the control gain has been computed using inequality (23), the required control force can exceed desired limits. In this case, the control gain is multiplied by a constant between the interval , conveniently adjusted during the control design. The control force is rewritten as
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
Figure 1 illustrates the model and its physical and geometric properties are presented in Table 1. Figures 4 and 5 show, respectively, the classical - and - diagrams for the linear flutter solution. According to these results flutter speed is equal to 12.7 m/s.
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 Figures 6 and 7 were obtained by extracting the eigenvalues from the matrix defined for different values of equivalent stiffness , that is, the control surface stiffness assuming values from zero to . The values of are shown in Figure 8, where the nonlinear function assumed a unitary value and the control force . The HBM also provided an estimate for the first harmonic of the system response according to Figure 9.
4.3. Controller Design
Using the methodology described in Section 4, three nonlinear functions were used to describe the aeroelastic system with freeplay (, ). These functions were computed assuming the maximum oscillation amplitude equal to five degrees; that is, . Figure 10 shows a comparison between their actual values and computed values by (12).
This leads to eight dynamic matrices that are written combining the maximum and minimum values of these functions, such that
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 (). Also, it was considered that the parameter is equal to for all cases. However, in particular for experimental applications where limitations exist on the actuator force, the parameter can be conveniently chosen between ].
For the cases in which the freeplay amplitude was , the initial conditions were and . Similarly, for the cases in which , the initial conditions were and . 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 case 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.
Case 2 (unstable LCO). In this case was considered 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.
For this second case the plunge and pitch degrees of freedom behavior can be observed in the phase plan shown in Figures 18 and 19. Finally, Figure 20 shows the control force for and similar results were obtained for the other cases. Note that the parameter can be chosen in order to decrease the required control force mainly for the first time steps.
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.
Appendix
The integration process involving Hénon’s technique can result in different time steps, as discussed in Section 2.2 and shown in Figure 21. However, in order to apply the Fourier transform to compute the system responses in frequency domain, it is necessary to keep a constant time step. In this case, after using the amount to identify the transition point shown in Figure 2, another time integration can be performed considering the time step . Then, the system response for the constant time step cannot include that preliminary point obtained using . This procedure was employed to compute the responses plotted in Figures 11 and 12.
Conflict of Interests
The authors declare that the research was conducted in the absence of any commercial or financial relationships that could be construed as a potential conflict of interests.