Analytical modeling of panel flutter and active control in supersonic variable stiffness composite laminates

Abstract This work focuses on panel flutter and active control in supersonic variable stiffness composite laminates, making progress on the analytical modeling and combined exploration of curvilinear fiber composites tailoring and piezoelectric sensors/actuators, as promising structural design technologies, for aeroelastic control. The Classical Laminated Plate Theory and the First-Order Piston Theory are used as structural and aerodynamic models, respectively. Flutter analyses are carried out for simply supported plates, either purely elastic laminates or piezoelectric composite laminates. The tailor-ability of curvilinear fiber composites and the effect of proportional control are discussed. Ultimately, the presented results provide a comprehensive benchmark for future assessments.


Introduction
The multidisciplinary subject of smart composite structures, integrating piezoelectric sensors and actuators, has experienced a remarkable growth in terms of research in the last three decades, underlying the development of the next generation of multifunctional and adaptive structural systems designed to achieve active vibration control, noise attenuation, shape control or structural heath monitoring [1][2][3][4][5]. Likewise, the design of fiber reinforced composites, per se, has evolved up until the current point where variable stiffness composite (VSC) laminates, with curvilinear fiber paths, emerge as a resourceful and highly promising structural technology to make progress on the development of advanced composite structures for cutting-edge applications. In light of the increasing interest in smart structures and VSC laminates, this work focuses on panel flutter and active aeroelastic control in supersonic composite laminates, exploring curvilinear fiber composites to further improve the aeroelastic response behavior of constant stiffness composite (CSC) laminates, reinforced by conventional straight fibers. The assessment is provided making use of the Rayleigh-Ritz method, pushing forward the analytical modeling of panel flutter in advanced supersonic composites, while contributing to the limited number of available literature on the combined application of curvilinear fiber composites tailoring and piezoelectric sensors and actuators for aeroelastic control.
Active vibration control of structural systems resorting to piezoelectric elements has been intensively studied throughout the years [6][7][8][9][10][11][12][13][14] (citing just a few noteworthy works). As far as smart piezoelectric plate/shell composite structures are concerned, they are typically made of a composite laminate core with surface bonded piezoelectric layers or patches, polarized in the transverse direction, i.e. in extension mode. Therefore, when a transverse electric field is applied on the piezoelectric face layers, they undergo inplane deformations, forcing the host composite structure to bend accordingly [15,16]. The direct piezoelectric effect converts mechanical energy into electrical energy and provides sensing capabilities to the structure, whereas the inverse piezoelectric effect provides the actuation capabilities by converting the input electrical energy into mechanical deformation. Moreover, the control law can be defined, for example, in terms of the electric potentials of the sensor and actuator layers, applying then proportional-derivative (PD) feedback control [6,11,12] or a linear-quadratic regulator (LQR) based controller [8,14].
The application of piezoelectric sensors and actuators has been also explored for active control of flow induced vibrations and aeroelastic instabilities such as flutter [17][18][19][20][21][22][23][24][25]. In particular, the active panel flutter control in supersonic air-vehicles (e.g. supersonic jets, operational unmanned aerial vehicles and space rockets) aims to mitigate the occurrence of self-excited oscillations in the skin panels exposed to supersonic flow, which arises from the coupled interaction of inertial, aerodynamic and elastic forces. Thus, minimizing structural fatigue and failure, while broadening the flight envelope of the vehicle. The improved aeroelastic response of supersonic laminated plates integrating active piezoelectric layers is corroborated by the works developed by Song and Li [26][27][28], including both proportional and LQR controllers, as well as thermal buckling. Nevertheless, the design optimization of such structural systems is crucial to ensure both effectiveness and efficiency, including, for instance, the best placement of distributed sensors and actuators in the form of piezoelectric patches [29,30] and the selection of the associated control gains [31].
On the other hand, from the point of view of the evergrowing advances in composite materials science and automated manufacturing technologies, the highly promising concept of variable stiffness composites, with curvilinear fiber paths, has gained an increased attention within aeroelastic tailoring of wings [32][33][34] and skin panels [35][36][37]. These curvilinear fiber composite laminates, also known as variable angle-tow (VAT) laminates, clearly improve the tailor-ability of the conventional straight fiber composites due to a broader design space of the fibers control angles. Therefore, Guimarães et al. [38] focused on the optimization of smart VSC panels for active flutter suppression, applying PD control through surface bonded piezoelectric patches. It is shown that the combination of curvilinear fiber composites and active control can be exploited to achieve improved aeroelastic responses in terms of supersonic flutter. Nonetheless, both manufacturing limitations and defects cannot be disregarded for a complete and consistent design of curvilinear fiber composites in terms of stresses, failure and buckling [39][40][41]. In fact, it is worth remarking that refined structural models are also mandatory for an accurate analysis of VSC laminates [42][43][44], especially when a globallocal characterization is intended in moderately thick structural components.
Despite the noteworthy contributions by Guimarães et al. [37,38], which are mainly devoted to design optimization, very few papers in the literature deal purposely with the analytical modeling of panel flutter in supersonic VSC laminates and application of active control to further improve its aeroelastic stability. The present work aims to take a step forward on these topics and provide a comprehensive benchmark for future assessments. In line with the leading works on buckling and in-plane response behavior of VAT laminates by Olmedo and G€ urdal [45,46], in the early 90 s, analytical solutions for VSC laminates are herein derived resorting to the Rayleigh-Ritz method, considering simply supported boundary conditions. The proposed formulation uses the Classical Laminated Plate Theory (CLPT) [47] as structural model, while describing the aerodynamic pressure due to the supersonic flow by the well-established First-Order Piston Theory [48][49][50][51]. Actually, this is the same approach followed in prior researches [26-28, 37, 38]. Nonetheless, this paper presents some novelties with respect to the available literature, namely by focusing, from the analytical modeling point of view, on the case of VSC panels as well as the combined exploration of curvilinear fiber composites tailoring and piezoelectric sensor and actuator layers for aeroelastic control.
Numerical applications are divided in three parts and provide not only a numerical validation of the proposed formulation, but most especially a comprehensive assessment and comparison of the aeroelastic response behavior of supersonic composite laminates with curvilinear fiber paths (as opposed to straight), including its performance when applied to smart panels for active flutter control. Firstly, for both CSC and VSC laminates, the model validation is carried out by comparison with available solutions in the literature, including finite element (FE) solutions [44,52] in free vibration analysis, as well as Rayleigh-Ritz solutions [37] in panel flutter analysis. Secondly, to provide a complete understating of the effect of both fibers orientations and stacking sequences on panel flutter response, further analyses of CSC and VSC laminates of interest are provided for three different flow angles. Lastly, the active aeroelastic flutter control of supersonic smart composite laminates is presented for plates fully covered by piezoelectric layers, resorting to proportional control. To assess a prior validation of the full aero-electro-elastic model with available literature, these final numerical applications are based on [28]. However, this work ends up adding much to the original benchmark by including novel configurations with curvilinear fiber composites, as well as the effect of various proportional control gains. Hence, the potential of combining VSC laminates and piezoelectric layers for improving flutter performance of supersonic panels is discussed and the effect of the proportional control gain is thoroughly characterized. In addition, the aeroelastic response behaviors of the controlled systems are compared with the uncontrolled ones, highlighting design considerations to take into account during the selection of the control gain in order to ensure the efficiency and integrity of the smart panels in real applications.

Analytical aeroelastic model
As intended by this work, the Rayleigh-Ritz method is used to obtain analytical aeroelastic flutter solutions for supersonic multilayered panels, including VSC laminates and smart laminates with piezoelectric face layers. A schematic representation of a smart composite laminate under supersonic flow is shown in Figure 1, assuming the case of a composite core fully covered by piezoelectric skins. An arbitrary controller K links the electric potentials across the actuator (a) and sensor (s) piezoelectric layers, which are perfectly bonded to the upper and lower surfaces of the composite core (c), respectively. In general, the upper surface of the plate is subjected to yawed supersonic flow, with an in-plane angle K. The proposed analytical formulation makes use of the CLPT and the First-Order Piston Theory as structural and aerodynamic models, respectively.
Additionally, the involved layers can be either piezoelectric layers poled through-thickness (i.e. extension mode) or purely elastic layers, including both straight and curvilinear fiber composites. It is assumed that the materials of the different layers can be pointwise orthotropic, at most. Hence, the plane stress (r zz ¼ 0) constitutive equations [47], written in the global coordinate system (x, y, z), are given as follows: where r ij represents the stress components, e ii the infinitesimal normal strains, c ij the engineering shear strains, D i the electric displacement components and E i the electric field components. The material properties consist of reduced elastic constants Q ij , reduced piezoelectric coefficients e ij and reduced dielectric constants ij [47,53]. The reduced elastic constants are given explicitly in Appendix A according to the typical in-plane rotation matrix [47], as well as to multiple angles [54] rather than the the previous powers of sines and cosines (the two methods are discussed ahead within the scope of the present formulation for the case of VSC laminates). For variable stiffness composites with shifted curvilinear fibers, lying at a fixed distance between each other, the fiber angle distribution is a continuos function of the in-plane coordinates, i.e. h ¼ hðx, yÞ: Therefore, the reduced elastic constants in the global coordinate system are continuously graded as well, i.e. Q ij ¼ Q ij ðx, yÞ: In the present work, it is assumed only linear fiber angle variations, including the case of a sole distribution along the x-axis, as well as distributions along both x-and y-axes. Following the leading works by Olmedo and G€ urdal [45,46], the linear fiber angle distribution along the x-axis is parametrized by: where the two control angles < T 0 , T 1 > are T 0 ¼ hða=2Þ and T 1 ¼ hð0Þ ¼ hðaÞ: On the other hand, the linear fiber angle variation along both x-and y-axes is given by: where the four control angles are defined as T 00 ¼ hð0, 0Þ, T 10 ¼ hða, 0Þ, T 01 ¼ hð0, bÞ and T 11 ¼ hða, bÞ, being interpolated by linear Lagrange functions. The control angles are ordered as < T 00 , T 10 , T 01 , T 11 >: Throughout an entire VSC laminate, the layers can have different control angles or instead the same control angles, but being rotated around the origin according to a prescribed stacking sequence as in conventional straight fiber laminates.
Under the assumption of infinitesimal strains, the straindisplacement relations are given by: where ðu 1 , u 2 , u 3 Þ ðu, v, wÞ are the displacement components in the x-, y-and z-axis, respectively. The engineering shear strains are given by c ij ¼ 2e ij : In light of the Kirchhoff hypotheses for thin plates, within the CLPT [47], the displacement field of the smart laminate is given as shown: uðx, y, z, tÞ ¼ u 0 ðx, y, tÞ À z @w 0 @x (5a) vðx, y, z, tÞ ¼ v 0 ðx, y, tÞ À z @w 0 @y (5b) wðx, y, z, tÞ ¼ w 0 ðx, y, tÞ where the subscript 0 stands for the mid-plane location (z ¼ 0). From Eqs. (4) and (5), one can verify the transverse inextensibility of the CLPT, as well as the assumption of null transverse shear strains, i.e. e zz ¼ 0 and c xz ¼ c yz ¼ 0, respectively. Furthermore, it is assumed that the electric field is applied on the transverse direction of the piezoelectric layers by surface electrodes, which means a linear through-thickness distribution of the electric potential. Thus, the in-plane electric field components are null, i.e. E x ¼ E y ¼ 0, and the transverse component has a constant value through-thickness of E z ¼ /=h p , where / is the electric potential difference across the surface electrodes and h p the thickness of the piezoelectric layer. In addition, since the transverse shear deformations are neglected in the present formulation, one concludes from Eq. (1) that: (i) the transverse shear stresses are null, i.e. r xz ¼ r yz ¼ 0, even for piezoelectric layers and (ii) the only non zero electric displacement component is in the transverse direction, i.e. D x ¼ D y ¼ 0: Since the bending terms are the most relevant for free vibration analysis of thin plates, the membrane components u 0 and v 0 are neglected in Eq. (5). Therefore, assuming that the panels are simply supported in all four edges, the transverse displacement w 0 is expanded in a bi-sinusoidal Fourier series to fulfill the essential boundary conditions. Superposing the modal coordinates q mn , the Rayleigh-Ritz approximation of the transverse displacement is given as follows: A more compact form of Eq. (6) is achieved resorting to vectorial notation as shown: where W ¼ fw 11 ::: w 1N w 21 ::: w 2N ::: w M1 ::: w MN g T is a vector containing the bi-sinusoidal functions, and q ¼ fq 11 ::: q 1N q 21 ::: q 2N ::: q M1 ::: q MN g T the vector of modal coordinates. The Hamilton's principle [47] is used to derive the aeroelastic equations of the laminate, assuming that the kinetic energy T, potential energy U and work done by the aerodynamic forces W are given by: where X k and q k are the volume and the density of the klayer, respectively, while S is the surface subjected to the distributed aerodynamic pressure Dp: Besides, the superscript dot denotes the time derivative.
In line with the First-Order Piston Theory [49,50,55], the pressure exerted by supersonic fluid on the flat plate is given as follows: Dp ¼ Àk @w @x cos K þ @w @y sin K À l @w @t (9) where the first and second terms are related to the aerodynamic stiffness and aerodynamic damping, respectively. The flow angle, also known as yaw angle, is denoted by K as previously shown in Figure 1. Moreover, the dynamic pressure parameter k and the aerodynamic damping parameter l are given as k ¼ 1 , q 1 , U 1 and M 1 are the dynamic pressure, density, speed and Mach number of the free stream, respectively (noting that for validity ffiffi ffi 2 p < M 1 < 5 [5]). On the basis of the aforementioned assumptions, one can make use of the Hamilton's principle to derive the following aeroelastic equilibrium equations: where M is the modal mass matrix, C a the aerodynamic damping matrix, K the modal (purely elastic) stiffness matrix, K a the aerodynamic stiffness matrix and K / the electromechanical coupling matrix, while € q and _ q are the second and first time derivatives of the modal coordinates q (i.e. acceleration and velocity, respectively). Moreover, the electric potential applied across the actuator layer is represented by / a , which varies only in time. The detailed expressions of the modal matrices are given in Appendix B.

Aeroelastic equilibrium equations with active control
As followed in prior investigations by Reddy [6], as well as by Moita et al. [10] and Song and Li [26][27][28][29]31], the closed circuit charge measured through the electrodes of the sensor layer, i.e. the sensed charge, is obtained according to Gauss law, assuming negligible converse piezoelectric effect and the absence of external applied charges. Since the only non zero electric displacement component is the transverse one, the charge developed in the sensor layer is given in terms of the modal coordinates by: where z s is the transverse coordinate of the mid surface of the sensor layer. Considering the bottom piezoelectric layer as sensor, one concludes from Figure 1 that z s ¼ Àðh c þ h s Þ=2: Hence, one obtains the output electric potential of the piezoelectric sensor layer as shown: where N ¼ ðh s =ð 33 S s ÞÞU, remarking that 33 is the dielectric constant of the piezoelectric material of the sensor layer.
For plates fully covered by piezoelectric layers, the area of the sensor is S s ¼ ab.
To provide active control to the structure, each piezoelectric layer acts as a plate capacitor connected to a control circuit, which supplies power to the actuator depending on the sensor output ( Figure 1). For proportional-derivative feedback control in terms of the electric potentials across the sensor and actuator layers [11,26], the control law is given by: where G p and G d are the proportional and derivative control gains, respectively. Introducing Eqs. (12) and (13) in Eq. (10), one obtains the modal equilibrium equations for an actively controlled aeroelastic system as follows: Besides the dissipative damping terms, the nonsymmetric aerodynamic stiffness matrix (Appendix B) implies, per se, complex solutions for the aeroelastic equations of motion, which is due to the non-conservative nature of the aerodynamic forces. Hence, the general solution of the homogeneous differential equation (14) is in the form of q ¼qe st , where s is the complex eigenvalue related to the vibration frequency andq is the associated eigenvector corresponding to the vibration mode. Considering the non-trivial solution, the eigenvalues are determined by the following quadratic eigenvalue problem: where the real and imaginary parts of the eigenvalues are given by RðsÞ ¼ Ànx n and IðsÞ ¼ x ¼ x n ffiffiffiffiffiffiffiffiffiffiffiffi ffi 1 À n 2 p , denoting here the damped and undamped frequencies as x and x n , respectively, while n is the damping factor.
As the dynamic pressure parameter increases, flutter occurs when the damping factor of one mode turns from positive to negative [51] and therefore the system becomes dynamically unstable above the flutter dynamic pressure parameter k F . Furthermore, the eigenvalues emerge as complex conjugated pairs after reaching the flutter bound, while the frequencies of the modes involved in flutter coalesce. Since the aerodynamic damping plays a stabilizing role in the occurrence of supersonic panel flutter, leading to higher flutter bounds with respect to the case without including it [35], the exclusion of the aerodynamic damping contribution results in safer designs. In the absence of dissipative terms, the damping factor becomes negative at the same precise value of critical dynamic pressure parameter for which the modes coalesce.
For the upcoming numerical applications, the aerodynamic damping is neglected in order to assess a more conservative preliminary study, making solely use of proportional control (which is more effective than the derivative one in panel flutter control [26]). As a result, all damping terms are neglected and Eq. (15) is then reduced to: The aeroelastic equilibrium equations are coded resorting to symbolic computation within the commercial software MATLAB and the eigenvalues are extracted with the native eigensolver command. Moreover, the surface integrals are obtained analytically, except the ones regarding the modal stiffness matrix K (Appendix B) for variables stiffness composites with linear Lagrange distribution of the fiber angle, along both x-and y-axis, as given in Eq. (3), which are carried out numerically. To reduce the computational effort on the calculation of the modal stiffness matrix for VSC laminates-either analytically or numerically-the reduced elastic constants in the global coordinate system are calculated according to the method introduced by Tsai and Pagano [54]. In line with the expressions presented in Appendix A, this approach allows to reduce significantly the number of dependencies on the sines and cosines of the fiber angle with respect to the expressions derived from the rotation matrix and therefore the integration in-plane is achieved more easily.
It is further emphasized that the electromechanical coupling matrix K / vanishes for purely elastic laminates, as well as the dynamic pressure parameter for the case of free vibration in vacuum, i.e. k ¼ 0.

Numerical applications
As intended in this work, the aeroelastic flutter response of simply supported composite laminates in supersonic flow is investigated through a comprehensive study, including both straight and curvilinear fiber composites, as well as piezoelectric layers for active control. The numerical applications are divided in three subsections, encompassing the validation of the present model as well as further research on panel flutter and active control of supersonic VSC laminates. Firstly, the developed analytical formulation is validated in free vibration and panel flutter analysis by comparison with available solutions in the literature, including both CSC and VSC laminates. Secondly, to highlight the effect of the fibers orientations and stacking sequences on the aeroelastic flutter stability of square panels, the flutter bounds of various three-layer composite plates of interest are assessed. For the sake of completeness, these benchmarks are provided for three different yaw angles: K ¼ 0 and 90 for normal flow along the x-and y-axis, respectively, as well as K ¼ 45 for yawed flow. Lastly, the active flutter control of smart composite laminates fully covered by piezoelectric layers is presented for plates subjected to normal supersonic flow along the x-axis. The stabilizing effect of the proportional feedback control is discussed for four different laminates, involving either straight or curvilinear fiber composite layers. Hence, the uncontrolled and controlled aeroelastic responses are compared, while underlining the contributions of curvilinear fiber composites in active panel flutter control.
It is further emphasized that since all damping terms are neglected, the flutter bounds of the panels are estimated by the occurrence of the first coalescence of natural frequencies, as the dynamic pressure parameter k increases monotonically. The modes involved in the occurrence of flutter are explicitly indicated for each case, in addition to providing graphic representation for some noteworthy illustrative examples.

Model validation: free vibration and panel flutter analysis
The formulation and computational implementation of the proposed analytical model are firstly validated in free vibration analysis of simply supported CSC and VSC laminates by a comparison with FE solutions [44,52]. The test cases consist of three-layer square plates, with a ¼ b ¼ 1 m and side-to-thickness ratio of a=h ¼ 100, including a cross-ply laminate (0/90/0) as well as three VSC laminates with linear fiber angle distributions along the x-axis, in line with Eq.
(2). The material properties of the equal thickness composite layers are the following: E 1 ¼ 173 GPa, E 2 ¼ 7:2 GPa, G 12 ¼ 3:26 GPa, 12 ¼ 0:29 and q ¼ 1540 kg/m 3 . For each case, Table 1 shows the first eight predicted natural frequencies alongside converged FE solutions, increasing the number of bi-sinusoidal terms in the Fourier expansion (keeping the maximum number of terms in each in-plane direction equal, i.e. M ¼ N). In particular, the FE solutions are obtained with models based on shear deformation theories, namely Reddy's TSDT in [52] and layerwise FSDT in [44].
For the CSC laminate, each vibration mode corresponds to a sole bi-sinusoidal term (e.g. m ¼ n ¼ 1 for the fundamental mode) and therefore the convergence depends only if all first eight (m, n) modes are included. In the present case, Table 1 shows that with three terms in each direction, the seventh and eighth natural frequencies do not correspond to the expected ones. On the other hand, for VSC laminates, each natural mode shape is a linear combination of the assumed bi-sinusoidal trial functions. This difference is clearly highlighted by the fundamental modes represented in Figure 2, where one can notice that the curvilinear fiber composites, due to the distributed stiffness, exhibit a significant change in the mode shape with respect to the conventional cross-ply laminate. Hence, the assumed number of trial functions plays a direct role on the evaluation of all natural frequencies and mode shapes.
In view of the convergence analysis shown in Table 1, at least six terms in each in-plane direction are needed to provide converged solutions for the lowest natural frequencies of VSC laminates. Moreover, since the present formulation is free of transverse shear, in view of the CLPT, Table 1 also reveals that the natural frequencies are slightly overestimated with respect to the reference solutions, which are based on shear deformation theories. However, since thin plates are considered (i.e. a=h ¼ 100), the converged analytical solutions are, at most, 5% higher than the FE solutions. Hence, from a purely practical standpoint, the solutions can be considered in agreement. Nonetheless, one additional remark may be relevant, at this point, regarding the fact that for laminates with flexural anisotropy, involving bending-twisting coupling, the proposed bi-sinusoidal expansion functions may overestimate the stiffness [56], leading to higher values of natural frequencies. As a result, the accurate modeling of elastic coupling effects is considered of the utmost importance for the proper design and analysis of curvilinear fiber composites in further ensuing research.
Regarding the code validation on panel flutter analysis of both CSC and VSC laminated plates, two optimally designed configurations proposed by Guimarães et al. [37] are here considered. The simply supported panels consist of eightlayer rectangular plates, with a ¼ 0.6 m and b ¼ 0.4 m, having  [37], the supersonic flow is assumed along the x-axis and six sinusoidal terms are applied in each in-plane direction to approximate the transverse displacement. Figure 3 presents the evolution of the first two natural frequencies with the nondimensionalized dynamic pressure parameterk ¼ kða=hÞ 3 =G 12 : Additionally, the variation of the damping factor is also included to show that it is null below the flutter bound, turning negative as the natural frequencies coalesce. It is worth remarking that the present analytical formulation is exactly the same as considered by Guimarães et al. [37], applying bi-sinusoidal trial functions and neglecting the aerodynamic damping as well. Therefore, the predicted flutter bounds match precisely the values presented by Guimarães et al. [37]:k F ¼ 779 for the CSC laminate andk F ¼ 796 for the VSC laminate. As a result, the developed code is validated for panel flutter analysis of both straight and curvilinear fiber composite laminates.

Flutter analysis of VSC panels
To provide a complete understanding of the effect of both fibers orientations and stacking sequences on the supersonic flutter of square composite plates, as well as a comprehensive assessment of analytical aeroelastic flutter solutions for benchmarking, various CSC and VSC laminates are investigated in this subsection. It is assumed that the simply supported panels are subjected to normal supersonic flowaligned along either the x-or the y-axis, i.e. K ¼ 0 or 90, respectively-as well as to yawed flow, with K ¼ 45 , to further characterize and compare the aeroelastic response behaviors of the composite laminates.
The intended panels consist of three-layer thin plates, with a ¼ b ¼ 1 m and a=h ¼ 250, including both CSC and VSC layers. The equal thickness composite layers have the same material properties considered for the laminates investigated in free vibration analysis. In fact, the VSC laminates previously considered in free vibration analysis are also investigated here in terms of supersonic flutter, including three more VSC laminates (with linear fiber angle variations along the x-axis). Additionally, the VSC configuration with linear fiber angle variation along the x-and y-axes, considered for the model validation in flutter analysis, is included as well (however, in a three-layer form). As in the previous flutter analysis of rectangular plates, six terms in each inplane direction are used to approximate the transverse displacement, i.e. M ¼ N ¼ 6, ensuring converged solutions. Table 2 shows the first three predicted natural frequencies in the absence of flow, as well as the nondimensionalized flutter parametersk F ¼ k F ða=hÞ 3 =G 12 and the associated frequencies x F of each laminate. For normal flow along the x-axis, among all laminates considered, the (0/90/ 0) shows the higher flutter bound, followed closely by the VSC laminate with fiber angle distributions in both in-plane axes. Moreover, the symmetry of the problem is clearly shown by comparing the flutter margins of the (0/90/0) with respect to (90/0/90), for the three flow directions. More  precisely, the flutter bound of the (90/0/90) for K ¼ 90 is exactly the same of the (0/90/0) for K ¼ 0 , coinciding if K ¼ 45 : Overall, one concludes from the results presented in Table 2 that aligning the fibers of the outer layers with the flow direction increases significantly the flutter resistance of the panels. It is worth remarking that VSC laminates demonstrate a superior tailor-ability to achieve tradeoff solutions for different flow directions than the CSC laminates, including higher natural frequencies as well, e.g. the (<0,45>/ < À45,À60>/<0,45>) laminate.
Furthermore, for the (45/À45/45) and (<0,45>/< À45, À60>/<0,45>) composite laminates, Figure 4 presents the evolution of the first six natural frequencies with the nondimensionalized dynamic pressure parameter, considering supersonic flow along the x-axis. In addition, Figure 5 focuses on the case of the VSC laminate with fiber angle distributions in both in-plane axes, for K ¼ 0 and 45 : It is shown that flutter typically occurs with the coalescence of the first two modes, even thought it can appear among the second and third modes. Despite the existence of higher order flutter modes, i.e. later coalescences of natural frequencies, it is expected that the structural integrity might be severally compromised before reaching such extreme operational conditions. As perceived from the trends shown in Figure 5, the yaw angle has a great impact on the evolution of the frequencies with the dynamic pressure parameter, leading to completely different aeroelastic response behaviors in terms of flutter bounds and modes involved in the occurrence of flutter.

Active flutter control analysis of VSC panels
Based on the model formulation, the active aeroelastic flutter control analysis of both CSC and VSC laminates is carried out for simply supported panels fully covered by piezoelectric layers. It is assumed solely the case of normal flow along the x-axis (K ¼ 0 ), as well as variable stiffness composites with linear fiber angle variations along the same direction. The test cases consist of four piezoelectric composite laminates, which include two cases involving straight fibers and two other addressing curvilinear fibers. The complete stacking sequences, representing the piezoelectric layers by the letter p, are the following: Table 2. Nondimensionalized flutter dynamic pressure parameterk F and frequency x F (Hz) of simply supported CSC and VSC square plates, with a=h ¼ 250, under normal and yawed supersonic flow.   These numerical applications are based on the original benchmark by Li and Song [28], where the smart angle-ply laminate with composite core ½45= À 45 s is investigated (being here purposely considered for validation with available literature). The cross-ply laminate, with core ½0=90 s , is here included for improved flutter performance and further comparison. On the other hand, the two novel variable stiffness configurations are designed making use of the conclusions drawn in the previous analysis of purely elastic laminates, i.e. aligning the fibers of the outer layers with the flow direction (0 ) to enhance flutter resistance. It is worth remarking that the first VSC laminate has the same layup of the angle-ply laminate at the center of the panel, but with the fibers aligned with flow direction on the left and right sides. On the contrary, the second VSC laminate keeps the angle-ply layup on the left and right sides of the panel, aligning the fibers with the flow at the center. The comprehensive evaluation of the aeroelastic response of these new configurations ends up adding much to the original benchmark [28] by considering not only the emerging and highly promising curvilinear fiber composites, but also the effect of various proportional control gains (as opposed to solely unitary ones).
In line with Li and Song [28], the square plates have a side a ¼ b ¼ 0:1 m and a total thickness of h ¼ 0.0012 m. Each piezoelectric layer has thickness h p ¼ 0:0002 m and the composite core is divided in four layers of equal thickness. The material properties of the composite layers are the following: E 1 ¼ 150 GPa, E 2 ¼ 9 GPa, G 12 ¼ 7:1 GPa, 12 ¼ 0:3 and q ¼ 1600 kg/m 3 . Likewise, the material properties of the transversely isotropic and thickness-poled piezoelectric layers are the following: 12 ¼ 0:3, e 31 ¼ e 32 ¼ 22:05 C/m 2 , 33 = 0 ¼ 1695 and q ¼ 7600 kg/m 3 , where 0 ¼ 8:85 Â 10 12 F/m is the vacuum dielectric constant. Following the nondimensionalized form used in [28], the dimensionless dynamic pressure parameter is given byk ¼ ka 3 =D 0 , with D 0 ¼ 6:4 Nm. Firstly, Table 3 shows the convergence analysis of the aeroelectro-elastic model on the evaluation of the nondimensionalized flutter parameter of the smart angle-ply and the first smart VSC laminate, considering: (i) the uncontrolled system with null control gain G p ¼ 0; (ii) the controlled system with unitary positive feedback gain G p ¼ 1 and (iii) the controlled system with negative feedback gains G p ¼ À1, À 2 and À3: The results lead to the conclusion that negative feedback increases the flutter bound of the uncontrolled systems, as also shown in detail by Figure 6. Actually, the unitary negative control gain provides an increase in 10% to 15% of the flutter dynamic pressure parameter with respect to the uncontrolled system. On the contrary, the positive feedback control decreases the flutter resistance of the panels. For the smart angle-ply laminate, it is emphasized that the results shown in Figure 6 are in agreement with the solutions reported by Li and Song [28], where the same graph can be found.
Moreover, Table 3 reveals that the convergence is much slower as the value of the negative feedback gain increases in magnitude. This tendency is further highlighted in Figure  7, assuming À7 G p 0 and varying the maximum number of terms in the expansion of the transverse displacement. In view of the solutions convergence, the number of terms in the flow direction plays a more pronounced role on the evaluation of the flutter bound, i.e. M > N. Taking into account the results shown in Table 3 and Figure 7, ten terms in the x-direction and five in the y-direction-M ¼ 10 and N ¼ 5, respectively-ensure converged solutions within the range of control gains of interest.  Table 3. Convergence analysis of the nondimensionalized flutter parameterk F of smart composite laminates, under supersonic flow with K ¼ 0 , for various control gains. Furthermore, Table 4 presents the nondimensionalized flutter parameters of the intended smart composite laminates, applying M ¼ 10 and N ¼ 5. A graphic representation of the numerical results is also provided in Figure 8. The effectiveness of the active proportional control on the augmentation of the flutter bound is easily perceived. In fact, for the present cases, the active control can increase the flutter bound of the uncontrolled laminates up to 80%. However, it is worth noting that in the present configuration of the test cases, each piezoelectric layer has nearly 17% of the total volume of the laminate, which represents approximately 35% of the total mass. Therefore, a careful selection of the piezoelectric materials for the sensors and actuators, as well as their geometry and most efficient placement, in the form of distributed patches (instead of using an entire single layer), are mandatory to ensure feasible applications of active panel flutter control with a good compromise between the added weight of the active control system and its contribution to the aeroelastic stability augmentation.
Comparing the laminates, the smart angle-ply exhibits the lowest flutter bound for G p ¼ 0, as expected from the misalignment of the fibers with respect to the flow direction. More specifically, from the comparison of the two variable stiffness configurations, one concludes that aligning the fibers with the flow direction at the leading and trailing edge supports leads to higher flutter bounds than just steering the fibers with the flow at the center, especially when low proportional control gains are applied. Among the selected control gains, the VSC laminates show the maximum flutter resistance for G p ¼ À3, outperforming the conventional CSC laminates despite of the applied control gain. It is further emphasized that an optimal value of control gain is neither too low nor too high, as perceived from Figure 8. Actually, for low control gains-say G p < 2-flutter still occurs with the coalescence of the first two natural frequencies. Once the control gain is high enough, flutter appears at considerably higher values of the dynamic pressure parameter, being now due to the coalescence of the third and fourth modes. Nevertheless, above a certain limit value of control gain, the control system losses efficiency in the sense that with a lower gain, one can achieve either the same or a later coalescence of the third and fourth modes (minimizing as well the possibility of depolarization of the actuator layer under high voltage/electric field).
As a result, to ensure the efficiency and operability of the smart panels during different operational conditions, the control gain must be properly selected. However, it is pointed out that due to the theoretical nature of this work, which is based on a set of assumptions and simplifications, further refined analyses and experiments are necessary to provide a fully accurate assessment of the control authority in real design applications. As a limitation, the present formulation does not allow a complete characterization of both electromechanical and purely elastic coupling effects, in addition to neglecting the contribution of transverse shear deformations. Even so, some useful design considerations can still be drawn from the provided results, holding true regardless of the proposed analytical model, namely on the combined application of curvilinear fiber composites and proportional control for improving the aeroelastic performance of smart laminates in active panel flutter control. To be precise, the results lead to the conclusion that the variable stiffness can be tailored to enhance the effect of a given range of control gains (or a specific one of interest, instead), whereas the negative proportional control can be exploited to yield substantial improvements on the panel flutter stability of both CSC and VSC laminates. Hence, the simultaneous integration of properly tailored curvilinear fiber composites and optimized piezoelectric elements can make progress on the aeroelastic design of advanced supersonic panels, especially suited to engineer the most innovative and cutting-edge air-vehicles.
To further illustrate the effect of the proportional control gain on the aeroelastic response of supersonic panels, Figure  9 shows the variation of the first six natural frequencies of the [p/<45,0>/< À45,0>/< À45,0>/<45,0>/p] smart VSC laminate with the nondimensionalized dynamic pressure parameter, considering the uncontrolled system as well as the controlled one with increasing values of negative feedback gain. A careful examination of Figure 9 reveals that the negative control firstly delays the coalescence of the lowest two natural frequencies. Increasing the magnitude of the control gain, the first two modes uncouple and flutter occurs with the coalescence of the third and fourth modes for higher values of dynamic pressure.
Additionally, it is also shown in Figure 9 that a zero-frequency condition emerges if the control gain is too high. This condition represents a static instability, known as divergence within the scope of the aeroelastic instabilities. For G p ¼ À5, divergence appears under low dynamic pressure, but disappears as the pressure increases. Although not shown, the divergence zone increases for G p ¼ À6, while the coalescence of the third and fourth modes is anticipated. The same holds for G p ¼ À7, however, the gap in dynamic pressure between the disappearing of divergence and the occurrence of flutter is now very reduced (thus shortening the stable region where the structure can operate). In effect, these results further underline the paramount importance of a careful selection of the most suitable control gain for each operational condition, taking also into account appropriated modeling methodologies for an accurate analysis of smart panels and its aeroelastic design optimization.

Conclusions
In this work, an assessment of panel flutter and active control in supersonic variable stiffness composite panels is presented through a Rayleigh-Ritz model, taking a step forward on the analytical modeling and combined exploration of curvilinear fiber composites tailoring and piezoelectric sensor and actuator layers, as two emerging and highly promising structural design technologies, for aeroelastic control and suitable design of advanced supersonic panels. As a result, the present paper ends up adding much to thus far available in the literature by further exploring and benchmarking, from the analytical modeling point of view, the aeroelastic response behavior of supersonic VSC laminates and its active panel flutter control through piezoelectric elements.
The proposed analytical formulation makes use of the CLPT to describe the displacement field, while assuming the First-Order Piston Theory to provide the aerodynamic pressure exerted by the supersonic flow. More specifically, only bending deformations are retained and the transverse displacement is expanded in a Fourier series of bi-sinusoidal terms to fulfill the simply supported boundary conditions. The control law provides active stiffness to the structure by means of proportional feedback control. Numerical applications consider the panel flutter analysis of purely elastic composite laminates and smart composite laminates fully covered by piezoelectric layers. As intended, for each case, both straight and curvilinear composites are investigated and compared. It is worth remarking that the developed code is firstly validated by comparison with free vibration and panel flutter solutions available in the literature.
In terms of the composites tailoring, it is concluded that aligning the fibers of the outer layers with the flow direction provides a significant improvement on the flutter response of supersonic multilayered panels. Additionally, VSC composites with curvilinear fiber paths show a superior tailor-ability to achieve tradeoff solutions for different flow directions, as well as higher natural frequencies. Regarding the active flutter control analysis, the numerical results lead to the conclusion that variable stiffness technology can be exploited to improve the aeroelastic response of smart composite laminates in supersonic flow. In fact, even with linear fiber angle variations along the flow direction, the stiffness can be distributed either to increase the flutter bound for a range of control gains or to enhance the effect of a particular control gain of interest. On the other hand, the proportional control strategy involving negative feedback gains has indeed a great potential for increasing the flutter resistance of supersonic panels. Moreover, it is emphasized that the selection of the control gain plays a key role on the static and dynamic aeroelastic responses of supersonic smart panels. In particular, for low control gains, flutter still occurs with the coalescence of the first two natural frequencies in the present cases. Increasing the magnitude of the (negative) proportional control gain, the coalescence of the first two natural frequencies can be avoided and flutter occurs with a later coalescence of the third and fourth modes, pushing forward the flutter bound of the uncontrolled systems up to 80%. However, for even higher values of control gain, the smart panels can be statically unstable at low dynamic pressures and the coalescence of the third and fourth modes occurs slightly sooner, i.e. the flutter dynamic pressures are now closer to the ones obtained using a low control gain. Therefore, a proper selection of the control gain is crucial to ensure the efficiency and operability of the smart panels during the different operational conditions of the flight mission profile. In effect, the combination of purposely tailored curvilinear fiber paths and optimized piezoelectric sensors and actuators can be extremely valuable for the aeroelastic control of advanced and multifunctional supersonic panels, offering a broad range of feasible design variables within the underlying multiobjective design optimization problems.
Ultimately, the presented results provide a comprehensive benchmark for future assessments, including both purely elastic and piezoelectric composite laminates, either with straight or curvilinear fiber composite layers. Hence, this work may allow further ensuing research, namely on the refined modeling (analytical and numerical) of smart composite panels to render highly accurate active aeroelastic flutter control analyses and its design optimization in terms of curvilinear fibers and piezoelectric elements.

Disclosure statement
The authors declare that they have no known competing financial interests or personal relationships that could have appeared to influence the work reported in this paper.

Data availability
The data that supports the findings of this study is available from the corresponding author upon reasonable request.