Control Surface Freeplay Effects Investigation on Airfoil’s Aero-Elastic Behavior in the Sub-Sonic Regime

Abstract

One of the main limitations of linearity assumptions in airfoil’s aero-elastic problems is the inability to predict the system behavior after starting the instability. In reality, nonlinearities may prevent the amplitudes from going to infinity. This paper presents a methodological approach for predicting airfoil aero-elastic behavior to investigate the control surface freeplay effects on the state responses and the flutter speed. For this purpose, the airfoil structural model is firstly developed while using the Lagrange’s method. The aerodynamic model is then generated by utilizing the Theodorsen approach for lift and moment calculation and Jones approximation with P-method for unstable aerodynamic modelling. After that, the aero-elastic model is developed by combination of structural and aerodynamic models and a numerical integration method is used to extract the time responses in the state space. The flutter analysis has been completed by utilizing the P-method for the system without freeplay and by the time response approach for the system with freeplay. The results that were obtained from simulations confirm the effectiveness of the proposed method to predict the aero-elastic behavior and stability condition of a two-dimensional airfoil as well as to estimate the flutter speed with reasonable accuracy and low computational effort. Moreover, a sensitivity analysis of freeplay degree on time response results has been done and the results are discussed in detail. It is also showed that the control surface freeplay decreases the flutter speed. The results of the paper are also validated against publicly available data.

1. Introduction

The fluid-structure interactions (FSI) study is one of the most important steps in the design procedure of deformable flying objects. These interactions could be stable or oscillatory. In oscillatory interactions, the strain that is induced in the solid structure causes it to move, such that the source of strain is reduced, and the structure returns to its former state only for the process to repeat. The aeroelastic problems did not occur if aircraft structures were rigid. However, modern aircraft structures are flexible with respect to weight considerations. Aeroelastic instabilities involve aerodynamic forces that are generated due to the motion of the structure, and may cause structural deflection. The aerodynamic effects, such as buffeting, refer to the fluctuating nature of the wing (and its interaction with the structure) and sometimes cause aerodynamic instability. When buffeting, flutter instability, and etc. occurred, aerodynamic and aeroelastic effects play a major role in the design of flexible structures [1]. The aerodynamic and aeroelastic effects will cause nonlinearities in the system that result in different dynamic behaviors, like limit cycle oscillations (LCO), infinite increase in the amplitude, and chaotic behavior. Therefore, the study of these effects plays a vital role in the design and analysis of the flying objects. The effect of laminate lay-up on flutter speed of composite missile stabilizers for the subsonic flight conditions is comprehensively investigated in [2] based on a two-dimensional (2D) flutter model. The authors have used a finite element approach based on the Lanczos method to achieve optimum results in respect to solution accuracy and computational cost.

One of the main limitations of the linearity assumption is the inability to predict the system behavior after starting the instability. In reality, some kind of nonlinearities may prevent the amplitudes from going to infinity. In nonlinear systems, vibration with limited amplitude or chaotic motion occurs in upper or lower of the flutter speed—the speed that the response of domain does not change with time [3]. Despite the fact that these two types of motions do not lead to destructive failure in structures, these can lead to fatigue failure or destructive effect on control systems. Therefore, the effect of nonlinear phenomena on aerodynamics behavior should be considered in the design process of flying objects.

Generally, the instability in aero-elastic systems can be traced from aerodynamic or structure sources. The nonlinear terms of aerodynamic can be the effects of shock motion, high angle of attack, flow separation, transonic flow, to name but a few. Besides, the nonlinear terms of structures can be inclusive of freeplay between components or control systems, large deformations, material nonlinearity, the effects of friction in control cables or bearings, the limitations in kinetic energy of control surfaces, and concentrated non-linear terms that may be related to elastic deformation in connections.

Kholodar explained why a certain amount of the preload or angle of attack could be the main reason for the unexpected higher-frequency oscillations. This study shows that a certain preload is required to excite limit cycle oscillations [4]. Marsden and Price have done an experimental analysis to show that, when freeplay nonlinearity is introduced in the pitch degree of freedom, the flutter speed that is estimated by extrapolation of the damping curve decreases with increasing the freeplay size [5]. A free vibration analysis is performed while using the finite element, and the fictitious mass methods in [6] in which aero-elastic responses included LCO, unstable LCO, and periodic motion are observed for air speeds below the linear flutter boundary. Moosazadeh, et al. assumed a fully nonlinear structure with nonlinear third order piston theory to show that the stability increases with freeplay in high speed and increasing the velocity will decrease the damping effect in post flutter behavior [7]. The challenges and characteristics of the flow past the body at critical angles of attack based on modern aerodynamics explained in [8]. The author showed that the separation starts from the trailing edge and moves towards the leading edge, proportionally to the increase of the angle of attack, although critical angles of attack depends on the type of the airfoil, the local characteristics for the defined conditions of the flow, and the dynamics of the angle of attack. The best actual lift-curve slope of the conventional symmetrical NACA 0012 airfoil according to the Mach and Reynolds flow numbers is obtained in [9]. Rasuo analyzed the principal factors that influence the accuracy of two-dimensional wind tunnel test results [10], and the influences of Reynolds number and Mach number on two-dimensional wind-tunnel testing [11]. Ocokoljić, et al. investigated the requirements for the accuracy of measurements in a wind tunnel, such as calibration of the wind tunnel test section [12]. Other relevant studies have also focused on stall-induced vibration analysis [13,14], bifurcation and chaos analysis [15,16,17], non-linear dynamics [18], and stability analysis [19]. It could be concluded that one of the main objectives of all studies in this field is to develop a reliable, precise, and fast method for aero-elastic behavior, stability, and the flutter speed analysis while avoiding numerical integration issues of instability and inaccuracy.

Table 1. Research milestones in airfoil aero-elastic behavior analysis.

	Milestone	Topic/Main Achievement
1	Theodorsen, 1934 [20]	Closed-form solution for the problem of an unsteady aerodynamic load on an oscillating airfoil.
2	Yang, et al., 1980 [21]	Flutter analysis of a NACA 64A006 airfoil in small disturbance transonic flow
3	Leishman and Nguyen, 1989 [22] Kurniawan, 2013 [23]	State-space representation of unsteady airfoil behavior
4	Menon and Mittal, 2018 [24]	Computational modeling and analysis of aero-elastic wing flutter

lists the main milestones in this field with their advantages and disadvantages.

Gap Analysis

• The first step by Theodorsen was limited by this assumption that the harmonic oscillations in inviscid and incompressible flow is subject to small disturbances.
• Woolston also just studied a simple aero-elastic system by investigating the effect of freeplay on flutter velocity [25].
• Air Force Wright Aeronautical Laboratories focused on flutter analysis of a NACA 64A006 airfoil (as a case study) pitching and plunging in small-disturbance, un-steady transonic flow.
• The steady-state formulations has been used to tackle the indicial lift on a thick-aerofoil response method in the literature. It is practical but still limited.
• The method presented by Menon and Mittal is precise and comprehensive (they have used ViCar3D code [26,27]) but the huge computational effort is the weakness of this study.

The above analysis confirms the lack of a fast and precise method for predicting and analyzing the effects of freeplay on control surface on the aero-elastic behavior of airfoil. In this paper, a new methodological approach is proposed to predict the aero-elastic behavior, stability, and the flutter speed of an airfoil with very low computational effort and reasonable accuracy in comparison with the published studies. For this purpose, the structural model is developed first while using Lagrange’s equation. Subsequently, using Theodorsen approach, the aerodynamic model is generated. The unstable aerodynamic is also modeled while using Jones approximation. Next, the aero-elastic model is developed by combination of structural and aerodynamic models. After that, the numerical results are generated by numerical integration in state space. Finally, the results analysis, including the flutter speed prediction, stability analysis, sensitivity analysis on freeplay degree, and the effect of freeplay on the time responses are presented and discussed.

2. Equations of Motion and the Proposed Method

For airfoils with high aspect ratio, a two-dimensional modelling approach with equivalent section properly simulates the aero-elastic behavior of 3-D airfoil. The structure is modeled with mass-spring and the aerodynamic is modeled with two dimensional assumptions. Mass-spring structural models have proper accuracy in the most two-dimensional nonlinear structural problems.

The general form of the equations of motion for an elastic system in Laplace domain can be represented, as follows [28]:

$$\left\{\frac{U^2}{b^2}{Mp}^2+\frac{U}{b}Cp+K+qQ\left(p\right)\right\}\tilde{q}=0$$

(1)

where $U$ is free stream velocity, $q=\raisebox{1ex}{${\rho U}^2$}\!\left/ \!\raisebox{-1ex}{$2$}\right.$ is dynamic pressure, and b is reference length that equal to half of chord of two dimensional airfoil. $M$, $C$, and $K$ are mass, damping, and stiffness matrixes, and $\tilde{q}$ is generalized coordinates vector. The $Q\left(p\right)$ matrix in the generalized aerodynamic forces is defined, as follows:

$$f_{\mathrm{aero}}=q_{\infty }Q\left(p\right)\tilde{q}$$

(2)

And $p$ is below dimensionless variable Laplace transform:

$$p=g+\mathrm{ik}$$

(3)

where $g$ is equal to $\gamma k$, $k$ is equal to $\raisebox{1ex}{$b\omega $}\!\left/ \!\raisebox{-1ex}{$U$}\right.$, $\omega$ is oscillatory frequency, and $\gamma$ is the rate of decline in transient motion.

There are different methods to obtain eigenvalues of Equation (1) and analyze the instability of system behaviors. In this study, the unsteady aerodynamic and p-method is used because of having more general applications than other method.

Table 2. Proposed method steps.

1	Structural modelling
	• Lagrange’s Equation • Freeplay Model
2	Aerodynamic modelling
	• Theodorsen approach for lift, aerodynamic moment and restoring moment calculation • Jones approximation, Fourier transformation, Wagner function, and Duhamel’s integral for unstable aerodynamic
3	Aero-elastic modelling
	• Combination of structural and aerodynamic models
4	Numerical solution
	• Numerical integral using Runge-Kutta technique • Approach validation with flutter speed concept
5	Freeplay investigation
	• Numerical simulation in the presence of freeplay • Sensitivity analysis of the freeplay degree • Validation with published studies

shows the itemized steps of the proposed method. These steps are described in detail in the next sections.

2.1. Two-Dimensional Structural Modeling of Airfoil

In this paper, the Lagrange’s equations are used in the modeling of structural system. These equations are the form of equations of motion that includes kinetic and potential energy, and have many applications, especially in the analysis of dynamical systems. The general form of Lagrange equation is as follows:

$$-\frac{d}{dt}\frac{\partial (T-U)}{\partial {\dot{q}}_i}+\frac{\partial (T-U)}{\partial q_i}+Q_i=0i=1\dots n$$

(4)

where $T=T({\dot{q}}_i,q_i,t)$ is kinetic energy, $U=U({\dot{q}}_i,q_i,t)$ is potential energy, $q_i$ are independent generalized coordinates of system, $Q_i$ are the vectors of generalized forces that include conservative and non-conservative forces, whereas this forces achieved by principle of virtual work, and n is the number of independent degrees of freedom of the system. The difference between kinetic and potential energy is called the system’s Lagrangian. Kinetic and potential energy parts of the Equation (4) will both be calculated in the next sections.

2.1.1. Kinetic Energy of Two Dimensional Airfoil with Control Surface

First, the coordinate of each point of airfoil is explained to calculate the velocity. According to Figure 1, the location of each point in the X-Z coordinates system that origin of this coordinate system located in elastic axis is as follows:

$$\overrightarrow{r}=z\alpha (t)\overrightarrow{i}+[h(t)-x\alpha (t)]\overrightarrow{k}$$

(5)

Additionally, coordinates of each point of the control surface is calculated, as follows:

$$\overrightarrow{r^{\prime }}=[(z-z_{\mathrm{hing}})\beta (t)+z\alpha (t)]\overrightarrow{i}+[-x\alpha (t)-(x-x_{\mathrm{hing}})\beta (t)+h(t)]\overrightarrow{k}$$

(6)

By taking the first time derivative of Equation (6), the velocity of each point is obtained, so the kinetic energy of airfoil with control surface given, as follows:

$$\begin{array}{ll}K& =\frac{1}{2}{\displaystyle \underset{S_a}{\iint }\rho [{(\frac{dh}{dt}-x\frac{d\alpha }{dt})}^2+z^2{(\frac{d\alpha }{dt})}^2]dxdz}\\ & +\frac{1}{2}{\displaystyle \underset{S_{cs}}{\iint }\rho [{\left[(z-z_{\mathrm{hing}})\frac{d\beta }{dt}+z\frac{d\alpha }{dt}\right]}^2}\\ & +{\left[-x\frac{d\alpha }{dt}-(x-x_{\mathrm{hing}})\frac{d\beta }{dt}+\frac{dh}{dt}\right]}^2]dxdz\end{array}$$

(7)

2.1.2. Potential Energy of Two Dimensional Airfoil with Control Surface

According to Figure 1, based on [29], the potential energy of airfoil is composed of two parts: The first part is derived from displacement in Z-direction that modeled with spring $K_h$, and the second part is torsional potential energy of airfoil that modeled with torsional spring $K_{\alpha }$. Additionally, the potential energy of control surface derived from torsion of control surface about hinge axis that is stored in torsional spring $K_{\beta }$. Accordingly, the potential energy of airfoil with control surface is calculated, as follows:

$$U=\frac{1}{2}(K_h{h}^2+K_{\alpha }{\alpha }^2+K_{\beta }{\beta }^2)$$

(8)

2.2. Control Surface Freeplay Modelling

It is assumed that the control surface has freeplay. This eventuates to this assumption that the torsional spring (connection of the control surface to the airfoil) cannot carry forces for some angles of control surface’s displacement. Hence, spring stiffness would be a nonlinear function of angle of torsion. This non-linear behavior is shown in Figure 2 that approximated into multi-sections.

Hence, restoring moment of control surface $M_{S\beta }$ has been shown in (9):

$$M_{S\beta }=\left\{\begin{array}{cc}{k}_{\beta }(\beta -\delta )& \left|\beta \right|\ge \delta ,\beta >0\\ k_{\beta }(\beta +\delta )& \left|\beta \right|\ge \delta ,\beta <0\\ 0& \mathrm{otherwise}\end{array}\right\}$$

(9)

3. The Aerodynamics Model

The aerodynamic forces that are caused by unsteady flow for airfoil are calculated by Theodorsen in [20]. For airfoil displacement modelling in Z-direction, airfoil torsion, and control surface torsion three springs $K_h$, $K_{\alpha }$, and $K_{\beta }$ are used, as shown in Figure 1. In this model, chord length of airfoil is 2 b and elastic axis located at point that airfoil rotates about it. The center of mass from this airfoil is located at $x_{\alpha }b$ distance from elastic axis and the hinge axis of control surface is located at $\left(a+c\right)b$ distance from elastic axis. In addition, the center of mass of control surface considered at $x_{\beta }b$ distance from hinge axis.

Theodorsen calculated lift $L$, aerodynamic moment $M_{\alpha }$, and restoring moment $M_{\beta }$ with the assumption of complete oscillatory motion, as follows:

$$\begin{array}{ll}L=& -\rho b^2(v\pi \dot{\alpha }+\pi \ddot{h}-\pi ba\ddot{\alpha }-v{T}_4\dot{\beta }-T_{1}b\ddot{\beta })-\\ & 2\pi \rho vbC(k)\left[v\alpha +\dot{h}+b(\frac{1}{2}-a)\dot{\alpha }+\frac{1}{\pi }{T}_{10}v\beta +b\frac{1}{2\pi }{T}_{11}\dot{\beta }\right]\end{array}$$

(10)

$$\begin{array}{l}{M}_{\alpha }=-\rho b^2(\pi (\frac{1}{2}-a)vb\dot{\alpha }+\pi b^2(\frac{1}{8}+a^2)\ddot{\alpha }\\ +(T_4+T_{10})v^2\beta +\left[T_1-T_8-(c-a)T_4+\frac{1}{2}{T}_{11}\right]vb\dot{\beta }\\ -\left[T_7+\left(c-a\right)T_1\right]b^2\ddot{\beta }-a\pi b\ddot{h})\\ +2\rho v{b}^2\pi (a+\frac{1}{2})C(k)[v\alpha +\dot{h}+b(\frac{1}{2}-a)\dot{\alpha }+\frac{1}{\pi }{T}_{10}v\beta +b\frac{1}{2\pi }{T}_{11}\dot{\beta }]\end{array}$$

(11)

$$\begin{array}{l}{M}_{\beta }=-\rho b^2([-2T_9-T_1+T_4(a-\frac{1}{2})]vb\dot{\alpha }\\ +2T_{13}{b}^2\ddot{\alpha }+\frac{1}{\pi }{v}^2\beta (T_5-T_4{T}_{10})\\ -\frac{1}{2\pi }vb\dot{\beta }{T}_4{T}_{11}-\frac{1}{\pi }{T}_3{b}^2\ddot{\beta }-T_{1}b\ddot{h})\\ -\rho v{b}^2{T}_{12}C(k)[v\alpha +\dot{h}+b(\frac{1}{2}-a)\dot{\alpha }+\frac{1}{\pi }{T}_{10}v\beta +b\frac{1}{2\pi }{T}_{11}\dot{\beta }]\end{array}$$

(12)

where ρ is flow density and $v$ is flow velocity. Additionally, h, $\alpha$, and $\beta$ are displacement in Z-direction, airfoil pitch, and pitch degree of control surface. Positive directions of X, Z, $\alpha$, and $\beta$ are shown in Figure 1. The origin of coordinate system XYZ located at elastic axis. The constants $T_1$ to $T_{14}$ that are used in Equations (10)–(12) are as follows:

$$\begin{array}{l}{T}_1=-\frac{1}{3}\sqrt{1-c^2}(2+c^2)+c{\cos }^{-1}c\\ T_2=c(1-c^2)-\sqrt{1-c^2}(1+c^2){\cos }^{-1}c+c({\cos }^{-1}c{)}^2\\ T_3=-(\frac{1}{8}+c^2)({\cos }^{-1}c{)}^2+\frac{1}{4}(c\sqrt{1-c^2}{\cos }^{-1}c(7+2c^2)-\frac{1}{8}(1-c^2)(5c^2+4)\\ T_4=-{\cos }^{-1}c+c\sqrt{1-c^2}\\ T_5=-(1-c^2)-({\cos }^{-1}c{)}^2+2c\sqrt{1-c^2}{\cos }^{-1}c\\ T_6=T_2\\ T_7=-(\frac{1}{8}+c^2){\cos }^{-1}c+\frac{1}{8}c\sqrt{1-c^2}(7+2c^2)\\ T_8=-\frac{1}{3}\sqrt{1-c^2}(2c^2+1)+c{\cos }^{-1}c\\ T_9=\frac{1}{2}\left[\frac{1}{3}{\left(\sqrt{1-c^2}\right)}^3+a{T}_4\right]\\ T_{10}=\sqrt{1-c^2}+{\cos }^{-1}c\\ T_{11}={\cos }^{-1}c(1-2c)+\sqrt{1-c^2}(2-c)\\ T_{12}=\sqrt{1-c^2}(2+c)-{\cos }^{-1}c(2c+1)\\ T_{13}=\frac{1}{2}[-T_1-(c-a)T_1]\\ T_{14}=\frac{1}{16}+\frac{1}{2}ac\end{array}$$

(13)

Additionally, $C\left(k\right)$ in Equations (10)–(12) is called Theodorsen function, given as follow:

$$C(k)=\frac{H_1^{(2)}(k)}{\left[H_1^{(2)}(k)+i{H}_0^{(2)}(k)\right]}$$

(14)

where $H_v^{(2)}$ is called the Hankel function of the second kind introduced based on Bessel functions, as follows:

$$H_v^{(2)}=\mathrm{Bessel}{J}_v-I\mathrm{Bessel}{Y}_v$$

(15)

Modeling of Unstable Aerodynamic Using the Jones Approximation

Using the assumption of oscillatory motion with frequency of ω, it can be written for system coordinate:

$$\alpha (t)=\overline{\alpha }{e}^{wt},h(t)=\overline{h}{e}^{wt},\beta (t)=\overline{\beta }{e}^{wt}$$

(16)

As a result, by substituting Equation (16) in Equations (10)–(12), aerodynamic force and moment in Laplacian domain are obtained, and, by using the inverse Fourier transform, Wagner function $\varphi (s)$ and Duhamel’s integral, aerodynamic forces in time domain can be achieved with the oscillatory motion assumption for the system coordinate:

$$\begin{array}{ll}L=& -\rho b^2\left[v\pi \dot{\alpha }+\pi \ddot{h}-\pi ba\ddot{\alpha }-v{T}_4\dot{\beta }-T_{1}b\ddot{\beta }\right]-\\ & 2\pi \rho vb\{\left[v\alpha (0)+\dot{h}(0)+b(\frac{1}{2}-a)\dot{\alpha }(0)+\frac{1}{\pi }{T}_{10}v\beta (0)+b\frac{1}{2\pi }{T}_{11}\dot{\beta }(0)\right]\varphi (t)\\ & +{\displaystyle {\int }_0^t\varphi (t-\sigma )\left[v\dot{\alpha }+\ddot{h}+b(\frac{1}{2}-a)\ddot{\alpha }+\frac{1}{\pi }{T}_{10}v\dot{\beta }+b\frac{1}{2\pi }{T}_{11}\ddot{\beta }\right]}d\sigma \}\end{array}$$

(17)

$$\begin{array}{ll}{M}_{\alpha }& =-\rho b^2\left[\begin{array}{l}\pi \left(\frac{1}{2}-a\right)vb\dot{\alpha }+\pi b^2\left(\frac{1}{8}+a^2\right)\ddot{\alpha }+\left(T_4+T_{10}\right)v^2\beta +\\ \left(T_1-T_8-(c-a)T_4+\frac{1}{2}{T}_{11}\right)vb\dot{\beta }-\\ \left(T_7+(c-a)T_1\right)b^2\ddot{\beta }-a\pi b\ddot{h}\end{array}\right]+\\ & 2\rho v{b}^2\pi \left(a+\frac{1}{2}\right)\{[v\alpha (0)+\dot{h}(0)+b(\frac{1}{2}-a)\dot{\alpha }(0)+\frac{1}{\pi }{T}_{10}v\beta (0)\\ & +b\frac{1}{2\pi }{T}_{11}\dot{\beta }(0)]\varphi (t)\\ & +{\displaystyle {\int }_0^t\varphi (t-\sigma )\left[v\dot{\alpha }+\ddot{h}+b(\frac{1}{2}-a)\ddot{\alpha }+\frac{1}{\pi }{T}_{10}v\dot{\beta }+b\frac{1}{2\pi }{T}_{11}\ddot{\beta }\right]}d\sigma \}\end{array}$$

(18)

$$\begin{array}{ll}{M}_{\beta }& =-\rho b^2\left[\begin{array}{l}\left\{-2T_9-T_1+T_4(a-\frac{1}{2})\right\}vb\dot{\alpha }+2T_{13}{b}^2\ddot{\alpha }+\\ \frac{1}{\pi }{v}^2\beta \left(T_5-T_4{T}_{10}\right)-\frac{1}{2\pi }vb\dot{\beta }{T}_4{T}_{11}-\frac{1}{\pi }{T}_3{b}^2\ddot{\beta }-T_{1}b\ddot{h}\end{array}\right]\\ & -\rho v{b}^2{T}_{12}\{[v\alpha (0)+\dot{h}(0)+b(\frac{1}{2}-a)\dot{\alpha }(0)+\frac{1}{\pi }{T}_{10}v\beta (0)\\ & +b\frac{1}{2\pi }{T}_{11}\dot{\beta }(0)]\varphi (t)\\ & +{\displaystyle {\int }_0^t\varphi (t-\sigma )\left[v\dot{\alpha }+\ddot{h}+b(\frac{1}{2}-a)\ddot{\alpha }+\frac{1}{\pi }{T}_{10}v\dot{\beta }+b\frac{1}{2\pi }{T}_{11}\ddot{\beta }\right]}d\sigma \}\end{array}$$

(19)

Hence, we can obtain unstable aerodynamic equations by substituting proper Wagner function. Reference [31] represented a proper approximation for Wagner function, as follows:

$$\begin{array}{l}\varphi (t)=1-c_1{e}^{-{\epsilon }_{1}t}-c_2{e}^{-{\epsilon }_{2}t}\\ c_1=0.165,c_2=0.335\\ {\epsilon }_1=0.0455\frac{v}{b},{\epsilon }_2=0.3\frac{v}{b}\end{array}$$

(20)

The amount of downwash is calculated in ¾ of chord in order to use the Jones approximation. For two-dimensional airfoil with control surface downwash has been obtained, as follows:

$$w_{0.75c}(t)=v\alpha +\dot{h}+b(\frac{1}{2}-a)\dot{\alpha }++\frac{1}{\pi }{T}_{10}v\beta +b\frac{1}{2\pi }{T}_{11}\dot{\beta }$$

(21)

The aerodynamic relations for two-dimensional airfoil with control surface are given in Equations (22)–(24) by using the Jones approximation:

$$\begin{array}{ll}L=& -\rho b^2\left[v\pi \dot{\alpha }+\pi \ddot{h}-\pi ba\ddot{\alpha }-v{T}_4\dot{\beta }-T_{1}b\ddot{\beta }\right]-\\ & 2\pi \rho vb\{\left[v\alpha (0)+\dot{h}(0)+b(\frac{1}{2}-a)\dot{\alpha }(0)+\frac{1}{\pi }{T}_{10}v\beta (0)+b\frac{1}{2\pi }{T}_{11}\dot{\beta }(0)\right]\varphi (t)\\ & +w_{0.75c}(t)-{\displaystyle \sum _{i=1}^n{\gamma }_i{B}_i}\}\end{array}$$

(22)

$$\begin{array}{ll}{M}_{\alpha }& =-\rho b^2\left[\begin{array}{l}\pi \left(\frac{1}{2}-a\right)vb\dot{\alpha }+\pi b^2\left(\frac{1}{8}+a^2\right)\ddot{\alpha }+\left(T_4+T_{10}\right)v^2\beta +\\ \left(T_1-T_8-(c-a)T_4+\frac{1}{2}{T}_{11}\right)vb\dot{\beta }-\\ \left(T_7+(c-a)T_1\right)b^2\ddot{\beta }-a\pi b\ddot{h}\end{array}\right]+\\ & 2\rho v{b}^2\pi \left(a+\frac{1}{2}\right)\{[v\alpha (0)+\dot{h}(0)+b(\frac{1}{2}-a)\dot{\alpha }(0)+\frac{1}{\pi }{T}_{10}v\beta (0)\\ & +b\frac{1}{2\pi }{T}_{11}\dot{\beta }(0)]\varphi (t)+w_{0.75c}(t)-{\displaystyle \sum _{i=1}^n{\gamma }_i{B}_i}\}\end{array}$$

(23)

$$\begin{array}{ll}{M}_{\beta }& =-\rho b^2\left[\begin{array}{l}\left\{-2T_9-T_1+T_4(a-\frac{1}{2})\right\}vb\dot{\alpha }+2T_{13}{b}^2\ddot{\alpha }+\\ \frac{1}{\pi }{v}^2\beta \left(T_5-T_4{T}_{10}\right)-\frac{1}{2\pi }vb\dot{\beta }{T}_4{T}_{11}-\frac{1}{\pi }{T}_3{b}^2\ddot{\beta }-T_{1}b\ddot{h}\end{array}\right]\\ & -\rho v{b}^2{T}_{12}\{[v\alpha (0)+\dot{h}(0)+b(\frac{1}{2}-a)\dot{\alpha }(0)+\frac{1}{\pi }{T}_{10}v\beta (0)\\ & +b\frac{1}{2\pi }{T}_{11}\dot{\beta }(0)]\varphi (t)+w_{0.75c}(t)-{\displaystyle \sum _{i=1}^n{\gamma }_i{B}_i}\}\end{array}$$

(24)

where:

$$k=\frac{\omega b}{v},s=\frac{vt}{b},\varphi (s)=\frac{1}{2\pi i}{\displaystyle {\int }_{+\infty }^{-\infty }\frac{C(k)}{k}{e}^{iks}dk},{\gamma }_i=\frac{v}{b}{\alpha }_i{\beta }_i$$

(25)

Additionally, the below equation should be satisfied by Equations (22)–(24):

$${\dot{B}}_i+({\beta }_i\frac{v}{b})B_i=w_{0.75c}(t)$$

(26)

4. The Aero-Elastic Model

By combination of Equations (22)–(24), (5), and (8) for unstable aerodynamic and applying the Lagrange equations, the aero-elastic equations for two dimensional airfoil with control surface are obtained as:

$$\begin{array}{l}\ddot{\alpha }\left[r_{\alpha }{}^2+\kappa (\frac{1}{8}+a^2)\right]+\dot{\alpha }\left[\frac{v}{b}\kappa (\frac{1}{2}-a)\right]+\alpha \frac{C_{\alpha }}{M{b}^2}\\ +\ddot{\beta }\left[r_{\beta }{}^2+(c-a)x_{\beta }-\frac{T_7}{\pi }\kappa -(c-a)\frac{T_1}{\pi }\kappa \right]\\ +\frac{1}{\pi }\dot{\beta }\kappa \frac{v}{b}\left[-2p-\left(\frac{1}{2}-a\right)T_4\right]+\beta \kappa \frac{v^2}{b^2}\frac{1}{\pi }(T_4+T_{10})+\ddot{h}\left(x_{\alpha }-a\kappa \right)\frac{1}{b}-\\ 2\kappa \left(a+\frac{1}{2}\right)\frac{v}{b}\{\left[v\alpha (0)+\dot{h}(0)+b(\frac{1}{2}-a)\dot{\alpha }(0)+\frac{1}{\pi }{T}_{10}v\beta (0)+b\frac{1}{2\pi }{T}_{11}\dot{\beta }(0)\right]\varphi (t)\\ +w_{0.75c}(t)-{\displaystyle \sum _{i=1}^n{\gamma }_i{B}_i}\}=0\end{array}$$

(27)

$$\begin{array}{l}\ddot{\alpha }\left[r_{\beta }{}^2+(c-a)x_{\beta }-\frac{T_7}{\pi }\kappa -(c-a)\frac{T_1}{\pi }\kappa \right]+\dot{\alpha }\left(p-T_1-\frac{1}{2}{T}_4\right)\frac{v}{b}\frac{\kappa }{\pi }+\\ \ddot{\beta }\left(r_{\beta }{}^2-\frac{1}{{\pi }^2}\kappa T_3\right)-\frac{1}{2{\pi }^2}\dot{\beta }\left[T_4{T}_{11}\frac{v}{b}\kappa \right]+\beta \left[\frac{1}{{\pi }^2}\frac{v^2}{b^2}\kappa \left(T_5-T_4{T}_{10}\right)\right]+\frac{M_{S\beta }}{M{b}^2}\\ +\ddot{h}\left(x_{\beta }-\frac{1}{\pi }\kappa T_1\right)\frac{1}{b}\\ +\frac{T_{12}}{\pi }\kappa \frac{v}{b}\{\left[v\alpha (0)+\dot{h}(0)+b(\frac{1}{2}-a)\dot{\alpha }(0)+\frac{1}{\pi }{T}_{10}v\beta (0)+b\frac{1}{2\pi }{T}_{11}\dot{\beta }(0)\right]\varphi (t)\\ +w_{0.75c}(t)-{\displaystyle \sum _{i=1}^n{\gamma }_i{B}_i}\}=0\end{array}$$

(28)

$$\begin{array}{l}\ddot{\alpha }\left[x_{\alpha }-\kappa a\right]+\dot{\alpha }\frac{v}{b}\kappa +\ddot{\beta }\left(x_{\beta }-\frac{T_1}{\pi }\kappa \right)-\dot{\beta }\frac{v}{b}{T}_4\kappa \frac{1}{\pi }+\ddot{h}\left(1+\kappa \right)\frac{1}{b}+h\frac{C_h}{M}\frac{1}{b}\\ +2\kappa \frac{v}{b}\{\left[v\alpha (0)+\dot{h}(0)+b(\frac{1}{2}-a)\dot{\alpha }(0)+\frac{1}{\pi }{T}_{10}v\beta (0)+b\frac{1}{2\pi }{T}_{11}\dot{\beta }(0)\right]\varphi (t)\\ +w_{0.75c}(t)-{\displaystyle \sum _{i=1}^n{\gamma }_i{B}_i}\}=0\end{array}$$

(29)

5. Numerical Solution Result

No general analytical method has been provided yet for nonlinear problems; and, complex solving methods are being used instead.

In this paper, a new simple approach is proposed to decrease the computational efforts of the current complex approaches. This method is based on the state space integration approach. The initial conditions for airfoil in small increments of time have been calculated by an estimation technique, and these increments are added together during the solution. This procedure is repeated until the solution reaches steady-state or converges to infinite value. The accuracy of this method depends on the used estimation technique and the interval of time increments both of them are related to the type and level of the system nonlinearity.

There are two common estimation methods:

• the finite difference method; and,
• the Runge-Kutta method.

The main advantage of the Runge-Kutta algorithms is that this method eliminates the need for partial differential calculations.

In the next subsections, the solution methods and results of two-dimensional airfoil with unsteady aerodynamic with and without freeplay on freeplay is presented and discussed.

5.1. Method of Solution

For two-dimensional airfoil and with using of unsteady aerodynamics, the governing equations on aero-elastic of the system can be presented, as follows:

$$M\ddot{q}+C\dot{q}+Kq=F$$

(30)

where:

$$q=\left[\begin{array}{c}\alpha \\ \beta \\ h\end{array}\right]$$

To solve this equation using the Runge–Kutta method, it should be rewritten in state space domain. The state variables are defined, as follows:

$$\begin{array}{ccc}\begin{array}{c}{x}_1=\alpha \\ x_2=\dot{\alpha }\\ {\dot{x}}_2=\ddot{\alpha }\end{array}& \begin{array}{c}{x}_3=\beta \\ x_4=\dot{\beta }\\ {\dot{x}}_4=\ddot{\beta }\end{array}& \begin{array}{c}{x}_5=\frac{h}{b}\\ x_6=\frac{\dot{h}}{b}\\ {\dot{x}}_6=\frac{\ddot{h}}{b}\end{array}\end{array}$$

(31)

$${A\dot{Q}+BQ=F}_1$$

(32)

Accordingly, we have:

$${\dot{Q}=A}^{-1}{F}_1{-A}^{-1}BQ$$

(33)

where:

$$\begin{array}{l}A=\left[\begin{array}{ccc}\begin{array}{ccc}0& 0& 0\\ 0& 0& 0\\ 0& 0& 0\end{array}& \begin{array}{ccc}{b}_1& b_4& b_7\\ c_1& c_3& c_7\\ a_1& a_3& a_5\end{array}& \begin{array}{cc}\begin{array}{c}0\\ 0\\ 0\end{array}& \begin{array}{c}0\\ 0\\ 0\end{array}\end{array}\\ \begin{array}{ccc}1& 0& 0\\ 0& 1& 0\\ 0& 0& 1\end{array}& \begin{array}{ccc}0& 0& 0\\ 0& 0& 0\\ 0& 0& 0\end{array}& \begin{array}{cc}\begin{array}{c}0\\ 0\\ 0\end{array}& \begin{array}{c}0\\ 0\\ 0\end{array}\end{array}\\ \begin{array}{ccc}\begin{array}{c}0\\ 0\end{array}& \begin{array}{c}0\\ 0\end{array}& \begin{array}{c}0\\ 0\end{array}\end{array}& \begin{array}{ccc}\begin{array}{c}0\\ 0\end{array}& \begin{array}{c}0\\ 0\end{array}& \begin{array}{c}0\\ 0\end{array}\end{array}& \begin{array}{cc}1& 0\\ 0& 1\end{array}\end{array}\right],B=\left[\begin{array}{ccc}\begin{array}{ccc}{b}_9+b_3& b_6+b_{11}& 0\\ c_9& c_5+\underset{\_}{M_{s\beta }(\beta )}+c_{11}& 0\\ a_8& a_{10}& a_6\end{array}& \begin{array}{ccc}{b}_2+b_{10}& b_5+b_{12}& b_8\\ c_2+c_{10}& c_4+c_{12}& c_8\\ a_2+a_9& a_4+a_{11}& a_7\end{array}& \begin{array}{cc}\begin{array}{c}{b}_{13}\\ c_{13}\\ a_{12}\end{array}& \begin{array}{c}{b}_{14}\\ b_{14}\\ a_{13}\end{array}\end{array}\\ \begin{array}{ccc}0& 0& 0\\ 0& 0& 0\\ 0& 0& 0\end{array}& \begin{array}{ccc}-1& 0& 0\\ 0& -1& 0\\ 0& 0& -1\end{array}& \begin{array}{cc}\begin{array}{c}0\\ 0\\ 0\end{array}& \begin{array}{c}0\\ 0\\ 0\end{array}\end{array}\\ \begin{array}{ccc}\begin{array}{c}-v\\ -v\end{array}& \begin{array}{c}-\frac{T_{10}v}{\pi }\\ -\frac{T_{10}v}{\pi }\end{array}& \begin{array}{c}0\\ 0\end{array}\end{array}& \begin{array}{ccc}\begin{array}{c}-(\frac{1}{2}-a)b\\ -(\frac{1}{2}-a)b\end{array}& \begin{array}{c}-\frac{T_{11}b}{2\pi }\\ -\frac{T_{11}b}{2\pi }\end{array}& \begin{array}{c}-1\\ -1\end{array}\end{array}& \begin{array}{cc}\frac{{\beta }_{1}v}{b}& 0\\ 0& \frac{{\beta }_{2}v}{b}\end{array}\end{array}\right]\\ F_1=\left[\begin{array}{c}\begin{array}{c}\frac{2\kappa v}{b}(\frac{1}{2}+a)\left[v\alpha (0)+\dot{h}(0)+b(\frac{1}{2}-a)\dot{\alpha }(0)+\frac{1}{\pi }{T}_{10}v\beta (0)+b\frac{1}{2\pi }{T}_{11}\dot{\beta }(0)\right]\varphi (t)\\ \frac{-T_{12}}{\pi }\frac{\kappa v}{b}\left[v\alpha (0)+\dot{h}(0)+b(\frac{1}{2}-a)\dot{\alpha }(0)+\frac{1}{\pi }{T}_{10}v\beta (0)+b\frac{1}{2\pi }{T}_{11}\dot{\beta }(0)\right]\varphi (t)-\underset{\_}{M_{S{\beta }_0}(\beta )}\\ \frac{-2\kappa v}{b}\left[v\alpha (0)+\dot{h}(0)+b(\frac{1}{2}-a)\dot{\alpha }(0)+\frac{1}{\pi }{T}_{10}v\beta (0)+b\frac{1}{2\pi }{T}_{11}\dot{\beta }(0)\right]\varphi (t)\\ 0\end{array}\\ \begin{array}{c}0\\ 0\\ 0\\ 0\end{array}\end{array}\right],Q=\left[\begin{array}{c}\begin{array}{c}{x}_1\\ x_3\\ x_5\\ x_2\end{array}\\ \begin{array}{c}{x}_4\\ x_6\\ B_1\\ B_2\end{array}\end{array}\right]\\ \begin{array}{cc}{M}_{s\beta }(\beta )=\{\begin{array}{l}{k}_{\beta }\\ 0\\ k_{\beta }\end{array}& ,\end{array}{M}_{s{\beta }_0}(\beta )=\{\begin{array}{l}\begin{array}{cc}-k_{\beta }\delta & \beta >\delta \end{array}\\ \begin{array}{cc}0& -\delta <\beta <\delta \end{array}\\ \begin{array}{cc}{k}_{\beta }\delta & \beta <-\delta \end{array}\end{array}\end{array}$$

The values of $a_1$ to $a_{13}$, $b_1$ to $b_{14}$, and $c_1$ to $c_{14}$ are geometric parameters that their relationships are given in the appendix.

Having the initial conditions, Equation (32) or Equation (33), which are nonlinear first order ordinary equations, could be solved by the numerical integration method. The initial conditions of the system are:

$$\alpha (0),\beta (0),h(0),\dot{\alpha }(0),\dot{\beta }(0),\dot{h}(0),B_1(0),B_2(0)$$

(34)

The initial conditions $\alpha (0)=5\mathrm{deg}$, $\beta (0)=5\mathrm{deg}$, and $\dot{h}(0)=0.03$ m have been assumed in solving procedure. The initial conditions for derivatives of these parameters are achieved based on the desired speed for each mode. $B_1(0)$ and $B_2(0)$ are assumed zero value.

Moreover, the geometric parameters, inertia, stiffness, and mass of the airfoil and its control surface are shown in

Table 3. Mass, Inertial, Stiffness, and Geometrical parameters.

Parameter	Value
Mass/length of wing-aileron	$1.558\mathrm{kg}/\mathrm{m}$
Mass of wing	$0.62868\mathrm{kg}$
$S_{\alpha }$ (per span)	$0.08587\mathrm{kg}\mathrm{m}$
$S_{\beta }$ (per span)	$0.00395\mathrm{kg}\mathrm{m}$
$K_{\alpha }$	$14861/{\mathrm{s}}^2$
$K_{\beta }$	$1551/{\mathrm{s}}^2$
$K_h$	$180911/{\mathrm{s}}^2$
Span	$0.52\mathrm{m}$
Chord	$0.254\mathrm{m}$

.

Figure 3 shows the flowchart of the flutter speed analysis procedure that is used in this paper with and without freeplay. As it is shown in this figure, in the absence of the freeplay, the equations are linear and, therefore, the P-method could be used. However, adding the freeplay results in nonlinear equations that could not be dealt with the P-method anymore. In this case, time response analysis is being used to find the flutter speed for the system. The numerical results for both cases are presented here.

5.2. Numerical Results without Freeplay Effect

The results of numerical solution without control surface freeplay are described in this section. In this case, the flutter happened when one of the branches in the damping figure goes above zero. It means that the instability is started and the amplitude will be increased. It is observed that the flutter speed is equal to 23.9 m/s by plotting the time response of the system and the speed variations in generalized coordinates (Figure 4).

To confirm this conclusion, the time response of α, β, and h in generalized coordinates are shown in Figure 5, Figure 6, Figure 7, Figure 8, Figure 9, Figure 10, Figure 11, Figure 12 and Figure 13. Figure 5, Figure 6 and Figure 7 show the results for the speed of 20 m/s, as shown in these figures the results are converging for all three variables. Figure 8, Figure 9 and Figure 10 show the oscillatory time responses with constant amplitude (which occurs in the flutter speed) for the speed of 23.9 m/s. Figure 11, Figure 12 and Figure 13 confirm that the time response results of all three variables diverged for the speed of 35 m/s.

It is worthwhile to mention that the amount of amplitude of α or torsional mode of airfoil is predominated toward the bending mode or torsional mode of control surface. The results of this section confirm the validity of the used method in modeling and aero-elastic behavior prediction of the airfoil using the concept of flutter speed. Accordingly, the effects of the control surface freeplay will be discussed in the next section.

5.3. The effect of Control Surface Freeplay

The results of numerical simulation considering the freeplay effects are presented in this section. For the flutter analysis with respect to the freeplay effect, the system response time is plotted. As expected, the surface freeplay decreases the flutter speed due to decreasing the equivalent strength of the flap. The flutter speed for two-dimensional airfoil with the freeplay effect became lower with increased LCO up to 23.9 m/s.

The graphs of time responses of α, β and h have been shown below to confirm the above-mentioned conclusion. (Figure 14, Figure 15 and Figure 16 for the speed lower than the flutter, Figure 17, Figure 18 and Figure 19 for the flutter speed, and Figure 20, Figure 21 and Figure 22 for the upper flutter speed).

More results are shown in Figure 23 for different degrees of freeplay to investigate the sensitivity of the results to freeplay degrees. This figure shows that:

• the amplitude of oscillation of response in case of freeplay is a bit more than normal;
• if the degree of freeplay is more than the response of the primary domain (4 deg for the in hand problem), the time responses change dramatically and create large domain; and,
• for freeplay degrees of less than primary domain of response of control surface the time response results are close together, and all of them converged in the lower flutter speed.

5.4. Stability Analysis

The state space representation is used, as shown in Figure 24, Figure 25, Figure 26, Figure 27, Figure 28, Figure 29, Figure 30, Figure 31 and Figure 32, which presents the responses of α, β, and h for lower, upper, and in the flutter speed in order to study the stability of time responses (without freeplay and with $2^{\circ }$ freeplay). In a lower speed of flutter (Figure 24, Figure 25 and Figure 26), the responses started from the initial value, and converged in the equilibrium point. In the flutter speed (Figure 27, Figure 28 and Figure 29), the response of β has been converged, but the responses of α and h became like rings. This phenomenon occurred, because, in this case, the responses have oscillation, but the domain of oscillations has remained constant. Additionally, in upper speed of flutter (Figure 30, Figure 31 and Figure 32), the responses diverged over time and trend to infinite.

5.5. Validation of the Results

This section compares the results that were obtained from the current work with those of other published studies in order to confirm the validity of the proposed approach.

Table 4. Comparison of flutter speed for two-dimensional airfoil.

Flutter Speed (m/s)	-
23.9	Present work
23.9	Numerical [32]
20.6	Experimental [32]

shows the summary of this comparison. This table clearly confirms that the proposed approach can predict the aero-elastic behavior of complex geometry in incompressible flow with reasonable accuracy in comparison with other huge and complex algorithms and experimental methods.

6. Conclusions

This paper focused on a new low-cost method for prediction of aero-elastic behavior of complex geometry in incompressible flow with a reasonable level of accuracy. First, a two-dimensional aero-elastic airfoil model with control surface freeplay (based on Theodorsen method, Lagrange method, and Jones approximation) was developed and the generated equations were solved while using numerical integration. After checking the validity of the generated model (converging behavior in lower flutter speed, oscillatory response with constant amplitude at flutter speed, and diverging at speeds upper than flutter), the freeplay was added to the equations. In the absence of the freeplay, the P-method could be used to find the flutter speed, as the equations are linear. In the case of the freeplay effect, the time response approach is used for flutter analysis. The results showed that the flutter speed for two-dimensional airfoil without freeplay is 23.9 m/s, which is in a very good agreement with publicly available studies. In the presence of freeplay, the oscillations start at a lower speed and the amplitude of the oscillation for a particular initial condition is strongly dependent on the amount of free-play. It is also worthwhile to mention that the system response could be flutter, LCO, or chaotic, depending on the geometrical property, system inertia, and the initial velocity of the system.