Open AccessArticle

Ply Optimization of Composite Laminates for Processing-Induced Deformation and Buckling Eigenvalues Based on Improved Genetic Algorithm

Qingchuan Liu

Xiaodong Wang

Zhidong Guan

Zengshan Li

and

Lingxiao Yang

School of Aeronautic Science and Engineering, Beihang University, Beijing 100191, China

Author to whom correspondence should be addressed.

Materials 2025, 18(2), 345; https://doi.org/10.3390/ma18020345

Submission received: 11 December 2024 / Revised: 3 January 2025 / Accepted: 6 January 2025 / Published: 14 January 2025

Download

Browse Figures

Versions Notes

Abstract

The structure of thermoset composite laminated plates is made by stacking layers of plies with different fiber orientations. Similarly, the stiffened panel structure is assembled from components with varying ply configurations, resulting in thermal residual stresses and processing-induced deformations (PIDs) during manufacturing. To mitigate the residual stresses caused by the geometric features of corner structures and the mismatch between the stiffener-skin ply orientations, which lead to PIDs in composite-stiffened panels, this study proposes a multi-objective stacking optimization strategy based on an improved adaptive genetic algorithm (IAGA). The viscoelastic constitutive model was employed to describe the modulus variation during the curing process to ensure computational accuracy. In this study, the IAGA was proposed to optimize the ply-stacking sequence of L-shaped stiffeners in composite laminated structures. The results demonstrate a reduction in the spring-in angle to 0.12°, a 50% improvement compared to symmetric balanced stacking designs, while the buckling eigenvalues were improved by 20%. Additionally, the IAGA outperformed the traditional non-dominated sorting genetic algorithm (NSGA), achieving a threefold increase in the Pareto solution diversity under identical constraints and reducing the convergence time by 70%. These findings validate the effectiveness of asymmetric ply design and provide a robust framework for enhancing the structural performance and manufacturability of composite laminates.

Keywords:

processing-induced deformation; finite element analysis; process simulation; genetic algorithm

1. Introduction

Composite materials are widely used in aerospace, automotive, and other fields due to their lightweight and high-strength characteristics. However, during the manufacturing process, composite materials inevitably generate thermal residual stresses, which affect the shape and dimensional accuracy of the structure after demolding, severely impacting its performance and causing assembly issues [1]. This problem is particularly significant for large and complex structures, where curing-induced deformation generates substantial assembly stresses, reducing the strength and fatigue life of the structure and potentially leading to structural failure [2]. The generation of curing deformation in composite laminated plates is primarily associated with temperature changes during the curing process, material anisotropy, and resin curing shrinkage. These factors interact with each other during the curing process, resulting in interlaminar residual stresses and causing complex stress fields, spring-in, and distortion deformation [3,4]. The residual stresses generated during the curing and cooling stages arise from both chemical and thermal effects [5]. Chemical shrinkage leads to resin contraction [6], while the laminate’s anisotropic mechanical properties due to the stacking sequence introduce differences in the in-plane and out-of-plane thermal expansion coefficients and chemical shrinkage coefficients of individual plies with different fiber orientations. Under the intense temperature fluctuations and complex chemical reactions during the curing process, in-plane and interlaminar residual stresses are generated, ultimately leading to curing deformation [7].

When the temperature decreases during the curing process of composite laminates, the layers will shrink in different ways due to inconsistencies in fiber orientation [8]. However, at this point, all the layers are bonded together, resulting in the generation of thermal residual stresses. For flat plate structures, the predominant deformation mode is distortion [9,10]; for structures with curvature features, such as L-shaped structures, the primary deformation mode is spring-in deformation [11]. However, for complex structures, their geometric characteristics exhibit both flat plate and curvature features. For example, a stiffened panel structure is a combination of flat plate and L-shaped structures, making the mechanism of thermal residual stresses even more complex.

To address the challenge of curing deformation in complex structures, researchers have proposed various methods, such as optimizing the curing process parameters using a manufacturer-recommended cure cycle (MRCC) [12], employing tools with thermal expansion coefficients similar to those of the prepreg material [13], and utilizing a ply design to fundamentally reduce the generation of residual stresses [1]. The utilization of a ply design is considered an effective strategy for actively controlling curing deformation [14]. By optimizing ply-stacking sequences and fiber orientations, especially through the application of asymmetric ply designs, residual stresses and deformations during the curing process can be significantly reduced, resulting in the improved overall performance of the composite material [15].

In the research and application of composite materials, ply optimization has garnered increasing attention as an effective design strategy. Ply optimization not only enhances the mechanical properties of composite materials, but also effectively addresses issues arising from curing deformation. Li employed asymmetric ply configurations to design laminates that resist hygrothermal shear deformation, achieving asymmetric ply designs for laminates comprising 8–14 layers [16]. Bruno et al. demonstrated that, by using “double–double“ layup (DD) stacking and adjusting fiber orientations, structural warping deformations can be reduced, proving that the depth of fiber orientation influences deformation [17]. Pelletier improved the strength and stiffness of materials while reducing the weight by optimizing the laminate stacking sequence and thickness [18]. Matsuzaki et al. utilized a finite element analysis and other numerical simulation methods to investigate the distribution of temperature and stress fields during the curing process, further achieving a reduction in curing deformation through an optimized ply design [19]. Cameron et al. developed a multi-objective optimization model that simultaneously optimizes the mechanical properties of materials and the deformation during the curing process, achieving more ideal design results [14].

In recent years, with the development of computer technology, the application of genetic algorithms (Gas) has gradually deepened in various fields. In the optimization of composite laminate stacking, genetic algorithms, with their strong global search capability, adaptability to discrete variables, and flexibility for multi-objective optimization, are highly effective in solving nonlinear and complex constraint optimization problems [20]. A GA is an optimization technique inspired by the principles of natural selection. As a population-based search method, it leverages the concept of “survival of the fittest” to iteratively refine solutions [21]. Some researchers [22] have applied genetic algorithms to optimize the stacking sequence of stiffened panels, improving their maximum buckling load. They have also discussed the influence of genetic parameters when using genetic algorithms to optimize stacking sequences for enhancing the buckling load. Choudhary et al. improved the buckling strength of composite cylindrical shells by 94% using genetic algorithms and a finite element analysis to optimize stacking sequences [23]. Huynh and Lee enhanced the laminate thickness design efficiency through a two-stage approach combining neural networks and genetic algorithms to predict the buckling factors [24]. Although genetic algorithms are commonly used for optimization, traditional Gas can easily become trapped in local extrema when solving more complex optimization problems [25], leading to unstable results and a low accuracy. Practical experience shows that traditional genetic algorithms often struggle to converge to a global optimum. As a result, many researchers have addressed these issues by improving the algorithm, and the introduction of the NSGA-II has greatly enhanced the stability of genetic algorithms. Srinivas and Deb [26] proposed the NSGA, a foundational method in multi-objective optimization. However, it lacks elitism, requires a sharing parameter, and has a high computational complexity. To address these issues, Deb et al. [27] introduced NSGA-II, which incorporates elitism and improves efficiency. Despite these advancements, NSGA-II struggles with diversity preservation in many-objective problems. To mitigate this, Luo et al. [28] proposed a dynamic crowding distance mechanism. Similarly, Coello and Pulido [29] developed the multi-objective micro-genetic algorithm (micro-GA), employing an archive for non-dominated solutions. Nonetheless, Pareto-based approaches often face challenges in many-objective optimization tasks [30].

Yuan et al. [31,32] proposed a multi-objective optimization method to reduce manufacturing defects and improve the efficiency in thick composite parts. Their method optimizes the curing process using Gas, combining the non-dominated sorting NSGA-II with a surrogate model, effectively reducing the curing time, the maximum temperature overshoot, and the maximum gradient of the degree of cure (DoC) simultaneously. Gao et al. [33] introduced a multi-objective genetic algorithm-based optimization strategy to address the multi-objective optimization problem in the curing process of ultra-thick composite laminates. By balancing the curing uniformity and residual stresses, their strategy achieved efficient and reliable process optimization. Fan et al. used a GA to optimize the ply layout and reduction strategy, balancing the structural weight, strength, and manufacturing performance, thereby improving the manufacturing efficiency and reliability of composite structures with ply-drop designs [34]. Kumpati et al. [35] presented a multi-objective GA for optimizing the stacking sequence of lightweight composite structures, with a focus on ensuring compliance with engineering design guidelines specific to the stacking sequence design.

The objective of this study was to solve the composite laminate stacking optimization problem through an improved adaptive GA (IAGA). First, the complex multi-field coupling processes involved in composite material manufacturing, including thermo-chemical, flow-compaction, and stress-deformation models, were analyzed, and a finite element analysis (FEA) method for the composite curing process was established. Then, based on the theoretical foundation of Gas, the role of the fitness function and genetic operators is discussed. Finally, considering the practical optimization problem—specifically, the characteristics of composite laminate stacking distributions—and the limitations of the original GA’s fixed crossover and mutation probabilities, as well as the constraints of NSGA-II, an adaptive GA was proposed. This improvement ensures a faster convergence and avoids local optima, while minimizing PIDs and residual stresses in composite components under the condition of maintaining the structural buckling stiffness, offering a more efficient and effective method for composite laminate design.

The Innovations of this paper are as follows: A high-fidelity numerical model of composite laminate PIDs, considering a simplified viscoelastic constitutive model, was established as the data input for stacking optimization, achieving a balance between data accuracy and computational efficiency. Based on theoretical derivations and the characteristics of the practical problem, an improved strategy for the GA was proposed to overcome issues such as convergence difficulties and susceptibility to local optima, which are inherent in traditional algorithms. Using an asymmetric and unbalanced ply multi-objective optimization approach, multi-objective optimization for both PIDs and the mechanical performance of the L-shaped stiffener was achieved. A characteristic index describing the influence of ply symmetry was summarized and validated with data, advancing the exploration of boundaries in composite material structure design.

In this paper, Section 2.1 introduces the relevant theories and mechanisms of PIDs, discussing the sources of PIDs, deriving the calculation of buckling eigenvalues, and providing related descriptions of the FEA. Section 2.2 analyzes the limitations of a traditional GA and proposes specific improvement strategies. In Section 3, the optimization results are analyzed and discussed, demonstrating the superior capabilities of the IAGA model. Finally, Section 4 presents the conclusion, containing clear data and summarizing the achievements of this study.

2. Theory and Methods

2.1. Curing Process Theoretical Analysis

Processing-induced deformation in composites is primarily composed of thermo-chemical, flow-compaction, and stress-deformation models.

2.1.1. Thermo-Chemical Models

The temperature and cure degree distribution during the curing process of composite materials is essentially a heat conduction problem with a nonlinear internal heat source, where the internal heat source originates from the exothermic chemical reactions of the matrix resin. In this study, the mathematical model for this problem is established using Fourier’s law of heat conduction and the energy balance relationship. Based on Fourier’s law of heat conduction and the curing kinetics equation, the coupled heat conduction and resin curing dynamics lead to the heat conduction–curing analysis control equations [9,36].

\frac{\partial}{\partial x} (K_{x} \frac{\partial T}{\partial x}) + \frac{\partial}{\partial y} (K_{y} \frac{\partial T}{\partial y}) + \frac{\partial}{\partial z} (K_{z} \frac{\partial T}{\partial z}) + \dot{q} = ρ c \frac{\partial T}{\partial t}

(1)

In these equations,

ρ

c

T

, and

t

represent the density, specific heat, temperature, and time, respectively;

K_{x}

K_{y}

, and

K_{z}

are the thermal conductivity coefficients of the material in the global coordinates; and

\dot{q}

is the internal heat generation rate per unit volume.

The curing kinetics part uses a phenomenological kinetic model and applies the Di-Benedetto equation to establish the relationship between the glass transition temperature and the cure degree of the QY9611 resin, as well as the glass transition temperature itself [37]. The equation is shown below, with the corresponding parameters listed in Table 1.

\frac{d α}{d t} = A \exp (\frac{- Δ E}{R T}) {(1 - α)}^{n}

(2)

T_{g} = T_{g, 0} + \frac{(T_{g, \infty} - T_{g, 0}) λ α}{1 - (1 - λ) α}

(3)

Here,

T_{g}

is the current glass transition temperature,

T_{g, 0}

is the initial glass transition temperature (uncured state), and

T_{g, \infty}

is the glass transition temperature at full cure (α = 1); λ is a parameter related to the curing conditions; α is the degree of cure; and R is the universal gas constant.

2.1.2. Resin Flow-Compaction Process

Based on Gutowski’s mechanistic model of resin flow/fiber interaction, prepreg laminates can be considered as a porous medium with nonlinear elastic deformation properties, filled with a viscous fluid. The excess resin that flows out of the composite material can be regarded as a seepage problem in a saturated porous medium. The flow behavior of the resin under curing pressure is described using the seepage theory of porous media. When the single-phase fluid in a porous medium is in a saturated state, assuming that the pressure during the autoclave curing process is constant, the total stress

σ_{i j}

can be decomposed into two parts:

σ_{i j} = {\underline{σ}}_{i j}^{f} - δ_{i j} P_{r}

(4)

where

{\underline{σ}}_{i j}^{f}

represents the effective stress in the fiber;

δ_{i j}

is the Kronecker delta function (where

δ_{i j}

= 1 when

i = j

and

δ_{i j}

= 0 when

i \neq j

and

P_{r}

is the resin pressure).

According to Darcy’s law, the compaction equation for the composite material is given by the following:

\frac{K_{x x}}{V_{f}} \cdot \frac{\partial^{2} P_{r}}{\partial x^{2}} + \frac{K_{y y}}{V_{f}} \cdot \frac{\partial^{2} P_{r}}{\partial y^{2}} + \frac{1}{V_{0}^{2}} \cdot \frac{\partial}{\partial z} (V_{f} K_{z z} \frac{\partial P_{r}}{\partial z}) = μ \frac{\partial}{\partial t} (\frac{1 - V_{f}}{V_{f}})

(5)

In the equation,

V_{0}

represents the fiber volume fraction when the load is zero;

V_{f}

is the fiber volume fraction;

μ

is the resin viscosity; and

K_{i j}

is the permeability of the composite material. The permeability is related to the fiber volume fraction, fiber diameter, and fiber structure, with an empirical formula given by the following:

K_{i j} = \frac{r_{f}^{2}}{4 K_{0}} \cdot \frac{{(1 - V_{f})}^{3}}{V_{f}^{2}}

(6)

where

r_{f}

is the fiber radius and

K_{0}

is the Kozeny constant, which varies depending on the fiber structure and resin flow direction. According to the Kozeny–Garman theory, the material porosity is expressed as follows:

e = \frac{V o l - V_{o l}}{V o l_{f}} = \frac{1}{V_{f}} - 1

(7)

where

V o l

and

V o l_{f}

represent the total volume of the composite material and the volume occupied by the fibers, respectively. In this study, the fiber volume fraction was set to 0.626. The process of resin flow in curing is shown in Figure 1. Figure 1 illustrates the curing process of a composite laminate under pressure. Initially, the resin is unevenly distributed with voids present. Upon applying curing pressure, the resin flows to fill the voids, creating resin-filled and unfilled areas. In the final stage, the pressure ensures complete resin penetration, resulting in a fully cured, compact, and void-free laminate.

2.1.3. Stress-Deformation Process

Considering both computational accuracy and efficiency, this study used a simplified resin constitutive model similar to the viscoelastic constitutive law [38], with the following expression:

E_{m} = \{\begin{matrix} E_{m}^{0} & T > T_{g} (α) \\ {[\frac{T_{c} - T}{T_{c} - T_{g} (α)}]}^{2} E_{m}^{0} + {[\frac{T - T_{g} (α)}{T_{c} - T_{g} (α)}]}^{2} E_{m}^{1} & T \leq T_{g} (α) a n d T > T_{c} \\ E_{m}^{1} & T \leq T_{g} (α) a n d T < T_{c} \end{matrix}

(8)

T_{c}

is the characteristic temperature. Equation (8) defines

E_{m}

as a piecewise function of the temperature

T

, describing the material’s modulus across different thermal regions. For

T > T_{g} (α), E_{m}

equals the rubbery modulus

E_{m}^{0}

. For

T \leq T_{g} (α)

and

T > T_{c}, E_{m}

is a weighted combination of the rubbery modulus

E_{m}^{0}

and glassy modulus

E_{m}^{1}

. For

T \leq T_{g} (α)

and

T \leq T_{c}

E_{m}

equals the glassy modulus

E_{m}^{1}

. This captures the transition between the glassy and rubbery states.

This model differs from CHILE in that, based on the experimental results [39,40], it is assumed that

T_{g} (α)

marks the onset of modulus increase, and no parameter

T_{c 1}

is required. According to experimental data [40], the resin modulus increases rapidly at first, and then slows down as the increase completes. Therefore, the modulus does not change linearly within the ranges

T \leq T_{g} (α)

and

T > T_{c}

, but rather exponentially. Based on the fitting of the RVE model calculation results, the exponent is taken as 2. The characteristic temperature

T_{c}

is calculated from the viscoelastic test data. The reaction kinetics and glass transition temperature of QY9611 resin are described by the following [37]:

\frac{d α}{d t} = 2.28 \times 10^{7} e x p (\frac{- 11114.9}{T}) (1 - α)^{1.485}

(9)

T_{g} = - 264.04 + \frac{120.46 α}{1 - 0.5176 α}

(10)

Figure 2 illustrates the development of residual stresses during the curing process due to the mismatch in thermal and chemical shrinkage between plies with different fiber orientations. This mismatch causes interlayer residual stresses, which manifest as curing-induced deformation in the laminate.

2.1.4. Calculation of Buckling Eigenvalues

The composition and stacking sequence of laminate ply ratios significantly affect the structural load-bearing performance. Given that stiffened panels exhibit a higher sensitivity to compressive loads, their mechanical performance is evaluated by calculating the buckling eigenvalues, assessing the structural stability under optimized curing deformation.

The modification of laminate stacking sequences has a significant impact on the load-bearing performance of the structure. Since stiffened panels are highly sensitive to compressive loading, their mechanical behavior is assessed by calculating the buckling eigenvalues, which serve as a crucial indicator of the structural stability.

The stiffness matrix

[D]

of the laminated plate is calculated based on classical laminate plate theory (CLPT):

[D] = \sum_{k - 1}^{N} [Q]^{k} (h_{k} - h_{k - 1})

(11)

The matrix

[Q]^{k}

represents the reduced stiffness of the

k

-th layer, where

h_{k}

and

h_{k - 1}

denote the positions of the top and bottom surfaces, respectively.

N

is the total number of layers.

Under plane stress conditions, the buckling differential equation of laminated plates is shown in Equation (12).

D_{11} \frac{\partial^{4} w}{\partial x^{4}} + 2 (D_{12} + 2 D_{66}) \frac{\partial^{4} w}{\partial x^{2} \partial y^{2}} + D_{22} \frac{\partial^{4} w}{\partial y^{4}} + N_{x} \frac{\partial^{2} w}{\partial x^{2}} + N_{y} \frac{\partial^{2} w}{\partial y^{2}} + 2 N_{x y} \frac{\partial^{2} w}{\partial x \partial y} = 0

(12)

The boundary condition is fixed at one end.

\begin{matrix} w (0, y) = 0 \\ \frac{\partial w}{\partial x} (0, y) = 0 \end{matrix}

(13)

In Equation (12),

w

is the deflection function,

D_{i j}

represent elements of the stiffness matrix of the laminated plate, and

N_{x}, N_{y}, {a n d N}_{x y}

are the in-plane stresses.

To solve differential Equation (12), the deflection function is set as a double sine function:

w (x, y) = W s i n (\frac{m π x}{a}) s i n (\frac{n π y}{b})

(14)

In Equation (14),

W

represents the amplitude;

m, n

are half-wave numbers; and

a, b

are the length and width of the laminated plate.

The deflection function can be substituted into the buckling differential equation to obtain the eigenvalue equation for the critical load:

[D_{11} {(\frac{m π}{a})}^{4} + 2 (D_{12} + 2 D_{66}) {(\frac{m π}{a})}^{2} {(\frac{n π}{b})}^{2} + D_{22} {(\frac{n π}{b})}^{4}] + N_{x} {(\frac{m π}{a})}^{2} + N_{y} {(\frac{n π}{b})}^{2} = 0

(15)

N_{x, c r} = λ_{b u c k l i n g} \times N_{r e f}

(16)

The symbols

D_{11}, D_{12}, D_{22}

, and

D_{66}

represent the elements of the laminate’s bending stiffness matrix. Specifically,

D_{11}

characterizes the bending stiffness of the plate primarily in the

x

-direction and

D_{22}

does so in the

y

-direction, while

D_{12}

and

D_{66}

describe the coupling effects between different bending modes and the in-plane shear stiffness, respectively. Solving the above system of equations yields the buckling eigenvalue when

N_{r e f} = 1

2.1.5. FEA Model and Mechanical Parameters

Considering that the primary source of deformation in stiffened panels, widely used in aircraft skin structures, is the spring-in deformation of the stiffeners causing overall structural deformations, and that stiffeners play a crucial role in bearing compressive loads [41], this study used the ABAQUS 2019 software to simulate the spring-in angle and buckling eigenvalues of a typical L-shaped laminate stiffener structure and built a GA platform using SIMULIA Isight software. Due to the limitations of shell elements in capturing the mechanical behavior in the thickness direction of the material, three-dimensional solid elements were used to model and compute the PIDs of the composite laminate [42]. The material properties were implemented through the user-defined subroutines UMAT and UEXPAN. The laminate structure was modeled using C3D20R elements, while the tool was modeled using C3D8R elements. The FEA model used in this study was validated in a previous paper published by other researchers from our research group [38]. The FEA model in this study was established with the same material parameters, boundary conditions, and loading settings, with the different dimensions as shown in Figure 3. The curing process and DoC curve are shown in Figure 4.

The mechanical properties and CTEs of ZT7H/QY9611 in a glassy state were tested, as shown in Table 2.

The FEA models are shown in Figure 5 and Figure 6. For the FEA model of PIDs, the applied boundary conditions included the curing pressure on the surface of the L-shaped structure and frictional contact with the mold. The loading temperature curve is shown in Figure 4. After the heat preservation stage, the mold was removed using the “model change” function in Abaqus, while the L-shaped structure was fixed to allow free spring-in.

For the buckling analysis, the first-order buckling mode can clearly assess the buckling behavior of the L-shaped structure. Therefore, the buckling eigenvalue of this mode was used as the evaluation criterion. As analyzed earlier, the buckling analysis had a clear physical significance, and its accuracy was relatively higher compared to the simulation of the complex nonlinear process of curing deformation. Therefore, the trend of change in the buckling eigenvalues during the optimization process was used as the optimization objective.

2.2. Improved Genetic Algorithm Based on NSGA-II

The aforementioned FEA established an accurate mapping relationship between the ply angles, stacking sequences, PIDs, and buckling stiffness, providing a numerical foundation for the development of subsequent optimization models. A set of 16 plies with specific angles and stacking sequences was utilized as the inputs for the optimization model, with the PIDs and buckling stiffness defined as the optimization objectives.

2.2.1. Limitations of NSGA-II in Ply Optimization

The primary distinction between the NSGA-II algorithm and a standard genetic algorithm lies in the hierarchical approach applied in NSGA-II before executing the selection operator based on the dominance relationships among individuals [43,44]. By employing a non-dominated sorting method, NSGA-II ensures that superior individuals have a higher probability of being passed on to the subsequent generation. Additionally, NSGA-II incorporates an elitism strategy, combining the parent and offspring populations to facilitate competition and produce the next generation. The standard NSGA-II process is shown in Figure 7. This figure illustrates the core workflow of NSGA-II. The parent population

P_{t}

and offspring population

Q_{t}

are merged into

R_{t}

followed by non-dominated sorting to divide

R_{t}

into ranked solution sets (

L_{1}, L_{2}, \dots

). Solutions are selected based on rank and, if necessary, crowding distance to ensure diversity, forming the next-generation population

P_{t + 1}

. This process balances the solution quality and diversity in multi-objective optimization. The core features are as follows: (1) the introduction of a fast non-dominated sorting algorithm; (2) the utilization of crowding distance and crowding comparison operators to maintain diversity; and (3) the implementation of an elitism strategy to preserve the best solutions.

For the stacking optimization problem addressed in this study, further modifications were made to the algorithm. These improvements mitigated the convergence challenges caused by the high input dimensionality and prevented the algorithm from becoming trapped in local optima, which can result from the highly nonlinear characteristics of a PID analysis.

2.2.2. Adaptive Non-Dominated Sorting

In traditional NSGA-II, the handling of non-dominated sorting is static. However, during ply optimization, as the optimization progresses, the distribution of the population composed of 16 layers with different angles and the numerical optimization of PIDs and

λ_{b u c k l i n g}

continuously change. Therefore, this study introduces adaptive non-dominated sorting, dynamically adjusting the tightness of sorting based on the generational progress and changes in the objective space. This adjustment balances the maintenance of diversity and the convergence of solutions.

In the early stages, to preserve diversity, a relatively relaxed dominance relationship among solutions is allowed. In the later stages, as the solutions approach the Pareto front, the constraints on dominance relationships are gradually tightened, promoting the convergence of the solution set toward the Pareto front. This improves the algorithm’s flexibility in handling multiple ply input conditions for laminates. The factor tightness (t) increases progressively with the number of generations t, making the sorting increasingly strict:

T i g h t n e s s (t) = \frac{1}{1 + e x p (γ \cdot (t - t_{m a x}))}

(17)

Here,

γ

is a constant that controls the rate of increase in tightness, and

t_{m a x}

is the maximum number of generations. When

t = 0

, the

t i g h t n e s s (t)

is relatively small (approaching 0), indicating that, in the early stages, the sorting rules are more relaxed. As the number of generations increases, the

t i g h t n e s s (t)

gradually grows, making the sorting rules stricter and promoting the convergence of individuals toward the Pareto solution set. The crowding distance calculation formula is modified as follows after introducing the tightness factor:

d_{i} = \sum_{k} (\frac{f_{k} (x_{i + 1}) - f_{k} (x_{i - 1})}{f_{m a x} (k) - f_{m i n} (k)}) \cdot T i g h t n e s s (t)

(18)

Here,

f_{k} (x)

represents the value of the objective function, while

f_{m a x}

and

f_{m a x}

are the maximum and minimum values of the objective function

f_{k}

, respectively.

2.2.3. Disturbance-Driven Crossover

Traditional crossover operators tend to become conservative when handling complex nonlinear problems, making them prone to becoming trapped in local optima. By introducing a perturbation mechanism, the randomness of the crossover process is increased, maintaining the diversity of the ply population and enhancing the algorithm’s exploration capability. Additionally, the perturbation intensity decreases linearly with the number of generations to ensure convergence.

For two parent individuals

p_{1}

p_{2}

, and

p_{3}

, the new individual after crossover can be adjusted by incorporating perturbation as follows:

p_{3} (i) = (ξ p_{1} (i) + (1 - ξ) p_{2} (i)) + η \cdot N (0, σ^{2})

(19)

Here,

η

is the perturbation intensity,

N (0, σ^{2})

represents standard normal distribution noise, and

ξ

is the crossover coefficient, ranging between 0.7 and 0.9.

2.2.4. Dynamic Gaussian Mutation Operator

In multi-objective optimization problems, the algorithm needs to balance exploring the solution space (diversity) and exploiting known high-quality solutions (convergence). Traditional mutation operators, particularly those with fixed Gaussian mutation amplitudes, may lead to premature convergence (lack of diversity) or the insufficient exploration of the solution space (excessive mutation), thereby reducing the algorithm’s performance.

By replacing the mutation amplitude’s relationship with the generation t with an exponential decay function, the algorithm adopts a larger mutation amplitude in the early stages, which gradually decreases over time. The variation in the mutation amplitude is expressed by Equation (20).

σ (t) = σ_{0} \cdot {(1 - \frac{t}{T})}^{β}

(20)

Here,

σ_{0}

represents the initial mutation amplitude,

T

is the maximum number of generations, and

β

is the decay coefficient. As the number of generations increases, the mutation amplitude gradually decreases, reducing the disturbances to already optimized solutions.

During the optimization process, the mutation amplitude is dynamically adjusted based on the convergence status of the objectives. If a particular objective has already converged, its mutation amplitude is reduced to avoid the excessive exploration of that objective.

σ_{k} (t) = σ_{0} \cdot {(1 - \frac{t}{T})}^{β} \cdot (1 + α \cdot V a r (f_{k}))

(21)

Here,

V a r

represents the variance of objective

k

within the population, and

α

is the adjustment coefficient. The Gaussian mutation formula for generating a new solution

x^{'}

is as follows:

x_{i}^{'} = x_{i} + N (0, σ_{k} (t)^{2})

(22)

N (0, σ_{k} (t)^{2})

represents a normal distribution with a mean of 0 and a variance equal to the dynamically adjusted mutation amplitude

σ_{k} (t)

The above strategy not only improves the global search capability of NSGA-II, but also enhances the local convergence speed, thereby increasing the efficiency and accuracy of the algorithm when handling complex multi-objective problems like ply optimization. The process flowchart of the IAGA model is shown in Figure 8.

2.2.5. Optimization Objectives and Constraints

The stacking design constraints for laminated plates in engineering have been summarized by researchers [45]. In addition to the symmetry typically required in traditional laminate designs, the 10% rule for ply distribution must also be maintained. Specifically, for the 0°, ±45°, and 90° directions, at least 10% of the layers must be oriented in each direction.

Considering that the primary optimization objective is curing deformation, appropriate weights are assigned to focus on minimizing the curing deformation. A target aggregation method is applied to combine multiple objectives into a single comprehensive objective.

F_{a g g} (x) = \sum_{k = 1}^{M} w_{k} f_{k} (x)

(23)

In this context,

M

denotes the total number of objectives and

w_{k}

represents the weight assigned to the

k_{t h}

objective. We assigned a weight of 0.7 to the PID and 0.3 to

λ_{b u c k l i n g}

\{\begin{matrix} F i n d X = (P_{1}, P_{2}, P_{3}, \dots, P_{16}) \\ M i n {P I D, \frac{1}{λ_{b u c k l i n g}}} \\ S . T . Count = ({C T}_{1}, {C T}_{2}, {C T}_{3}, {C T}_{4}) \end{matrix}

(24)

In this constraint,

P

represents the ply angles,

P I D

refers to processing-induced deformation,

λ_{b u c k l i n g}

denotes the maximum buckling eigenvalue, and

C T

indicates the statistical count of each ply type (0/±45/90). The schematic diagram of the optimization process is shown in Figure 9.

In CLPT, the B matrix within the ABD stiffness matrix characterizes the coupling between bending and tensile in composite laminates. This matrix plays a critical role in governing the spring-in deformation behavior. Consequently, this study employed the Frobenius norm of the B matrix as a quantitative metric to evaluate the degree of asymmetry in laminate stacking sequences, integrating it into the IAGA.

The constitutive matrix Q for each layer can be calculated using the uniaxial material stiffness of the composite, such as

E_{1}

E_{2}

G_{12}

, and

v_{12}

. For an isotropic material, the Q matrix has the following relationships:

Q = [\begin{matrix} \frac{E_{1}}{1 - v_{12} v_{21}} & \frac{{v_{12} E}_{2}}{1 - v_{12} v_{21}} & 0 \\ \frac{{v_{12} E}_{2}}{1 - v_{12} v_{21}} & \frac{E_{2}}{1 - v_{12} v_{21}} & 0 \\ 0 & 0 & G_{12} \end{matrix}]

(25)

where

v_{21}

is the other component of Poisson’s ratio, which is typically calculated as follows:

v_{21} = v_{12} \cdot \frac{E_{2}}{E_{1}}

(26)

In a laminate, the

B

matrix is the sum of the contributions of all layers to the total stiffness matrix. For each layer

k

, the contribution of the

B

matrix is as follows:

B_{k} = Q_{k} \cdot \frac{h_{k + 1}^{2} - h_{k}^{2}}{2}

(27)

where

h_{k}

and

h_{k + 1}

represent the positions of the top and bottom surfaces of two adjacent layers, respectively.

For a matrix

B

, its Frobenius norm

{‖B‖}_{F}

is defined as follows:

{‖B‖}_{F} = \sqrt{\sum_{i = 1}^{m} \sum_{j = 1}^{n} {|B_{i j}|}^{2}}

(28)

where

B = [B_{i j}]

is an

m \times n

matrix and

B_{i j}

is the element in the

i

-th row and

j

-th column of matrix

B

. The relationship between

{‖B‖}_{F}

and PID is further discussed in Appendix A.

3. Results and Discussion

3.1. Original NSGA-II Optimization Results

The original NSGA-II algorithm was used to perform stacking optimization for a 16-layer L-shaped laminate. After a genetic process with a population size of 12 (p = 12) and 20 generations (t = 20), the optimal individual (ply stack) along with the corresponding PID and

λ_{b u c k l i n g}

were obtained, taking the symmetric and balanced stacking sequence [45/90/−45/0]_2s as a reference. The genetic curve of the optimization process is shown in Figure 10, where the horizontal axis represents the number of iterations and the vertical axis represents the optimization objective.

The gray line in the figure represents the PID and

λ_{b u c k l i n g}

results after each iteration using the FEA model, while the red line represents the improved solutions selected from the best individual in each generation, with the selection method also based on non-dominated sorting. Reference stacking is a commonly used symmetric balanced stacking, with a curing deformation of 0.354° and a buckling eigenvalue of 42,321. The model converged after 175 iterations, with a final curing deformation of 0.28° and a corresponding buckling eigenvalue of approximately 47,000. From the evolution curve, it can be observed that most individuals in the population tend to converge to the same extremum, while multiple near-optimal solutions are eliminated, causing the algorithm to become trapped in a local optimum. This results in an imperfect evolutionary process, which is due to the limitations of the NSGA-II algorithm in handling high-dimensional nonlinear problems.

3.2. Optimization Results of IAGA Model

After 240 iterations, the current model parameters were able to identify only three Pareto points, as shown in Figure 11. This result is insufficient to meet the design requirements and does not achieve improvements in the PID. Consequently, the IAGA was developed based on the proposed strategies: adaptive non-dominated sorting, disturbance-driven crossover, and the dynamic Gaussian mutation operator. Using the same population size and genetic algorithm parameters, the results obtained are presented in Figure 12 and Figure 13.

By comparing the result curves of the two models, it is evident that the IAGA proposed in this study is significantly better than NSGA-II in handling ply optimization problems. With the introduction of adaptive non-dominated sorting, the IAGA adopts a more relaxed dominance relationship between solutions, enabling it to quickly explore high-quality ply configurations with a PID value of approximately 0.12° and

λ_{b u c k l i n g}

of 44,000 in the initial stages. Additionally, the inclusion of a perturbation mechanism allows for the exploration of the surrounding solution space, ensuring that the identified solutions remain globally optimal rather than locally optimal. This led to the discovery of a series of ply configurations with a PID value of approximately 0.125° and

λ_{b u c k l i n g}

of around 50,000. The improved dynamic Gaussian mutation operator further enhanced the algorithm’s ability to explore the solution space, resulting in a greater number of Pareto points. Compared to the PID value of 0.28°,

λ_{b u c k l i n g}

of 46,766, and the limited number of Pareto points obtained by the traditional NSGA-II, the IAGA demonstrates not only a faster convergence, but also significantly superior optimization results, the results obtained are presented in Figure 14.

3.3. IAGA Incorporates Asymmetry

In this study, the Frobenius norm of the B matrix was used as an indicator to measure the asymmetry of the laminate stacking sequence and was incorporated into the framework of the genetic algorithm for optimization. The optimization results are shown in Figure 15.

The results indicate that the introduction of the asymmetry indicator

{‖B‖}_{F}

had a positive impact on the optimization problem. Compared to the model without

{‖B‖}_{F}

in the optimization objectives, the optimization capability for PID improved by 20%, while the optimization capability for

λ_{b u c k l i n g}

remained consistent. However, the optimization efficiency slightly decreased, which is related to the overall complexity of the model. Figure 16 shows results of IAGA incorporates

{‖B‖}_{F}

. Figure 16 shows the result of NSGA-II, IAGA and IAGA-B. Figure 17 shows the Pareto front under the three optimization objectives.

Comparing the PID and λ calculation results before and after ply optimization against the symmetric balanced stacking allows for the verification of the optimization results, as shown in Figure 18, Figure 19 and Figure 20. The deformation results after ply optimization show some differences compared to the symmetric balanced deformation pattern. The optimized results exhibit a slight twisting deformation, but it is minimal.

The buckling mode also exhibits differences, with the buckling eigenvalue of the optimized ply significantly increased.

To highlight the excellent capabilities of the IAGA algorithm for ply optimization problems, optimization calculations were performed using the multi-island genetic algorithm (MIGA), neighborhood cultivation genetic algorithm (NCGA), and archive-based micro-genetic algorithm (AMGA). The parameters for these models were set according to the default recommended settings in Isight, and the results are shown in Figure 21.

3.4. Influence of Population Size

In this study, while using the improved NSGA-II genetic algorithm for composite ply optimization design, it was observed that the population size is a critical factor influencing the algorithm’s convergence and accuracy. Since the population size determines the number of candidate solutions in each generation, it directly affects the exploration capability of the search space and the quality of the final solution. A smaller population size reduces the computational cost, but is prone to becoming trapped in local optima, leading to a decreased robustness of the algorithm. Increasing the population size enhances the diversity, allows for broader coverage of the solution space, and provides more candidate solutions to explore the search space, thereby increasing the likelihood of finding the global optimal solution. However, a larger population size requires more fitness evaluations and selection operations in each generation, resulting in higher computational costs and the need for more computational resources and time. The comparison of computational efficiency between the two models with a population size of 64 is shown in Figure 22.

3.5. Computation Time

In industrial applications, feasible optimization solutions are inevitably constrained by their computation time and CPU requirements. Therefore, it is relevant to briefly mention the model size and computational demands. All work in this study was conducted on a workstation with Intel 12,700 K (16 logical processors) and 32 GB of RAM. The model consists of 3900 elements, with 10 CPU cores dedicated to calculating the PID and the remaining cores used for computing λ, both running concurrently. The runtime for a single iteration, including model INP file generation, solver execution, and result data transfer, was 4 min. A comparison of the computation times between the NSGA-II and IAGA models is shown in Figure 23. The computation time here refers to the total time elapsed until convergence.

The PID optimization capability refers to the percentage decrease in the optimized values of the spring-in angle and λ compared to the results of the reference stacking. A cross-comparison of the results shown in Table 3 from the NSGA-II, IAGA, asymmetry-enhanced IAGA-B, MIGA, NCGA, and AMGA revealed that the IAGA significantly reduced the computation time while markedly improving the optimization performance. After introducing asymmetry as an optimization objective, the convergence speed slightly decreased; however, the exploration efficiency of the model improved. This indicates that asymmetry plays an important role in PID optimization.

The improved IAGA model significantly enhanced the optimization efficiency and result quality, effectively reducing processing-induced deformations (PIDs) and improving the buckling eigenvalues. Additionally, its flexible adaptability and efficient convergence provide new insights into addressing complex multi-objective optimization problems. For laminated plates with other geometric features, based on the previous analysis, the mechanism of spring-in remains the same. Therefore, this method is equally applicable to other complex structures with corner features. However, the method relies on high-fidelity finite element analysis models, presenting limitations in terms of computational resource demands. Future research could explore the application of surrogate models, extend the approach to more complex geometries and loading conditions, and incorporate additional manufacturing constraints to enhance its practicality.

4. Conclusions

This study focused on multi-objective ply optimization for the PIDs and buckling eigenvalues of composite laminated structures. A high-fidelity FEA model was established, and an optimization framework was developed to address the PID of L-shaped stiffeners. Starting from the concept of asymmetric ply design, this research proposes several improvements to the traditional NSGA-II genetic algorithm, including adaptive non-dominated sorting, disturbance-driven crossover, and a dynamic Gaussian mutation operator. These proposals were derived from an in-depth analysis of the physical mechanisms underlying the curing process. The improved IAGA model significantly enhances the optimization efficiency compared to the traditional model, while also further strengthening its multi-objective optimization capabilities. It mitigates the shortcomings of the GA in handling high-dimensional nonlinear problems, such as becoming trapped in local optima, and successfully establishes an efficient and accurate improved genetic algorithm.

The improved IAGA model significantly enhances the optimization efficiency, reducing the curing deformation (PID) by 50% compared to symmetric balanced stacking, and improving the λ by 20%. Compared to the NSGA-II algorithm, the IAGA achieved a 3.5-fold improvement in convergence efficiency while delivering quantitatively superior optimization outcomes, such as identifying configurations with a PID of 0.12° and a buckling eigenvalue of over 50,000. Under the same number of iterations and constraints, the IAGA model successfully identified a Pareto solution set with a capacity of 10, representing a threefold improvement in exploration capability compared to the traditional model, which could only generate approximately three Pareto solutions. The incorporation of asymmetry as an optimization objective further enhanced the optimization capability by 20% for PIDs while maintaining the structural stability. These results validate the feasibility and advantage of applying non-symmetric ply designs in composite laminates and provide valuable insights for practical engineering applications.

Our ongoing research focuses on utilizing surrogate models to replace the FEA for PID computations, with the aim of significantly enhancing the optimization efficiency and enabling rapid ply optimization under constrained conditions. Concurrently, the feasibility of extending this methodology to accommodate more complex geometries and stacking configurations is being systematically investigated.

Author Contributions

Conceptualization, Q.L. and X.W.; methodology, Q.L.; software, X.W. and L.Y.; validation, Q.L., X.W. and Z.G.; formal analysis, Z.L.; investigation, Q.L. and L.Y.; resources, X.W.; data curation, X.W.; writing—original draft preparation, Q.L.; writing—review and editing, X.W.; visualization, Z.L. and L.Y.; supervision, Z.G.; project administration, Z.G.; funding acquisition, X.W. All authors have read and agreed to the published version of the manuscript.

Funding

This research was funded by the National Key Research and Development Program of China, grant number 2021YFB3401700; the Young Elite Scientists Sponsorship Program by CAST, grant number 2021QNRC001; and the Aeronautical Science Foundation of China, grant number 2022Z039051001.

Institutional Review Board Statement

Not applicable.

Informed Consent Statement

Not applicable.

Data Availability Statement

The raw data supporting the conclusions of this article will be made available by the authors upon request.

Conflicts of Interest

The authors declare no conflicts of interest.

Nomenclature

AMGA	archive-based micro-genetic algorithm
CLPT	classical laminate plate theory
CTE	coefficient of thermal expansion
DD	“double–double” layup
DoC	degree of cure
FEA	finite element analysis
GA	genetic algorithm
IAGA	improved adaptive genetic algorithm
MIGA	multi-island genetic algorithm
MRCC	manufacturer-recommended cure cycle
NCGA	neighborhood cultivation genetic algorithm
NSGA	non-dominated sorting genetic algorithm
PID	processing-induced deformation
RVE	representative volume element
$ρ$	density
$c$	specific heat
$T$	temperature
$t$	time
$\dot{q}$	internal heat generation rate per unit volume
$K_{x}$	thermal conductivity coefficients of x-direction
$K_{y}$	thermal conductivity coefficients of y-direction
$K_{z}$	thermal conductivity coefficients of z-direction
$T_{g}$	current glass transition temperature
$T_{g, 0}$	initial glass transition temperature (uncured state)
$T_{g, \infty}$	$glass transition temperature at full cure (α = 1$ )
$T_{c}$	characteristic temperature
$λ$	parameter of curing conditions
$α$	degree of cure
$σ_{i j}$	total stress
${\underline{σ}}_{i j}^{f}$	effective stress in the fiber
$δ_{i j}$	Kronecker delta function
$P_{r}$	resin pressure
$V_{0}$	fiber volume fraction when the load is zero
$V_{f}$	fiber volume fraction
$μ$	resin viscosity
$K_{i j}$	permeability
$r_{f}$	fiber radius
$K_{0}$	Kozeny constant
$V o l$	total volume
$V o l_{f}$	volume occupied by the fibers
$E_{m}^{0}$	rubbery modulus
$E_{m}^{1}$	glassy modulus
$[Q]^{k}$	$reduced stiffness of the k$ -th layer
$h_{k}$	positions of ply
$N$	total number of layers
$D_{i j}$	elements of the stiffness matrix of the laminated plate
$N_{x}$	in-plane stresses of x-direction
$N_{y}$	in-plane stresses of y-direction
$N_{x y}$	in-plane stresses of z-direction
$γ$	rate of increase in tightness
$f_{k} (x)$	value of the objective function
$η$	perturbation intensity
$N (0, σ^{2})$	standard normal distribution noise
$ξ$	crossover coefficient
$σ_{k} (t)$	mutation amplitude
$M$	total number of objectives
$λ_{b u c k l i n g}$	buckling eigenvalue
$C T$	statistical count
${‖B‖}_{F}$	Frobenius norm of matrix B

Appendix A

To demonstrate the definitive influence of the

{‖B‖}_{F}

on the PID in L-shaped structures, the PID corresponding to the

{‖B‖}_{F}

of generated stacking sequences was calculated. The stacking sequence list is shown in Table A1, and the PID results are displayed in Figure A1.

Table A1. Ply sequence.

Num	Ply Sequence	${‖B‖}_{F}$	PID
P1	[45/90/-45/-45/90/-45/-45/-45/0/-45/-45/-45/0/-45/45/90]	71.17	0.414
P2	[-45/-45/45/90/-45/90/45/90/45/45/45/90/-45/90/-45/0]	33.19	0.131
P3	[0/45/45/-45/45/-45/45/90/90/45/45/90/-45/-45/45/45]	37.86	0.286
P4	[0/45/-45/90/0/-45/45/45/45/90/-45/-45/0/45/90/-45]	35.34	0.133
P5	[90/0/-45/45/-45/45/45/0/-45/45/-45/45/-45/-45/0/90]	60.64	0.369
P6	[90/45/45/0/45/45/45/90/45/0/45/0/90/45/45/-45]	27.41	0.124

Figure A1. PID corresponding to

{‖B‖}_{F}

Figure A1. PID corresponding to

{‖B‖}_{F}

References

Wang, B.; Fan, S.; Chen, J.; Yang, W.; Liu, W.; Li, Y. A Review on Prediction and Control of Curing Process-Induced Deformation of Continuous Fiber-Reinforced Thermosetting Composite Structures. Compos. Part A Appl. Sci. Manuf. 2023, 165, 107321. [Google Scholar] [CrossRef]
Amann, C.; Kreissl, S.; Grass, H.; Meinhardt, J. A Review on Process-Induced Distortions of Carbon Fiber Reinforced Thermosets for Large-Scale Production. Prod. Eng. 2017, 11, 665–675. [Google Scholar] [CrossRef]
Ding, A.; Li, S.; Sun, J.; Wang, J.; Zu, L. A Comparison of Process-Induced Residual Stresses and Distortions in Composite Structures with Different Constitutive Laws. J. Reinf. Plast. Compos. 2016, 35, 807–823. [Google Scholar] [CrossRef]
Shokrieh, M.M.; Safarabadi, M. Understanding Residual Stresses in Polymer Matrix Composites; Elsevier Ltd.: Amsterdam, The Nethterlands, 2021; ISBN 978-0-12-818817-0. [Google Scholar]
Goncalves, P.T.; Arteiro, A.; Rocha, N. Micro-Mechanical Analysis of the Effect of Ply Thickness on Curing Micro-Residual Stresses in a Carbon/Epoxy Composite Laminate. Compos. Struct. 2023, 319, 117158. [Google Scholar] [CrossRef]
Callister, W.D. Materials Science and Engineering: An Introduction; Wiley: New York, NY, USA, 1994; ISBN 978-0-471-58128-4. [Google Scholar]
Ding, A.; Li, S.; Wang, J.; Zu, L. A Three-Dimensional Thermo-Viscoelastic Analysis of Process-Induced Residual Stress in Composite Laminates. Compos. Struct. 2015, 129, 60–69. [Google Scholar] [CrossRef]
Abouhamzeh, M.; Sinke, J.; Benedictus, R. Prediction Models for Distortions and Residual Stresses in Thermoset Polymer Laminates: An Overview. J. Manuf. Mater. Process. 2019, 3, 87. [Google Scholar] [CrossRef]
Magnus Svanberg, J.; Anders Holmberg, J. Prediction of Shape Distortions Part I. FE-Implementation of a Path Dependent Constitutive Model. Compos. Part A Appl. Sci. Manuf. 2004, 35, 711–721. [Google Scholar] [CrossRef]
Svanberg, J.M.; Holmberg, J.A. Prediction of Shape Distortions. Part II. Experimental Validation and Analysis of Boundary Conditions. Compos. Part A Appl. Sci. Manuf. 2004, 35, 723–734. [Google Scholar] [CrossRef]
Liu, C.; Shi, Y. A Thermo-Viscoelastic Analytical Model for Residual Stresses and Spring-in Angles of Multilayered Thin-Walled Curved Composite Parts. Thin-Walled Struct. 2020, 152, 106758. [Google Scholar] [CrossRef]
Hou, J.; You, B.; Xu, J.; Fu, T.; Wang, T. Numerical Simulation and Multi-Objective Optimization for Curing Process of Thermosetting Prepreg. Appl. Compos. Mater. 2022, 29, 1409–1429. [Google Scholar] [CrossRef]
Khairul Alam, M.; Anghelescu, M.S. Analysis of Deformation and Residual Stresses in Composites Processed on a Carbon Foam Tooling. J. Compos. Mater. 2009, 43, 2057–2070. [Google Scholar] [CrossRef]
Cameron, C.; Hozić, D.; Stig, F.; van der Veen, S. A Method for Optimization against Cure-Induced Distortion in Composite Parts. Struct. Multidisc. Optim. 2023, 66, 51. [Google Scholar] [CrossRef]
Chwal, M.; Muc, A.; Stawiarski, A.; Barski, M. Residual stresses in multilayered composites—General overview. Compos. Theory Pract. 2016, 16, 132–138. [Google Scholar]
Li, J.; Li, D. Extension-Shear Coupled Laminates with Immunity to Hygro-Thermal Shearing Distortion. Compos. Struct. 2015, 123, 401–407. [Google Scholar] [CrossRef]
Vermes, B.; Tsai, S.W.; Massard, T.; Springer, G.S.; Czigany, T. Design of Laminates by a Novel “Double–Double” Layup. Thin-Walled Struct. 2021, 165, 107954. [Google Scholar] [CrossRef]
Pelletier, J.L.; Vel, S.S. Multi-Objective Optimization of Fiber Reinforced Composite Laminates for Strength, Stiffness and Minimal Mass. Comput. Struct. 2006, 84, 2065–2080. [Google Scholar] [CrossRef]
Matsuzaki, R.; Yokoyama, R.; Kobara, T.; Tachikawa, T. Multi-Objective Curing Optimization of Carbon Fiber Composite Materials Using Data Assimilation and Localized Heating. Compos. Part A Appl. Sci. Manuf. 2019, 119, 61–72. [Google Scholar] [CrossRef]
Haftka, R.T. Genetic Algorithms for Optimization of Composite Laminates. In Mechanics of Composite Materials and Structure; Soares, C.A.M., Soares, C.M.M., Freitas, M.J.M., Eds.; Springer: Dordrecht, The Nethterlands, 1999; pp. 431–442. ISBN 978-94-011-4489-6. [Google Scholar]
Michalewicz, Z. Genetic Algorithms + Data Structures = Evolution Programs; Springer: Berlin/Heidelberg, Germany, 1996; ISBN 978-3-642-08233-7. [Google Scholar]
Riche, R.L.; Haftka, R.T. Optimization of Laminate Stacking Sequence for Buckling Load Maximization by Genetic Algorithm. AIAA J. 1993, 31, 951–956. [Google Scholar] [CrossRef]
Choudhary, P.K. Optimal Design of Composite Cylindrical Shells Subject to Compression Buckling Strength. Multidiscip. Model. Mater. Struct. 2023, 19, 829–847. [Google Scholar] [CrossRef]
Huynh, T.N.; Lee, J. Optimal Thickness Distribution Design for Blending Composite Laminates Using Buckling Factor Prediction. Compos. Struct. 2024, 327, 117693. [Google Scholar] [CrossRef]
Zang, W.; Ren, L.; Zhang, W.; Liu, X. A Cloud Model Based DNA Genetic Algorithm for Numerical Optimization Problems. Future Gener. Comput. Syst. 2018, 81, 465–477. [Google Scholar] [CrossRef]
Srinivas, N.; Deb, K. Muiltiobjective Optimization Using Nondominated Sorting in Genetic Algorithms. Evol. Comput. 1994, 2, 221–248. [Google Scholar] [CrossRef]
Deb, K.; Pratap, A.; Agarwal, S.; Meyarivan, T. A Fast and Elitist Multiobjective Genetic Algorithm: NSGA-II. IEEE Trans. Evol. Computat. 2002, 6, 182–197. [Google Scholar] [CrossRef]
Luo, B.; Zheng, J.; Xie, J.; Wu, J. Dynamic Crowding Distance? A New Diversity Maintenance Strategy for MOEAs. In Proceedings of the 2008 Fourth International Conference on Natural Computation, Jinan, China, 18–20 October 2008; Volume 1, pp. 580–585. [Google Scholar]
Coello, C.A.; Toscano Pulido, G. A Micro-Genetic Algorithm for Multiobjective Optimization. In Proceedings of the Evolutionary Multi-Criterion Optimization; Zitzler, E., Thiele, L., Deb, K., Coello Coello, C.A., Corne, D., Eds.; Springer: Berlin/Heidelberg, Germany, 2001; pp. 126–140. [Google Scholar]
Emmerich, M.T.M.; Deutz, A.H. A Tutorial on Multiobjective Optimization: Fundamentals and Evolutionary Methods. Nat. Comput. 2018, 17, 585–609. [Google Scholar] [CrossRef] [PubMed]
Yuan, Z.; Kong, L.; Gao, D.; Tong, X.; Feng, Y.; Yang, G.; Yang, Z.; Li, S. Multi-Objective Approach to Optimize Cure Process for Thick Composite Based on Multi-Field Coupled Model with RBF Surrogate Model. Compos. Commun. 2021, 24, 100671. [Google Scholar] [CrossRef]
Yuan, Z.; Tong, X.; Yang, G.; Yang, Z.; Song, D.; Li, S.; Li, Y. Curing Cycle Optimization for Thick Composite Laminates Using the Multi-Physics Coupling Model. Appl. Compos. Mater. 2020, 27, 839–860. [Google Scholar] [CrossRef]
Gao, Y.; Ye, J.; Yuan, Z.; Ling, Z.; Zhou, Y.; Lin, Z.; Dong, J.; Wang, H.; Peng, H.X. Optimization Strategy for Curing Ultra-Thick Composite Laminates Based on Multi-Objective Genetic Algorithm. Compos. Commun. 2022, 31, 101115. [Google Scholar] [CrossRef]
Fan, H.-T.; Wang, H.; Chen, X.-H. An Optimization Method for Composite Structures with Ply-Drops. Compos. Struct. 2016, 136, 650–661. [Google Scholar] [CrossRef]
Kumpati, R.; Skarka, W.; Skarka, M.; Brojan, M. Enhanced Optimization of Composite Laminates: Multi-Objective Genetic Algorithms with Improved Ply-Stacking Sequences. Materials 2024, 17, 887. [Google Scholar] [CrossRef]
Yan, X. Finite Element Modeling of Curing of Epoxy Matrix Composites. J. Appl. Polymer. Sci. 2007, 103, 2310–2319. [Google Scholar] [CrossRef]
Wen, Y. Study on the curing reaction kinetics of modified bismaleimide resin. China Plast. Ind. 2018, 46, 22–26+32. [Google Scholar]
Liu, X.; Guan, Z.; Wang, X.; Jiang, T.; Geng, K.; Li, Z. Study on Cure-Induced Residual Stresses and Spring-in Deformation of L-Shaped Composite Laminates Using a Simplified Constitutive Model Considering Stress Relaxation. Compos. Struct. 2021, 272, 114203. [Google Scholar] [CrossRef]
Minakuchi, S.; Niwa, S.; Tagakaki, K.; Takeda, N. Composite Cure Simulation Scheme Fully Integrating Internal Strain Measurement. Compos. Part A Appl. Sci. Manuf. 2016, 84, 53–63. [Google Scholar] [CrossRef]
Khoun, L.; Centea, T.; Hubert, P. Characterization Methodology of Thermoset Resins for the Processing of Composite Materials—Case Study: CYCOM 890RTM Epoxy Resin. J. Compos. Mater. 2010, 44, 1397–1415. [Google Scholar] [CrossRef]
Yue, G.; Zhang, B.; Dai, F.; Du, S. Three-Dimensional Cure Simulation of Stiffened Thermosetting Composite Panels. J. Mater. Sci. Technol. 2010, 26, 467–471. [Google Scholar] [CrossRef]
Zeng, X.; Raghavan, J. Role of Tool-Part Interaction in Process-Induced Warpage of Autoclave-Manufactured Composite Structures. Compos. Part A Appl. Sci. Manuf. 2010, 41, 1174–1183. [Google Scholar] [CrossRef]
Tian, M.; Gong, X.; Yin, L.; Li, H.; Ming, W.; Zhang, Z.; Chen, J. Multi-Objective Optimization of Injection Molding Process Parameters in Two Stages for Multiple Quality Characteristics and Energy Efficiency Using Taguchi Method and NSGA-II. Int. J. Adv. Manuf. Technol. 2017, 89, 241–254. [Google Scholar] [CrossRef]
Rong, Y.; Zhou, Q.; Huang, Y.; Chang, Y.; Zhang, G.; Shao, X. Multi-Objective Optimization of Laser Brazing with the Crimping Joint Using ANN and NSGA-II. Int. J. Adv. Manuf. Technol. 2016, 85, 1239–1247. [Google Scholar] [CrossRef]
Optimal Design of Laminated Composite Structures with Ply Drops Using Stacking Sequence Tables. Compos. Struct. 2014, 107, 559–569. [CrossRef]

Figure 1. Resin flow-compaction procedure.

Figure 2. Deformation differences in ply layers with different angles.

Figure 3. Dimensions of FEA model.

Figure 4. Curing process curve and degree of cure.

Figure 5. FEA model of PIDs.

Figure 6. FEA model of

λ_{b u c k l i n g}

Figure 6. FEA model of

λ_{b u c k l i n g}

Figure 7. Standard NSGA-II process.

Figure 8. IAGA model.

Figure 9. Ply optimization with FEA and IAGA.

Figure 10. Original NSGA-II optimization curve of (a) PID and (b)

λ_{b u c k l i n g}

Figure 10. Original NSGA-II optimization curve of (a) PID and (b)

λ_{b u c k l i n g}

Figure 11. Pareto points of original NSGA-II optimization.

Figure 12. Improved model optimization curve of (a) PID and (b)

λ_{b u c k l i n g}

Figure 12. Improved model optimization curve of (a) PID and (b)

λ_{b u c k l i n g}

Figure 13. Pareto points of improved model optimization.

Figure 14. Comparison of optimization capabilities between IAGA and NSGA-II. (a) PID; (b)

λ_{b u c k l i n g}

Figure 14. Comparison of optimization capabilities between IAGA and NSGA-II. (a) PID; (b)

λ_{b u c k l i n g}

Figure 15. IAGA incorporates

{‖B‖}_{F}

. (a) PID; (b)

λ_{b u c k l i n g}

Figure 15. IAGA incorporates

{‖B‖}_{F}

. (a) PID; (b)

λ_{b u c k l i n g}

Figure 16. Results comparison of NSGA-II, IAGA and IAGA incorporates

{‖B‖}_{F}

. (a) PID; (b)

λ_{b u c k l i n g}

Figure 16. Results comparison of NSGA-II, IAGA and IAGA incorporates

{‖B‖}_{F}

. (a) PID; (b)

λ_{b u c k l i n g}

Figure 17. Three-dimensional Pareto front.

Figure 18. Comparison of PID before and after optimization (whole model). (a) Before optimization; (b) After optimization.

Figure 19. Comparison of PID before and after optimization (L-shape). (a) Before optimization; (b) After optimization.

Figure 20. Comparison of λ before and after optimization. (a) Before optimization; (b) After optimization.

Figure 21. Comparison of optimization capabilities between different GA models. (a) PID; (b)

λ_{b u c k l i n g}

Figure 21. Comparison of optimization capabilities between different GA models. (a) PID; (b)

λ_{b u c k l i n g}

Figure 22. Population size of 64. (a) PID; (b)

λ_{b u c k l i n g}

Figure 22. Population size of 64. (a) PID; (b)

λ_{b u c k l i n g}

Figure 23. Comparison of computation times.

Table 1. Curing kinetics parameters and total heat of reaction for QY9611 resin.

Material	A/(s⁻¹)	∆E/(J/mol)	n	H_u/(J/kg)	$T_{g, 0}$ /(K)	$T_{g, \infty}$ /(K)	$λ$
QY9611	2.28 × 10⁷	92410	1.485	3.0032 × 10⁵	264.04	513.75	0.4824

Table 2. Properties of ZT7H/QY9611.

Property	Rubbery	Glassy
$E_{1} / G P a$	140	144
$E_{2} = E_{3 f} / G P a$	0.16	10.2
$v_{12} = v_{13 f}$	0.3	0.3
$v_{23}$	0.6	0.4
$G_{12} = G_{13 f} / G P a$	0.05	6.0
$G_{23} / G P a$	0.04	3.0
$α_{1} / (μ ε \cdot K^{- 1})$	0.2
$α_{2} = α_{3} / (μ ε \cdot K^{- 1})$	40.9
$β_{1} / μ ε$	−167
$β_{2} = β_{3} / μ ε$	−8810

Table 3. Comparison of optimization capability.

Model	Iterations to Convergence	PID Optimization Capability	λ Optimization Capability	Computation Time
NSGA-II	175	20.90%	10.50%	11.7 h
IAGA	50	64.69%	17.35%	3.3 h
IAGA-B	80	68.93%	14.80%	5.3 h
MIGA	144	2.54%	8.75%	9.6 h
NCGA	196	−55.65%	9.64%	13.1 h
AMGA	205	−26.55%	11.40%	13.7 h

Disclaimer/Publisher’s Note: The statements, opinions and data contained in all publications are solely those of the individual author(s) and contributor(s) and not of MDPI and/or the editor(s). MDPI and/or the editor(s) disclaim responsibility for any injury to people or property resulting from any ideas, methods, instructions or products referred to in the content.

© 2025 by the authors. Licensee MDPI, Basel, Switzerland. This article is an open access article distributed under the terms and conditions of the Creative Commons Attribution (CC BY) license (https://creativecommons.org/licenses/by/4.0/).

Share and Cite

MDPI and ACS Style

Liu, Q.; Wang, X.; Guan, Z.; Li, Z.; Yang, L. Ply Optimization of Composite Laminates for Processing-Induced Deformation and Buckling Eigenvalues Based on Improved Genetic Algorithm. Materials 2025, 18, 345. https://doi.org/10.3390/ma18020345

AMA Style

Liu Q, Wang X, Guan Z, Li Z, Yang L. Ply Optimization of Composite Laminates for Processing-Induced Deformation and Buckling Eigenvalues Based on Improved Genetic Algorithm. Materials. 2025; 18(2):345. https://doi.org/10.3390/ma18020345

Chicago/Turabian Style

Liu, Qingchuan, Xiaodong Wang, Zhidong Guan, Zengshan Li, and Lingxiao Yang. 2025. "Ply Optimization of Composite Laminates for Processing-Induced Deformation and Buckling Eigenvalues Based on Improved Genetic Algorithm" Materials 18, no. 2: 345. https://doi.org/10.3390/ma18020345

APA Style

Liu, Q., Wang, X., Guan, Z., Li, Z., & Yang, L. (2025). Ply Optimization of Composite Laminates for Processing-Induced Deformation and Buckling Eigenvalues Based on Improved Genetic Algorithm. Materials, 18(2), 345. https://doi.org/10.3390/ma18020345

Note that from the first issue of 2016, this journal uses article numbers instead of page numbers. See further details here.

Article Menu

Ply Optimization of Composite Laminates for Processing-Induced Deformation and Buckling Eigenvalues Based on Improved Genetic Algorithm

Abstract

1. Introduction

2. Theory and Methods

2.1. Curing Process Theoretical Analysis

2.1.1. Thermo-Chemical Models

2.1.2. Resin Flow-Compaction Process

2.1.3. Stress-Deformation Process

2.1.4. Calculation of Buckling Eigenvalues

2.1.5. FEA Model and Mechanical Parameters

2.2. Improved Genetic Algorithm Based on NSGA-II

2.2.1. Limitations of NSGA-II in Ply Optimization

2.2.2. Adaptive Non-Dominated Sorting

2.2.3. Disturbance-Driven Crossover

2.2.4. Dynamic Gaussian Mutation Operator

2.2.5. Optimization Objectives and Constraints

3. Results and Discussion

3.1. Original NSGA-II Optimization Results

3.2. Optimization Results of IAGA Model

3.3. IAGA Incorporates Asymmetry

3.4. Influence of Population Size

3.5. Computation Time

4. Conclusions

Author Contributions

Funding

Institutional Review Board Statement

Informed Consent Statement

Data Availability Statement

Conflicts of Interest

Nomenclature

Appendix A

References

Share and Cite

Article Metrics

Article Access Statistics

Further Information

Guidelines

MDPI Initiatives

Follow MDPI