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

Abstract

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 macro-scale 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.

Key words： robust topology optimization; lattice structures; multi-material; material uncertainty; load uncertainty; univariate dimension reduction

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 [¹–³], acoustic manipulation [⁴,⁵] and energy absorption [⁶]. Mesoscale lattice structures have also been embedded into global structures to achieve multifunctional heat exchangers [⁷,⁸], and parts with high energy absorption [⁹–¹²] and low weight [¹³,¹⁴]. In particular, lattice structures are exceptionally popular in the biomedical field due to their porosity, which allows biocompability with organic tissue, high performance-to-weight ratio, and ease of customization [¹⁵–¹⁷]. By exploiting the advancements in multi-material AM, greater performance improvements [¹²], like graded stiffness [²] and a larger range of stiffness [¹⁸], are possible.

This work focuses on periodic, mesoscale, multi-material 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 [¹⁹] and multiphase phase-field [²⁰] 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 [²¹] and pseudo-sensitivity [²²] approaches exist, but are prohibitive for large problems. Zhang et al. [²³] 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. [²³], 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 Bendsoe and Sigmund [²⁴]. 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) [²⁵]. 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 [²⁶] 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. [²⁷] 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 [²⁸] 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 [²⁹], 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 [³⁰,³¹]. Park et al. [³²] 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 [³³–³⁵] and material [³⁶,³⁷] uncertainties exist for single material continuous problems. There is relatively more research for random loads, such as Refs. [³⁸–⁴⁰]. 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) [⁴¹] 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 [²⁴]

(1)E_e=E({\rho }_e)={\rho }_e^{p}E,

where {\rho }_e\in (\mathrm{0},\mathrm{1}] is the artificial density variable of element e such that it is greater than some small number {\rho }_{\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 [²⁵]:

(2)E_e=E_e({\rho }_e)={\rho }_e^p(E_{\mathrm{1}}-E_{\mathrm{0}})+E_{\mathrm{0}},

where E_{\mathrm{1}} is the elastic modulus of the solid material, E_{\mathrm{0}} is a small value that represents the void phase, and p is the penalty coefficient to discourage intermediate densities within [\mathrm{0},\mathrm{1}]. This improved formulation not only allows E_{\mathrm{0}} to be independent of p but also facilitates filtering schemes such as the minimum length scale filter.

Gibiansky and Sigmund [⁴²] showed the following 3-phase (two distinct solid materials and one void phase) scheme for composites that includes a second penalization parameter, q:

(3)E({\rho }_{e\cdot \mathrm{1}},{\rho }_{e\cdot \mathrm{2}})={\rho }_{e\cdot \mathrm{1}}^p[(\mathrm{1}-{\rho }_{e\cdot \mathrm{2}}^q)E_{\mathrm{1}}+{\rho }_{e\cdot \mathrm{2}}^q{E}_{\mathrm{2}}],

where E_{\mathrm{2}} is the modulus of the second solid material, and the design variables are {\rho }_{e\cdot \mathrm{1}},{\rho }_{e\cdot \mathrm{2}}\in [\mathrm{0},\mathrm{1}]. This model was tested by Stegmann and Lund in Ref. [²⁶], 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:

(4)\begin{array}{ll}{E}_e({\rho }_{e\cdot \mathrm{1}},{\rho }_{e\cdot \mathrm{2}})\hfill & ={\rho }_{e\cdot \mathrm{1}}^p[{\rho }_{e\cdot \mathrm{2}}^p{E}_{\mathrm{2}}+(\mathrm{1}-{\rho }_{e\cdot \mathrm{1}}^p)E_{\mathrm{1}}]+(\mathrm{1}-{\rho }_{e\cdot \mathrm{1}}^p)E_{\mathrm{0}}\hfill \\ \hfill & ={\rho }_{e\cdot \mathrm{1}}^p{\rho }_{e\cdot \mathrm{2}}^p(E_{\mathrm{2}}-E_{\mathrm{1}})+{\rho }_{e\cdot \mathrm{1}}^p(E_{\mathrm{1}}-E_{\mathrm{0}})+E_{\mathrm{0}}.\hfill \end{array}

From the second line, it is clear that the first set of variables, {\rho }_{e\cdot \mathrm{1}}, determines the optimal topology (solid or void) of the overall structure, while the second set, {\rho }_{e\cdot \mathrm{2}}, selects the material at each solid element. This idea is demonstrated for a 3-phase scenario in Table 1.

Tab.1 Expansion of the proposed interpolation scheme for three phases (M = 2)

Density variables	Young’s modulus
{\rho }_{e\cdot \mathrm{1}}=\mathrm{0}, {\rho }_{e\cdot \mathrm{2}}=\mathrm{0}	E_{\mathrm{0}}
{\rho }_{e\cdot \mathrm{1}}=\mathrm{1}, {\rho }_{e\cdot \mathrm{2}}=\mathrm{0}	E_1
{\rho }_{e\cdot \mathrm{1}}=\mathrm{1}, {\rho }_{e\cdot \mathrm{2}}=\mathrm{1}	E_2

This work further proposes extending Eq. (4) to a general material interpolation scheme that can handle any number of materials:

(5)E_e({\rho }_{e\cdot \mathrm{1}},{\rho }_{e\cdot \mathrm{2}},\dots ,{\rho }_{e\cdot M})={\displaystyle \sum _{j=\mathrm{1}}^M\left[(E_j-E_{j-\mathrm{1}}){\displaystyle \prod _{k=\mathrm{1}}^j{\rho }_{e\cdot k}^p}\right]}+E_{\mathrm{0}},

where M is the number of distinct, non-void materials. When M=\mathrm{1} and M=\mathrm{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:

(6)$$\begin{gathered} \underset{\rho }{\min \mathrm{i}\mathrm{m}\mathrm{i}\mathrm{z}\mathrm{e}}{J}_0=c\left(\rho \right)={\mathit{u}}^{\mathrm{T}}\mathit{K}\left(\rho \right)\mathit{u} \\ \mathrm{s}\mathrm{u}\mathrm{b}\mathrm{j}\mathrm{e}\mathrm{c}\mathrm{t}\mathrm{t}\mathrm{o}\mathit{K}\left(\rho \right)\mathit{u}=\mathit{f}, \\ {J}_1=\underset{e=1}{\overset{P}{\mathrm{\Sigma }}}\left({\rho }_{e\cdot 1}{v}_e\right)-{\overline{J}}_1\le 0, \\ {J}_2=\underset{e=1}{\overset{P}{\mathrm{\Sigma }}}\left({\rho }_{e\cdot 1}{\rho }_{e\cdot 2}{v}_e\right)-{\overline{J}}_2\le 0, \\ ⋮ \\ {J}_M=\underset{e=1}{\overset{P}{\mathrm{\Sigma }}}\left({\rho }_{e\cdot 1}{\rho }_{e\cdot 2}\mathrm{...}{\rho }_{e\cdot M}{v}_e\right)-{\overline{J}}_M\le 0, \\ 0\le {\rho }_{e\cdot i}\le 1,e=1,\mathrm{2...},P,i=1,2,\mathrm{...}M \end{gathered}$$

where \rho is the vector containing all design variables {\rho }_{e\cdot 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(\rho )={\displaystyle \sum _{e=\mathrm{1}}^P{E}_e{K}_e}, with the unit local stiffness matrix K_e. The volume ratios of the solid materials in each element are constrained by equations J_{\mathrm{1}},J_{\mathrm{2}},\dots ,J_M and the upper bounds of those ratios, {\bar{J}}_{\mathrm{1}},{\bar{J}}_{\mathrm{2}},\dots ,{\bar{J}}_M.

Since the topology optimization problem is solved with gradient-based methods, the derivative of the compliance c(\rho ) with respect to the density variables, {\rho }_{e\cdot i}^{}, is derived using the adjoint variable method [²⁴] as

(7)\frac{\partial c\left(\rho \right)}{{\partial \rho }_{e\cdot i}}=-u^{\mathrm{T}}\frac{\partial \mathit{K}\left(\rho \right)}{{\partial \rho }_{e\cdot i}}\mathit{u}=-p{\rho }_{e\cdot i}^{p-1}\underset{j=i}{\overset{M}{\mathrm{\Sigma }}}\left[(E_j-E_{j-1})\underset{k=1,k\ne i}{\overset{j}{\mathrm{\Pi }}}{\rho }_{e\cdot i}^{p-1}\right]{\mathit{u}}_e^{\mathrm{T}}{\mathit{K}}_e{\mathit{u}}_e,

where u_e is the element displacement vector of element e.

Robust topology optimization

To capture uncertainty mathematically, the random variables {\omega }_m\in {\Theta }^M and {\omega }_f\in {\Theta }^F are introduced, where {\Theta }^M and {\Theta }^F are the random sample spaces [⁴³] 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.

(8)$$\begin{gathered} \underset{\rho }{\min \mathrm{i}\mathrm{m}\mathrm{i}\mathrm{z}\mathrm{e}}{J}_0=\mu \left(c(\mathit{\rho },{\mathit{\omega }}_m,{\mathit{\omega }}_f)\right)+\beta \sigma \left(c(\mathit{\rho },{\mathit{\omega }}_m,{\mathit{\omega }}_f)\right) \\ \mathrm{s}\mathrm{u}\mathrm{b}\mathrm{j}\mathrm{e}\mathrm{c}\mathrm{t}\mathrm{t}\mathrm{o}\mathit{K}(\rho ,{\mathit{\omega }}_m)\mathit{u}=\mathit{f}\left({\mathit{\omega }}_f\right), \\ {J}_1=\underset{e=1}{\overset{P}{\mathrm{\Sigma }}}\left({\rho }_{e\cdot 1}{v}_e\right)-{\overline{J}}_1\le 0, \\ {J}_2=\underset{e=1}{\overset{P}{\mathrm{\Sigma }}}\left({\rho }_{e\cdot 1}{\rho }_{e\cdot 2}{v}_e\right)-{\overline{J}}_2\le 0, \\ ⋮ \\ {J}_M=\underset{e=1}{\overset{P}{\mathrm{\Sigma }}}\left({\rho }_{e\cdot 1}{\rho }_{e\cdot 2}\mathrm{...}{\rho }_{e\cdot M}{v}_e\right)-{\overline{J}}_M\le 0, \\ 0\le {\rho }_{e\cdot i}\le 1,e=1,\mathrm{2...},P,i=1,2,\mathrm{...}M \end{gathered}$$

where c(\rho ,{\omega }_m,{\omega }_f) is the compliance under uncertainty, characterized by its mean and standard deviation, \mu (c) and \sigma (c). A weighted sum of these two statistical moments using a weight constant, \beta, results in a multi-objective optimization problem, which allows \beta to be increased in order to put more emphasis on minimizing the variation in compliance of the final design.

Uncertainty quantification and propagation

Multi-material interpolation scheme under uncertainty

When material uncertainty exists, an extended version of the proposed multi-material scheme (Eq. (5)) can be used to calculate the Young’s modulus of each element, E_e:

(9)E_e({\rho }_{e\cdot \mathrm{1}},{\rho }_{e\cdot \mathrm{2}}^{},\dots ,{\rho }_{e\cdot M}^{})\hspace{0.17em}=\hspace{0.17em}{\displaystyle \sum _{j=\mathrm{1}}^M\left[(E_j({\omega }_j)-E_{j-1}({\omega }_{j-\mathrm{1}})){\displaystyle \prod _{k=\mathrm{1}}^j{\rho }_{e\cdot k}^p}\right]}\hspace{0.17em}+E_0,

where each modulus E_j({\omega }_j) is a realization of its corresponding random variable, {\omega }_j\in {\omega }_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 {\omega }_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 [⁴¹]. For the method to be accurate, all random variables should be chosen such that there are no strong interactions between them [⁴⁴,⁴⁵]. Furthermore, it is assumed that all of the random variables are mutually independent, and for notational simplicity, they are not separated into {\omega }_m and {\omega }_f. Instead, they are lumped into \omega =[{\omega }_m,{\omega }_f]\in {\Theta }^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 \mu (c) and {\sigma }^{\mathrm{2}}(c) can be expressed as

(10)\mu (c(\rho ,\omega ))={\displaystyle {\int }_{{\Theta }^Z}c(\rho ,\omega )}\hspace{0.17em}{f}_{\omega }(\omega )d\omega ,

(11){\sigma }^{\mathrm{2}}(c(\rho ,\omega ))={{\int }_{{\Theta }^Z}[c(\rho ,\omega )-\mu (c(\rho ,\omega ))]}^{\mathrm{2}}{f}_{\omega }(\omega )d\omega ,

where f_{\omega }(\omega ) 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, {\omega }_j, such that

(12)\begin{array}{l}c(\rho ,\omega )\hspace{0.17em}\cong \hspace{0.17em}{\displaystyle \sum _{j=\mathrm{1}}^{Z}c(\rho ,\mu ({\omega }_{\mathrm{1}}),\dots ,{\omega }_j,\dots ,\mu ({\omega }_Z))}-(Z-\mathrm{1})c(\rho ,\mu ({\omega }_{\mathrm{1}}),\dots ,\mu ({\omega }_Z))\\ \hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}={\displaystyle \sum _{j=\mathrm{1}}^{Z}c(\rho ,{\mu |}_{{\mu }_j={\omega }_j})}-(Z-\mathrm{1})c(\rho ,\mu ).\end{array}

Note that the univariate functions c(\rho ,{\mu |}_{{\mu }_j={\omega }_j}) are independent since the random variables are independent, and therefore Eqs. (10) and (11) can be approximated as follows [⁴⁶]:

(13)\mu (c(\rho ,\omega ))={\displaystyle \sum _{j=\mathrm{1}}^Z\mu (c(\rho ,{\mu |}_{{\mu }_j={\omega }_j}))}-(Z-\mathrm{1})c(\rho ,\mu ),

(14){\sigma }^{\mathrm{2}}(c(\rho ,\omega ))={\displaystyle \sum _{j=\mathrm{1}}^Z{\sigma }^{\mathrm{2}}(c(\rho ,{\mu |}_{{\mu }_j={\omega }_j}))}.

Then, utilizing one-dimensional Gauss-type quadrature to calculate the univariate moments in Eqs. (13) and (14), the final expressions are derived as

(15)\mu (c(\rho ,\omega ))\cong {\displaystyle \sum _{j=\mathrm{1}}^Z\left[{\displaystyle \sum _{k=\mathrm{1}}^N{w}_{j\cdot k}}c(\rho ,{\mu |}_{{\mu }_j=l_{j\cdot k}})\right]}-(Z-\mathrm{1})c(\rho ,\mu ),

(16){\sigma }^{\mathrm{2}}(c(\rho ,\omega ))\cong {\displaystyle \sum _{j=\mathrm{1}}^Z{\displaystyle \sum _{k=\mathrm{1}}^N{w}_{j\cdot k}{\left[c(\rho ,{\mu |}_{{\mu }_j=l_{j\cdot k}})-{\displaystyle \sum _{n=\mathrm{1}}^N{w}_{j\cdot n}\hspace{0.17em}c(\rho ,{\mu |}_{{\mu }_j=l_{j\cdot n}})}\right]}^{\mathrm{2}}}}

where N is the number of quadrature nodes, l_{j\cdot k} is the kth node of the random variable {\omega }_j, and w_{j\cdot k} is the corresponding weight.

Using these methods, the statistical moments of compliance can be estimated with Z\cdot N+\mathrm{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 [³³].

Sensitivity analysis

For gradient-based optimization of the robust problem, the derivative of the weighted objective function with respect to each artificial density variable, {\rho }_{e\cdot i}, is

(17){\displaystyle \frac{\partial J_{\mathrm{0}}}{\partial {\rho }_{e\cdot i}}}={\displaystyle \frac{\partial \mu (c(\rho ,\omega ))}{\partial {\rho }_{e\cdot i}}}+\beta {\displaystyle \frac{\partial \sigma (c(\rho ,\omega ))}{\partial {\rho }_{e\cdot i}}}\hspace{0.17em}={\displaystyle \frac{\partial \mu (c(\rho ,\omega ))}{\partial {\rho }_{e\cdot i}}}+{\displaystyle \frac{\beta }{\mathrm{2}\sigma (c(\rho ,\omega ))}}{\displaystyle \frac{\partial {\sigma }^{\mathrm{2}}(c(\rho ,\omega ))}{\partial {\rho }_{e\cdot i}}}

where e=\mathrm{1},\mathrm{2},\mathrm{...},P, and i=\mathrm{1},\mathrm{2},\mathrm{...},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 {\omega }_j, which is sampled using Gauss-type quadrature. Thus,

(18){\displaystyle \frac{\partial \mu (c(\rho ,\omega ))}{\partial {\rho }_{e\cdot i}}}\cong {\displaystyle \sum _{j=\mathrm{1}}^Z\left[{\displaystyle \sum _{k=\mathrm{1}}^N{w}_{j\cdot k}}{\displaystyle \frac{\partial c(\rho ,{\mu |}_{{\mu }_j=l_{j\cdot k}})}{\partial {\rho }_{e\cdot i}}}\right]}-(Z-\mathrm{1}){\displaystyle \frac{\partial c(\rho ,\mu )}{\partial {\rho }_{e\cdot i}}},

(19)\begin{array}{r}{\displaystyle \frac{\partial {\sigma }^{\mathrm{2}}(c(\rho ,\omega ))}{\partial {\rho }_{e\cdot i}}}\cong \mathrm{2}{\displaystyle \sum _{j=\mathrm{1}}^Z{\displaystyle \sum _{k=\mathrm{1}}^N{w}_{j\cdot k}}[c(\rho ,{\mu |}_{{\mu }_j=l_{j\cdot k}})-\mu (c(\rho ,{\mu |}_{{\mu }_j=l_{j\cdot k}}))]}\hspace{0.17em}\\ \hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\times \hspace{0.17em}\left[{\displaystyle \frac{\partial c(\rho ,{\mu |}_{{\mu }_j=l_{j\cdot k}})}{\partial {\rho }_{e\cdot i}}}-{\displaystyle \sum _{n=\mathrm{1}}^N{w}_{j\cdot n}}{\displaystyle \frac{\partial c(\rho ,{\mu |}_{{\mu }_j=l_{j\cdot n}})}{\partial {\rho }_{e\cdot i}}}\right].\end{array}

Similar to Eq. (7) but with the addition of random variables, the derivative of compliance under uncertainty is

(20)\frac{\partial c(\rho ,\omega )}{{\partial \rho }_{e\cdot i}}=-u^{\mathrm{T}}\frac{\partial \mathit{K}(\rho ,\omega )}{{\partial \rho }_{e\cdot i}}\mathit{u}=-p{\rho }_{e\cdot i}^{p-1}\underset{q=i}{\overset{M}{\mathrm{\Sigma }}}\left[(E_q\left({\omega }_q\right)-E_{q-1}\left({\omega }_{q-1}\right))\underset{k=1,k\ne i}{\overset{q}{\mathrm{\Pi }}}{\rho }_{e\cdot k}^p\right]{\mathit{u}}_e^{\mathrm{T}}{\mathit{K}}_e{\mathit{u}}_e,

where {\omega }_q\in {\omega }_m. This derivative can be calculated after solving K(\rho ,{\omega }_m)u=f({\omega }_f), which is deterministic due to UDR since {\omega }_m and {\omega }_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 {\rho }_{e\cdot \mathrm{1}}=\mathrm{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).

Fig.1 (a) Boundary conditions; (b) initial lattice unit cells

Full size|PPT slide

Tab.2 Material properties of the 4-phase examples

Material	Color	Mean Young’s modulus/MPa	Poisson’s ratio	Volume ratio constraint
Void	‒	1.0e‒8	0.3	‒
Soft	Green	500	0.3	10%
Medium	Blue	1000	0.3	10%
Hard	Yellow	2000	0.3	10%

The optimization problems are solved via the method of moving asymptotes (MMA) [⁴⁷], which is efficient despite the large number of design variables and volume constraints. The design variables {\rho }_{e\cdot i} are initially set equal to their corresponding volume ratio limits, e.g., {\rho }_{e\cdot \mathrm{1}}={\bar{J}}_{\mathrm{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 [⁴⁸]. The optimization is stopped when the relative difference between the objective function values between two iterations is less than {\mathrm{10}}^{-\mathrm{6}}.

For UDR, all examples use N=\mathrm{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_{\mathrm{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).

Fig.2 Optimization history of the objective function under material uncertainty when b = 10

Full size|PPT slide

Fig.3 Pareto curve under material uncertainty with different objective function weights b

Full size|PPT slide

Fig.4 Optimal solutions under material uncertainty when (a) b = 0 (deterministic) and (b) b = 10 (robust)

Full size|PPT slide

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.

Tab.3 Mean and standard deviation of the compliance (N∙mm) results under material uncertainty using UDR with 8 quadrature nodes and Monte Carlo with 10000 samples

Solution type	Statistical moment	UDR/quadrature	Monte Carlo	Difference/%
Deterministic (b = 0)	Mean	0.43036	0.43027	0.02
	Std. Dev.	0.02063	0.02087	1.13
Robust (b = 10)	Mean	0.46244	0.46262	0.04
	Std. Dev.	0.02018	0.02072	2.61

Std. Dev.: Standard deviation

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).

Fig.5 Boundary conditions with the random external load, which varies in both magnitude (f) and angle (q)

Full size|PPT slide

Fig.6 Optimization history for the objective function under load uncertainty when b = 10

Full size|PPT slide

Fig.7 Pareto curve for different objective function weights under load uncertainty (a) including b = 0 (deterministic), and (b) showing robust solutions only (b = 1 to b = 10)

Full size|PPT slide

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), 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).

Fig.8 Optimal solution under load uncertainty when b = 0 (deterministic)

Full size|PPT slide

Fig.9 Optimal solutions under load uncertainty when (a) b = 1 and (b) b = 10

Full size|PPT slide

Tab.4 Compliance (N∙mm) results under load uncertainty using UDR with 8 quadrature nodes and Monte Carlo with 10000 samples

Solution type	Statistical moment	UDR/quadrature	Monte Carlo	Difference/%
Deterministic (b = 0)	Mean	1.02900	1.04167	1.22
	Std. Dev.	0.82306	0.92341	10.87
Robust (b = 1)	Mean	0.74000	0.74949	1.27
	Std. Dev.	0.42108	0.49094	14.23
Robust (b = 10)	Mean	0.76359	0.77067	0.92
	Std. Dev.	0.38892	0.44047	11.70

Std. Dev.: Standard deviation

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).

Fig.10 Optimization history for the objective function under both material and load uncertainty when b = 10

Full size|PPT slide

Fig.11 Pareto curve for different objective function weights under both material and load uncertainty (a) including b = 0 (deterministic), and (b) showing robust solutions only (b = 1 to b = 10)

Full size|PPT slide

Fig.12 Optimal solution under both material and load uncertainty when b = 0 (deterministic)

Full size|PPT slide

Fig.13 Optimal solutions under both material and load uncertainty when (a) b = 1 and (b) b = 10

Full size|PPT slide

Tab.5 Compliance (N∙mm) results under material and load uncertainties using UDR with 8 quadrature nodes and Monte Carlo with 10000 samples

Solution type	Statistical moment	UDR/quadrature	Monte Carlo	Difference/%
Deterministic (b = 0)	Mean	1.03223	1.03533	0.30
	Std. Dev.	0.82331	0.92008	10.52
Robust (b = 1)	Mean	0.74380	0.76006	2.14
	Std. Dev.	0.42198	0.49447	14.66
Robust (b = 10)	Mean	0.77624	0.78349	0.93
	Std. Dev.	0.39135	0.45541	14.07

Std. Dev.: Standard deviation

Conclusions

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.

Acknowledgements

Ms. Yu-Chin Chan would like to acknowledge support from the Digital Manufacturing and Design Innovation Institute (DMDII) through award number 15-07-07. Additionally, this material is based upon work supported by the National Science Foundation Graduate Research Fellowship Program under Grant No. DGE-1842165. Any opinions, findings, and conclusions or recommendations expressed in this material are those of the author(s) and do not necessarily reflect the views of the National Science Foundation.

Open Access

This article is licensed under a Creative Commons Attribution 4.0 International License, which permits use, sharing, adaptation, distribution and reproduction in any medium or format, as long as you give appropriate credit to the original author(s) and the source, provide a link to the Creative Commons licence, and indicate if changes were made.

The images or other third party material in this article are included in the article’s Creative Commons licence, unless indicated otherwise in a credit line to the material. If material is not included in the article’s Creative 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.

To view a copy of this licence, visit http://creativecommons.org/licenses/by/4.0/.

References

1 Chen Y, Li T, Scarpa F, Lattice metamaterials with mechanically tunable Poisson’s ratio for vibration control. Physical Review Applied, 2017, 7(2): 024012

2 Chen D, Zheng X. Multi-material additive manufacturing of metamaterials with giant, tailorable negative Poisson’s ratios. Scientific Reports, 2018, 8(1): 9139

3 Yuan S, Shen F, Bai J, 3D soft auxetic lattice structures fabricated by selective laser sintering: TPU powder evaluation and process optimization. Materials & Design, 2017, 120: 317–327

4 Yu K, Fang N X, Huang G, Magnetoactive acoustic metamaterials. Advanced Materials, 2018, 30(21): 1706348

5 Bilal O R, Ballagi D, Daraio C. Architected lattices for simultaneous broadband attenuation of airborne sound and mechanical vibrations in all directions. 2018, arXiv:1809.01252v1

6 Shan S, Kang S H, Raney J R, Multistable architected materials for trapping elastic strain energy. Advanced Materials, 2015, 27(29): 4296–4301

7 Maloney K J, Fink K D, Schaedler T A, Multifunctional heat exchangers derived from three-dimensional micro-lattice structures. International Journal of Heat and Mass Transfer, 2012, 55(9–10): 2486–2493

8 Son K N, Weibel J A, Kumaresan V, Design of multifunctional lattice-frame materials for compact heat exchangers. International Journal of Heat and Mass Transfer, 2017, 115(Part A): 619–629

9 Brennan-Craddock J, Brackett D, Wildman R, The design of impact absorbing structures for additive manufacture. Journal of Physics: Conference Series, 2012, 382: 012042

10 Ozdemir Z, Tyas A, Goodall R, Energy absorption in lattice structures in dynamics: Nonlinear FE simulations. International Journal of Impact Engineering, 2017, 102: 1–15

11 Habib F N, Iovenitti P, Masood S H, Fabrication of polymeric lattice structures for optimum energy absorption using multi jet fusion technology. Materials & Design, 2018, 155: 86–98

12 Stankovic T, Mueller J, Egan P, A generalized optimality criteria method for optimization of additively manufactured multimaterial lattice structures. Journal of Mechanical Design, 2015, 137(11): 111405

13 Moon S K, Tan Y E, Hwang J, Application of 3D printing technology for designing light-weight unmanned aerial vehicle wing structures. International Journal of Precision Engineering and Manufacturing-Green Technology, 2014, 1(3): 223–228

14 Terriault P, Brailovski V. Modeling and simulation of large, conformal, porosity-graded and lightweight lattice structures made by additive manufacturing. Finite Elements in Analysis and Design, 2018, 138: 1–11

15 Murr L E, Gaytan S M, Medina F, Next-generation biomedical implants using additive manufacturing of complex, cellular and functional mesh arrays. Philosophical Transactions of the Royal Society A: Mathematical, Physical and Engineering Sciences, 2010, 368(1917): 1999–2032

16 Parthasarathy J, Starly B, Raman S. A design for the additive manufacture of functionally graded porous structures with tailored mechanical properties for biomedical applications. Journal of Manufacturing Processes, 2011, 13(2): 160–170

17 Zhang X Z, Leary M, Tang H P, Selective electron beam manufactured Ti-6Al-4V lattice structures for orthopedic implant applications: Current status and outstanding challenges. Current Opinion in Solid State and Materials Science, 2018, 22(3): 75–99

18 Saxena K K, Das R, Calius E P. 3D printable multimaterial cellular auxetics with tunable stiffness. 2017, arXiv:1707.04486

19 Wang M Y, Wang X. “Color” level sets: A multi-phase method for structural topology optimization with multiple materials. Computer Methods in Applied Mechanics and Engineering, 2004, 193(6–8): 469–496

20 Zhou S, Wang M Y. Multimaterial structural topology optimization with a generalized Cahn-Hilliard model of multiphase transition. Structural and Multidisciplinary Optimization, 2007, 38: 89–111

21 Hiller J D, Lipson H. Multi material topological optimization of structures and mechanisms. In: Proceedings of the 11th Annual Conference on Genetic and Evolutionary Computation. Montreal: ACM, 2009, 1521–1528

22 Ramani A. Multi-material topology optimization with strength constraints. Structural and Multidisciplinary Optimization, 2011, 43(5): 597–615

23 Zhang X S, Paulino G H, Ramos A S. Multi-material topology optimization with multiple volume constraints: A general approach applied to ground structures with material nonlinearity. Structural and Multidisciplinary Optimization, 2018, 57(1): 161–182

24 Bendsøe M P, Sigmund O. Topology Optimization: Theory, Methods, and Applications. Berlin: Springer, 2004

25 Sigmund O. Morphology-based black and white filters for topology optimization. Structural and Multidisciplinary Optimization, 2007, 33(4–5): 401–424

26 Stegmann J, Lund E. Discrete material optimization of general composite shell structures. International Journal for Numerical Methods in Engineering, 2005, 62(14): 2009–2027

27 Gaynor A, Meisel N A, Williams C B, Multiple-material topology optimization of compliant mechanisms created via polyjet 3D printing. Journal of Manufacturing Science and Engineering, 2014, 136(6): 061015

28 Yin L, Ananthasuresh G K. Topology optimization of compliant mechanisms with multiple materials using a peak function material interpolation scheme. Structural and Multidisciplinary Optimization, 2001, 23(1): 49–62

29 Zuo W, Saitou K. Multi-material topology optimization using ordered SIMP interpolation. Structural and Multidisciplinary Optimization, 2017, 55(2): 477–491

30 Sing S L, Wiria F E, Yeong W Y. Selective laser melting of lattice structures: A statistical approach to manufacturability and mechanical behavior. Robotics and Computer-Integrated Manufacturing, 2018, 49: 170–180

31 Ahmadi S M, Hedayati R, Ashok Kumar Jain R K, Effects of laser processing parameters on the mechanical properties, topology, and microstructure of additively manufactured porous metallic biomaterials: A vector-based approach. Materials & Design, 2017, 134: 234–243

32 Park S I, Rosen D W, Choi S, Effective mechanical properties of lattice material fabricated by material extrusion additive manufacturing. Additive Manufacturing, 2014, 1–4: 12–23

33 Chen S, Chen W. A new level-set based approach to shape and topology optimization under geometric uncertainty. Structural and Multidisciplinary Optimization, 2011, 44(1): 1–18

34 Jansen M, Lombaert G, Diehl M, et al. Robust topology optimization accounting for misplacement of material. Structural and Multidisciplinary Optimization, 2013, 47(3): 317–333

35 Schevenels M, Lazarov B S, Sigmund O. Robust topology optimization accounting for spatially varying manufacturing errors. Computer Methods in Applied Mechanics and Engineering, 2011, 200(49–52): 3613–3627

36 Jalalpour M, Tootkaboni M. An efficient approach to reliability-based topology optimization for continua under material uncertainty. Structural and Multidisciplinary Optimization, 2016, 53(4): 759–772

37 Changizi N, Jalalpour M. Robust topology optimization of frame structures under geometric or material properties uncertainties. Structural and Multidisciplinary Optimization, 2017, 56(4): 791–807

38 Alvarez F, Carrasco M. Minimization of the expected compliance as an alternative approach to multiload truss optimization. Structural and Multidisciplinary Optimization, 2005, 29(6): 470–476

39 Dunning P D, Kim H A, Mullineux G. Introducing loading uncertainty in topology optimization. AIAA Journal, 2011, 49(4): 760–768

40 Deng J, Chen W. Concurrent topology optimization of multiscale structures with multiple porous materials under random field loading uncertainty. Structural and Multidisciplinary Optimization, 2017, 56(1): 1–19

41 Rahman S, Xu H. A univariate dimension-reduction method for multi-dimensional integration in stochastic mechanics. Probabilistic Engineering Mechanics, 2004, 19(4): 393–408

42 Gibiansky L V, Sigmund O. Multiphase composites with extremal bulk modulus. Journal of the Mechanics and Physics of Solids, 2000, 48(3): 461–498

43 Yin X, Lee S, Chen W, et al. Efficient random field uncertainty propagation in design using multiscale analysis. Journal of Mechanical Design, 2009, 131(2): 021006

44 Lee S H, Chen W. A comparative study of uncertainty propagation methods for black-box-type problems. Structural and Multidisciplinary Optimization, 2009, 37(3): 239–253

45 Lee S H, Chen W, Kwak B M. Robust design with arbitrary distributions using Gauss-type quadrature formula. Structural and Multidisciplinary Optimization, 2009, 39(3): 227–243

46 Zhao Y G, Ono T. Moment methods for structural reliability. Structural Safety, 2001, 23(1): 47–75

47 Svanberg K. The method of moving asymptotes—A new method for structural optimization. International Journal for Numerical Methods in Engineering, 1987, 24(2): 359–373

48 Svanberg K. MMA and GCMMA—Two Methods for Nonlinear Optimization. Technical Report. 2007