Smart Wing Flutter Suppression

: In this work, it has been shown the effect of a piezoelectric material on postponing the ﬂutter phenomenon and even removing it completely on a regular wing. First, the system response of a smart wing with only plunge DOF and pitch DOF are presented. Using an efﬁcient piezopatch can effectively decay the oscillations of the smart wing in a very short time. In addition, implementing one piezopatch in the plunge DOF of a regular wing with three DOF can postpone the ﬂutter speed by 81.41%, which is a considerable increase in the ﬂutter speed. We then present the effect of adding one more piezopatch to a smart wing in the pitch DOF to further postpone the ﬂutter phenomenon. The ﬂutter speed in a smart wing can be postponed by 115.96%, which is a very considerable value. Finally, adding one more piezopatch on a smart wing in the control DOF can completely remove the ﬂutter phenomenon from the wing, which represents a great achievement in the dynamic aeroelectic behavior of a wing.


Introduction
Aeroelastic analysis of a modern wing with high flexibility is crucial. The ability to control the aeroelastic instability due to high flexibility is very important to achieve the desired aerodynamic performance in a wing design [1,2]. One important aeroelastic analysis is flutter resulting from the merging of two or more vibration modes during flight. The flutter phenomenon can reduce the flight envelope or require a redesign of the wing. Appearing flutter can compromise the long-term durability of the wing structure and the flight performance, operational safety, and energy efficiency of the aircraft. Hence, preventing flutter is crucial for the modern airplane [3][4][5][6][7].
Smart materials have been used in wing structures for many years. Although there are different smart materials, piezoelectric materials have received the most attention. Considering the direct and inverse effects of piezoelectric materials, they can perform as sensors and/or actuators on a wing. In fact, they can be used as actuators and dampers to manage the aeroelastic behaviour of the wing. One effective way to avoid redesign the wing is to use piezoelectric materials to significantly delay the flutter [8]. Piezoelectric actuators have been used in active aeroelastic control of an adaptive wing [9]. They have also been implemented in honeycomb material for a flapping wing [10]. In addition, optimal piezoelectric material has been used to control a plate subjected to time-dependent boundary moments and forcing function by vibration damping [11]. Moreover, piezoelectric patches have been studied on passive flutter control of damaged composite laminates by employing the finite element method [12]. Recently, piezoelectric layers have been used on thick porous plates to study the aeroelastic flutter [13]. Furthermore, aeroelastic optimization has been investigated for materials with piezoelectric actuators and sensors [14]. Adding a shunt circuit to a piezoelectric material can create a piezopatch to effectively modify the wing's aeroelastic behaviour. Previously, there were practical limits in the low frequency range like the one typically existing in aeroelastic phenomena due to the large required inductance in passive aeroelastic control. However, it is now possible to have a small inductor integrated into a piezopatch dedicated to aeroelastic control [15]. Since standard inductors usually have too much internal resistance for resonant shunt application, they are not a practical component to integrate into a piezopatch. Implementing closed magnetic circuits with high-permeability materials allows the design of large-inductance inductors with high quality factors.
The use of shunted piezopatch permits damping augmentation in the wing structure without causing any instability. In addition, they need little to no power and are simple to apply. Their necessary hardware includes the piezoelectrics and a simple electric circuit including a capacitor, inductor, and resistor. The shunted piezopatch can control aeroelastic vibration of the wing by consuming the energy created from wing vibrations. In fact, it can reduce the vibrations of specific modes and frequencies.
In this paper, the effect of piezoelectric material on increasing the flutter speed is investigated in detail by considering a simple aeroelastic system. The system is a 2D wing with a piezoelectric patch which has plunge, pitch, and control rotation degrees of freedom (DOF) as well as unsteady aerodynamic forces. The objective of this work is to represent the role of piezoelectric patches, which can substantially influence a simple aeropizoelastic system.
In Section 2, the equations of motion of a smart wing with only plunge DOF are represented, and we explain how to solve those equations to obtain the plunge velocity, displacement, electrical current, and electric charge. Then, the fixed points of the system and its stability around those points are investigated to show the system response. Example 1 indicates the effective decay in the oscillation of a smart wing in comparison to a regular wing. In addition, a smart wing with only pitch DOF is presented to obtain its equations of motion in free vibrations. By solving the system of equations, the angular velocity, pitch angle, electrical current, and electric charge can be obtained. The stability of the system is then investigated around the fixed point. Example 2 indicates the system response of a smart wing with only pitch DOF and how effectively the oscillation can be decayed in a smart wing in comparison to that of a regular wing.
In Section 3, a two-dimensional smart wing with the plunge, pitch, and control DOF and a piezopatch in the plunge DOF is represented to obtain the equations of motion under unsteady aerodynamic loads. Solving the system of equations produces the plunge velocity, displacement, electrical current, and electric charge as well as the pitching velocity, rotation, electrical current, and electric charge. Later, by finding the flutter speed, we show how adding one piezopatch can effectively postpone the flutter. Section 4 represents a smart wing with the plunge, pitch, and control DOF and piezopatches in the plunge and pitch DOF. We show that the flutter speed can be further increased by having two piezopatches. Finally, in Section 5, a smart wing with three DOF and three piezopatches in the plunge, pitch, and control DOF is presented to indicate how it can remove the flutter.

Aeroelastic Analysis of Smart Wing
Before studying an aeroelastic smart wing, it is necessary to study the stability of aeroelastic smart wings. The time response of the aeroelastic system is given by [16] where v i is the smart wing spatial deformation, e λ i t is the smart wing temporal deformation, and b i is the eigenvector. It is a good idea to investigate the character of the fixed point of single DOF smart wing in the plunge and pitch motions separately.

A Smart Wing with Only Plunge DOF
Consider a smart wing which has just the plunge DOF, as shown in Figure 1.

A Smart Wing with Only Plunge DOF
Consider a smart wing which has just the plunge DOF, as shown in Figure 1. The equations of motion for a smart wing with plunge DOF in free vibrations are as follows where is the mass of the smart wing, is the structural damping of the smart wing, is the structural stiffness, ℎ is the smart wing's instantaneous displacement, is the plunge electromechanical coupling, is the plunge capacitance of piezoelectric material, is the plunge inductance of piezoelectric material, is the plunge resistance of piezoelectric material, and is the plunge electric charge. As mentioned before, the plunge electromechanical coupling can be obtained as ⁄ , where is the plunge coupling coefficient. Considering ℎ , , ℎ, and , Equation (2) can be written as first-order differential equations Defining and , Equation (3) can be written as where represents linear functions, and , , , and are the smart wing states and denote the system's velocity, displacement, electrical current, and electric charge responses, respectively. The single DOF aeroelastic smart wing system has four eigenvalues that describe the stability of the fixed point. The fixed points, or static solutions, of the system are obtained from the solutions of where m is the mass of the smart wing, C h is the structural damping of the smart wing, K h is the structural stiffness, h(t) is the smart wing's instantaneous displacement, β h is the plunge electromechanical coupling, C ph is the plunge capacitance of piezoelectric material, L h is the plunge inductance of piezoelectric material, R h is the plunge resistance of piezoelectric material, and q h is the plunge electric charge. As mentioned before, the plunge electromechanical coupling can be obtained as β h = e h /C ph , where e h is the plunge h, x 2 = . q h , x 3 = h, and x 4 = q h , Equation (2) can be written as first-order differential equations Defining Equation (3) can be written as where f represents linear functions, and x 1 , x 2 , x 3 , and x 4 are the smart wing states and denote the system's velocity, displacement, electrical current, and electric charge responses, respectively. The single DOF aeroelastic smart wing system has four eigenvalues that describe the stability of the fixed point. The fixed points, or static solutions, of the system are obtained from the solutions of f(x, q) = 0 (5) or, equivalently, By considering Equation (4), Equation (6) can be written as The solution of Equation (7) can be obtained [2] x where v i is the ith eigenvector of A, λ i is the ith eigenvalue of A, and b i is the ith element of b = V −1 x 0 , where V is the eigenvector of A and x 0 is the initial condition.

Example 1.
A smart wing with plunge DOF in system response.
As the first example, a smart wing with only plunge DOF ( Figure 1) has been considered. It has the following characteristics: m = 0.3872 Kg, C h = 0.3237 Ns/m, K h = 13, 380 N/m, e h = 7.55 × 10 −3 C/m, C ph = 268 nF, L h = 106 H, R h = 4050 Ω, and the initial conditions x 1 (0) = 0 m/s, x 2 (0) = 0 A, x 3 (0) = 0.1 m, and x 4 (0) = 0 C. The system response is plotted in Figure 2. The solid line indicates the displacement of the smart wing and the dashed line depicts the displacement of the regular wing. As shown in Figure 2, the piezoelectric patch very effectively decays the vibrations. Both system responses are oscillatory, but their amplitudes decay with time towards zero. This is known as a damped response. However, the amplitude of the smart wing response decays much faster than the one of the regular wing response. The smart wing oscillation decays in almost 0.6 s, while the one of the regular wing takes around 12 s to decay. , (5) or, equivalently, By considering Equation (4), Equation (6) can be written as (7) where The solution of Equation (7) can be obtained [2] (9) where is the ith eigenvector of , is the ith eigenvalue of , and is the ith element of , where is the eigenvector of and is the initial condition.
Example 1 A smart wing with plunge DOF in system response.
As the first example, a smart wing with only plunge DOF ( Figure 1) has been considered. It has the following characteristics: and 0 0 C. The system response is plotted in Figure 2. The solid line indicates the displacement of the smart wing and the dashed line depicts the displacement of the regular wing. As shown in Figure 2, the piezoelectric patch very effectively decays the vibrations. Both system responses are oscillatory, but their amplitudes decay with time towards zero. This is known as a damped response. However, the amplitude of the smart wing response decays much faster than the one of the regular wing response. The smart wing oscillation decays in almost 0.6 s, while the one of the regular wing takes around 12 s to decay.   Furthermore, the phase plane plot for the velocity and displacement indicates that the point (0, 0) evokes the system trajectory, as shown in Figure 3. The smart wing trajectory starts from the initial displacement and velocity at the far right and turns to the centre of the phase plane where (0, 0) is the fixed point, x F = 0. In fact, the phase plane plot shows that the fixed point attracts the smart wing trajectory. Furthermore, the phase plane plot for the velocity and displacement indicates tha the point 0,0 evokes the system trajectory, as shown in Figure 3. The smart wing tra jectory starts from the initial displacement and velocity at the far right and turns to th centre of the phase plane where 0,0 is the fixed point, . In fact, the phase plane plot shows that the fixed point attracts the smart wing trajectory. Likewise, the phase plane for the electrical current and charge begins at the initia conditions for the electric charge and current, which are zeros, and it twists out counter clockwise until reaching the maximum values. The trajectory then turns towards the star point 0,0 , as depicted in Figure 4.

A Smart Wing with only Pitch DOF
As the second investigation for the fixed point, we considered a smart wing with only pitch DOF, as depicted in Figure 5.  Likewise, the phase plane for the electrical current and charge begins at the initial conditions for the electric charge and current, which are zeros, and it twists out counterclockwise until reaching the maximum values. The trajectory then turns towards the start point (0, 0), as depicted in Figure 4.
Furthermore, the phase plane plot for the velocity and displacement indicates that the point 0,0 evokes the system trajectory, as shown in Figure 3. The smart wing trajectory starts from the initial displacement and velocity at the far right and turns to the centre of the phase plane where 0,0 is the fixed point, . In fact, the phase plane plot shows that the fixed point attracts the smart wing trajectory. Likewise, the phase plane for the electrical current and charge begins at the initial conditions for the electric charge and current, which are zeros, and it twists out counterclockwise until reaching the maximum values. The trajectory then turns towards the start point 0,0 , as depicted in Figure 4.

A Smart Wing with only Pitch DOF
As the second investigation for the fixed point, we considered a smart wing with only pitch DOF, as depicted in Figure 5.

A Smart Wing with Only Pitch DOF
As the second investigation for the fixed point, we considered a smart wing with only pitch DOF, as depicted in Figure 5. The smart wing equations of motion in free vibrations are given as follows: where is the mass moment of inertia, is the pitching acceleration, is the pitching structural damping, is the pitching velocity, is the pitching structural stiffness, is the pitching angle, is the pitch electromechanical coupling, is the pitch electric charge, is the pitch inductance of piezoelectric material, is the rate of the pitch electrical current, is the pitch resistance of piezoelectric material, is the pitch electrical current, is the pitch capacitance of piezoelectric material, is the elastic axis, and is the piezoelectric axis, as depicted in Figure 5. Equation (10) can be written as first-order differential equations by assuming , , , and Equation (11) can be rewritten by considering and , as follows: , where are linear functions and , , , and are the smart wing states and denote the system's pitching velocity, pitching angle, pitch electrical current, and pitch electric charge responses, respectively. There are four eigenvalues for the single DOF aeroelastic system that indicate the stability of the fixed point. One can obtain the fixed point, or static solution, of the system by solving , (13) or, equivalently, The smart wing equations of motion in free vibrations are given as follows: where I α is the mass moment of inertia, ..
α is the pitching acceleration, C α is the pitching structural damping, . α is the pitching velocity, K α is the pitching structural stiffness, α is the pitching angle, β α is the pitch electromechanical coupling, q α is the pitch electric charge, L α is the pitch inductance of piezoelectric material, .. q α is the rate of the pitch electrical current, R α is the pitch resistance of piezoelectric material, . q α is the pitch electrical current, C pα is the pitch capacitance of piezoelectric material, x f is the elastic axis, and x p is the piezoelectric axis, as depicted in Figure 5. Equation (10) can be written as first-order differential equations by assuming Equation (11) can be rewritten by considering q = [I α C α K α β α L α C pα R α ] T and T , as follows: where f are linear functions and x 1 , x 2 , x 3 , and x 4 . are the smart wing states and denote the system's pitching velocity, pitching angle, pitch electrical current, and pitch electric charge responses, respectively. There are four eigenvalues for the single DOF aeroelastic system that indicate the stability of the fixed point. One can obtain the fixed point, or static solution, of the system by solving f(x, q) = 0 (13) or, equivalently, Using Equation (12), it is possible to write Equation (14) as Designs 2022, 6, 29 The solution of Equation (15) can be given as [9] x where where V is eigenvector of A and x 0 is the initial condition.

Example 2.
A smart wing with pitch DOF in system response.
In the second example, a smart wing with only pitch DOF has been given, as shown in Figure 6. The smart wing characteristics are m = 0.3872 Kg, C α = 0.1 Nms/rad, 3c, x f = 0.4c, and c = 0.25 m. In addition, the initial conditions are given as x 1 (0) = 0 rad/s, x 2 (0) = 0 A, x 3 (0) = 0.1 rad, and x 4 (0) = 0 C. The mass moment of inertia of the smart wing by assuming the uniform thickness and mass distribution can be calculated as The system response of the smart wing is depicted in Figure 6. The solid line represents the pitching angle of the smart wing and the dashed line shows the pitching angle of the corresponding regular wing. It is clear that piezoelectric patch decays effectively the pitching vibrations, as shown in Figure 6. The oscillation of the smart wing decays almost 0.0198 s however, the vibration of the corresponding regular wing takes 0.3 s to decay.
Using Equation (12), it is possible to write Equation (14) as (15) where The solution of Equation (15) can be given as [9] (17) where is the ith eigenvector of , is the ith eigenvalue of , and is the ith element of , where is eigenvector of and is the initial condition.
Example 2 A smart wing with pitch DOF in system response.
In the second example, a smart wing with only pitch DOF has been given, as shown in Figure   In addition, the phase plane plot for the pitching velocity and angle has been shown in Figure 7, where the point 0,0 evokes the system trajectory. The initial pitching angle and velocity is the start point of the smart wing trajectory at the far right and the trajectory is turning to the fixed point, , where the center of the phase plane is 0,0 . In addition, the phase plane plot for the pitching velocity and angle has been shown in Figure 7, where the point (0, 0) evokes the system trajectory. The initial pitching angle and velocity is the start point of the smart wing trajectory at the far right and the trajectory is turning to the fixed point, x F = 0, where the center of the phase plane is (0, 0).  The initial conditions for the electric charge and current, which are zeros, is the start point for the phase plane for the electrical current and charge. The trajectory twists out counterclockwise until reaching its maximum values; after that it turns towards the start point 0,0 , as indicated in Figure 8. The effects of different values on the system response have been shown in Figure  9. Low values of mean the piezoelectric patch located very close to the leading edge can decay the pitching oscillation quicker than the high-values one.  The initial conditions for the electric charge and current, which are zeros, is the start point for the phase plane for the electrical current and charge. The trajectory twists out counterclockwise until reaching its maximum values; after that it turns towards the start point (0, 0), as indicated in Figure 8.  The initial conditions for the electric charge and current, which are zeros, is the start point for the phase plane for the electrical current and charge. The trajectory twists out counterclockwise until reaching its maximum values; after that it turns towards the start point 0,0 , as indicated in Figure 8. The effects of different values on the system response have been shown in Figure  9. Low values of mean the piezoelectric patch located very close to the leading edge can decay the pitching oscillation quicker than the high-values one.  The effects of different x p values on the system response have been shown in Figure 9. Low values of x p mean the piezoelectric patch located very close to the leading edge can decay the pitching oscillation quicker than the high-values one.
In Section 3, we illustrate how a smart wing with three DOF, plunge (including a piezopatch), pitch, and control, called a piezo-plunge-wing, can behave in comparison to a regular wing in aeroelastic analysis, and how adding one piezopatch can effectively postpone the flutter phenomenon on a smart wing. The effects of different values on the system response have been shown in Figure  9. Low values of mean the piezoelectric patch located very close to the leading edge can decay the pitching oscillation quicker than the high-values one. In Section 3, we illustrate how a smart wing with three DOF, plunge (including a piezopatch), pitch, and control, called a piezo-plunge-wing, can behave in comparison to a regular wing in aeroelastic analysis, and how adding one piezopatch can effectively postpone the flutter phenomenon on a smart wing.

Two-Dimensional Unsteady Plunge-Pitch-Control Smart Wing
A linear 2D smart wing which has plunge, pitch, and control degrees of freedom is shown in Figure 10. The model contains an airfoil with a piezoelectric patch in the plunge DOF, called piezo-plunge-wing. The system has the plunge, pitch, and control DOF indi-

Two-Dimensional Unsteady Plunge-Pitch-Control Smart Wing
A linear 2D smart wing which has plunge, pitch, and control degrees of freedom is shown in Figure 10. The model contains an airfoil with a piezoelectric patch in the plunge DOF, called piezo-plunge-wing. The system has the plunge, pitch, and control DOF indicated by h, α, and β, respectively. The DOF β represents the angle of the control surface around its hinge, located at distance x h from the leading edge, and K β denotes the stiffness of the control surface. The equations of motion can be obtained by using the Lagrange's equations and the Kirchhoff's law as [16,17] where is the mass per unit length of the wing, is the static mass moment of the wing around the pitch axis , is the mass moment of inertia around the pitch axis , is the static mass moment of the control surface around the hinge axis , is the control surface moment of inertia around the hinge axis, is the product of inertia of the wing and control surface, is the lift, is the pitching moment of the wing around the pitch axis , is the pitching moment of the control surface around the hinge axis , is the plunge electromechanical coupling, is the plunge capacitance of piezoelectric material, is the plunge inductance of piezoelectric material, is the plunge resistance of piezoelectric material, and is the plunge electric charge. The electromechanical coupling, , depends on the plunge coupling coefficient, , and the plunge capacitance, , and it can be calculated by ⁄ . Considering unsteady aerodynamics, the lift and moments can be written as follows [18,19]  The equations of motion can be obtained by using the Lagrange's equations and the Kirchhoff's law as [16,17] where m is the mass per unit length of the wing, S αh is the static mass moment of the wing around the pitch axis x f , I α is the mass moment of inertia around the pitch axis x f , S β is the static mass moment of the control surface around the hinge axis x h , I β is the control surface moment of inertia around the hinge axis, I αβ is the product of inertia of the wing and control surface, L is the lift, M x f is the pitching moment of the wing around the pitch axis x f , M xh is the pitching moment of the control surface around the hinge axis x h , β h is the plunge electromechanical coupling, C ph is the plunge capacitance of piezoelectric material, L h is the plunge inductance of piezoelectric material, R h is the plunge resistance of piezoelectric material, and q h is the plunge electric charge. The electromechanical coupling, β h , depends on the plunge coupling coefficient, e h , and the plunge capacitance, C ph , and it can be calculated by β h = e h /C ph . Considering unsteady aerodynamics, the lift and moments can be written as follows [18,19]  Substituting Equations (19) to (21) into Equation (18) gives a set of equations of motion which is only time dependent and can be solved numerically by using the backward finite difference scheme for numerical integration [19]. However, implementing the exponential form of the Wagner function's approximation, the equations of motion can be rewritten as ordinary differential equations which can be solved analytically rather than numerically, which is much more practical [20,21]. The Wagner function's approximation is as follows The full unsteady aeroelastic equations of motion can be written as ..
is the aerodynamic states vector, Φ(t) is Wagner's function, A is the structural mass and inductance matrix, B is the aerodynamic mass matrix, E is the structural stiffness and resistance matrix, F is the aerodynamic stiffness matrix, W is the aerodynamic state influence matrix, g is the initial condition excitation vector, and W 1 and W 2 are the aerodynamic state equation matrices. Equation (23) can be written in purely first order ordinary differential equations form by where T is the 14 × 1 state vector, M = A + ρB, I 4×4 is a 4 × 4 unit matrix, 0 4×4 is a 4 × 4 matrix of zeros, 0 4×6 is a 4 × 6 matrix of zeros, 0 6×4 is a 6 × 4 matrix of zeros, and 0 10×1 is a 10 × 1 vector of zeros. The initial conditions are x(0) = x 0 . The initial condition g . Φ(t) is an excitation whose effect decays exponentially. In order to obtain steady-state solutions, the initial condition is ignored in this paper; therefore, Equation (24) becomes Example 3. A smart wing with plunge, pitch, and control DOF and a piezopatch in plunge DOF.
After running the simulation, the flutter speed was 74.2973 m/s, which showed an 81.41% increase in the flutter speed of a regular wing with the same characteristics and without the piezoelectric patch. Figure 11 depicts the variation of damping ratios of a regular wing and smart wing with respect to the airflow velocity or airspeed. It can be seen that having a piezoelectric patch on the wing can effectively increase the flutter speed.
Furthermore, due to the piezoelectric effect, there was no flutter in the plunge mode; however, flutter happens in the pitch mode, as shown in Figure 11.
Designs 2022, 6, x FOR PEER REVIEW 12 of 25 After running the simulation, the flutter speed was 74.2973 m s ⁄ , which showed an 81.41% increase in the flutter speed of a regular wing with the same characteristics and without the piezoelectric patch. Figure 11 depicts the variation of damping ratios of a regular wing and smart wing with respect to the airflow velocity or airspeed. It can be seen that having a piezoelectric patch on the wing can effectively increase the flutter speed. Furthermore, due to the piezoelectric effect, there was no flutter in the plunge mode; however, flutter happens in the pitch mode, as shown in Figure 11. In addition, the real part of eigenvalues versus the freestream velocity is shown in Figure 12. Again, Figure 12b indicates flutter appears in the pitch mode. There is an effective increase in the flutter speed of the smart wing in comparison to the regular wing. Furthermore, the imaginary part of eigenvalues versus the freestream velocity is depicted in Figure 13. Figure 13b shows that flutter appeared in the pitch mode and there was an effective increase in the flutter speed of the smart wing in comparison to the regular wing. In addition, the real part of eigenvalues versus the freestream velocity is shown in Figure 12. Again, Figure 12b indicates flutter appears in the pitch mode. There is an effective increase in the flutter speed of the smart wing in comparison to the regular wing. After running the simulation, the flutter speed was 74.2973 m s ⁄ , which showed an 81.41% increase in the flutter speed of a regular wing with the same characteristics and without the piezoelectric patch. Figure 11 depicts the variation of damping ratios of a regular wing and smart wing with respect to the airflow velocity or airspeed. It can be seen that having a piezoelectric patch on the wing can effectively increase the flutter speed. Furthermore, due to the piezoelectric effect, there was no flutter in the plunge mode; however, flutter happens in the pitch mode, as shown in Figure 11. In addition, the real part of eigenvalues versus the freestream velocity is shown in Figure 12. Again, Figure 12b indicates flutter appears in the pitch mode. There is an effective increase in the flutter speed of the smart wing in comparison to the regular wing. Furthermore, the imaginary part of eigenvalues versus the freestream velocity is depicted in Figure 13. Figure 13b shows that flutter appeared in the pitch mode and there was an effective increase in the flutter speed of the smart wing in comparison to the regular wing. Furthermore, the imaginary part of eigenvalues versus the freestream velocity is depicted in Figure 13. Figure 13b shows that flutter appeared in the pitch mode and there was an effective increase in the flutter speed of the smart wing in comparison to the regular wing. By using Equation (8), the matrix Q can be formed and its eigenvalues and eigenvectors can be calculated for two different airspeeds, U = 10 m/s and the flutter speed, U = 74.2973 m/s. There are six complex eigenvalues which represent the structural state dynamics of the smart wing. These complex eigenvalues are conjugates of those of the regular wing. There are six real eigenvalues for the aerodynamics state dynamics. Furthermore, there are two real eigenvalues representing the piezoelectric state dynamics. In each eigenvector, the first three elements give structural velocities, the next three correspond to structural displacements, the next six elements represent aerodynamic state displacements, and finally, the last two are for piezoelectric electric charge.
The smart wing eigenvalues for the three structural modes at U = 10 m/s are as follows λ 1 = −1.3460 ± 42.7410i, λ 2 = −6.2904 ± 110.9803i, λ 3 = −5.4720 ± 205.9954i and its corresponding eigenvectors which present the smart wing structural mode shapes are where, in each mode shape, the first element presents plunge displacement, the second one indicates pitch angle, and the last one gives control surface angle. Generally, in aeroelastic systems, the degrees of freedom are coupled to each other and cannot occur independently. Mostly, control surface and pitch displacements happen in mode two and three. Despite the regular wing, mode one contains significant pitch angle. The smart wing deformation in the three modes has been depicted in Figure 14. There is almost similarity in pitch and control in modes two and three; however, there is significant pitch in mode one. In addition, the smart wing eigenvalues at airspeed U = 74.2973 m/s are λ 1 = −21.2035 ± 13.2734i, λ 2 = 0.0000 ± 97.5068i, λ 3 = −6.0159 ± 204.6810i and its corresponding mode shapes are In comparison to eigenvalues at airspeed U = 10 m/s, the real parts of λ 1 are much more negative, and the real part of λ 2 is almost zero. In addition, at U = 74.2973 m/s, the control components of mode shapes ϕ 2 and ϕ 3 are very close together and there is a significant pitch in mode one, as shown in Figure 15. In Section 4, a smart wing with three DOF and two piezopatches in the plunge and pitch DOF, called piezo-plunge-pitch-wing, is presented to illustrate its behavior in comparison to a regular wing in aeroelastic analysis and how implementing two piezopatches can further postpone the flutter phenomenon on a smart wing.

A Smart Wing with Plunge, Pitch, and Control DOF and Piezopatches in Plunge and Pitch DOF
A 2D smart wing with plunge, pitch, and control DOF which has two piezopatches, one in the plunge DOF and the other in the pitch DOF, is considered, as indicated in Figure  16; it can be called a piezo-plunge-pitch-wing. The system has the same characteristics of the piezo-plunge-wing in section three. In Section 4, a smart wing with three DOF and two piezopatches in the plunge and pitch DOF, called piezo-plunge-pitch-wing, is presented to illustrate its behavior in comparison to a regular wing in aeroelastic analysis and how implementing two piezopatches can further postpone the flutter phenomenon on a smart wing.

A Smart Wing with Plunge, Pitch, and Control DOF and Piezopatches in Plunge and Pitch DOF
A 2D smart wing with plunge, pitch, and control DOF which has two piezopatches, one in the plunge DOF and the other in the pitch DOF, is considered, as indicated in Figure 16; it can be called a piezo-plunge-pitch-wing. The system has the same characteristics of the piezo-plunge-wing in Section 3.

A Smart Wing with Plunge, Pitch, and Control DOF and Piezopatches in Plunge and Pitch DOF
A 2D smart wing with plunge, pitch, and control DOF which has two piezopatches, one in the plunge DOF and the other in the pitch DOF, is considered, as indicated in Figure  16; it can be called a piezo-plunge-pitch-wing. The system has the same characteristics of the piezo-plunge-wing in section three.  The equations of motion of the smart wing can be obtained by using the Lagrange's equations and the Kirchhoff's law as , and x p are defined as in Equation (1), L α is the pitch inductance of piezoelectric material, R α is the pitch resistance of piezoelectric material, C pα is the pitch capacitance of piezoelectric material, β α is the pitch electromechanical coupling, and q α is the pitch electric charge. The pitch electromechanical coupling, β α , depends on the pitch coupling coefficient, e α , and the pitch capacitance, C pα , and it can be calculated by β α = e α /C pα . Similar to Section 3, the full unsteady aeroelastic equations of motion can be written as ..
where y = h α β q h q α T is the displacement and charge vector.
Equation (29) can be written in purely first order ordinary differential equations form by q α h α β q h q α w 1 · · · w 6 T is the 16 × 1 state vector, M = A + ρB, I 5×5 is a 5 × 5 unit matrix, 0 5×5 is a 5 × 5 matrix of zeros, 0 5×6 is a 5 × 6 matrix of zeros, 0 6×5 is a 6 × 5 matrix of zeros, and 0 11×1 is a 11 × 1 vector of zeros. The initial conditions are x(0) = x 0 . The initial condition g . Φ(t) is an excitation whose effect decays exponentially. In order to obtain steady-state solutions, the initial condition is ignored in this paper therefore Equation (30) becomes Example 4. A smart wing with plunge, pitch, and control DOF and a piezopatch in plunge and pitch DOF in aeroelastic analysis.
In the fourth example, an additional piezopatch was used to control vibrations in pitch DOF in the smart wing in example three. Therefore, a smart wing with plunge, pitch, and control DOF and two piezopatches, one in the plunge and the other in the pitch DOF, as shown in Figure 16, was considered with the following parameters. It assumed the same characteristics as the smart wing existing in example three with the pitch piezopatch parameters as the pitch coupling coefficient e α = 0.019 C/m, the pitch capacitance of piezoelectric material C pα = 1450 nF, the pitch inductance of piezoelectric material L α = 103 H, and the pitch resistance of piezoelectric material R α = 2674 Ω.
Simulation results indicated that adding an additional piezopatch in the pitch DOF can remove the flutter phenomenon in the pitch mode, as depicted in Figure 17. Hence, by having two piezopatches, one in the plunge and the other in the pitch DOF, it is possible to avoid flutter in both DOF. However, there was flutter in the control DOF.  Figure 17 shows that the smart wing with two piezopatches has flutter at 88.4353 m s ⁄ in the control DOF, which indicates a 115.96% increase in the flutter speed of a regular wing with the same characteristics without the piezopatch and a 19.03% increase in the flutter speed of a smart wing with the same characteristics that has only one piezopatch in the plunge DOF. Clearly, using two piezopatches can remove the flutter phenomenon in the plunge and pitch modes, but flutter happens in the control mode, as shown in Figure 17b.
Furthermore, the real part of eigenvalues versus the freestream velocity has been depicted in Figure 18. In Figure 18b, it can be seen that there is no flutter in the plunge and pitch modes; however, flutter appears in the control mode. There is also a considerable increase in the flutter speed of the smart wing with two piezopatches in comparison to that of the smart wing with only one piezopatch.  Figure 17 shows that the smart wing with two piezopatches has flutter at 88.4353 m/s in the control DOF, which indicates a 115.96% increase in the flutter speed of a regular wing with the same characteristics without the piezopatch and a 19.03% increase in the flutter speed of a smart wing with the same characteristics that has only one piezopatch in the plunge DOF. Clearly, using two piezopatches can remove the flutter phenomenon in the plunge and pitch modes, but flutter happens in the control mode, as shown in Figure 17b.
Furthermore, the real part of eigenvalues versus the freestream velocity has been depicted in Figure 18. In Figure 18b, it can be seen that there is no flutter in the plunge and pitch modes; however, flutter appears in the control mode. There is also a considerable increase in the flutter speed of the smart wing with two piezopatches in comparison to that of the smart wing with only one piezopatch. phenomenon in the plunge and pitch modes, but flutter happens in the control mode, as shown in Figure 17b.
Furthermore, the real part of eigenvalues versus the freestream velocity has been depicted in Figure 18. In Figure 18b, it can be seen that there is no flutter in the plunge and pitch modes; however, flutter appears in the control mode. There is also a considerable increase in the flutter speed of the smart wing with two piezopatches in comparison to that of the smart wing with only one piezopatch. In addition, the imaginary part of eigenvalues versus the freestream velocity is indicated in Figure 19. Figure 19b shows that flutter appears in the control mode and there is an effective increase in the flutter speed of the smart wing in comparison to the regular wing. In addition, the imaginary part of eigenvalues versus the freestream velocity is indicated in Figure 19. Figure 19b shows that flutter appears in the control mode and there is an effective increase in the flutter speed of the smart wing in comparison to the regular wing. By using Equation (31), the matrix can be formed, and its eigenvalues and eigenvectors can be calculated for two different airspeeds, = 10 m s ⁄ and the flutter speed, = 88.4353 m s ⁄ . There are six complex eigenvalues which represent the structural state dynamics of the smart wing. These complex eigenvalues are conjugates of those of the regular wing. There are six real eigenvalues for the aerodynamics state dynamics. Furthermore, there are four real eigenvalues representing the piezoelectric state dynamics. In each eigenvector, the first three elements give structural velocities, the next three correspond to structural displacements, the next six elements represent aerodynamic state displacements, and finally, the last four are for piezoelectric electric charges.
The smart wing eigenvalues for the three structural modes at = 10 m s ⁄ are as follows where, in each mode shape, the first element presents plunge displacement, the second indicates pitch angle, and the last element provides control surface angle. Generally, in aeroelastic systems, the degrees of freedom are coupled to each other and cannot occur independently. Mostly, control surface and pitch displacements happen in mode two and three. In contrast to the regular wing, mode one contains significant pitch angle. The By using Equation (31), the matrix Q can be formed, and its eigenvalues and eigenvectors can be calculated for two different airspeeds, U = 10 m/s and the flutter speed, U = 88.4353 m/s. There are six complex eigenvalues which represent the structural state dynamics of the smart wing. These complex eigenvalues are conjugates of those of the regular wing. There are six real eigenvalues for the aerodynamics state dynamics. Furthermore, there are four real eigenvalues representing the piezoelectric state dynamics. In each eigenvector, the first three elements give structural velocities, the next three correspond to structural displacements, the next six elements represent aerodynamic state displacements, and finally, the last four are for piezoelectric electric charges.
The smart wing eigenvalues for the three structural modes at U = 10 m/s are as follows where, in each mode shape, the first element presents plunge displacement, the second indicates pitch angle, and the last element provides control surface angle. Generally, in aeroelastic systems, the degrees of freedom are coupled to each other and cannot occur independently. Mostly, control surface and pitch displacements happen in mode two and three. In contrast to the regular wing, mode one contains significant pitch angle. The piezo-plunge-pitch-wing deformation in the three modes has been depicted in Figure 20. It is clear that there is slight similarity in pitch and control in mode two and three; however, there is significant pitch in mode one.
The smart wing eigenvalues for the three structural modes at = 10 m s ⁄ are as follows where, in each mode shape, the first element presents plunge displacement, the second indicates pitch angle, and the last element provides control surface angle. Generally, in aeroelastic systems, the degrees of freedom are coupled to each other and cannot occur independently. Mostly, control surface and pitch displacements happen in mode two and three. In contrast to the regular wing, mode one contains significant pitch angle. The piezo-plunge-pitch-wing deformation in the three modes has been depicted in Figure 20. It is clear that there is slight similarity in pitch and control in mode two and three; however, there is significant pitch in mode one. In addition, the smart wing eigenvalues at airspeed U = 88.4353 m/s are λ 1 = −39.2529 ± 11.3444i, λ 2 = −11.7376 ± 80.3532i, λ 3 = 0.0000 ± 96.3908i and its corresponding mode shapes are In comparison to eigenvalues at airspeed U = 10 m/s, the real parts of λ 1 and λ 2 are much closer, the real part of λ 1 is much more negative, and the real part of λ 3 is almost zero. In addition, at U = 88.4353 m/s, the control components of mode shapes ϕ 2 and ϕ 3 are almost symmetrical, and there is a significant pitch in mode one, as shown in Figure 21.
Next section presents a smart wing with three DOF and three piezopatches in the plunge, pitch, and control DOF, called a piezo-plunge-pitch-control-wing, to illustrate its behaviour in comparison to a regular wing in aeroelastic analysis, and how using three piezopatches can completely remove the flutter phenomenon on a smart wing. In comparison to eigenvalues at airspeed = 10 m s ⁄ , the real parts of 1 and 2 are much closer, the real part of 1 is much more negative, and the real part of 3 is almost zero. In addition, at = 88.4353 m s ⁄ , the control components of mode shapes 2 and 3 are almost symmetrical, and there is a significant pitch in mode one, as shown in Figure 21. Next section presents a smart wing with three DOF and three piezopatches in the plunge, pitch, and control DOF, called a piezo-plunge-pitch-control-wing, to illustrate its

A Smart Wing with Plunge, Pitch, and Control DOF and Piezopatches in the Plunge, Pitch, and Control DOF
A 2D smart wing with plunge, pitch, and control DOF which has three piezopatches, one in the plunge DOF, one in the pitch DOF, and the other one in the control DOF, is considered as indicated in Figure 22; it can be called a piezo-plunge-pitch-control-wing. The system has the same characteristics of the smart wing in the previous section. behaviour in comparison to a regular wing in aeroelastic analysis, and how using three piezopatches can completely remove the flutter phenomenon on a smart wing.

A Smart Wing with Plunge, Pitch, and Control DOF and Piezopatches in the Plunge, Pitch, and Control DOF.
A 2D smart wing with plunge, pitch, and control DOF which has three piezopatches, one in the plunge DOF, one in the pitch DOF, and the other one in the control DOF, is considered as indicated in Figure 22; it can be called a piezo-plunge-pitch-control-wing. The system has the same characteristics of the smart wing in the previous section. The equations of motion of the smart wing can be obtained by using the Lagrange's equations and the Kirchhoff's law as x p , L α , R α , C pα , and β α are defined as in Equation (28), L β is the control inductance of piezoelectric material, R β is the control resistance of piezoelectric material, C pβ is the control capacitance of piezoelectric material, β β is the control electromechanical coupling, and q β is the control electric charge. The control electromechanical coupling, β β , depends on the control coupling coefficient, e β , and the control capacitance, C pβ , and it can be calculated by β β = e β /C pβ . Similar to Section 3, the full unsteady aeroelastic equations of motion can be written as ..
where y = h α β q h q α q β T is the displacement and charge vector.
Equation (35) can be written in purely first order ordinary differential equations form by where q β h α β q h q α q β w 1 · · · w 6 T is the 18 × 1 state vector, M = A + ρB, I 6×6 is a 6 × 6 unit matrix, 0 6×6 is a 6 × 6 matrix of zeros, 0 6×6 is a 6 × 6 matrix of zeros, and 0 12×1 is a 12 × 1 vector of zeros. The initial conditions are x(0) = x 0 . The initial condition g . Φ(t) is an excitation whose effect decays exponentially. In order to obtain steady-state solutions, the initial condition is ignored in this paper; therefore, Equation (36) becomes Example 5. A smart wing with plunge, pitch, and control DOF and a piezopatch in plunge, pitch, and control DOF in aeroelastic analysis.
In the fifth example, an additional piezopatch was used to control vibrations in control DOF in the smart wing in example four. Therefore, a smart wing with plunge, pitch, and control DOF and three piezopatches, one in plunge, one in pitch, and the other in control DOF, as shown in Figure 22, was considered with the following parameters. It assumed the same characteristics for smart wing existing in example four with the control piezopatch parameters as the control coupling coefficient e β = 0.019 C/m, the control capacitance of piezoelectric material C pβ = 26.8 nF, the control inductance of piezoelectric material L β = 103 H, and the control resistance of piezoelectric material R β = 1274 Ω.
Simulation results indicated that adding an additional piezopatch in the control DOF can remove the flutter phenomenon in the control mode, as depicted in Figure 23. Hence, by having three piezopatches, one in the plunge, one in the pitch, and the other in the control DOF, it is possible to avoid flutter in all DOF, eradicating all flutter.
control DOF in the smart wing in example four. Therefore, a smart wing with plunge, pitch, and control DOF and three piezopatches, one in plunge, one in pitch, and the other in control DOF, as shown in Figure 22, was considered with the following parameters. It assumed the same characteristics for smart wing existing in example four with the control piezopatch parameters as the control coupling coefficient = 0.019 C m ⁄ , the control capacitance of piezoelectric material = 26.8 nF , the control inductance of piezoelectric material = 103 H, and the control resistance of piezoelectric material = 1274 Ω.
Simulation results indicated that adding an additional piezopatch in the control DOF can remove the flutter phenomenon in the control mode, as depicted in Figure 23. Hence, by having three piezopatches, one in the plunge, one in the pitch, and the other in the control DOF, it is possible to avoid flutter in all DOF, eradicating all flutter.  Figure 23b shows that the smart wing with three piezopatches has no flutter. Using three piezopatches can remove the flutter phenomenon completely in the plunge, pitch, and control modes, as shown in Figure 23b.
Furthermore, the real part of the eigenvalues versus the freestream velocity has been depicted in Figure 24. In Figure 24b, it is clear that there is no flutter in the plunge, pitch, or control modes.  Figure 23b shows that the smart wing with three piezopatches has no flutter. Using three piezopatches can remove the flutter phenomenon completely in the plunge, pitch, and control modes, as shown in Figure 23b.
Furthermore, the real part of the eigenvalues versus the freestream velocity has been depicted in Figure 24. In Figure 24b, it is clear that there is no flutter in the plunge, pitch, or control modes. Moreover, the imaginary part of eigenvalues versus the freestream velocity has been shown in Figure 25.  Moreover, the imaginary part of eigenvalues versus the freestream velocity has been shown in Figure 25. Figure 25b indicates flutter does not appear in the smart wing.
By using Equation (37), the matrix Q can be formed and its eigenvalues and eigenvectors can be calculated for two different airspeeds, U = 10 m/s and the flutter speed, U = 100 m/s. There are six complex eigenvalues which represent the structural state dynamics of the smart wing. These complex eigenvalues are conjugates of those of the regular wing. There are six real eigenvalues for the aerodynamics state dynamics. Furthermore, there are six real eigenvalues representing the piezoelectric state dynamics. In each eigenvector, the first three elements give structural velocities, the next three correspond to structural displacements, the next six elements represent aerodynamic state displacements, and finally, the last six are for piezoelectric electric charges. Moreover, the imaginary part of eigenvalues versus the freestream velocity has been shown in Figure 25. Figure 25b indicates flutter does not appear in the smart wing. By using Equation (37), the matrix can be formed and its eigenvalues and eigenvectors can be calculated for two different airspeeds, = 10 m s ⁄ and the flutter speed, = 100 m s ⁄ . There are six complex eigenvalues which represent the structural state dynamics of the smart wing. These complex eigenvalues are conjugates of those of the regular wing. There are six real eigenvalues for the aerodynamics state dynamics. Furthermore, there are six real eigenvalues representing the piezoelectric state dynamics. In each eigenvector, the first three elements give structural velocities, the next three correspond to structural displacements, the next six elements represent aerodynamic state displacements, and finally, the last six are for piezoelectric electric charges.
The smart wing eigenvalues for the three structural modes at = 10 m s ⁄ are as follows The smart wing eigenvalues for the three structural modes at U = 10 m/s are as follows λ 1 = −2.3739 ± 46.5181i, λ 2 = −12.1637 ± 82.6875i, λ 3 = −5.4879 ± 204.1061i and its corresponding eigenvectors which present the smart wing structural mode shapes are where, in each mode shape, the first element presents plunge displacement, the second indicates pitch angle, and the last element gives control surface angle. Generally, in aeroelastic systems, the degrees of freedom are coupled to each other and cannot occur independently. Mostly, control surface and pitch displacements happen in mode two and three. In contrast to the regular wing, mode one contains significant pitch angle. The smart wing deformation in the three modes has been depicted in Figure 26. It is clear that there is almost symmetry in pitch in modes one and two; however, their values are significant.
In addition, the smart wing eigenvalues at airspeed U = 100 m/s are λ 1 = −40.2408 ± 12.0910i, λ 2 = −12.8883 ± 81.8056i, λ 3 = −5.5620 ± 204.1067i and its corresponding mode shapes are In comparison to eigenvalues at airspeed U = 10 m/s, the real parts of λ 1 and λ 2 are much closer, the real part of λ 1 is much more negative, and λ 2 and λ 3 are close to those for U = 10 m/s. In addition, at U = 100 m/s, the plunge component of mode shape ϕ 1 is much higher than that of U = 10 m/s. Furthermore, there are significant pitches in modes one and two, as shown in Figure 27. 0.1348 −0.3931 0.8565 where, in each mode shape, the first element presents plunge displacement, the second indicates pitch angle, and the last element gives control surface angle. Generally, in aeroelastic systems, the degrees of freedom are coupled to each other and cannot occur independently. Mostly, control surface and pitch displacements happen in mode two and three. In contrast to the regular wing, mode one contains significant pitch angle. The smart wing deformation in the three modes has been depicted in Figure 26. It is clear that there is almost symmetry in pitch in modes one and two; however, their values are significant. In comparison to eigenvalues at airspeed = 10 m s ⁄ , the real parts of 1 and 2 are much closer, the real part of 1 is much more negative, and 2 and 3 are close to those for = 10 m s ⁄ . In addition, at = 100 m s ⁄ , the plunge component of mode shape 1 is much higher than that of = 10 m s ⁄ . Furthermore, there are significant pitches in modes one and two, as shown in Figure 27.

Limitation of the Model
In this study, in order to simplify the problem, the structure was modelled by considering linearity in the structural behaviour. However, there is nonlinearity in the aerodynamic part. Therefore, although there is no nonlinearity in the structure, the whole model is nonlinear.

Limitation of the Model
In this study, in order to simplify the problem, the structure was modelled by considering linearity in the structural behaviour. However, there is nonlinearity in the aerodynamic part. Therefore, although there is no nonlinearity in the structure, the whole model is nonlinear.

Conclusions
In this paper, we showed how the flutter phenomenon can be postponed and even completely removed on a wing by using piezoelectric patches. The main contribution or finding that we've provided in this paper is a practical way to suppress flutter on a wing by implementing a passive aeroelastic control including piezoelectric patches and shunt circuits in which the size of required inductance is small. In Section 2, system response of a smart wing with only plunge DOF and pitch DOF has been presented. It is clear that using an efficient piezopatch can effectively decay the oscillations of the smart wing in a very short time. The smart wing vibration with only plunge DOF decayed in almost 0.6 s; however, the regular wing without a piezoelectric patch took around 12 s to decay. In addition, the oscillation of the smart wing with only pitch DOF decayed in almost 0.0198 s; however, the vibration of the corresponding regular wing took 0.3 s to decay. As illustrated in Section 3, implementing one piezopatch in the plunge DOF of a regular wing with three DOF can postpone the flutter speed by 81.41%, which is a considerable increase in the flutter speed. In addition, we show how the flutter phenomenon can shift from the plunge mode in a regular wing to the pitch mode in a smart wing. We also present the effect of adding one more piezopatch to a smart wing in the pitch DOF to further postpone the flutter phenomenon. The flutter speed in a smart wing can be postponed by 115.96%, which is a very considerable value. Furthermore, the flutter phenomenon disappears from the pitch mode; however, it remains in the control mode. Finally, adding one more piezopatch on a smart wing in the control DOF can completely remove the flutter phenomenon from the wing, which represents a great achievement in dynamic aeroelectic behavior of a wing.