Optimizing Topology and Fiber Orientations With Minimum Length Scale Control in Laminated Composites

Topology optimization has gained wide acceptance as a powerful design tool for identifying highly efficient structural solutions. The goal is to optimize the distribution of a given material across a design domain with the resulting connectivity defining the structural topology [1,2]. The vast majority of research in topology optimization has assumed that the material used to manufacture the resulting structural designs is isotropic. Although appropriate for many materials and fabrication processes, the assumption is generally invalid for the case of laminated composites. A typical fiber-reinforced laminate composite is composed of a stack of anisotropic laminate sheets that are joined together. A common goal in laminate design is thus to determine the optimal orientation of each laminate sheet within a stack (the stacking sequence) and across the structural domain.

Ideally, laminate orientations would be tailored at each point in space within the structural domain, allowing the use of classical homogenization-based topology optimization methods that include an orientation angle as an independent design variable [3]. Taking such an approach, however, would likely result in optimized orientation distributions that are incredibly complex and detailed. At a larger scale, one can consider the case of the precise arrangement of continuous fibers using modern advanced fiber placement machines, motivating the development of associated continuous fiber placement design methods [4,5]. In contrast, engineers in aerospace and wind turbine industries often prescribe a discrete set of allowable orientations for each laminate sheet (e.g., 0 deg, ±45 deg, 90 deg) and a minimum sheet size over which this orientation is held constant. The corresponding design optimization problem is a large scale, discrete-linked variable problem, a challenging class of problems to solve [6].

The discrete material optimization (DMO) method [7] has arisen as an effective approach to handling the discrete material selection problem in topology optimization. The DMO approach relaxes the discrete condition and introduces a continuous independent design variable for each candidate material orientation at each location the material can be placed (typically each finite element in the discretized problem). Elemental stiffness is then expressed as a combination of the stiffnesses corresponding to candidate orientations, with the solid isotropic material with penalization (SIMP) method used to drive design variables to 0–1 magnitudes, and local constraints on each element to ensure only one orientation is used at a location [7]. Voids can also be accounted for by considering them as a candidate orientation having minimal stiffness [8].

As the resulting DMO problem is continuous, gradient-based optimizers can be used to identify optimized solutions as in classical topology optimization. However, optimizing discrete orientation at the finite element level can lead to rapidly varying orientation distributions that may still be challenging to manufacture, as well as introduce issues of solution mesh dependency, both of which are well-known challenges in topology optimization [9]. The original DMO approach circumvented these issues by predefining patches composed of multiple finite elements, each required to have the same orientation. These patch discretizations are fixed during the course of the optimization, meaning that the designer is required to define the patch sizes, shape, and locations a priori. Noting this limitation, Sørensen and Lund [10] recently proposed implementing projection filters within the DMO framework, essentially allowing the optimizer to influence these decisions.

Projection-based methods in topology optimization essentially express physical design variables, typically relative densities, as a continuous function of independent optimization variables [11]. This function is defined for each finite element over a geometric space, called the neighborhood, using some form of a regularized Heaviside step function. This step function is key to achieving clear (near-discrete) distributions of relative density, offering a distinct advantage over linear filtering techniques [12–14] that lead to blurry topological boundaries or interfaces between adjacent materials. It is the combination of the neighborhood geometry and regularization function that provide the designer with control over geometric features, which have been manipulated to achieve the minimum length of designed structural features [11], void features [15,16], or both [15–17], machining constraints [18], additive manufacturing overhang constraints [19,20], and even discrete objects such as fibers [21–24], among others. An important property of the regularized Heaviside functions in these works is that all have either positive or negative curvatures through the entire range of the variable magnitudes, an important detail that guarantees the geometric control being pursued is rigorous and quantifiable.

An alternative to these Heaviside projection formulations is so-called thresholding projection functions, referred to as volume-preserving filters in the original work [25]. These filters typically use hyperbolic tangent functions as the approximation to the step function and tend to produce clear (near-binary) representations of topology, as in standard projection. The hyperbolic tangent functions are also quite convenient for perturbing a given topology, such as in robust topology optimization approaches [26]. Sørensen and Lund [10] successfully implemented the threshold projection methodology in the context of DMO and demonstrated that orientations tended to become grouped together using this approach, without the need to predefine patches of constant orientation. They also observed that the length scale of the designed features were less than the projection radius and that mesh independence of solutions was not guaranteed. Wu et al. [27] extended this idea to consider a variation on the thresholding function and fabricated designed samples, but mesh independent solutions likewise cannot be guaranteed in the approach. These observations are not unique to DMO implementations and are ultimately a result of the properties of thresholding projection. Despite its advantages, thresholding projection alone does not guarantee strict geometric control, but rather may only be said to influence geometric properties. Although often numerical results show a relation between the radius of the projection neighborhood and the length scale of designed features, it can be readily shown that features at the scale of a single element can be designed, regardless of the projection radius, when using only the thresholding projection approach.

Herein, we look to extend the novel projection-based DMO approach presented by Sørensen and Lund [10] by replacing the thresholding projection scheme with the original projection step function of Guest et al. [11], with the aim of guaranteeing the minimum length scale of orientation distributions in a physically meaningful and convenient manner. Specifically, the designer may specify the minimum patch size, defined as a physical dimension, as well as the corresponding shape of this minimum patch. In this work, we consider projected shapes that are square instead of the widely used radial projection patches, thereby allowing tight packing of the projected orientation patches and providing a more natural connection with manufactured sheets and panels that are then assembled to form a structure. Depending on the scale of the considered structure, presented solutions could likely also be created by tow steering or other advanced fiber placement technologies where the minimum length scale would also need to be constrained. We clearly acknowledge, however, that the proposed approach does not capture the full freedom provided by such technologies. Instead, we are focusing here on structures that are fabricated through the assembly of patches having a constant orientation. We note that by using different projection radii for each orientation, the designer may prescribe different minimum sizes for holes and cavities and/or for these composite patches. Finally, we note that guaranteeing the minimum length scale also provides a path to ensuring that DMO solutions are mesh independent, as observed in classical topology optimization [9,11].

The original [7] and updated DMO formulations and algorithms are presented in Sec. 2, followed by results found when solving the benchmark minimum compliance design problems previously studied in Refs. [8,10,28].

We consider the case of optimizing the distribution of orthotropic materials with different possible strong axis fiber orientations across a structural design domain. When the void phase is also considered, the problem becomes one of simultaneously optimizing topology and fiber orientations across the topology.

In order to simultaneously optimize topology and fiber orientation, an additional variable controlling the topology could be used as in typical multi-material topology optimization problems [12,31]. Alternatively, the linear equality constraints (3) could be relaxed to linear inequality constraints to allow no candidate to be chosen [29]. Here, we instead consider the void phase as one of the candidate materials [8] such that the proposed projection algorithm can conveniently and effectively control the minimum length scale of both solid material phases and the void phase. The elasticity tensor of the void phase (e.g., jth candidate material) is set as C_j = C₀ in Eq. (5), thereby ensuring that voids do not add additional stiffness to the structure.

It is generally well known from studies in classical topology optimization that the uniform weighting function is better suited for obtaining discrete solutions than the linear weighting function.

The shape of the projection domain in HPM, often called the manufacturing building block [21] or manufacturing primitive [17,24,33], becomes very important when attempting to mimic manufacturing processes and achieve strict control of minimum length scale with multiple materials or orientations. In order for the approach to be effective, projected material shapes must lie tangent, without gaps between them and without overlaps of different material projections. This is quite similar and thoroughly explained in past projection works involving multiple materials with radial projection (Fig. 1), including two-phase length scale control [16,17] and discrete object projection [21] in density-based topology optimization. To facilitate tangency, as well as utilize a manufacturing primitive that is more consistent with composite panels, we utilize a square projected shape of side length 2r_min as shown in Fig. 3.

To ensure the minimum length scale through the entire design area, special care should be taken when using the projection method at the boundaries of the design domain. As in Ref. [21], the independent (nodal) design variables are removed or set to inactive near the boundaries of the design domain within a distance of r_min, while the elemental design variables in these regions do receive the projection from nodal design variables in neighboring regions.

As the resulting optimization problems are continuous, all problems utilize SIMP formulations and are solved using a gradient-based optimizer with sensitivities computed using the adjoint method. The results presented herein were found using the method of moving asymptotes (MMA) [34,35] with the objective function scaled such that the initial objective function value is equal to 100, following the objective function range of between 1 and 100 recommended in Ref. [35]. It should be noted that the use of MMA as given in Refs. [34,35] for solving problems with a large number of active constraints is generally inefficient but can be improved using alternate methods to solve the convex subproblem, such as with interior point methods [36]. Other gradient-based solvers can also be readily used with the DMO approach (e.g., Refs. [7,8,10]). The finite element analysis, sensitivities calculation, and the optimization are all implemented with matlab code.

A continuation method is applied to the SIMP exponent to help avoid convergence to undesirable local optima. The exponent is initially set to a magnitude of p = 1.5 and is increased by 0.5 if the optimizer has converged or a maximum number of iterations has been reached. The convergence criterion used here was a change in the objective function of 0.01%, and the maximum number of iterations before a continuation step is taken was set to 50. The maximum total number of iterations was set to 250.

The projection implementation utilized a constant β approach [37], and therefore, no continuation method was used for the magnitude of β. As discussed in detail in Ref. [37], the magnitude of β was selected as a function of a design problem r_min/h ratio, where h is the size of the finite element, in order to achieve near-discrete solutions, resulting in magnitudes of β = 10, 40, and 90 for the considered r_min to h ratios. The authors also solved the same problems using a continuation method on β without observing meaningful change in solution quality, and therefore, it is recommended that a constant β approach be used.

Finally, the use of penalty function (16) with large enough magnitude of the penalty coefficient α may be sufficient to achieve discrete solutions without the use of SIMP penalization. This idea was used in Ref. [39] with a discreteness constraint and was tested here with good success. A more rigorous and broad study is required before recommending this for the algorithm proposed here.

The proposed method is demonstrated on three benchmark maximum stiffness design problems, considering both (1) material orientation selection optimization in a solid structure and (2) simultaneous topology and material selection optimization. It should be noted that the first two presented benchmark problems assume orientation at a point is constant through the thickness, and thus, only requires one layer of design variables, while the third example considers multiple layers. The material properties used for all problems are shown in Table 1, which provides Young’s moduli in the principal directions (E11, E22, and E33), shear moduli (G12, G23, and G13), and Poisson’s ratio (υ₁₂). We note the material properties used for P.1 represent a generic orthotropic material [28], while the properties for P.2 and P.3 represent a glass fiber-reinforced polymer [8,10].

Problem P.1 is academic in nature and follows the short cantilever plate solved in Refs. [7,28] that has L = 6 m, H = 2 m, t = 0.1 m, and a uniformly distributed load P = 10 KN/m applied on the top edge of the plate (see Fig. 4(a)). Following Refs. [7,28], six different allowable fiber orientations of (0 deg, 30 deg, −30 deg, 60 deg, −60 deg, 90 deg) are permitted. The finite element analysis is solved using four-node bilinear quadrilateral plane stress elements.

Problem P.2 is a pinned plate subjected to two independent load cases ([8]; Fig. 4(b)). The dimension of the plate is L = 4 m, H = 2 m, and t = 0.5 × 10⁻³ m. The applied loads F₁ = 100 N and F₂ = 100 N are two independent unit forces acting in opposite directions as shown. The objective function is an averaged compliance of the structure under the two independent load cases. In this case, the void phase and volume constraints of the solid materials are taken into consideration, and the topology and orientation within the topology are simultaneously optimized. The volume constraint for this case is set as 70% of the full design domain. The candidate materials are orthotropic with four different allowable orientations of (0 deg, 45 deg, −45 deg, 90 deg) degrees, as well as a void phase. The finite element analysis is also solved using four-node bilinear quadrilateral plane stress elements.

Problem P.3 is a clamped laminated composite shell that has L = 10 m, H = 10 m, and t = 0.05 m (Ref. [10]; Fig. 4(c)). The applied uniformly distributed load is P = 0.1 KN/m². The candidate materials for all layers are an orthotropic with four different allowable orientations of (0 deg, 45 deg, −45 deg, 90 deg) degrees and a void phase. The finite element analysis is solved using an eight-node shell element with five degrees-of-freedom per node developed specifically for laminated composite shell structures (see Refs. [41,42] for details).

Figure 5 shows the orientation distributions achieved when solving problem P.1 using a finite element mesh of 60 × 20, totaling 7200 elemental orientation design variables, with and without projection. As can be seen from Fig. 5(a), where no projection is used, the design variable orientations are permitted to vary rapidly from element to element due to the elemental orientation being an independent design variable. This not only leads to mesh dependent solutions due to the ill-posedness of the minimum compliance problem but also produces regions of rapidly varying material orientations known as checkerboard patterns (Fig. 5(a)), a well-known instability in topology optimization when using low-order finite elements [43,44]. Using higher-order finite elements may alleviate this problem but comes at the cost of increased computational expense and further does not solve the mesh dependency issue. We, therefore, turn to imposing a minimum length scale on material orientation phases to overcome these issues and provide additional manufacturability benefits.

Figure 5(b) displays the solution using Heaviside projection, and it is clear that checkerboard patterns have been eliminated and the overall orientation layouts are significantly simplified compared with the case without projection. It is observed that the minimum length scale, which is selected as r = 0.1 m for this case and is indicated by the black square, appears to be violated in the bottom right corner of the beam in Fig. 5(b). However, it is theoretically and numerically not possible to violate the definition of minimum length scale using HPM if solutions are fully discrete (M_nd = 0%). Thus, the appearance of length scale violation is due to the fact that we are plotting only the dominate material orientation phase in each finite element. In this corner of the beam, the element orientation design variables are not discrete and have not changed substantially from the initial guess, ultimately converging to multiple intermediate magnitudes per element. This is indicated by the global measure of M_nd = 4.75% for this design and confirmed by examining the distributions of elements in this region. For example, the single element in the bottom right corner has converged to a solution with all six candidates having essentially the same intermediate magnitude for x_iel (five candidates having 0.16708, and one candidate having 0.16707).

Although the relaxed DMO formulation enables discrete solutions, it relies on the exponent penalty in SIMP to motivate convergence to discrete distributions. Strain energies in the bottom right corner of the beam are negligible, and thus, the SIMP penalization has limited effectiveness in this region, particularly since material must be present (voids are not one of the candidate phases). It is theoretically possible that increasing the exponent to extreme values would help, but we instead penalize the existence of intermediate magnitudes directly in this case by using the penalty term of Eq. (16). Using this additional penalty drives the solution to the nearly discrete solution shown in Fig. 5(c) having M_nd = 0.37%, with only a minor increase in compliance. Additionally, we note that the minimum length scale is now satisfied at all locations. Due to the existence of this negligible strain energy region in the bottom right corner of the beam, and the fact that we do not allow voids as a material choice, we only use this quadratic penalty for all solutions related to this problem P.1.

The convergence history for the solution of Fig. 5(c) is shown in Fig. 6. This plot illustrates that the normalized structural compliance and discreteness consistently decline during the optimization process. Examining this plot closely reveals a small jump in the objective function at the 50th iteration, as well as less noticeable jumps in later iterations. These small jumps are due to the continuation steps where the SIMP penalty coefficients are increased. It is shown that the non-discreteness measure of the design consistently decreases from the initial guess, eventually converging to zero.

The effect of changing the minimum allowable length scale of patches of like orientation is illustrated in Fig. 7. It is clearly seen that the minimum length scale is satisfied in all cases and that increasing the minimum length scale leads to simpler orientation distributions. For the largest length scale case (Fig. 7(c)), only three of the candidate six orientations are used in the final design. The “cost” of this simplified design is an increase in compliance. The designer is thus trading improved manufacturability for a loss in structural stiffness.

Finally, we note that one can very easily assign different minimum allowable length scales to the different candidate materials by simply changing the length scale and/or shape of the materials’ projection. This is illustrated in Fig. 8 where two different sets of minimum length scales are prescribed to subsets of candidate orientations as indicated by the black squares. The minimum features for all fiber orientation phases are strictly satisfied in the optimized design. Although this likely has limited use in the setting of square orthotropic panels, it is very useful if the different candidate materials are manufactured using different processes [16], or, in the case of the next example, if inserted holes and cutouts have different minimum size requirements than structural features.

Figure 9 shows the optimized topology and orientation distributions when solving problem P.2 using a finite element mesh of 80 × 40, totaling 16,000 elemental design variables, for various magnitudes of minimum allowable length scales. This figure clearly indicates that the minimum length scale of all orientation patches and holes satisfies the prescribed minimum length scales indicated on each topology. The latter is due to the treatment of the void phase as a zero-stiffness candidate material in the DMO formulation. As Heaviside projection is applied to each of these candidate materials independently, we have direct minimum length scale control over all material phases, including voids, and orientations. Comparing Figs. 9(a) and 9(b), it is clear that the small orientation patches of 0 and 90 deg in Fig. 9(a) are eliminated when the required minimum length scale of material is increased. Additionally, the feature in the center of the structure (having 90 deg orientation) changes shape as the length is increased from Figs. 9(a) to 9(c), in order to achieve the minimum length scale in the patch as well as the patches above and below this feature. Figures 9(d) through 9(f) illustrate that changing the minimum length scale of the void phase (holes) leads to changes in the structure shape, while the general topology remains consistent in this example.

Figure 9 confirms that the minimum length scale is satisfied on the structural material as well as the voids and that structural topology and orientation can be simultaneously optimized. As expected, compliance of the optimized topologies increases as the required minimum length scale increases, again illustrating the tradeoff between structural performance and manufacturability. Comparing Figs. 9(a) and 9(b), as well as Figs. 9(d) and 9(b) where the same magnitude of β = 40 is used, it is seen that the non-discreteness measure of the optimized designs also increases with the length scale of the projection. This is a well-understood effect and can be mitigated by increasing the projection parameter as described in Ref. [37].

The convergence behavior for these examples is consistent with Fig. 6, which contains the convergence history for the solution of Fig. 9(d). Similar to the case of problem P.1, the jumps in structural compliance are observed due to the continuation steps on the SIMP penalty and become smaller and less noticeable as the non-discreteness decreases. It is shown that the non-discreteness measure of the design monotonically decreases from the initial guess, eventually converging to a small magnitude.

To illustrate the distribution of design variables more clearly, the magnitude of each candidate material orientation and void phase of Fig. 9(b) are separately plotted in Fig. 10 using a black and white contour. As can be seen in this figure, each element is clearly assigned only a single orientation or void phase, and there is no overlap between the individual candidate materials. It is also observed that the elements are either black or white, with very little gray appearing, indicating a (near-)discrete distribution.

As the proposed approach enforces a minimum length scale on designed features, including void, structural material, and orientation patches, the approach should be capable of circumventing the fundamental instability of solution mesh dependence. This is in contrast to using the thresholding filter, which does not provide mesh independence, as noted in Ref. [10]. Problem P.2 is now solved on coarse and fine meshes using consistent magnitudes of the allowable minimum length scale. Results are shown in Fig. 11, where it is seen that the topologies are consistent with only small changes at the material boundaries.

To see the effects of the structural volume constraint, optimized designs are displayed in Fig. 12 for problem P.2 for the case of r = 0.1 m, R = 0.05 m, and volume fractions of 100%, 80%, 60%, and 40%. It is shown that all four fiber orientation phases are involved in the 100% and 80% volume fraction solutions, while the 0 deg fiber orientation is eliminated from smaller volume fraction designs. Internal holes also develop as volume fraction decreases, eventually approaching a truss-like design for the 40% volume fraction solution (Fig. 12(d)). The minimum allowable length scales for void phase and fiber orientation phases are again strictly satisfied through the design domain.

The preceding examples use a single global volume constraint. It is very straightforward to apply individual volume constraints to each candidate material, which may be desirable for some applications, particularly if using different base materials. The optimized fiber orientation distributions for multiple volume constraints are shown in Fig. 13. For example, the global volume constraint 70% and an additional volume constraint 20% for the fiber of 0 deg would lead to a total volume constraint of 14% on this fiber orientation phase. The resulting topology in this case is different than Fig. 11(b) and offers a larger compliance (less stiff). This is fully expected, as applying individual volume constraints is more restrictive unless the constraint magnitudes exactly match the volume fractions found without using local volume constraints. The non-discreteness of all the optimized designs is very small, and the minimum allowable length scales are strictly satisfied through the design domain. It should be mentioned that not all equality volume constraints are strictly satisfied for the final design as a result of dependency between the volume constraints and minimum length scale. Observed violations, however, have been quite small. For example, the obtained volume fractions for discrete fiber orientations (0 deg, 45 deg, −45 deg, 90 deg) in Fig. 13 is (13.9%, 17.4%, 17.4%, 20.6%), which only slightly varies from the applied volume constraint of (14%, 17.5%, 17.5%, 21%).

The design of laminated fiber-reinforced structures often considers the case of multiple plies stacked together to form a composite structure. In this case, the fiber orientation and topology may vary with each stacked ply, leading to a changing design through the thickness of the composite. The optimization formulation presented in this paper is now applied to the clamped laminated composite shell of problem P.3 (Fig. 4(c)) where the topology and orientation in each ply are simultaneously optimized. It is important to note that the method as presented assumes that the thickness of each ply has been assigned a priori. Additionally, void material appearing on the inside of the structure in this context typically represents the use of a light foam or balsa wood to serve as permanent scaffolding during manufacturing, while offering little contribution to mass and stiffness [10,30]. Figure 14 displays the solution for a five-ply composite with a maximum allowable global volume fraction of 70% considering a minimum allowable length scale on the voids of R = 0.5 m and two different sets of minimum allowable length scales on the orientation patches: r = 0.50 m and r = 0.75 m, as indicated in the figure. The designs in each case are symmetric through the thickness and fiber orientations are also symmetrically distributed in each ply, which is expected as symmetric loading and boundary conditions are applied on the square structure. It is observed that the outer layers (#1 and #5) are completely or nearly solid, while the middle layer (#3) is entirely “void” (permanent scaffold material) in both cases, as expected to maximize bending stiffness. The minimum length scale of the orientation and “voids” are satisfied in all layers, and we note that increasing the designer-prescribed minimum allowable length scale leads to simpler but less stiff topologies. In this example, the compliance increases from 13.22 Nm to 13.66 Nm when increasing the minimum length scale r from 0.50 m to 0.75 m.

The convergence plots in terms of compliance and non-discreteness for the solutions of Figs. 14(a,1)–14(e,1) are plotted in Fig. 6. Again, a large jump in structural compliance occurs at the iteration of 50 due to the increment of SIMP penalty, and the jumps due to continuation become smaller and more unnoticeable, as the non-discreteness decreases. It also shows that non-discreteness measure of the design consistently decreases from the initial guess, eventually converging to small magnitude.

The design for an optimized shell structure having 10 layers of uniform thickness and a global volume fraction of 50% is shown in Fig. 15. It is noted that the middle two layers 5 and 6, shown in Figs. 15(e) and 15(j), are completely void and would again utilize a low mass material to serve as scaffolding, agreeing well with designs in Ref. [10]. The design topologies are also noted to be nearly symmetric around the midplane of the laminate due to the symmetry of the problem. As can be seen from the figures, the fiber orientation phase fraction also gradually increases from the middle layer to the face layer to maximize the structural stiffness. The minimum allowable length scale of fiber orientation phases is strictly satisfied in every layer and through the design domain. Interestingly, one could solve the same problem with 100% volume fraction (not allowing holes), and the solution improves compliance only 9% despite doubling the amount of material used. This confirms that the topology optimized solution of Fig. 15 does well at reducing mass with minimized loss in stiffness.

In this work, projection functions consistent with the original HPM [11] are used in the framework of DMO for simultaneously optimizing the structural topology and fiber orientation across the topology in the design of laminated composites. The work builds directly on the work of Sørensen and Lund [10], with the HPM projection functions now providing control over the minimum allowable length scale of all features, including patches of like orientation and holes, with the capability to apply different minimum length scales to these feature types.

Numerical examples show consistent trends with results from the literature for single-ply and multiple-ply benchmark design problems [8,10,28] and are shown to strictly satisfy designer-prescribed length scales and achieve near-discrete distributions of orientations selected from a discrete set of orientation options. Larger prescribed minimum length scales lead to simpler orientation layouts and topologies, and in some cases, reduce the number of selected orientations within an optimized solution. Ultimately, this provides the designer a path toward exploring tradeoffs in manufacturability and structural performance, suggesting that the projection-based method can play an important role in the design of composite laminates.

It is emphasized that the current work considers the case where structures are constructed from an optimized assembly of orthotropic panels (and voids) and that the interface of which is not explicitly represented. Although many of the solutions can be fabricated using fiber placement technologies such as tow steering, the proposed approach does not fully leverage the design freedom provided by such technologies. Additionally, we note that only structural stiffness is considered as the design objective and volume fractions (mass) as the functional constraints. Of course, there are many other considerations in the design of stacked ply laminates, including (for example) symmetry, balance, varying orientation between adjacent plies, and ply drop rate [45], as well as other structural properties including strength. Existing literature has provided a path toward topology optimization considering many of these properties [30,46], and thus, the presented approach should be readily extendable to many of these cases.

This work was supported by the National Aeronautics and Space Administration (NASA) Space Technology Research Institute (STRI) for Ultra-Strong Composites by Computational Design (US-COMP) (Grant No. NNX17AJ32G). Any opinions, findings, and conclusions or recommendations expressed in this article are those of the author(s) and do not necessarily reflect the views of NASA.

The datasets generated and supporting the findings of this article are obtainable from the corresponding author upon reasonable request. The authors attest that all data for this study are included in the paper. Data provided by a third party are listed in Acknowledgment.

c =

structural compliance

e =

index for finite elements

f =

nodal applied loads

h =

finite element size

i =

index for material phases

l =

index for layers

p =

SIMP penalty parameter

q =

quadratic penalty term

u =

nodal displacements

w =

weighting function

x =

elemental design variables after projection

C =

elastic constitutive matrix

E =

Young’s modulus

G =

shear modulus

K =

global stiffness matrix

V =

elemental volume

X =

coordinate vector (of node or element centroid)

$V¯$ =

allowable maximum volume

r_min =

filter radius and minimum allowable length scale

M_nd =

measure of non-discreteness

n^c =

number of candidate material phases

nel =

number of finite elements

nl =

number of layers

β =

parameter for the regularized Heaviside function

μ =

elemental magnitude of filtered (weighted averaged) independent design variables

ρ =

elemental volume fraction

υ =

Poisson’s ratio

ϕ =

independent design variables