Optimal Design of Laminates Orientation in a Composite Material by Genetics Algorithm and Simulated Annealing

1. Introduction

For the design of a composite material laminate, it is essential to consider several factors to satisfy the engineering and financial requirements. The first decision drivers are the raw materials with which the composite will be manufactured. Reinforcement and matrix will determine together the end material's physical, chemical, and mechanical properties. However, the number of plies, orientation, and stacking sequence will define the performance of the final structure to resist the applied loads [1].

One of the unique characteristics of composites is that mechanical properties can be tailored to resist specific load conditions. The fibers are stronger in their longitudinal axis. Therefore, the design aims to orient as many layers as possible in the direction of the applied stress. As layers are stacked through the thickness at different orientations, hundreds of combinations can exist that could fulfill the mechanical purpose [2,3].

In the case of uniaxial or biaxial stress, selecting plies' orientation can be easy, which can be accomplished by experimentation or computer simulation. However, as the stress state becomes triaxial due to combined loads and moments, it becomes a real challenge requiring a more exhaustive analysis of the orientation and stack-up sequence of the constituent plies. The common practice on composite structures design suggests the implementation of quasi-isotropic configuration, [0/90/+-45], with [50-10-40] volume proportions, respectively. The goal is to have an isotropic and balanced material which could fulfill the most common load cases but is not optimized.

The optimization of composite laminates contemplates two main goals: 1) minimization of weight [4,5] and 2) maximization of the load-carrying capacity [1,6]. The optimization of the specific strength is translated into a reduction of raw materials, manufacturing tasks, and laps, and, therefore, cost improvement [7].

These two goals have motivated the use of several optimization methods (gradient-based, direct search, heuristics, and hybrid), turning the composite's structural design into a more systematic and well-defined task, less dependent on the designer and more focused on the 5material performance [6,8,9].

For composite laminate structures, two stages need to be accomplished. First, in constant stiffness consideration, the composite part is treated as a single element, and the goal is to find an optimal stacking sequence uniform for the entire structure. Second, in variable stiffness consideration, the composite part is partitioned into several elements, each with a particular fiber orientation and stacking sequence [10,11].

For the first stage, direct search methods have been employed due to their advantage of not requiring information on the gradient and for not being sensitive to the initial search point; founding their results by using function values of preceding iterations. These methods strongly portray non-linear and noisy search spaces [12,13,14,15]. Also, the methods cover more uniformly ample search space; due to their inherent stochasticity, it escapes more easily from local optima. In this category, two powerful tools are found, genetic algorithms (GA) [16,17,18,19,20,21] and simulated annealing (SA) [16,22,23,24,25,26]. Newer and more sophisticated optimization algorithms are reported elsewhere [27,28]. However, our goal is to present a rigorous procedure based on strain and safety criteria to support the design of laminated material.

The laminates' orientation affects the composite material's mechanical behavior in several ways. Laminates oriented in the same direction as the applied load will experience more significant stiffness and strength in that direction. In contrast, laminates oriented perpendicular to the applied load will experience higher flexibility and toughness in that direction. Since the optimal orientation of laminates in a composite material depends on the specific application and the desired mechanical properties, this work aims to evaluate the optimization methodologies based on Genetic Algorithm and Simulated Annealing to support the design of fiber-reinforced laminates and define rigorous criteria to decide if the standard practices' orientations guarantee at the same time the maximum strain and Tsai-Wu safety factor criteria.

This paper presents a numerical investigation to evaluate two different optimization algorithms applied to composite structures. The model description takes place in the first section by defining the corresponding physical and mathematical considerations. Then, in the second section, the problem formulation is set to calculate the optimal orientations design in a stacking sequence which fulfills at the same time both the maximum strain criterion and the Tsai-Wu failure criterion. The third section shows the results for three different study cases. Finally, conclusions are cited regarding the performance of the Genetic Algorithm (GA) and Simulated Annealing (SA) as optimization techniques for calculating fiber-reinforced composites loaded under consideration of classic laminate theory.

2. Methodology

2.1. Model description

The optimization is carried out by assuming considerations that allow formulating the corresponding mathematical model which describes the physical process. Such suppositions are listed below:

• All laminas have an orthotropic behavior, which implies that the stiffness matrix is reduced to a 6 X 6 matrix system with 21 independent elastic constants.
• The laminate thickness is small compared to the other dimensions, implying that all laminas are restricted to an in-plane stress state, so stresses in the z-direction are negligible. This means that the stiffness matrix becomes a 3 X 3 matrix system with 4 independent elastic constants.
• There is a perfect bond between different laminas. Here it is assumed that all of them are subjected to the same midplane strain.
• Strain distribution in the laminate thickness is linear so that the response of each lamina depends on its relative location regarding the midplane of the laminate.
• All laminas are macroscopically homogeneous and have a linear-elastic behavior. This means that mechanical properties do not depend on the spatial location.
• The laminate works under a constant stiffness principle. This means that the stiffness is the same in the entire domain of the laminate and no discontinuity of each ply is allowed.
• Residual stresses derived from the curing process are negligible.

2.2. Theory of laminates

This classical theory is based on the assumptions deployed above for a thin-shell structure so that the mechanical properties of a composite material are evaluated by considering in-plane stress conditions. The generalized Hooke's equation for an orthotropic, three-dimensional composite material gives a linear relationship between stress and strain as follows:

$$\left(\begin{array}{c}{\sigma }_x\\ {\sigma }_y\\ {\sigma }_z\\ {\tau }_{xz}\\ {\tau }_{yz}\\ {\tau }_{xy}\end{array}\right)=\left(\begin{array}{cccccc}{C}_{11}& C_{12}& C_{13}& C_{14}& C_{15}& C_{16}\\ C_{21}& C_{22}& C_{23}& C_{24}& C_{25}& C_{26}\\ C_{31}& C_{32}& C_{33}& C_{34}& C_{35}& C_{36}\\ C_{41}& C_{42}& C_{43}& C_{44}& C_{45}& C_{46}\\ C_{51}& C_{52}& C_{53}& C_{54}& C_{55}& C_{56}\\ C_{61}& C_{62}& C_{63}& C_{64}& C_{65}& C_{66}\end{array}\right)\left(\begin{array}{c}{\epsilon }_x\\ {\epsilon }_y\\ {\epsilon }_z\\ {\gamma }_{xz}\\ {\gamma }_{yz}\\ {\gamma }_{xy}\end{array}\right)$$

(1)

and, according to the assumptions defined above, the equation that describe the behavior of a lamina in an orthotropic composite material is given by:

$$\left(\begin{array}{c}{\sigma }_x\\ {\sigma }_y\\ {\tau }_{xy}\end{array}\right)=\left(\begin{array}{ccc}{C}_{11}& C_{12}& 0\\ C_{12}& C_{22}& 0\\ 0& 0& C_{66}\end{array}\right)\left(\begin{array}{c}{\epsilon }_x\\ {\epsilon }_y\\ {\gamma }_{xy}\end{array}\right)$$

(2)

Where $C_{ij}$ elements constitute the material's stiffness matrix and depend on the material properties (E₁, E₂, G₁₂, and ν₁₂) in the principal directions.

A laminate composed of $k$ laminas of orthotropic composite material (Figure 1) is expected to be loaded with in-plane forces and moments that will produce strains and curvatures, respectively [29,30]. The resultant of forces and moments acting on the laminate can be obtained by integration of stresses along the thickness (z) according to the equations (Eqs. (3) and (4), respectively.

$$\left[N_x{N}_y{N}_{xy}\right]={\displaystyle {\int }_{-t/2}^{t/2}\left[{\sigma }_x{\sigma }_y{\sigma }_{xy}\right]}dz={{\displaystyle \sum _{k=1}^N{\displaystyle {\int }_{z_{k-1}}^{z_k}\left[{\sigma }_x{\sigma }_y{\sigma }_{xy}\right]}}}_{k}dz$$

(3)

$$\left[M_x{M}_y{M}_{xy}\right]={\displaystyle {\int }_{-t/2}^{t/2}\left[{\sigma }_x{\sigma }_y{\sigma }_{xy}\right]}zdz={{\displaystyle \sum _{k=1}^N{\displaystyle {\int }_{z_{k-1}}^{z_k}\left[{\sigma }_x{\sigma }_y{\sigma }_{xy}\right]}}}_{k}zdz$$

(4)

Note that $N_i$ (i = x, y, xy) are the forces and $M_i$ the moments, both, per unit length. These terms are related to the strains and curvatures by means of the stiffness of the laminate.

The representation of such properties in terms of the ABD matrix is:

$$\left(\begin{array}{c}{N}_x\\ N_y\\ N_{xy}\\ M_x\\ M_y\\ M_{xy}\end{array}\right)=\left(\begin{array}{cccccc}{A}_{11}& A_{12}& A_{16}& B_{11}& B_{12}& B_{16}\\ A_{12}& A_{22}& A_{26}& B_{12}& B_{22}& B_{26}\\ A_{16}& A_{26}& A_{66}& B_{16}& B_{26}& B_{66}\\ B_{11}& B_{12}& B_{16}& D_{11}& D_{12}& D_{16}\\ B_{12}& B_{22}& B_{26}& D_{12}& D_{22}& D_{26}\\ B_{16}& B_{26}& B_{66}& D_{16}& D_{26}& D_{66}\end{array}\right)\left(\begin{array}{c}{\epsilon }_x^0\\ {\epsilon }_y^0\\ {\gamma }_{xy}\\ {\kappa }_x\\ {\kappa }_y\\ {\kappa }_{xy}\end{array}\right)$$

(5)

Coefficients related to the extensional stiffness elements are evaluated according to:

$$A_{ij}={\displaystyle \sum _{k=1}^N{C}_{ij}\left({\theta }_k\right)\left(z_k-z_{k-1}\right)}$$

(6)

the coupling extension bending components are calculated as:

$$B_{ij}={\displaystyle \sum _{k=1}^N{C}_{ij}\left({\theta }_k\right)\left(z_k^2-z_{k-1}^2\right)}$$

(7)

and the bending stiffness elements are evaluated as

$$D_{ij}=={\displaystyle \sum _{k=1}^N{C}_{ij}\left({\theta }_k\right)\left(z_k^3-z_{k-1}^3\right)}$$

(8)

where $C_{ij}\left({\theta }_k\right)$ are the elements of the stiffness matrix of a lamina expressed in the laminate's coordinate axes, $z_k$ indicates the relative position of a layer respect to the midplane of the laminate, and ${\theta }_k$ is the orientations angle associated to the stacking sequence.

2.3. Problem formulation

The goal is to obtain the optimal orientations design in a stacking sequence which fulfills at the same time both the maximum strain criterion and the Tsai-Wu failure criterion. Thus, an objective function is proposed, based on the strains and the stresses that involves a minimization process to deformations and a maximization to the safety factor. The optimization process is carried out on a symmetric laminate, so it is enough to consider half of the plies, and by the nature of this kind of laminate, values $B_{ij}$ will be zero.

As a benchmark, a weakest-link failure approach was adopted regarding the first ply of the laminate, which is subjected to the highest stresses. The Tsai-Wu safety factor is maximized. The design variables are the orientations of the layers, and the number of plies and ply thickness are kept constant. To assure the static strength of the laminate, the safety factor of all plies is also calculated.

The Tsai-Wu failure criterion indicates when a ply subjected to an in-plane stress condition will fail. It happens if and only if the following expression is satisfied:

$$F_1{\sigma }_1+F_2{\sigma }_2+F_{11}{\sigma }_1^2+F_{22}{\sigma }_2^2+F_{66}{\tau }_{12}^2+F_{12}{\sigma }_1{\sigma }_2\ge 1$$

(9)

where ${\sigma }_1$, ${\sigma }_2$, and ${\tau }_{12}$ are the first ply stresses obtained from Eq. (2) The terms $F_i$ can be evaluated as: $F_1=\frac{1}{{\sigma }_{\mathrm{l}\max }^T}-\frac{1}{{\sigma }_{\mathrm{l}\max }^C}$,$F_2=\frac{1}{{\sigma }_{\mathrm{t}\max }^T}-\frac{1}{{\sigma }_{\mathrm{t}\max }^C}$, $F_{11}=\frac{1}{{\sigma }_{\mathrm{l}\max }^T{\sigma }_{\mathrm{l}\max }^C}$, $F_{22}=\frac{1}{{\sigma }_{t\max }^T{\sigma }_{\mathrm{t}\max }^C}$, $F_{66}=\frac{1}{{\tau }_{lt\max }^2}$, and $F_{12}=-\frac{1}{2}\sqrt{F_{11}{F}_{22}}$.

Specifically, the Tsai-Wu safety factor is defined as:

$$SF\left(\theta \right)=\frac{-b\pm \sqrt{b^2-4ac}}{2a}$$

(10)

The coefficients of Eq. (10) are evaluated as: $a=F_{11}{\sigma }_1^2+F_{22}{\sigma }_2^2+F_{66}{\tau }_{12}^2+F_{12}{\sigma }_1{\sigma }_2$, $b=F_1{\sigma }_1+F_2{\sigma }_2$, and $c=-1$.

Notice that there is an implicit dependence of $SF$ on the orientation of a ply as the stresses depend on the stiffness matrix which is angle and thickness-dependent.

In a composite material of $N$ different plies with non-periodic orientations and located at equidistant positions from the middle, the objective function is constituted by the contributions of the safety factor of all plies in the stacking sequence according to:

$$OF\left(\theta ;p\right)={\displaystyle \sum _{k=1}^{N}SF\left({\theta }_k\right)}$$

(11)

To maximize the objective function, it is necessary to include the constraints that complement the safety factor criteria. The optimization problem can be announced as:

$$\begin{array}{l}MaximizeOF\left(\theta ;p\right)\\ s.t.\\ SF\left({\theta }_k\right)-1>0,\hspace{1em}\forall k=1,2,\dots ,N\end{array}$$

(12)

where $p$ is the vector of constant parameters, such as the material properties which do not change during the optimization process, ${\theta }_k$ is the orientation angle of each ply that is allowed to vary in the range $-180°\le {\theta }_k\le 180$.

3. Results and discussion

3.1. Application cases

3.1.1. In-plane loads

The study case 1, with the aim to verify the robustness of both optimization algorithms was a composite laminate under in-plane loads. The shear load per unit length $N_{xy}$ is higher than the axial loads $N_x$ and $N_y$. Shear load is twice the axial loads per unit length. This load case is considered as reference laminate the following configuration: [±45/0/90/±45]_s. This layup is intended to have most of the plies with a 𝜃= 45º to increase the shear strength and add a couple of 0º and 90º layers to improve axial stiffness. The reference laminate has SF = 1.57 at its weakest ply and an objective function value of -3.5362 (see Table 1).

Figure 2 and Figure 3 show, respectively, i) the results for objective function values and ii) number of evaluations of objective function. From Figure 2, it can be highlighted that GA is able to find solutions with the same value of objective function, -7.2200, which is a clear indicator of the good stability of GA with respect to the input parameters (see Table 2). The results for SA fall in a range between -5.7274 and -6.9932, revealing a lower stability of the algorithm and its drawback to escape from local optimum (see Table 3). Nonetheless, solutions found by SA are clearly superior to the initial solution.

Comparing both algorithms, the value of the objective function by SA is lower than GA, but has a high variability. GA average function value is -7.2200 and a variability coefficient of 7.14% (see Table 2). The SA has an average function value of -6.3996 and a variability of 5.8% (see Table 3). The dispersion for the GA is despicable compared with the ones obtained by SA. Although SA solutions are satisfactory, these results show lower performance for the in-plane load condition than GA solutions.

Execution time required to obtain the optimal solution using SA is lower than GA. SA time gaining is almost twice fast as GA. By nature, GA needs more computational time to compute optimal solutions.

There are no significant differences between the layup configuration for the worst and best solution for GA solutions. Those solutions keep the same proportion of 8 plies with orientation near = 45º and 4 plies with near = 135º (see Table 2). These results show the high stability of the algorithm. A similar trend is found for SA solutions, but it is probable that this algorithm got stuck in a local optimal before it can improve the orientation values. In this study case, the initial solution is in the form of [±45/0/90/±45]_s, and the best solution obtained is [(45/135)₂/45₂]_s, which means that plies at 0º and 90º are not necessary, even when there are axial loads applied (see Table 3). Due to the magnitude and ratio of loading case, mainly shear, a ±45° based configuration is sufficient to grant the laminate strength. Solutions obtained by both algorithms, according to the optimization criterion and based on the SF of the weakest ply, represent better configurations than the reference laminate, however GA solutions remain superior.

3.1.2. Combined moments

The study case 2, composite laminate was charged combining bending and torsion moments. In this study case, the torsion moment has the highest magnitude. Initial solution of the laminate has a layup sequence as follows: [-45/45/45/-45/ 0/90]_s. Layup sequence has ±45 plies at the bottom and the top of the laminate to increase the in-plane shear stiffness as shown in Table 4. This configuration has a SF of 0.7240 in the weakest layer and an objective function value of -2.5567. Because one of the layers has a SF that does not satisfy the Tsai-Wu criterion, the objective function value, with this layup sequence, was penalized by diminishing its value.

From Figure 4, GA found 5 local optimal solutions with the following objective function values: -5.4950, -5.5206, -5.6228, -5.6230 and -5.6532. It can be highlighted that only one execution results in a value of -5.4950, meanwhile 16 executions have objective functions values below -5.6.

The response of GA to input parameters is stable, by obtaining solutions that asymptotically approach one objective function value. The average value of the objective function in this study case is -5.0672 with a variation coefficient of 0.8824%.

Different objective function values were computed in all solutions for the SA algorithm. This trend showsSA's difficulty in escaping from local optima and converge. SA objective values oscillate between -3.3396 and -5.3112, with an average value of -4.4810 and a coefficient of variation of 11.43%. Despite most SA solutions being satisfactory, statistical data show a lower stability and performance compared to GA.

From Figure 5, the GA required ten times more evaluations to find the optimal solution than SA. The GA needed to process at least 17200 evaluations, meanwhile SA could afford it with 3152 evaluations. Comparing the computational cost, the SA obtained a solution in an average time of 21.60 s, while GA took 35.51 s to find a solution.

As shown in Table 5 the best and worst solution from the optimization algorithms are proposed with 45º and 135º (or -45º) lamina. The position of the (-45º) laminate has a considerable impact, because although most of the solutions satisfy the SF, the objective function value depends on the position of this particular orientation.

With the results shown in Table 6, solutions proposed by the SA are better than the reference laminate. However, SA solutions have a lower value of objective function, which is noticeable on the weakest lamina showing a lower SF. The worst solution proposed by SA has a lamina with a SF < 1, which does not fulfill the Tsai-Wu criterion. Even when SA solutions trend to lamina with = 45° and = -45º, it is possible that for those solutions the algorithm got stuck in a local optimum.

In this study case, the reference solution has [-45/45/45/-45/0/90]_s as a layup sequence. The best solution obtained is of the type [(45/-45)₂/45₂]_s. For this load case, it is also unnecessary to add plies at = 0° and = 90º, even when there are axial loads applied. This load case shows that optimization algorithms can discard "common practices layups" in which plies with orientation of 0° and 90° are always added in the laminate configuration

3.1.3. In-plane loads and combined moments

In the study case 3, composite laminate was loaded with in-plane loads and combined moments. Typically, this is the loading case for which most composite plate structures are designed. For this case, the reference layup sequence is as follows: [0/45/-45/45/0/90]_s, as depicted in Table 7. The layup sequence was selected by placing the 0º and 45º plies near the laminate's surface to increase the bending and torsion stiffness. Plies at 0º and 90º are near the neutral to strengthen the imposed membrane loads. This composite has an average objective function value of -0.2887 and a SF of 0.6940 on its weakest ply. As it can be seen, the value of the objective function for this solution is significantly reduced because one of its plies does not comply with the Tsai-Wu criteria.

Considering the evolution of the objective function for the optimal solution, all GA solutions have a higher value, nonetheless, SA solutions are very close to them. Both algorithms obtained better quality solutions than the reference one, however, SA showed more variability. The average value for the objective function by GA is -2.0271 with a coefficient of variation of 3.30%. The average value for the objective function by SA is -1.7695 with a coefficient of variation of 17.69%.

Similarly to case 1 and 2, the GA required a higher number of evaluations of the objective function compared to SA. For case 3, GA required at least 17800 function evaluations with a maximum of 122400. In contrast, SA required at least 3197 evaluations with a maximum of 12503. The number of evaluations for GA is almost one order of magnitude bigger than SA. The number of evaluations needed to find the optimal solution is reflected in the computational time. In this sense the SA had a better solution time, from 8.43 seconds to 29.98 seconds, while the GA spent from 21.67 seconds to 126.65 seconds. The difference in computational time is evident.

In this study case, GA solutions promote layers near = 22° and = 24° at the surface, as it can be seen in Table 8. Plies from 3^rd to 6^th have complementary angles. For example, the 3rd layer of the best solution has = 145º, and the same layer of the worst solution has = 27º.

It can be noticed in Table 8, that the GA solution with the worst value of objective function has at its less resistant layer (1^st ply) a SF superior to the less resistant layer of the best solution (ply 10^th), which is at its strength limit. Nonetheless, this last solution is considered of better quality for the optimization criterion employed.

The 3^rd and 4^th layers of the solution with the highest value of objective function, are those that cause this increment. Both solutions satisfy the Tsai-Wu criterion and represent better options than the reference one (see Table 8).

On the other hand, it can be noticed in Table 9 that SA proposes similar orientation angles for the first four layers in its best solution compared with GA. The worst solution by SA has the lowest value of the objective function. This solution presents a ply with a SF = 0.93, which does not comply with the Tsai-Wu criterion (see Table 9). It can be noticed that the orientation angles for the best and worst SA solutions are very different. It can be implied that the SA algorithm could be stuck in a local optimum, being unable to improve the orientation angles, raising the value of the objective function and, therefore, being penalized.

Although the lower performance to find solutions for the considered optimization criterion, the SA algorithm fulfills to obtain on 18 of 20 executions, solutions that were not penalized and which are considered better than the reference layup.

In this study case the original layup sequence is [0/45/-45/45/0/90]_s, meanwhile the best solution proposed by the optimization algorithms is [21/24/145/140/45/49]_s. An interesting observation is that if the layers are regrouped in pairs, each orientation has a gap of 5° . In this study case, the optimal solution shows that the inclusion of plies with = 0º or 90º is not mandatory. Each load case can have a related layup sequence, where the correct angle orientation can significantly improve the performance of the composite structure.

Figure 6. Objective function values for GA and SA for study case 3.

Figure 6. Objective function values for GA and SA for study case 3.

Figure 7. Number of evaluations for GA and SA for study case 3.

Figure 7. Number of evaluations for GA and SA for study case 3.

4. Conclusions

This paper studied the performance of Genetic Algorithm (GA) and Simulated Annealing (SA) as optimization techniques for calculating fiber-reinforced composites loaded under consideration of classic laminate theory. The main objective is to obtain the optimal orientation layup, which fulfills at the same time both the maximum strain criterion and the Tsai-Wu failure criterion for different loading cases.

GA presented a good performance but with a higher computational cost concerning SA. In all study cases, GA found solutions that overpass the expectations of the initial (reference) solution. GA also showed great stability, any input data do not influence it, and in most cases, it was able to find at different executions, convergence solutions.

SA was able to find optimal solutions with a good value of objective function. Nonetheless, it is strongly influenced by the input parameters. The variability is higher than GA, but it can obtain solutions in a very short time. SA optimized solutions are also always better than the proposed ones.

GA solutions have a higher value of objective function than SA. SA has a higher variability in the results with respect to GA, which suggests a more limited capability of the algorithm to escape from local optima.

For the first load case, in-plane loads, the initial solution is in the form of [±45/0/90/±45]_s, and the best solution obtained by both optimization algorithms is in the form of [(45/135)₂/45₂]_s. For the second load case, combined moments, [-45/45/45/-45/0/90]_s is the initial layup sequence, then the best solution obtained provided by both optimization algorithms is in the form of [(45/-45)₂/45₂]_s. For the third study case, combined moments, the original layup sequence is [0/45/-45/45/0/90]_s, meanwhile the best solution calculated by the optimization algorithms is [21/24/145/140/45/49]_s. In all three cases, plies at 0º and 90º are unnecessary even when applied axial loads.

Both techniques, GA and SA, have proven advantages and disadvantages. However, both algorithms found better solutions than the reference ones.

It is also highlighted that layup solutions rarely contain plies oriented near the common practice angles in all three study cases. With this study, a slight variation in the orientation angles could significantly impact mechanical performance. Exploration of these non-common angles would be challenging without tools such as GA or SA optimization techniques.

The use of optimization techniques can improve design time and cost. Further work will consider an experimental comparison between the reference layup versus the optimal layup calculated.