An enhanced curvature-constrained design method for manufacturable variable stiffness composite laminates

In this paper, design strategies are developed to explore better approaches of enforcing local layer-wise curvature constraints in the optimization of variable stiffness laminates in order to ensure the manufacturability of optimized designs based on the limitations of automated fiber placement. The methods developed here aim to improve an existing approach of imposing the curvature constraint directly on the fiber angles (i.e., direct control method) and are suitable for a design framework that uses lamination parameters as primary design variables. One approach developed here, termed the indirect control method, enforces the curvature constraint indirectly with better computational efficiency through the spatial gradient of the lamination parameters. It is shown that the curvature constraint on the actual fiber angles can also be satisfied with a sufficiently stringent upper bound albeit it produces overly conservative designs. Alternatively, an enhanced approach, termed the hybrid control method, is developed by combining the direct method and a relaxed version of the indirect control method. The case studies of minimum compliance design indicate that it provides the best manufacturable design among the three methods in the context of variable stiffness laminates using lamination parameters.


Introduction
Robotic-driven manufacturing techniques for composite materials, such as Automated Fiber Placement (AFP), allow to manufacture laminates that have non-homogeneous stiffness, commonly known as variable stiffness laminates (VSLs) [1]. With the advent of these advanced manufacturing techniques, VSL composites were soon recognized as a feasible way to improve the structural response (see, e.g., [2]). Hence, in parallel with new improvements of the manufacturing technique, companion design methods are also being actively developed to exploit the material design freedom enabled by VSLs. Some examples are methods based on Genetic Algorithms to optimize for vibration response [3] and for strength [4]. For gradient-based methods, Stegmann and Lund [5] developed the Discrete Material Optimization Method to maximize the buckling load of a wind turbine. Later, this method has been applied in topology optimization with, among others, buckling, eigenfrequency or displacement constraints [6]. The aforementioned methods use a parametrization of the design space based directly on the fiber paths (using either continuously varying fiber angles or a limited set of angles).
Alternatively, other design methods start from a broader design space without specifying or constraining a priori the number of layers in the laminate, but require post-processing to translate the optimal designs into a specific number of layers and actual fiber paths within the layers. Correspondingly, these design strategies can be classified as multistep methods. In a first step, an optimal design is identified within a broad design space typically parametrized by homogenized properties such as lamination parameters [7][8][9][10][11] or polar formalism [12][13][14]. Homogenization refers here to the through-the-thickness effective properties in a relatively thin laminate as commonly encountered in aerospace structures. Subsequently, in a second step, the optimal homogenized parameters are used to retrieve physical and/or geometrical information about the optimal laminate. These multistep methods are usually efficient in terms of exploring a large design space at a minimal computational cost.
Experimental tests confirm that VSLs can indeed improve the performance of components. For example, a rectangular specimen with a cut-out made of a VSL was found to have an increase in 13% of its buckling capacity over a conventional straight fiber laminate [15,16]. More recently, Khani et al. [17] implemented an experimental test to validate that the ultimate tensile failure load of a VSL wing lower-skin with a large access hole is 35% higher than that of the composite with homogeneous (constant) stiffness.
Despite recent progress in manufacturing techniques, AFP still has limitations in terms of the variable path designs that can actually be produced. In particular, fiber paths that have sharp turns result in defects such as wrinkles, gaps and/or overlaps, which can decrease the load carrying capacity of the composite or may become a location where delamination nucleates. In order to guarantee that the minimum turning radius of the design (reciprocal of curvature) remains above a critical value (manufacturing limit), a constraint can be imposed on the fiber path curvature. However, it has been challenging to impose this type of constraint in multistep methods without a significant performance loss. This is particularly problematic in methods where the retrieval step is carried out by matching the physical and geometrical properties of the laminate to the optimal homogenized parameters without considering the mechanical performance. Since the net gain in the optimization process can be jeopardized by these type of losses, it is important to address these issues within the overall design process.
Two basic types of methods are considered here to impose the curvature constraint, namely direct and indirect control methods. In the direct control methods, the curvature constraint is imposed directly on the fiber angles [18][19][20]. In contrast, methods where a constraint is imposed on the homogenized parameters can be called ''indirect control methods" since the constraint only restricts the curvature of the fiber angles indirectly through an implicit relation between them and the homogenized properties (see [21] for a recent example). The drawback of the indirect method is that in general it cannot guarantee a priori the satisfaction of the local fiber path curvature constraint.
In the context of a multi-step method [22], the existing direct control method [20] induces large performance loss between the consecutive steps by imposing local curvature constraints on fiber angles in the intermediate step. Simultaneously, its computational cost is relatively heavy. In the present work, a new indirect control method is developed in the first step of this multi-step method. By imposing an upper bound on the spatial gradient of the lamination parameters, the effect of the curvature of the fiber paths can be considered in the beginning of the sequential optimization procedure. Its computational efficiency is relatively high since the local curvature constraints on fiber angles in multiple layers are smeared into one gradient constraint through the lamination parameters. Moreover, the performance loss induced by the local curvature constraints in the intermediate steps in the multi-step method can be effectively reduced with this method. However, this method has to be very conservative in order to successfully constrain the curvature of the fiber paths on the individual layers. As a result, an enhanced design strategy, termed the hybrid control method, is proposed by utilizing a relaxed version of the indirect control method in conjunction with the existing direct control method. The three methods (direct, indirect and hybrid) are tested with compliance minimization problems and a comparative analysis is carried out to assess their performance. The numerical results confirm that the hybrid method produce the lowest compliance with a decent computational efficiency when manufacturing constraints are considered. For simplicity, the formulation in the present work is restricted to planar composite laminates (plates) but the general strategy can be extended to a three-dimensional context (thin shells).
The structure of this work is as follows: the formulation of the optimization problem is described in Section 2 and the curvature constraint on lamination parameters for manufacturable design is introduced in Section 3. The three aforementioned methods to control the steering of the fiber path are presented in Section 4. In Section 5, a basic verification problem, namely a rectangular plate with a point load, is used to assess the numerical performance of the methods. A comparative analysis of the three different methods to apply curvature constraints is shown in Section 6. A more complex numerical test, namely a square plate with a circular hole, is included in Section 7 in order to further assess the performance of the methods under more challenging conditions. Finally, concluding remarks are given in Section 8.

Lamination parameters
The current work focuses on the compliance optimization of planar structures (plate laminates), which are described by the in-plane stiffness. One commonly-used way to homogenize the through-the-thickness stiffness properties of a composite laminate is by using Lamination Parameters (LPs) [23]. The in-plane stiffness A is parameterized using four lamination parameters where h is the thickness of the laminate, and the lamination parameters are obtained by In Eq. (2), z ¼ z=h is the normalized position in the thickness direction and h z ð Þ is the in-plane fiber angle at z. The detailed expressions for the five constant matrices (C 0 ; C 1 ; C 2 ; C 3 ; C 4 ) appearing in Eq. (1) can be found in [7]. The stiffness matrix can be expressed as a function of the lamination parameters and the thickness of the layers. Assuming that the thickness is given (e.g., constant thickness throughout the plate), then the four lamination parameters at each point may be used as the design variables for a VSL. Consistent with Eq. (2), the lamination parameters are confined within a feasible convex domain, which is given by [24] 2V 2 For a symmetric and balanced laminate, V 2 and V 4 are zero, which implies that the feasible domain simplifies to [25] 2V 2 1 À V 3 À 1 6 0; ð4aÞ For simplicity, the design space is henceforth limited to symmetric and balanced laminates, however the design strategy may be extended to a more general case.

Three-step optimization for variable stiffness composite
A multistep approach to VSL optimization, which has proven to be very versatile, is the three-step optimization method presented in [22]. The steps, and their connection to the manufacturing constraints, are as follows: Step 1 Optimize a chosen objective functional, such as compliance or buckling, with the lamination parameters as design variables. Typically, the objective functional expressed in terms of the lamination parameters is convex [10,26,27], which enables the use of efficient optimization algorithms to find a global optimum. This is the primary optimization step. However, since the actual fiber angles/paths are not immediately available from the optimal lamination parameters, manufacturability cannot be guaranteed in this step.
Step 2 Determine the fiber angles of each ply (layer) in a specific number of plies that closely match the optimal lamination parameters obtained in Step 1. This step, referred to as the angle retrieval step in the present work, is formally an inverse problem with an objective functional meant to minimize the difference between the optimal homogenized properties from Step 1 and the actual homogenized properties corresponding to a finite number of plies. Alternatively, it can be implemented to optimize the objective functional in Step 1 with respect to the fiber angles to approach the optimal solution solution obtained in Step 1. This is refered to as angle optimization (see [27]). Manufacturing constraints can be added in this step (see [19]) with an additional loss in performance compared to unconstrained properties.
Step 3 Construct the actual fiber paths using the optimal fiber angles obtained in Step 2. This is essentially a post processing step that uses methods similar to those employed to construct streamlines in computational fluid mechanics ( [28]). The final design to be manufactured by AFP is obtained after this step. In this step it is typically not required to impose a curvature constraint since this has already been enforced in the previous steps. Correspondingly, emphasis is placed on the first two steps of the method.

Two-level approximation for the VSL
The objective functionals of the optimization problems in Step 1 (primary performance objective such as structural compliance) and Step 2 (objective functional for error minimization) can be both formally expressed using the same approach based on a twolevel approximation. Correspondingly, a common general form of the two-level approximation for the VSL is introduced in this section in anticipation of the optimization strategy. Since the problem is solved numerically using a finite element discretization (see [27]), it is also convenient to introduce a discretized version of the two-level approximation for VSLs. To this end, consider a finite element discretization of a laminate which is characterized by a collection of nodes associated to the finite elements.
The first level approximation f I ð Þ for an objective function f (for either Step 1 or 2) with respect to the in-plane global stiffness A is expressed as where n represents the total number of nodes of the finite element model, A i is the 3 Â 3 in-plane stiffness matrix at the i th node (with A À1 i denoting its inverse), U i and W i are 3 Â 3 coefficient matrices of the objective functional and: denotes the Frobenius product (inner product). The in-plane stiffness depends on the laminate layout (e.g. material elastic properties, number of layers, thicknesses, orientation), all or some of which can be viewed as design variables. In the present work, attention is focused on the fiber orientations as design variables for fixed thickness and elastic properties. Symbolically, all the necessary data required to specify a given design can be collected in a vector x, the so-called design variable, and the inplane stiffness at every node i can be expressed as a function of the design, i.e., Þ is the value of first level approximation f I ð Þ at x 0 . In the previous equation, the vector gj x 0 is the gradient of f I ð Þ with respect to x evaluated at x 0 ; Hj x 0 is the Gauss-Newton part of the Hessian matrix at x 0 to ensure the positive semi-definiteness of the approximation, Dx ¼ x À x 0 , the superscript T indicates the transpose and the single dot indicates composition. More details can be found in [19].

Optimization strategy
The strategy regarding Step 1 and Step 2 of the three-step method is as follows: the overall framework employed to solve the optimization problem is Svanberg's conservative convex separable approximations (CCSA) [29]. The approach is illustrated in Fig. 1, which shows the sub-problems used in the iterative method. Mehrotra's interior point method (IPM) is used to address the Karush-Kuhn-Tucker (KKT) conditions of the sub-problems [30]. The combination of the CCSA and Mehrotra's IPM integrates the robustness of the CCSA in convergence and the rapid convergence rate of this IPM due to its predictor-corrector steps. The element used in the FEM model in Step 1 of the optimization is an 8noded serendipity element [31], where the lamination parameters are defined on the vertices. In the remainder of the text, the nodal value correlated to the design variables refers to the value at the vertices. Since Step 2 does not require FEM for analysis, the discretization is carried out using triangular elements (see e.g., [32]). Formally, the sub-problem to be solved in both Step 1 and Step 2 is as follows: subject to x 2 Feasible domain: ð7bÞ For Step 1, for the minimum compliance problem using symmetric and balanced laminates, the design variables are collected in a vec- where the local values of the lamination parameters V 1 and V 3 at a node i are collected in the vectors V 1 and V 3 , i.e., The global feasible domain for x is defined upon applying Eq. (4) at each node i.
Step 1 is the second level approximation of the compliance. Further details, including the computation of the first level approximation f I ð Þ , can be found in [27]. The stopping criterion of the CCSA in this work is as follows: the inner loop iteration stops when the duality gap of the subproblem in the IPM is equal or less than 10 À10 . The outer loop stops when the relative variation of the compliance from Eq. (5) with respect to the compliance from the finite element method is equal or less than 10 À3 . For Step 2, the objective functional is formulated in accordance with an inverse problem to retrieve the fiber angles. The corresponding first level approximation is computed from Eq. (5) using , and A Ã i is the optimal in-plane stiffness at node i obtained during Step 1. In this step, the design variable x represents the fiber angles at layer l associated to node i, i.e., , where n d is the number of design layers.
Since the final laminate is to be symmetric and balanced, each design layer represents four layers in the actual laminate: a negative of the design layer is right next to the layer, and the complete stack is symmetric. For example, if two design layers are expressed as h 1 =h 2 ½ , the actual physical laminate is h 1 = À h 1 =h 2 = À h 2 ½ s , with the subscript s indicating symmetry. Therefore, the number of design layers n d is a quarter of the total layers n l . The feasible domain for each fiber angle at each location is Àp; p ½ . The objective function f II ð Þ x ð Þ for the sub-problem in Step 2 is the second level approximation of the error between the optimized properties from Step 1 and the actual properties of a finite number of layers. For implementation purposes, the vector x in Step 2 (and also in Step 1) can be vectorized with a single index.
Additional manufacturing constraints, which affect the actual feasible domain, will be specified in more detail in Section 4.2. Details about the streamline analogy in Step 3, which pertains to the continuous fiber paths, can be found in [28] and will not be discussed here. In the next section, gradient constraints on the lamination parameters, which aim to indirectly impose a curvature constraint, will be described.

Relationship between curvature constraints and gradient constraint on the lamination parameters
The curvature of a fiber path in the l th layer of a VSL can be expressed in terms of the norm of the gradient of the corresponding fiber angle, i.e., krh l x; y ð Þk, where r is the gradient with respect to the in-plane coordinates x and y on the surface of the model. In order to control the minimum turning radius in the l th layer of a VSL, the direct implementation of the curvature constraints on the fiber angles is as follows: krh l x; y ð Þk 6 1 r min where r min is the minimum allowable turning radius. However, this constraint requires a layer-by-layer approach, which is computationally expensive for complex structures with a relatively large number of layers. Alternatively, one can implement the gradient constraints on the lamination parameters to implicitly constrain the gradient of the fiber angle. One way to achieve this is to link the constraints on the lamination parameters with that on the fiber angles through the chain rule. Consider for example the lamination parameter V 1 , which in view of Eq. (2) can be expressed as where h l ¼ z lþ1 ð Þ À z l is the normalized thickness of the l th layer.
From this relation, the gradient of V 1 x; y ð Þ can be obtained as Taking the norm of the vectors in Eq. (10) and using the Cauchy-Schwarz inequality, it follows that Inequality (11) indicates that krV 1 x; y ð Þk can act as a lower bound for the sum of krh l x; y ð Þk. To directly control the curvature in each Step 1 (e.g., minimization of structural compliance) and Step 2 (matching parameters) are solved using the conservative convex separable approximations (CCSA) with Mehrotra's predictor-corrector interior point method (IPM). The third step (construction of fibers paths) is solved using a streamline analogy.
layer, an upper bound needs to be imposed on krh l x; y ð Þk as indicated in Eq. (8). However, a method to indirectly control the minimum turning radius is to impose a sufficiently stringent constraint on the sum of krh l x; y ð Þk in order to limit the layer-wise contributions to the sum. To avoid using explicitly a parametrization that requires the fiber angles during this step, an upper bound constraint is imposed on krV 1 x; y ð Þk. To this end, the condition Eq. (8) can be used in (11), i.e., an artificially stringent upper bound on krV 1 x; y ð Þkcan be imposed with the purpose of constraining krh l x; y ð Þk as follows: where h ¼ P n l l¼1 h l ¼ 1 and d > 0 being an adjustable parameter, called the upper bound factor. The upper bound factor d is a scalar smaller than 1 to smear the effect of sin 2h l x; y ð Þ ð Þ j j in each layer in Eq. (12) into one parameter. The parameter d is chosen such that it allows to indirectly impose an upper bound on krh l x; y ð Þk at all points and all layers through a bound on the gradient of the lamination parameters.
The same upper bound factor d is used for krV 3 x; y ð Þk to simulate the effect of sin 4h l x; y ð Þ ð Þ j j in each layer. Correspondingly, the following constraint is imposed on the gradient of In general, it is not possible to choose a priori a value of d such that Eq. (8) is automatically satisfied during Step 1. However, a suitable value of d can be identified a posteriori as shown in A. A discussion on the choice of d will be provided in the sequel but for the foregoing analysis it is sufficient to remark that it is feasible to find a suitable value.

Numerical implementation of gradient constraints on the lamination parameters
As indicated in the previous section, the indirect control of the steering of the fiber path in Step 1 of the three-step optimization method requires the gradients of the lamination parameters (i.e., rV 1 x; y ð Þ and rV 3 x; y ð Þ). In the present implementation, a fournoded quadrilateral element with four Gauss points is used. For subsequent use, the basic notation and interpolations are recorded here. In particular, the value of the lamination parameters V a ; a ¼ 1; 3, at a Gauss point g is calculated by where N i;g is the shape function of the i th node evaluated at Gauss point g, and V a;i is the value of V a at the i th node. Consequently, the gradient of V a at Gauss point g is calculated by Finally, the calculation of the norm of the gradient at each Gauss point in an element is given by where w g ¼ det J is the determinant of the Jacobian matrix J of the mapping between the physical domain and the master element at the g th Gauss point. The constraints Eqs. (14) and (15) are imposed element-wise using Eq. (18) with the corresponding bound. Three different methods to apply curvature constraints will be described in the next section.

Methods to apply curvature constraints for variable stiffness laminates
The curvature constraint has been implemented in previous works directly in Step 2 in an element-wise fashion in the form of steering constraints (see e.g., [18][19][20]). This approach will be henceforth called the direct control method. The method has the advantage to guarantee that the curvature constraint is satisfied locally at each point of each layer in a VSL composite plate. However, since the constraint is imposed during the angle retrieval step (Step 2 of the three-step method), it is uncoupled from the primary objective functional of the first step. Consequently, one drawback of this approach is that it may lead to a significant loss in performance between Step 1 and Step 2 in terms of the primary objective functional.
To address this loss in performance using the direct control method, two new strategies to apply the curvature constraints are proposed in this section. The first strategy is called the indirect control method, where the gradient constraints are imposed on the lamination parameters in Step 1 only. The second strategy, referred to as the hybrid control method, combines the use of the gradient constraints on the lamination parameters in Step 1 and the steering constraints on the fiber angles in Step 2.

General three-step framework with curvature constraints
Formally, the three methods (direct, indirect and hybrid) may be expressed in a general three-step framework with curvature constraints. Each method corresponds to a distinct and characteristic set of parameters. In this framework, the Step 1 corresponds to the primary objective functional, with or without gradient constraints on the lamination parameters depending on the method. For definiteness, the primary objective in Step 1 is chosen as the compliance C of a VSL, normalized by a referential compliance C 0 that may correspond to an initial or a benchmark design. In principle the methods used here can be applied to other types of objective functionals as long as they are parametrized with lamination parameters as design variables. In Step 2 (angle retrieval step) the objective function D represents the sum over all nodes i of the difference between the optimal stiffness matrix A Ã i obtained in Step 1 and the stiffness matrix A i for a VSL with a finite number of layers and a given set of fiber angles at the same node i, i.e., The optimization problem in Step 2 is to minimize D with respect to the fiber angles on each layer and for all nodes with or without steering constraints. If the steering constraints are taken into account, these may be imposed locally in each element and each design layer as f 2 e;l 6 f 2 max , where the value of f e;l is the fiber path curvature in the e th element and the l th design layer, and f max ¼ 1=r min being the upper bound based on the minimum turning radius r min of the AFP (more details can be found in [19]). Hence, for all methods, Step 1 and 2 are formally expressed as follows: 1. The minimum compliance problem (primary optimization problem) is expressed as subject to 2V 2 1;i À V 3;i À 1 6 0 ð20bÞ 2. The angle retrieval step The distinction between the methods is through the choice of the key parameters d and f max as indicated below. In order to distinguish the different types of constraints, the term ''gradient constraint" is primarily used for a gradient constraint on the lamination parameters in Step 1, the term ''steering constraint" refers to a constraint on the gradient of the fiber angles per layer in Step 2 and the expression ''curvature constraint" is used interchangeably for both types of constraints.

Direct control method
In Step 1 of the direct control method, the normalized compliance C=C 0 is minimized with respect to the lamination parameters without constraints on their local gradients. However, rather than suppressing the constraints, it is formally equivalent to relax the gradient constraints on the lamination parameters by setting d ! 1 so that they will not be active, hence the direct control method is characterized by a sufficiently large value of d. This approach has the advantage of employing the same code as other strategies at an acceptable computational cost since the constraints do not become active. In Step 2 (angle retrieval step) the objective function D as given in Eq. (19) is minimized with the parameter f max ¼ 1=r min chosen as an upper bound on the steering constraints, i.e., based on the actual minimum turning radius r min .

Indirect control method
In the proposed indirect control method, the normalized compliance C=C 0 is minimized in Step 1 with respect to the lamination parameters with constraints on their local gradients. In this case, it is not possible to guarantee a priori that the actual fiber path curvature f remains below a critical value f max (or, equivalently, that the radius of curvature r remains above a critical value r min ). Hence, it is necessary to choose a stringent upper bound factor d 2 0; d Ã ½ where the limit value d Ã is the largest upper bound that indirectly restricts the curvatures for a manufacturable design in Step 1. The advantage of the indirect control method is that the number of constraints is relatively small to achieve manufacturable designs, i.e., one per element for the gradient constraints on the lamination parameters. Simultaneously, the steering constraints in Step 2 are not needed. However, this value d Ã is not known a priori hence the actual upper bound factor d used requires calibration. In Step 2, the difference D is minimized without a steering constraint, which formally can be seen as taking a sufficiently large upper bound f max ! 1. In this case, the steering constraints in Step 2 are not active.

Hybrid control method
The proposed hybrid control method consists of minimizing the normalized compliance C=C 0 of the VSL using the gradient constraints on the lamination parameters in Step 1 but with a relaxed upper bound factor d 2 d Ã ; 1 ½ Þand, subsequently, minimizing the difference D with local steering constraints on the fiber angles in Step 2 using f max ¼ 1=r min . Observe that d in the hybrid control method is chosen in a different range than for the indirect control method, hence the intermediate design obtained after Step 1 (primary optimization) will in general be different for the indirect and hybrid methods. Furthermore, because of this relaxed value for d, the hybrid control method requires that the steering constraints in Step 2 should be based on the actual minimum turning radius (just like the direct control method) in order to guarantee that the manufacturing requirements are satisfied.
For clarity, the formulations of the three different curvature constraint methods are summarized in Table 1. The three methods are applied to a benchmark case with the purpose of assessing the performance associated to each method.

Test case 1: rectangular plate with point load
To study the distinct methodologies to impose curvature constraints, a simple test case of a cantilever rectangular plate is considered. The dimensions of the plate are a ¼ 1:2 m, b ¼ 0:4 m. The left edge is clamped and a unit point load F is applied at the right bottom corner. The linear elastic properties of each ply are taken as E 1 ¼ 148 GPa, E 2 ¼ 9:65 GPa, G 12 ¼ 4:55 GPa, m 12 ¼ 0:3 with 1 denoting the local fiber direction. The quasi-isotropic layout is used as a reference to compare the performance of the nonhomogeneous optimal design, hence the objective is to minimize C=C QI with C being the compliance of a non-homogeneous design and C QI being the compliance of the quasi-isotropic layout.
It is worth pointing out that a verification of the finite element implementation and the (unconstrained) optimization algorithm was done by comparing the results of the current work with those presented in [25], which were obtained from an independentlydeveloped code and algorithm. The compliance obtained from the current code had a relative difference of less that 3% compared to the results in [25] using the same geometrical, material and design data. Both solutions provided similar distributions of V 1 and V 3 and small differences may be attributed to the fact that the work of [25] uses isogeometric analysis whereas the present code uses a more traditional finite element implementation.

Primary test on direct, indirect and hybrid control method
To illustrate in detail the effect of the curvature constraints with the direct, indirect and hybrid control methods, one specific case for each method is presented in this section. All the cases presented in this section are solved with n d ¼ 6 design layers using a symmetric and balanced laminate, hence with n l ¼ 24 layers. The thickness of each layer is 0:6 mm. The initial lamination parameters at each node in Step 1 are obtained with design layers oriented as 40 ; 70 ; 40 ; 10 ; 70 ; 70 ½ , the initial fiber angles in Table 1 Parameters for three different manufacturing constraint schemes.

Method Gradient constraints in Step 1
Steering constraints in Step 2 Step 2 are 10 for each layer at each node. The r min is chosen as r min ¼ 0:8 m in this test case.

Direct control method with representative minimum turning radius
The distributions of V 1 and V 3 of the optimal solution from Step 1 and Step 2 in the direct control method are shown in Fig. 2. From the figure, the distributions of both V 1 and V 3 after Step 1 have relatively steep gradients with values changing from À1 to 1 throughout a small region. After retrieving the fiber angles in Step 2, the corresponding lamination parameters vary more gradually. In terms of the optimal compliance obtained in the two steps, the optimal normalized compliance from Step 1 is 0:482. The compliance for the ''manufacturable" design obtained after Step 2 is 0:593, which implies a relative loss of performance of 23%. This is due to the fact that the design space in Step 1 is larger (unconstrained) than the design space in Step 2 (constrained).
The fiber paths of the first design layer (i.e., the outer layer) are also shown in Fig. 2 (Step 3). In general, the resulting design follows the pattern that is expected for a minimum compliance laminate whereby the fibers optimally oppose the deformation by transmitting the applied load on the right to the clamped edge on the left following a classical arch-like layout. The resulting design complies with the manufacturing constraint of a minimum turning radius, which is in fact attained locally in this design after Step 2.

Indirect control method with implicit representative minimum turning radius
In the indirect control method, it was found through numerical experimentation, as shown in Appendix A.1, that the value of d ¼ 0:08 allows to satisfy the curvature constraints at every layer without explicitly imposing a steering constraint.
The distributions of the optimal V 1 and V 3 from Steps 1 and 2 are shown in Fig. 3. The design itself is nearly homogeneous due to the effect of a small value d in Step 1. Indeed, the normalized optimal compliance in Step 1 is 0:636 whereas the normalized compliance of the ''manufacturable" design after Step 2 is 0:647, with only a small increase of 1:7%. This small difference, which may be also observed in the contour plots of V 3 between Step 1 and 2, can be ascribed to the fact that not all lamination parameters can be matched exactly with only 6 design layers (see [33]).
The fiber paths obtained in Step 3 illustrates that the indirect control method can achieve a manufacturable design, but it is in general inefficient in terms of exploiting the local curvature. Indeed, the fiber paths of the first layer (and also in other layers) are relatively straight and the smallest turning radius in the design is 1:34 m, which implies that the steering capability is not being used to its fullest extent.

Hybrid control method with representative minimum turning radius
In the hybrid control method, partially enforcing the curvature constraints in Step 1 has the advantage of coupling geometrical requirements to the primary objective functional but leaving sufficient flexibility in the design space to prevent the negative effects of the indirect control method.
The optimal V 1 and V 3 in Step 1 of the hybrid control method, with relaxed gradient constraints on the lamination parameters by d ¼ 0:8 is shown in Fig. 4. In this case, the local steering constraints in Step 2 only induce relatively small changes during the angle retrieval step. Indeed, the normalized optimal compliance from Step 1 is 0:523 whereas the normalized compliance from Step 2 is 0:549, indicating a relatively small increase of about 5%. From the fiber paths of the first layer shown in Fig. 4 (Step 3), it can be seen that the design obtained from the hybrid method also follows a classical arch-like approach for minimum compliance.

Comparison of the optimal designs for one representative minimum turning radius
The numerical results of the primary test in Section 5.1 are summarized in Table 2. From the table, the best design (smallest compliance) in Step 1 is obtained from the direct control method. However, the order of performance changes for the final design obtained after Step 2, where the hybrid control method provides the best design. The fact that the hybrid control method delivers the best design among the three methods can be ascribed to the coupling between the curvature constraints (via the gradient of the lamination parameters) and the primary objective functional. Indeed, the partial loss in performance obtained from the hybrid method after Step 1 is compensated by a less significant loss in performance in Step 2 compared to the direct control method, which enforces the geometrical curva-ture constraint in an uncoupled fashion. The indirect control method behaves overly conservative that, as shown in the table, the smallest turning radius in the design remains significantly above the r min . From this perspective, the hybrid control method can be interpreted as a ''pre-conditioner" of the design in Step 1 that mollifies the loss of performance in Step 2 while it still takes full advantage of the steering of the fiber paths.

Comparative analysis of test case 1 for distinct values of model parameters and minimum turning radius
In the previous section, it is shown that the hybrid control method provides the best design considering a curvature constraint for a selected value of r min . In this section, the comparative analysis is extended to a wide range of values with r min ¼ 0:4; 0:8; 1:2; 1:6; 2:0 f g m, n d ¼ 6 design layers and d ¼ 0:05; 0:08; 0:8 f g . An overall assessment of the performance of the direct, indirect and hybrid control methods is achieved in terms of the compliance and computational efficiency.

Comparison of the optimal compliance in Step 2 of the three methods
The optimal normalized compliance for different values of r min obtained from the three methods which satisfies the curvature constraint at the end of Step 2 is shown in Fig. 5. As a reference, the optimum compliance without manufacturing constraints, which coincides with the performance obtained after Step 1 of the direct control method, is indicated in the figure as a dashed line. A signifi- Fig. 3. Optimal lamination parameters V1 and V3 for minimum compliance after Steps 1 and 2 obtained using the indirect control method with d ¼ 0:08; r min ! 0 and 6 design layers. For illustration purposes, the fiber paths in the first layer and its balanced counterpart are also shown (Step 3). cant loss in performance between Step 1 and Step 2 for the direct control method can be observed from this figure. The ratio of the increase in compliance ranges from about 14% for r min ¼ 0:4 m to about 35% for r min ¼ 2 m. Moreover, the indirect control method provides the worst design of the three methods (highest optimal compliance) for the whole range of r min . The hybrid control method with d ¼ 0:8 generates the best design of the three methods.
The results confirm that the indirect control method is inefficient in optimizing the compliance with manufacturable designs since its constraints are overly stringent. The hybrid control method provides the best design given that the value of the upper bound factor d is appropriately chosen (as will be shown in subsequent section).

Comparison of the computational cost for the curvature constraint methods
Imposing curvature constraints in a multistep method increases the computational cost, hence it is relevant to quantify this effect within the same computational framework. From Fig. 6, it can be observed that the indirect control method requires the least amount of time to converge. The reason is that the number of gradient constraints on the lamination parameters is relatively small comparing to the layer-wised steering constraints. Also, the local steering constraints and their sensitivities, which are computationally costly, are not needed in this method.
The CPU time for the hybrid control method with d ¼ 0:8 increases further compared to the indirect control method due to the large number of local steering constraints added in Step 2. For the cases of different values of d around 0:8 in this method, the CPU time cost is affected in a non-monotonic way. However, its variation range is not sensitive to d and is bounded above by that of the direct control method. The CPU time required for the direct control method is the highest among the three methods in general because more steering constraints become active in Step 2. In addition, it was found that the number of iterations in Step 1 for the direct control method is actually larger than for the other methods, which implies that imposing gradient constraints in Step 1 actually helps to accelerate the convergence. The results in these tests indicate that the hybrid control method not only provides the best optimal compliance and but it is also more efficient in terms of computational cost compared to the existing direct control method.

Optimal upper bound factor d in the hybrid control method
In order to find the optimum range of d in the hybrid control method for the general case, five sets of numerical tests are imple-  mented here, with the value of d uniformly varied from 0:4 to 1:6 for each set. The normalized compliance obtained after Step 2 is shown in Fig. 7. From this figure, it can be observed that for the larger values of r min considered (namely r min ¼ 1:6 m; 2:0 m), the compliance obtained only changes slightly, indicating a weak dependency on d for the range of values considered in the parametric analysis. However, if there is more freedom for the steering of the fiber path, which corresponds to smaller values of r min such as 0:4 m and 0:8 m, the results become more sensitive with respect to d. The optimum value is around d ¼ 0:8 for r min ¼ 0:4 m and around d ¼ 1 for r min ¼ 0:8 m. In general, the results suggest that the optimal value of d increases with increasing values of r min . However, it has less effect on the optimal solution in this situation, because a large value of r min limits the freedom to optimize the design. Therefore, the optimum value of d should be focused on the cases with small r min . Through numerical experimentation, a value of d around 1 is recommended.  7. Test case 2: square plate with a hole under distributed load A second design case is analyzed for a square plate with a circular hole clamped on one edge and subjected to a uniformly distributed tangential load t ¼ 20 k N m À1 . The intention is to study a more complex geometry that is commonly encountered in aerospace structures. The side length a ¼ 0:4 m and the radius of the cut-out is R ¼ 0:1 m. The material properties of the laminate are the same as those in test case 1. The r min in this case is r min ¼ 0:2 m. The value of r min is chosen in accordance with the panel's dimension to allow for a more complex design with more fiber steering.

Direct control method
With the direct control method, the optimal values for the lamination parameters V 1 and V 3 from the unconstrained Step 1 and the constrained Step 2 are shown in Fig. 8. Similar to the previous test case, the values of the lamination parameters after Step 1 have relatively steep gradients; the values are then forced to have a more gradual change after imposing the steering constraints in Step 2.
For illustration purposes, the fiber path orientations of the first layer are also shown in Fig. 8 as part of Step 3 of the three-step method. In this design, the local steering constraint was found to be active in all design layers. From the fiber paths, it is possible to recognize that the load is optimally transmitted from the edge where it is applied to the rigid support on the opposite side, adapting the path to the traction-free surfaces on the cut-out and the corresponding areas of stress concentrations along the cut-out. The effect of the steering constraints in Step 2 is that the perfor-mance of the design is reduced since the normalized compliance increases from 0:613 in Step 1 to 0:698 in Step 2, which represents a change of nearly 14 %.

Hybrid control method
The indirect control method is also tested for this case, whereas the contour plots are not shown here to cut the description. As in the test case 1, the optimal normalized compliance obtained from the indirect control method is still the highest due to its over conservativeness with small d ¼ 0:12 in Step 1.
In the hybrid control method, the d is chosen as d ¼ 1:0 due to a lower r min in Step 1. The corresponding optimal values of V 1 and V 3 in Steps 1 and 2 are shown in Fig. 9 with the first design layer as a visual reference of the optimal design. In this case, the normalized optimal compliance from Step 1 is 0:632, and the compliance corresponding to the retrieved angles in Step 2 is 0:679, hence the performance loss between the two steps is relatively moderate (increase of 7:4% in compliance). The smallest turning radius found in all layers coincides with r min ¼ 0:2 m.

Comparison of the optimal designs for one representative minimum turning radius
In this section, a comparison is made between the designs for the square plate with a central hole using the three methods. The results are summarized in Table 3 for the specific case of r min ¼ 0:2 m. Similar to the results for the rectangular plate analyzed in Section 5.2, the hybrid control method generates the design with the best performance at the end of Step 2 for test case 2. Although in this case the differences between the direct and hybrid control method are smaller, the design obtained from the hybrid control method is superior in terms of performance, which again indicates that the hybrid control method is the best option for imposing curvature constraints in the three-step optimization method. For the sake of conciseness, no parametric analysis is reported here for test case 2 (i.e., only the results for r min ¼ 0:2 m are shown), but the results for other values of r min indicate a similar ranking in terms of the performance of the designs as those in Test case 1.

Conclusion
In this paper, the three-step optimization design method is enhanced to reduce the performance loss between Step 1 (optimization with respect to lamination parameters) and Step 2 (angle retrieval step) when curvature constraints are considered. A new indirect control method is proposed to constrain the fiber path curvature through the use of gradient constraints on the lamination parameters in Step 1 while carrying out Step 2 without steering constraints on the fiber curvature. A second method termed the hybrid control method, which combines a relaxed version of the indirect method with an existing method (direct control method), is also proposed to satisfy curvature constraints in designs of composite structures for AFP technology.
Through numerical experimentation in two representative examples, it is shown that the indirect control method can indeed produce feasible designs but it severely limits the design space and prevents taking full advantage of the fiber steering capacity of AFP machines. Due to its detrimental effect on the optimization capacity, the indirect control method is not recommended as a standalone strategy. However, and more importantly, numerical results demonstrate that designs obtained from the new hybrid control method outperform individual designs obtained from the direct and indirect control methods and, in addition, converge faster to a solution compared to the existing direct control method.
The improved performance of the designs obtained from the hybrid control method can be ascribed to a couple of factors, namely (1) the global coupling between curvature constraints and primary objective in Step 1 allows to partially increase the curvature throughout the structure while accounting for the transmission of loads (global redistribution of loads) and (2) the relaxation of the constraints in Step 1 provides sufficient flexibility to locally adjust the fiber path in Step 2 up to the minimum allowable turning radius, hence fully employing the fiber steering capacity. In view of these findings, the hybrid control method becomes the recommended strategy to enforce curvature constraints, at least in the context of the minimum compliance problem. Applications of this method to other commonly-used objectives, such as buckling or strength, are to be implemented in the future work.

Acknowledgements
Zhi Hong would like to thank the financial support from the China Scholarship Council CSC (No. 201406020095) for this research.

Appendix A. Feasibility of the indirect control method
In this appendix, the feasibility of the indirect control method is tested systematically by changing the value of the upper bound factor d, the minimum allowable turning radius r min and the number of design layers in the angle retrieval step n d . Discrete values in the following ranges are used: r min 2 0:4; 2:0 ½ m, d 2 0:04; 0:14 ½ and n d ¼ 1; . . . ; 6. The problem solved is the same as the one described in test case 1 with the same material properties and thicknesses; the only difference are the values of d; r min and n d . Hence, the upper bound used in Step 1 for each case is according to Eqs. (14) and (15).
To determine whether a design obtained from the indirect control method is actually feasible or not, the curvature constraint is verified a posteriori, i.e., the smallest curvature found by searching all local values of the curvature of the design obtained after Step 2 and the corresponding Step 3 are compared to the minimum allowable turning radius. The upper bounds in Eqs. (14) and (15) depend on both r min and d (for a given thickness), hence it is convenient to present the results in two distinct formats, namely feasible/infeasible regions as a function of d for one fixed value of r min and feasible/ infeasible regions as a function of r min for one fixed value of d. In both formats the results are given for various values of n d .

A.1. Feasible/infeasible range for distinct upper bound factors d in indirect control method
The range of the upper bound factor d, where the minimum turning radius can be constrained with the indirect control method is illustrated in Fig. A.10 for r min ¼ 0:8 m for distinct values of the number of design layers n d . From the figure, it can be observed that, for example, for n d ¼ 6 design layers r min ¼ 0:8, the smallest    turning radius is below the critical value when d ¼ 0:1 (infeasible) but it is above when d ¼ 0:08 (feasible). As expected, the results indicate that the smallest turning radius decreases as d increases since the curvature constraint is progressively relaxed. However, the smallest turning radius in a design does not vary monotonically as a function of the number of design layers, which reflects the non-uniqueness of the angle retrieval process (inverse problem). It can be observed from the figure that the designs obtained with only 1 or 2 design layers (n d ¼ 1; 2) are clearly within the admissible design region, i.e., the corresponding designs tend to have relatively small curvatures (large radius of curvature). As the number of design layers increases, so does the curvature (i.e., the smallest radius of curvature decreases) and eventually the design may become infeasible.
A.2. Feasible/infeasible range for distinct minimum turning radius r min in indirect control method The range of values of the minimum turning radius r min such that the design can be constrained with the indirect control method is illustrated in Fig. A.11 for d ¼ 0:05 for distinct values of the number of design layers n d . As may be observed from the figure, it turns out that, for all values considered, a feasible design can be obtain. This reflects that the value d ¼ 0:05 is sufficiently small to generate feasible designs. The corresponding design, however, has relatively small curvature (i.e., limited steering), which indicates that imposing the curvature constraint indirectly via the gradient of the lamination parameters is plausible but at the expense of severely limiting the design capacity.