*Article* **Fundamental Frequency Optimization of Variable Angle Tow Laminates with Embedded Gap Defects**

**João Carvalho 1, Abdolrasoul Sohouli 2,\* and Afzal Suleman 1,2**


**Abstract:** Variable stiffness composite laminates can improve the structural performance of composite structures by expanding the design space. This work explores the application of variable stiffness laminated composite structures to maximize the fundamental frequency by optimizing the tow angle. To this end, an optimization framework is developed to design the fiber angle for each layer based on the maximization of the fundamental frequency. It is assumed that the design process includes the manufacturing constraints encountered in the automated fiber placement process and a linear fiber angle variation. The current study improves existing results by considering embedded gap defects within the optimization framework. The plates are assumed symmetric, with clamped and simply supported boundary conditions. The optimal results and a comparison between the non-steered and steered plates with and without gaps are presented. Results show that, although gaps deteriorate the structural performance, fiber steering can still lead to an increase in the fundamental frequency depending on the plate's geometry and boundary conditions.

**Keywords:** variable stiffness composite; variable angle tow; automated fiber placement; Defect Layer Method; natural frequency; genetic algorithm optimization

### **1. Introduction**

Nowadays, the use of composites materials is increasing in a wide range of sectors such as the aerospace, automotive, naval and others. This is the main motivation to develop innovative and cost effective composite manufacturing techniques which allow the production of composite laminates meeting certain requirements in a more effective way. In particular, the development of the Automated Fiber Placement (AFP) technique opens the possibility of manufacturing Variable Stiffness Composite Laminates (VSCLs), which allow the stiffness to vary spatially. VSCLs include, for instance, those made with variable fiber spacing (VFS) and those in which the fibers are deposited following curvilinear paths, denominated as variable-angle tow (VAT) laminates. They allow new design possibilities by enlarging the design space when compared to conventional composite laminates. Hence, several authors have been studying them with the objective of improving composite laminates characteristics. For instance, the studies presented in [1–3] show the influence of these kind of laminates in the buckling behaviour of laminated panels.

Although most studies for vibration analysis are performed on composite laminated plates with constant stiffness, some authors also considered the use of VAT laminates. Most of these studies aimed to study the influence of the curvature of the fibers on the fundamental frequency achieved by a plate. Honda et al. [4] considered splines to represent the fiber paths of curvilinear fibers. Honda et al. [5] developed a multi-objective approach in order to maximize both fundamental frequencies and in-plane strengths of VAT plates, while Pereira et al. [6] performed a multi-objective optimization considering the fundamental frequency and the first mode specific damping capacity (SDC). Serhat

**Citation:** Carvalho, J.; Sohouli, A.; Suleman, A. Fundamental Frequency Optimization of Variable Angle Tow Laminates with Embedded Gap Defects. *J. Compos. Sci.* **2022**, *6*, 64. https://doi.org/10.3390/jcs6020064

Academic Editors: Stelios K. Georgantzinos and Francesco Tornabene

Received: 31 December 2021 Accepted: 17 February 2022 Published: 20 February 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/).

et al. [7], also maximized the fundamental frequency of variable stiffness panels but proposing a more efficient method to design variable stiffness composite panels using lamination parameters. Although the steering radius was not directily constrained during the optimization process, the design space is controlled to ensure a smooth changes in the layer angles. The design methodology efficiency developed in [7] is improved by Serhat et al. [8] study by combining the lamination parameters formulation with a spectral Tchebychev modeling approach. In this study, the fundamental frequencies of the variable stiffness designs were maximized for different boundary conditions and aspect ratios. It was verified that optimal variable stiffness designs provide higher frequencies (up to 10%) when compared to constant stiffness ones. More recently, Farsadi et al. [9] included plates with different aspect ratios, boundary conditions and skew angles. The fundamental frequency was also studied and optimized together with other composite characteristics.

The most common manufacturing constraint considered when designing curvilinear composite laminates is the maximum curvature constraint [10–14]. Akhavan and Ribeiro [10], for example, considered this constraint in the study of the natural frequencies and mode shapes of VAT with a linear curvilinear fiber path. They used a p-version finite element and a Third-order Shear Deformation Theory (TSDT). It was verified that the use of curvilinear fibers can be advantageous to adjust frequencies and mode shapes. Demir et al. [14] developed a lamination parameter optimization algorithm using a penalty parameter to enforce a maximum value for the fiber curvatures. The proposed methodology is applied to design for the minimum compliance in-plane and out-of-plane problems. In a recent study, Rashed and Demir [15] maximized the fundamental frequency of variable stiffness plates. Manufacturing constraints were considered and the method employed uses lamination parameters as the design variables. The optimum lamination parameters distribution found is in good agreement with other studies present in the literature.

The previous studies mentioned assumed that the VAT laminated plate is ideal, which means that the models used to obtain the natural frequencies and mode shapes did not take into account the manufacturing process, in particular, the AFP and the respective manufacturing induced imperfections. As a result of that, some models were developed in order to take defects into consideration and better predict the structural behaviour of VAT laminates. Blom [16] developed a methodology where each ply of the laminate is modeled including gaps and overlaps generated by tow drops. Afterwards, a 2-D mesh is generated on top of each ply and the centroids of all the elements in all the plies are checked for the tow drop region. If it is the tow drop region considered as a gap, then at that position resin properties are assigned to the respective element stack. In case of an overlap, the thickness assigned to the element stack is twice the thickness. It is important to emphasize that a tow drop defect is only identified if the respective element centroid is located over it. In consequence, a very refined mesh is required to have a good precision which is associated to a higher computational cost. A similar strategy was developed by Fayazbakhsh [17,18]. However, each element properties are assigned in accordance with its volume fraction of gaps and overlaps using micro-mechanics. In this way, it is not required such a refined mesh and computational cost is saved. Several authors considered this modeling methodology to include the effect of both gaps and overlaps [19–21]. Marouene et al. [19] study considered this last methodology to model variable stiffness plates with the goal of maximizing the in-plane stiffness and the buckling load. The numerical results obtained were then validated with experimental studies and a good accordance between them was found. Akbarzadeh et al. [20] study was focused on evaluating the impact that AFP manufacturing defects, namely, gaps and overlaps, have on the laminate properties of variable stiffness plates. Both static bending, buckling and free vibration have been studied. It also examined the role of shear deformation on those plates. On the one hand, for very thin plates, all equivalent single layer theories used provided similar results. On the other hand, the results of moderately thick plates showed some discrepancies, especially between those obtained using the Classical Laminate Theory (CLPT) and the Third Order Shear Deformation Theory (TSDT). It was verified an increase of the buckling load and

fundamental frequency when compared to the steered plate with no defects, if overlaps were considered. When gaps were the plate manufacturing defect, a decrease in both mentioned parameters was verified.

In addition to the studies presented in the previous paragraph, other novel methods to include manufacturing defects were developed in [22,23]. Brooks and Martins [22] treated tow paths as the streamlines of a unit-vector field. In the developed mathematical formulation the gap/overlap formation is related with the divergence of a 2D vector field. A positive divergence indicates a gap formation and a negative value the occurrence of an overlap. Mishra et al. [23] explored the effect of gaps on the stiffness and buckling load of variable stiffness panels. The developed methodology used to account the effect of gaps showed to be significantly more efficient than the Defect Layer Method described in [17,18].

Other authors considered instead more detailed and comprehensive 3D analysis [24]. Marrouzé et al. [25], Heinecke et al. [26] and Willie et al. [27] approach models gaps and overlaps in both laminate, lamina and micro structure level. Li et al. [28] developed 3D meshing tools to automatically generate ply by ply models with gaps and overlaps. It automatically inserts intra and inter-ply cohesive elements in order to capture the influence of splitting and delamination. As a result, out-of-plane waviness and ply thickness variations are automatically modeled. These detailed 3D finite element models can capture effects like load redistribution and stress localization. Therefore, these type of analysis could be justifiable if higher levels of detail are necessary.

Regarding experimental studies, in Antunes et. al. study [29], which is the a continuation of Rodrigues et al. study [30], experimental modal analysis were performed on a rectangular plate considering various boundary conditions including one edge clamped and the remaining edges free (CFFF), two opposite edges clamped and the other two free (CFCF) and all the edges free (FFFF). Moreover, modal damping ratios were also identified in addiction to both natural frequencies and mode shapes. In this article, the theoretical model is based on two p-version type finite element models, one based on the Classical Plate Theory (CPT) and another based on the First-order Shear Deformation Theory (FSDT). Similar results were obtained using the two models, which can be justified by the use of a very thin plate. The authors reported that the discrepancies found between the experimental and numerical results were due to the occurrence of gaps and overlaps in the plate, which the FEM model did not take into account.

In this work, an optimization framework is developed with the objective of designing composite plates with curvilinear fibers for maximum fundamental frequency. Square and rectangular plates under fully clamped and fully simply supported boundary conditions are optimized using a Genetic Algorithm. The optimization is done first considering only straight fiber and then using curvilinear fibers with the purpose of evaluating the difference between the maximum frequencies obtained. The optimization of curvilinear fiber plates is performed considering both the absence and different values of the maximum curvature constraint of the AFP machine, thereby allowing to assess its effect on both maximum frequency and optimal designs. Moreover, the optimization of curvilinear fiber plates is performed considering the absence and the presence of embedded gap defects to evaluate its effect on both maximum frequencies and optimal designs. Gap imperfections are taken into consideration using the Defect Layer Method [17]. The authors consider that the consideration of defects within the optimization framework is the principal contribution of this work.

The remainder of this article is organized as follows: first, the variable stiffness plates considered are described in Section 2. This Section also includes the methodology for modeling the VATs and for including the effect of gaps in the FE model. The optimization problem is formulated in Section 3. Section 4 presents the developed framework used to solve the optimization problem. Section 5 discusses and analyses the results obtained. Finally, Section 6 presents the main conclusions drawn of this work.

### **2. Materials and Methods**

### *2.1. Variable Angle Tow Laminates*

In the AFP process, a band of fibers, denominated course, is composed by individual units named tows. The courses are deposited following a reference fiber path which can be curvilinear. Different curvilinear fiber paths have been studied by various authors. In this work, the fiber path presented in the following section is the one considered to solve the optimization problem.

### 2.1.1. Fiber Path Definition

The reference course fiber path considered was first introduced in [31] in 1993. In this formulation, the fiber angle (*θ*), measured with respect to the *x*∗ axis, varies linearly along the axis *x*∗. The respective fiber angle variation is presented in Equation (1). The fiber path trajectory is presented in Equation (2) [32,33]. The variable *T*0, represents the value of *θ* in the middle of the plate, while the variable *T*<sup>1</sup> represents the value of *θ* when *x*<sup>∗</sup> = ±*a*/2. The variable *φ* represents the orientation of *x*\* relative to the global axis *x*. This fiber path is represented in Figure 1a.

In this work, *φ* is considered as zero. Thus, both axis *x* and *x*<sup>∗</sup> coincide and *T*<sup>1</sup> represents the fiber angle at both edges in *x* of a plate. The middle and edge fiber angles are written as < *T*0, *T*<sup>1</sup> >. As an example, the fiber path correspondent to < 60◦, 15◦ > is presented in Figure 1b for a plate with length *a* in the *x* axis.

$$
\theta(\mathbf{x}^\*) = \phi + (T\_1 - T\_0) 2 \frac{|\mathbf{x}^\*|}{a} + T\_0 \tag{1}
$$

$$y^\*(\mathbf{x}^\*) = \begin{cases} \frac{a}{2(T\_1 - T\_0)} [-\ln(\cos(T\_0)) + \\ \ln\left(\cos\left(T\_0 + \frac{2(T\_0 - T\_1)}{a} \mathbf{x}^\*\right)\right)], -\frac{a}{2} \le \mathbf{x}^\* \le 0 \\ \frac{a}{2(T\_0 - T\_1)} [-\ln(\cos(T\_0)) + \\ \ln\left(\cos\left(T\_0 + \frac{2(T\_1 - T\_0)}{a} \mathbf{x}^\*\right)\right)], 0 \le \mathbf{x}^\* \le \frac{a}{2} \end{cases} \tag{2}$$

**Figure 1.** Reference fiber path: (**a**) Reference fiber path, with *φ* = 0◦; (**b**) Reference fiber path < *T*0, *T*<sup>1</sup> >=< 60◦, 15◦>, with *φ* = 0◦.

### 2.1.2. Induced Defects

During the deposition of the course, the AFP machine's head is perpendicular to the local fiber angle. Thus, each point along the AFP head has the same orientation as the one corresponding to the reference path. Therefore, the points that form each tow of the course can be calculated using Equation (3). In this equation, *nt* represents the number of tows in a course, *tw* is the tow width and *i* is the tow index which range increases by one from −*nt*/2 to *nt*/2. Moreover, the superscript *ref* is used to denote the points in the reference path.

$$\begin{aligned} \mathbf{x} &= \mathbf{x}^{r\varepsilon f} - i \cdot t\_{\text{tr}} \sin \left( \theta^{r\varepsilon f} \right) \\ \mathbf{y} &= \mathbf{y}^{r\varepsilon f} + i \cdot t\_{\text{tr}} \cos \left( \theta^{r\varepsilon f} \right) \end{aligned} \tag{3}$$

The fiber paths of the adjacent courses are obtained by moving the centerline fiber path in the *y* direction with the shift distance, which leads to the formation of manufacturing defects within the laminate. The shift distance (*ds*) considered is the minimum vertical distance between the top and bottom boundaries of the reference course and it is calculated as presented in Equation (4). The points of the shifted course (*ys*, *xs*) are obtained as presented in Equation (5).

$$d\_s = \min(\frac{n\_t t\_w}{\cos\left(\theta^{ref}\right)})\tag{4}$$

$$\begin{aligned} x\_s &= x^{ref} - i \cdot t\_w \sin\left(\theta^{ref}\right) \\ y\_s &= y^{ref} + i \cdot t\_w \cos\left(\theta^{ref}\right) + d\_s \end{aligned} \tag{5}$$

The deposition of both reference and adjacent courses is presented in Figure 2a considering a course with 8 tows. The reference course and its reference curve are represented in black and green, respectively. The adjacent course and its reference curve are represented in blue and red, respectively. In Figure 2a the intersection points between the reference and adjacent course are also visible. These intersections lead to the formation of defects which can be either gaps or overlaps depending on the way the tows are cut. If the entire course is cut or restarted when there is an intersection with the adjacent course, a large gap or overlap area is formed. So, in order to solve this problem and minimize the defect area, the tow cut and restart capability of the AFP machine is used. It allows to individually control each tow, to cut and restart its deposition and leads to a reduced defect area. In this work, an one-sided cut strategy is employed, which means that the tows are only cut on one side of each course. It is also used the zero percent coverage parameter in the tow drop strategy, which implies that the cut is performed as soon as one edge of a tow intersects the boundary of the adjacent course, also known as complete gap strategy or 0% coverage parameter. The gap formation which results from this cut strategy is presented in Figure 2b, where the gap set is represented in a gray colour.

**Figure 2.** (**a**) Tow cut and restart positions; (**b**) Gap pattern between consecutive courses.

The distance between two consecutive gap sets is the shift distance (Equation (4). So, the gap set is translated vertically with the shift distance across a ply. The gap pattern over an entire ply is represented in Figure 3a for a square plate, but considering that the AFP machine deposits 32 tows, each one with a width of 3.175 mm. The gap distribution is then intersected with the ply mesh in order to calculate the gap area fraction (*Am*) within each element as represented in Figure 3b. Each element equivalent properties are then calculated using a 'modified' rule of mixtures (Equations (6)–(10)). In these equations the subscript *c* represents the non-defective composite and the subscript *m* the matrix or gaps. The composite area fraction, *Ac*, and the gap area fraction, *Am*, are related by *Am* = 1 − *Ac*. It is important to note that both area fractions are equivalent to volume fractions because it is assumed a constant thickness across the ply, with or without gaps.

$$E\_1 = A\_c E\_{1c} + A\_m E\_m \tag{6}$$

$$E\_2 = \frac{E\_{2\mathcal{L}} E\_m}{A\_m E\_{2\mathcal{L}} + A\_{\mathcal{L}} E\_m} \tag{7}$$

$$G\_{12} = \frac{G\_{12c} G\_m}{A\_m G\_{12c} + A\_c G\_m} \tag{8}$$

$$G\_{23} = A\_c G\_{23c} + A\_m G\_m \tag{9}$$

$$\mathbf{v}\_{12} = A\_{\mathcal{E}} \mathbf{v}\_{12\mathcal{E}} + A\_{\mathcal{W}} \mathbf{v}\_{\mathcal{W}} \tag{10}$$

**Figure 3.** (**a**) Gap distribution on a plate (**b**) Gap area fraction in each element.

### 2.1.3. Manufacturing Constraints

The maximum curvature constraint or inverse, minimum curvature radius, depends on the AFP machine and limits the maximum value of curvature of the reference fiber path of each course. Typical curvature constraints are presented in Table 1. This constraint limits the design space by limiting the fiber angle distribution that each lamina possesses. This constraint is dependent on the AFP machine and must be evaluated before the production stage in order to ensure that the composite can be manufactured. For that reason, it is most often used as an optimization constraint when VAT laminates are optimized to maximize or minimize determined characteristics.


**Table 1.** Curvature Constraint, adapted from [34].

The curvature equation was analytically obtained using the fiber reference path presented in Equation (1) and also the definition of curvature of a function with a a single variable presented in Equation (11).

$$K = \frac{\frac{d^2f(x)}{dx^2}}{(1 + (\frac{df(x)}{dx})^2)^{3/2}} \tag{11}$$

*f*(*x*) was substituted in Equation (11) by *y*∗(*x*∗) of the reference fiber path (Equation (2)) treated here as *y*(*x*), because *φ* = 0. The first and second derivatives of *y* with respect to *x* are calculated in the following equations:

$$\frac{dy}{dx} = \tan(\theta) \tag{12}$$

$$\frac{d^2y}{dx^2} = \frac{d\tan\theta}{dx} = \frac{1}{\cos^2\theta} \frac{2(T\_1 - T\_0)}{a} \tag{13}$$

The trigonometric relation presented in Equation (14) is also considered:

$$1 + \tan^2 \theta = \frac{1}{\cos^2 \theta} \tag{14}$$

After substituting Equations (12)–(14) in Equation (11), the curvature equation for the reference fiber path considered is obtained (Equation (15)). This last equation was also considered in references [10,22].

$$K = 2\frac{T\_1 - T\_0}{a} \cos\left( (T\_1 - T\_0)2\frac{|\mathbf{x}^\*|}{a} + T\_0 \right) \tag{15}$$

### **3. Optimization**

The optimization problem that this work intends to solve is defined in Equation (16). The goal is to maximize the fundamental frequency, treated here as *f1*, of a laminated plate. It is assumed that the laminate is symmetric and that the total number of plies (*Np*) is equal to eight. The problem is solved considering different cases where the vector of design variables (**x**) depends on whether only rectilinear or curvilinear fibers (whose orientation is presented in Equation (1) are considered. The cases that only consider rectilinear fibers are treated as non-steered cases and the cases that allow curvilinear fibers are treated as steered cases. In the cases where only rectilinear fibers are considered (*T*<sup>0</sup> = *T*1), the design variables correspond to the orientation of each layer, *T*<sup>0</sup> (Equations (17) and (18)) and, therefore, the number of design variables is four (one per ply for a half laminate). In the cases where steered fibers are considered, the design variables correspond to the angles *T*<sup>0</sup> and *T*<sup>1</sup> of each layer (Equations (17) and (18)). Consequently, there is a total of eight design variables (two per ply for a half laminate). Each design variable is allowed to continuously vary between −90◦ and 90◦, except when gaps were considered. In those cases the domain of design variables is allowed to vary between −89◦ and 89◦, due to singularities of the *Matlab* subroutine used. In some of the steered cases, the maximum curvature constraint, displayed in Equation (15), is imposed in order to guarantee that the laminate can be

manufactured using an AFP machine. Thus, in some steered cases the curvature of each ply is constrained by the maximum value of curvature allowed (*Kmax*). It is not necessary to impose the maximum curvature constraint in the non steered cases, because the fibers are rectilinear and therefore the curvature (*K*) is always zero. Thus, the optimization of the non-steered cases is an unconstrained problem. The boundary conditions considered are the fully clamped and fully simply supported boundary conditions.

$$\begin{aligned} \text{maximize} & \quad f\_1(\mathbf{x}) \\ \text{by varying} & \quad \mathbf{x} \in \mathbb{R}^n \\ \text{subject to} & \quad K\_n(\mathbf{x}) \le K\_{\text{max}}, \quad \text{if stored} \\ & n = 1, 2, \dots, \dots, N\_p \end{aligned} \tag{16}$$

$$\text{given} \quad \mathbf{x} = [x^{(1)}, \dots, x^{(i)}, \dots, x^{(N\_p)}]^T \tag{17}$$

such that <sup>−</sup> <sup>90</sup>◦ <sup>&</sup>lt; *<sup>x</sup>*(*i*) <sup>&</sup>lt; <sup>90</sup>◦ (17)

$$\begin{array}{ll}\text{where} & \mathbf{x}^{(i)} = \begin{cases} [T\_0^{(i)}]\_\prime & \text{if non-steered} \\ [T\_0^{(i)}, T\_1^{(i)}]\_\prime & \text{if steeered} \end{cases} \end{array} \tag{18}$$

The notation used to distinguish the different cases is orientation distribution boundary condition—unconstrained/constrained. The orientation distribution notation can be NS (non-steered), if only non-steered plies are considered or LS (linearly steered), if the fiber path orientation varies linearly across *x*∗ (Equation (1)). The boundary conditions fully clamped (CCCC) and fully simply suported (SSSS) are treated here as C and S, respectively. The unconstrained and constrained cases are denominated UN or CON, respectively. In the constrained cases, it is considered a maximum curvature constraint (*Kmax*) with two different values: 1.57 m−<sup>1</sup> [16], treated here as constraint A, and 3.28 m−<sup>1</sup> [32], denominated here as constraint B. For instance, LS-C-CON-A is a plate with a linear fiber angle variation, fully clamped and with a maximum curvature constraint of 1.57 m<sup>−</sup>1.

The plates considered are symmetric with a total of eight plies, each one with a thickness of 0.159 mm. The material properties considered are displayed in Table 2. The length in the *x*∗ axis, *a*, is kept constant and equal to 0.5 m throughout all the plates. A square plate and two rectangular plates with aspect ratios (*a*/*b*) equal to 0.5 and 2.0 are considered.


**Table 2.** Carbon epoxy Cytec® G40-800/5276-1 and resin properties [17].

### **4. Implementation**

The developed framework can be divided into two main components: the *Matlab* and the Finite Element Model (*Python*/ABAQUS) environments as it is represented in Figure 4. The optimization process starts in the *Matlab* where the optimizer, *Matlab*'s built-in genetic algorithm, is implemented and where each generation's population is created. After that, the population information is passed to a *Python* script where the geometry of the plate, property assignment, boundary and loading conditions are defined before running the structural analysis using *ABAQUS*. The input data generated depend if the case considered only non-steered, steered without imperfections and steered plates with imperfections. Following the structural analysis, its results are then again passed to the *Matlab* environment where they are post-processed and the resulting fitness is evaluated by the genetic algorithm optimizer. If there is convergence, the optimization stops, otherwise another population is created and the optimization process continues in the next generation.

**Figure 4.** Framework scheme.

### *4.1. Finite Element Model*

The geometry creation, element property assignment as well as the structural analysis is defined in *Python* scripts and performed by *ABAQUS*. The orientation of the fibers and material properties can differ in each iteration and it is generated as the input data by *Matlab*. In the non-steered cases all the elements of a ply have the same orientation, so the orientation is defined layer by layer and therefore *Matlab* only passes each layer orientation to the *Python Script*. In the steered cases, the element orientation in a ply is calculated in accordance with Equation (1), where the fiber orientation varies linearly along the *x*∗. As a consequence, the orientation can no longer be defined layer by layer and has to be defined element by element in each layer which leads to an increase of computational cost. *Matlab* generates each ply *T*<sup>0</sup> and *T*<sup>1</sup> for half laminate (due to its symmetry) which are read by the *Python* script. Then, each element orientation in a ply is assigned by substituting in Equation (1) each element centroid position in *x*∗ and its respective ply *T*<sup>0</sup> and *T*1. If defects are being considered in the steered case, each element property is dependent on its resin percentage. So, not only the orientation is defined element by element in each layer, but also the material properties. For that reason, the steered cases considering defects are those that require a higher computational cost.

It is chosen a mesh with a total of 2500, 5000, 1250, FSDT based shell S4R mesh elements [35] for the square plates and plates with *a*/*b* = 0.5 and *a*/*b* = 2.0, respectively. This choice was made after a mesh convergence study with the objective of simultaneously assuring a good compromise between the model convergence and the computational cost necessary to run a structural analysis.

### *4.2. Matlab Environment*

The main script starts by generating a population. The problem is constrained and the generated population should satisfy the optimization constraints which are defined in a different script. After this, each individual is checked by a script that saves the optimization history. That history has saved all the individuals already analysed and their respective fitness evaluation. So, in case a certain individual has already been evaluated, then there is no need to run the structural analysis and this step can be skipped. This saves computational time during the optimization process. If the individual has not yet been evaluated, the respective input data needed to run the Python Script and ABAQUS CAE is created.

The creation of input data depends on whether the plate considered is non-steered or steered, and if defects are or are not being considered. It is important to note that only symmetrical composites were considered in all cases, so the input data is only referred to half of the laminate. When defects are considered, namely gaps, they are modeled in a *Matlab* subroutine in the same way as in Fayabazh et al. [18]. Afterwards, each element resin percentage is calculated and passed to the Python Script where the equivalent element properties are calculated [18].

Following the input data creation the next steps that include the generation of the plate geometry, property and fiber path assignment, boundary and loading conditions definition and structural analysis are performed by the *Python* script which used *ABAQUS* to run the structural analysis. After the structural analysis is finished, the output data, in this case, the plate natural frequency, is passed to the *Matlab* script that defines the objective function in order to evaluate its fitness. Finally, after all the individuals of a generation are evaluated, the optimizer checks if the convergence has yet been achieved. If it meets the convergence criterion, the optimization process is terminated. Otherwise, the optimizer generates the next population using genetic operators such as selection, crossover and mutation and the optimization process continues. A population size of 80, a function tolerance of 10−5, a maximum number of stall generations equal to 25 and a constraint tolerance of 10−<sup>6</sup> are considered. This choice was made after a study where these parameters were varied with the objective of minimizing the computational cost necessary during the optimization process while assuring a very high probability of finding the optimal solution.

### *4.3. Validation*

The finite element model was verified by comparing its results with the linear natural frequencies obtained in the study performed by Akhavan and Ribeiro [10] and by Akbarzadeh et al. [20].

In the first study, imperfections are not considered and the fiberpath considered is also the one presented in Equations (1) and (2). The respective plate dimensions and material properties are presented in Table 3. The results, namely, the first three natural frequencies are shown in Table 4, where a fully clamped boundary condition was considered. In the same table it is possible to verify that the results obtained using the finite element model developed in this work are in accordance with those obtained by Akhavan and Ribeiro [10].

The plate properties considered in the second study are presented in Table 5, where square plates (*a* = *b* = 1 m) are considered. The comparison of results is presented in Table 6, where the dimensionless fundamental frequency (*ω*¯ = *<sup>ω</sup>a*<sup>2</sup> *h ρ E*2 ) calculated using the developed framework and obtained by Akbarzadeh et al. [20] are displayed and compared. It is worth mentioning that the results were compared considering a defect free plate and the same one with complete gap manufacturing defects. It is also important to note that Akbarzadeh et al. [20] considered a plate with a constant curvature fiber path, which can be the cause of the slight difference in the results displayed in Table 6.

**Table 3.** Characteristics of plate 1 [10].



**Table 4.** Natural frequencies [10].

**Table 5.** Prepreg and resin properties [20].


**Table 6.** Dimensionless fundamental frequency—*ω*¯ .


### **5. Results and Discussions**

This section is dedicated to present and discuss the optimization results which aim to maximize the fundamental frequency.

### *5.1. Non Steered Plates*

This section has the goal of presenting the optimal results of plates that only consider non steered plies. The first natural frequency is presented in Table 7 as well as the optimal orientation of each ply (*T*∗ *0i*). The computational cost is represented by the number of generations (Gen.) required to perform the optimization.

Focusing first on the square plates, it is observable in Table 7 that the fully clamped supported case achieves a higher first natural frequency (51.601 Hz) than the one obtained in the simply supported case (31.223 Hz). It is also noticeable that the optimal orientation of each ply is different in both cases. On the one hand, in the NS-C, the optimal orientation alternates between approximately −90◦ and nearly 0◦. On the other hand, in the NS-S, the exterior ply and its respective symmetric one have an orientation of approximately 45◦, while all the interior ones have an orientation near −45◦. That again proves that the optimal orientation of each ply is dependent on the imposed boundary conditions. Comparing now both cases, but focusing on the last column of Table 7, it is possible to verify that both required the same number of generations to perform the optimization.

The non steered rectangular optimal results are also presented in Table 7. When the aspect ratio is equal to 0.5, the optimal orientations are near zero degrees for all plies and for both boundary conditions. A similar behaviour is observed for an aspect ratio of 2, where all plies have an orientation near 90◦ for both boundary conditions. So, for non steered plies, these two cases are equivalent, because all plies are oriented parallel to the smaller side of the plate for both boundary conditions considered. Regarding the computational cost, it is visible that all the non steered cases have a similar number of function evaluations due to

the fact that all of them have the same number of design variables. The small discrepancy is explained by the randomness associated to the genetic algorithm used.


**Table 7.** Non Steered Results.

### *5.2. Defect Free Steered Optimization Results*

The optimal results for the ideal square steered plates are presented in Table 8 and shown in Figure 5. The optimization is performed considering unconstrained and constrained cases with different maximum curvature constraints (*Kmax* = 1.57 m−<sup>1</sup> (A) and *Kmax* = 3.28 m−<sup>1</sup> (B)) in order to assess the curvature effect on the final solution. All these cases are optimized for fully clamped and fully simply supported boundary conditions. Table 8 displays the optimal *T*<sup>0</sup> and *T*<sup>1</sup> for the first four plies (because the composite laminate is symmetric) and the respective first three natural frequencies for the cases considered. It is also visible the number of generations (Gen.) required to perform the optimization.

The optimization results for a square plate considering the fully clamped boundary condition are presented in columns four to six of Table 8. The unconstrained optimization (no maximum curvature considered), the constrained case A (*Kmax* = 1.57 m−1) and the constrained case B (*Kmax* = 3.28 m−1) are presented in the fourth, fifth and sixth columns, respectively. An observation that can be made is that the unconstrained case (LS-C-UN) is the one that achieves the highest first natural frequency (59.471 Hz) which represents an increase of 15.25% with respect to the square NS-C case (Table 7). However, six plies of this case have a maximum curvature (*Kmax* = 6.42 m−<sup>1</sup> for ply 1 and 8, *Kmax* = 6.52 m−<sup>1</sup> for ply 2 and 7, *Kmax* = 6.53 m−<sup>1</sup> for ply 3 and 6) higher than the *Kmax* values considered here. Consequently, this laminate can not be manufactured by an AFP machine. In spite of that, the result of the unconstrained case is useful not only to show where the optimal solution is, but also as, with the development of new manufacturing techniques that enable to produce plates with higher fiber curvatures, it could become a feasible design in the close future.

Focusing now in both constrained cases, LS-C-CON-A (*Kmax* = 1.57 m−1) and LS-C-CON-B (*Kmax* = 3.28 m−1), Table 8 shows that LS-C-CON-B is the case with the highest fist natural frequency (56.814 Hz) between the two. That constitutes an increase of 10.10% with respect to the square NS-C case (Table 7), while the LS-C-CON-A increase is only 4.58%. LS-C-CON-B is also the constrained case that allows a higher maximum curvature, which, again, reinforces the influence of this parameter in the frequency obtained. It is also worth noting that for both cases, *Kmax* is an active constraint for the first two layers and their respective symmetric ones. Therefore, it can be argued that lowering the maximum curvature allowed tends for the frequency to decrease. Regarding plies four and five of LS-C-CON-B, their maximum curvature approaches zero similarly to the unconstrained case. That is also the case for plies three and six in the LS-C-CON-A case.

The optimization results considering the fully simply supported boundary condition case are presented in columns seven to nine of Table 8. The unconstrained optimization (no maximum curvature considered), the constrained case A (*Kmax* = 1.57 m−1) and the constrained case B (*Kmax* = 3.28 m−1) are presented in the seventh, eighth and ninth columns, respectively. The unconstrained optimal solution (LS-S-UN) has a first natural frequency equal to 31.229 Hz which is very similar to the one achieved by both constrained cases. The same phenomena is verified with respect to the second and third natural frequencies. So, it is no surprise, that all the cases have similar *T*<sup>0</sup> and *T*<sup>1</sup> combinations and, consequently, also similar maximum fiber curvatures (*K*∗ *max*) in each ply. Regarding the *K*∗ *max* in each ply, it is important to note that its value in all cases and in all plies is near zero, which implies that the maximum curvature is not an active constraint. It also means that the unconstrained case is possible to manufacture using the AFP machine (because the maximum curvature of all plies is lower than both *Kmax* considered here) and that the first maximum natural frequency achieved is near the one obtained considering only rectilinear fibers. This last affirmation can be proved by comparing the first natural frequency values obtained in Table 8 with the one considering only non steered fibers in Table 7 for a fully simply supported boundary condition. The respective increase of the first natural frequency caused by the introduction of steered fibers is only 0.02%. So, it can be concluded that the use of non steered fibers is almost as efficient as steered fibers in order to maximize the first natural frequency of a square plate with a fully simply supported boundary condition.

**Figure 5.** Square plate optimal reference fiber paths layouts. (**a**)LS-C-CON-A defect free; (**b**) LS-S-CON-A defect free. (**c**) LS-C-CON-A with complete gap. (**d**) LS-S-CON-A with complete gap.

Focusing now on the plates with *a*/*b* = 0.5 (visible in Figure 6), columns two and three of Table 8, it is observable that the orientation of each ply is near 0◦ for the fully clamped boundary condition. Therefore, the optimal solution found is similar to the one obtained considering only rectilinear fibers (Table 7). So, it appears that the use of curvilinear fibers is not justifiable in this case. On the contrary, for the simply supported boundary condition, there is a slight improvement of 1.49% on the optimal solution with the introduction of curvilinear fibers. The influence of the curvature of the fibers is especially relevant because the maximum curvature is an active constraint in the first two plies and their symmetric ones (*Kmax* = 1.57 m<sup>−</sup>1). Regarding the plates with *a*/*b* = 2.0 (last two columns of Table 8), the maximum curvature is also an active constraint for the fully simply supported boundary condition. For this case, the orientation in the middle of all plies (*T*0) is near ±90◦, which is also the optimal orientation of all plies for the correspondent non steered case (Table 7). The orientation at the edge of the plies is clearly restricted by the maximum curvature constraint, that, as aforementioned, is an active constraint in all plies. Analyzing now the case with the fully clamped boundary condition, it is visible that all plies have almost zero curvature and an orientation near ±90◦. So, the optimal solution obtained is very similar to the one achieved in the non steered case.

**Figure 6.** Rectangular plate (*a*/*b* = 0.5) optimal reference fiber paths optimal layouts. (**a**) LS-C-CON-A defect free; (**b**) LS-S-CON-A defect free; (**c**) LS-C-CON-A with complete gap; (**d**) LS-S-CON-A with complete gap.

**Table 8.** Optimal defect free steered results.


### *5.3. Defect Free Steered Results with Complete Gap*

The previous section results do not consider the occurrence of any type of defect inside the plate, which is a limitation of the model used. Therefore, the results would be different from those shown. In order to better estimate that difference, the three first natural frequencies of the previous constrained optimal layouts are calculated considering the existence of gaps caused by the one sided cut with a zero percent coverage parameter of an AFP machine. The results are shown in Table 9 for square and rectangular plates considering the deposition of courses with 8, 16 and 32 tows. Each tow has a width of 3.175 mm. It is also visible the volumetric gap fraction (*υG*[%]) of each analysed case. As previously mentioned, in a one sided cut strategy the tows of each course are only cut and restarted on one side of the course, where the tow drop areas are formed. By fixing each tow width, if the number of tows increases, so does the total course width. In consequence, the number of course intersections inside the plate decreases, which leads to a reduction of the tow drop area. So, as expected, it is clear that the volumetric gap fraction increases in Table 9 with the decrease of the number of tows.

Another important aspect to notice in Table 9 is that all the first natural frequencies decreased when compared with those obtained considering defect free steered plates (Table 8). It is also observable that with the increase of the volumetric gap fraction, the fundamental frequency diminishes. This could imply that the existence of gaps is responsible for the decrease of the first natural frequency. The stiffness reduction in the areas where gaps are located can be the explanation for that frequencies' decrease. It is also worth mentioning that, in spite of that, the use of steered fibers in a square plate with a fully clamped boundary condition considering 32 or 16 tows can still achieve a higher first natural frequency when compared to a plate with non steered fibers (Table 7) for both *Kmax*. However, if the boundary condition used is a fully simply supported one, the results show that the opposite happens. The highest first natural frequency is instead achieved by the plate that uses non steered fibers. This may be explained by the fact that when considering ideally steered square plates, the obtained results were very similar to those obtained for non steered plates (only 0.02% increase). When further adding the effect of gaps, the first natural frequency decreases past the one obtained for non steered square plates. The same argument could be used regarding the fully clamped rectangular plate cases, where the achieved frequencies are lower than those obtained in the respective optimal non steered plates as expected. In the rectangular plates simply supported cases, it is observed that with the decrease of the number of tows, the gap area increases, and the fundamental frequency decreases until it is no longer advantageous to have curvilinear fibers considering this manufacturing process.


**Table 9.** Optimal defect free steered results with complete gaps introduced.

### *5.4. Complete Gap Steered Optimization Results*

It was shown in the latest section that the occurrence of gaps leads to a decrease of the fundamental frequency, which can be correlated to the reduction of the elastic properties in the areas where gaps are located. It was also shown that the higher the gap percentage is, or the lower the number of tows used, the lower the fundamental frequency is. However, the optimization has not yet been performed considering the Defect Layer Method included in the finite element model built, which is interesting in order to verify not only how the defects affect the fundamental frequency, but also to check how they affect the optimal laminated layout. With that in mind, the results presented in the following sections assume that the AFP machine deposits 32 tows, each one with a width of 3.175 mm for both fully clamped and fully simply supported boundary conditions. Only constrained cases are presented in order to assure that the plates could be manufactured.

The square plate results are displayed in Table 10 for both boundary conditions and maximum curvature constraints values considered. As expected, all the optimal fundamental frequencies achieved are lower than the respective optimal ones obtained with the absence of defects. For instance, regarding the square plate cases with a fully clamped boundary condition, fourth and fifth columns of Table 10, it is observable that there is a decrease of approximately 0.83% and 2.65% in cases LS-C-CON-A (Figure 5) and B, respectively, relative to the equivalent ones but without considering defects (Table 8). Nonetheless, it is still verified an increase of 3.71% and 7.19%, respectively, with respect to the NS-C solution (Table 7). Therefore, it can be argued that although a decrease of the frequency is verified by the occurrence of gaps, it is still advantageous to use curvilinear fibers to maximize the frequency for square plates with this boundary condition. It is also important to note that the optimal layout for both cases differs from the one obtained when no defects are considered. This evidence can be verified by comparing the optimal designs from Tables 8 and 10. Nevertheless, the maximum curvature constraint is still an active one in some of the plies and, when it is not, the maximum curvature obtained is near it. The reason why the optimal layout has changed when gaps are considered can be explained by comparing the results from Tables 9 and 10 when considering 32 tows. In Table 9 it is visible that LS-C-CON-A and B have a gap volumetric fraction (*υG*[%]) of 3.51% and 3.14%, respectively, while in Table 10 1.88% and 2.21% are the respective gap volume fractions. The reduction in the occurrence of gaps is one of the factors that explains that the optimal solution considering gaps is superior to the one found considering ideal plates with gaps inserted afterwards. So, when the layout is optimized with the occurence of gaps included it may be said that in order to maximize the fundamental frequency there is a trade-off between having the desired fiber curvature and the minimum gap volume fraction.

Regarding the square plate, but focusing now on the fully simply supported boundary condition (sixth and seventh columns of Table 10), it is again visible that the maximum curvature of all plies is near zero. In consequence, the maximum curvature is not an active constraint in both cases A and B. This phenomenon has been previously observed when analysing the same cases, but without defects being considered (Table 8). The plate layout is also similar to the one previously presented with the orientation of each layer tending to ±45◦, so the frequencies obtained are very similar to those obtained considering only non steered fibers (NS-S, Table 7) and steered fiber but without gaps included (Table 8). Since the optimal solutions found in this section still present some curvature, although tending to zero, some gaps occur. In consequence, the optimal fundamental frequency is 0.56% and 0.46% lower in cases LS-S-CON A and B, respectively, when compared with the respective non steered optimal solution (Table 7). So, because of the previously stated reasons, it can be argued that it is more advantageous to use only rectilinear fibers for square plates with this boundary condition.

The rectangular plate results are visible in Table 10 considering only *Kmax* = 1.57 m−<sup>1</sup> for both boundary conditions and aspect ratios. As it could be expected, the optimal fundamental frequencies achieved are lower than those obtained when defects are not considered. Focusing first on the plates with *a*/*b* = 0.5 (visible in Figure 6), there is a decrease of 0.31% and 0.74% for the LS-C-CON-A and LS-S-CON-A, respectively, with respect to their equivalent cases where complete gaps are not included (Table 8). On the one hand, in the fully clamped case, it is observable that the maximum curvature of each ply tends to zero. On the other hand, in the simply supported boundary condition, the constraint is active in some plies, as in the respective defect free cases. It is also worth noting that there is a decrease for the fully clamped boundary condition of approximately 0.31% when compared to the optimal solution considering only non steered plies (Table 7). As a result is can be said that it is preferable to use rectilinear fibers for this case. On the contrary, when the BC is simply supported there is an increase near 0.74% when compared


of gaps.

to the respective optimal NS case (Table 7), even taking into consideration the occurrence

**Table 10.** Optimal steered results with complete gaps.

A similar phenomenon occurs for plates with *a*/*b* = 2.0, however in these cases the frequencies achieved are much higher. The plate's maximum curvature tends to zero in all plies for the fully clamped boundary condition and it is near the maximum curvature constraint in some plies for the simply supported boundary condition (Table 10). It is verified a reduction in frequency of 0.73% and 2.21% for LS-C-CON-A and LS-S-CON-A, respectively, in relation to their optimal equivalent defect free cases (Table 8). Since there is a reduction of 0.73% in the LS-C-CON-A with respect to its respective optimal non steered case (Table 7), the non steered optimal solution is preferable for this boundary condition just like in the previous case. Focusing now on the simply supported boundary condition, there is a slight decrease of 0.14% with respect to its corresponding non steered optimal solution (Table 7). In addition, the maximum curvatures especially in plies 2 and 4, and their respective symmetric ones, are near the maximum curvature constraint. Therefore, it can also be said that it is also preferable to use non steered plies for this boundary condition as opposed to the previously case with an *a*/*b* = 0.5. Since the maximum curvature of the plies does not tend to zero, it is possible to speculate that a better solution could be found if the maximum curvature constraint imposed was higher.

In a similar way as with square plates, the optimal layout when complete gap imperfections are considered during the optimization differs from the one obtained for ideal plates. These differences in the orientations can be observed by comparing the optimal layouts in Tables 8 and 10. When adding the complete gap imperfections to the ideal optimal solutions (Table 9), it is visible for a plate with *a*/*b* = 0.5 that the *υG*[%] obtained is lower than the one obtained when defects are taken into consideration within the optimization framework (Table 10). On the contrary, for a plate with *a*/*b* = 2.0, *υG*[%] is lower in Table 10 than on Table 9 considering a course with 32 tows. Therefore, when the layout is optimized with the occurrence of gaps included it is once again visible that there is a trade-off between having the desired fiber path and the minimum gap volume fraction depending on the plate's geometry and boundary condition.

### **6. Conclusions**

An optimization framework has been implemented with the objective of tailoring the performance of laminated composite structures by maximizing the fundamental natural frequency. The optimization is performed for both the non-steered and steered fibers in order to evaluate the effects of curvilinear fibers on the structural performance. Furthermore, manufacturing issues such as tow drop gaps are also taken into account using the Defect Layer Method with the objective of assessing its effect on the fundamental frequency.

The steered plates were first optimized considering no manufacturing embedded defects. It was found that the fibers curvature effect on the fundamental frequency depends not only on the boundary condition considered, but also on the plate geometry. For square plates, it was observed that the use of steered fibers is advantageous when a fully clamped boundary condition is imposed, while the optimal fiber orientation is rectilinear when a fully simply supported boundary condition is imposed. However, the use of curvilinear fibers was found to be advantageous for rectangular plates with simply supported boundary conditions.

Before considering embedded defects within the optimization framework, complete gap defects were added to the optimal ideal steered plate designs. When tow drop defects such as gaps were considered, results show that these rich resin areas lead to a reduction in the fundamental frequency, which appears to be a consequence of a reduction of the material equivalent elastic properties. It was also verified that reducing the number of tows leads to an increase in the occurrence of rich resin areas and therefore to a higher decrease in the fundamental frequency.

The plates were then optimized considering embedded defects within the optimization framework. The results revealed that the optimal designs that consider ideal plates can differ from the ones where manufacturing defects are considered. It can also be verified that not considering defects within the optimization framework leads to a different (often worse) gap area percentage. In spite of the occurrence of complete gaps, it was observed that the use of steered fibers is still advantageous when a fully clamped boundary condition is imposed in square plates. Moreover, the optimal fiber orientation is also rectilinear when a fully simply supported boundary condition is imposed in square plates. Regarding rectangular plates, results show that curvilinear fibers are only advantageous for a rectangular plate with *a*/*b* = 0.5 and with a simply supported boundary condition.

It has been shown that the use of curvilinear fibers could be an effective way to tailor the natural frequency of laminated composite plates. Nonetheless, variable angle tow laminates are still a relatively recent field, therefore experimental validation and verification would significantly contribute to the improvement of the numerical models proposed in this paper. Future work would include critical buckling load and to account for the minimum cut length constraint of the AFP machines.

**Author Contributions:** J.C.: Methodology, Software, Validation, Formal Analysis, Writing—original draft. A.S. (Abdolrasoul Sohouli): Methodology, Software, Validation, Formal Analysis, Writing original draft. A.S. (Afzal Suleman): Supervision, Conceptualization, Resources, Writing—review & editing. All authors have read and agreed to the published version of the manuscript.

**Funding:** This research received no external funding.

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

**Informed Consent Statement:** Not applicable.

**Data Availability Statement:** Not applicable.

**Acknowledgments:** The authors acknowledge the Fundação para a CiênciaeaTecnologia (FCT), through IDMEC, under LAETA, project UIDB/50022/2020. A.S. (Afzal Suleman) also acknowledges the NSERC Canada Research Chair and Discovery funding programs.

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

### **Abbreviations**

The following abbreviations are used in this manuscript:


### **References**

