*Article* **FEM-Validated Optimal Design of Laminate Process Parameters Based on Improved Genetic Algorithm**

**Xing Mou, Zhiqiang Shen, Honghao Liu, Hui Xv, Xianzhao Xia and Shijun Chen \***

School of Mechanical Engineering, Hefei University of Technology, Hefei 230000, China; mouxing2021@163.com (X.M.); shenzhiqiang158@foxmail.com (Z.S.); 13700588738@163.com (H.L.); xuhuicslg@163.com (H.X.); xiaxz4557@126.com (X.X.) **\*** Correspondence: newbonchen@hotmail.com; Tel.: +86-136-85-519-410

**Abstract:** In tape placement process, the laying angle and laying sequence of laminates have proven their significant effects on the mechanical properties of carbon fibre reinforced composite material, specifically, laminates. In order to optimise these process parameters, an optimisation algorithm is developed based on the principles of genetic algorithms for improving the precision of traditional genetic algorithms and resolving the premature phenomenon in the optimisation process. Taking multi-layer symmetrically laid carbon fibre laminates as the research object, this algorithm adopts binary coding to conduct the optimisation of process parameters and mechanical analysis with the laying angle as the design variable and the strength ratio R as the response variable. A case study was conducted and its results were validated by the finite element analyses. The results show that the stresses before and after optimisation are 116.0 MPa and 100.9 MPa, respectively, with a decrease of strength ratio by 13.02%. The results comparison indicates that, in the iterative process, the search range is reduced by determining the code and location of important genes, thereby reducing the computational workload by 21.03% in terms of time consumed. Through multiple calculations, it validates that "gene mutation" is an indispensable part of the genetic algorithm in the iterative process.

**Keywords:** composite materials; genetic algorithm; important genes; strength ratio; finite element analysis

### **1. Introduction**

Carbon fibre reinforced composite materials (CFRP) are widely used in the fields of aviation equipment structures [1], fuselages and wings manufacturing due to their high corrosion resistance, high strength and lightweight [2,3]. As a typical representative of the practical application of composite materials, laminates have a direct impact on the mechanical properties of CFRPs, which depend highly on their laying angle and laying sequence. Hence, many investigations have been conducted in-depth on exploring the influence of these two variables and searching for their optimal combinations.

A genetic algorithm (GA) is a random search method abstracted by simulating gene selection and genetic mechanisms in the natural environment. It can not only solve linear problems but also can be used to solve nonlinear problems; and genetic algorithm is sensitive to the problem of discretization, so it has attracted the attention of scholars [4], and it has been applied to many occasions. At present, scholars worldwide favor the use of genetic algorithms to optimise the design of the mechanical properties of laminates, and certain research results have been obtained. Wang et al. [5] applied a multi-island genetic algorithm to optimise the damage of different loads to the laminate during low-velocity impact and optimised the laminate sequence. Strain energy was taken as the optimisation objective and the laying angle as the design variable. Their results show that the impacted area was reduced by 42.1%, and the residual compression strength increased by 10.79%. It is notable that the laying angle was optimised as a continuous variable rather than discrete

**Citation:** Mou, X.; Shen, Z.; Liu, H.; Xv, H.; Xia, X.; Chen, S. FEM-Validated Optimal Design of Laminate Process Parameters Based on Improved Genetic Algorithm. *J. Compos. Sci.* **2022**, *6*, 21. https:// doi.org/10.3390/jcs6010021

Academic Editors: Stelios K. Georgantzinos and Francesco Tornabene

Received: 17 December 2021 Accepted: 4 January 2022 Published: 7 January 2022

**Publisher's Note:** MDPI stays neutral with regard to jurisdictional claims in published maps and institutional affiliations.

**Copyright:** © 2022 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/).

as in the other investigations mentioned above. Chen et al. [6] aimed at optimising the fibre laying angle to achieve the minimum stress of the laminate, and used a genetic algorithm to obtain the best laying angle and minimum stress of the laminate. Their method realised the optimisation of the laying angle of two-layer laminates and three-layer laminates by taking the design variable as a continuous variable.

In order to improve the reliability, some optimisation strategies were proposed by Xiu and Cui [7], which took the buckling load of laminate as the optimisation objective, and the laying angle of laminate as constraint condition, and used the neural network model and genetic algorithm to optimise the ratio of laying angle to the buckling load of the laminate. The optimisation method proposed in their investigation is capable of producing a better laying sequence with the same material weight under the condition that the laying angles are given and the number of layers is a constant; nevertheless, the critical process parameter, laying angle, cannot be optimised by this method and is just taken as a given constant. Feng et al. [8] took the strength of laminate as the optimisation objective, and the laying angle of laminate as the constraint condition. In order to obtain a reasonable laying angle and ply number, Jin et al. [9] proposed a two-step optimisation strategy and developed a three-level optimisation model with embedded genetic algorithm and Tsai-Wu failure criterion. The concept of sublayer was introduced into the optimisation method for the optimisation of the location of layers, geometric dimensions, number of layers, and laying sequence. The authors claimed that the introduction of this concept improved the manufacturability of laminates, however, the laying angle was not optimised in this work. By embedding the Tsai-Hill failure criterion, Park et al. [10] used genetic algorithms to apply different loads to the symmetrical composite laminates under different boundary conditions to perform mechanical optimisation design. Tournament selection and uniform crossover methods were adopted in the optimisation process. The laying angle was taken as the design variable as a continuous variable, and a low failure index was obtained. The results showed that the modified genetic algorithm was capable of identifying a global optimal solution for the parameter's optimisation of laid laminates.

Furthermore, new encoding approaches and operators have been developed for genetic algorithms to accelerate the convergence speed and computational accuracy. A genetic algorithm was used to calculate the fibre laying angle and validated that the calculation results can effectively improve the strength of the laminate. This work improved the genetic algorithm by modifying the encoding approach and the integer encoding was embedded in the iteration process for the laminate parameters optimisation. The advantage of this modification is that integer encoding has the potential to avoid unnecessary iterations of floating number, thereby improving the computational efficiency theoretically. The laying angle in this optimisation was specified to 0◦, ±45◦, and 90◦, leading to weakening the designability of laminates to some extent. By designing a new mutation operator and crossover operator, Wang et al. [11] adopted a symbol encoding method for the development of a self-adaptive genetic algorithm with the bending stiffness parameter as the fitness function. They deeply explored the influence of the adaptively changing crossover operator and the mutation operator on the genetic algorithm and further optimised the laying order. Compared with the traditional genetic algorithm, their modification enhanced the reliability, convergence speed, and computational efficiency of the optimisation process. Another self-adaptive genetic algorithm was proposed by Yang et al. [12] through the automatic change of the mutation rate and crossover rate according to fitness value. This algorithm took the maximum deformation of laminate as optimisation objective, the laminate angle as constraint condition, and called ABAQUS to analyse the stress of the laminate. The results show that this algorithm is easier to get convergence than the traditional genetic algorithm. Hwang et al. [13] designed an elite comparison operator for the genetic algorithm to accelerate the iteration process for the optimisation of laminates. This elite comparison operator was used to identify the differences in the design variables from two different generations, thereby maintaining the similar ones and changing the different ones to find their possible values. Two optimisation problems were solved by this modified genetic algorithm and the results

indicate that the designed operator has the potential to speed up the convergence of the optimisation process.

The above investigations mainly focus on the change of optimisation parameters, design of encoding methods, and improvement of operators in the application of traditional genetic algorithms. However, in order to improve the precision of traditional genetic algorithms and eliminating the premature phenomenon in the optimisation process, this article improves the genetic algorithm (improved genetic algorithm, referred as IGA below) on the basis of the principles of traditional genetic algorithms through the application of binary coding to explore important gene codes and their locations. In this IGA, the laying angle is set as the constraint condition, and the minimum stress as the optimisation objective, which customised this algorithm for the optimisation of process parameters of laminates tape placement. Based on the Tsai-Wu failure criterion [9], the laying angle and laying sequence were optimised; and the obtained stress was compared with that of the traditional genetic algorithm to validate the IGA. In addition, the optimal solutions were further validated by a finite element analysis of a laminate with the same geometry, loading, and boundary conditions.

### **2. Mathematical Model of Laminate and the IGA**

*2.1. On- and off-Axis Stress Conversion of Laminates*

Generally speaking, a multi-layer laminate can be formed by stacking multiple monolayer boards at different laying angles. According to the designability of composite materials and the theoretical basis of elastic mechanics, the on-axis direction of the monolayer boards can be laid at different angles and in different orders according to the loading conditions to meet the requirements of strength and stiffness design of the structural parts under the load conditions.

Research for the classical laminate theory (CLT), the establishment and derivation of formulae are mainly based on the following assumptions:


For simplification of the formula derivation process, the median plane of the laminate is taken as the reference plane, viz. the plane perpendicular to the z-direction at the middle point of the thickness of the laminate. According to the above assumptions, let the stress function *x* and *y* be *ϕ*, then the corresponding stress components *σx*, *σy*, *τxy*, are given in Equation (1) as below:

$$
\sigma\_x = \frac{\partial \varphi^2}{\partial y^2}, \sigma\_y = \frac{\partial \varphi^2}{\partial x^2}, \quad \tau\_{xy} = \frac{-\partial \varphi^2}{\partial x \partial y} \tag{1}
$$

Since the classic laminate theory only considers plane stress and ignores the existence of interlayer stress in the iterative process, the stiffness coefficient of each layer of the laminate is different so that the internal force of the laminate can only be integrated layer by layer. Assuming that the stress of the *kth* layer of the laminate is *σk* . , the internal force of the laminate can be expressed as:

$$\left\{ \begin{array}{c} N\_{\mathbf{x}} \\ N\_{\mathbf{y}} \\ N\_{\mathbf{xy}} \end{array} \right\} = \sum\_{k=1}^{n} \int\_{z\_{k-1}}^{z\_{k}} \left\{ \begin{array}{c} \sigma\_{\mathbf{x}}^{k} \\ \sigma\_{\mathbf{y}}^{k} \\ \tau\_{\mathbf{x}}^{k} \end{array} \right\} dz\_{\mathbf{z}} \; \left\{ \begin{array}{c} M\_{\mathbf{x}} \\ M\_{\mathbf{y}} \\ M\_{\mathbf{xy}} \end{array} \right\} = \sum\_{k=1}^{n} \int\_{z\_{k-1}}^{z\_{k}} \left\{ \begin{array}{c} \sigma\_{\mathbf{x}}^{k} \\ \sigma\_{\mathbf{y}}^{k} \\ \tau\_{\mathbf{x}}^{k} \end{array} \right\} z dz \tag{2}$$

The stress–strain relationship of the *kth* layer of the laminate is as follows:

$$\left\{ \begin{array}{c} \sigma\_{\mathbf{x}}^{k} \\ \sigma\_{\mathbf{y}}^{k} \\ \tau\_{\mathbf{xy}}^{k} \end{array} \right\} = \left[ \begin{array}{c} \overline{Q\_{11}^{k}} & \overline{Q\_{12}^{k}} & \overline{Q\_{16}^{k}} \\ \overline{Q\_{21}^{k}} & \overline{Q\_{22}^{k}} & \overline{Q\_{26}^{k}} \\ \overline{Q\_{61}^{k}} & \overline{Q\_{62}^{k}} & \overline{Q\_{66}^{k}} \end{array} \right] \left\{ \begin{array}{c} \varepsilon\_{\mathbf{x}}^{k} \\ \varepsilon\_{\mathbf{y}}^{k} \\ \gamma\_{\mathbf{xy}}^{k} \end{array} \right\} \tag{3}$$

Substituting Equation (2) into Equation (1), and through integration, the stress–strain relationship of the laminate can be obtained as follows:

$$\begin{Bmatrix} \begin{Bmatrix} N\_{\mathcal{X}} \\ N\_{\mathcal{Y}} \\ N\_{\mathcal{X}\mathcal{y}} \\ M\_{\mathcal{X}} \\ M\_{\mathcal{Y}} \\ M\_{\mathcal{X}\mathcal{y}} \end{Bmatrix} = \begin{bmatrix} A\_{11} & A\_{12} & A\_{16} & B\_{11} & B\_{12} & B\_{16} \\ A\_{21} & A\_{22} & A\_{26} & B\_{21} & B\_{22} & B\_{26} \\ A\_{61} & A\_{62} & A\_{66} & B\_{61} & B\_{62} & B\_{66} \\ B\_{11} & B\_{12} & B\_{16} & D\_{11} & D\_{12} & D\_{16} \\ B\_{21} & B\_{22} & B\_{26} & D\_{21} & D\_{22} & D\_{26} \\ B\_{61} & B\_{62} & B\_{66} & D\_{61} & D\_{62} & D\_{66} \end{Bmatrix} \begin{Bmatrix} \varepsilon\_{\mathcal{X}}^{0} \\ \varepsilon\_{\mathcal{Y}}^{0} \\ \gamma\_{xy}^{0} \\ k\_{x} \\ k\_{xy} \end{Bmatrix} \tag{4}$$

where, *Aij*, *Bij*, *Dij* are

$$\begin{aligned} A\_{ij} &= \sum\_{k=1}^{n} \int\_{z\_{k-1}}^{z\_k} \overline{Q\_{ij}^k} dz = \sum\_{k=1}^{n} \overline{Q\_{ij}^k} (z\_k - z\_{k-1}) \\ B\_{ij} &= \sum\_{k=1}^{n} \int\_{z\_{k-1}}^{z\_k} \overline{Q\_{ij}^k} dz = \frac{1}{2} \sum\_{k=1}^{n} \overline{Q\_{ij}^k} (z^2\_k - z^2\_{k-1}) \\ D\_{ij} &= \sum\_{k=1}^{n} \int\_{z\_{k-1}}^{z\_k} \overline{Q\_{ij}^k} z^2 dz = \frac{1}{3} \sum\_{k=1}^{n} \overline{Q\_{ij}^k} (z^3\_k - z^3\_{k-1}) \end{aligned} \tag{5}$$

where, *i*, *j* = 1, 2, 6; *Aij* is the in-plane stiffness coefficient, *Bij* is the coupling stiffness coefficient, and *Dij* is the bending stiffness coefficient.

The flexibility matrix of the laminate is

$$S' = \begin{bmatrix} A' & B' \\ H' & D' \end{bmatrix} \tag{6}$$

where, *<sup>A</sup>* = *<sup>A</sup>*−<sup>1</sup> + *<sup>A</sup>*−1*B*(−*BA*−1*<sup>B</sup>* + *<sup>D</sup>*) −1 *BA*<sup>−</sup>1; *<sup>B</sup>* = *<sup>B</sup>*−<sup>1</sup> − *<sup>A</sup>*−1*BD*; *<sup>H</sup>* = *<sup>B</sup>*−<sup>1</sup> − *<sup>D</sup>*−1*BA*<sup>−</sup>1; *<sup>D</sup>* = (−*BA*−1*<sup>B</sup>* + *<sup>D</sup>*) −1 .

As the presence of the matrix [B] indicated, the laminate is prone to the coupling of bending deformation and tension and compression deformation when subjected to force; in addition, in-plane deformation is relatively easy to occur when subject to bending stress, and further leads to the warpage of the laminate after curing (that is, *Bij* = 0, which means it is hard for the calculation of laminates). For this reason, the composite material laminates generally used in engineering should be symmetrically laid (*Bij* = 0) to avoid the calculation difficulties caused by the coupling stiffness [14]. Note that symmetry is assumed in this study to make sure the laminate only subject to tensile strain but no flexure [14], though it is difficult to always keep symmetry exactly in industry. Therefore, the flexibility matrix S of the laminate at this time can be simplified as:

$$S' = \begin{bmatrix} A^{-1} & 0 \\ 0 & D^{-1} \end{bmatrix} \tag{7}$$

At the same time, through the stress conversion matrix [*Tσ*], the off-axis stress of each layer can be converted into the on-axis stress. The specific conversion formula is shown in Equation (8):

$$\left\{ \begin{array}{c} \sigma\_1\\ \sigma\_2\\ \tau\_{12} \end{array} \right\} = \left[ \begin{array}{cc} m^2 & m^2 & 2mn\\ n^2 & m^2 & -2mn\\ -mn & mn & m^2 - n^2 \end{array} \right] \left\{ \begin{array}{c} \sigma\_x\\ \sigma\_y\\ \tau\_{xy} \end{array} \right\} \tag{8}$$

where, *m* = cos *θ*, *n* = sin *θ*, and *θ* is the laying angle of each layer of fibres.

### *2.2. Failure Criterion for Laminates*

After comprehensively considering the advantages and disadvantages of multiple strength criteria, Tsai and Wu [15] proposed a new strength criterion based on the tensor form. That is, the matrix breaks based on Tsai-Wu criterion, as below:

$$F\_{11}\sigma\_1^2 + F\_{22}\sigma\_2^2 + F\_{66}\sigma\_6^2 + 2F\_{12}\sigma\_1\sigma\_2 + F\_1\sigma\_1 + F\_2\sigma\_2 = 1\tag{9}$$

where, *σ*<sup>1</sup> is the longitudinal stress, *σ*<sup>2</sup> is the transverse stress, *σ*<sup>6</sup> is the shear stress, *XT* is the longitudinal tensile strength, *XC* is the longitudinal compressive strength, *YT* is the transverse tensile strength, *YC* is the transverse compressive strength, S is the shear strength; *F*1, *F*11, *F*2, *F*22, and *F*<sup>66</sup> are strength parameters as shown below in Equation (10).

$$\begin{aligned} F\_1 &= \frac{1}{X\_{\overline{T}}} - \frac{1}{X\_{\overline{C}}} & F\_{11} &= \frac{1}{X\_{\overline{C}} X\_{\overline{T}}} \\ F\_2 &= \frac{1}{Y\_{\overline{T}}} - \frac{1}{Y\_{\overline{C}}} & F\_{22} &= \frac{1}{Y\_{\overline{T}} Y\_{\overline{C}}} \\ F\_{66} &= \frac{1}{S^2} \end{aligned} \tag{10}$$

The value of *F*<sup>12</sup> here has little effect on the calculation of the strength criterion, and this is validated by the investigation cited here [16]. Then let *F*<sup>12</sup> = 0, the calculation result of Equation (9) and that of the Hoffman strength criterion are relatively close with an error of 10%. Hence, Equation (9) is improved as below in Equation (11).

$$\left| F\_{11} \sigma\_1^2 + F\_{22} \sigma\_2^2 + F\_{66} \sigma\_6^2 + 2F\_{12} \sigma\_1 \sigma\_2 + 2F\_{16} \sigma\_1 |\sigma\_6| + 2F\_{26} \sigma\_2 |\sigma\_6| + F\_1 \sigma\_1 + F\_2 \sigma\_2 = 1 \right. \tag{11}$$

The value of *F*16, *F*<sup>26</sup> here is obtained by a semi-empirical formula [17], and the value of *F*<sup>12</sup> is the value when the Tsai-Wu strength criterion and the Hoffman criterion are the same. See Equation (12) for details.

$$F\_{16} = \frac{1}{\frac{1}{5S\sqrt{X\_T X\_C}}} \quad F\_{26} = \frac{1}{\frac{1}{5S\sqrt{Y\_T Y\_C}}} \quad \left\{ \begin{array}{c} \\ \\ \\ \end{array} \right\} \tag{12}$$

$$F\_{12} = -\frac{1}{2} F\_{11} \end{array} \tag{12}$$

Note that the value of *<sup>F</sup>*<sup>12</sup> here is not <sup>−</sup>1/2√*F*11*F*<sup>22</sup> because:


$$\underbrace{\begin{matrix} \cdot \cdot \cdot \cdot \cdot \cdot \cdot \cdot \cdot \cdot \cdot \cdot \cdot \cdot \cdot \cdot \cdot \\ \cdot \cdot \cdot \cdot \cdot \cdot \cdot \cdot \cdot \cdot \cdot \cdot \cdot \cdot \cdot \cdot \\ \cdot \cdot \cdot \cdot \cdot \cdot \cdot \cdot \cdot \cdot \cdot \cdot \cdot \end{matrix}}\_{\cdot \cdot 1}$$

**Figure 1.** Strength curve of Tsai-Wu failure criterion [15].

From Equation (11), when the calculated value on the left side of the equation equates 1, it indicates that the material has just reached the critical failure state; when the left side is greater than 1, it indicates that damage has begun; when the left side is less than 1, it indicates that the material is in a safe state, and there is still a certain safety margin.

In order to obtain better calculation speed and expected values in genetic algorithms, a suitable. fitness value is generally specified. Here, the strength ratio is applied to the laminate, namely

$$R = \frac{\sigma\_{\bar{l}}}{\sigma} \tag{13}$$

where *R* is the ratio of a certain component of the ultimate stress to the corresponding component of loading stress, and *i* = 1, 2, 6.

The meaning of "correspondence" here [18] assumes that the loading is linear, that is, the stress components increase simultaneously in a certain proportion, and this situation is in line with the actual application in industrial production. Transforming Equation (13) to get *Rσ* = *σi*, substitute this into Equation (11) and after simplification to produce Equation (14) below.

$$\begin{aligned} \left(F\_{11}\sigma\_1^2 + F\_{22}\sigma\_2^2 + F\_{66}\sigma\_6^2\right)\mathcal{R}^2 + 2(F\_{12}\sigma\_1\sigma\_2 + F\_{16}\sigma\_1|\sigma\_6| + F\_{26}\sigma\_2|\sigma\_6|)\mathcal{R}^2\\ + (F\_1\sigma\_1 + F\_2\sigma\_2)\mathcal{R} = 1 \end{aligned} \tag{14}$$

This is a quadratic equation in one variable. It is easy to get Equation (15) as below.

$$R = \frac{-b + \sqrt{b^2 + 4a}}{2a} \tag{15}$$

where:

$$\begin{aligned} a &= F\_{11}\sigma\_1^2 + F\_{22}\sigma\_2^2 + F\_{66}\sigma\_6^2 + 2(F\_{12}\sigma\_1\sigma\_2 + F\_{16}\sigma\_1|\sigma\_6| + F\_{26}\sigma\_2|\sigma\_6|) \\ b &= F\_1\sigma\_1 + F\_2\sigma\_2 \end{aligned} \tag{16}$$

Equation (15) is called the strength ratio equation. From the definition of the strength ratio, it can be seen that when *R* > 1, the material has not failed; when the applied stress is increased to (*R* − 1) times, the material fails. R cannot be less than 1 because this is mathematically meaningless; on contrast, *R* < 1 in engineering applications still has a certain reference value, that is, it indicates that the applied stress must be reduced or the size of the structure must be increased to ensure that the structure is not damaged.

### *2.3. Genetic Algorithm and Its Improvement*

Genetic algorithm is based on the natural law of survival of the fittest by simulating the theory of genetics and biological evolution proposed by Darwin, thereby retaining the excellent genes and individuals [19].

The core content of genetic algorithm mainly includes four parts, namely selection, crossover, mutation, and fitness evaluation. The optimisation process is shown in Figure 2.

**Figure 2.** Iteration process of genetic algorithm.

At present, the improvement of genetic algorithm mainly focuses on self-adaptive genetic operators, that is, generating a certain number of individuals, improving the crossover rate and genetic probability between individuals during genetic operations, and using average fitness to improve calculation efficiency. The method proposed in this paper is to start from the population; when generating the population, only one individual is generated. Then the genes of this individual mutate later in the evolution process, and the gene mutation means gene transfer from one vertex to another according to the definition of the hyperplane. The calculation is performed by gene mutation; the better genes and positions are selected through fitness, and negative feedback adjustment is used to generate individuals with important genes, and then identify the corresponding laying angle through decoding. The basic mathematical theory mentioned above is:

For the full set space composed of individuals whose gene length is *l* [20], that is, the individual space is given in Equation (17)

$$S\_{sp} = \{0, 1\}^{l} \tag{17}$$

where we regard this as the spatial vertex set with the dimension *l*. When *l* = 3, it becomes a three-dimensional cube.

As Figure 3 indicated, for three-dimensional individual space, the genes at each vertex can be connected by 12 straight lines. It can also be seen that in a straight line, we only need to change the value of one of the genes, then the gene can be converted into the gene at the other end of the same line; in any plane, as long as the values of two genes are changed, the genes at the four vertices can be switched to each other. Therefore, in a three-dimensional individual space, the straight line and the plane are called the "hyperplane", and the interconnection between them constitutes the "subspace" in the space.

**Figure 3.** Three-dimensional individual space.

Therefore, according to the above definition of "hyperplane", for the *i-th* gene, if

$$f\left(\left(\mathbf{x}\_{1\prime}\cdots\mathbf{x}\_{i-1\prime},\mathbf{l},\mathbf{x}\_{i+1\prime}\cdots\mathbf{x}\_l\right)\right)\geq f\left(\left(\mathbf{x}\_{1\prime}\cdots\mathbf{x}\_{i-1\prime},0,\mathbf{x}\_{i+1\prime}\cdots\mathbf{x}\_l\right)\right)\tag{18}$$

for any *xk* = 0, 1(*k* = *i*) holds, or

$$f\left(\left(\mathbf{x}\_{1\prime}\cdots\mathbf{x}\_{i-1\prime},0,\mathbf{x}\_{i+1\prime}\cdots\mathbf{x}\_l\right)\right)\geq f\left(\left(\mathbf{x}\_{1\prime}\cdots\mathbf{x}\_{i-1\prime},1,\mathbf{x}\_{i+1\prime}\cdots\mathbf{x}\_l\right)\right)\tag{19}$$

for any *xk* = 0, 1(*k* = *i*) holds, the gene *i* is taken as an important gene. When Equation (18) is established, important gene 1 is the better choice; while when Equation (19) is established, important gene 0 is the better choice.

Once the "hyperplane" is determined, it is possible to locate the gene value of the corresponding position after confirming that the gene value of some gene position is valid, thereby reducing the search range and further reducing the computational workload. The specific process of the IGA is shown in Figure 4.

**Figure 4.** Flow chart of the IGA.

### **3. IGA Applied to Composite Laminate Design**

*3.1. Ply Design Requirements for Laminates*

For the ply design of laminates, the coupling effect caused by stretching and bending during the curing process should be avoided since it causes resin cracking, warping deformation and even delamination for laminates. Therefore, a design generally complies with the following process design rules [21].


### *3.2. Genetic Algorithm Coding*

In the process of searching optimal solutions, the genetic algorithm generally uses a finite length string code to define the solution to be solved, such as a binary code. The layup of a laminate is a matter of discrete variables, so the choice of variables affects the length of the code string to some extent. Since the other encoding schemes, such as integer encoding and floating-point number encoding, are not unique in the subsequent mutations, they are not conducive to genetic operations such as crossover and mutation. Therefore, traditional binary encoding is selected hereof. The concrete corresponding method is:

$$\text{00: } 0^\circ; 01; \text{ } -45^\circ; 10; \text{ } 45^\circ; 11; \text{ } 90^\circ$$

Since the initial population generates randomly, there must be a certain number of initial individuals that do not meet the requirements of the laminate layup design; these individuals must be eliminated, and the others that meet the requirements should remain. The subsequent operations conduct all the above operations for each subsequent generation, thereby ensuring that the population individuals are suitable. The operation is like a "human intervention" in nature.

### *3.3. Fitness Function*

In the genetic algorithm, fitness is used to evaluate the goodness of individuals adapting to the "environment", thereby further weighing the degree to which individuals can reach or approach the optimal solution in the process of genetic evolution [22]. Therefore, the fitness function selected for this calculation is

$$R\_s = R - (5n\_1 - 10n\_2) \quad \text{ } i = 1, \text{ } 2 \dots n \tag{20}$$

where, *n* is the total number of laminates, *R*(*i*) is the strength ratio of each layer of fibre, *n*<sup>1</sup> is the number of adjacent four layers with the same laying angle, *n*<sup>2</sup> is the number of layers that do not meet the requirements of the layup ratio, 5*n*<sup>1</sup> + 10*n*<sup>2</sup> is the penalty factor. Because the initial individuals generated may not all meet the design requirements of the laminate, as mentioned above, a certain high penalty is required for the fitness factor in the searching process. Hence, through reduction of the fitness function value of unsuitable. Individuals, the probability of being selected in selection replication is further reduced. In actual operations, there may be a phenomenon that a high penalty makes *Rs* a negative value, and such a value is not conducive to the calculation of the subsequent results. Therefore, when it is a negative value, set *Rs* = 0.05, thereby forcing it to become a small positive value by "human interference" to facilitate subsequent calculations. According to the equations above, the optimisation process of the strength of laminate by the IGA algorithm is the process of analysing and calculating the best strength ratio of the single layer in the population.

### *3.4. Relevant Parameters Selection*

In this paper, T300/5228 epoxy carbon fibre with good specific strength is selected, and the material parameters are shown in Table 1.


**Table 1.** Mechanical properties of T300 epoxy resin carbon fibre.

In an optimisation process, the communication of individuals between populations is not a complete fusion, but on the basis of absorbing each other's excellent genes and simultaneously ensuring that their own excellent genes are preserved. Therefore, individuals need to use different crossover rates *PC* and mutation rates *Pm*. Furthermore, the values of *PC* and *Pm* also affect the global and local search capabilities and lead to different optimisation results; empirically, researchers recommend to adopt a large value of *PC* and small value of *Pm*. Here in this article, the two probability factors are specified as 0.95 and 0.005, respectively. The specific parameter settings related to the genetic algorithm are detailed in Table 2.

**Table 2.** Parameter settings of genetic algorithm.


At the same time, in order to speed up the optimisation process, the evaluation parameters of the adaptive genetic operator are introduced [23], namely

$$
\Omega I = \frac{R\_{\text{max}} - R\_{\text{avg}}}{R\_{\text{avg}} - R\_{\text{min}} + \varepsilon} \tag{21}
$$

where, *Rmax* is the maximum fitness value in each generation of individuals; *Ravg* is the average fitness value in each generation of the population; *Rmin* is the minimum fitness value in each generation; *ε* is an infinitesimal with the main purpose to prevent the extreme case of 0 in the denominator.

Through Equation (21), the value of *PC* and *Pm* in the two cases of *U* > 1 and *U* < 1 is further derived so as to speed up the calculation speed.

### **4. A Case Study**

### *4.1. Relevant Parameters Selection*

Given that a rectangular composite laminate beam with four sides simply supported is designed, the size of the laminate is, length *b* = *2l* with the value of 1300 mm, width *a* = 132 mm, thickness *h* = 2.25 mm, and the bearing pressure perpendicular to the surface *q* = 0.5 MPa. The design adopts a symmetrical 18-layer layup scheme. The geometric model of the laminate is illustrated in Figure 5.

**Figure 5.** Structure and geometric dimensions of the laminate model.

The schematic diagram of the simply supported laminate beam is shown in Figure 6, as below.

**Figure 6.** Loading diagram of simply supported laminate beam.

Since *σ<sup>y</sup>* is caused by the lateral load *q*, and *q* is a fixed value, which means that *σ<sup>y</sup>* is a function of *y*, here we can set (See the detailed derivation process of the equations in this section at the Appendix A.)

$$
\sigma\_y = f(y) = \frac{\partial^2 \varphi}{\partial x^2} \tag{22}
$$

Then integrating Equation (22) it becomes

$$
\sigma\_y = \frac{1}{2} \mathbf{x}^2 f(\mathbf{x}) + \mathbf{x} f\_1(\mathbf{x}) + f\_2(y). \tag{23}
$$

Substituting Equation (23) into compatibility Equation (26) it becomes

$$\begin{split} \frac{d^2}{2} \frac{d^4 f(y)}{dy^4} S\_{11} + \left( \frac{d^4 f\_1(y)}{dy^4} S\_{11} - 2 \frac{d^3 f(y)}{dy^3} S\_{16} \right) \mathbf{x} + \left( 2S\_{11} + S\_{66} \right) \frac{d^2 f(y)}{dy^2} - 2 \frac{d^2 f\_1(y)}{dy^2} S\_{16} \\ + \frac{d^4 f\_2(y)}{dy^4} S\_{11} = 0 \end{split} \tag{24}$$

where, the balance equation and compatibility equation are shown in Equations (25) and (26), and the stress function *ϕ* is shown in Equation (1) above

$$\begin{cases} \frac{\partial \tau\_x}{\partial x} + \frac{\partial \tau\_{xy}}{\partial y} = 0 \\\\ \frac{\partial \tau\_y}{\partial y} + \frac{\partial \tau\_{xy}}{\partial x} = 0 \end{cases} \tag{25}$$

$$\begin{cases} \frac{\partial \tau\_x}{\partial x} + \frac{\partial \tau\_{xy}}{\partial y} = 0 \\\\ \frac{\partial \tau\_y}{\partial y} + \frac{\partial \tau\_{xy}}{\partial x} = 0 \end{cases} \tag{26}$$

If for any *x*, Equation (23) holds, then all the coefficients of Equation (24) need to be 0, that is

$$\frac{S\_{11}}{2} \frac{d^4 f(y)}{dy^4} = 0$$

$$\frac{d^4 f(y)}{dy^4} S\_{11} - 2 \frac{d^3 f(y)}{dy^3} S\_{16} = 0 \tag{27}$$

$$(2S\_{11} + S\_{66}) \frac{d^2 f(y)}{dy^2} - 2 \frac{d^2 f\_1(y)}{dy^2} S\_{16} + \frac{d^4 f\_2(y)}{dy^4} S\_{11} = 0.$$

Hence, it can set *f*(*y*) as

$$f(y) = \frac{A\_1}{6}y^3 + \frac{A\_2}{2}y^2 + A\_3y + A\_4\tag{28}$$

According to the principle of Saint-Venant it gives

$$\begin{aligned} \tau\_{xy} &= 0 & y &= \pm \frac{h}{2} \\ \sigma\_y &= q & y &= + \frac{h}{2} \\ \sigma\_y &= 0 & y &= -\frac{h}{2} \end{aligned} \quad \left\{ \begin{array}{c} \\ \end{array} \right\} . \tag{29}$$

At the same time, the component force according to the shear stress is *Q*, namely

$$Q = -ql \quad \text{x} = l \quad . \tag{30}$$

Substituting Equation (28) into boundary condition Equations (29) and (30), and according to the axial component force in the axial direction is 0, the stress expression of *σx*, *σy*, *τxy* can be derived as [24]

$$
\sigma\_{\mathbf{x}} = \frac{12q}{h^3} \left( \frac{2S\_{12} + S\_{66}}{6S\_{11}} y^3 - \frac{1}{2} \mathbf{x}^2 y \right) + \left( \frac{6ql^2}{h^3} - \frac{3q}{10h} \frac{2S\_{12} + S\_{66}}{S\_{11}} \right) y
$$

$$
\sigma\_y = -\frac{2q}{h^3} y^3 + \frac{3q}{2h} y + \frac{q}{2} \tag{31}
$$

$$
\tau\_{xy} = \mathbf{x} \left( \frac{6q}{h^3} y^2 - \frac{3q}{2h} \right)
$$

In order to conduct optimisation, R is first obtained by combining Equations (31) and (15) according to the above requirements, and then it is substituted into Equation (20) to produce the penalty fitness function; and the MATLAB script was written. For speeding up the convergence process, here the adaptive operator of Equation (21) was introduced, then the size of the crossover rate *PC* was calculated and so was the mutation rate *Pm* between individuals, thereby further obtaining the target value. The compiled calculation program was conducted on a PC with the CPUs (AMD Ryzen 71700 Eight-Core Processor, 2.99 GHz); and the important gene was 0 with 2 important gene positions, which were the 1st and 5th genes, respectively. In order to avoid contingency caused by one or two cycles, the program for finding important genes was run repeatedly and the important genes generated after each cycle were counted and sorted, and then through counting the number of occurrences and position numbers, finally it got the important gene positions as the 1st and 5th genes. Then injecting the obtained important genes into new individuals randomly generated, the strength ratio R was obtained during the evolution process. It needs to run multiple times to prevent local convergence that may occur accidentally and fail to reach the global convergence value.

### *4.2. Results Discussion*

The optimisation diagram of final strength ratio R is shown in Figure 7. It suggests that after multiple runs, almost all optimisation curves had fluctuated in the fitness value in the early stage. This is because it is necessary to determine whether the population individual meets the requirements of the layup design after each iteration; if not, eliminate them and introduce new individuals. Therefore, in the continuous genetic evolution process, individuals that did not meet the requirements were eliminated, and new individuals that meet the requirements were added in thereby ensuring that the original design layup criteria were always met.

This value always changed in the optimisation process in the later stage. As indicated by the first, fourth, fifth, and seventh optimisation curves in Figure 7, the best strength ratio appeared at a certain genetic stage, and the genetic generation number was kept short. The "best value" is produced by "gene mutation" and is not the global optimal solution. Simultaneously, the 4th and 7th optimisation curves gradually increase with the number of iterations. Although the optimal strength ratio was found due to the "gene mutation", the final optimisation curve converges locally, which also indicates that the genetic algorithm has the disadvantage of local convergence, and thus, unable to find the global optimal solution smoothly.

**Figure 7.** Ratio of number of iterations of the IGA (*n*) to strength ratio (*R*).

In order to determine whether the curve is reasonable and whether the calculation result meets the initial design requirements, the curve was processed through polynomial data fitting. The seven strength ratios were averaged in the fitting process, and the final fitting result is shown in Figure 8.

**Figure 8.** The fitted curve of the number of iterations of the IGA (*n*) to strength ratio (*R*).

The expression of the red curve in Figure 8 is given below as,

$$R = A\_4 n^4 + A\_3 n^3 + A\_2 n^2 + A\_1 n + A\_0 \tag{32}$$

where

$$\begin{aligned} A\_4 &= -6.5604\varepsilon - 13 & A\_3 &= -7.3372\varepsilon - 10 \\ A\_2 &= -2.4327\varepsilon - 7 & A\_1 &= 3.0651\varepsilon - 5 \\ A\_0 &= 3.8198 \end{aligned} \tag{33}$$

According to statistics, Goodness of fitting *T*<sup>2</sup> is the indicator to evaluate the degree of fit [25]. It is expressed in Equation (34) below.

$$T^2 = \frac{\sum\_{i=1}^{N} \left(\hat{R}\_i - \overline{R}\right)^2}{\sum\_{i=1}^{N} \left(R\_i - \overline{R}\right)^2} \tag{34}$$

where, *R*ˆ *<sup>i</sup>* is the fitted value; *R* is the average value of the initial data; *Ri* is the *i-th* initial data. Putting the average of the strength ratios of the IGA algorithm into Equation (34), it becomes *T*<sup>2</sup> = 0.943. The value range of *T*<sup>2</sup> is 0 to 1. The closer to 1, the better the fitting effect is. In order to further determine the reliability of this improved algorithm, the number of iterations was increased up to 550 generations, and the optimisation iteration was performed three times with their results averaged. At the same time, the calculation result of the fitted curve was substituted for comparison, as shown in Figure 9.

**Figure 9.** The curve of reliability of the IGA.

Figure 9 shows that after 550 generations with the maximum difference ratio of 0.1%, the minimum difference ratio of −0.18%, the maximum difference range of 0.208%, and *<sup>T</sup>*<sup>2</sup> of 0.943 at the same time; hence, it can be considered that the IGA algorithm and the fitted curve are reliable.

It can be seen from Tables 3 and 4 and Figure 7 that when the number of iterations evolves to 52 generations at the earliest stage, the strength ratio R reaches the maximum value of 3.825, suggesting a significant improvement with the increase rate of 48.83%, and further, in the iterative process, the laying angle ±45◦ gradually approaches the outermost layer, which can also effectively improve the impact resistance and stability of the laminate. At the same time, through decoding the optimised layering scheme inversely, it suggests that the number of genes 0 is gradually increasing. This gradual increase also shows that

0 is an important gene, which is consistent with the initial calculation results. Compared with the traditional genetic algorithm, in order to ensure a single variable, the population initially formed by the IGA was retained, and simultaneously ensuring that other initial settings remained unchanged; then optimisation iterations were performed. The calculation results are shown in Figure 10 and Table 5.



**Table 4.** Comparison before and after layup optimisation.


\* The small "s" next to the right bracket denotes symmetric layup hereof.

**Figure 10.** The curve of traditional genetic algorithm number of iterations N-strength ratio, R.


**Table 5.** Calculation results of traditional genetic algorithm.

\* "[9]" denotes there are nine layers with a laying angle of 0◦.

Figure 10 presents that the strength ratio R obtained by the GA algorithm mainly shows an increasing trend at the beginning, and the maximum strength ratio *R* is 8.819. At the same time, Table 5 shows that the laying angles are all 0◦ at this stage, and then the strength ratio *R* fluctuates. This can be attributed to the "gene mutation" due to mutation during the evolution process. From Table 4, it is known that the "gene mutation" has produced Gene 1; however, as indicated by the layup scheme in Table 4, the laying angle does not meet the original design criteria, which led to the phenomenon of "prematurity".

At the same time, through the comparison between Figures 7 and 10 it results that:


### **5. FEM Validation**

In order to validate the proposed IGA, a finite element model of a laminate was developed with a length of 1200 mm, a width of 130 mm, and thickness of 2.25 mm. The laminate has a total of 18 monolayers with a thickness of each layer of 0.125 mm. The Laminate plate was constrained by fixed the x, y, and z directions on all four sides (U1 = U2 = U3 = 0). That is to say, the laminate plate was set as a simply supported beam. Then the plate was meshed with hexagonal shell elements with eight nodes (SC8R) and reduced integration. A uniform load of 0.5 MPa was applied perpendicular to the *x*–*y* plane of the laminate plane. Finally, the finite element analysis was conducted through the commercial software ABAQUS and the simulation results before and after optimisation were obtained, as shown in Figure 11.

As indicated by Figure 11a,b, the maximum stress in the fibre direction before optimisation is 116.0 MPa, and the maximum stress after optimisation is 100.9 MPa. The stress value is greatly reduced, and the reduction ratio reaches 13.02%, indicating that the optimisation results meet the initial requirements. The stresses at the other two directions, namely s22 and s12 changed from 8.57 MPa and 8.80 MPa to 5.464 MPa and 11.11 MPa, respectively. These two directions are not the directions mainly subject to applied forces and the corresponding stresses are also relatively small to s11. At the same time, it can be found from the figures that the maximum stress is mainly concentrated on the two long sides of the laminate. This is because the use of four-sided simple support leads to stress concentration.

**Figure 11.** Stress contour of the laminate before and after optimisation: (**a**) Stress contour before optimisation; (**b**) Stress contour after optimisation.

### **6. Conclusions**

Through the development of an IGA, the laying angle and laying sequence of laminates were optimised and the optimal results were compared with those of traditional genetic algorithms. The reliability of the optimisation was further validated by comparison with the finite element analysis results and the conclusions are drawn as below.

• Based on the genetic algorithm, this paper optimised the layering angle of laminate, and adopted the identification of important genes to determine the important genes and their positions and reduced the query range so that the calculation cost in terms of time was reduced by 21.03%.


**Author Contributions:** Conceptualization, X.M., Z.S. and S.C.; methodology, X.M. and H.L.; software, X.M.; validation, X.M., H.X. and X.X.; formal analysis, X.M.; investigation, X.M. and Z.S.; resources, X.X.; data curation, H.L. and H.X.; writing—original draft preparation, X.M. and Z.S.; writing review and editing, X.M., Z.S. and S.C.; visualization, X.M.; supervision, S.C.; project administration, S.C.; funding acquisition, S.C. All authors have read and agreed to the published version of the manuscript.

**Funding:** This research is financially supported by the Fundamental Research Funds for the Central Universities (No.: JZ2021HGQB0274).

**Institutional Review Board Statement:** Not applicable.

**Informed Consent Statement:** Not applicable.

**Data Availability Statement:** Not applicable.

**Conflicts of Interest:** The authors declare no conflict of interest.

### **Nomenclature**


### **Appendix A**

Considering a composite beam with simply supported ends, it has

$$\begin{cases} \begin{array}{c} \sigma\_{\overline{z}} = 0 \\ \tau\_{yz} = \tau\_{zy} = \tau\_{zx} = \tau\_{xz} = 0 \end{array} \end{cases} \tag{A1}$$

Here, only *σx*, *σy*, *τxy* are not zero. Then according to statics, it obtains the following expressions.

• Balance equation (ignoring body stress)

$$\begin{cases} \frac{\partial \tau\_x}{\partial \mathbf{x}} + \frac{\partial \tau\_{xy}}{\partial y} = 0\\ \frac{\partial \tau\_y}{\partial \mathbf{y}} + \frac{\partial \tau\_{xy}}{\partial \mathbf{x}} = 0 \end{cases} \tag{A2}$$

Let stress ∅ be a function of *x*, *y*, the following relationship holds:

$$
\sigma\_x = \frac{\partial^2 \mathcal{Q}}{\partial y^2}, \sigma\_y = \frac{\partial^2 \mathcal{Q}}{\partial x^2}, \tau\_{xy} = -\frac{\partial^2 \mathcal{Q}}{\partial x \partial y} \tag{A3}
$$

• Geometric compatibility equations

$$\frac{\partial^2 \varepsilon\_x}{\partial y^2} + \frac{\partial^2 \varepsilon\_y}{\partial x^2} = \frac{\partial^2 \gamma\_{xy}}{\partial x \partial y} \tag{A4}$$

• Physical equations

$$S = A^{-1} \tag{A5}$$

$$\begin{cases} \varepsilon\_{\mathcal{X}} = S\_{11}' \sigma\_{\mathcal{X}} + S\_{12}' \sigma\_{\mathcal{Y}} + S\_{16}' \tau\_{\mathcal{X}\mathcal{Y}} \\ \quad \varepsilon\_{\mathcal{Y}} = S\_{12}' \sigma\_{\mathcal{X}} + S\_{22}' \sigma\_{\mathcal{Y}} + S\_{26}' \tau\_{\mathcal{X}\mathcal{Y}} \\ \quad \gamma\_{\mathcal{X}\mathcal{Y}} = S\_{16}' \sigma\_{\mathcal{X}} + S\_{26}' \sigma\_{\mathcal{Y}} + S\_{66}' \tau\_{\mathcal{X}\mathcal{Y}} \end{cases} \tag{A6}$$

• Substituting (A3) into (A6) and then into (A4), it can obtain the equation below,

$$S\_{22}'\frac{\partial^4 \mathcal{Q}}{\partial x^4} - 2S\_{26}'\frac{\partial^4 \mathcal{Q}}{\partial x^3 \partial y} + \left(2S\_{12}' + S\_{66}'\right)\frac{\partial^4 \mathcal{Q}}{\partial x^2 y^2} - 2S\_{16}'\frac{\partial^4 \mathcal{Q}}{\partial x \partial y^3} + S\_{11}'\frac{\partial^4 \mathcal{Q}}{\partial y^4} = 0 \tag{A7}$$

The expressions of *σx*, *σy*, *τxy* from (A7) as below

<sup>1</sup> Since *σ<sup>y</sup>* is caused by the constant shear force q, and its value does not depend on *x*, it is a function of *y*. Therefore, we can set *σ<sup>y</sup>* = *f*(*y*)

$$\frac{\partial^2 \mathcal{Q}}{\partial \mathbf{x}^2} = \sigma\_y = f(y) \tag{A8}$$

Through integral formula (A8), it can produce:

$$\mathcal{Q} = \frac{1}{2}\mathbf{x}^2 f(y) + \mathbf{x}f\_1(y) + f\_2(y) \tag{A9}$$

<sup>2</sup> Therefore, substituting (A9) into the compatible equation, each item of (A7) can be obtained as follows: *∂*4∅

$$\begin{aligned} \frac{\partial^4 \mathcal{D}}{\partial x^4} &= 0, \frac{\partial^4 \mathcal{D}}{\partial x^3 \partial y} = 0 \\\\ \frac{\partial^4 \mathcal{D}}{\partial x^2 y^2} &= \frac{d^2 f(y)}{dy^2} \end{aligned} \tag{A10}$$

$$\begin{aligned} \frac{\partial^4 \mathcal{D}}{\partial x \partial y^3} &= x \frac{d^3 f(y)}{dy^3} + \frac{d^3 f\_1(y)}{dy^3} \\\\ \frac{\partial^4 \mathcal{D}}{\partial y^4} &= x^2 \frac{d^4 f(y)}{2 dy^4} + x \frac{d^4 f\_1(y)}{dy^4} + \frac{d^4 f\_2(y)}{dy^4} \end{aligned} \tag{A10}$$

Substituting into the compatible Equation (A4), it becomes

*∂*4∅

$$\left(2S\_{12}' + S\_{66}'\right)\frac{d^2f(y)}{dy^2} - 2S\_{16}'\left(\mathbf{x}\frac{d^3f(y)}{dy^3} + \frac{d^3f\_1(y)}{dy^3}\right) + S\_{11}'\left(\mathbf{x}^2\frac{d^4f(y)}{2dy^4} + \mathbf{x}\frac{d^4f\_1(y)}{dy^4} + \frac{d^4f\_2(y)}{dy^4}\right) = 0\tag{A11}$$

If it holds for any x of the above formulas, then it gets

$$\frac{S\_{11}'}{2} \frac{d^4 f(y)}{dy^4} = 0 \tag{A12a}$$

$$S\_{11}' \frac{d^4 f\_1(y)}{dy^4} - 2S\_{16}' \frac{d^3 f(y)}{dy^3} = 0 \tag{A12b}$$

$$2\left(2S\_{12}' + S\_{66}'\right)\frac{d^2f(y)}{dy^2} - 2S\_{16}'\frac{d^3f\_1(y)}{dy^3} + S\_{11}'\frac{d^4f\_2(y)}{dy^4} = 0\tag{A12c}$$

If (a) holds, then

$$f(y) = \frac{A\_1 y^3}{6} + \frac{A\_2 y^2}{2} + A\_3 y + A\_4 \tag{A13}$$

If (b) holds, then

$$f\_1(y) = \frac{A\_1}{12} \frac{S\_{16}'}{S\_{11}'} y^4 + \frac{B\_1}{6} y^3 + \frac{B\_2}{2} y^2 + B\_3 y + B\_4 \tag{A14}$$

Because *B*<sup>4</sup> is not used in next equation, it can be omitted. If (c) holds, then

$$f\_{\Omega}(y) = \left[\frac{1}{30}A\_1\left(\frac{S\_{16}'}{S\_{11}'}\right)^2 - \frac{A\_1}{120}\frac{\left(2S\_{12}' + S\_{66}'\right)}{S\_{11}'}\right]y^5 + \left[\frac{1}{12}B\_1\frac{S\_{16}'}{S\_{11}'} - \frac{A\_2}{24}\frac{\left(2S\_{12}' + S\_{66}'\right)}{S\_{11}'}\right]y^4 + \frac{C\_1}{6}y^3 + \frac{C\_2}{2}y^2 + C\_3y + C\_4 \tag{A15}$$

Because (*C*<sup>3</sup> *y* + *C*4) is not used in next equation, it can be omitted.

Finally, through combining the Equations (A13)–(A15), it obtains stress function as below,

$$\begin{aligned} \mathcal{Q} &= \frac{1}{2} \mathbf{x}^2 f(\mathbf{y}) + x f\_1(\mathbf{y}) + f\_2(\mathbf{y}) = \frac{1}{2} \mathbf{x}^2 \left( \frac{A\_1 \mathbf{y}^3}{6} + \frac{A\_2 \mathbf{y}^2}{2} + A\_3 \mathbf{y} + A\_4 \right) + \\ &\ge \left( \frac{A\_1}{12} \frac{S\_{16}'}{S\_{11}'} \mathbf{y}^4 + \frac{B\_1}{6} \mathbf{y}^3 + \frac{B\_2}{2} \mathbf{y}^2 + B\_3 \mathbf{y} \right) + \left[ \frac{1}{30} A\_1 \left( \frac{S\_{16}'}{S\_{11}'} \right)^2 - \frac{A\_1}{120} \frac{\left( 2S\_{12}' + S\_{66}' \right)}{S\_{11}'} \right] \mathbf{y}^5 + \\ &\ge \left[ \frac{1}{12} B\_1 \frac{S\_{16}'}{S\_{11}'} - \frac{A\_2}{24} \frac{\left( 2S\_{12}' + S\_{66}' \right)}{S\_{11}'} \right] \mathbf{y}^4 + \frac{C\_1}{6} \mathbf{y}^3 + \frac{C\_2}{2} \mathbf{y}^2 \end{aligned} \tag{A16}$$

Through the expression of ∅, here it obtains its stress component σ*x*, σ*y*, τ*xy*

$$
\sigma\_x = \frac{\partial^2 \sigma}{\partial y^2} = \left[\frac{2}{3} A\_1 \left(\frac{S\_{1b}'}{S\_{11}'}\right)^2 - \frac{A\_1}{6} \frac{\left(2S\_{12}' + S\_{6b}'\right)}{S\_{11}'}\right] y^3 + \left[A\_1 \frac{S\_{1b}'}{S\_{11}'} \mathbf{x} + B\_1 \frac{S\_{1b}'}{S\_{11}'} - \frac{A\_2}{2} \frac{\left(2S\_{12}' + S\_{6b}'\right)}{S\_{11}'}\right] y^2 + \dotsb \tag{A17}
$$

$$
\left[\frac{A\_1}{2} \mathbf{x}^2 + B\_1 \mathbf{x} + C\_1\right] y + \left[\frac{A\_2}{2} \mathbf{x}^2 + B\_2 \mathbf{x} + C\_2\right]
$$

$$
\sigma\_y = f(y) = \frac{A\_1 y^3}{6} + \frac{A\_2 y^2}{2} + A\_3 y + A\_4
$$

$$
\tau\_{xy} = -\left(\frac{\partial^2 \sigma}{\partial x \partial y}\right) = -\left[\mathbf{x}\left(\frac{A\_1}{2} y^2 + A\_2 y + A\_3\right) + \frac{A\_1}{3} \frac{S\_{1b}'}{S\_{11}'} y^3 + \frac{B\_1}{2} y^2 + B\_2 y + B\_3\right]
$$

1. When the upper and lower boundaries *y* = ±h/2, *τxy* = 0 (this holds for any *x*), through the above equation, the relational expressions of *A*2, *A*3, *B*1, *B*2, *B*3, and *A*<sup>1</sup> can be derived.

$$\begin{cases} \begin{array}{c} \text{x} \left( \frac{A\_1}{2} \frac{h^2}{4} + A\_2 \frac{h}{2} + A\_3 \right) + \frac{A\_1}{3} \frac{S\_{16}'}{S\_{11}} \frac{h^3}{8} + \frac{B\_1}{2} \frac{h^2}{4} + B\_2 \frac{h}{2} + B\_3 = 0 \\\\ \text{x} \left( \frac{A\_1}{2} \frac{h^2}{4} - A\_2 \frac{h}{2} + A\_3 \right) - \frac{A\_1}{3} \frac{S\_{16}'}{S\_{11}} \frac{h^3}{8} + \frac{B\_1}{2} \frac{h^2}{4} - B\_2 \frac{h}{2} + B\_3 = 0 \end{array} \end{cases} \tag{A18}$$

Because *τxy* = 0 holds for any *x*, the first items of the Equation (A18) can produce:

$$\begin{cases} \frac{A\_1}{2} \frac{h^2}{4} + A\_2 \frac{h}{2} + A\_3 = 0\\ \frac{A\_1}{2} \frac{h^2}{4} - A\_2 \frac{h}{2} + A\_3 = 0 \end{cases} \tag{A19}$$

Solving the Equation (A19), *<sup>A</sup>*<sup>2</sup> <sup>=</sup> <sup>0</sup> *<sup>A</sup>*<sup>3</sup> <sup>=</sup> <sup>−</sup>*h*<sup>2</sup> <sup>8</sup> *A*<sup>1</sup> can be obtained. The other items of the Equation (A18) can produce:

$$\begin{cases} \begin{array}{c} \frac{A\_1}{3} \frac{S\_{16}'}{S\_{11}'} \frac{h^3}{8} + \frac{B\_1}{2} \frac{h^2}{4} + B\_2 \frac{h}{2} + B\_3 = 0 \\\\ -\frac{A\_1}{3} \frac{S\_{16}'}{S\_{11}'} \frac{h^3}{8} + \frac{B\_1}{2} \frac{h^2}{4} - B\_2 \frac{h}{2} + B\_3 = 0 \end{array} \end{cases} \tag{A20}$$

Solving the Equation (A20), *B*<sup>3</sup> *B*<sup>2</sup> can be obtained

$$\begin{cases} \begin{array}{c} B\_3 = -\frac{h^2}{8} B\_1 \\\\ B\_2 = -\frac{h^2}{12} \frac{S\_{16}'}{S\_{11}'} A\_1 = \frac{S\_{16}'}{S\_{11}'} \frac{q}{h} \end{array} \tag{A21}$$

When the upper and lower boundaries the *y* = ±h/2, the equations of *y* = +*<sup>h</sup>* <sup>2</sup> , *σ<sup>y</sup>* = *q <sup>y</sup>* <sup>=</sup> <sup>−</sup>*<sup>h</sup>* <sup>2</sup> , *<sup>σ</sup><sup>y</sup>* <sup>=</sup> <sup>0</sup> ,*σ<sup>y</sup>* <sup>=</sup> *<sup>f</sup>*(*y*) <sup>=</sup> *<sup>A</sup>*1*y*<sup>3</sup> <sup>6</sup> <sup>+</sup> *<sup>A</sup>*2*y*<sup>2</sup> <sup>2</sup> + *A*3*y* + *A*<sup>4</sup> hold. Then the relational expression of *A*2, *A*3, *A*<sup>4</sup> and *A*<sup>1</sup> can be derived. The derivation process is shown as below

$$\begin{cases} \quad \frac{A\_1}{6} \frac{h^3}{8} + \frac{A\_2}{2} \frac{h^2}{4} + A\_3 \frac{h}{2} + A\_4 = \mathbf{q} \\\ \quad - \frac{A\_1}{6} \frac{h^3}{8} + \frac{A\_2}{2} \frac{h^2}{4} - A\_3 \frac{h}{2} + A\_4 = 0 \end{cases} \tag{A22}$$

Solving the Equation (A22), *A*<sup>3</sup> *A*<sup>4</sup> can be obtained

$$\begin{cases} \quad A\_3 = \frac{q}{l} - \frac{h^2}{24} A\_1 \\ \quad A\_4 = \frac{q}{2} - \frac{h^2}{8} A\_2 = \frac{q}{2} \end{cases} \tag{A23}$$

From the Equations (A19) and (A23), *A*<sup>1</sup> *A*<sup>3</sup> can be obtained

$$\begin{cases} \ A\_1 = -12 \frac{q}{\hbar^3} \\ \ A\_3 = \frac{3q}{2\hbar} \end{cases} \tag{A24}$$

2. With the boundary of *x* = *l*, the shear force Q = −*ql*,

$$\int\_{-\frac{b}{2}}^{\frac{b}{2}} \left\{-\left[x\left(\frac{A\_1}{2}y^2 + A\_2y + A\_3\right) + \frac{A\_1}{3}\frac{S\_{16}'}{S\_{11}'}y^3 + \frac{B\_1}{2}y^2 + B\_2y + B\_3\right]\right\} dy = -ql \tag{A25}$$

Then the relational expression of *A*2, *A*3, *B*1, *B*2, *B*<sup>3</sup> and *A*<sup>1</sup> can be derived. The derivation process is shown, as below

$$\left. \left( -\frac{A\_1}{2} l \frac{1}{3} y^3 - A\_2 l \frac{1}{2} y^2 - A\_3 l y - \frac{A\_1}{3} \frac{S\_{16}'}{S\_{11}'} \frac{1}{4} y^4 - \frac{B\_1}{2} \frac{1}{3} y^3 - B\_2 \frac{1}{2} y^2 - B\_3 y \right) \right|\_{-\frac{h}{2}}^{\frac{h}{2}} = -ql \quad \text{(A26)}$$

Integration of the above formula, it produces

$$lA\_1 + B\_1 = \frac{24ql - 24hlA\_3 - 24lhB\_3}{h^3} \tag{A27}$$

Because *<sup>A</sup>*<sup>3</sup> = <sup>3</sup>*<sup>q</sup>* <sup>2</sup>*<sup>h</sup>* as shown in Equation (A24) and *<sup>B</sup>*<sup>3</sup> <sup>=</sup> <sup>−</sup>*h*<sup>2</sup> <sup>8</sup> *B*<sup>1</sup> in Equation (A21), *B*<sup>1</sup> *and B*<sup>3</sup> can be obtained here as below.

$$\begin{cases} \begin{array}{c} B\_1 = 0\\ B\_3 = -\frac{h^2}{8} B\_1 = 0 \end{array} \tag{A28}$$

3. Since the *σ<sup>x</sup>* at the x-axis direction is zero; the moment is also zero. Hence, it gives

$$\int \begin{array}{cc} \int\_{-\frac{l}{2}}^{\frac{l}{2}} (\sigma\_x)\_{x=l} dy = 0 \end{array} \tag{A29a}$$

$$\int \!\_{-\frac{b}{2}}^{\frac{b}{2}} (\sigma\_{\mathbf{x}})\_{\mathbf{x}=l} y dy = 0 \tag{A29b}$$

Through this Equation (A29), the relational expression of *A*2, *B*1, *B*2, *B*3, *C*1, *C*2, and *A*<sup>1</sup> can be derived. From (a) of (A29), it gets *C*<sup>2</sup> = 0; and from (b) of (A29), it gets *C*<sup>1</sup> = 6 5 *q h S* 16 *S* 11 2 <sup>+</sup> <sup>6</sup> *ql*<sup>2</sup> *<sup>h</sup>*<sup>3</sup> <sup>−</sup> <sup>3</sup> 10 *q h* (2*S* 12+*S* 66) *S* 11 Finally, here the coefficients, *A*1, *A*2, *A*3, *A*4, *B*1, *B*2, *B*3, *C*1, *C*<sup>2</sup> are all been obtained, the following equation holds:

$$\begin{split} \mathcal{Q} &= \frac{1}{2} \mathbf{x}^2 \left( -\frac{2q}{\hbar} y^3 + \frac{3q}{2\hbar} y + \frac{q}{2} \right) + \mathbf{x} \left( -\frac{q}{\hbar} \frac{\mathbf{s}\_{1\mathbf{d}}'}{\mathbf{s}\_{1\mathbf{1}}} y^4 + \frac{\mathbf{s}\_{1\mathbf{q}}' q}{2\boldsymbol{S}\_{1\mathbf{1}} \hbar} y^2 \right) - \frac{2q}{\hbar^3} \left[ \left( \frac{\mathbf{s}\_{1\mathbf{d}}'}{\mathbf{s}\_{1\mathbf{1}}'} \right)^2 - \frac{\left( 2\mathbf{s}\_{1\mathbf{2}}' + \mathbf{s}\_{\mathbf{d}}' \right)}{4\mathbf{s}\_{1\mathbf{1}}'} \right] y^5 + \\ &\frac{1}{6} \left[ \frac{6q}{\hbar^3} \left( \frac{\mathbf{s}\_{1\mathbf{d}}'}{\mathbf{s}\_{1\mathbf{1}}'} \right)^2 + \frac{6q^2}{\hbar^3} - \frac{3q}{10\hbar} \frac{2\mathbf{s}\_{1\mathbf{2}}' + \mathbf{s}\_{\mathbf{d}}'}{\mathbf{s}\_{1\mathbf{1}}'} \right] y^3 \end{split} \tag{A30}$$

Finally, *σx*, *σy*, *τxy* can be obtained through Equation (A30) as below.

$$
\sigma\_x = \frac{-12q}{h^3} \left\{ \left[ \frac{2}{3} \left( \frac{S\_{\rm fb}'}{S\_{\rm f1}'} \right)^2 - \frac{2S\_{\rm f2}' + S\_{\rm f6}'}{6S\_{\rm f1}'} \right] y^3 + \frac{S\_{\rm f}'}{S\_{\rm f1}'} \text{xy}^2 + \frac{1}{2} \mathbf{x}^2 y \right\} + 
$$

$$
\left[ \frac{6q}{5h} \left( \frac{S\_{\rm fb}'}{S\_{\rm f1}'} \right)^2 + \frac{6ql^2}{h^3} - \frac{3q}{10h} \frac{2S\_{\rm f2}' + S\_{\rm f6}'}{S\_{\rm f1}'} \right] y + \frac{S\_{\rm f6}'}{S\_{\rm f1}'} \frac{q}{h} \mathbf{x} \tag{A31}
$$

$$
\sigma\_y = -\frac{2q}{h^3} y^3 + \frac{3q}{2h} y + \frac{q}{2}
$$

$$
\tau\_{xy} = -\left[ \mathbf{x} \left( -\frac{6q}{h} y^2 + \frac{3q}{2h} \right) - \frac{4q}{h^3} \frac{S\_{\rm f6}'}{S\_{\rm f1}'} y^3 + \frac{S\_{\rm f6}'}{S\_{\rm f1}'} \frac{q}{2} y \right]
$$

### **References**

