Research Article

Parametric Analysis of Position and Direction of Laminated Composite C-Spar on Aeroelastic Flutter in Aircraft Tail

Ahmad Reza Ghasemi ^a* , Hamid Rabieyan-Najafabadi ^b,
Hossein Nejatbakhsh ^a, Amin Gharaei ^c

^a Faculty of Mechanical Engineering, University of Kashan, Kashan, 87317-53153, Iran

^b Faculty of New Sciences and Technologies, University of Tehran, Tehran, Iran

^c Faculty of Engineering, Yazd University, Yazd, Iran

1.     Introduction

The aerodynamic loads on the wings, tails, and blades are often tolerated by the main spar structure, which is manufactured using laminated composite materials. Nowadays, C-spars are used in the tails to withstand the loads and moments. The spar web is made of laminated composite material with a high degree of multiaxial layups, and the spar flanges contain a high number of unidirectional laminates. The special airfoil section shape creates lift and drag aerodynamic loads, while bending and torsional moments are borne by the spars and skin, respectively. The shape, materials, layups, position, and direction of the spar have decisive effects on the performance, vibration, and static/dynamic behaviors of the wings, tails, and blades. The mutual interaction of aerodynamic, structural, and inertial forces (dynamic aeroelasticity) can cause catastrophic failure of the aircraft, known as flutter phenomena. As the wind speed increases in an aircraft, there may be a point at which the structural damping is insufficient to dampen out the increasing motions caused by the addition of aerodynamic energy to the structure. This vibration can lead to structural failure, and therefore, considering flutter characteristics is an essential part of designing aircraft parts.

Numerical, analytical, and experimental studies for better perception of the flutter phenomena of the aircraft tail have been presented in recent years [1]. Optimized the wing structure with aeroelastic constraints [2-3], and aeroelasticity tailored and behaviors of the wing box and blades [4-5] are some of the analytical and numerical investigations for a better understanding of this phenomenon in the blades. A number of studies focused on the effects of the initial in-plane and out-of-plane curvature [6-8] on the aeroelasticity instability (flutter) velocity, and some others used aviation standards for design and distinguish flutter speed point of the blades and wings [9-11]. The aeroelasticity (flutter)  analysis of swept aircraft wings [12-13] shows that concentrated mass has a complicated influence on the flutter boundaries.

The unsteady air loads affect the aeroelastic analyses using the doublet lattice method investigated by Van Zyl and Mathews [14]. The computational aeroelasticity methods for analysis of T-tail stability (flutter) boundaries for a free-flying aircraft in a transonic regime using a CFD formulation investigated by Attorni et al. [15]. Yu and Hu [16] have used an ultrasonic motor as an actuator in active flutter control of the airfoil section in the wind tunnel test. The aeroelasticity response of a bending–torsional coupling of wind turbine blade section with three degrees of freedom has been investigated by Stäblein et al. [15]. The structural stiffness of the airfoil section has assumed linear and bending–torsional coupling considered by the coupling coefficient in the stiffness matrix in their research [17]. Li and Ekici [18] using a one-shot method showed that the flutter boundaries can be very accurately predicted by prescribing a very small pitching amplitude.

Tang and Dowell [19] developed a numerical code for aeroelasticity analysis of the horizontal tail and mentioned that for a horizontal tail with flexible bending and torsion stiffness, a stall aerodynamic model and wing structural nonlinearity need to be considered. Zaki et al. [20] optimized the laminated composite arrangements to achieve maximum flutter with a minimum weight penalty. This research was one of the few studies in which the wing consisted of two spars, and the location of the spars was considered in the wing chord length for flutter investigation. Latif et al. [21] showed that changing the spar thickness contributes most significantly to the flutter speeds, whereas increasing the rib thickness decreases the flutter speed at high thickness values. Kumar et al. [22] modeled the geometric, structural, and aerodynamic parameters of the airfoil as independent variables. Furthermore, a review of the flutter of T-tail configurations [1] and a review of the nonlinear aeroelasticity of high aspect-ratio wings [23] showed that further work still needs to be conducted to identify aeroelastic instabilities that may occur in the tails and control surfaces.

Farsadi and Javanshir conducted a calculation of Flutter and Dynamic Behavior of Advanced Composite Swept Wings with Tapered cross-sections in Unsteady Incompressible Flow [24]. Touraj Farsadi et al. [25] studied the geometrically nonlinear aeroelastic behavior of pre-twisted composite wings modeled as thin-walled beams. Prasant Kumar Swain et al. [26] performed an aeroelastic analysis of a laminated composite plate with material uncertainty. Narayan Sharma et al. [27] conducted a stochastic frequency analysis of a laminated composite plate with curvilinear fiber. Narayan Sharma et al. [28] performed a free vibration analysis of a functionally graded porous plate using a 3-D degenerated shell element. Amin Gharaei et al. [29] proposed an analytical approach for the aeroelastic analysis of tail flutter.

The literature review of the flutter phenomena, especially the flutter speed of the aircraft tail, shows that many of the research studies modeled the airfoil skin and did not concentrate on the effects of internal components such as the spar. This may be because of the importance of the skin in torsional frequency, but the effect of the spar must also be considered. In this research, the aeroelastic behavior, including flutter and divergence speed, of an aircraft tail section consisting of C-spars, is investigated. The main purpose of this manuscript is to study the influence of the different locations and directions of the C-spar. The elastic section modulus (S=I/c) and the mass of the laminated composite spar (M_s) are considered constant. Geometrical and physical parameters are determined using the ANSYS commercial software, and a two-degree-of-freedom classical flutter analysis method is used in the analytical process. Furthermore, by changing the position of the C-spar from the leading edge of the airfoil to the mid-chord, and also by changing the direction of the spar toward and backward to the leading edge, the aerodynamic instability, including flutter and divergence speeds of the wing, is analyzed.

2.     Problem Descriptions

In the analysis of aeroelastic flutter, an aerodynamic model (incompressible 2D subsonic model) is developed to capture the complex airflow using computational fluid dynamics or wind tunnel testing. The structural model incorporates the mechanical properties of the aircraft tail. The solution to aeroelastic instability involves coupling the aerodynamic and structural models and employing numerical methods such as computational structural dynamics. By analyzing the coupled system under various conditions, critical flutter speeds and the stability of the tail can be determined. It is important to note that the specific models and numerical techniques employed may vary based on the complexity of the problem and available resources, with computational aeroelasticity being an advanced approach for capturing the intricate interactions between aerodynamics and structural dynamics.

In this research, the aeroelastic behavior of a V-tail aircraft with a C-spar has been investigated. The airfoil section with the position (a) and direction (b) of the C-spar beam is shown in Figure 1. The common position of the spar is changed from 10% to 50% of the chord length from the leading edge of the airfoil section, as shown in Figure 1(a). The direction of the C-spar, represented by the orientation of the spar towards the left hand (LH) and backward right hand (RH) to the leading edge, can be seen in Figure 1(b).

In the airfoil section with two degrees of freedom, as shown in Fig.2, the geometrical parameters are detected. The chord length is C and the half of the chord is determined as  in this Fig.. The points AC, SC, and CG refer to the aerodynamic center, shear center(elastic center) or reference point in the flutter analysis, and the center of gravity, respectively. The location of the points SC and CG are determined using dimensionless parameters  and  with respect to the mid-chord, as shown in Fig.2, which . When these parameters are zero, the point lies on the mid-chord, and when they are positive/negative, the points lie toward the trailing/leading edges. In the structural design process, the skin is designed with 4 composite layers with 45-degree layups. These layups will be constant in all cases.

Fig. 2. Tail section with two degrees of freedom and bending/torsional spring stiffnesses

The wing structural bending and torsional stiffnesses are modeled using discrete linear springs as shown in Fig. 2 [6-7, 30]. In the case of a V-tail with a smooth taper in the aircraft, the length of the mean chord at the mid-length of the tail is chosen for analysis. The bending and torsional stiffnesses are then determined based on the material properties and thickness of the skin.

3.     Governing Equations

For the flutter analysis of the wing, calculation of torsional stiffness and bending stiffness is necessary. Bending stiffness was determined by equation1 as follows [31]:

in the above equation, E, I_x, and L are Young’s modulus, a moment of inertia about the x-axis and length of the wing, respectively.

To determine the torsional stiffness, equation 2 is used [31]. In this equation, torsional stiffness is calculated by dividing the torsion of the wing, T by the torsional angle, .

The torsional angle of the wing without spar found by equation3 [31]:

A, G, and q are denoted to the airfoil area, shear modulus, and shear flow of the wing respectively. The  is circumferential segment lengths and t is the skin thickness.

Since the wing with spar is divided into two cells and each cell has a torsional angle, these angles must be equal to ensure structural integrity. Equation 4 is used to determine the angles of torsion in each cell [31].

the terms ,  represent summations around the entire perimeters of cell 1 and cell 2, respectively.  indicates the value of the interior web and these terms are calculated by equation 5:

The terms q1 and q2 are determined by equation6 as follows [29]:

By substituting the previously mentioned values into equation 2, the torsional stiffness can be determined. Similarly, the bending stiffness can be obtained using Equation 1. It's important to note that the equations mentioned are applicable to isotropic materials. In the case of orthotropic materials, an equal modulus approach must be used.

To find them, an orthotropic transformed compliance matrix  by φ angle ply for spar and skin layups must be evaluated. By and from compliance matrix, tension, and shear modulus ( ) can be found in equation 7 and equation 8. Substituting to equation1 and equation2, the required stiffness will be found [30].

The natural bending and torsional frequency are found by equation9 [22]:

4.     Aeroelastic Eigenvalue Problem Solution

To simplify the problem, we consider the following four dimensionless variables as follows [30]:

where r, σ, μ, and U are the dimensionless radius of gyration of the airfoil section about the shear center SC, the ratio of the uncoupled bending to torsional frequencies, the model mass to the mass of the air affected by the model, and reduced velocity as the dimensionless free stream speed of the air, respectively.

From Lagrange equations, the aeroelastic stiffness matrix can be written, that the determinant of this matrix must be equal to zero (equation11)

There are two complex conjugate pairs of roots, that . For a given airfoil section, the behavior of the real and imaginary parts of the complex roots as functions of U will be calculated and discussed [30]. In the next section, this formulation has been applied to a V-tail hollow airfoil section and a section with spar. The divergence and flutter speeds are investigated for different spar locations. The obtained results for this specific airfoil section are discussed and analyzed.

5.     Results and Discussions

5.1.  Geometrical Parameters Analysis

In this research, various configurations of spars are investigated. The suggested thickness of the skin is 1.0 mm. The chord length of the airfoil is C = 431 mm, and the aerodynamic parameters are determined and presented in Table 1 using the ANSYS commercial software. To determine the bending and torsional stiffnesses using the finite element method, a new approach is employed. The flexural wing section value is defined as the ratio of the force applied at the bending section to the deflection of that section relative to the root section. Similarly, the torsional stiffness of the wing is calculated as the ratio of the moment to the angular rotation of the airfoil section. The torsional moment is produced by a pair of vertical forces applied to a frame fixed at the measuring section [31].

Mechanical properties of the unidirectional and woven carbon fiber reinforced polymers (CFRP) that are used in the tail, have been mentioned in Table 2. The E, ν, and ρ denoted the elastic modulus, Poisson's ratio, and density, respectively; and subscripts 1, 2, and 12 denoted the on-axis directions of the laminated composite materials.

The bidirectional (BD) woven carbon fiber reinforced polymers (CFRP) with ±45° layups are used in skin and spar-web. In the spar cap, unidirectional (UD) and ±45° BD are used, simultaneously. The difference among cap layups is because of preventing from changing elastic section modulus (S=I/c). Also, cap width changes to be constant in the elastic section modulus S.

Table 1. Geometric dimensions of wing with and without spar, based on different spar locations

Table 2. Mechanical properties of CFRP laminated composite materials

In Fig. 3, a two-dimensional reference wing section with a spar is depicted. Bending and torsional loading were applied to determine the bending and torsional stiffnesses of the tail [33-34]. For the bending stiffness, the value of the tail stiffness is defined as the ratio of the applied force  at the bending section to the deflection ( ) relative to the root, as illustrated in Fig. 3(b). Similarly, for the torsional stiffness, the value is determined as the ratio of the bending moment  applied at the root section to the angle of rotation ( ) of the reference section relative to the root section, as shown in Fig. 3(c). Therefore, the bending/torsional stiffnesses can be written as:

In the finite element simulation depicted in Fig. 3, various spar locations and directions were considered, and the bending and torsional stiffnesses were obtained. This method, which involves using finite element analysis to determine the geometric properties and stiffnesses of the section, is referred to as the numerical method. In contrast, the analytical solution for these properties and stiffnesses is discussed in Section 3 of the research.

Fig. 3. Finite element simulation analysis to obtain a) Bending stiffnesses b) Torsional stiffnesses analyze

Table 3. Bending and torsional characterization (from theory) of airfoil sections for various spar configurations

Finite element models involve 40040 S4R elements. Bending and torsion loads are applied to each side of the tail symmetrically as boundary conditions. The static general finite element solver is employed for this purpose. A mesh study is conducted by using a mesh size ranging from 2 to 10 millimeters and results indicate that the mesh size has a minimal impact on the stresses with less than 8 percent.

The stiffnesses of the wing can be calculated analytically. Using this analytical approach, the frequencies and the ratio of the uncoupled bending to torsional frequencies have been determined for various types of sections. These results, obtained through the analytical method, are presented in Table 3. In addition to the analytical results, the finite element simulation results obtained using commercial software are also included in Table 3. A comparison is made between the numerical results and the analytical results, and the differences between these two methods are mentioned. It is stated that the observed differences between the results of these two methods are acceptable and within an acceptable range.

5.2.  Flutter Analysis of Tail Section

In the flutter analysis of orthotropic airfoil sections, the goal is to compare different spar configurations, locations, and positions. The steady-state method and thin airfoil theory for two degrees of freedom are utilized for this purpose. The flutter speeds in incompressible flow for different configurations are calculated and plotted using the governing equations. In this section, a specific section of CFRP (Carbon Fiber Reinforced Polymer) V-tail is considered for the flutter analysis. By substituting the CFRP parameters into the governing equations, the behavior of the complex roots as functions of the flight speed (V) is investigated to determine the divergent and flutter speeds. The imaginary part of the modal frequencies solution provides the primary graph of the flutter speed. The flutter point is located at the intersection of the torsional and bending (plunging) modal frequencies, where they coalesce. This point signifies the critical speed at which flutter instability occurs in the airfoil section.

5.2.1.          Flutter of Various Spar Locations

By substituting the mechanical properties of CFRP, as mentioned in Table 3, into the governing equations for different spar positions and locations as specified in Table 4, the divergences and flutter speeds of the composite sections are obtained. The plots of the imaginary parts of the roots versus dimensionless speed for various spar configurations, as well as without a spar, are shown in Figure 4 and Fig. 5. In these figures, it can be observed that the flutter speed occurs at the intersection between pitching and plunging oscillations, while the divergence occurs at zero frequency. As the distance from the leading edge increases, both the dimensionless divergence and the flutter speed decrease. These figures provide a visual representation of the effects of spar positions and locations on the stability characteristics of the composite sections. It demonstrates how different configurations and the absence of a spar can impact the occurrence of divergence and flutter.

Fig. 4. The imaginary part of the modal solution, indicates the modal frequency versus V for the right spar web

Fig. 5. The imaginary part of the modal solution, indicates the modal frequency versus V for the left spar web

Table 4. Flutter and divergence speed (Km/h) for various spar cap and web locations

In Fig. 5, as the distance from the leading edge increases, the dimensionless divergence and flutter speed decrease. This means that with a greater distance from the leading edge, the airfoil section becomes more stable, exhibiting a higher resistance to divergence and flutter. Regarding the conversion of dimensionless speeds to real speeds, multiplying the dimensionless speeds by the product  (where b represents a characteristic length and ω_θ represents the dimensional angular frequency) can provide an estimation of the real speeds. However, it is important to note that the behavior of dimensionless speeds may not directly correspond to the behavior of real speeds. The dimensionless speeds, when plotted, may exhibit a rational distance among them, allowing for easier comparison and analysis. However, this behavior may not be directly observed in the real speeds. The relationship between dimensionless speeds and real speeds is influenced by various complex factors, including aerodynamics and structural characteristics. Therefore, it is crucial to consider the dimensional values and carefully analyze the specific context and factors affecting the airfoil's behavior when interpreting and comparing dimensionless and real speeds.

5.2.2.          Discussion and Comparisons

Except for the Mach number, the effective parameters in flutter speed are x_θ, r, σ, and μ.

The frequency ratio (σ) is indeed influenced by the torsional frequency. In strength structural design, a straight wing is often assumed to behave like a beam subjected to aerodynamic bending loads. The spar is designed to withstand these aerodynamic loads while maintaining constant flexural stiffness. When changing the spar position, it is important to keep the flexural stiffness constant. This is achieved by maintaining a constant elastic section modulus (S=I/c) in all sections with a spar. As a result, the values of kh (the flexural stiffness factor) and ωh (the bending frequency) do not change.

1. a) Ratio of the uncoupled bending to torsional frequencies

1. b) Dimensionless radius of gyration of the airfoil
   section about the shear center SC

1. c) The model mass to the mass of the air
   affected by the model,

1. d) Reduced velocity as the dimensionless
   free stream speed of the air

Fig. 6 The imaginary Effective parameters versus spar position

However, moving the spar will affect the torsional rigidity. As the torsional frequency (ω_θ) increases, the frequency ratio (σ) decreases. This is because σ is defined as the ratio of the torsional frequency to the bending frequency. There are other parameters that can influence flutter, and these are plotted in Fig. 6. The position of the spar affects x_θ, which is the main factor influencing the inertia coupling stiffness (I_SC). Consequently, ISC affects r² (the square of the radius of gyration) and σ. It is important to consider these factors and their interactions when designing and analyzing the wing structure to ensure stability and prevent flutter.

In Fig. 7, as the spar is moved from the leading edge to the trailing edge, several observations can be made:

• The flutter speed (V_F) and the divergence speed (V_D) will decrease. This is due to the decrease in torsional stiffness resulting from the change in spar position.
• The distance between the flutter speed (V_F) and the divergence speed (V_D) will decrease. This reduction in the separation between VF and VD is generally not desirable, as it indicates a reduced margin of stability.
• The left-hand spar-web position (closer to the leading edge) shows a better condition to prevent flutter. This is contrary to what is shown in Figure 4 and Figure 5, where the results are presented in a dimensionless form. It is important to note that in wing design, the actual dimensional values should be considered to accurately assess the stability and flutter characteristics.

Regarding the design requirements, such as those outlined in the FAR23.629 standard, it is recommended to have a margin of at least 20% between the maximum aircraft design speed and the flutter speed. Additionally, it is common in design references to consider a minimum margin of 1.5 in the ratio between torsional and bending frequencies. In this research, with a 20% margin from the flutter speed, all the ratios are above 3, indicating a satisfactory design margin in accordance with these criteria.

It is crucial to carefully analyze the dimensional values and consider the specific design requirements and standards when interpreting the results and making design decisions to ensure the stability and safety of the wing structure.

Fig. 7 The flutter and speed versus the spar position

5.3.  Effect of Spar Direction

In strength structural design, the sizing of spar caps is an important consideration. The direction of the C-spar (the direction of the spar caps) may depend on the manufacturing process, but in structural design, it is typically not a critical factor. The position of the spar web (the vertical member connecting the spar caps) will affect the size of the torsional cell and shear flows within the structure. However, it does not directly impact the sizing of the skin thickness, which is primarily determined by buckling considerations. It is important for the designer to be cautious because the position of the spar web can significantly influence the aeroelastic characteristics of the wing structure. In this research, the effect of the spar web position is investigated. In Fig. 8, when the spar web is moved from the left-hand direction to the right-hand direction, the distance between the center of gravity and the elastic center will be considerably increased. This leads to an increase in the moment of inertia (ISC), which in turn decreases the torsional frequency (ωθ) and the frequency ratio (σ). Ultimately, this increase in moment of inertia and decrease in frequency ratio leads to an increase in the flutter speed (as shown in Fig. 6 and Fig. 8). These findings highlight the importance of considering the position and direction of the spar web in aeroelasticity considerations. The position of the spar web can have a significant impact on the torsional characteristics of the wing structure, affecting its stability and flutter behavior.

Fig. 8 Spars in one position (40%c) with different directions (different web positions)

6.     Validation

For validation of the method, an old approximate equation is used based on the bending–torsional coupling [34] that refers to NACA RM L7G02. It is rewritten according to this research symbols as follows:

where  is the ratio of the mass of air (for diameter cylinder equal to the chord of the wing) to the mass of the wing, and both of them are taken for the length equals the span of the wing. Also,  is the ratio of the mass-radius of gyration referred to shear center to the half-chord of the airfoil section.

Table 5. The discrepancy of flutter speeds for an aluminum section between results obtained from this research and the primary equation (Eq. 13)

The calculated flutter speed obtained in this research can be compared to the results obtained from Eq. 14. These results for an aluminum section for various thicknesses based on primary Eq. 14 and results of the coupling-torsional coupling in this research are presented in Table 5. These results show that this approach has good agreement with other methods, and the percent discrepancy between these two methods is less than 10 %. As the primary Eq. 14 is based on torsional stiffnesses and disregarding the bending stiffnesses, it seems that this discrepancy between results is not out of reality and is acceptable for validating this research results.

7.     Conclusions

Based on the results and conclusions obtained from the study on the V-tail aircraft using the steady method in incompressible flow, the following key findings were observed:

• This study presents a novel approach to designing and analyzing aircraft structures by incorporating composite lay-ups, isotropic relations, and FEM analysis. By considering these factors, designers can optimize the design and achieve the necessary levels of structural integrity.
• Adding a spar to a hollow section of the V-tail does not significantly improve the aeroelastic behavior. Only a modest increase of approximately 20% in flutter speed was observed. The primary effect of the spar is on the bending stiffness, which does not lead to a substantial increase in flutter speed.
• The position of the spar cap and the location of the spar web were investigated. It was found that locating the spar web on the left-hand side of the caps is better for improving the aeroelastic behavior. This finding is crucial for structural designers when determining the optimal spar configuration.
• Although flutter occurs before divergence, the distance between the flutter and divergence speeds can be considerable. Moving the spar from the leading edge to the mid-chord can decrease this margin, potentially leading to a reduced margin of stability.
• The study noted that the variations among flutter speeds for different spar configurations can be ignored to some extent. It was concluded that changing the spar configuration does not have a deterministic effect on flutter speed, with approximately 20% variation observed among the speeds. The calculations were based on conservative steady assumptions, and considering unsteady assumptions could result in considerably higher flutter speeds.
• The best configuration for achieving the maximum flutter speed was found to be when the spar cap and web spar were positioned at the maximum distance from the mid-chord (with the cap at 10% and the web on the left-hand side). Designers should take this into consideration when optimizing the flutter characteristics of the V-tail aircraft.

These conclusions highlight the importance of spar configuration and positioning in the aeroelastic behavior of the V-tail aircraft, providing valuable insights for structural designers to enhance stability and prevent flutter.

Overall, this article presents an applicable problem for aircraft designers who need to achieve standard levels for the structure to meet flutter FAR standard permissions.

Funding Statement

This research did not receive any specific grant from funding agencies in the public, commercial, or not-for-profit sectors.

Conflicts of Interest

The author declares that there is no conflict of interest regarding the publication of this article.

References

[1]   Murua, J., Martínez, P., Climent, H., van Zyl, L., & Palacios, R., 2014. T-tail flutter: Potential-flow modelling, experimental validation and flight tests. Progress in Aerospace Sciences, 71, 54-84.

[2]   Mastroddi, F., Tozzi, M., & Capannolo, V., 2011. On the use of geometry design variables in the MDO analysis of wing structures with aeroelastic constraints on stability and response. Aerospace Science and Technology, 15(3), 196-206.

[3]   Dillinger, J. K. S., Klimmek, T., Abdalla, M. M., & Gürdal, Z., 2013. Stiffness optimization of composite wings with aeroelastic constraints. Journal of Aircraft, 50(4), 1159-1168.

[4]   Stodieck, O., Cooper, J. E., Weaver, P. M., & Kealy, P., 2017. Aeroelastic tailoring of a representative wing box using tow-steered composites. AIAA journal, 55(4), 1425-1439.

[5]   Li, Z., Wen, B., Dong, X., Peng, Z., Qu, Y., & Zhang, W., 2020. Aerodynamic and aeroelastic characteristics of flexible wind turbine blades under periodic unsteady inflows. Journal of Wind Engineering and Industrial Aerodynamics, 197, 104057.

[6]   Piovan, M. T., Domini, S., & Ramirez, J. M., 2012. In-plane and out-of-plane dynamics and buckling of functionally graded circular curved beams. Composite Structures, 94(11), 3194-3206.

[7]   Amoozgar, M. R., & Shahverdi, H., 2019. Aeroelastic stability analysis of curved composite blades in hover using fully intrinsic equations. International Journal of Aeronautical and Space Sciences, 20(3), 653-663.

[8]   Amoozgar, M. R., Fazelzadeh, S. A., Khodaparast, H. H., Friswell, M. I., & Cooper, J. E., 2020. Aeroelastic stability analysis of aircraft wings with initial curvature. Aerospace Science and Technology, 107, 106241.

[9]   Shokrieh, M. M., & Behrooz, F. T., 2001. Wing instability of a full composite aircraft. Composite structures, 54(2-3), 335-340.

[10] Ghasemi, A. R., Jahanshir, A., & Tarighat, M. H., 2014. Numerical and analytical study of aeroelastic characteristics of wind turbine composite blades. Wind and Structures, 18(2), 103-116.

[11] Ghasemi, A. R., & Mohandes, M., 2016. Composite blades of wind turbine: Design, stress analysis, aeroelasticity, and fatigue. Wind turbines-design, control and applications, 1-26.

[12] Mazidi, A., & Fazelzadeh, S. A., 2010. Flutter of a swept aircraft wing with a powered engine. Journal of Aerospace Engineering, 23(4), 243-250.

[13] Mazidi, A., & Fazelzadeh, S. A., 2013. Aeroelastic modeling and flutter prediction of swept wings carrying twin powered engines. Journal of Aerospace Engineering, 26(3), 586-593.

[14] Van Zyl, L. H., & Mathews, E. H., 2011. Aeroelastic analysis of T-tails using an enhanced Doublet Lattice Method. Journal of Aircraft, 48(3), 823-831.

[15] Attorni, A., Cavagna, L., & Quaranta, G., 2011. Aircraft T-tail flutter predictions using computational fluid dynamics. Journal of Fluids and Structures, 27(2), 161-174.

[16] Yu, M., & Hu, H., 2012. Flutter control based on ultrasonic motor for a two-dimensional airfoil section. Journal of fluids and structures, 28, 89-102.

[17] Stäblein, A. R., Hansen, M. H., & Pirrung, G., 2017. Fundamental aeroelastic properties of a bend–twist coupled blade section. Journal of Fluids and Structures, 68, 72-89.

[18] Li, H., & Ekici, K., 2018. A novel approach for flutter prediction of pitch–plunge airfoils using an efficient one-shot method. Journal of Fluids and Structures, 82, 651-671.

[19] Tang, D., & Dowell, E. H., 2013. Computational/experimental aeroelastic study for a horizontal-tail model with free play. AIAA journal, 51(2), 341-352.

[20] Zaki, N. Z., Mazaha, F. A., Abdul-Latif, A., Mansor, S., Ab Wahid, M., Dahalan, M. N., & Nasir, M. N. M., 2019. Effect of Skin and Spar Laminate Orientations on Flutter of Composite UAV Wing. Journal of Aeronautics, Astronautics and Aviation, 51(2), 201-211.

[21] Latif, R. F., Khan, M. K. A., Javed, A., Shah, S. I. A., & Rizvi, S. T. I., 2020. A semi-analytical approach for flutter analysis of a high-aspect-ratio wing. The Aeronautical Journal, 125 1284 , 410-429.

[22] Kumar, S., Onkar, A. K., & Manjuprasad, M., 2020. Stochastic Modeling and Reliability Analysis of Wing Flutter. Journal of Aerospace Engineering, 33(5), 04020044.

[23] Afonso, F., Vale, J., Oliveira, É., Lau, F., & Suleman, A., 2017. A review on non-linear aeroelasticity of high aspect-ratio wings. Progress in Aerospace Sciences, 89, 40-57.

[24] Farsadi T., Javanshir J., 2017. Calculation of Flutter and Dynamic Behavior ofAdvanced Composite Swept Wings with Tapered Cross Section in Unsteady Incompressible Flow. Mechanics of Advanced Materials and Structures 1537-6494.

[25] Farsadi T., Rahmanian M., Kayran A., 2018. Geometrically nonlinear aeroelastic behavior of pretwisted composite wings modeled as thin walled beams. Journal of Fluids and Structures 83 (2018) 259–292.

[26] Swain P. K., Sharma N., Maiti D. K., Singh B. N., 2019. Aeroelastic Analysis of Laminated Composite Plate with Material Uncertainty. American Society of Civil Engineers, 2020, 33(1): 04019111.

[27] Narayan Sharma , Prasant Kumar Swain , D. K. Maiti & B. N. Singh, 2020. Stochastic frequency analysis of laminated composite plate with curvilinear fiber. Mechanics of Advanced Materials and Structures, 11 Aug 2020.

[28] Sharma N., Tiwari P., Maiti D. K., Maity D., 2021. Free vibration analysis of functionally graded porous plate using 3-D degenerated shell element. Composites Part C: Open Access 6, 2021, 100208.

[29] Gharaei A., Rabieyan-Najafabadi H., Nejatbakhsh H., Ghasemi A. R., 2022. An analytical approach for aeroelastic analysis of tail flutt. Advances in Computational Design, Vol. 7, No. 1, 2022, pp. 69-79.

[30] Hodges, D. H., & Pierce, G. A., 2011. Introduction to structural dynamics and aeroelasticity (Vol. 15). Cambridge university press.

[31] Peery, David J., 1983, Aircraft structures. Dover Publications.

[32] Kaw, A. K., 2005. Mechanics of composite materials. CRC press.

[33] Aleksandrowicz, R., & Lucjanek, W. Sailplane stiffness measurements. OSTIV Publications, 5.

[34] Budiansky, B., Kotanchik, J. N., & Chiarito, P. T., 1947. A Torsional Stiffness Criterion for Preventing Flutter of Wings of Supersonic Missiles. National Advisory Committee for Aeronautics Langley Field Va Langley Aeronautical Laboratory.