Stress-Constrained Design of Functionally Graded Lattice Structures With Spline-Based Dimensionality Reduction

Topology optimization of continuum structures and functionally graded materials typically focuses on maximizing stiffness or minimizing compliance, but the design of optimal material distribution with stress constraints has gained momentum in recent years. The existing literature has addressed several inherent challenges in stress-constrained optimization that are not present in purely volume-constrained compliance minimization. Examples include the highly nonlinear and nonconvex nature of the design problem, often associated with singular topologies (cf. [1,2]), and the locality of the stress constraints, which requires a balance between efficiency and accuracy (cf. [3–5]). Most of these investigations are concerned with the solid-void topology design of homogeneous structures. When designing functionally graded lattice structures with stress constraints, however, a reliable and physically meaningful prediction of yield is needed for lattice unit cells with intermediate densities between 0 and 1. While homogenization enables the interpretation and realization of intermediate material densities in terms of the effective stiffness of specific mesostructures, the homogenized or macroscale stress underpredicts the peak stress in the underlying mesostructures (here, the term mesostructure refers to the arrangement of the material into unit cells that comprise a lattice structure). The point of initial failure in the architected material depends on the detailed local geometry and material composition of the unit cell, which are not retained during the homogenization process.

This work focuses on stress-constrained multiscale design problems. It differs from the existing work in the literature which focuses on stress-constrained designs at the mesostructural level (cf. [6–8]). These studies aim to design architected materials with reduced stress concentrations to delay the point of failure in the mesostructure given any macroscopic loads. However, they do not enable the efficient prediction of when failure occurs in the mesostructure as a function of the homogenized strains, which is of primary interest to this paper. The ability to do so is needed in the multiscale design of functionally graded lattices composed of these mesostructures.

Various methods have been proposed in the literature to define yielding in materials with intermediate densities for the purpose of structural optimization. In the context of designing solid-void topologies, power-law relationships between stress and material density have been introduced as an interpolation scheme that is designed to alleviate singularities as density approaches zero [2,4]. Physics-based stress–density relationships have been derived analytically and used to optimize special classes of materials, such as rank-2 sequential laminates in two dimensions (2D) by Duysinx and Bendsøe [9] and then laminates of any rank by Allaire et al. [10]. Lipton and Stuebner [11] optimized functionally graded lattices with multiscale stress constraints that connect macroscale to local-scale stresses for fiber-reinforced structures in 2D. Stump et al. [12] obtained stress-constrained designs of functionally graded lattices with arbitrary material compositions using the averaged stress within each phase, but they did not account for the geometry of the material interfaces, which are especially important for architected materials. A generalizable yet accurate way to describe yield at intermediate densities for arbitrary mesostructural architectures has not emerged until very recently. Pasini et al. [13] used on-the-fly homogenization nested within the optimization loop to directly determine the effective stiffness and local peak stresses when designing functionally graded lattices in 2D. A more computationally scalable approach applicable to three dimensions (3D) was proposed by Cheng et al. [14]. The yield envelopes of metallic mesostructures at intermediate densities were approximated by Hill’s criterion. Surrogate models describing the material constants in Hill’s criterion as a function of density were constructed from plastic simulations of the unit cell.

An alternative way to predict yield at the mesoscale is proposed here for the design of functionally graded lattices with any unit cell architectures consisting of solid and void. The proposed measure of yield differs from Ref. [14] in that it can be determined by analyzing the available simulation data for homogenization, thus avoiding additional nonlinear simulations of the unit cell. Furthermore, it is guaranteed to be conservative to the accuracy of the unit cell simulations. It is based on the worst-case stress in the mesostructure, which was introduced in Ref. [6] in the context of designing architected materials with improved strength performance. Here, it is used to provide a conservative estimate of mesoscale yield when designing at the macrostructural level.

The stress prediction capability is part of a computationally tractable framework for designing lattice structures for stiffness and strength. Surrogate models capturing the behavior of the mesostructure as a function of density are generated to avoid the difficult problem of analyzing detailed fine-scale features when designing the coarse-scale structure. Furthermore, the coarse-scale density field is parametrized by B-splines. B-splines have been used previously to represent density in topology optimization where the authors demonstrated the ability of splines to provide automatic filtering and length scale control, produce mesh-independent results for multi-resolution design, and handle design-dependent loads [15–17]. The use of splines also provides control over the smoothness of the density field, which could help mitigate the possibility of large stress concentrations at transitions from nearly solid elements to low volume fraction elements with small strut sizes. Smooth density variation could also aid in manufacturing, where sharp changes in the structure could lead to discontinuities and high-temperature gradients in laser-based processes and a corresponding degradation in performance. Finally, operating on spline control points can result in orders of magnitude fewer design parameters than element-wise design approaches. This compact representation enables tractable numerical differentiation for gradient-based optimization. This feature is leveraged in this work to illustrate the flexibility of the approach to be easily implemented even when using the commercial analysis software that may make it tedious or infeasible for engineers to obtain analytical gradients of the simulation. However, it should be noted that analytical gradients are available for the spline representation itself, so whenever analytical gradients are readily available, they can be incorporated into the approach to realize computational savings during the optimization process. A supplemental document (available in the Supplemental Materials on the ASME Digital Collection) provides additional details on deriving the sensitivities of the spline representation and using them within the context of the types of optimization problems described in this work.

The procedure for obtaining surrogate models of stiffness and strength is presented in Sec. 2. The details of the B-spline representation are given in Sec. 3. Sections 4 and 5 describe the optimization problem formulation and results. Finally, Sec. 6 concludes with a discussion of the demonstrated results and possibilities for future work to address the weaknesses of the proposed method and make further progress in multiscale design.

The development of surrogate or reduced-order models (ROMs) to describe the homogenized mechanical response of the mesostructures as a function of mesostructural parameters is necessary for computationally tractable multiscale design [18–20]. The presence of fine-scale features in a functionally graded lattice requires an exponential increase in meshing and simulation costs if the detailed mesostructural geometries are analyzed directly in traditional approaches to computer-aided design [21]. The approach presented here mitigates the computational burden by precomputing the structural response of defined lattice structure unit cell geometries of interest and capturing the response using simple surrogate models. The constitutive properties for each element in the finite-element mesh of the macroscale lattice can then be defined by querying the cheap surrogate models for the homogenized mechanical properties of the unit cell of interest at the prescribed macroscale relative density. This approach maintains reasonable physical accuracy with much cheaper finite-element simulations.

Training data for the surrogate models are obtained via finite-element analyses (FEAs) of the unit cells, the geometries of which are represented using the geometric projection approach described in Refs. [20,22–25]. Rod primitives are used to define the struts of the unit cell by the 3D coordinates of each endpoint and the rod diameter, thus providing mesh-independent means of defining the unit cell designs. To analyze the structural performance of the unit cell, the rods are projected into an orphan mesh of 80³ 3D continuum finite elements, yielding a total of 512,000 elements. In the finite-element projection, the geometry is described implicitly by assigning appropriate material properties to each element. This is done by sampling around the centroid of each element in a spherical volume that circumscribes the cubic brick element in order to determine the fraction of the element that lies within the boundaries of rod primitives. Elements that lie completely within the rods describing the unit cell geometry are assigned a volume fraction of 1 and Young’s modulus of the constituent material of the rod. Elements that lie completely outside of the geometry are considered to be void space and assigned Young’s modulus many orders of magnitude below the constituent material, although still non-zero to maintain numerical stability of the finite-element simulation. Elements that lie partly within the geometry are assigned Young’s modulus that linearly interpolates Young’s modulus of the constituent material proportionally to the volume fraction of the element. The analysis presented here uses the same approximate interpolation approach as Ref. [23] where elements are separated into discretized volume fraction bins to make writing the analysis input file in the commercial FEA software abaqus more tractable. The accuracy of the design discretization through geometric projection and the resulting finite-element simulations increase with a finer mesh resolution, but must be balanced against the associated computational expense. The mesh resolution used here was deemed a reasonable tradeoff between accuracy and expense in light of the previous geometric projection work such as Ref. [22]. However, adaptive mesh refinement and exact interpolation approaches would yield further increases in the accuracy of the analyses and the resulting surrogate models of homogenized unit cell properties. The use of the geometric projection method to represent designs incurs an additional computational cost compared to voxelized design representations, but the drastic reduction in design parameters and enforcement of simple geometric primitives that are known to be manufacturable make it suitable for our application.

Several different rod diameters were sampled to span a range of relative densities from approximately 2% up to 100% or fully solid. The homogenized constitutive tensor was obtained at each relative density. The data were then fit using a weighted least-squares approach where the relative error of each entry was minimized instead of the absolute error to avoid washing out the results at low volume fractions. Each of the three independent elements of the cubically symmetric constitutive tensor for the unit cell geometries of interest was fit using independent third-order polynomial functions. These surrogate models were then included in the abaqus user-defined material model (UMAT) to define the material properties for each unit cell geometry.

The unit cell is assumed to yield when the peak von Mises stress in the unit cell exceeds the yield stress of the constituent material. The yield surface of a specific unit cell geometry is a function of the relative density as well as the magnitude and direction of the homogenized stress state. An isotropic yield model is proposed here based on the worst-case stress in the microstructure. This approach has the advantage of being simple yet conservative by avoiding the difficulty associated with capturing the complex anisotropic shape of the yield surface in a high-dimensional stress space. At the same time, it can be obtained from analyzing the linear elasticity simulation data available from the homogenization process discussed in Sec. 2.2.

Using Eq. (3), it is possible to show that the worst-case stress λ in fact points to the largest eigenvalue of $Tv$⁠. The readers are welcome to Ref. [6] for additional details.

This means that the local von Mises stress at any point y^e in the mesostructure is at most λ times the magnitude $‖SM‖$ of the homogenized stress S_M. A single worst-case stress measure is then determined for the entire unit cell by computing λ_Max as the largest λ(y^e) over all y^e. For a given relative density, λ_Max indicates the largest von Mises stress that can possibly be present anywhere in the unit cell as an amplification factor of the homogenized stress magnitude. It is independent of the direction of S_M. The resulting yield model is therefore isotropic. However, the worst-case stress has merit in its simplicity in that it can be described by a scalar function of the relative density. More importantly, it provides a conservative estimate of the homogenized stress magnitude, below which designers can be certain that the peak von Mises stress in the mesostructure is below the yield stress of the constituent material.

To construct a ROM for the worst-case stress, let $ρ¯$ denote the relative density. Then, for each $ρ¯j$ where there is simulation data for the constitutive ROM, the worst-case stress measure is obtained in the following manner: the matrix $Tv(ye)$ is assembled for every point y^e on the unit cell mesh; the largest eigenvalue λ(y^e) of $Tv(ye)$ is then computed; $λMax(ρ¯j)$⁠, which is the maximum λ(y^e) over all y^e on the mesostructure, is finally determined. Once all the $(ρ¯j,λMax(ρ¯j))$ pairs are gathered, the weighted least-squares procedure from Sec. 2.2 is used to fit the data. The ROM consists of a third-order polynomial and an additional basis proportional to $1/ρ¯$⁠, which is introduced to capture the rapid growth of λ_Max as $ρ¯$ goes to zero. The ROM is obtained by fitting a fourth-order polynomial to $(ρ¯j,ρ¯jλMax(ρ¯j))$ and lowering the degree in each term by 1. Figure 1 shows that the proposed form of the ROM improves the quality of the fit compared to a fourth-order polynomial ROM.

To motivate the use of B-splines, a simple 2D test was conducted to study the effect of removing the artificial penalty in the traditional solid isotropic material with penalization (SIMP) approaches to topology optimization. The artificial penalty in SIMP makes intermediate material densities inefficient, thus encouraging designs to converge to solid and void regions. This penalty is imposed because traditional subtractive manufacturing techniques produce only fully solid material. However, advanced additive manufacturing techniques have enabled the production of architected materials that can be manufactured at intermediate densities. By removing the SIMP penalty, the 2D test aimed to investigate the effect of allowing intermediate densities. A simple cantilever beam example was set up with a 20 × 32 element mesh and a downward load applied in the middle of the free end of the beam. One optimization was conducted with the typical SIMP penalty exponent of 3, and another was conducted with a linear scaling (or a penalty exponent of 1). Figure 2 shows the intuitive results that intermediate material densities are structurally efficient when linear elasticity scaling is used without an artificial penalty. The ability to effectively distribute loads over a wider area while still obeying total mass constraints is one of the primary potential benefits of multiscale design of lattice structures. This result is similar to what has previously been demonstrated in related works on variable thickness sheet design (cf. [30,31]) and free material design (cf. [32,33]). The primary distinction between the cited works and the present one lies in the use of surrogate models for the stress-constrained design of lattices, as well as the use of splines for reduced dimensionality of the design problem.

The presence of smoother material distributions is encouraging for the use of B-splines as the design representation since splines can naturally represent smooth topologies with a much smaller number of parameters than element-based approaches. Another study was conducted with an identical loading case to Fig. 2(b), but with a 10 × 10 field of control points as design parameters. Figure 3 illustrates the result of this investigation. The two designs are comparable both visually and in terms of the resulting structural compliance. For the finite-element mesh resolution of 20 × 32 used here, this translates to a reduction by more than a factor of 6 in the number of design parameters. Even larger dimensionality reduction would be achieved in cases where the analysis requires a finer mesh resolution. However, lattice unit cell geometries will all scale differently with relative density. Depending on the unit cell selected and the loading case applied, the optimal design may still have sharp changes in topology. In these cases, increasing the control point resolution for the B-spline surfaces would be necessary to converge toward the true optimal design. A control mesh refinement study is included in Sec. 5.1 to determine the appropriate number of design parameters.

Equation (13a) imposes a physically realistic upper bound of 100% on the relative density and a lower bound of 2% based on manufacturability considerations. Equation (13b) defines a material budget for the optimization where V_Max is the maximum volume fraction of the full structure. Equation (13c) is the stress constraint defined based on the worst-case stress ROM proposed in Sec. 2.3. Equation (13d) is the partial differential equation (PDE) constraint given by the structural equilibrium equation in terms of the homogenized stress and strain. Equation (13e) defines the homogenized constitutive relationship based on the ROM in Sec. 2.2. The readers are reminded that both λ_Max and ℂ_H are spatially varying quantities, which are functions of $ρ¯$⁠, which is in turn a function of P. In this paper, the material field distribution to be optimized is 2D, and this planar design field is extruded to allow the analysis of full 3D lattice structures that are of constant density in the third dimension.

The optimization problem is solved numerically using the interior-point algorithm embedded within the fmincon function in matlab. Thanks to the reduced dimensionality of the design problem, gradients of the objective and any nonlinear constraints with respect to P can be computed via finite-difference without incurring a significant cost. This feature simplifies the use of the off-the-shelf commercial software in the optimization pipeline. For the results presented here, abaqus is used to enforce the PDE constraint in Eq. (13d) by performing a structural analysis at each design iteration. Equation (13e) is implemented in UMAT for user materials. During a finite-element simulation, the values of $ρ¯$ are needed at the Gauss quadrature points of each element, where the constitutive relationship is used to compute the stress from the strain. At each quadrature point, $ρ¯$ is first determined from Eq. (8) as a function of its location and the B-spline control point values. The constitutive ROM in Sec. 2.2 is then evaluated to determine the local stress–strain relationship for the given value of $ρ¯$⁠. Likewise, for Eq. (13c), the value of $ρ¯$ at each Gauss quadrature point is used to evaluate λ_Max at the same point according to the worst-case stress ROM from Sec. 2.3. This evaluation is implemented in the post-processing user subroutine URDFIL, in which the value of $ρ¯$ and the magnitude of the homogenized stress, $‖SM‖$⁠, are queried for each quadrature point to predict the worst-case stress $λMax‖SM‖$ in the underlying mesostructure. The subsequent sections provide more detailed discussions on the numerical treatment of the density constraint in Eq. (13a) and the volume fraction calculation in Eqs. (12b) and (13b), the stress constraint in Eq. (13c), and the beam and L-bracket geometries on which the material distribution is optimized.

This section describes the implementation of Eqs. (13a) and (13b) as linear inequality constraints. For the density constraints in Eq. (13a), simple box constraints on the design parameters P are insufficient to bound the pointwise values of $ρ¯$ due to the implicit relationship between $ρ¯$ and P through the B-spline parameterization. It must therefore be enforced as an additional set of constraints.

The smoothness of $R$ at the hinge and its slope for ξ < 0 are respectively controlled by parameters α and β. The hinge function, as shown in Fig. 4, takes on a small value when Eq. (13c) is satisfied at the point p in region q. As a result, the value of S_tol in Eq. (18) controls the amount of stress constraint violation in each region. Suitable values of α, β, and S_tol, along with the number and size of the N_R regions dividing the domain, determine how tightly the stress constraint is satisfied over the structure. For the results presented in this paper, α = 600 and S_tol = β = 4 × 10⁻⁶ are chosen with a suitably large N_R such that the maximum worst-case stress does not exceed the yield stress by more than 1–2%, while ensuring smooth convergence of the optimization. It should be noted that the worst-case stress is inherently conservative and that a safety factor could be added to compensate for the small constraint violation. Alternatively, a tighter control of the local constraint satisfaction could be achieved by adaptively tightening S_tol over the course of the optimization as proposed in Ref. [34].

For this work, the regions are chosen to be those formed by intersecting the knot intervals of the B-spline representation of the material field. This is illustrated in Fig. 5. For the choice of B-spline parameterization described in Sec. 3 and more specifically Eqs. (9) and (10), densities everywhere within each region indexed by s, t are determined by the set of 3 × 3 control points with indices from s to s + 2 and from t to t + 2. As examples, Fig. 5 highlights two regions along with the control points responsible for the density variation within each region. Dividing the structure in this manner makes it easier to satisfy the stress constraint globally during optimization.

The design workflow is demonstrated with an extruded cantilever beam in 3D. The beam, shown at the top of Fig. 6, is of rectangular shape with a dimension of 30 × 5 × 2.4 m. The constituent material for the graded lattice is assumed to be aluminum, which has Young’s modulus of 68.9 GPa, Poisson’s ratio of 0.33, and a yield strength of 276 MPa. A surface traction of 25 MPa is applied over the middle 16% of the beam cross-sectional area, which is highlighted on the right side of Fig. 6. All optimization studies of the cantilever beam in this paper are subject to a volume constraint, V_Max, of 40%.

The domain and the loading exhibit anti-symmetry about the mid-plane in y, which can be leveraged to reduce the costs of the structural simulation and design. Anti-symmetry in y means that the y-displacements are fully symmetric in the top and bottom halves of the domain, whereas the x- and z-displacements are, respectively, symmetric in magnitudes but opposite in signs. The simulation of the full domain can therefore be reduced to the top half of the geometry above the symmetry plane, as long as appropriate boundary conditions are imposed on the symmetry face. The resulting domain with the applied boundary conditions is shown near the bottom of Fig. 6. On the yz-face at x = 0, standard clamped boundary conditions (ENCASTRE BC in abaqus) are applied where all displacement components are fixed at zero. On the opposite face at the free end of the beam, a downward surface traction is applied over the bottom 8% of the halved beam cross section. The anti-symmetry condition (YASYMM in abaqus) is enforced on the xz-face at the bottom of the halved beam, where the x- and z-components of the displacements are set to zero.

Symmetry in the material field is additionally enforced. Let there be i = 1, …, m control points along the length of the beam and j = 1, …, n along the height. Let j = 1, 2 further denote control points closest to the symmetry plane. Then, it can be shown by manipulating Eqs. (9) and (10) that a zero gradient of $ρ¯$ across the symmetry plane can be ensured by setting $ρ¯i,1=ρ¯i,2$ for all i = 1, …, m. This assumption also results in a total of m(n − 1) design parameters. Section 5.2 validates the anti-symmetry condition on a uniform and a compliance-minimized density distribution, respectively.

The present workflow is additionally demonstrated on an L-bracket geometry that is commonly used to benchmark stress-constrained structural topology optimization methods. In this case, the material distribution on the L-bracket is optimized without modifying the geometry, providing a comparison between the present approach and the classical solid-void design approaches. The geometry of the L-bracket is illustrated in Fig. 7 along with its dimensions and boundary conditions. The L-bracket is clamped at the top (ENCASTRE BC in abaqus). A surface traction of 100 MPa in the downward or negative y direction is applied to the area highlighted on the bottom right, which is 1 × 1 m. The geometry is discretized by elements of size 1 × 1 × 1 for structural analyses. The constituent material for the lattice is aluminum as in the cantilever beam problem.

Figure 8 illustrates the B-spline parameterization of the material field. A B-spline grid of resolution 10 × 10 is used, which overlays the L-bracket geometry. Control points in gray do not contribute to the relative density in any part of the structural domain. Thus, they are excluded from the set of design parameters. Figure 8 also shows the N_R regions for the aggregated stress constraints in the refined grid on the bottom left of the figure (also shown in purple in the electronic version). The original regions formed by the knot intervals, as shown in Fig. 5, are refined in both x and y to ensure that the size of each region is small enough to achieve tight control of the largest stress in each region. There is a total of 84 design parameters and 175 stress constraints.

To demonstrate that the proposed framework can be used to design for different mesostructures, two types of lattice unit cells are considered for the optimization problem outlined in Sec. 4: the octet truss [35] and an alternative mesostructural architecture designed to have a higher normal stiffness but lower shear stiffness compared to the octet truss [23]. Geometries of the unit cells are illustrated in Fig. 9. Expressions for the constitutive and worst-case stress ROMs are listed in Tables 1 and 2 for the two unit cell types, respectively. Figure 10 further provides a visual comparison between the properties of the two unit cell types by plotting the ROMs from Tables 1 and 2 with the relative density on the x-axis. In all cases, material constants in the homogenized constitutive relationships are normalized by Young’s modulus of the constituent material. Results in this section are organized into four parts. Section 5.1 investigates the appropriate number of control points to be used in subsequent studies. Section 5.2 verifies the anti-symmetry assumption. Section 5.3 includes the stress-constrained compliance minimization results for the octet truss. These results are then compared with compliance-minimized designs using the alternative unit cell type with and without a stress constraint in Sec. 5.4. Finally, Sec. 5.5 presents two optimized lattice designs for the L-bracket based on stress-constrained volume minimization and volume-constrained compliance minimization, respectively.

This section investigates the effect of control mesh refinement on the optimized compliance of a functionally graded structure made of octet truss unit cells. The stress constraint is not considered in this study. The optimized material distributions for three combinations of control mesh resolution are shown in Fig. 11. The control point resolutions are superimposed on the top half of the domain, on which the optimization is performed. Recall that for the quadratic B-spline material representation used in this paper, the number of knot intervals is two less than the number of control points in each direction.

It can be observed that the lowest resolution of 6 × 6 is able to capture the general regions of high and low relative densities, with optimized compliance within 8% of the optimized compliance using 14 × 8 control points. At the intermediate refinement level of 10 × 7 control points, the difference in optimized compliance is less than 1.6% of that achieved with 14 × 8 control points. The compliance of an unoptimized structure with a uniform distribution of $ρ¯=0.4$ is 3.171 × 10⁶ Nm. Therefore, in all three cases, the optimization has improved the performance of the structure by 70%. On the other hand, a larger number of design parameters results in more optimization iterations due to the need to explore a larger design parameter space. In this case, in which the design gradients are approximated with finite-differencing, increasing the number of design parameters results in a further increase in the computational cost. Using more control points to parametrize the material field thus offers diminishing returns in designing these functionally graded structures. Based on the findings of this investigation, subsequent design studies will be conducted with a control point resolution of 10 × 7 to achieve a smooth gradation in density with enough fidelity to improve structural performance while keeping the dimensionality of the design space low.

This section verifies the anti-symmetry assumption by comparing the x- and z-displacement fields on the halved domain with those obtained via simulations on the full beam geometry. This study is conducted on a uniform density distribution of 0.4 with results shown in Fig. 12 and on the compliance-minimized density distribution with 10 × 7 control points with results shown in Fig. 13. Contours in the first column are displacements on the xy-face at z = 2.4, whereas contours in the second column are displacements on the yz-face at x = 30. For both density fields, the x- and z-displacements are equal in magnitude and opposite in sign, thus validating the anti-symmetry in the solution fields. Furthermore, the zero-contour of the displacement fields on the full domain is aligned with the location of the symmetry plane in all cases. Therefore, the anti-symmetric boundary conditions enforced on the symmetry face of the halved domain are also validated. Finally, a node-by-node comparison of all components of the displacement shows that the simulation results on the full domain can be reproduced using those on the halved domain with appropriate sign changes to numerical precision.

This section investigates various stress-constrained designs using the octet truss unit cell and compares their performance to compliance-minimized designs without the stress constraint. The optimal material distribution for minimal compliance is independent of the magnitude of the applied load at the free end of the beam. This independence is due to the linear scaling between the load and the displacements over the linear elastic domain. In contrast, a stress-constrained design is load-dependent. For this reason, two different load cases are considered: an end load of 25 MPa and 50 MPa, respectively. Figure 14 shows the worst-case stress distribution, normalized by the yield stress, on the compliance-minimized design from Fig. 11(b). In Fig. 14(a), there is a small region of constraint violation where the symmetry plane intersects the free end of the beam, and the worst-case stress in the mesostructure is 1.188 times greater than the yield stress of the constituent material. When the applied load is increased to 50 MPa, regions near the top and bottom corners of the fixed end are also failing, with the worst-case stress in the mesostructure up to 2.376 times that of the yield stress.

Figures 15(a) and 16(a) show the compliance-minimized density distributions from problem formulations that include the stress constraint for the two load cases. The corresponding normalized worst-case stress distributions can be found in Figs. 15(b) and 16(b). Figure 17 additionally plots the changes in the largest aggregated stress constraint over the first 100 optimization iterations for both the 25 MPa and 50 MPa load cases. It shows that the aggregated stress constraint is satisfied within a reasonable number of iterations in both cases. It furthermore shows that the constraint stops decreasing once it reaches the threshold, S_tol, in Eq. (18), which is shown as a dashed line in Fig. 17. The optimization algorithm considers the design feasible when the value of the aggregated stress constraint is less than the threshold, which explains why the worst-case stress slightly exceeds the yield stress in certain regions in Figs. 15(b) and 16(b). However, in both cases, the local constraint violation is small, and the feasibility of the structure has been improved by incorporating the aggregated stress constraint as indicated by the reduction in maximum worst-case stress.

For the 25 MPa load case, the largest worst-case stress is reduced to roughly 1.011 times the yield stress of the constituent material. The mild stress constraint violation is due to the choice of the aggregation function discussed in Sec. 4.2. It is possible to enforce the stress constraint more tightly by varying the user-specified parameters. Doing so, however, makes the optimization problem more nonlinear and more challenging to solve. A constraint violation of 1% is reasonable given that the worst-case stress is a conservative estimate and that a factor of safety is often added in practice. The reduction in the maximum worst-case stress is accompanied by an insignificant 0.02% increase in compliance. This result is expected because only a very small portion of the structure is failing in Fig. 14(a) without the stress constraint. It is therefore not surprising that a small change in the density distribution is sufficient to ensure the structural integrity without sacrificing much of the overall structural stiffness. In the case of the 50 MPa load, the optimized material distribution in Fig. 16(a) is more distinctly different from that in Fig. 11(b). The largest worst-case stress is reduced to roughly 1.017 times the yield stress of the material. This is accompanied by a 13% increase in compliance compared to the non-stress-constrained result. Both load cases demonstrate the importance of considering structural failure in the design of functionally graded lattice structures, as well as the ability of the B-spline material parameterization to produce designs satisfying different functional requirements.

This section investigates functionally graded lattice designs with the alternative unit cell type described in Fig. 6(b) and Table 2. For a lattice beam with a uniform density distribution of 0.4, the compliance is 1.956 × 10⁶ Nm, which outperforms an octet truss lattice with uniform density by almost 40%. Figure 18 shows the optimized density and normalized worst-case stress distribution without enforcing the stress constraint. The optimized compliance is 8.9% higher than that of the octet truss shown in Fig. 11(b), and it is evident that the worst-case stress in many regions of the structure is above the yield stress of the constituent material. This result is not unexpected given the large bending moment experienced by the high aspect-ratio cantilever beam. Figure 19 illustrates the stress-constrained design using the alternative unit cell in terms of the density and normalized worst-case stress distributions. Figure 19(a) shows a density distribution that varies more gradually than that in Fig. 18(a). This change reduces the largest worst-case stress in Fig. 19(a) on the domain to 1.013 times the yield stress. High stress areas indicated in black in Fig. 18(b) have also been substantially reduced. This is accompanied by a 22.5% increase in the structural compliance compared to Fig. 18(a). A comparison between Figs. 11(b) and 18(a), as well as between Figs. 15(a), 16(a), and 19(a) further shows that the proposed B-spline parameterization is able to produce visually distinct material distributions and that the proposed design framework can accommodate mesostructures with different mechanical properties.

Two design studies are conducted on the L-bracket geometry with the octet truss as the mesostructure. The first is a stress-constrained volume minimization problem. Equation (12b) is used as the objective, subject to all constraints in Eq. (13) except for (13b). The design is initialized with a uniform density distribution of 40%. Figure 20(a) shows the optimized density field. The resulting volume fraction is 37.5%, the compliance is 3.348 × 10⁷ Nm, and the maximum stress constraint violation is approximately 0.66% of the yield stress. Distribution of the normalized worst-case stress in Fig. 20(b) shows that the stress constraint is satisfied in most parts of the domain. The present study of the L-bracket differs from ones in the literature (cf. [4,36]) in that changes in topology due to complete removal of materials in a region are not permitted due to the lower bound in the relative density. For this reason, the shape of the reentrant corner is unaltered and is instead reinforced via adding materials to that region. The use of B-spline parameterization further produces a smoother material field distribution compared to designing with the densities of individual elements. Despite these differences, the overall distribution of materials over the L-bracket as well as the optimized volume fraction are comparable with the results reported by similar studies in the literature.

The second study minimized the compliance with a volume constraint given by Eq. (13b), where $VMax=37.5%$ is the minimized volume fraction from the previous study. In contrast with the first study where there is no consideration of compliance, the second study does not consider the maximum stress in the objective or in the constraints. The optimized density distribution is shown in Fig. 21(a), accompanied by the normalized worst-case stress distribution in Fig. 21(b). The optimized compliance is 3.050 × 10⁷ Nm, which is 8.9% lower than the compliance from the first study as expected. The decrease in compliance comes at the cost of multiple regions at risk of structural failure, which are indicated in black in Fig. 21(b). The largest worst-case stress on the L-bracket is almost twice the yield stress of the constituent material.

The two studies in this section show that the present framework produces designs that are comparable with those from the literature. Furthermore, they demonstrate that the B-spline parameterization is applicable to material field design on non-rectangular domains.

This paper describes a computational framework for designing functionally graded lattices for stiffness and strength. The lattices are composed of unit cell mesostructures. The analysis and optimization of such multiscale structures are made computationally tractable with the use of surrogate models. More specifically, stress-constrained design optimization is enabled by ROMs characterizing yield at the mesostructural level based on a worst-case stress measure. A compact design representation that also allows for control over the smoothness of the material field is achieved via B-splines. The proposed methodology is demonstrated by optimizing functionally graded lattices made of two types of unit cell structures. It is shown that the present framework allows for the design of structurally efficient lattices able to withstand different design loads. The B-spline parameterization offers sufficient fidelity to achieve significant improvements in structural compliance and to produce visually distinct material distributions for different unit cells and for design domains of different shapes.

In terms of future work, the proposed methodology can be applied to explore multiscale lattice designs with a wider range of unit cell types and in different application settings. Extending the B-spline parameterization from 2D to a fully 3D design parameter field would enable greater control over the design domain and yield potential improvements in structural performance. Implementing a workflow with analytical gradients of the material field would also provide computational savings and potentially improved accuracy in the results. The stress prediction capability can be improved by incorporating anisotropy into the yield model such that failure of the mesostructure depends on the direction of the homogenized stress tensor. Accuracy of the yield model may also benefit from taking into account different smoothed joint geometries, which could have an impact on the point of initial failure in the mesostructure. Finally, additional improvements in the structural efficiency of the optimized lattice structures may be possible with a computational framework that concurrently designs the material field and the shape/topology with suitable manufacturability constraints.

The authors would like to thank Daniel A. Tortorelli from the Lawrence Livermore National Laboratory for his suggestions throughout this research. This work was performed under the auspices of the U.S. Department of Energy by Lawrence Livermore National Laboratory under contract DE-AC52-07NA27344 and with funding from the Defense Advanced Research Projects Agency (DARPA). The views, opinions, and/or findings expressed are those of the authors and should not be interpreted as representing the official views or policies of the Department of Defense or the U.S. government. LLNL-PROC-768641.