Optimizing Variable-Axial Fiber-Reinforced Composite Laminates: The Direct Fiber Path Optimization Concept

The concept of aligning reinforcing fibers in arbitrary directions offers a new perception of exploiting the anisotropic characteristic of the carbon fiber-reinforced polymer (CFRP) composites. Complementary to the design concept of multiaxial composites, a laminate reinforced with curvilinear fibers is called variable-axial (also known as variable stiffness and variable angle tow). The Tailored Fiber Placement (TFP) technology is well capable of manufacturing textile preforming with a variable-axial fiber design by using adapted embroidery machines. This work introduces a novel concept for simulation and optimization of curvilinear fiber-reinforced composites, where the novelty relies on the local optimization of both fiber angle and intrinsic thickness build-up concomitantly. This framework is called Direct Fiber Path Optimization (DFPO). Besides the description of DFPO, its capabilities are exemplified by optimizing a CFRP open-hole tensile specimen. Key results show a clear improvement compared to the current often used approach of applying principal stress trajectories for a variable-axial reinforcement pattern.


Introduction
Recently, the demand for energy efficient systems leveraged the use of CFRP lightweight composites in structural components. These materials are increasingly being employed in aeronautical, aerospace, and automotive applications. Due to the high cost of carbon fibers, their efficient usage becomes essential [1]. By employing a variable-axial (VA) fiber design, stiffness and strength properties may be improved when comparing to classical CFRP designs [2]. Thereby, the term VA means varying the fiber orientation at the ply level. The desired performance of CFRP composites is achieved by guiding the loads almost exclusively along the fiber orientation and thus minimizing the shear load of the matrix. For a technical realization, TFP technology, which was developed at Leibniz-Institut für Polymerforschung Dresden (Germany), is well suited. Basics and some applications of TFP technology are described in [3,4]. The placement of carbon fibers is usually carried out by stitching dry rovings, as shown in Figure 1. The roving is guided through a rotatable roving pipe onto a base material, where a sewing thread applied in the zig-zag-pattern holds it in place.
Several approaches have been developed to optimize VA composites. An extensive overview of curvilinear fiberreinforced composites was recently performed by Ribeiro et al. [5]. Under the name variable angle tow steering, Weaver et al. [6] improved the postbuckling performance of composite panels with a VA layout, whereas Panesar and Weaver [7] optimized blended bistable laminates suitable for morphing flap applications. Duvaut et al. [8] implemented a varying fiber density in order to consider local stress intensity. For a similar purpose, the local layer thickness was varied by Parnas et al. [9] as an additional design parameter. Groh and Weaver [10] proposed a minimum-mass design of a typical aircraft wing panel under end-compression. Khani et al. [11] developed a mathematical optimization algorithm for variable stiffness panels using lamination parameters. Van Campen et al. [12] proposed a methodology to convert known lamination parameters distribution for a VA composite laminate into realistic fiber angles, with minimum loss of structural performance. Cho and Rowlands [13] reduced stress concentrations in an open-hole laminate with a genetic algorithm.
In contrast to optimization procedures, principal stress criterion has been often used for deriving curvilinear fiber paths, e.g., [14][15][16][17]. Kelly et al. [18], Waldmann et al. [19], and Malakhov and Polilov [20] designed the curvilinear fiber path based on the concept of aligning fibers following the load path by placing fibers along principal stresses. Both approaches were assumed to be optimization criteria, although no mathematical optimization process was explicitly carried out and only an a priori design criterion was applied. Furthermore, gradient based numerical optimization processes [9,[21][22][23] and optimization approaches based on evolutionary algorithms [24][25][26][27][28] to design VA composites have been also developed.
However, due to the often chosen approach of varying angles of single quasi-isotropic finite elements (FE), the number of design variables is considerably high, also increasing computational costs. Therefore, most examples have been computed by using FE models with a limited number of finite elements. Increased design freedom causes increased design problem complexity. For example, Conti et al. [17] found that using fiber angles as design variables inevitably leads to an ill-behaved objective function with many local minima.
Usually, a VA fiber pattern can be implied in a varying density of fibers, which causes a nonuniform thickness of the dry preforms. Hence, this heterogeneous thickness build-up is extremely complex to be accounted for in the analysis. Given this complexity, current approaches neglect the thickness accumulation and consider only the fiber angle variation.
Thus, the major criticism of many state-of-the-art optimization approaches is that there is no mathematical optimization procedure to begin or to operate without necessary information on the manufacturing process due to the lack of an appropriate modeling procedure. However, the knowledge of the thickness distribution and the local fiber orientation corresponding to an arbitrary fiber layout, which can be produced by TFP, is essential for the part design process. Spickenheuer et al. [1,29] and Albers et al. [30] made initial attempts to separate the optimization process of a curvilinear fiber-reinforced composite manufactured via TFP from the actual numerical models in order to limit the number of required design variables, making them independent of the applied FE mesh resolution. Thus, once a sufficiently accurate modeling of VA fiber layouts is established, optimization techniques can be applied to the fiber pattern. This allows fully utilizing the high degrees of freedom in the design process and the maximization of the anisotropic material characteristics of CFRPs.
Given the identified gaps in the current state-of-the-art in properly modeling VA composites, this work presents a novel design procedure for illustrating the capability for generating a VA pattern for an open-hole tensile specimen, where an optimal fiber pattern cannot be easily derived. Hence, the novel optimization approach, called Direct Fiber Path Optimization (DFPO) for VA composites, will be introduced and numerically evaluated on the example of an open-hole tensile specimen.

Model Setup.
According to the state-of-the-art, modeling of composite structures is mostly limited to stacking layers with a constant thickness and a constant fiber angle within each layer. Models for structural analysis of uniform spirals and single curved tapes of parallel fibers and constant thickness have been successfully applied additionally by an increase in modeling efforts [31]. In this case, an analytic description for local preform thickness and the fiber orientation is known, which can be used to build appropriate FE models for structural simulation. However, the existing model limitations are too strong if one plans to apply optimization strategies to fully exploit the potential of CFRP manufactured by TFP and to adapt production requirements. Although there are many approaches of employing the mathematical description of the optimization algorithm to deduce the numerical model, e.g., the geometry of each iteration step, this work is going to describe the modeling independent of the optimization and thus as a generic module for any VA structure with a similar placement characteristic, strictly following the manufacturing characteristics of TFP. The objective of the modeling is the elastic description of the laminate compliance and the prediction of initial failures based on a physically based failure criterion. A mesoscaled model is used to evaluate the specific properties of the TFP process following the recommendations raised by Uhlig et al. [32].
To generate continuous layers with TFP, the rovings have to be placed with a slight overlap to avoid gaps between them. If placed with constant thickness and fiber orientation, parallel laminates can be produced for a small range of distances between neighboring rovings. The thickness is calculated according to the following (see Figure 2): where A rov is the roving cross-section area, the distance between neighboring rovings, the roving fineness, the fiber volume fraction, and the fiber density.
For arbitrary nonparallel roving placement, the thickness evaluation becomes more complex. As a starting point for  (1). Figure 3: Placement pattern as a sequence of straight lines.
the analytical description of such a preform the placement path is used, which is the basis for the fiber placement with a TFP machine. This path or more generally a sequence of paths will be referenced as the design pattern. The simplest mathematical description is a sequence of straight lines in two dimensions. Curved placement paths, e.g., containing primitives such as arcs or splines, will be approximated with a sequence of short straight lines within the accuracy of production. A straight line is defined by the starting point ( 0 , 0 ) and end point ( 1 , 1 ) inside the placement plane. All points in between are described by where is the parameterization variable ranging from zero to one. Note that either the total sequence of straight lines can be connected, in the case if there is just one fiber path, or at least some succeeding lines are not connected, which represents completely separate fiber paths, as can be seen in Figure 3. This type of design pattern contains only the information of the fiber placement paths including the length of rovings but no information about the width or cross-section area. A formal extension of the straight path information combined with the cross-section area V is the line thickness distribution line for line segments with length : Here, the concept of the Dirac delta distribution was used to define the density function. This line thickness distribution function is only an intermediate step as the total fiber volume is concentrated along the infinitely thin lines and an infinitely high thickness is obtained on the lines and zero elsewhere. However, the function already fulfills the normalization condition. By integrating over the total design space or any area containing all line segments the total fiber volume is obtained: For practical purposes this thickness distribution is not very useful as it lacks the information about the width of a typical roving which is placed by an embroidery machine. This width usually depends on the type of rovings used, e.g., the number of filaments, material density, and most importantly, on a machine parameter, the width of the zigzag stitch used to fix the roving on the base material.
By convolution with different smoothing functions, the information about the width of the roving can be added. A very convenient approach is the coarse-graining by convolution with a Gaussian centered at ( , ) of width determined by the placement width to obtain the Gaussian thickness distribution: The coarse-graining is done by integrating the function in the whole plane of , . By using the definition of the line thickness distribution , a solution for this convolution can be expressed in terms of error functions. This makes a numerical implementation very fast. This Gaussian thickness distribution represents a Gaussian weighted average of the roving volume density in the area around the point at which the thickness needs to be computed. For single straight rovings, this thickness distribution leads to a Gaussian crosssection area, which roughly approximates the real crosssection areas for the TFP process, as shown in microsections by Uhlig et al. [32]. Other smoothing functions such as a cylindrical average approximate the cross-section of a single roving to a closer degree. However, the resulting laminates exhibit many discontinuities, which negatively influence the convergence of the modeling. With the Gaussian thickness distribution the laminate boundary needs to be defined by a cut-off thickness, as the Gaussian is nowhere exactly zero.
The main challenge for numerical modeling is to obtain the geometry and the fiber orientation based on the placement pattern of a single layer. Successive layers are stacked on top of each other without regard to draping behavior, which is fine as long as thickness gradients of the lower layers are small enough. The description is restricted to layers of noncrossing rovings or at least to roving placements where overlapping rovings cross at small angles, such that an element wise average of fiber orientation is meaningful. Note that, for many examples which contain a self-crossing roving path, the layer can be split into smaller layers with noncrossing rovings. In Figure 4, a schematic description of the modeling procedure is shown. Based on a twodimensional (2D) mesh of the planar design space a threedimensional finite element model is derived using localized information of the Gaussian thickness distribution and the averaged fiber orientation. The thickness is evaluated at each corner node and the fiber orientation at the center of each element. The fiber orientation is well defined for linear line segments. The elemental fiber orientation is averaged by a thickness weighted average of all linear line segments, which contribute to the total thickness at the center point of each element.
Successive layers can be stacked on top of each other. The resulting three-dimensional (3D) FE model represents a piece-wise linearization per element of the locally averaged characteristics, namely, thickness and fiber angle. Alternatively, the thickness and fiber angle can be combined at the center of the FE into a 2D layered shell element description to obtain a model for the same fiber layout with less computational cost. The main difference arises from neglecting the out-of-plane component of the fiber orientation and thickness gradients within an element.
Next, two numerical examples are considered by using DFPO. For both cases, the following parameters are used: fiber volume fraction ( ) of 58%; roving fineness ( ) of 400 tex; density ( ) of 1.76 g.cm −3 ; and width smoothing parameter ( ) of 1 mm.

Case 1: Open-Hole Tensile
Specimen. An open-hole specimen under tensile loading is chosen to demonstrate the modeling capabilities and the optimization of VA laminates by employing the DFPO approach. Their specimen geometry and dimensions are presented in Figure 5(a). In order to directly evaluate the capabilities of the proposed optimization framework, the specimen comprises two layers, achieved by stacking a carbon fiber TFP layer (layer to be optimized) on top of the base material (±45 ∘ woven fabric with area weight of 256 / 2 ). Figure 6 shows in detail the two-layer openhole specimen in study.
Based on a 2D meshing of the supporting plane, the local thickness is evaluated at each node for each laminate layer along the FE mesh, as it is shown in blue color scale (Figure 7(b)). In addition, the elemental fiber orientation (Figure 7(a)) is set as the averaged fiber orientation at the center of each element. In areas where the current considered fiber pattern places no rovings, the thickness computation yields effectively zero. However, to provide a continuous mesh in this case, a very small thickness of 0.001 mm is set at the corresponding nodes and the corresponding element material properties are set to resin properties (blue elements in Figure 6). The corresponding FE model additionally incorporates at the bottom of the laminate a layer of constant thickness (0.24 mm) of base material, as Figure 6 depicts.
Symmetrical boundary conditions are applied along all axes of the specimen. The load is applied at the top-edge of the specimen. These details can be seen in Figure 6. Finite element simulations are carried out in ANSYS APDL using quadratic SOLID186 and linear SOLID185 elements (ANSYS library reference).

Case 2: Narrow-Middle Tensile Specimen.
In order to provide another example for the applicability of the proposed DFPO framework, a sample under the same loading conditions has been considered. For that, a narrow-middle specimen under tensile loading is analyzed and optimized. Details on the geometry and dimensions of the narrowmiddle tensile specimen are shown in Figure 5(b). In order to evaluate the capabilities of DFPO, similarly to the openhole specimen, the sample consists of two layers, attained by stacking a carbon fiber TFP layer (layer to be optimized) on top of the base material (± 45 ∘ carbon fiber woven fabric with areal weight of 256 . −2 ). The material properties of both UD carbon fiber/epoxy TFP layer and the carbon fiber/epoxy woven fabric laminated composites used in the FE models and optimizations are presented in Table 1.

Optimization Process
The optimization problem for the fiber path is described by the minimization of an objective function, in which compliance minimization is the objective function, which analogously stands for stiffness maximization ( ) under variation of each roving placement path C i . Within the context of the actual optimization, two compliances are aimed to be minimized as follows: where minimizing the maximum of the displacement indirection is the objective function for stiffness optimization, whereas minimization of the maximum of MIA (mode interaction parameter) is the objective function for strength optimization. This MIA parameter is related to the physically based failure mode concept developed by Cuntze [33].   Concept (FMC) is based on the stress and strengths quantities, which means that MIA (failure parameter) is calculated based on the stress state of the laminate at each interaction along the analysis. In other words, if ≥ 1, then the laminate fails; analogously if < 1, the laminate is safe. Additionally, all failure modes can be combined into a single numerical value suitable for optimization with the mode interaction (MIA) quantity. Since the whole formulation of Cuntze's FMC is very extensive, its full description can be seen in [33].
Mathematically, the dimensionality of the optimization problem of even a single roving path is infinite. However, due to limited production accuracy, the placement path can be modeled using a finite set of parameters within some placement path representation.
The optimization flowchart is implemented and presented in Figure 8. The parameterized fiber layout is represented by a finite set of coefficients, e.g., spline control points. The 2D fiber path is computed which in turn is analyzed by the 3D modeling tool to generate the finite element model. The local thickness and fiber orientation are taken into account. Loads and boundary conditions are applied and then the model is solved. Based on this solution, the target optimization value (compliance minimization or stiffness  maximization) is derived. The optimization value is the sole input value for gradient free optimization algorithms, such as BOBYQA (Bound Optimization BY Quadratic Approximation) by Powell [34], which can modify the fiber path parameters within predefined boundaries to achieve a minimal displacement value. As long as no gradients are derived only gradient free optimization algorithms can be used. BOBYQA provides a fast converging algorithm for smooth optimization functions due to its quadratic approximation also implementing box constraints that can be used to restrict the fiber pattern to within reasonable locations. Details on the optimization parameters are given in Section 3.2. In general, other optimization values, such as failure stress, can be applied. However, the convergence to overall good solutions is much better for stiffness optimization in comparison to strength optimization. Thus, for a strength optimization, a stiffness optimized layout is used as an initial layout.

Convergence Study.
For the use in optimization procedures, the numerical model must be sufficiently stable and free of mesh dependence, once otherwise numerical fluctuations lead to nonconverging behavior in the optimization algorithm.
For layers that fully cover the design space such that neighboring rovings overlap, i.e., maximum displacement in -direction max( ) (Figure 9(a)) and maximum MIA (Figure 9(b)), the simulation converges or stabilizes with increasing the number of elements ( ), as Figure depicts. Regarding stiffness optimization (Figure 9(a)), the FE model composed of quadratic elements (SOLID 186) easily converges for any element size, whereas for the FE model with linear elements (SOLID 185), the model converges well with a minimum number of 20,000 elements. On the other hand, for strength optimization (Figure 9(b)), both element types take a while to converge, but for a mesh density of 200,000 elements, the FE model converges for both linear and quadratic elements. In this way, for the stiffness objective function, the mesh with 20,000 elements has been employed in all further optimization and FE analyses.
The convergence is only achieved if the boundaries of the rovings overlap the previous and next rovings, thus forming continuous layers without gaps. If the rovings do not fill the whole mesh, this base mesh elements need to be aligned along the bounding contour of the fiber layers to allow a realistic material description per element.

Open-Hole and Narrow-Middle Specimens Optimization.
In this section, the parameterization of the fiber layout is described in more detail. For both examples, only the 0 ∘ layer is optimized. However, in general, multiple layers can be parameterized in a similar way and the collective parameter sets are combined to form a single optimization parameter vector. A basis or an initial fiber layout is chosen and the parameterization describes only modifications of this layout. For the 0 ∘ layer of both examples, a layout of equidistant straight and parallel fibers is chosen as an initial layout. Deviations from this layout are restricted to shifts indirection (see coordinate system in Figure 6), which limits possible layouts to angles of less than 90 ∘ between fiber orientation and the load which is parallel to the -direction. In addition, closed loops cannot be described with such an approach. The angle limitation is useful especially if multiple layers are considered where fiber layers are assigned to specific "tasks", which should not be exchanged between layers during the optimization. (Closed loops and abruptly ending fibers within the part are also impractical for production with TFP.) Similar to Nagendra et al. [35] the fiber path is modeled based on 2D cubic B-splines. However, only deviations from the initial path are described, the straight and parallel fiber layout in this case, with the spline functions. The x-coordinate of the placement path is given by where are spline basis functions and the control points define the optimization parameters. The linear scaling factors and determine the total length scale. An equidistant set of defines the different rovings next to each other in xdirection and the total set of curves for each roving path along the y-direction is given by variation of V: at the beginning and increase up to 112 obtained by node insertion after BOBYQA algorithm converges for a lower resolution. In principle, BOBYQA algorithm converges even for larger number of optimization parameters of several hundreds of parameters. However, the manufacturing precision limits meaningful increase of the resolution. The optimization is considered to be converged if the control points do not change by more than 0.005 mm between successive iterations. The initial resolution of 16 parameters converges in about 60 iterations and takes about 10 min in a typical workstation. Figure 10 shows the various layouts of the open-hole specimen. The reference layout with equidistant and parallel fibers is given in Figure 10(a)), the stiffness optimization result is in Figure 10(b)), and for comparison the result of a principal stress orientation of fibers is given in Figure 10(c)). The optimization results provide a different solution when compared to previously optimized fiber pattern for openhole tensile specimen, as can be seen in [9,14,28], where they employed the principal stress criterion. Not surprisingly, DFPO achieves better improvement than those ones. The disturbance of fibers reaches much farther away from the hole, such that globally straighter fibers with overall similar length are obtained.

Results and Discussion
In addition to the open-hole specimen, another example is provided to demonstrate the potential of the DFPO framework for another case. Then, a tensile specimen with a narrow section in the middle is considered, where the ratio of the narrow section to the full width is 50%. Due to the smooth transition region of the narrowed section, the principal stress layout (Figure 11(c)) works very well in this case and more fiber rovings divert from the straight path (Figure 11(a)). The DFPO solution is qualitatively similar to the open-hole solution but with stronger fiber concentration due to the stronger narrowing of the defect (Figure 11(b)). In addition, in this case, the effect of the optimization using DFPO is much more "global" compared to the principal stress layout. Figure 12 presents the stiffness and strength increase of the optimized fiber layouts relative to the reference design ( Figure 10(a)). The principal stress oriented layout (Figure 10(c)) yields to a 5% increase in stiffness (Figure 11(a)) (20% for the second example) and about 139% increase in strength in terms of Cuntze fiber failure mode interaction max(MIA) (Figure 12(b)) (237% for the second example), whereas the DFPO-optimized layout (Figure 10(b)) results in about 9% increase of stiffness (25% for the second example) and 197% increase in strength (275% for the second example). Please note that the boundary conditions of the optimizations were such that the total number of rovings next to each other was fixed and thus the volume and mass change for different fiber layouts. However, the increase in volume of 0.6% for principal stress and 1.6% for DFPO (5.9% and 6.3%, respectively, for the second example) is smaller than the gain in both stiffness and strength.
Mathematical Problems in Engineering In contrast to the principal stress design, DFPO represents a real optimization procedure and consequently takes global and not just local features of the specimen into account. The thickness distribution is nonuniform in both cases and a thickness concentration near the defect of the structure is observed. In the DFPO case this thickness concentration extends further from the defect area than in the principal stress layout. The fiber length of single rovings is much more uniform along each family of specimen for the DFPOoptimized, such that the load balance of all rovings under tensile load is better. Compared to other optimization techniques where elemental fiber orientations and thickness values are optimized without correlations induced by endless fibers, in DFPO, each fiber layout considered in every optimization iteration is already manufacturable and no subsequent adaptation is necessary. Thus, these gains obtained by the optimization can be fully transferred to the application.

Conclusions
The key objective of this investigation was to present a novel methodology for optimizing the fiber path with a variable-axial fiber reinforcement design by employing a novel optimization methodology, called Direct Fiber Path Optimization (DFPO). The main achievement is the local optimization of both fiber angle and thickness at each finite element along the base mesh in order to reach global optimum. DFPO demonstrated its capabilities on the optimization of both open-hole and narrow-middle examples under uniaxial tension. For both cases, the results show a clear increase in both stiffness and strength compared to a reference design with equidistant straight fiber-reinforced parallel fibers, as well as compared to the principal stress oriented layouts.

Data Availability
The data used to support the findings of this study are available from the corresponding and first authors (bittrich-lars@ipfdd.de) upon request.