A new method for simultaneous material and topology optimization of composite laminate structures using Hyperbolic Function Parametrization

This paper presents a new discrete parametrization method for simultaneous topology and material optimization of composite laminate structures, referred to as Hyperbolic Function Parametrization (HFP). The novelty of HFP is the way the candidate materials are parametrized in the optimization problem. In HFP, a filtering technique based on hyperbolic functions is used, such that only one design variable is used for any given number of material candidates. Compared to state-of-the-art methods such Discrete Material and Topology Optimization (DMTO) and Shape Function with Penalization (SFP), HFP has much fewer optimization variables and constraints but introduces additional non-linearity in the optimization problems. A comparative analysis of HFP, DMTO and SFP are performed based on the problem of maximizing the stiffness of composite plates under a total volume constraint and multiple manufacturing constraints using various loads, boundary conditions and input parameters. The comparison shows that all three methods are highly sensitive to the choice of input parameters for the optimization problem, although the performance of HFP is overall more consistent. HFP method performs similarly to DMTO and SFP in terms of the designs obtained and computational cost. However, HFP obtains similar or better objective function values compared to the DMTO and SFP methods.


Introduction
Composite materials are known to have superior strength-to-weight and stiffness-to-weight properties compared to e.g.metallic materials.These properties, along with the inherent ability to enable tailoring of mechanical behavior of structures, have been the primary motivation for choosing composites in manufacturing of products and structural components.This has historically been driven by the needs of the aerospace industry.The past two decades have seen a rapid development of new and innovative composite manufacturing technologies and an increased availability of diverse and cost-effective composite materials.This has resulted in a substantial growth in application of composites within industrial sectors such as aerospace, energy generation, infrastructure, heavy industry and automotive.
The mechanical performance of composites structures is governed by a large number of design and manufacturing parameters due to the material characteristics of composites, and complexity of the manufacturing processes.Structural optimization is often utilized in an attempt to handle all parameters, requirements and limitations in the design of composite structures.This requires an efficient parametrization scheme of the relevant variables to generate optimal and manufacturable designs.This has resulted in a significant body of work concerning composite parametrization techniques performed within the field of structural optimization over the years.
Composite parametrization techniques can be roughly sorted into three categories: Continuous Parametrization (CP), Indirect Parametrization (IP) and Discrete Parametrization (DP).
In CP, the design variables (fiber angles, ply thickness, material distribution, etc.), are treated directly as continuous variables in the optimization problem.Continuous fiber angle optimization (CFAO) uses fiber angles θ as design variables in the problem formulation [1,2].CFAO leads to a non-convex design space with multiple stationary points.Solutions generated with CFAO will be highly dependent on the initial configuration of the design variables and will likely be suboptimal both from an optimization aspect, with a high probability of obtaining a local minima, and from an engineering aspect as the optimized designs may be difficult to realize in manufacturing.Other examples of CP methods are the multiple phase optimization approach by Zhou et al. [3] and homogenized material optimization proposed by Hozić et al. [4].In both cases the initial optimization applies a smearing/ homogenization of the material properties for a fixed set of candidate materials across the design domain, with the aim of finding the optimal material distribution.Due to the effects of homogenization these methods are not well suited to handle out-of-plane loads and in general overestimate the mechanical properties of the designs.In order to generate manufacturable composite design, one or more additional optimization steps are also required for both methods.
IP uses lamination parameters, first introduced by Tsai et al. [5,6], as intermediate design variables in the optimization problem.The stiffness of a composite laminate is fully parametrized using a set of, at most, 12 linearly interdependent lamination parameters.Grenestedt and Gudmundson [7] showed that the feasible regions for the lamination parameters are convex, leading to a convex design space of the optimization problem, thus making it possible to obtain a globally optimal solution.Later, Bloomfield et al. [8] developed a method which, for a predefined set of fiber angles, generates explicit closed-form relations defining the feasible regions for the full set of 12 lamination parameters.The main disadvantages with IP is that they are in general limited to consider one candidate material in the optimization problem, and don't generate a direct description of the laminate data for the design.This makes it difficult to include e.g.composite failure criteria and design and manufacturing constraints, and as a result, IP requires additional optimization steps to generate manufacturable designs.Despite the shortcomings, IP has successfully been used in design problem considering various engineering challenges.Some recent works include Wu et al. [9], who studied the post-buckling behavior of composite plates subjected to axial compression loads where lamination parameters were used to generate the optimal design of the plate, and Stanford et al. [10] who implemented IP for mass minimization of a wing structure under aeroelastic constraints.
Discrete parametrization has gained most interest within structural optimization for parametrization of composites.With DP, many of the limitations with CP and IP techniques are overcome.In DP, the design variables are treated as discrete variables which gives a good representation of the true nature of composites.This enables the use of a wide range of design variablesincluding fiber angles, ply density, ply thickness, and candidate materialsin the optimization problem, and makes it straightforward to include failure criteria and design and manufacturing constraints, as well as extending the optimization problems to multi-material formulations.However, the main disadvantage of DP techniques is that they result in large-scale integer optimization problems.A common approach to solve such problems is to use metaheuristic optimization algorithms, some recent examples include Deveci et al. [11], Reguera and Cortínez [12], Sreehari and Maiti [13].Although metaheuristic algorithms are frequently used for stacking sequence optimization of laminates in a variety of engineering problems, they tend to be computationally expensive even for moderate sized problems and are not well suited in applications involving topology optimization as elaborated by Sigmund [14].
A better approach, which enables the use of efficient gradient-based optimization algorithms, to avoid large-scale integer optimization problems is to relax the discreteness requirement of DP and allow a continuous variation but penalize non-discrete (intermediate) values.The first example of such an approach was the Discrete Material Optimization (DMO) method introduced by Stegmann and Lund [15,16] and further improved by Hvejsel and Lund [17,18].
DMO is based on expressing the constitutive matrix of a ply as a weighted sum of candidate materials using a multi-material interpolation function for the constitutive properties, with penalization applied to design variables associated to each candidate.Bruyneel [19] presented the Shape Function with Penalization (SFP) method as an alternative to DMO.SFP uses an interpolation scheme based on shape functions of a four-node finite element, such that a fixed set of four candidates only require two design variables; in contrast to DMO which requires one design variable per candidate.The cost for reducing the number of variables is additional non-linearity in the optimization problem.A generalized version of SFP, labelled Bi-valued Coding Parametrization (BCP), was then introduced by Gao et al. [20].This method has no upper limit on the number of candidates that can be parametrized in the optimization problem.Later, Sørensen et al. [21] presented an extension to the DMO method called Discrete Material and Thickness Optimization (DMTO) including also a ply density variable in the interpolation function of DMO, thus allowing for simultaneous material and topology optimization.New versions of DMTO have been presented by Sørensen and Lund [22] and most recently Sjølund et al. [23] in which the authors present further improvements to how the ply density variables are parametrized in DMTO.Allaire and Delgado [24] presented a level-set method for shape and stacking sequence optimization of composite laminates.This method can be seen as version of the DMTO, but the optimization of the shape (topology) and the stacking sequence (material choice) of the plies is treated separately using a staggered solution approach.Kiyono et al. [25] presented a method named Normal Distribution Function Optimization (NDFO).This approach uses a normal distribution function to parametrize a discrete set of candidates using only one design variable when considering continuously varying fiber path problems.
For a more in-depth overview of the work within this area of research, the authors recommend review articles by Ghiasi et al. [26,27] and more recent articles by Xu et al. [28], Albazzan et al. [29] and Nikbakt et al. [30].
In the present work we introduce a new DP method, referred to as Hyperbolic Function Parametrization (HFP).HFP is somewhat related to NDFO by Kiyono et al. [25] in the sense that the number of design variables necessary to parametrize any given number of candidate materials is reduced to only one.A detailed description of the DMTO, SFP and HFP methods below is followed by a numerical comparison of these three parametrization methods.The comparison is based on stiffness maximization under a total volume constraint and composite manufacturing constraints of a multi-layered composite laminate plate.Such a comparison has to the best of the authors knowledge not yet been done in the literature.

Discrete parametrization of composites
Consider a composite laminate plate with a design domain Ω made up of n l plies and divided into n e finite elements (FEs).The domain Ω el is the volume occupied by the l:th layer in the e:th element.The thickness distribution of the plate is sought using a density-based approach with an element-wise constant discretization described by a variable ρ ∈ R n e .The material in each layer of each element is assumed to be linearly elastic and chosen from a discrete set with n c candidates.For example, such a set can be fiber angles {0 • , 45 • , − 45 • , 90 • } of a unidirectional (UD) composite material or any other given combination of candidate materials, i.e. metals, various types of composite materials, etc., Fig. 1 illustrates the general idea of parametrization of the candidate materials through the interpolation function (1).The material candidates are given by a variable x.The discretization of the candidate materials is typically taken to be patch-wise constant, with each element in x giving the candidates of a given layer in a number of elements aggregated together to form a patch, a description of the implemented patch formulation is given in Section 3.1.
The parametrization of the laminate is performed on the composite ply level, for which the constitutive matrix of ply l in element e is given by where matrix E c represents the constitutive properties of a candidate material.Matrix E 0 ≈ 0 is introduced to ensure that E el is positive definite for every admissible design (strictly speaking it suffices to have a non-zero E 0 in just one of the layers).Note that both E c and E 0 are symmetric and positive definite matrices.The functions υ l and ω elc in (1) are interpolation functions for the ply densities and candidate materials, respectively.
The interpolation function ( 1) is a generalized version of the function introduced by Sørensen and Lund [31] for the DMTO method.
Depending on how ω elc is defined in (1), we obtain the DMTO, SFP or HFP method described in Sections 2.2.1, 2.2.2 and 2.2.3, respectively.
Similarly, depending on how υ el is defined in (1) different ply density penalization functions are obtained.In this work we have chosen to implement the same penalization function υ el (ρ), which is described in Section 2.1, for all three parametrization methods.

Ply density interpolation
In the present work we have implemented the version of the ply density penalization function proposed by Sørensen and Lund [22], so υ el (ρ) in ( 1) is expressed as , ∀(e, l) ( using the Rational Approximation of Material Properties (RAMP) penalization scheme [32], in which the parameter q > 1 is a penalization factor. 1 The variable ρel is referred to as the physical ply density and is obtained by applying two filters to the non-physical (optimization) design variable ρ e : first an in-plane filter applied on ρ e to generate a set of intermediate design variables ρ e ; then an a out-of-plane filter applied to ρ e to generate ρel .These filters are described below.

In-plane density filter
The in-plane density filter gives a density ρ = ρ(ρ) and is taken as a standard linear density filter, see [34,35], defined by where in which R is the radius of the filter, d(e, i) the Euclidean distance between the centroids of element e and i, and N e = {j | d(e, j)⩽R} is the filtering neighbourhood of element e.The use of the in-plane density filter herein serves a dual purpose.Firstly it ensures FE mesh-independent designs.Secondly it is used to control the thickness variation of the design, a frequent design requirement imposed on composite designs.The thickness variation, which implicitly states the ply drop-off rate in a composite design, is adjusted using the filter radius R. Increasing the value of R will reduce abrupt changes in thickness variations across the design domain.However, the in-plane density filter only provides an upper limit on the thickness variation, and it is not possible to use the in-plane density filter to specify and maintain a fixed rate of thickness variation.This motivates the addition of an explicit constraint on the thickness variation in the design-problem (21) defined below.

Out-of-plane density filter
The variables ρ and ρ represent the total material distribution within the design domain and do not contain any composite ply data.To generate discrete ply data from the material distribution given by ρ, Sørensen and Lund [22] proposed the application of a filtering method previously introduced by Gersborg and Andreasen [36] for filter-based casting constraints in general topology optimization.
In the present work, Thickness filter Nr.2 (TF2) in [22] is implemented as the out-of-plane density filter.In TF2, hyperbolic functions are used as filtering operators on the intermediate variables ρ.The main purpose of the filter is to determine if a physical ply density variable, ρel , should take on the value 0 or 1 based on the value of ρ e .The out-of-plane density filter is defined as ρel = H ρ (ρ e ; s l , β) = 1 − tanh(βρ e ) + tanh(β(s l − ρ e )) tanh(βρ e ) + tanh(β(1 − ρ e )) , ∀(e, l).
(4) 1 Whether or not q > 1 gives the desired penalization hinges of course on the optimization problem formulation, see [33].

D. Hozić et al.
Here β is a parameter that controls the filter width, and are coordinate points in a normalized through-the-thickness coordinate system.Definition (5) assumes that all plies have the same thickness, so each coordinate point s l defines the center position of a ply l in the laminate.The coordinate points s l are used in H ρ as solid/void boundary threshold values for the ply density variables such that Fig. 2 illustrates H ρ applied at l = 3 for an arbitrary laminate consisting of n l = 4 plies.The vertical axis represents the values of the continuous intermediate variable ρ e , the horizontal axis represents the discrete ply density variables ρel , while the green diamond markers are the solid/void threshold values s l for each ply given by (5).The dashed blue curve represents the filtering response of ρel for a low value of β = 20, which provide a wider filter width, resulting in a higher measure of non-discreteness of the ply density variables ρel .Higher values of β = 160, represented by the solid red curve, provide a narrow filter width thus giving a more discrete filtering response of ρel as H ρ in (4) behaves more as the unit-step function in (6).Overall, this shows that β has a large impact on the non-discreteness level of the ply density variables ρel where higher values of β are desirable as H ρ in (4) give ρel closer to either 0 or 1.
Unfortunately, initializing β to a high value leads to numerical difficulties such a convergence to poor local minima in the optimization process.Therefore a continuation strategy in which β is gradually increased during the optimization process is introduced.The implementation of the continuation strategy and penalization on the physical ply density variables ρel is discussed is Section 3.4.

Candidate material interpolation 2.2.1. DMTO -Discrete Material and Topology Optimization
Multiple versions of DMTO have been proposed, the main difference being the definition of the ply density penalization function υ el in (1).
However common for all DMTO versions, as well as the preceding DMO method, is the definition of the candidate material penalization function 1), which is expressed using the RAMP penalization scheme as , ∀(e, l, c).(7) Here x elc are the candidate material design variables, each of which is associated to one candidate in the set {1, …, n c }.As the constitutive function (1) shows, the basic idea of DMTO is to express the constitutive properties for a given ply in a laminate as a weighted sum of the constitutive properties from the set of candidates with the penalization functions ( 7) used as weights.To ensure manufacturable designs with DMTO, the design variables should satisfy Condition (8a) and (8b) require all design variables in the problem formulation to take on values between 0 and 1. Condition (8c) states that only one candidate material is allowed per ply and is included in DMTObased problem formulations as explicit linear constraints to reduce the risk of obtaining non-manufacturable mixed-material designs.

SFP -Shape Function Parametrization
Bruyneel [19] proposed SFP as an alternative method to DMTO/ DMO.SFP uses shape functions of a (four-noded) Lagrange finite element to parametrize the candidate material penalization function ω elc .The SFP version of the interpolation function ( 1) is where the candidate material interpolation functions is given as 2 , ∀(e, l, c), (10) where w c are shape functions applied as weights such that each vertex of a reference element represents one candidate material, as illustrated in Fig. 3.
The four shape functions w c in (10) are defined as (a) Reference FE-element (b) SFP parametrization Fig. 3.The SFP parametrization using mapping of the Quad4 reference element. 2 Bruyneel [19] proposed that a SIMP penalization scheme should be applied to parametrize the ply density and candidate material penalization functions in SFP; however for the sake of comparison with DMTO and HFP, we have here chosen to apply the RAMP scheme.
where 1] are natural coordinates of the reference element that are used as candidate material design variables in SFP.
As formulated in [19], a disadvantage of SFP is its limitation to a set of four candidate materials, i.e. n c = 4.However, for n c = 4, SFP only requires two candidate material design variables per ply in (9) to parametrize all four candidate materials.In comparison, DMTO requires one design variable per candidate material, i.e. twice the number of design variables.In addition, the use of shape functions in the formulation of (10) means that, implicitly, the conditions set on ω elc in (8c) are satisfied for SFP.The constraint (8c) needs therefore not to be enforced explicitly in the problem formulation.Compared to DMTO, design problems based on SFP therefore require fewer design variables and constraints.

HFP -Hyperbolic Function Parametrization
In HFP a filtering technique is applied to ω elc in (1) such that, on a ply level, any given number of candidate materials n c are parametrized with only one candidate material design variable in the optimization problem.The HFP version of the interpolation function ( 1) is where ω elc (x) = ω elc (x(x(x) ) ) is the HFP version of the candidate material penalization function.Here x is the physical candidate material design variables used to evaluate E el while x and x are the intermediate and numerical variable counterparts of x.
The candidate material penalization function used in ( 12) is expressed using the RAMP penalization scheme as Here the physical variables xelc are generated by first applying a candidate material filter to the numerical variable x el to generate a set of filtered design variables x elc , which are then transformed into xelc by implementing a normalization scheme on x elc .

Candidate material filter.
The candidate material filter is based on the same principle as the out-of-plane density filter in Section 2.1.2,in that, hyperbolic functions are used as filtering operators on the numerical candidate material variables x el .The filter is implemented on a ply level to select one candidate material from a predefined set of candidates n c .The candidate material filter is defined as (14) where x elc are the layerwise filtered candidate material design variables, α controls the filter width, and each point corresponds to a candidate material.For large values of α, H x in (14) approximates a Dirac delta function such that layer-wise filtered candidate material design variables are generated according to Fig. 4 illustrates H x in ( 14) applied to a composite ply with n c = 4 candidate materials.The vertical axis represents layerwise filtered physical variables x elc , the horizontal axis represents the layerwise numerical design variables x el , and the green markers along the horizontal axis represent the points s c associated to the set of candidates {1, …, n c }.
The two curves in Fig. 4 shows x elc for two different values of α.Lower values of α, as shown by the dashed blue curve for α = 20, provide a wide filter around each point s c and can thus result in a higher degree of non-discreteness of x elc .In addition, small α-values will induce filtering interference around the points s c .This can be seen in Fig. 4 for α = 20 where e.g.x el2 = H x (x el = s 1 , s 2 , 20) ∕ = 0, leading to a nonmanufacturable mixed material design, i.e. more than one candidate is selected.Higher values of α, represented by the solid red curve for α = 160, generate a narrow filter around the points s c thus providing a more discrete filtering of x el with little filter interference.The main concern with high values of α, as indicated by the red curve in Fig. 4, is that large sections of the feasible domain for x el do not generate any filtering response by H x , i.e. x elc = 0 for all c, which can again lead to a nonfeasible filter result as none of the candidate are selected.This shows that α has a significant impact on both the non-discreteness of the filtered variables x elc and the numerical conditioning of H x in general.
To avoid numerical difficulties for large α-values in the candidate material filter, a continuation strategy is applied to gradually increase α.In addition, a normalization scheme, as proposed by Stegmann and Lund [15], is applied to the filtered variables x elc , giving the physical (normalized) variables Introducing the normalization scheme in (17) reduces the effect of filtering interference for low values of α and the overall non-discreteness level of the physical design variables.The implementation of the continuation strategy and penalization for the physical candidate material variables xelc is discussed in Section 3.4.

Optimization Problem
In the present work the HFP method in Section 2.2.3 is demonstrated on a minimum compliance problem with a number of design and composite manufacturing constraints.The same problem formulation has also been solved using the DMTO and SFP methods such that a comparison between results of the three methods can be performed.
The compliance is given by where the global displacement vector u(ρ, x) is obtained from the linear elastic finite element state equation in which F is the global load vector.For a design domain divided into n p patches (see Section 3.1), the global stiffness matrix is given by where A represents the assembly operator for positioning of local element stiffness matrices K e into the global matrix, B e is the strain-displacement matrix and E el (ρ, x) represents the effective constitutive matrix of a ply l in finite element e given by ( 12) for HFP, and ( 1) and ( 9) for DMTO and SFP, respectively.The minimum compliance problem can now be stated as Design constraint 1 (DC1), in which a e and h ply represent the elemental area and ply thickness, respectively, is a limit on the total amount of material in the design domain.Included in the problem is also two composite manufacturing constraints (MCs).MC1 is an explicit constraint on the thickness variation of the laminate presented by Sørensen and Lund [31] which is introduced to gain explicit control over this variation not provided by the in-plane density-filter (c.f.Section 2.1.1).The constraint is here formulated for a uniform mesh with equalsized square elements and the set M e gives the (at most) four immediate neighbouring of element e.To the best of the authors knowledge, this constraint has not been formulated in a mesh-independent way; doing so is however beyond the scope of the present article.MC2 limits the number of consecutive plies (n CL ) in the laminate that can have the same candidate material.All constraints in ( 21) are used for HFP, SFP and DMTO but their functional form differ.In addition, DMTO also requires condition (8c) to be included in the problem formulation as explicit constraints.For all three methods, DC1 and MC1 are linear constraints, using only the intermediate density variable ρ e , while for HFP and SFP, MC2 is a nonlinear constraint due to the candidate material parametrization functions ω elc .
More details on the design and manufacturing constraints are given in [31,22].

Patches
Sørensen et al. [21] proposed a division of the design domain into socalled patches such that all elements in a patch retain the same properties.Patches can be defined independently on both the density variables ρ and candidate material variables x in (21) thus limiting the complexity of the topology and/or candidate material distribution within the design domain; see Fig. 5.
Applying patches will significantly affect the numerical performance and solutions of the optimization problem, therefore it is important they should reflect the limitations in composite manufacturing such that the optimal results can be realized.
In the present work, patches are implemented for the candidate material design variables x in (21).The design domain Ω is divided into a number of non-overlapping patches n p ⩽n e such that where Ω el is the design domain of ply l in element e, and P p contains a unique set of element indices that belong to patches p = 1, …, n p .With the use of patches according to (22), candidate material variables on element level x el will share the candidate material design variables on patch level x pl in (21) in the sense that x el = x pl for all e ∈ P p and all l.

Design Sensitivity
The optimization problem ( 21) is solved using gradient-based optimization solvers that require design sensitivity data for both the objective function and constraints.The sensitivity data is obtained analytically by differentiation of the objective and constraint functions in (21) with respect to the density and candidate material design variables, respectively.The derivative of the objective function (18) with respect to the density design variable is given, using the state Eq. ( 19), as

P a t c h 2 P a t c h 1
where, using the chain rule, Using ( 3), ( 4) and ( 20), the partial derivatives in (24) are expressed as In (25a) the partial derivative of the interpolation function E ml is expressed using ( 12) and ( 2) as ∂υ ml (ρ) ∂ρ ml In similar fashion the derivative of the objective function (18) with respect to the candidate material design variable is given as where ∂K e (ρ, x) ∂x elc ∂x elc ∂x el (28) Using ( 14), ( 17) and ( 20) the partial derivatives in (28) become In (29a) the partial derivative of the interpolation function E el is expressed using (12) and (13) as As for the constraint functions in (21), the only non-trivial derivative is that of MC2.This constraint only depends on the design variables x and not on the displacements u(ρ, x) and is obtained straightforwardly using the chain rule.

Optimization Solver
The design problem ( 21) is solved using a variant of the Globally Convergent Methods of Moving Asymptotes (GCMMA) [37,38].To describe the method, let y = (ρ, x).Solutions to problem (21) (and its DMTO and SFP versions) are sought by solving a sequence, k = 0, 1, …, of strictly convex subproblems of the form where g (k,ℓ) 0 is a strictly convex approximation of the compliance, and g (k,ℓ) i ,i = 1,…,m, are convex approximations of the non-linear functions defining MC2 in (21).The set Y (k) is defined by the box constraints of the original problems and so-called move limits.
In order to have a globally convergent method we require that the approximations be conservative in the sense that To achieve conservatism, one or more inner iterations may be needed, i. e. ℓ > 1.Note that the linear constraints (DC1 and MC1) can be kept as they are in the subproblems without changing the global convergence properties of the method since exact approximations are conservative [39].
The approximation g (k,ℓ) 0 of the compliance is taken as the usual MMA approximation , where the coefficients p (k,ℓ) 0j and q (k,ℓ) 0j are given by (here i = 0) in which σ j depends on the asymptotes [39].For the convex approximation of the non-linear constraint MC2 we choose, for ease of implementation, an MMA approximation: Here however, the formulas (33) would lead to a dense gradient of g (k,ℓ) i (y), making it impossible to solve large-scale instances of our design problem.The reason for this is the addition of the terms κ i σ (k,ℓ) j /4 which contribute with non-zeros to the gradient of g (k,ℓ) i even when ∂f i /∂x j is zero.To preserve the sparsity of the original constraint we therefore define p (k,ℓ) ij (and similarly for q ) as where δ ij = 1 if ∂f i /∂x j can take on non-zero values and δ ij = 0 otherwise.
With this choice of coefficients, the sparsity in the gradient of MC2 is preserved.

Continuation Strategy
To ensure manufacturable designs, the physical design variables ρ and x should take on values in {0, 1} n when solving (21).Applying gradient-based optimization solvers requires the numerical variables ρ and x in (21) to be treated as continuous variables however, as stated in Sections 2 and 2.2.3 the use of the presented filtering techniques can not guarantee that discrete values of the physical design variables are obtained and therefore the RAMP penalization is applied.
To achieve an acceptable level of discreteness for the physical design variables ρ and x, a continuation strategy, in which a step-wise increase of the values for the filter parameters α and β along with the RAMP parameters p and q is implemented.The increase in penalization is performed on the basis of the non-discreteness measure of the physical design variables calculated at each converged solution of problem in (21).Following Sørensen et al. [21] we use the following measures of non-discreteness for ρ and x: where V el is the volume of ply l in element e, M dnd represents the ply density non-discreteness and M cnd the candidate material nondiscreteness.The convergence criteria of the continuation strategy is a modified version of the one suggested in [21], obtained by adding a secondary convergence condition to M dnd so that where ∊ 0/1 is the maximum allowed value of non-discreteness, index z is the current iteration step of the continuation strategy and M Δz dnd is the secondary condition for ply density non-discreteness given by The secondary condition considers the relative change of M z dnd over three consecutive continuation iterations z, and if the maximum change of M Δz dnd is below the limit ∊ 0/1 it is assumed that M z dnd has converged to its minimum value.Note that both the conditions in (36) requires M z cnd to retain an absolute value below the non-discreteness limit ∊ 0/1 to ensure that non-manufacturable multi-material designs are not generated.Algorithm 1 shows a pseudocode for the continuation strategy implemented for HFP, here M dnd controls the values of filter parameter β and RAMP parameter q associated to the density variables, while M cnd control the values of α and p associated with the candidate material variables.(continued ) In Algorithm 1, N β , N q , N α , and N p , represent the predefined sets of values for the parameters β, q, α, and p, respectively, controlled by the continuation strategy.The values of the parameters can be adjusted separately at each new continuation step, thus enabling each sequence to stop if the corresponding design variables fulfil the non-discreteness requirement or maximum level of penalization is achieved, e.g.β j = β max or p m = p max .The convergence requirement (36) of the continuation strategy is the final requirement of the optimization problem in ( 21) that needs to be satisfied in order for a design to be accepted as an optimal solution.If (36) is not satisfied a new continuation step is initiated and a new solution of problem ( 21) is sought with updated parameters according to Algorithm 1.

Numerical Examples
This section describes the numerical examples used to evaluate the HFP method.The examples are based on two geometries that are subjected to different types of boundary conditions and loads, along with different settings for the optimization problem in (21).The numerical examples are also solved using the DMTO and SFP methods described in Sections 2.2.1 and 2.2.2, respectively, so that a comparison of the results obtained by HFP, DMTO and SFP can be performed.For the numerical examples we use a unidirectional Glass Fiber Reinforced Polymer (GFRP-UD) material with material properties provided in [22] and summarized in Table 1.We consider two sets of candidate materials, where for both sets, each candidate represents a fiber orientation as given in Table 2.
The filter and penalization parameter sets used in the continuation strategy are given in Table 3.
The parameter sets for HFP were determined by conducting a parameter study with primary focus on obtaining best results.The parameter sets for DMTO are taken from [22], in which the authors proposed a set of optimal parameter values.For SFP the parameter sets were taken to be the same as for HFP.This decision was based on the observation that the performance of SFP was better using the parameter sets for HFP than those for DMTO.Note that for HFP, the initial starting value for parameter set N α in Table 3 is increased for runs with candidate material set Θ 2 to account for the increased number of candidates within the design space, see Fig. 4. The starting values of the design variables for all optimization runs is summarized in Table 4.
The convergence tolerances used by the optimization algorithm are listed in Table 5, in which ∊ 0/1 is the tolerance for non-discreteness of the design variables described in Section 3.4, while ∊ dv , ∊ obj , ∊ const are the tolerances for the design variables, objective function and constraints, respectively, used by the GCMMA solver described in Section 3.3.The sub-problems (31) are solved using the Interior Point OPTimizer (IPOPT) solver [40].
For the FE discretizations we use rectangular fully integrated secondorder, nine-noded iso-parametric Mindlin-Reissner plate elements with five degrees-of-freedom per node.For more details regarding the finite element formulation, see e.g.[41,42].

Numerical example 1 (NE1)
The first example, NE1, is based on the setup form [22] and is illustrated in Fig. 6, which shows a square plate with clamped edges subjected to a uniform normal pressure of 1.0 kPa across the upper surface.
The plate is modelled as a composite laminate with eight plies, n l = 8.A single design patch, n p = 1, is used for the entire design domain for the candidate material design variables.The limits in the manufacturing constraints in (21) are set to restrict the thickness variation between adjacent elements to n S = 1 for MC1, and for MC2 to restrict the number of consecutive plies with the same candidate material to n CL = 1.The composite plate is required to have one full ply span the entire design domain, this is achieved by setting the lower bound of the density variable box constraint to ρ e = 1/n l in (21).Two values for the in-plane density filter radius, R = 0 and R = 0.3 [m], are considered, in order to observe the effect of this filter on the final designs.A total of 20   optimization runs are performed for NE1, where the problem ( 21) is solved for both values of the in-plane density filter radius using all three parametrization methods and for two sets of candidates in Table 2 using HFP and DMTO as SFP is limited to just four candidates.The settings for the optimization runs are listed in Table 6, where NE1.01-12 are runs with n c = 4 candidates given by Θ 1 for all three parametrization method while NE1.13-20 are runs for n c = 6 candidates given by Θ 2 for HFP and DMTO.

Numerical example 2 (NE2)
The setup for NE2 is illustrated in Fig. 7.Here we consider a cantilever rectangular plate of length l = 1.0 [m] and width w = 0.2 [m] where an out-of-plane torque is applied to the free edge such that the plate is subjected to a torsion load case.The magnitude of the applied torque is The plate is modelled as a laminate with n l = 16 plies.The FE mesh consists of 200 × 20 elements.Three values for the in-plane density filter are considered, R = 0, 0.5w and w, along with two sets of candidates, thus in total fifteen (15) optimization runs are performed.In this example there is no requirement for the laminate to retain one full ply across the design domain, so ρ e = 0 in (21).Table 7 lists the optimization settings used by the optimization runs for NE2, where NE2.01-09 are runs with n c = 4 candidate materials given by Θ 1 in Table 2 for all three parametrization method while NE2.10-15 are corresponding runs for n c = 6 candidates given by Θ 2 in Table 2 for HFP and DMTO.

Results
In this section we present and discuss data from optimization runs based on the two numerical examples.The optimization problems involve a large number of parameters which may affect the performance of the parametrization methods.In our numerical experiments we have noted that the in-plane filter radius can have a fairly large impact and we have therefore chosen to study variation of that parameter.
All data presented in the following are from optimization runs which have converged to the tolerance limits for the GCMMA solver given in Table 5.

Numerical Example 1
The results for NE1 are presented in Table 8 for the 20 × 20 mesh and Table 9 for the 30 × 30 mesh, and visualized in Fig. 8 and 9 for n c = 4 and Fig. 10 and 11 for n c = 6, respectively.Tables 8 and 9 show that optimization runs performed with the same setting of the filter radius R and mesh retain similar compliance values for HFP and SFP.The same observation can be made when comparing runs with n c = 4 and n c = 6 candidates for DMTO and HFP, respectively.
For the 20 × 20 mesh, Table 8 show that HFP and SFP retains 7% higher compliance compared to DMTO for runs with R = 0, while for runs with R = 0.3, HFP and SFP have 2.4% lower compliance values compared to DMTO.Table 9 show that for all runs with the 30 × 30 mesh, HFP and SFP retain a 5% − 20% lower compliance compared to DMTO, depending on parameter settings and number of candidates.This shows that overall, HFP and SFP obtain better objective values than DMTO, indicating DMTO is more sensitive to the input parameters.
The data in Tables 8 and 9 show that in terms of function evaluations per design iteration, HFP require an average ratio of 3.7 which is slightly higher compared to DMTO and SFP that have an average ratio of 2.8.For the runs with n c = 4, SFP requires the least amount of function evaluations to obtain a solution.For 20 × 20 mesh, SFP requires 1.9 − 3.5 times less function evaluations than DMTO and 1.4 − 1.9 than HFP.For the 30 × 30 mesh, the number of function evaluations is more similar, with SFP requiring 1.2 times fewer than DMTO and 1.2 − 1.9 times fewer  9, where NE1.08 and NE1.09 using HFP and SFP require 1.2 − 1.4 times more function evaluations than corresponding run NE1.07 using DMTO.For the runs with n c = 6 we note that HFP on average require 1.4 times more function evaluations than DMTO to obtain a solution.However, with the exception of NE1.13, runs with HFP generate better designs than DMTO.We also note that the total number of function evaluations and design iterations required to generate an optimized design increase when the in-plane density filter is used, i.e. when R > 0. For all runs in Tables 8 and 9, the number of function evaluations increase between 1.4 and 3.7 times when in-plane filter is applied.
The candidate material non-discreteness measure (35b) retain very low values, ranging from 0% to 10 − 6 % for all cases, so in Tables 8 and 9 we approximate M cnd = 0.0%.These values are reached within two continuation iterations for all runs except NE1.01 in Table 8 and NE1.07 in Table 9.This shows that all three parametrization methods can produce almost fully discrete candidate material design variables with no artificial mixing of materials.The non-discreteness measure M dnd (35a) for the density variables retain values below the limit ∊ 0/1 ≤ 0.1% for optimization runs with R = 0.For R = 0.3, M dnd does not reach the tolerance limit set in Table 5 but end up between 1.0% and 1.2%, depending on the parametrization method,.The same behavior is observed for the larger mesh in Table 9 where M dnd reaches a value between 0.9% and 1.6%, indicating that the use of in-plane density filter together with manufacturing constraint MC1 has a large effect on the density non-discreteness of the designs regardless of mesh density.
The optimized solutions in Tables 8 and 9 for runs with n c = 4 are cross-ply laminates consisting of 0 • and 90 • plies, with a maximum number of plies ranging from 6 to 7 for runs with R = 0.3 to 8 plies for runs with R = 0 as specified in Table 10.Fig. 8 and 9 show that all methods yields a similar design when R = 0 and that all designs for the 30 × 30 mesh are refined versions of the designs in Fig. 8 for the 20 × 20 mesh.We also note that for R = 0.3, similar designs are generated using HFP (Fig. 8d, 9d) and SFP (Fig. 8f, 9f) while a different design is generated with DMTO (Fig. 8b, 9b).
For runs with n c = 6 candidates, Fig. 10 and 11 show that for R = 0, both HFP (Fig. 8c, 10c, 9c, 11c) and DMTO (Fig. 8a, 10a, 9a, 11a) give the same optimized solutions as for corresponding runs for n c = 4 candidates in Fig. 8 and Fig. 9.However, for runs with R = 0.3, DMTO (Fig. 8b, 10b, 9b, 11b) does not generate the same optimized solutions, while this is the case for HFP (Fig. 8d, 10d, 9d, 11d), indicating that HFP is more consistent and is able to find the same designs even when number of candidates are increased.
The results in Tables 8 and 9 shows that all optimization runs converge to a solution within four continuation iterations, except for runs NE1.01 and NE1.12 which require 22 and 13 continuation iterations, respectively.In both cases this is attributed to the optimization solver initially settling for a mixed-material design using the DMTO method which then required additional continuation iterations to find a fully discrete candidate material solution.For all the runs, the material distribution and candidate material choice of the design is essentially set during the two initial continuation iterations and subsequent iterations only yield small changes to the design in order to fulfill the nondiscreteness criteria (36) and (37).

Numerical Example 2
The results for NE2 are presented in Table 11, with the optimal stacking sequences for all optimization runs presented in Table 12.The corresponding optimized designs are visualized in Fig. 12 for runs with n c = 4 and in Fig. 13 for runs with n c = 6.
For NE2, the same overall behaviour is observed as for NE1 in Section 5.1.Like for NE1, objective function values for the optimization runs with HFP and SFP in Table 11 retain similar values for a given value or R, with an exception for runs with R = 0 where NE2.03 using SFP gives a 5% − 10% lower compliance compared to runs NE2.02 and NE2.01 using HFP and DMTO, respectively.The lowest objective function values are consistently obtained by HFP and SFP methods, which generate designs with 4% − 12% lower compliance compared DMTO for R = 0 and R = w and n c = 4 candidates.For n c = 6, corresponding runs show that HFP generates 4% − 5% lower compliance compared to DMTO.For runs with    11 show that all three parametrization method generate design with similar objective function values for both n c = 4 and n c = 6.The results in Table 11 show that the number of function evaluations per design iteration is similar for SFP and DMTO, with an average value of 2.8 for both methods, which is the same as observed for NE1 in Section 5.1, while for HFP the average value is 4.5 which is slightly higher than the corresponding value for observed in NE1 for HFP.For runs with n c = 4, Table 11 show that for NE2, non of the methods have a clear advantage in terms of computational efficiency, as the total number of design iterations and function evaluations required to find a solution changes significantly for all methods depending on the filter radius R.This is clearly illustrated by comparing NE2.05 with NE2.06 for R = 0.5w and NE2.08 with NE2.09 for R = w, respectively.Both sets of optimization runs generate designs with similar compliance, but NE2.05 using HFP requires about 2.5 times more design iterations and function evaluations than NE2.06 using SFP.The opposite is observed for NE2.09 using SFP which requires about two times more iterations and 1.3 times more function evaluations compared with NE2.08 using HFP.For runs with n c = 6, the results show an advantage for DMTO when considering computational efficiency as DMTO requires fewer function evaluations.
The non-discreteness measures for NE2 in Table 11 show that all optimization runs give fully discrete candidate design variables as M cnd ≈ 0. M dnd follows the same pattern as for NE1 where all runs with R = 0, except NE2.11, have values below the non-discreteness limit ∊ 0/1 ≤ 0.1%.For runs with R > 0 the primary tolerance limit is not fulfilled as M dnd settles on a value between 1.26% and 1.47% depending on optimization run.All runs converge to a solution within 4 continuation iterations.
Table 12 show the optimized stacking sequences for all optimization runs in Table 11.All stacking sequences contain predominantly ±45 • plies, which is expected considering the applied torque load case for NE2  Comparing the results presented in Fig. 12 we can see that all three methods generate a similar material distribution for each setting of the in-plane density radius R and number of candidate materials.The only major difference between the designs is that some runs in Fig. 12 and Fig. 13 generate mirrored designs.

Conclusions
A new discrete parametrization method labeled HFP for optimization of composite laminate structures has been proposed.The method applies a filter-based technique using hyperbolic functions to parametrize a predefined set of candidate materials.This approach enables HFP to parametrize any given number of candidates using only one candidate material design variable per ply.Compared to existing state-of-art parametrization methods such as DMTO and SFP, HFP reduces the size of the optimization problem both in terms of number of design variable and constraints, although this comes at the cost of introducing additional non-linearity in the optimization problem.
A numerical comparison based on two numerical examples performed for sets of four and six candidate materials, respectively, shows that all three methods obtain similar optimal design in terms of topology, choice of candidate material and design-variable non-discreteness.Considering objective values, HFP generates similar or better values compared to the DMTO and SFP, depending on the input parameter.Overall, the computational cost is on the same order of magnitude for all three methods, with DMTO and SFP requiring similar number of function evaluations per design iteration while requires HFP somewhat more function evaluations per design iteration compared to SFP and DMTO.
Overall, the results indicate that all three methods are highly sensitive to the choice of input parameters for the optimization problem, such as filter-and penalization parameters, number of candidate materials, and initial starting points, although the performance of HFP is more consistent compared to DMTO and SFP.

Declaration of competing interest
The authors declare that they have no known competing financial interests or personal relationships that could have appeared to influence the work reported in this paper.

Fig. 4 .
Fig. 4. Candidate material filter H x applied for n c = 4 candidate materials and α = 20, 160, the green markers represent points s c associated to the set of candidates {1, …, n c }.

Algorithm 1 :
Pseudocode of convergence requirements for the continuation strategy implemented for HFP if Mnd⩽∊ 0/1 Convergence = TRUE (continued on next page)

Fig. 6 .
Fig.6.Geometry and boundary conditions for NE1, a clamped square plate with edge length of 1.0 m, subjected to a uniform pressure of 1.0 kPa.The plate is modelled as a monolithic laminate with n l = 8 plies.

Fig.
Fig. Geometry and boundary conditions for NE2, a rectangular plate with length l = 1 m and width w = 0.2 m, subjected to a out-of-plane torque load of M τ = 1.0 kNm at the free edge.The plate is modelled as a monolithic laminate with n l = 16 plies.

Fig. 8 .
Fig. 8. Optimized results for NE1 using the 20 × 20 mesh for n c = 4. Left column represent design without in-plane density filter and right column represents designs with applied in-plane density filter.The rows represent one of the three parametrization methods; DMTO, HFP and SFP respectively.

Fig. 9 .
Fig. 9. Optimized results for NE1 using the 30x30 mesh for n c = 4. Left column represent design without in-plane density filter and right column represents designs with applied in-plane density filter.The rows represent one of the parametrization methods; DMTO, HFP and SFP respectively.

Fig. 10 .
Fig. 10.Optimized results for NE1 using the 20 × 20 mesh for n c = 6.Left column represent design without in-plane density filter and right column represents designs with applied in-plane density filter.The rows represent the DMTO and HFP parametrization methods, respectively.

Fig. 11 .
Fig. 11.Optimized results for NE1 using the 30 × 30 mesh for n c = 6.Left column represent design without in-plane density filter and right column represents designs with applied in-plane density filter.The rows represent the DMTO and HFP parametrization methods, respectively.

Fig. 12 .
Fig. 12. Optimized results for NE2 using 200 × 20 element mesh for n c = 4. Left column represent design without in-plane density filter, central column represents design with R = 0.5w and right column represents designs with R = 1.0w.

Table 1
Material properties of GRFP-UD used for the numerical examples.

Table 3
Filter and penalization parameter sets used in the continuation strategy.

Table 4
Starting values for design variables, used for all optimization runs.

Table 5
Convergence tolerances for Continuation and GCMMA.

Table 6
Optimization settings for NE1.

Table 7
Optimization settings for NE2.

Table 8
NE1 -Results for the 20 × 20 mesh, visualization of results is given in Fig.8.nC.It is the number of continuation steps.nIt is the number of outer GCMMA iterations.nF is the number of objective and constraint function evaluations.

Table 9
NE1 -Results for the 30 × 30 mesh, visualization of results is given in Fig.9.

Table 11 NE2
-Details of results visualized in Fig.12.R = 0.5w the results in Table