A Multiobjective Optimization Framework for Strength and Stress Concentration in Variable Axial Composite Shells: A metaheuristic approach

Abstract A metaheuristic approach for variable axial composites considering multiobjective optimization is investigated. The proposed methodology is based on the combination of three main parts: a methodology for defining the orientation of the fibers in the laminate, a structural analysis program (based on the Finite Element Method) and an optimization algorithm. It is important to highlight that a radial basis function (RBF), which describes a smooth fiber pattern, is generated using control points. The novelties of the present methodology consist of a proposal for a generalized parameterization technique, which allows the investigation of mechanical strength and stress concentration of variable axial composites. Thus, NSGA-II multiobjective genetic algorithm is used as optimization tool to define the fiber orientations. Besides, ax metaheuristic approach is used in situations when it is desirable to simultaneously minimize the stress concentration factor ( K t ) and a failure criterion index ( F I or Φ ). Two case studies are investigated: a double notched plate and a tube with a transverse hole. Graphical Abstract


INTRODUCTION
Composite materials have been applied in fields where a high strength-weight ratio is required.They replace conventional materials in applications when such substitutions are feasible and show advantages like the improved strength, stiffness, fatigue and resistance to impact, (Kaw 2006).
It is common to find composite materials as laminate plates or shells, which are produced by using a lamination process.During the lamination is possible to tailor the directional dependence of strength and stiffness of a laminate to match the design loadings, which actuate in the composite structures (Jones 1999) Usually, the laminas are reinforced by unidirectional fibers, which underestimate the potential of these materials in structural applications.However, due to the developments in robotics, new manufacturing processes such as TFP (Tailored Fiber Placement), AFP (Automated Fiber Placement), and CTS (Continuous Tow Shearing) have become a way to produce more optimized composite structures (Sobhani Aragh et al. 2021).These materials will be referred to as Variable Axial (VA) composites in the current paper, and their advantage is to steer the fibers to modify the stiffness and strength in situ.This allows, for example, adjusting the stiffness of a composite panel in order to increase its natural frequency (Guimarães et al. 2019).The CTS process is significantly outperformed by the TFP and AFP processes, where the similarity between the intended and produced components is better.The AFP process is a good choice for like-3D components, whereas the TFP has cheaper implementation cost and greater productivity.The main benefit of CTS is the layer deposition when compared to TFP.It is important to highlight that each kind of process can impose different restrictions in terms of fiber curvature that can be produced.However, in the present work, it was not considered any specific manufacturing process.
The efficiency of the resultant material is determined by the appropriate fiber trajectory, which must be established by the loads and requirements of the mechanical project.Depending on the field defined by the objective function, this optimal placement can be fairly complex.The current design approach relies upon a fiber direction parameterization technique in conjunction with an optimization algorithm.This methodology is found in many works in the literature since the ideal fiber trajectory is typically not obtained in a straightforward manner.
A common approach to parameterize the fiber placement is the use of a hypersurface to describe its angle at a given position.Honda et al. (2013) obtained the path of the fibers projecting over the finite elements mesh the tangential directions from the contour lines of a polynomial and use the mean curvature in conjunction with a failure criterion or the first modal frequency as objective functions in a multiobjective optimization framework.van Zanten et al. (2019) used a coarser manufacturing mesh as a hypersurface that has the angles of the fibers defined on its nodes.Through an inverse parametric transformation, these angles are interpolated to the centroids of elements, which describe the refined calculation mesh where displacement and stress components are evaluated.However, depending on the structure to be analyzed and on the complexity of the hypersurface, the optimization can present a high computational cost.Setoodeh et al. (2006) presents an alternative approach that uses lamination parameters to define the optimization problem with few design variables, which result in a convex space solution that allows to obtain a gradient optimization with fast convergence.Stodieck et al. (2017) also-applied lamination parameters as a base for their methodology using B-splines built from control points to describe the stiffness variation (via lamination parameters).Thus, a gradient-based optimization method (GCMMA) is used to identify the best value for buckling load.Finally, a genetic algorithm (GA) was used to extract curvilinear fiber paths from the target lamination parameters.Another technique to reduce the computational cost consists of replacing a burdensome structural analysis (as Finite Element Method -FEM) with a faster trained metamodel.Wang et al. (2020) developed a framework for axial buckling load optimization of VA composite cylinders using a reliability-based design optimization (RBDO) method, which takes the manufacturing parameters of the filament winding into account.The objective function of the optimization was the critical buckling load that is obtained through a Kriging metamodel trained with the results from finite element analysis.Almeida et al. (2019Almeida et al. ( , 2021) ) evaluated the performance of VA cylinders (manufactured by filament winding) for buckling load, obtaining significant improvements in the implemented objective function, and discussing the effects of imperfections, thickness build-up and non-symmetric laminates on the buckling behavior.
The above-presented frameworks usually have as objective function some failure criterion, natural frequencies, buckling load considering a manufacturing constraint, but not approaching stress concentration.Regarding unidirectional composites, Ueng and Lin (1987) proposed solutions for plates with two elliptical holes under normal loads by combining analytical solutions for a plate with a hole with the superposition method applied to the stress potential functions.One of the findings of this research is that, in contrast to an isotropic material, stress patterns remain the same, but stress magnitude changes.Weixing and Xinlu (1991) showed solutions for stress concentration on unidirectional laminated composite plates using FEM, and the results were compared to the analytical solutions obtained by Chiang (1998) and Structures, 2023, 20(6), e497 3/16 However, for stress concentration analysis in Variable Axial (VA) composites, only few works are found in the literature.
Vijayachandran and Waas (2022) studied the influence of steered fibers on the stress concentration around a hole in a laminate plate.The authors considered manufacturing parameters, and mechanical properties of the Robotic Automatic Fiber Placement (RAFP).A genetic algorithm was used as an optimizer, and reductions up to 33% in the stress concentration factor were obtained compared to a quasi-isotropic reference material.
Based on the aspects pointed out, in the present work, an investigation made by Santana et al. ( 2023) using a mono-objective optimization for mechanical strength is improved to a multiobjective optimization using the Non-Dominated Sorted Genetic Algorithm II (NSGA-II).It is important to highlight that, although this work presents some similarity to the previous one, the parameterization methods are essentially different.Control points (  ) are dispersed over the finite element mesh, and angle orientation values are assigned to them.These angles are then interpolated to the integration points using Radial Basis Functions (RBF).The Cartesian coordinates of integration points are obtained by the use of the finite element shape functions, which allow change from the natural coordinate system to the Cartesian coordinate system.The objective functions are the stress concentration factor (Kt) and a failure criterion (Failure Index -FI).One case study consists of a notched plate that has been studied analytically by Weixing and Xinlu (1991) and numerically by (Chiang 1998) considering conventional (non-VA) composites reinforced by unidirectional fibers.In another case study, a composite tube is analyzed, observing the capability of the VA composite to reduce the stress concentration while at the same time keeping safety values for failure criteria.Both Tsai-Wu and LARC03 failure criterion are used in the analyses.In summary, the main contribution of the present work is to show the effectiveness of the proposed generalized fiber parameterization methodology, which is a metaheuristic approach, highlighting the concurrence of mechanical strength and stress concentration in VA composite laminates performance.In addition, it is shown the importance to consider both parameters as objective functions for a multiobjective optimization algorithm for VA composite structures.

Proposed Approach: Finite Element Formulation, Fiber Orientation Parameterization and Failure Criteria
In the following subsections, the finite element formulation, the fiber parameterization framework and the failure criteria are presented.

Finite Element Formulation
An 8-node curved finite element is used in the structural analysis of the variable-stiffness laminated composite shell studied in this work.Ahmad et al. (1970) showed the formulation of this degenerated shell element, which considers first-order shear deformation theory, and Kumar and Palaninathan (1997) extended it to orthotropic materials.This analysis is performed with an in-house software developed in MATLAB (MATLAB 2010).The steps to implement the finite element model can be found in the previously cited references.In this subsection, only topics that are used to construct the fiber parameterization are covered.Figure 1 illustrates a finite element mesh with the nodes, the global coordinate system   , the curvilinear natural system   and the lamina coordinate system   , which is assembled at each integration point.The last system is used as a reference for positioning the fibers (the fiber angle at the integration point (  ) is a measured relative to the  1 component).and Structures, 2023, 20(6), e497  4/16 Since the coordinate systems are essential for the construction of the proposed framework, it is important to show some definitions.The   component is the unit normal vector of the element mid-surface located at a given integration point   whose coordinates are obtained using the shape functions of the finite element   (  ) as   =     (  ).This component is calculated as shown in Eq. 1: , , where, following the procedure developed by Hughes (2000),  , 1 and  , 2 are tangent vectors to the curvilinear natural coordinates  1 and  2 .These intermediate components are obtained as . (2) The remaining components of the coordinate system of the lamina are obtained by a procedure that ensures consistency between the finite elements.This approach certifies that the fiber path is consistent throughout the entire studied structure.The directions defined by the unitary vectors   and   are independent of the node numbering order of the elements and   is parallel to the plane composed by  1 and  3 and oriented towards the positive direction of the axis  1 (  .  ) > 0, except for a mid-surface perpendicular to the global axis  1 .
The unit vector   is oriented in the same direction as the global axis  2 (  =   ), if   is aligned to the global axis (  ‖  .  ‖ = ).For any other case,   is obtained as described in Eq. 3. (3)

Fiber Orientation Parameterization
The proposed framework consists of in assembling an interpolation function from control points, which are distributed over the area of interest of the structure that has angle values previously assigned angle values.The selected interpolation method is the Radial Basis Function (RBF), which is expressed by where ℎ�  � is a polynomial corrector that ensures the uniqueness of the fit and θ  is the angle interpolated at a node .
The variable γ  is a weight that is associated with the control point  for the RBF ϕ(. ), which accumulates the interactions with all surrounding control points as shown by Biancolini (2017).
It is assumed that ϕ(. ) is conditionally positive definite and ℎ is a linear polynomial.These assumptions allow to obtain the weights vector .The coefficients of the polynomial corrector  are obtained in terms of the angles assigned to the control points by the following system of equations where  is the interpolation matrix and  is the constraint matrix described in the literature (Biancolini, 2017).The angles at the integration points θ  are obtained by the interpolation of the nodal angles of a given finite element by using its shape functions as Composite Shells: A metaheuristic approach Pedro Bührer Santana et al.
Latin American Journal of Solids and Structures, 2023, 20(6), e497 5/16 1 ( , ) The selected RBF kernel is the thin-plate spline (Eq.7).In the present work, this kernel was selected since it has not any parameter to be adjusted and it has a variation-diminishing propriety that guarantees some smoothness of the approximations (Buhmann, 2003).Besides, the software developed by Chirokov (2023) was used to evaluate the radial basis functions obtained in the present work.

Failure Criteria
Both Tsai-Wu and LARC03 were used as failure criteria in the proposed approach.It is well known that Tsai-Wu failure criterion is very famous compared to LARC03.Therefore, the authors decided to show only the equations for LARC03, keeping in mind the equations of Tsai-Wu.Davila et al. (2005) considered two main failure modes in the composite material: matrix failures and fiber failures.For matrix failures, a set of criteria based on the concepts proposed by Hashin and Rotem (1973) and the fracture plane proposed by Puck (1998) are applied.
For the matrix under tension and shear stresses (Eq.8), the fracture planes (α) are normal to the plane of the plies and parallel to the fiber direction.(1 ) where    and    are the in-situ matrix strength in tension and in-situ longitudinal shear strength presented in Eq (10).In this work, the investigated case studies have thick plies (h > 0.7 mm) so  is obtained by and If the matrix is under compression and the fiber is under tension, the failure criterion (FI) for the matrix is given by 2 2 22 12 22 cos (sin cos ) where cos ( cos ) For this situation, the fracture planes are not necessarily normal to the ply (Equation 11).In this case, the Mohr-Coulomb (Salençon, 2001) effective stresses are used to calculate the angles of the fracture plane.In fact, in Eq. 11, α is the angle of the fracture plane, which is obtained interactively by the maximization of   and   , which is given by the transverse shear strength presented in Eq. 12.The terms   and   are the coefficients of transversal and longitudinal influence described in Equation 13.  and Structures, 2023, 20(6), e497 6/16 (13) The term   in Eq. 14 refers to the typical fracture angle for a laminated composite under transverse compression loading which is 53° while   and   are the transversal compression strength and   is the longitudinal shear strength.The superscript  refers to imperfections in the fiber alignment that are idealized as a local region of weaviness.The modified stresses for this situation are: where  is the sum of a small imperfection angles and rotational components due to shear loading.It is defined as and   is the total misalignment angle for the case of pure axial compression loading, described as For fiber in tension, the failure criterion is based on the effective strain acting along the fibers.The noninteracting maximum allowable strain for the LARC03 criterion is defined by: where  11 is the local strain in direction 1 and  1  is the strain strength in this same direction.Finally, for fiber and matrix in compression, the failure criterion yields while for fiber in compression and matrix into tension, the failure criterion is given by: (1 ) In summary, Figure 2 shows the main aspects of the proposed optimization procedure and how these topics are linked together.The design variables of the NSGA-II algorithm are angles assigned to the control points (  ).These angles are used as a framework to build the radial basis function, which is used to interpolate angles at the finite element nodes.Finally, the shape functions are used to interpolate the nodal angles to the integration points as shown in Figure 3.The stress values and the failure criterion indices were extrapolated to the nodes using the approach proposed by Hinton and Campbell (1974).The NSGA II software implemented by Song Lin ( 2023) was used, and a population of 200 individuals and 400 generations were assumed.In Figure 6, Pareto fronts show that the objective functions are concurrent and are suitable to have a multiobjective optimization process.
The obtained numerical results and the values of the fibers angles (or directions of the axis) at the control points are organized in Table 2.In addition, the fibers patterns, maximum principal stress and LARC03 failure criterion indices (FI) for the Scotchply 1002 and T300/N5208 laminate are presented in Figure 6.For a unidirectional laminate with the fibers aligned with the prescribed displacement, the Scotchply 1002 laminate presents a  of 0.1183 and a   of 4.4619, while the T300/N5208 laminate presented a  of 0.4239 and a   of 5.8469.A maximum reduction of 45.6% and 64.4% in LARC03 failure criterion indices is obtained for Design I using Scotchply 1002 and T300/N5208 laminates, respectively.The maximum reduction in   for both laminates is given in Design III, being 66.3% for Scotchply 1002 and 94.2% for T300/N5208.
The same laminate structure is optimized using unidirectional fibers for the same objective function, allowing the comparisons between the obtained results of the resultant VA designs and unidirectional ones.Three additional cases are proposed here: (a) a one-layer plate, (b) a symmetric laminate plate with stack sequence 1 2 2 1 and (c) a laminate plate with stack sequence 1 2 3 4 / / / θ θ θ θ     .The generated Pareto fronts obtained from the optimization processes executed for these plates are presented in Figure 7, considering that the lateral constraints for the angles are ∓  2 .
For the one layered plate, the minimum FI was obtained for [90] (all fibers aligned to the y axis, i.e. aligned to the load) and the minimum   was obtained for [0].For the symmetric plate, the minimum FI was obtained with a stack sequence [90/90] s and the minimum   with a stacking sequence [3.94/79.12]s .The last proposed stack sequence presents [90/90/90/90] for minimum FI and [0/49.78/90/88.60]for minimum   .
The three proposed unidirectional laminates resulted in the same maximum mechanical strength, which was already compared to the obtained VA Design.Based on Figure 7, the single-layered laminate (a) and the symmetric one (b) show similar behavior up to a   of 2.67.From this point, the symmetric laminate presents more flexibility, allowing it to reach a   of 1.46.Although this Design presents less stress concentration than Design III, it shows a higher penalization in the mechanical strength.The 1 2 3 4 / / / θ θ θ θ     shows higher flexibility than the one-layered plate and the symmetric laminate (case b).For reduction in the stress concentration factor, from a   of 1.67, it is more worthwhile to use the unidirectional fibers since it maintains lower FI and is easier to manufacture.In Figure 8, the improvement in the strength for Design I is explained by the fact that the fibers reinforce the structure in  direction near the circular cutout, which is the critical region for failure.At some distance from the circular cutout, the fibers converge to the plate edge, transferring a portion of the load to the matrix, as well as reducing the stress and the   .Design II and Design III show that as smaller the   , as more acute is the angle of the fibers since they reach the notch.This behavior increases the actuation of the matrix, decreasing the stresses in the notch region while penalizing the mechanical strength.

Tube with a transverse hole
This study case is based on the case study presented in Bi et al. (2020).It consists of a tube with a through hole.The tube has a length of 0.8 m with diameter of 0.03 m and the holes have a diameter of 0.01 m.The Scotchply 1002 laminate structure is made of five layers with 1 mm each.The geometry, mesh details and prescribed displacement are shown in Figure 9.In Figure 10, it is possible to observe twelve control points (in blue), which are distributed around the area of interest (VA region).However, as the structure has symmetry, the control points on the lower part of the structure have the same values as the upper one, then for the optimization process, only nine control points remain.The proposed approach was used for the optimization, considering a population of 200 individuals and 200 generations during the analysis.
In addition to the LARC03 criterion, in this case, the application of the Tsai-Wu failure criterion described in Eq. 22 was also investigated.
( ) The resulted Pareto fronts for both criteria are presented in Figure 11, and the coefficients generated for the selected designs, as well as their numerical results are presented in Table 3.It is possible to note that both fronts have very similar shapes.For a unidirectional composite with the fibers aligned to the direction of the displacement, it was obtained a Tsai-Wu failure index of 0.1406 and   of 4.8654.Design I is not able to offer a significant improvement (only 7.5%) over the baseline unidirectional laminate showing that, for mechanical strength value, the unidirectional configuration is the best for this laminate Increasing the number of layers with different fiber paths or studying other control points distributions could expand the solution space and lead to better results.In Figure 12, for Design I, it is noticeable that the fibers are almost straight near the hole, reinforcing the hypothesis that the straight unidirectional baseline design is near the optimal configuration for mechanical strength.In Design II and Design III, the same behavior noted in the previous case was observed with the fibers assuming a certain angle near the hole reducing the portion of the load supported by the fibers and, consequently, decreasing the stress and the mechanical strength in the cutout region.In Design I, the Tsai-Wu failure index increases by 104.8%, and the   decreases by 34.8%, while in Design II, there is an increment of 344% in the Tsai-Wu failure index, and the   decreases by 61.1%.The fibers patterns in this case follow a similar behavior that the one described in the plate, thus it is necessary to have a balance between mechanical strength value (FI) and stress concentration factor (  ).For LARC03 failure criterion, the unidirectional baseline design has a failure criterion value of 0.0788, and Design I show an improvement of 14% in the mechanical strength value and 12.6% in   .This more significant improvement in the failure criterion value is attributed to the accuracy of LARC03 in detecting the failure of the laminate since there are slight differences between fiber paths obtained using each of the failure criteria, as shown in Figure 13.

CONCLUSION
The present paper proposed a framework for multiobjective optimization of VA composites.The stress concentration factor (  ) and mechanical strength values, or failure criterion (FI), were then explored as objective functions.It was demonstrated that they are suitable (concurrent) for the multiobjective optimization process, resulting in a Pareto's front.
A double notched plate and a tube with holes were investigated.Referring to the plate, when compared to a unidirectional composite laminate where the fibers are aligned to the prescribed displacement, both objective functions are improved for all points of the Pareto's front for the T300/N5208 laminate.Regarding the tube, there is no significant improvement in terms of strength to the unidirectional design, but it was able to reduce the stress concentration close to the holes.

Figure 1 :
Figure 1: Coordinate systems of the finite element model.
Finally, T eff τ and L eff τ are effective stresses in such plane.Composite Shells: A metaheuristic approach Pedro Bührer Santana et al.Latin American Journal of Solids

Figure 2 :
Figure 2: Flowchart of the proposed approach.

Figure 3 :
Figure 3: Representation of fiber directions interpolated to the integration points.

Figure 8 :
Figure 8: Obtained fiber paths, max stress and failure criterion on the notched region.

Figure 9 :
Figure 9: (a) Geometry and FE mesh details.(b) Details of the region of the circular cutout.(c)Lateral view.

Figure 10 :
Figure 10: (a) Side view of the cylinder with control points (b) top view of the tube with control points and its coordinates.

Figure 12 :
Figure 12: Obtained fiber paths, max stress and failure criterion on the circular cutout region.

Figure 13 :
Figure 13: Design I and Design III fiber path for (a,b) LARC03 and (c,d) Tsai-Wu.
Latin American Journal of Solids

Table 2 .
Design angles (radians) for the obtained fiber paths and results obtained for the double notched plate.

Table 3 .
Design angles (radians) for the obtained fiber paths and results obtained for the tube with a transverse hole.