1. Introduction
The aerospace industry continually seeks advancements in materials and design techniques to enhance the performance, safety, and efficiency of aircraft and unmanned aerial vehicles (UAVs). One of the most critical challenges in this domain is the phenomenon of flutter, an aeroelastic instability that can have catastrophic consequences if not properly managed. Flutter occurs when aerodynamic forces interact with a structure’s natural modes of vibration, leading to self-excited oscillations that can rapidly grow in amplitude. This phenomenon can cause significant structural damage or even complete structural failure, posing severe risks to the safety and operational capability of flying vehicles. Consequently, maximizing flutter speed—the critical speed beyond which flutter occurs—is paramount in the design of aerospace structures, particularly for control surfaces such as fins and wings.
Composite materials, particularly carbon/epoxy laminates, have emerged as the materials of choice in aerospace applications due to their exceptional strength-to-weight ratio, durability, and the ability to tailor their properties for specific performance criteria. The customization of fiber orientations and stacking sequences within laminates allows engineers to optimize mechanical properties such as stiffness and strength, which directly influence the flutter characteristics of the structure. By manipulating these parameters, it is possible to design structures that not only meet strength and stiffness requirements but also exhibit enhanced aeroelastic performance.
Minguet and Dugundji (1988) [
1] employed genetic algorithms to optimize composite wing structures for flutter suppression. Similarly, Liu and Atluri (2000) [
2] utilized a gradient-based optimization method to suppress flutter in laminated composite plates. Additionally, Park and Chai (1999) [
3] investigated the use of genetic algorithms for the optimal design of composite panels with constraints, while Matsumoto and Tanaka (2005) [
4] focused on gradient-based techniques for flutter suppression in composite laminates. Reddy and Rozvany (2001) [
5] conducted an algorithmic investigation of the flutter suppression of composite wings. Zhang, Zheng, and Liu (2017) [
6] conducted an optimization of composite wing structures using genetic algorithms. Guo (2005) [
7] conducted an optimization of composite wing structures for maximum flutter speed. These studies collectively demonstrate the application of various optimization algorithms to enhance flutter characteristics in aerospace structures.
The classical laminate theory (CLT) and the finite element method (FEM) are widely used to model and analyze the behavior of composite laminates. FEM, particularly when employing algorithms like the Lanczos algorithm, enables the precise calculation of modal frequencies in bending and torsion. These modal frequencies are critical for assessing flutter performance, as they determine the structure’s response to aerodynamic loads. However, identifying the optimal laminate configurations that maximize flutter speed is a complex task that requires a detailed understanding of the interplay between material properties, geometric configurations, and aerodynamic forces [
8,
9].
To address this challenge, this study integrates advanced regression algorithms with FEM analysis to systematically optimize laminate lay-ups. The process begins with the generation of an initial seed dataset of laminate configurations, which serves as the foundation for further analysis. This dataset is progressively expanded using bespoke software designed to interface with FEA modal tools, allowing for the creation and augmentation of a comprehensive design space. The primary objective is to maximize the Delta function, a parameter that correlates with flutter speed, by exploring various laminate scenarios.
A critical aspect of this optimization process is the evaluation of different regression algorithms to predict the Delta function values accurately. The algorithms considered in this study include Fast Forest Regression, Fast Tree Regression, Sdca Regression, and Lbfgs Poisson Regression. Among these, the Fast Tree Regression algorithm demonstrated the highest accuracy and computational efficiency, making it the algorithm of choice for subsequent analyses [
10,
11,
12,
13].
Using the Fast Tree Regression algorithm, a design space of nearly 2000 potential laminate candidates was explored. The focus was on symmetric lay-ups, as symmetric laminates are crucial for avoiding undesirable coupling between bending and torsion. This coupling can lead to complex deformation modes and reduced structural stability, which are particularly unfavorable in the design of control surfaces for UAVs and missiles.
The results from this extensive numerical simulation process are summarized in the final optimized laminate configurations, which exhibit the highest Delta function values and, consequently, the highest expected flutter speeds. These findings underscore the effectiveness of combining tailored composite materials with advanced optimization algorithms to achieve specific aerodynamic performance goals.
In summary, this study provides a robust and innovative methodology for optimizing the flutter performance of composite laminates. By leveraging the strengths of FEM analysis and machine learning techniques, we have demonstrated the potential to significantly enhance the aerodynamic performance of aerospace structures. The insights gained from this research can be directly applied to the design of safer, more efficient, and higher-performing aerospace vehicles, contributing to advancements in aerospace technology and engineering.
The research presented in this paper is critically important for advancing the design and performance of unmanned aerial vehicles (UAVs). UAVs are indispensable in modern aerospace applications, ranging from military surveillance to environmental monitoring, and thus require exceptional aerodynamic efficiency and structural robustness to operate effectively under diverse and often challenging conditions. Flutter, a dangerous aeroelastic instability that can lead to the rapid failure of structural components, represents a significant threat to UAVs’ safety and operational integrity. Addressing this issue is paramount to ensuring the reliability and longevity of these vehicles.
This study employs the Fast Tree algorithm to optimize the fiber arrangements and lay-up configurations of carbon/epoxy composite plates, widely used in UAV construction due to their high strength-to-weight ratio and versatility. The advanced optimization techniques deployed in this research enable the development of composite structures with significantly improved flutter characteristics. Specifically, the optimized configurations exhibit higher Delta function values, indicating increased flutter speeds and enhanced structural stability. These improvements mean that UAVs designed with these optimized components can fly at higher speeds and endure more strenuous conditions without succumbing to flutter-induced failures. Moreover, the integration of finite element analysis (FEA) and comprehensive computational flutter analysis (CFA) ensures that the optimized designs are not only theoretically sound but also practically viable. The use of commercial software like Nastran for validation further underscores the robustness of the proposed methodology. By pushing the boundaries of flutter optimization, this research provides a crucial framework that can be applied to the design of UAVs, making them safer, more efficient, and more capable of performing complex missions. In conclusion, the findings of this study are poised to make a substantial impact on the aerospace industry, particularly in the realm of UAVs. The enhanced flutter performance achieved through this innovative approach underscores the potential for significant advancements in UAV design, promising greater reliability and expanded operational capabilities in various aerospace applications.
2. Materials and Methods
Aeroelastic analysis employs a range of methods and tools, including numerical simulations, experimental testing, and optimization techniques, to assess and optimize the dynamic behavior of aerospace structures under aerodynamic loads. These methods play a critical role in ensuring the reliability, safety, and performance of aircraft, spacecraft, and other flexible systems operating in demanding environments.
2.1. Flutter Model Used for Optimization
To analyze the aeroelastic flutter characteristics of aircraft structures, it is essential to identify appropriate structural theories and aerodynamic models to establish mechanical models for structural flutter systems. This involves establishing partial differential equations for the structures using energy methods like Hamilton’s principle or the virtual work principle. Subsequently, these partial differential equations are discretized using methods such as the Galerkin’s method or finite element method (FEM), resulting in the equations of motion for discrete parametric systems [
14,
15,
16,
17,
18,
19].
The governing equation of an aeroelastic system, based on Newton’s law of motion can be expressed in the matrix form, as follows:
where [
M], [
C], and [
K] represent mass, damping, and stiffness matrices of the aeroelastic system,
x(t) is the system deformation, and
F(t) represents the aerodynamic force applied on the structure (system) analyzed; see
Figure 1.
In quasi-steady flutter, the frequency of vibration related to the flapping and twisting of the structure is considerably lower compared to its linear flight speed. This phenomenon is more common at supersonic and hypersonic flight speeds and is less frequent at subsonic and transonic speeds. Despite its simplicity in formulation, quasi-steady flutter provides an ideal starting point for optimizing the orthotropic structure of the flutter.
In the quasi-steady approach, it is assumed that the aerodynamic force is expressed by the following relation:
In the previous equation, CL is the lift coefficient assumed to be the function of the airfoil shape and the angle of attack α, ρ is the density of air at flight altitude, V0 is the linear flight speed, and l is the chord length (typical chord).
This assumption renders the possibility to express the flutter speed (
VF) for a quasi-steady regime, in a closed form, as a function of bending (
ωh) and twisting (
ωα) natural frequencies in the following form (
Appendix B), convenient for use in future optimization tasks:
In relation (3),
rα is the radius of gyration about the mid-chord,
μ is the dimensionless plate airstream mass ratio, and
xa is the relative distance between the aerodynamic and shear center (expressed as a fraction of the chord length). The term
, in relation (3), is the lift gradient; in this research, rectangular composite plates are analyzed, (rectangular cross section), for which the lift gradient is 2π/rad, the location of the aerodynamic center is at 0.25 aft the leading edge, and the shear center at half the chord
l, hence for the structure analyzed, the value of
xa is 0.25 (
Figure 1). Taking this into account, relation (3) can be rewritten and expressed as
Further, factoring the previous relation (relation (4)), to separate terms directly dependent on the composite fibers arrangement in the laminas (
m is the plate mass, ρ
m is the composite material density, and
t is the plate thickness), relation (4) becomes:
The term A*, in relation (5), contains all the terms independent of fiber arrangement factored from relation (4).
Analyzing the previous equation, one may observe that to maximize the flutter speed of the structure (which is always desirable), one of the parameters that require maximization is the squared difference of natural frequencies in bending and torsion:
In this paper, ∆ is designated as the Delta function. The objective is to maximize this function by analyzing the two-dimensional fiber arrangement (alignment) within the laminas and their stacking sequence in the composite laminate.
2.2. Optimization Algorithm Analysis and Results
Traditional methods of optimizing fiber alignment involve time-consuming trial-and-error processes or computationally expensive simulations. However, recent advancements in optimization algorithms have led to the development of efficient techniques, such as the Fast Tree algorithm.
The comprehensive optimization analysis conducted in this study is illustrated in the following diagram; see
Figure 2.
Initially, a seed dataset was established. Employing numerical methods, FEM modal analysis was conducted on multiple orthotropic composite plates featuring stack-ups commonly utilized in applications such as composite fins and UAV components. The FEM approach relied on the Lanczos algorithm for extracting modal frequencies in both bending and torsion, which were essential for calculating the Delta function. Initially, experimental modal analysis was performed on several laminates to validate the numerically obtained modal results. The numerical modal analysis was executed using commercial Abaqus software, focusing on an AS-4/3501-6 composite system (uni-directional), with dimensions of 70 × 35 × 1 mm, with different lay-ups.
The term “seed dataset” denotes an initial collection of data utilized to initiate the training process, initialize algorithms, or expand data through various means. This dataset was progressively augmented until the coefficient of determination (R2) reached a minimum threshold of 0.95. During this phase, several regression algorithms were scrutinized in terms of performance and R2 values. The algorithms examined included Fast Forest Regression, Fast Tree Regression, Sdca Regression, and Lbfgs Poisson Regression. To accomplish this objective, custom software was developed to interface with the FEA modal analysis software. This software allows users to generate various scenarios for creating and expanding laminate datasets.
The process of potential algorithm analysis (intermediate phase) is depicted in
Figure 3, while the final outcomes for all analyzed algorithms, including their performance measured in terms of R2 and training duration, are detailed in
Figure 3 and summarized in
Table 1.
From the findings presented in
Table 1, the Fast Tree Regression algorithm emerged as the preferred choice for subsequent analyses. The outcomes associated with this algorithm were derived from computations conducted using the finite element method (FEM approach) for 58 distinct lay-ups, encompassing symmetric, quasi-isotropic, asymmetric, and balanced configurations. Complete seed laminate lay-ups can be found in Annex A. The selection of this specific number of seed laminate candidates ensured the attainment of the targeted coefficient of determination of 0.95 for the optimization problem under investigation.
The Fast Tree algorithm is a powerful optimization tool that utilizes tree-based data structures to efficiently search through the vast solution space of fiber orientations. Unlike traditional algorithms, which may struggle with high-dimensional optimization problems, the Fast Tree algorithm excels in handling complex design spaces with numerous variables.
One of the key advantages of the Fast Tree algorithm is its ability to rapidly converge toward optimal solutions. By intelligently exploring the design space and exploiting problem-specific characteristics, the algorithm can identify promising fiber alignments in a fraction of the time required by conventional methods.
Moreover, the Fast Tree algorithm accommodates various design objectives and constraints, allowing engineers to balance competing factors such as weight, cost, and performance requirements. This flexibility makes it an invaluable tool for optimizing fiber alignment in carbon composite laminates tailored to specific applications and performance targets.
In practical applications, the Fast Tree algorithm can be seamlessly integrated into existing design workflows, enabling engineers to explore a wide range of design options efficiently. By harnessing the power of optimization algorithms like the Fast Tree algorithm, engineers can unlock the full potential of carbon composite laminates and drive innovation in lightweight, high-performance structures.
Progressing further, the dataset was expanded to nearly 2000 potential laminate candidates, representing the design space for the flutter problem optimization at hand. As previously mentioned, this expansion was achieved using specially designed software developed for this purpose (
Figure 4).
3. Results and Discussion
After completing the optimization task (maximization of the Delta function), the obtained laminate lay-ups (25 best candidates, with maximal values of the Delta function, and hence, the expected highest flutter speeds) are summarized in
Table 2. It is worth noting that the constraint to this flutter optimization problem is that only symmetric laminate lay-ups are considered. This is due to the fact that the potential application of this methodology can be applied in the design of stabilizing and the control surfaces of UAVs and missiles where the coupling between bending and torsion is unfavorable. Relying on the main composite materials advantage, that they can be tailored to a specific application, they can be designed in such a way that there is no coupling between bending and torsion. This is only achievable if the laminate is symmetric. In the case of symmetric laminates, all elements in a coupling stiffness matrix are equal to zero, and hence, no coupling [
9].
This is seen by analyzing the complete ABD composite matrix and laminate constitutive equations:
In the previous relation, N is the resultant in-plane force vector, M is the resultant moment vector, ε0 is the mid-plane strain vector, and κ is the curvature vector. The A matrix represents the in-plane stiffness of the laminate, relating the in-plane forces to the in-plane strains and accounting for the membrane behavior of the laminate. The B matrix represents the coupling stiffness, relating the in-plane forces to the out-of-plane bending and twisting deformations, thereby describing the coupling between in-plane and bending responses. The D matrix represents the bending stiffness of the laminate, relating the bending moments to the curvatures and accounting for the flexural behavior of the laminate.
The coupling stiffness matrix (B) is defined by the following relation [
8]:
where
Qk is the reduced stiffness matrix for the
kth layer, and
zk and
zk−1 are the vertical positions of the top and bottom surfaces of the
kth layer, respectively. For the symmetric laminate lay-up, [
B] is zero, and relation (5) becomes:
To validate the findings, Aeroelastic Finite Element Method (AFEM) models were constructed, and computational flutter analysis (CFA) was conducted to determine flutter speeds for missile fins composed of carbon/epoxy laminates, comparing unoptimized and optimized lay-ups resulting from this study. The aeroelastic model is illustrated in
Figure 5, comprising a structural model and an aerodynamic model interconnected with splines. In this model, the rear fuselage, where the fin is attached, is assumed to be rigid. The entire domain of interest is discretized with plate elements of a laminate type, while the aerodynamic model is discretized with aero panels [
20,
21]. The enhanced P-K algorithm is utilized to ascertain flutter speeds for both the optimized and non-optimized cases analyzed. The structural model (as a part of the complete flutter model) consists of plate elements (laminate elements). The solver was Nastran SOL 145; the laminate elements are utilized to model composite materials, which consist of multiple layers (laminae) with varying orientations and material properties. These elements are particularly useful in aerospace, automotive, and other engineering applications where composite materials are common due to their high strength-to-weight ratios and customizable properties. The PCOMP card in Nastran defines the properties of a composite laminate. This includes the material ID, thickness, orientation angle, and the stacking sequence of each lamina within the laminate. The laminate properties defined by PCOMP can be assigned to shell elements such as CQUAD4 (four-node quadrilateral) or CTRIA3 (three-node triangular). These elements are used to model the structural behavior of the composite laminate under various loading conditions. MAT8 is used to define orthotropic material properties for each lamina. This includes properties like Young’s modulus, the shear modulus, Poisson’s ratio in the principal material directions, and density.
In our analysis, we began with a coarse mesh and progressively refined it until two consecutive mesh sizes yielded nearly identical results. This approach, known as mesh convergence, ensures that the solution is not significantly affected by further mesh refinement, providing confidence in the accuracy of the results.
Starting with a coarse mesh allows for a quicker initial assessment of the problem, helping to identify any major issues without the computational expense of a fine mesh. As we refine the mesh, the solution becomes more accurate, but the increase in computational cost must be justified by the improvement in results.
This iterative refinement process is a standard practice when dealing with new problems, as it systematically ensures that the mesh density is sufficient for accurate results. Alternatively, mesh sizes can sometimes be determined based on previous experience (as in this research; see comment 1) with similar problems. When prior experience is available, it can guide the initial mesh design, reducing the number of iterations required to achieve convergence.
By combining these approaches, we ensure that our mesh is both efficient and effective, balancing computational resources with the need for accurate, reliable results.
Regarding the constraints, and considering that the optimization of composite plate flutter is crucial, only the plate itself was analyzed while the rest of the structure was assumed to be rigid. The plate was connected to the rigid boom using fixed constraints at the nodes.
The material elastic coefficients required for a complete definition are assumed to be known, obtained from the OEM, and are provided in the material characteristics table in the main text. If these coefficients were not available, they could be determined using well-established micromechanics theories (requiring fiber data, matrix data, and volume fractions) or through experimentation. This process would yield the four independent elastic coefficients necessary to define an orthotropic (thin-walled) structure.
Table 3 provides the geometry, material, and lay-up data for both configurations under analysis.
The results for the modal frequencies of the bending and twisting modes, obtained from the modal analysis of the non-optimized composite plate, are depicted in the following figures (
Figure 6 and
Figure 7).
Table 4 presents the modal analysis results, including the calculated Delta function, for all analyzed cases.
Upon analyzing the Delta function values, it becomes evident that the optimized lay-up exhibits significantly higher Delta function values, suggesting the likelihood of higher flutter speeds for this configuration. With this consideration, comprehensive flutter analysis for both lay-ups is conducted using the PK(NL) method implemented in the commercial Nastran software (Aeroelasticity module 1) [
20]. The results of the computational flutter analysis are presented in
Figure 8 and
Figure 9, along with V-g and V-f graphs (
Figure 10 and
Figure 11).
It is postulated that flutter, indicative of plate instability, occurs when the structure’s damping value diminishes to zero. Hence, from the graphs depicted in
Figure 10 and
Figure 11, it is evident that the damping for the non-optimal lay-up reaches zero at speeds of 23 m/s. In contrast, for the optimal configuration obtained through the methodology described in this research, the flutter velocity is 30 m/s, marking a notable increase of over 23%. This enhancement, attributed to the optimal fiber arrangement for flutter across all laminas composing the composite plate, signifies a significant improvement in its structural integrity and performance.