Robust topology optimization of multi-material lattice structures under material and load uncertainties

Enabled by advancements in multi-material additive manufacturing, lightweight lattice structures consisting of networks of periodic unit cells have gained popularity due to their extraordinary performance and wide array of functions. This work proposes a density-based robust topology optimization method for meso- or macroscale multi-material lattice structures under any combination of material and load uncertainties. The method utilizes a new generalized material interpolation scheme for an arbitrary number of materials, and employs univariate dimension reduction and Gauss-type quadrature to quantify and propagate uncertainty. By formulating the objective function as a weighted sum of the mean and standard deviation of compliance, the tradeoff between optimality and robustness can be studied and controlled. Examples of a cantilever beam lattice structure under various material and load uncertainty cases exhibit the efficiency and flexibility of the approach. The accuracy of univariate dimension reduction is validated by comparing the results to the Monte Carlo approach.


Introduction
Lattice structures are periodic or aperiodic cellular structures with strut-based unit cells. Because they consist of many structural members and are naturally porous, they hold desirable properties, e.g., low weight and material cost, while exhibiting extraordinary mechanical, thermal or acoustic properties. In recent years, additive manufacturing (AM) has allowed easier manufacture of geometrically complex lattices, and as a result the structures have gained traction as both micro-architectured metamaterials and mesostructures. At the microscale, for example, lattice metamaterials have been designed for negative Poisson's ratio [1][2][3], acoustic manipulation [4,5] and energy absorption [6]. Mesoscale lattice structures have also been embedded into global structures to achieve multifunctional heat exchangers [7,8], and parts with high energy absorption [9][10][11][12] and low weight [13,14]. In particular, lattice structures are exceptionally popular in the biomedical field due to their porosity, which allows biocompability with organic tissue, high performance-toweight ratio, and ease of customization [15][16][17]. By exploiting the advancements in multi-material AM, greater performance improvements [12], like graded stiffness [2] and a larger range of stiffness [18], are possible.
This work focuses on periodic, mesoscale, multimaterial lattice structures, which can replace a portion or the entirety of a solid global part, and is also applicable to macroscale lattices. Due to the large number of lattice struts, or members, as well as the additional materials, the design space expands significantly. Thus, to design these structures, computational methods, e.g., topology optimization, which automatically lays out the materials, are highly desirable. While multi-material topology optimization has become more prevalent, most methods are not applicable to meso-and macroscale lattice structures. The color level-set [19] and multiphase phase-field [20] methods require solving the Hamilton-Jacobi and Cahn-Hilliard partial differential equations, respectively, on a continuous domain. Therefore, they cannot design lattice structures in which each strut is a discrete element. Techniques for discrete elements like the multi-material genetic algorithm [21] and pseudo-sensitivity [22] approaches exist, but are prohibitive for large problems. Zhang et al. [23] used the ground structure method to iteratively remove members, applying a new variable update method to effectively optimize multi-material trusses under multiple constraints.
Similar to Ref. [23], this work takes the strut member as an element and can either remove or assign materials to each. Departing from the above methods, however, a density-based approach is employed by calculating the elastic modulus based on the well-known 3-phase standard solid isotropic material with penalization (SIMP) scheme by Bendsøe and Sigmund [24]. It is also inspired by the modified version of SIMP that replaces the void elastic modulus with a small number to avoid singularity during finite element analysis (FEA) [25]. The formulation in this paper takes SIMP a step further by allowing multiple materials.
A few others have extended Sigmund's 3-phase interpolation schemes to M solid phases and can be broadly categorized into: 1) Schemes that add design variables for each new material, similar to standard and modified SIMP, and 2) schemes that limit the number of design variables in order to reduce cost. Stegmann and Lund [26] employed a model of the first category, initially allowing more design variables but later reducing the number using "patch" variables. They, however, followed the less stable standard SIMP formulation. Gaynor et al. [27] also proposed a first category scheme that, despite considering pratical PolyJet manufacturing constraints, unrealistically requires the difference between the elastic modulus of each material to be equal. From the second category, Yin and Ananthasuresh [28] used a peak function to transform one design variable into the modulus of multiple materials. Although it does not increase the design variables, the function contains horizontal slopes that may cause instability, and includes several curve parameters that must be meticulously selected by the user. Another method is Zuo and Saitou's ordered scheme [29], which avoids instability in its power functions but has a complicated formulation and fluctuations in the optimization history.
It should be noted that although using less design variables does reduce the time needed to evaluate the modulus and to update the variables, the vast majority of topology optimization, especially in large problems, is spent in FEA, which depends on the number of elements in the domain. As this number generally remains the same throughout optimization regardless of how many materials exist, it is not necessary to force a more complicated scheme in exchange for less design variables. Hence, this work proposes a multi-material interpolation scheme of the first category that is generalized for arbitrary numbers of materials, easy to implement and based on modified SIMP for more stability.
Furthermore, the uncertainty in the materials and loads are considered in this work. In AM, uncertainty is inherent in the material (e.g., variability in the powder of laser powder bed fusion) and the build process (which can lead to inconsistent geometry). In lattice structures, the strut dimensions and mechanical properties can vary greatly depending on the processing conditions [30,31]. Park et al. [32] quantified the variability in the strut geometry by calculating the probabilistic distribution of effective diameters, then deriving the effective elastic modulus. In this paper, a similar approach is taken, except here the elastic moduli of the materials are used directly as the uncertain parameters. By doing so, both the uncertainty in the material and the built geometry can be simultaneously considered by prescribing appropriate probability distributions. Load uncertainty that can arise during application is also added for more robust designs. Both material and load influence the part performance, which is used as the objective function in topology optimization; considering their variability can result in different and more robust designs.
Topology optimization methods that consider geometric [33][34][35] and material [36,37] uncertainties exist for single material continuous problems. There is relatively more research for random loads, such as Refs. [38][39][40]. Nonetheless, the majority of these methods use level-set topology optimization, projection filters to model uncertainty, or Monte Carlo to quantify uncertainty--all of which are difficult or impossible with multi-material lattice structures consisting of a large number of strut elements. In response, this work proposes utilizing a weighted sum of the mean and standard deviation of the performance criteria, e.g., compliance, as the objective function, where the statistical moments are calculated using univariate dimension reduction (UDR) [41] and Gauss-type quadrature sampling.
The new generalized multi-material interpolation scheme is proposed in Section 2. In Section 3, the deterministic formulation of the topology optimization using the scheme is introduced, followed by the robust formulation and details of uncertainty quantification in Section 4. Numerical examples of a mesoscale lattice structure with a global cantilever beam shape under different uncertainty cases is shown and discussed in Section 5. Finally, the conclusions are presented in Section 6.

Multi-material interpolation scheme
The proposed multi-material interpolation scheme follows the SIMP method, where the continuous design variables are the artificial densities of the elements. The material property of each element is expressed as a material interpolation function that penalizes intermediate values of the densities through an exponent, p. For example, the Young's modulus of each element e in a design containing one solid material can be described by the standard SIMP model as [24] where e 2ð0,1 is the artificial density variable of element e such that it is greater than some small number min to avoid singularity when solving the finite element problem, and E is the overall Young's modulus of the solid material.
Using the so-called modified SIMP scheme, singularity in the stiffness matrix can be circumvented by assigning a small positive number to the void phase as follows [25]: where E 1 is the elastic modulus of the solid material, E 0 is a small value that represents the void phase, and p is the penalty coefficient to discourage intermediate densities within ½0,1. This improved formulation not only allows E 0 to be independent of p but also facilitates filtering schemes such as the minimum length scale filter. Gibiansky and Sigmund [42] showed the following 3-phase (two distinct solid materials and one void phase) scheme for composites that includes a second penalization parameter, q: where E 2 is the modulus of the second solid material, and the design variables are e$1 , e$2 2 ½0,1. This model was tested by Stegmann and Lund in Ref. [26], where the two penalty parameters had to be carefully chosen in order to converge to a global optimum. Moreover, Eq. (3) follows the standard SIMP model instead of the more stable modified version in Eq. (2). In response, this work amalgamates Eq. (3) and the modified SIMP scheme for three phases, without the inclusion of another penalty parameter, as follows: From the second line, it is clear that the first set of variables, e$1 , determines the optimal topology (solid or void) of the overall structure, while the second set, e$2 , selects the material at each solid element. This idea is demonstrated for a 3-phase scenario in Table 1.
This work further proposes extending Eq. (4) to a general material interpolation scheme that can handle any number of materials: E e ð e$1 , e$2 ,:::, e$M Þ ¼ where M is the number of distinct, non-void materials. When M ¼ 1 and M ¼ 2, this scheme simplifies to the established expressions Eqs. (2) and (4). Using the proposed Eq. (5), M design variables per element are required, but the number of elements in the FEA simulation does not change when more materials are introduced. Since the majority of the computational expense in topology optimization is typically attributed to FEA, the cost of additional design variables due to multiple materials is relatively small. Based on the stable modified SIMP scheme, the proposed formulation is easily integrated into problems with any number of materials.

Deterministic topology optimization
The density-based deterministic optimization problem to minimize compliance is ð e$1 e$2 ::: e$M v e Þ -J M £0, 0£ e$i £1, e ¼ 1,2,:::,P, i ¼ 1,2,:::,M , where ρ is the vector containing all design variables e$i , c is the compliance of the structure, v e is the volume of element e, u is the displacement vector, f is the external load vector, and P is the total number of elements. The stiffness in each element, E e , is calculated using the proposed interpolation scheme (Eq. (5)), and the global stiffness matrix is KðρÞ ¼ Since the topology optimization problem is solved with gradient-based methods, the derivative of the compliance cðρÞ with respect to the density variables, e%i , is derived using the adjoint variable method [24] as Table 1 Expansion of the proposed interpolation scheme for three phases (M = 2)

Density variables
Young's modulus where u e is the element displacement vector of element e.

Robust topology optimization
To capture uncertainty mathematically, the random variables ω m 2Θ M and ω f 2Θ F are introduced, where Θ M and Θ F are the random sample spaces [43] corresponding to uncertainties in the material properties and loads, respectively. Here M is the number of uncertain solid phases as defined previously, while F is the number of uncertain external loads. The robust topology optimization is then formulated as follows: where cðρ,ω m ,ω f Þ is the compliance under uncertainty, characterized by its mean and standard deviation, ðcÞ and ðcÞ. A weighted sum of these two statistical moments using a weight constant, β, results in a multi-objective optimization problem, which allows β to be increased in order to put more emphasis on minimizing the variation in compliance of the final design.
where each modulus E j ðω j Þ is a realization of its corresponding random variable, ω j 2ω m . In the case when there is no material uncertainty (i.e., there is only load uncertainty, or the problem is deterministic), the above simply reverts to Eq. (5) due to the lack of ω m .

Calculation of statistical moments
The material and load uncertainties are propagated using UDR, which, when combined with Gauss-type quadrature, efficiently reduces the moments of multivariate probability distributions into a weighted sum of univariate functions [41]. For the method to be accurate, all random variables should be chosen such that there are no strong interactions between them [44,45]. Furthermore, it is assumed that all of the random variables are mutually independent, and for notational simplicity, they are not separated into ω m and ω f . Instead, they are lumped into ω¼½ω m ,ω f 2Θ Z , where Z is the total number of random variables. For example, if both material and load uncertainties are considered, Z ¼M þ F. The statistical moments ðcÞ and 2 ðcÞ can be expressed as 2 ðcðρ,ωÞÞ ¼ ! Θ Z ½cðρ,ωÞ -ðcðρ,ωÞÞ 2 f ω ðωÞdω, (11) where f ω ðωÞ is the joint probability density function of the random variables.
With UDR, the compliance can be approximated by reducing it to a sum of univariate functions. In each function, all random variables are held equal to their mean except for one, ω j , such that cðρ,ωÞ ffi Note that the univariate functions cðρ,μj j ¼ω j Þ are independent since the random variables are independent, and therefore Eqs. (10) and (11) can be approximated as follows [46]: ðcðρ,μj j ¼ω j ÞÞ -ðZ -1Þcðρ,μÞ, (13) 2 ðcðρ,ωÞÞ ¼ X Z j¼1 2 ðcðρ,μj j ¼ω j ÞÞ: Then, utilizing one-dimensional Gauss-type quadrature to calculate the univariate moments in Eqs. (13) and (14), the final expressions are derived as where N is the number of quadrature nodes, l j$k is the kth node of the random variable ω j , and w j$k is the corresponding weight. Using these methods, the statistical moments of compliance can be estimated with Z$N þ1 deterministic finite element simulations, a value that is, as discussed earlier, not dependent on the number of elements in the design domain. In comparison to tensor product quadrature, which is an alternative univariate moment estimation method requiring N Z FEA evaluations, UDR is much more efficient [33].

Sensitivity analysis
For gradient-based optimization of the robust problem, the derivative of the weighted objective function with respect to each artificial density variable, e$i , is where e ¼ 1,2,:::,P, and i ¼ 1,2,:::,M . As in the previous section, when applying UDR to estimate the sensitivities of the mean and variance, the random variables are fixed at their mean except for ω j , which is sampled using Gauss-type quadrature. Thus, Similar to Eq. (7) but with the addition of random variables, the derivative of compliance under uncertainty is where ω q 2ω m . This derivative can be calculated after solving Kðρ,ω m Þu ¼ f ðω f Þ, which is deterministic due to UDR since ω m and ω f are either held at the mean or equal to the quadrature node according to Eqs. (18) and (19).

Numerical examples
The robust multi-material formulation was demonstrated with a 4-phase 3D cantilever beam under different combinations of material and load uncertainty, where each strut of the lattice structure is one FEA element with three design variables. In all examples, the initial structure ( Fig. 1(a)) has dimensions 200 mmÂ100 mmÂ20 mm and a total of 19502 elements. It is created by repeating the cubic unit cell in Fig. 1(b), which is 10 mmÂ10 mmÂ10 mm and consists of 20 elements with 2 mm diameters, and removing overlapping struts. The uniformity of the unit cell is not guaranteed in the optimized result, however, as individual strut elements can be removed if e$1 ¼ 0.
The FEA at each iteration is performed in Altair HyperWorks OptiStruct, which has the capability to generate and simulate lattice structures, including differently oriented struts, as part of their lattice optimization module. Each lattice strut is modeled as a 1D CBEAM element, and has PBEAML properties with TYPE=ROD. The material properties and their volume ratio upper limits are listed in Table 2, and the boundary conditions are shown in Fig. 1(a).
The optimization problems are solved via the method of moving asymptotes (MMA) [47], which is efficient despite the large number of design variables and volume constraints. The design variables e$i are initially set equal to their corresponding volume ratio limits, e.g., e$1 ¼J 1 . To simplify the calculations, the volume of each element, v e , was assumed to equal one. In these examples, the generalized multi-material interpolation scheme (Eq. (9)) does not require high penalization, with p =2 sufficing to drive away most intermediate (gray) densities. However, it is necessary to significantly increase the MMA parameter asydecr from the default value to 0.97 in order to decrease oscillations in the optimization history. The asyincr parameter is lowered to 1.03 for faster convergence with a slight sacrifice to the compliance, the initial movelimit is 1.2, and the rest of the parameters are kept at the default values recommended by Svanberg [48]. The optimization is stopped when the relative difference between the objective function values between two iterations is less than 10 -6 .
For UDR, all examples use N ¼8 Gauss quadrature points. Increasing the number of points beyond this does not noticeably improve the accuracies of the estimated statistical moments in comparison to Monte Carlo with 10000 samples; instead it escalates the computational cost unncessarily. In addition, the weight b in the robust objective function J 0 (Eq. (8)) is varied to generate Pareto optimal fronts of each example, showing the effects of the applied uncertainties. When b = 0, the objective function reverts to the deterministic one (Eq. (6)) and the optimal topology in this case is essentially the deterministic one. As b is increased, more emphasis is put on minimizing the standard deviation of compliance, i.e., the uncertainty. Three cases of uncertainty are examined in the following sections: uncertainty in the 1) materials only, 2) load only, and 3) both materials and load.

Materials uncertainty
In the first case, the Young's moduli of the four material phases are taken as normally distributed random variables. Thus, each modulus is prescribed a mean and standard deviation. The means of the soft, medium and hard materials are listed in Table 2; the standard deviations are set to 100 MPa each, which is equivalent to 20%, 10%, and 5% of the mean values, respectively. The load is deterministic with a magnitude of 1.5 N ( Fig. 1(a)).
Despite applying a high weight on the standard deviation of compliance, b = 10, the history of the objective function decreases stably and monotonically overall to 10.2% of the original value in 38 iterations (Fig. 2). In Fig. 3, the tradeoff between the the mean and standard deviation of compliance can be observed. When b is decreased, the mean decreases while the standard deviation increases, i.e., more optimal but less robust solutions are found. This relationship can also be seen in Fig. 4, where the robust solution has more of the medium (blue) and hard (yellow) materials, which were assigned lower relative standard deviations (10% and 5% of their means, respectively) than the soft material (20% of its mean).
To validate the statistical moments that were estimated  using UDR, the final topologies were evaluated with 10000 Monte Carlo samples and compared in Table 3. The study shows that the UDR approximations are accurate within 3% of the Monte Carlo results.

Load uncertainty
For the example under load uncertainty, two normally distributed random variables are considered: The magnitude and the angle of the external load applied (Fig. 5). The mean of the load maginitude is 1.5 N and the standard deviation is 0.225 N, or 15% of the mean, while the mean and standard deviation of the angle are 0°and 10°, respectively. As before, the objective history for b = 10 is fairly smooth and stable, and the final value is 5.67% of the initial compliance after 37 iterations (Fig. 6). Unlike the case with only material uncertainty, the robust designs are dramatically different, with long horizontal bars of the hard (yellow) material added in order to brace the cantilever against loads with angles that are not 0° (Figs. 8 and 9). Consequently, the mean and standard deviation of compliance for the robust solutions are significantly lower than the deterministic one ( Fig. 7(a) and Table 4). Still, the Pareto optimal front considering only robust solutions shows a similar tradeoff as before ( Fig. 7(b)), and once again the topology with higher b consists of stronger materials (Fig. 9). When load uncertainty is considered, however, the standard deviation estimated using UDR is not as accurate as Monte Carlo sampling, with differences higher than 10% (Table 4).

Materials and load uncertainties
The final case demonstrates the proposed method's ability   to capture multiple sources of uncertainty by simultaneously considering material and load uncertainties. This example combines the two previous ones: The same normal distributions for the three solid materials as in Section 5.1, and the same for the load magnitude and angle as in Section 5.2, are applied. Again, the optimization when b = 10 converges well, achieving an objective value that is 5.71% of the initial in 34 iterations (Fig. 10). The Pareto frontier (Fig. 11) and topologies (Figs. 12 and 13) are akin to those when only load uncertainty is considered in that the mean and standard deviation of compliance for the robust solutions also decrease due to the formation of the hard material (yellow) bars. With the additional uncertainty in the materials, however, the compliance increases marginally and the accuracy of UDR deteriorates slightly (Table 5).    This work proposes a new generalized material interpolation scheme for an arbitrary number of materials, as well as the employment of UDR and Gauss-type quadrature to quantify and propagate any combination of uncertainty in the materials and loads. The multi-material robust topology optimization method is demonstrated with mesoscale lattice structures, although the methods can be easily transferred to macroscale lattices, such as trusses, as well as continuum structures.
The numerical examples of 4-phase cantilever beam lattice structures reveal that the method is efficient for large, multi-material problems with multiple sources of uncertainty. The optimization converges stably in under 40 iterations despite nearly 20000 elements and up to five random parameters. Using UDR and Gauss-type quadrature significantly cuts down the computational expense of calculating the statistical moments of the objective function, requiring only eight FEA evaluations per    iteration rather than, for example, 10000 using Monte Carlo simulation. Although there is some loss of accuracy with UDR, the difference compared to the Monte Carlo approach is between 0.02% and 2.14% for the mean of compliance, and 1.13%-14.66% for the standard deviation, which is acceptable in exchange for the computational savings. Furthermore, by setting the objective function as a weighted sum of mean and standard deviation, multiple designs along the Pareto frontier can be generated and compared to observe and regulate the effect of uncertainty on the problem.
With the rise of highly functional, multi-material lattice structures built by AM, this paper offers a solution to efficiently and robustly design these complex structures. Possible extensions of this work include but are not limited to: Applying objectives more complicated than compliance; assigning targeted properties to achieve functionally graded structures; incorporating the mesoscale design within a multiscale framework that also optimizes the global geometry; and optimizing multi-material conformal lattice-skin structures.   Commons licence and your intended use is not permitted by statutory regulation or exceeds the permitted use, you will need to obtain permission directly from the copyright holder.