Variable Sti ﬀ ness Composites: Optimal Design Studies

: This research work has two main objectives, being the ﬁrst related to the characterization of variable sti ﬀ ness composite plates’ behavior by carrying out a comprehensive set of analyses. The second objective aims at obtaining the optimal ﬁber paths, hence the characteristic angles associated to its deﬁnition, that yield maximum fundamental frequencies, maximum critical buckling loads, or minimum transverse deﬂections, both for a single ply and for a three-ply variable sti ﬀ ness composite. To these purposes one considered the use of the ﬁrst order shear deformation theory in connection to an adaptive single objective method. From the optimization studies performed it was possible to conclude that signiﬁcant behavior improvements may be achieved by using variable sti ﬀ ness composites. Hence, for simply supported three-ply laminates which were the cases where a major impact can be observed, it was possible to obtain a maximum transverse deﬂection decrease of 11.26%, a fundamental frequency increase of 5.61%, and a buckling load increase of 51.13% and 58.01% for the uniaxial and biaxial load respectively.


Introduction
Fiber-reinforced polymers have been used over the years in many different application areas [1], wherein the use of straight fibers with selected orientations within each ply deserved a significant research attention. More recently, with the technological advances in fiber placement [2] the manufacture of variable stiffness (VS) composites emerged, bringing the possibility to consider non rectilinear fiber paths within each ply, hence achieving coordinate dependent' fiber angles. These advances conferred enhanced tailoring ability to these composite materials, hence a greater capability to meet structures' requirements, which has been studied in recent years by a number of researchers. Among the published research work until this moment, it is, thus, relevant to give a brief overview on the diversity of the studies performed on these advanced composites subject.
Hence, one may refer the work due to Setoodeh et al. [3] which investigated the use of cellular automata (CA) on the optimization of the fiber angles to improve the in-plane stiffness since CA is appropriate for parallel implementations. Different situations were analyzed revealing a good convergence and robustness rate and significant gains in using VS. Following this study, Abdalla et al. [4] introduced the generalized reciprocal approximation to maximize the fundamental frequency. The design variables were linked to the nodes, instead of the elements, to assure the continuity of the distribution of the lamination parameters, being the maximization made at each node separately. The results showed 2 of 19 that higher frequencies could be achieved by using VS composites. Later, Setoodeh et al. [5] extended this method to buckling analysis, obtaining the same conclusions for the buckling load. In the context of buckling finite element analysis, Lopes et al. [6] studied the improvement of buckling load and first-ply failure of VS laminates taking into account tow-drops and overlaps. It was demonstrated that VS revealed higher buckling load and first-ply failure, but the results could be hindered if the resin-rich regions were taken into account. Gürdal et al. [7] used the classical lamination theory, the Rayleigh-Ritz method and the Trefftz criterion to develop an iterative analytical approach to analyze the buckling response on VS laminates. They concluded that it is possible to change either the buckling load or the in-plane stiffness while keeping the other constant. Hao et al. [8] studied also the buckling behavior of variable-stiffness panels but using a different approach based on isogeometric analysis. This study was performed for a linear fiber angle and a flow field fiber angle variation for different ply quantities, geometric parameters, boundary, and loading conditions being the results compared with the ones obtained by finite element analysis. The results obtained showed that the isogeometric analysis is well suited for the analysis of VS panels and that the geometric stiffness has a significant influence on the buckling load and on the corresponding buckling mode shapes.
Honda et al. [9] developed an analytical method, based on the Ritz method, to study the vibration of VS composite rectangular plates. They applied the method to a single-ply with parabolic shaped fibers and compared the results with straight fiber plates, verifying that the plate with the parabolic fibers can have higher frequencies and that they have a strong effect on the vibration modes. Houmat [10] used a three-dimensional elasticity theory and the finite element method to study the vibrations of VS rectangular plates and analyzed the inter-laminar modal flexural and the modal transverse shear stresses and the modal cross-sectional warping. The author concluded that the inter-laminar modal flexural stresses are discontinuous because of the difference of the fiber angle between adjacent layers, the sign of modal transverse shear stresses can change through the thickness due to the change of the fiber angle over in-plane dimensions and that the shape of modal cross-sectional warping is influenced by the mode number and stacking sequence of layers.
A p-version finite element based on a third-order shear deformation theory (TSDT) and considering geometric non-linearity was used by Akhavan et al. [11] to determine the deflections, the normal and the shear stresses of VS composite laminate under transverse loads. They concluded that due to the changes in the maximum stresses, the VS laminates are not always better than constant stiffness composite laminates.
Venkatachari et al. [12] studied the effects of the radius-to-thickness ratio, the thickness ratio, the temperature, the moisture, and the boundary conditions on the fundamental frequency for doubly curved, cylindrical, and spherical panels. They used a first-order and a higher order shear deformation theories and concluded that the pattern of frequency variation is similar for both cylindrical and spherical panels and that the fiber angles have a significant effect on the fundamental frequency. Sarvestani et al. [13] developed a semi-analytical methodology to conduct hygro-thermo-mechanical analysis on thin to relatively-thick fiber-steered composite, conical and cylindrical panels as well as circular plates. The authors registered up to 57% and 44% improvements in the buckling loads and fundamental frequencies respectively.
In the field of the optimization of these VS composites performance it is also possible to find some works although with a restricted scope. Among these works, Demir et al. [14] optimized the lamination parameters based on a least-squares finite element, considering a continuity constraint to find a manufacturable layup, having obtained good compliance values. Following this work, Shafighfard et al. [15] optimized the fiber angle direction for open-hole geometries by minimizing the compliance taking as design variables the lamination parameters, and afterward obtaining manufacturable designs. These authors achieved smaller compliances and maximum stress concentrations for the VS and determined that the optimal fiber angles did not match with the principal stress directions.
Hao et al. [16] developed a multi-stage design method based on lamination parameters. The method begins by firstly using B-spline surfaces to fit the lamination parameters followed by an isogeometric analysis along with a gradient-based method to attain the optimal stiffness. In short, the lamination parameters are converted into fiber angles using an evolutionary algorithm, and the fiber paths are obtained by using flow field functions, ensuring their manufacturability, and, lastly, the buckling load is taken as the objective function. The gradient-based method is used again for final optimization. This hybrid method showed good global optimization and high computational efficiency.
To improve the in-plane flexibility and the bending stiffness of wing skins, Murugan et al. [17] studied which material would be better suited, within a set of materials. The effects of boundary conditions and aspect ratio on the out-of-plane deflection were considered in the minimization of the in-plane stiffness and maximization of the bending stiffness, by spatially varying the volume fraction of the fibers. They concluded that the aspect ratio and the boundary conditions have a considerable effect on the out-of-plane deflection and the variable fiber volume can contribute to enhance the desired characteristics. Similarly, Murugan and Friswell [18] studied the effect of curved fiber paths on the in-plane flexibility and out-of-plane bending stiffness representative of a morphing wing skin to maximize the in-plane deformation and, simultaneously, minimize the out-of-plane deformation. They concluded that the curved fiber paths have a considerable influence in the in-plane flexibility and out-of-plane bending stiffness, varying with the aspect ratio. Furthermore, the ply's away from the neutral axis tend to be curved while the ones near tend to be straight.
Zhangming Wu et al. [19] proposed a formulation to represent non-linear varying fiber angles in which the coefficients of polynomials directly reflect the fiber angles at selected reference points. With that formula describing the fiber angles evolution, an optimization was performed to obtain the maximum buckling load, revealing similar results to ones in the literature using lamination parameters.
With this work one presents a study that aims at providing a transversal characterization of the optimal mechanical behavior of variable stiffness composites, from the linear statics perspective to the free vibrations and linear buckling perspectives, for a single layer and a three-layer configuration composite plates.
Hence, in the first part of this study, one presents the main aspects associated to the methodology used, followed by a set of verification case studies.
The application case studies are divided in two parts: The first one, wherein a set of static and buckling analyses aim at complementing a free vibrations analysis considered in the context of the verification case studies; and the second part and the main part of the present work, where a complete and transversal set of optimization studies are presented. In this last context, one performs a set of constrained optimization processes focused on a single layer and on three layers VS composite plates, where the objective functions are the minimization of the maximum transverse deflection, the maximization of the fundamental frequency, and the maximization of the critical buckling load. The design variables are the characteristic angles associated to the fibers' path in each ply.

Static Buckling, Free Vibrations, and Static Analyses
As mentioned, the present work aims at considering a set of optimization case studies that will be based on the free vibrations, static buckling, and also on the static analyses of variable stiffness composite plates. To these purposes the Lagrangean functional for plate elements is written as follows [20,21]: where the first and second terms correspond to the first and second order elastic strain energy, and represents the stress components associated to the initial state of stress previously determined via a linear static analysis. The third term denotes the kinetic energy and the last one the work performed by the external forces. By applying Hamilton principle and considering the whole discretized domain, one obtains the following equilibrium equation: which if one considers a linear static instability analysis yields the equation: In this equilibrium equation, K σ is the geometric stiffness matrix, K is the stiffness matrix, and q i is the eigenvector associated to the eigenvalue λ i . To the minor eigenvalue will correspond the critical buckling load.
On the other hand, if one intends to perform a dynamic free vibration analysis, and by assuming free harmonic vibrations, the equilibrium equation achieved will be: The eigenvalues of this equation correspond to the square of the natural frequencies associated to the vibration mode shape given by the eigenvector q i . To the minor eigenvalue square root will correspond the fundamental frequency of the plate.
If instead, one aimed at performing a static analysis, wherein the second order elastic strain energy and the kinetic energy will be null, then the equilibrium equation for each plate element will be given as: All these analyses were implemented through the finite element method, considering the use of the first order shear deformation theory (FSDT). Next to obtaining the element matrices and vectors and after carrying out their assembly to reproduce the whole discretized domain, the corresponding boundary conditions were imposed, so as to achieve the solutions for each specific combination of design variables.

Variable Stiffness Composites
The fiber within each ply follows a curvilinear path, defined according to the linearly varying angle along the x-coordinate, proposed by Gürdal and Olmedo [22]: where T 0 is the angle at the center of the plate and T 1 is the angle at the edges of the plate. The variable (a) denotes the length of the plate and θ x (x) is the angle at the x coordinate, being the notation < T0|T1 > used to represent such paths. For an easier visualization, Figure 1 depicts < 20|60 > and < 50| − 20 > fiber paths configurations, where x, y directions define the plate and each layer planes.
To avoid, as much as possible, kinks in the fibers the maximum curvature is submitted to the constraint [23]: where f (x) is the fibers' path function and ρ the curvature radius. According to this constraint, the manufacturing domain is represented in Figure 2, where, the region in white denotes the fiber angles combination domain that yield feasible manufacturing layers, while the grey region represents the fiber angles combination domain that will exceed the maximum allowed curvature.
where T is the angle at the center of the plate and T is the angle at the edges of the plate. The variable ( ) denotes the length of the plate and ( ) is the angle at the x coordinate, being the notation < 0| 1 > used to represent such paths. For an easier visualization, Figure 1   To avoid, as much as possible, kinks in the fibers the maximum curvature is submitted to the constraint [23]: where ( ) is the fibers' path function and the curvature radius. According to this constraint, the manufacturing domain is represented in Figure 2, where, the region in white denotes the fiber angles combination domain that yield feasible manufacturing layers, while the grey region represents the fiber angles combination domain that will exceed the maximum allowed curvature. To perform the finite element analyses, each plate finite element was assigned a fiber angle according to the element center x coordinate value, so the more elements along the x coordinate, the less the angle discretization error, when compared to the analytical definition, as shown in Figure 3.  To perform the finite element analyses, each plate finite element was assigned a fiber angle according to the element center x coordinate value, so the more elements along the x coordinate, the less the angle discretization error, when compared to the analytical definition, as shown in Figure 3. To avoid, as much as possible, kinks in the fibers the maximum curvature is submitted to the constraint [23]: where ( ) is the fibers' path function and the curvature radius. According to this constraint, the manufacturing domain is represented in Figure 2, where, the region in white denotes the fiber angles combination domain that yield feasible manufacturing layers, while the grey region represents the fiber angles combination domain that will exceed the maximum allowed curvature. To perform the finite element analyses, each plate finite element was assigned a fiber angle according to the element center x coordinate value, so the more elements along the x coordinate, the less the angle discretization error, when compared to the analytical definition, as shown in Figure 3.  In this figure one presents for a square plate, the angles associated to each element of 10 × 10 and 50 × 50 discretized domains. The solid line represents the exact angle and the dashed lines are the angles in each element.

First Order Shear Deformation Theory and Constitutive Relation
The displacement field associated to the first-order shear deformation theory [24,25], describes the plate kinematics as follows: where u(x, y, z), v(x, y, z), and w(x, y) represent the displacements of a plate arbitrary coordinate point, u 0 (x, y), v 0 (x, y), and w 0 (x, y) stand for the displacements of a mid-plane' point associated to it, and θ 0 x and θ 0 y denote the normal to the midplane' rotations, respectively. The generalized strains are obtained by considering the elasticity kinematical relations for small deformations, and by considering the assumptions associated to this theory and the material characteristics, one can write for each variable stiffness composite layer, the constitutive relation: The general expressions for the transformed elastic stiffness coefficients Q ij are given in literature, for constant stiffness layers. However, in the present case, the stiffness coefficients are not constant within each layer. They depend on the x coordinate, thus on θ 0 and θ 1 angles as denoted in the constitutive relation. Additionally, as according to this theory the equilibrium conditions are not satisfied at the plate's upper and lower surfaces, a shear correction factor k is needed [20,26,27].

Optimization
A typical minimization problem is commonly formalized as: where f (x) represents the objective function being x the design variables vector and g m (x) and h n (x) denote respectively behavioral inequality and equality constraints. The design variables are also subject to side constraints x L , x u . [28][29][30].
The optimization method used in the present work is an adaptive single objective method, which combines an optimal space-filling design of experiments, a kriging response surface, and mixed-integer sequential quadratic programming. This mathematical optimization method, based on a response surface, which enables in a lighter way the search for the global optimum, uses a minimum number of design points strictly necessary to allow building the kriging response surface [30,31].
In the present work, the optimization case studies are constrained optimization problems, focused on single objective functions, with continuous design variables. Concerning to the constraints, besides the side constraints one has also considered the constraint associated to the maximum curvature in a layer by layer basis. To perform the optimization studies, the finite element codes developed in Matlab [32] were linked to a hybrid adaptive optimization procedure [30].
This adaptive optimization technique can be described according to the present pseudocode, in a summarized form: This Algorithm 1 combines an optimal space-filling (OSF) design space, a kriging response surface and mixed-integer sequential quadratic programming (MISQP). Based on this hybrid approach the algorithm searches for the optimal solution based on the response surface obtained, which requires a minimum number of design points to build the kriging response surface [31]. This approach enables considering a reduced number of design points for the optimization process, which provides lower computation time requirements.
The process starts with an initial population with a number of individuals that will remain constant along the process. The design space is subdivided into a specified number of divisions which will also remain constant. After each domain reduction, done after each evaluation phase by assessing the most promising design space region, this new domain with closer delimitating borders will be again subdivided. In order to maintain the same number of design points, as some of the previous will be eliminated, the required number will be added.
The mixed-integer sequential quadratic programming technique will then search within the response surface generated for each search domain, considering different starting points. The identified optimum candidates are then evaluated, considering the kriging error predictor criterion. The stability of each solution will be then assessed by considering an additional refinement of the kriging's response surface, which after its validation a domain reduction centered on the accepted candidates will be performed. If the candidate is rejected, a new verification point is calculated and inserted in the current kriging response surface as a refinement point, restarting the MISQP process. The global optimum is achieved when the achieved solution is stable and is not susceptible to improve according to the parameters considered.

Case 1: Natural Frequencies of Isotropic Plates
An isotropic plate with a side ratio of α = a/b is considered along with four different boundary conditions, here generically denoted as XXYY, with the first and second X representing the sides in x = −a/2 and x = a/2, respectively, and the first and second Y the sides in y = −b/2 and y = b/2. The letter C will signify a clamped edge, S a simply supported edge, and F a free edge.
The results are shown in a non-dimensional frequency way, by considering In Table 1 one presents the results obtained using the classical theory and the first-order shear deformation theory. These results were compared to reference alternative solutions [26] which considered a classical approach. The relative deviations (dev) among the present classical theory results and [26] solutions are denoted within parenthesis and calculated as dev = res−res re f res re f × 100%. As it is possible to conclude from Table 1, there is a very good agreement between the present classical shear deformation theory model and the ones obtained by the reference. It is also possible to observe that the present first-order model also demonstrates a good performance. Table 2 presents a similar study but now for a side ratio α = 2.5. Although the reference does not consider the CCCC boundary conditions for this aspect ratio situation, in this study it was considered the inclusion of these results for completion purpose. These specific results are just compared among the present two models (CLPT and FSDT). Again, it is possible to understand the good performance of the models for the rectangular plates' non-dimensional frequencies. To note, in a more specific manner the good performance of the FSDT model.

Case 2: Buckling Critical Loads of Constant Stiffness Composite Plates
This second verification case considers a (0 • /90 • )s composite plate with an edge to thickness ratio of 10 (a/h = 10). The material properties relations are defined as G 12 /E 2 = 0.6 and ν 12 = 0.25. The aim of this case is the evaluation of the present models' buckling response for different degrees of orthotropy (E 1 /E 2 ), and the comparison of their results to an alternative solution [33]. The buckling response of the plates when submitted to a uniaxial compression load N x applied along the x direction, is presented in a non-dimensional form as: The results obtained are presented in Table 3, for the present FSDT model. In this case one has considered also for illustrative purposes, a situation where the shear correction factor assumes the value 5/6, and another situation where, as previously mentioned, this factor is calculated [27].
The relative deviations among the present results and the higher-order shear deformation theory alternative solutions [33] are denoted within parenthesis and are calculated as dev = res−res re f res re f × 100%. From a global perspective, the results obtained with the present first-order shear deformation model present a reasonable agreement with the higher-order model (HSDT) of the reference authors [33]. It is also possible to conclude that the calculation of the shear correction factor considering the specific characteristics of the laminate is able to provide results that are closer to the ones given by the higher-order theory, which would be expected.

Case 3: Natural Frequencies of Three-Layer Variable Stiffness Composite Plates
This last verification case considers a simply supported three-layer VS square plate, and its natural frequencies. The geometrical and material properties are described in Table 4. The results obtained with the present FSDT model are presented in Table 5.  Table 5 allows concluding on the very good agreement between the present FSDT results and the ones obtained by [34] that used a third-order shear deformation theory. The classical theory results obtained by [26] over predicts the natural frequencies as it would be expected.

Numerical Applications
The material properties used in this section cases studies are presented in Table 4. Each layer thickness is now set to 10/3mm. The stacking sequences of the variable stiffness laminates studied in the last verification study are now used to analyze the static deflection of the plate when subjected to a uniformly distributed pressure under different sets of boundary conditions. The pressure magnitude values used both in the static analysis and in the optimization study, are presented in Table 6.

Static and Buckling Analyses of Three-Layer Variable Stiffness Composite Plates
The maximum transverse displacements presented by the different plates, obtained with the FSDT model can be observed in Table 7. The boundary conditions in these cases correspond to the simply supported ones, already used in the verification case of Section 3.3. As we can observe, it is not possible to say that a specific stacking sequence is able to provide the best performance for all the boundary conditions considered. This is an expected behavior due to the variable stiffness character of these composites, being particularly relevant to note that this constitutes an additional design parameter to consider in the context of composite materials' tailor ability.
These stacking schemes were also considered in the context of linear buckling analysis, where two situations were studied; one related to a uniaxial compression state, along the x coordinate, and another consisting in a biaxial compressive load, with a ratio F x /F y = 1 and a magnitude of 100 kN. Table 8 demonstrates the buckling load multipliers (λ) for the first six modes. From the results obtained it is possible to conclude that as it could be expected, biaxial compression states yield lower buckling load multipliers regardless the instability mode one considers. It is also possible to observe that the first VS configuration provides the higher buckling loads, hence it would be a better solution considering the buckling problem.

Optimization
In this sub-section, it is presented a first set of optimization studies where a single VS layer is considered, in order to minimize its maximum transverse displacement. It is also considered the maximization of its fundamental frequency and its critical buckling load when submitted to a uniaxial and biaxial compression state. Next, one considers a three-layer configuration and one performs a similar set of optimization studies.

Single Layer Variable Stiffness Composite Plates
A single layer simply supported plate is considered in order to achieve the best configuration concerning the angles T 0 and T 1 (Figure 1) that will minimize the maximum transverse deflection (w max ) of the plate when subjected to a uniformly distributed pressure, with the magnitudes shown in Table 6, under a different set of boundary conditions, being the results shown in Tables 9 and 10. Concerning the remaining boundary conditions, SSSS, SFSF, and CFCF, although demonstrating better results for the VS plate, the improvement can be considered residual for the SFSF case while being more significant in SSSS boundary conditions' case. Hence in these last analyses' conditions, it can be concluded that in some cases there would be no particular advantages on the use of VS composites, contrarily to what was observed in buckling analyses and also in free vibration analyses.
One has considered now the maximization of the fundamental frequency (ω max ) and the critical buckling load multiplier (λ max ) when subjected to a uniaxial compression load along the x direction, and a biaxial compressive load, with F x /F y = 1, with a magnitude of 1 kN. In the optimization process, one has considered that the angles could vary in the interval [−90 • , 90 • ] considering additionally the curvature constraint already referred in Section 2.1.
Tables 11 and 12 present the results obtained for the fundamental frequency maximization and, for comparison purposes, the results for a constant stiffness (CS) unidirectional composite plate. As we may conclude, the optimized VS plate shows a better result for most cases when compared to the CS plate. In this case, it was possible to achieve maximum improvement of 8.164% for the SSSS case. Tables 13 and 14 present the optimal configurations obtained for the maximum critical buckling loads associated with the uniaxial compression state. As we may conclude, the optimized VS plate shows a better result for some cases when compared to the CS plate. In this case, it was possible to achieve a maximum improvement of 15.958% for the SSSS case.  Tables 15 and 16 presents the optimal configurations obtained for the maximum critical buckling loads associated with the biaxial compression state. As we may conclude, the optimized VS plate shows a better result for most cases when compared to the CS plate. It is particularly relevant to note that in the SSSS case, it was possible to achieve a maximum improvement of 43.448%.  The variation of the critical buckling load multiplier and the fundamental frequency with the change of T 0 and T 1 was further studied for the SSSS case. The surfaces were obtained by setting the angles from −90 • to 90 • with a step of 22.5 • . Figure 4 depicts the corresponding contour plots. The variation of the critical buckling load multiplier and the fundamental frequency with the change of T0 and T1 was further studied for the SSSS case. The surfaces were obtained by setting the angles from −90° to 90° with a step of 22.5°. Figure 4 depicts the corresponding contour plots. All these plots demonstrate to be antisymmetric in relation to the line T0 = T1. In all cases, although with less significance for the fundamental frequency, when the angle on the center of the plate (T0) is higher the results are lower.

Three-Layer Variable Stiffness Composite Plates
In the present subsection, the previously analyzed three-layer configuration was considered in the context of the optimization of its static, free vibration, and critical buckling load behavior. Hence, All these plots demonstrate to be antisymmetric in relation to the line T 0 = T 1 . In all cases, although with less significance for the fundamental frequency, when the angle on the center of the plate (T 0 ) is higher the results are lower.

Three-Layer Variable Stiffness Composite Plates
In the present subsection, the previously analyzed three-layer configuration was considered in the context of the optimization of its static, free vibration, and critical buckling load behavior. Hence, one has proceeded to the minimization of the maximum transverse deflection, and to the maximization of the fundamental frequency and critical buckling load.
Tables 17 and 18 present the results obtained for the minimization of the maximum transverse deflection of VS plates submitted to different boundary conditions. From these tables, it is possible to draw a similar conclusion to the one of a single-layer VS plate.
In the SSFF, CCCC, and CCFF case the CS composite with unidirectional fiber can be clearly used as they provide the same results as VS composites do. In the remaining boundary conditions' cases, the achieved optimal stacking sequence provides an effective contribution to lower the maximum transverse displacement of the plates, with more significance to the SSSS and CFCF cases.  Figure 5 depicts the color gradient contour plot of the optimal solution for the constant stiffness composite (a) and the variable stiffness composite (b) in the case of the simply supported boundary conditions.
As it is possible to observe, in the second case, the shape profile of the deflection surface is more balanced and the color gradient illustrates also the difference between the quantitative values presented in Table 17.
The maximization of the fundamental frequency of the simply supported three-layer composite was addressed next. The results obtained both for a constant stiffness and for a variable stiffness composite are presented in Table 19 and the corresponding mode shapes depicted in Figure 6. The maximization of the fundamental frequency of the simply supported three-layer composite was addressed next. The results obtained both for a constant stiffness and for a variable stiffness composite are presented in Table 19 and the corresponding mode shapes depicted in Figure 6.    As in the single layer composite previously optimized, it is possible to conclude on the improvement of the maximum frequency of the VS composite plate, although to a less extent. It is also visible from Figure 6 that the vibration mode shape for the optimal configuration VS case acquires a more balanced profile toward biaxial symmetry, quite different from the one presented by the CS composite. It can also be observed that VS plates present in most cases, the outer layers with similar angles. It can also be observed that for the cases analyzed the angles of the outer layers present amplitudes close to the ones corresponding to the single-layer optimal configuration. This is, with terminology adaptations, applicable to the CS plates, which is an expected result considering the higher contribution that these layers provide to the fundamental frequency that corresponds to a flexural mode.
The final optimization studies are related to the determination of the optimal layups that will provide the maximum critical buckling load multipliers both in uniaxial and biaxial compression load cases. The results obtained can be observed in Table 20 and the corresponding color gradient contour plots are visible in Figures 7 and 8. Table 20. Optimal angles of a three-layer VS plate for static buckling load multipliers. As in the single layer composite previously optimized, it is possible to conclude on the improvement of the maximum frequency of the VS composite plate, although to a less extent. It is also visible from Figure 6 that the vibration mode shape for the optimal configuration VS case acquires a more balanced profile toward biaxial symmetry, quite different from the one presented by the CS composite. It can also be observed that VS plates present in most cases, the outer layers with similar angles. It can also be observed that for the cases analyzed the angles of the outer layers present amplitudes close to the ones corresponding to the single-layer optimal configuration. This is, with terminology adaptations, applicable to the CS plates, which is an expected result considering the higher contribution that these layers provide to the fundamental frequency that corresponds to a flexural mode.
The final optimization studies are related to the determination of the optimal layups that will provide the maximum critical buckling load multipliers both in uniaxial and biaxial compression load cases. The results obtained can be observed in Table 20 and the corresponding color gradient contour plots are visible in Figures 7 and 8.  From Table 20 it can be concluded again that VS composites can have significant importance in the maximization of the critical buckling load corresponding to a uniaxial compression load case. This was already seen in the context of a single-layer. The aspect related to the magnitude of the angles presented by the outer layers is also here confirmed.  Table 20 it can be concluded again that VS composites can have significant importance in the maximization of the critical buckling load corresponding to a uniaxial compression load case. This was already seen in the context of a single-layer. The aspect related to the magnitude of the angles presented by the outer layers is also here confirmed. Concerning the corresponding critical mode shape, it can be observed that there is an elongation trend more visible in contrast to the one presented by the constant stiffness composite in the case of the uniaxial compression load. A similar elongation pattern is observed in Figure 8, for the biaxial critical buckling load, although more balanced when compared to the constant stiffness composite plate. Concerning the corresponding critical mode shape, it can be observed that there is an elongation trend more visible in contrast to the one presented by the constant stiffness composite in the case of the uniaxial compression load. A similar elongation pattern is observed in Figure 8, for the biaxial critical buckling load, although more balanced when compared to the constant stiffness composite plate.

Conclusions
Due to their known fiber curvilinear paths, VS composites can provide greater design flexibility when compared to constant stiffness laminates, as they may enable for an eventual adjustment to geometrical specificities, such as holes, and contribute for a better load redistribution. Moreover, the possibility of selecting the better configurations for specific boundary conditions can be a very important feature to improve structures' mechanical performance while considering functional or manufacturability constraints.
This work presents a study on the optimization of VS composite plates considering the minimization of the maximum transverse deflection, the maximization of the fundamental frequency, and of the critical buckling load both in a uniaxial and biaxial compression load condition.
A set of studies based on first-order shear deformation theory were developed in a preliminary phase of the work for verification purposes and also with the aim of providing a more complete set of analyses namely for the three-layer VS composite, prior to the optimization case studies.
The optimization studies, based on an adaptive hybrid technique, have considered single-layer and three-layer composite plates, having as design variables the characteristic fiber angles, T 0 and T 1 , within each layer. To guarantee manufacturability conditions, according to results obtained by other authors, the present optimization studies took into account that feasibility constraint.
From the results of the single layer case, it can be seen that the VS composite plates present a particularly improved behavior when they are simply supported on all their edges, having the maximum transverse deflection decreased by 14.30%, the fundamental frequency has increased by 8.16%, and the buckling load multipliers has increased by 15.96% and 43.45% for the uniaxial and biaxial load, respectively.
For the three-layer configuration, the VS composite plates also shown significant improvement over the CS composite ones, particularly for the simply supported boundary condition, having the maximum transverse deflection decreased by 11.26%, the frequency increased by 5.61%, and the buckling load multipliers increased by 51.13% and 58.01% for the uniaxial and biaxial load, respectively.
In a global appreciation, VS composite plates present better results when compared to CS composite ones. However, the results obtained allow concluding that in some situations that depend not only on the boundary conditions but also on the nature of the objective function to be assessed, the advantage of the former may be not so significant.
The achieved results provide a transversal overview on the potential advantages of variable stiffness composites, nevertheless for each specific optimization case, eventually ruled by other constraints, other specific studies may be required.