A new multi-scale topology optimization framework for optimal combinations of macro-layouts and local gradings of TPMS-based lattice structures

Abstract

In this work, optimal combinations of macro-layouts and local gradings of triply periodic minimal surface (TPMS)-based lattice structures are obtained by using multi-scale topology optimization. A new innovative framework is proposed by using two density variables in each finite element. The first variable is the local relative lattice density and it sets the effective orthotropic elastic properties of the element, which in turn are obtained by using numerical homogenization of representative volume elements of the particular TPMS-based lattice structure of interest. The second variable is a standard topology optimization macro density variable, which defines if the element should be treated as a void or contain the graded lattice structure by letting this variable be governed by the rational approximation of material properties (RAMP) model. By using such density variables for all elements, the compliance is minimized by separately constraining the volume of lattice structure and the volume of macro-layout by using two independent constraints. For benchmarks in 3 D, it is demonstrated that the stiffness is increased significantly by including local grading of the lattice structure compared to using a constant lattice density. It is also demonstrated how ultra-lightweight designs can be generated using the multi-scale formulation, and how the optimal multi-scale solutions easily can be realized to printable stl-files by using implicit based geometry modeling. Finally, the new multi-scale topology optimization framework is utilized to generate an optimal design combination of macro-layout and local grading of frame-based Gyroid structure for the established GE-bracket benchmark.

Keywords:

• Topology optimization
• functional grading
• triply periodic minimal surfaces
• numerical homogenization
• additive manufacturing constraints

1. Introduction

The use of lattice structures in design for additive manufacturing has quickly emerged as a popular and efficient design alternative for creating innovative multifunctional lightweight solutions (Seharing, Azman, and Abdullah 2020; Plocher and Panesar 2018). In particular, the family of triply periodic minimal surfaces (TPMS) studied in detail by Schoen (1970) for generating lattice structures seems extra promising. The Gyroid surface is probably the most well-known of Schoen’s 12 TPMS-based surfaces, but others exist such as e.g., the G-prime and the Schwarz-D surface. An extensive study of the mechanical properties of several TPMS-based lattices can be found in the work by Maskery et al. (2018). Typically, in the design process, bulk parts are replaced by these kind of lattice structures, which then are globally graded with finite element analysis in a trial-and-error manner in order to satisfy design requirements and additive manufacturing constraints.

A better mechanical design approach would be to simultaneously replace the bulk part with a new optimal macro-layout containing locally optimal graded lattice structures. In this paper, such an approach is presented by proposing a new multi-scale topology optimization framework for simultaneously generating optimal combinations of macro-layouts and local gradings of periodic lattice structures. Additive manufacturing constraints on the thickness of the lattices are also included in the formulation by constraining the relative density of the lattice with lower and upper limits. Recently, Jansen and Pierard (2020) presented similar ideas using a hybrid approach in order to simultaneously optimize the shape of the macroscale structure and the functionally graded lattices inside this design domain. Their approach was demonstrated by studying the well-known GE bracket. This benchmark was recently studied by Strömberg (2021a) by grading shell-based TPMS lattice structures using a single variable approach. The macro-layout in that work was obtained by forcing low densities below a prescribed limit toward zeros using a Sigmoid filter. In this work, the macro layout is obtained as one part of the optimal solution using a standard macro topology optimization density variable governed by RAMP (Stolpe and Svanberg 2001), which in turn is filled with graded lattice structure given as a second part of the optimal solution using an additional local lattice density variable governed by effective elastic properties of the lattice.

A review on current and future trends in topology optimization, covering grading of lattices using topology optimization, can be found in (Liu et al. 2018). Strategies for functionally graded lattice structures by using topology optimization were discussed by Panesar et al. (2018). A recent review on design, optimization and additive manufacturing of cellular structures was given by Nazir et al. (2019), considering both strut- and TPMS-based lattices structures. A most recent review on multi-scale topology optimization is presented by Wu, Sigmund, and Groen (2021). Most of the work on optimal grading of lattices structures using topology optimization is for strut-based lattices or on honeycomb lattices in 2 D. A mathematical rigourous work on topology optimization of periodic microstructures by the homogenization method is presented by Allaire, Geoffroy-Donders, and Pantz (2019). An early work on grading honeycomb lattices in 2 D without macro-layout was suggested and validated by Zhang et al. (2015). This work was extended in (Zhang, Liu, and To 2017) to include anisotropic elasticity. Anisotropic optimization of lattice structures was recently considered by Li et al. (2020). Daynes et al. (2017) graded struts by minimizing compliance and utilized the principle stress directions to define the layouts of the struts. Jin, Li, and Zhang (2018) performed sizing of struts by mapping topology optimization results onto strut-based lattices. A two-step sequential level-set approach for macro-layout and grading of struts in 2 D was presented in Wang et al. (2018a). Natural frequency optimization of strut lattices was performed in Cheng et al. (2018a) and of honeycomb structures in Wang et al. (2018b). Maximizing natural frequencies of periodic strut structures with Kriging-assisted topology optimization was proposed by Zhang, Gao, and Xiao (2020). Lattice structure topology optimization for heat conduction design was investigated by Cheng et al. (2018b). Stress constraints for lattice grading were considered in (Cheng, Bai, and To 2019). Most recently, thermal distortions from additive manufacturing processes were included by using the inherent strain method (Takezawa et al. 2020).

An early work on two-scale topology optimization is presented by Rodrigues, Guedes, and Bendsø (2002), where a hierarchical computational procedure for generating an optimal combination of macro-layout with local element material properties was suggested. A few years later, Liu, Yan, and Cheng (2008) proposed the PAMP-model for two-scale topology optimization of cellular material in order to design ultralight structures. Vicente et al. (2016) performed concurrent topology optimization for minimizing frequency response of two-level hierarchical structures. Furthermore, concurrent two-scale topological design of multiple unit cells using a combined level set and density model approach was presented by Wang and Kang (2019). Most recently, Liu, Kang, and Luo (2020) proposed a method to obtain connectable micro-structures when using two-scale concurrent topology optimization of lattice structures.

Studies on optimal graded TPMS-based lattice structures by topology optimization appear less frequently. Recently, an interesting approach for optimal graded TPMS-based lattices was proposed by Li, Dai, et al. (2019). Strömberg (2021a) graded shell-based Gyroid, G-prime and Schwarz-D lattice structures using a single density variable with a Sigmoid filter for forcing densities below a lower limit toward zero. In this work, optimal combinations of macro-layouts and local gradings of the frame-based Gyroid structure are generated using the new proposed innovative framework.

The outline of the paper is the following: in the next section the governing equations of the global stiffness matrix using the two density variables as well as different volume measures of the design domain, macro layout and lattice structure are given. In Section 3, the new innovative multi-scale topology optimization formulation is presented together with sensitivity analysis and adopted filters. In Section 4, the effective elastic properties of the frame-based Gyroid structure as a function of the relative lattice density are derived. In Section 5, the proposed framework is demonstrated for a 3 D L-shaped design domain subjected to both bending and torsion. The compliance is minimized for different constraints on the volume fractions of the macro layout and the lattice structure, while keeping the volume of bulk material fixed. In such manner, different optimal combinations of macro layouts and local gradings of TPMS-based lattice structures are generated. This is done for a single force and a weighted load case of three orthogonal forces. In addition, an ultra-lightweight design of the L-shaped benchmark is generated for the single load case. The design domain of the established GE bracket is also studied for a weighted load case. At the end of this section, by using implicit-based geometry, the optimal combinations of macro layouts and local gradings of lattices are transformed to printable stl-files. Finally, some concluding remarks are given.

2. Governing equations

A non-homogenous linear orthotropic elastic body with fixed displacements and prescribed external forces, which is discretized with linear finite elements is considered, see Figure 1. Two density variables, γ_e and ρ_e, are introduce for each element e; γ_e is the local relative density of lattice in the element and ρ_e is a standard macro density variable defining if the element should be treated as a void of filled with lattice structure. Such density variables for all elements are utilized to set up a global stiffness matrix for a design domain with local graded lattice structures according to the procedure presented below.

Figure 1. Each finite element e in the body is equipped with two density variables, ρ_e and γ_e. ρ_e is governed by RAMP and determines if an element is a void or is filled with lattice material. γ_e gives the relative density of lattice material for each element e and it is governed by effective elastic properties derived by numerical homogenization.

The effective elastic properties of each element e is obtained by numerical homogenization of representative volume elements (RVEs) of the TPMS-based lattice structure of interest such that the effective elasticity matrix in Voigt notation is given by (1) $${\mathit{C}}_e={\mathit{C}}_e({\gamma }_e)=\left[\begin{array}{cccccc}{f}_{11}{c}_{11}& f_{21}{c}_{21}& f_{31}{c}_{31}& 0& 0& 0\\ f_{21}{c}_{21}& f_{22}{c}_{11}& f_{32}{c}_{31}& 0& 0& 0\\ f_{31}{c}_{31}& f_{32}{c}_{31}& f_{33}{c}_{33}& 0& 0& 0\\ 0& 0& 0& f_{44}{c}_{44}& 0& 0\\ 0& 0& 0& 0& f_{55}{c}_{44}& 0\\ 0& 0& 0& 0& 0& f_{66}{c}_{66}\end{array}\right],$$(1) where $f_{\mathit{\text{ij}}}=f_{\mathit{\text{ij}}}({\gamma }_e)$ are material interpolation (regression) laws as function of the relative lattice density $0\le l_b\le {\gamma }_e\le u_b\le 1$ of element e, where l_b and u_b are prescribed lower and upper limits on γ_e, respectively.

Typically, the lower limit l_b represents the thinnest member that can be printed and u_b is utilized for establishing practical removal of powder. Furthermore, c_ij in (Eq. 1(1) $${\mathit{C}}_e={\mathit{C}}_e({\gamma }_e)=\left[\begin{array}{cccccc}{f}_{11}{c}_{11}& f_{21}{c}_{21}& f_{31}{c}_{31}& 0& 0& 0\\ f_{21}{c}_{21}& f_{22}{c}_{11}& f_{32}{c}_{31}& 0& 0& 0\\ f_{31}{c}_{31}& f_{32}{c}_{31}& f_{33}{c}_{33}& 0& 0& 0\\ 0& 0& 0& f_{44}{c}_{44}& 0& 0\\ 0& 0& 0& 0& f_{55}{c}_{44}& 0\\ 0& 0& 0& 0& 0& f_{66}{c}_{66}\end{array}\right],$$(1) ) are the elastic properties of the bulk material of the lattice structure, which are assumed to be governed by transversely isotropic elasticity, i.e., ${\mathit{C}}_e{(1)}^{-1}=$ (2) $$\left[\begin{array}{cccccc}1/E& -\nu /E& -{\nu }_{13}/E& 0& 0& 0\\ -\nu /E& 1/E& -{\nu }_{13}/E& 0& 0& 0\\ -{\nu }_{13}/E& -{\nu }_{13}/E& 1/E_{33}& 0& 0& 0\\ 0& 0& 0& 1/G_{23}& 0& 0\\ 0& 0& 0& 0& 1/G_{23}& 0\\ 0& 0& 0& 0& 0& 2(1+\nu )/E\end{array}\right],$$(2) where E and E₃₃ are Young’s moduli, ν and ν₁₃ are Poisson’s ratios and G₂₃ is the out-of-plane shear moduli. Thus, the 3-direction is aligned with the building direction of the additive manufacturing process of the lattice and it is assumed that isotropic properties are obtained in the build plane.

Thus, each element e has unique linear orthotropic elastic properties, defined by ${\mathit{C}}_e$ in (Eq. 1(1) $${\mathit{C}}_e={\mathit{C}}_e({\gamma }_e)=\left[\begin{array}{cccccc}{f}_{11}{c}_{11}& f_{21}{c}_{21}& f_{31}{c}_{31}& 0& 0& 0\\ f_{21}{c}_{21}& f_{22}{c}_{11}& f_{32}{c}_{31}& 0& 0& 0\\ f_{31}{c}_{31}& f_{32}{c}_{31}& f_{33}{c}_{33}& 0& 0& 0\\ 0& 0& 0& f_{44}{c}_{44}& 0& 0\\ 0& 0& 0& 0& f_{55}{c}_{44}& 0\\ 0& 0& 0& 0& 0& f_{66}{c}_{66}\end{array}\right],$$(1) ), which in turn depends on the relative lattice density γ_e. The stresses and strains, i.e., (3) $$\mathit{\sigma }=\left\{\begin{array}{c}{\sigma }_{11}\\ {\sigma }_{22}\\ {\sigma }_{33}\\ {\sigma }_{23}\\ {\sigma }_{13}\\ {\sigma }_{12}\end{array}\right\}\text{\hspace{0.5em}and\hspace{0.5em}}\mathit{ϵ}=\left\{\begin{array}{c}{ϵ}_{11}\\ {ϵ}_{22}\\ {ϵ}_{33}\\ 2{ϵ}_{23}\\ 2{ϵ}_{13}\\ 2{ϵ}_{12}\end{array}\right\},$$(3) are, of course, governed by $\mathit{\sigma }={\mathit{C}}_e({\gamma }_e)\mathit{ϵ},$ and the corresponding local effective finite element stiffness matrix for element e is denoted by (4) $${\mathit{k}}_e={\mathit{k}}_e({\gamma }_e).$$(4)

In the global assembly procedure, the RAMP model is utilized in order to define if the finite element e should be considered to be a void or filled with lattice structure. Thus, the global stiffness matrix is generated by (5) $$\mathit{K}=\mathit{K}(\mathit{\rho },\mathit{\gamma })=\underset{e}{\cap }\frac{{\rho }_e}{1+n(1-{\rho }_e)}{\mathit{k}}_e({\gamma }_e),$$(5) where, $\mathit{\rho }=\left\{{\rho }_e\right\},\hspace{0.17em}\mathit{\gamma }=\left\{{\gamma }_e\right\},\hspace{0.17em}\cap$ represents an assembly operator, n is the RAMP factor and $0\le {\rho }_e\le 1$ is the relative macro density of lattice structure for element e. Here, ${\rho }_e=1$ implies a completely filled element with lattice structure and ${\rho }_e=0$ represents a void. In the numerical implementation, the zero is replaced with a small number ϵ in order to avoid singular stiffness matrices. Instead of using the RAMP model, the established SIMP model can be adopted in (Eq. 5(5) $$\mathit{K}=\mathit{K}(\mathit{\rho },\mathit{\gamma })=\underset{e}{\cap }\frac{{\rho }_e}{1+n(1-{\rho }_e)}{\mathit{k}}_e({\gamma }_e),$$(5) ), which then reads (6) $$\mathit{K}=\mathit{K}(\mathit{\rho },\mathit{\gamma })=\underset{e}{\cap }{\rho }_e^n{\mathit{k}}_e({\gamma }_e).$$(6)

The total volume of bulk material generated in the assembly procedure is given by (7) $$V^{\text{bulk}}=V^{\text{bulk}}(\mathit{\rho },\mathit{\gamma })=\sum _e{\rho }_e{\gamma }_e{V}_e,$$(7) where V_e represents the total volume of each element e when ${\rho }_e={\gamma }_e=1.$ By using (Eq. 7(7) $$V^{\text{bulk}}=V^{\text{bulk}}(\mathit{\rho },\mathit{\gamma })=\sum _e{\rho }_e{\gamma }_e{V}_e,$$(7) ), the following different volume measures can be identified (8) $$\begin{array}{c}{V}^{\text{design}}=V^{\text{bulk}}(\mathbf{1},\mathbf{1}),\hfill \\ V^{\text{macro}}=V^{\text{macro}}(\mathit{\rho })=V^{\text{bulk}}(\mathit{\rho },\mathbf{1}),\hfill \\ V^{\text{lattice}}=V^{\text{lattice}}(\mathit{\gamma })=V^{\text{bulk}}(\mathbf{1},\mathit{\gamma }),\hfill \end{array}$$(8) where $V^{\text{design}}$ is the total volume of the design domain, $V^{\text{macro}}$ is the total volume of the macro layout filled with lattice structure, $V^{\text{lattice}}$ is the total volume of lattice in the design domain and $\mathbf{1}$ is a vector of ones. The volume of lattice in the macro layout of lattice structure $V^{\text{macro}}$ is given by (9) $$V_{\text{macro}}^{\text{lattice}}=V^{\text{lattice}}-(V^{\text{design}}-V^{\text{macro}})l_b,$$(9) which will converge toward the volume of bulk material $V^{\text{bulk}},$ because ${\rho }_e=ϵ$ in the void region. $V^{\text{macro}}$ and $V^{\text{lattice}}$ are constrained separately in the multi-scale topology optimization formulation presented in (Eq. 14(14) $$\begin{array}{c}\underset{(\mathit{\rho },\mathit{\gamma },\mathit{d})}{\mathrm{\text{min}}}c={\mathit{F}}^T\mathit{d}\hfill \\ \mathrm{s}.\mathrm{t}.\{\begin{array}{c}\mathit{K}(\mathit{\rho },\mathit{\gamma })\mathit{d}=\mathit{F},\hfill \\ {\mathrm{V}}^{\text{macro}}(\mathit{\rho })\le {\widehat{V}}^{\text{macro}},\hfill \\ {\mathrm{V}}^{\text{lattice}}(\mathit{\gamma })\le {\widehat{V}}^{\text{lattice}},\hfill \\ ϵ\mathbf{1}\le \mathit{\rho }\le \mathbf{1},\hfill \\ {\mathrm{l}}_{\mathrm{b}}\mathbf{1}\le \mathit{\gamma }\le u_b\mathbf{1},\hfill \end{array}\hfill \end{array}$$(14) ) in the next section. The different volume measures introduced above in (Eq. 8(8) $$\begin{array}{c}{V}^{\text{design}}=V^{\text{bulk}}(\mathbf{1},\mathbf{1}),\hfill \\ V^{\text{macro}}=V^{\text{macro}}(\mathit{\rho })=V^{\text{bulk}}(\mathit{\rho },\mathbf{1}),\hfill \\ V^{\text{lattice}}=V^{\text{lattice}}(\mathit{\gamma })=V^{\text{bulk}}(\mathbf{1},\mathit{\gamma }),\hfill \end{array}$$(8) ) and (Eq. 9(9) $$V_{\text{macro}}^{\text{lattice}}=V^{\text{lattice}}-(V^{\text{design}}-V^{\text{macro}})l_b,$$(9) ) using $V^{\text{bulk}}$ in (Eq. 7(7) $$V^{\text{bulk}}=V^{\text{bulk}}(\mathit{\rho },\mathit{\gamma })=\sum _e{\rho }_e{\gamma }_e{V}_e,$$(7) ) are illustrated in Figure 2 for a standard benchmark solved by (Eq. 14(14) $$\begin{array}{c}\underset{(\mathit{\rho },\mathit{\gamma },\mathit{d})}{\mathrm{\text{min}}}c={\mathit{F}}^T\mathit{d}\hfill \\ \mathrm{s}.\mathrm{t}.\{\begin{array}{c}\mathit{K}(\mathit{\rho },\mathit{\gamma })\mathit{d}=\mathit{F},\hfill \\ {\mathrm{V}}^{\text{macro}}(\mathit{\rho })\le {\widehat{V}}^{\text{macro}},\hfill \\ {\mathrm{V}}^{\text{lattice}}(\mathit{\gamma })\le {\widehat{V}}^{\text{lattice}},\hfill \\ ϵ\mathbf{1}\le \mathit{\rho }\le \mathbf{1},\hfill \\ {\mathrm{l}}_{\mathrm{b}}\mathbf{1}\le \mathit{\gamma }\le u_b\mathbf{1},\hfill \end{array}\hfill \end{array}$$(14) ) using $V^{\text{bulk}}=0.25V^{\text{design}},\hspace{0.17em}{V}^{\text{macro}}=0.5V^{\text{design}},\hspace{0.17em}{l}_b=0.2$ and $u_b=0.6.$ Thus, $V_{\text{macro}}^{\text{lattice}}=V^{\text{bulk}}$ and (Eq. 9(9) $$V_{\text{macro}}^{\text{lattice}}=V^{\text{lattice}}-(V^{\text{design}}-V^{\text{macro}})l_b,$$(9) ) yields (10) $$V^{\text{lattice}}=(0.25+(1-0.5)0.2)V^{\text{design}}=0.35V^{\text{design}}.$$(10)

Figure 2. Illustration of the different volume measures defined in (Eqs. 8(8) $$\begin{array}{c}{V}^{\text{design}}=V^{\text{bulk}}(\mathbf{1},\mathbf{1}),\hfill \\ V^{\text{macro}}=V^{\text{macro}}(\mathit{\rho })=V^{\text{bulk}}(\mathit{\rho },\mathbf{1}),\hfill \\ V^{\text{lattice}}=V^{\text{lattice}}(\mathit{\gamma })=V^{\text{bulk}}(\mathbf{1},\mathit{\gamma }),\hfill \end{array}$$(8) and 9(9) $$V_{\text{macro}}^{\text{lattice}}=V^{\text{lattice}}-(V^{\text{design}}-V^{\text{macro}})l_b,$$(9) ) by using the volume of bulk material in (Eq. 7(7) $$V^{\text{bulk}}=V^{\text{bulk}}(\mathit{\rho },\mathit{\gamma })=\sum _e{\rho }_e{\gamma }_e{V}_e,$$(7) ).

Consequently, (11) $$\begin{array}{c}{V}^{\text{macro}}(\mathit{\rho })\le 0.5V^{\text{design}},\hfill \\ V^{\text{lattice}}(\mathit{\gamma })\le 0.35V^{\text{design}}\hfill \end{array}$$(11) are used in (Eq. 14(14) $$\begin{array}{c}\underset{(\mathit{\rho },\mathit{\gamma },\mathit{d})}{\mathrm{\text{min}}}c={\mathit{F}}^T\mathit{d}\hfill \\ \mathrm{s}.\mathrm{t}.\{\begin{array}{c}\mathit{K}(\mathit{\rho },\mathit{\gamma })\mathit{d}=\mathit{F},\hfill \\ {\mathrm{V}}^{\text{macro}}(\mathit{\rho })\le {\widehat{V}}^{\text{macro}},\hfill \\ {\mathrm{V}}^{\text{lattice}}(\mathit{\gamma })\le {\widehat{V}}^{\text{lattice}},\hfill \\ ϵ\mathbf{1}\le \mathit{\rho }\le \mathbf{1},\hfill \\ {\mathrm{l}}_{\mathrm{b}}\mathbf{1}\le \mathit{\gamma }\le u_b\mathbf{1},\hfill \end{array}\hfill \end{array}$$(14) ).

From a engineering point of view, the following volume fractions are of interest, see also (Eq. 29(29) $$\begin{array}{c}{\widehat{\xi }}^{\text{bulk}}=\frac{V^{\text{bulk}}(\{{\rho }_e^{\text{heav}}>{\rho }_t\},\mathit{\gamma })}{V^{\text{bulk}}(\mathbf{1},\mathbf{1})},\hfill \\ {\widehat{\xi }}^{\text{macro}}=\frac{V^{\text{bulk}}(\{{\rho }_e^{\text{heav}}>{\rho }_t\},\mathbf{1})}{V^{\text{bulk}}(\mathbf{1},\mathbf{1})},\hfill \\ {\widehat{\xi }}_{\text{macro}}^{\text{lattice}}=\frac{V^{\text{bulk}}(\{{\rho }_e^{\text{heav}}>{\rho }_t\},\mathit{\gamma })}{V^{\text{bulk}}(\{{\rho }_e^{\text{heav}}>{\rho }_t\},\mathbf{1})}\mathrm{ }\text{for}\mathrm{ }{V}^{\text{bulk}}(\{{\rho }_e^{\text{heav}}>{\rho }_t\},\mathbf{1})>0.\hfill \end{array}$$(29) ): (12) $$\begin{array}{c}{\xi }^{\text{bulk}}=\frac{V^{\text{bulk}}}{V^{\text{design}}},\hfill \\ {\xi }^{\text{macro}}=\frac{V^{\text{macro}}}{V^{\text{design}}},\hfill \\ {\xi }_{\text{macro}}^{\text{lattice}}=\frac{V_{\text{macro}}^{\text{lattice}}}{V^{\text{macro}}},\hfill \end{array}$$(12) where ${\xi }^{\text{bulk}}$ is the fraction of bulk material in the design domain, ${\xi }^{\text{macro}}$ is the fraction of macro layout and ${\xi }_{\text{macro}}^{\text{lattice}}$ is the fraction of lattice in the macro layout. By starting from $V_{\text{macro}}^{\text{lattice}}=V^{\text{bulk}},$ inserting (Eq. 12(12) $$\begin{array}{c}{\xi }^{\text{bulk}}=\frac{V^{\text{bulk}}}{V^{\text{design}}},\hfill \\ {\xi }^{\text{macro}}=\frac{V^{\text{macro}}}{V^{\text{design}}},\hfill \\ {\xi }_{\text{macro}}^{\text{lattice}}=\frac{V_{\text{macro}}^{\text{lattice}}}{V^{\text{macro}}},\hfill \end{array}$$(12) ), one can derive (13) $${\xi }_{\text{macro}}^{\text{lattice}}{\xi }^{\text{macro}}={\xi }^{\text{bulk}}.$$(13)

This relationship is utilized in the formulation of the benchmarks in the numerical study in Section 5, such that an upper limit on the bulk material can be given relative to the upper limits on $V^{\text{macro}}$ and $V^{\text{lattice}}$ used in (Eq. 14(14) $$\begin{array}{c}\underset{(\mathit{\rho },\mathit{\gamma },\mathit{d})}{\mathrm{\text{min}}}c={\mathit{F}}^T\mathit{d}\hfill \\ \mathrm{s}.\mathrm{t}.\{\begin{array}{c}\mathit{K}(\mathit{\rho },\mathit{\gamma })\mathit{d}=\mathit{F},\hfill \\ {\mathrm{V}}^{\text{macro}}(\mathit{\rho })\le {\widehat{V}}^{\text{macro}},\hfill \\ {\mathrm{V}}^{\text{lattice}}(\mathit{\gamma })\le {\widehat{V}}^{\text{lattice}},\hfill \\ ϵ\mathbf{1}\le \mathit{\rho }\le \mathbf{1},\hfill \\ {\mathrm{l}}_{\mathrm{b}}\mathbf{1}\le \mathit{\gamma }\le u_b\mathbf{1},\hfill \end{array}\hfill \end{array}$$(14) ).

3. Topology optimization, sensitivity analysis and filters

In this section, a new innovative multi-scale topology optimization formulation is proposed for minimizing the compliance of the design domain with one constraint on the macro layout volume $V^{\text{macro}}$ and an additional constraint on the lattice volume $V^{\text{lattice}}.$

The new multi-scale formulation reads (14) $$\begin{array}{c}\underset{(\mathit{\rho },\mathit{\gamma },\mathit{d})}{\mathrm{\text{min}}}c={\mathit{F}}^T\mathit{d}\hfill \\ \mathrm{s}.\mathrm{t}.\{\begin{array}{c}\mathit{K}(\mathit{\rho },\mathit{\gamma })\mathit{d}=\mathit{F},\hfill \\ {\mathrm{V}}^{\text{macro}}(\mathit{\rho })\le {\widehat{V}}^{\text{macro}},\hfill \\ {\mathrm{V}}^{\text{lattice}}(\mathit{\gamma })\le {\widehat{V}}^{\text{lattice}},\hfill \\ ϵ\mathbf{1}\le \mathit{\rho }\le \mathbf{1},\hfill \\ {\mathrm{l}}_{\mathrm{b}}\mathbf{1}\le \mathit{\gamma }\le u_b\mathbf{1},\hfill \end{array}\hfill \end{array}$$(14) where $\mathit{d}$ is the displacement vector, $\mathit{F}$ contains the external forces, c is the compliance, $\mathit{\text{Kd}}=\mathit{F}$ is the static equilibrium equation, ${\widehat{V}}^{\text{macro}}$ and ${\widehat{V}}^{\text{lattice}}$ are the upper limits on the macro layout volume and the lattice volume, respectively. The standard compliance problem is recovered by setting u_b=1 and ${\widehat{V}}^{\text{lattice}}=V^{\text{design}},$ and the local grading formulation for the whole design domain using effective elasticity properties, as presented in Strömberg (2021a), is obtained by setting ${\widehat{V}}^{\text{macro}}=V^{\text{design}}.$

By dividing the volume constraints in (Eq. 14(14) $$\begin{array}{c}\underset{(\mathit{\rho },\mathit{\gamma },\mathit{d})}{\mathrm{\text{min}}}c={\mathit{F}}^T\mathit{d}\hfill \\ \mathrm{s}.\mathrm{t}.\{\begin{array}{c}\mathit{K}(\mathit{\rho },\mathit{\gamma })\mathit{d}=\mathit{F},\hfill \\ {\mathrm{V}}^{\text{macro}}(\mathit{\rho })\le {\widehat{V}}^{\text{macro}},\hfill \\ {\mathrm{V}}^{\text{lattice}}(\mathit{\gamma })\le {\widehat{V}}^{\text{lattice}},\hfill \\ ϵ\mathbf{1}\le \mathit{\rho }\le \mathbf{1},\hfill \\ {\mathrm{l}}_{\mathrm{b}}\mathbf{1}\le \mathit{\gamma }\le u_b\mathbf{1},\hfill \end{array}\hfill \end{array}$$(14) ) by $V^{\text{design}}$ and utilizing (Eqs. 9(9) $$V_{\text{macro}}^{\text{lattice}}=V^{\text{lattice}}-(V^{\text{design}}-V^{\text{macro}})l_b,$$(9) , 12–14) can be rephrased to (15) $$\begin{array}{c}\underset{(\mathit{\rho },\mathit{\gamma },\mathit{d})}{\mathrm{\text{min}}}c={\mathit{F}}^T\mathit{d}\hfill \\ \mathrm{s}.\mathrm{t}.\{\begin{array}{c}\mathit{K}(\mathit{\rho },\mathit{\gamma })\mathit{d}=\mathit{F},\hfill \\ {\mathrm{V}}^{\text{macro}}(\mathit{\rho }){/\mathrm{V}}^{\text{design}}\le {\xi }^{\text{macro}},\hfill \\ {\mathrm{V}}^{\text{lattice}}(\mathit{\gamma }){/\mathrm{V}}^{\text{design}}\le {\xi }_{\text{macro}}^{\text{lattice}}{\xi }^{\text{macro}}+(1-{\xi }^{\text{macro}}){\mathrm{l}}_{\mathrm{b}},\hfill \\ ϵ\mathbf{1}\le \mathit{\rho }\le \mathbf{1},\hfill \\ {\mathrm{l}}_{\mathrm{b}}\mathbf{1}\le \mathit{\gamma }\le u_b\mathbf{1}.\hfill \end{array}\hfill \end{array}$$(15)

The new multi-scale topology optimization problem in (Eq. 14(14) $$\begin{array}{c}\underset{(\mathit{\rho },\mathit{\gamma },\mathit{d})}{\mathrm{\text{min}}}c={\mathit{F}}^T\mathit{d}\hfill \\ \mathrm{s}.\mathrm{t}.\{\begin{array}{c}\mathit{K}(\mathit{\rho },\mathit{\gamma })\mathit{d}=\mathit{F},\hfill \\ {\mathrm{V}}^{\text{macro}}(\mathit{\rho })\le {\widehat{V}}^{\text{macro}},\hfill \\ {\mathrm{V}}^{\text{lattice}}(\mathit{\gamma })\le {\widehat{V}}^{\text{lattice}},\hfill \\ ϵ\mathbf{1}\le \mathit{\rho }\le \mathbf{1},\hfill \\ {\mathrm{l}}_{\mathrm{b}}\mathbf{1}\le \mathit{\gamma }\le u_b\mathbf{1},\hfill \end{array}\hfill \end{array}$$(14) ) or (Eq. 15(15) $$\begin{array}{c}\underset{(\mathit{\rho },\mathit{\gamma },\mathit{d})}{\mathrm{\text{min}}}c={\mathit{F}}^T\mathit{d}\hfill \\ \mathrm{s}.\mathrm{t}.\{\begin{array}{c}\mathit{K}(\mathit{\rho },\mathit{\gamma })\mathit{d}=\mathit{F},\hfill \\ {\mathrm{V}}^{\text{macro}}(\mathit{\rho }){/\mathrm{V}}^{\text{design}}\le {\xi }^{\text{macro}},\hfill \\ {\mathrm{V}}^{\text{lattice}}(\mathit{\gamma }){/\mathrm{V}}^{\text{design}}\le {\xi }_{\text{macro}}^{\text{lattice}}{\xi }^{\text{macro}}+(1-{\xi }^{\text{macro}}){\mathrm{l}}_{\mathrm{b}},\hfill \\ ϵ\mathbf{1}\le \mathit{\rho }\le \mathbf{1},\hfill \\ {\mathrm{l}}_{\mathrm{b}}\mathbf{1}\le \mathit{\gamma }\le u_b\mathbf{1}.\hfill \end{array}\hfill \end{array}$$(15) ) is solved using sequential linear programming (SLP) using the sensitivities and filters presented below. The sensitivities of the compliance c, using the adjoint method, are given by (16) $$\frac{\partial c}{\partial {\rho }_e}=-{\mathit{d}}^T\frac{\partial \mathit{K}}{\partial {\rho }_e}\mathit{d},$$(16) (17) $$\frac{\partial c}{\partial {\gamma }_e}=-{\mathit{d}}^T\frac{\partial \mathit{K}}{\partial {\gamma }_e}\mathit{d},$$(17) where (18) $$\frac{\partial \mathit{K}}{\partial {\rho }_e}=\frac{n+1}{{(1+n(1-{\rho }_e))}^2}{\mathit{k}}_e({\gamma }_e)$$(18) and (19) $$\frac{\partial \mathit{K}}{\partial {\gamma }_e}=\frac{{\rho }_e}{1+n(1-{\rho }_e)}\frac{\partial {\mathit{k}}_e}{\partial {\gamma }_e},$$(19) where $\partial f_{\mathit{\text{ij}}}/\partial {\gamma }_e$ are needed in (Eq. 19(19) $$\frac{\partial \mathit{K}}{\partial {\gamma }_e}=\frac{{\rho }_e}{1+n(1-{\rho }_e)}\frac{\partial {\mathit{k}}_e}{\partial {\gamma }_e},$$(19) ).

The sensitivities of the volume $V^{\text{macro}}$ and $V^{\text{lattice}}$ are, of course, (20) $$\frac{\partial V^{\text{macro}}}{\partial {\rho }_e}=\frac{\partial V^{\text{lattice}}}{\partial {\gamma }_e}=V_e.$$(20)

The sensitivities are treated by a density filter for both ρ_e and γ_e (Bourdin 2001), and, in addition, ρ_e is passing a smooth Heaviside function (Wang, Lazarov, and Sigmund 2011). Thus, ρ_e is treated by (21) $${\rho }_e^{\text{filt}}=\frac{\sum _{g=1}^{n_{\text{el}}}{\delta }_g{V}_g{\rho }_g}{\sum _{f=1}^{n_{\text{el}}}{\delta }_f{V}_f},$$(21) where $n_{\text{el}}$ is the number of finite elements, (22) $${\delta }_f={\delta }_f(e)={\left(r_{\text{min}}-\text{dist}(e,f)\right)}_{+},$$(22)

$\text{dist}(e,f)$ denotes the distance between the center of element e and f, and $r_{\text{min}}$ is the filter radius, which is set to 1.5 times the characteristic length of the finite elements in the numerical examples. (Eq. 21(21) $${\rho }_e^{\text{filt}}=\frac{\sum _{g=1}^{n_{\text{el}}}{\delta }_g{V}_g{\rho }_g}{\sum _{f=1}^{n_{\text{el}}}{\delta }_f{V}_f},$$(21) ) is also applied on γ_e. Furthermore, (23) $${\rho }_e^{\text{heav}}=\frac{\mathrm{\text{tanh}}(\beta \eta )+\mathrm{\text{tanh}}(\beta ({\rho }_e^{\text{filt}}-\eta ))}{\mathrm{\text{tanh}}(\beta \eta )+\mathrm{\text{tanh}}(\beta (1-\eta ))},$$(23) where η defines the threshold and β sets the slope of the smooth Heaviside filter. In the numerical examples, $\eta =0.55,$ β is ramped from 1 to 20, and the filter is activated after 50 SLP iterations. The Heaviside filter is depicted in Figure 3 for this setting of η and different values on β. The adoption of the smooth Heaviside filter implies a very thin layer of intermediate values between the void region with ${\rho }_e^{\text{heav}}=ϵ$ and the lattice filled region with ${\rho }_e^{\text{heav}}=1,$ see e.g., the thin layer of intermediate values in $V^{\text{macro}}$ presented in Figure 2, see also the histogram of ${\rho }_e^{\text{heav}}$ presented in Figure 8. In the post-processing made in this paper, intermediate values ${\rho }_e^{\text{heav}}>0.6$ are interpreted as a part of the region with TPMS-based lattice structure and ${\rho }_e^{\text{heav}}\le 0.6$ are treated as voids. This could of course be modified by using a threshold value closer to 1, but the post-processed results will look very similar because almost all intermediate values are efficiently removed by the Heaviside filter. The latter is also demonstrated in Figure 8 by studying the measure of discreteness defined in (Eq. 30(30) $$M_{\text{nd}}=\sum _{e=1}^{n_{\text{el}}}\frac{4{\rho }_e^{\text{heav}}(1-{\rho }_e^{\text{heav}})}{n_{\text{el}}}\times 100\%.$$(30) ).

Figure 3. The smooth Heaviside filter in (Eq. 23(23) $${\rho }_e^{\text{heav}}=\frac{\mathrm{\text{tanh}}(\beta \eta )+\mathrm{\text{tanh}}(\beta ({\rho }_e^{\text{filt}}-\eta ))}{\mathrm{\text{tanh}}(\beta \eta )+\mathrm{\text{tanh}}(\beta (1-\eta ))},$$(23) ) for $\eta =0.55$ and different values on β.

Thus, in problem (Eq. 14(14) $$\begin{array}{c}\underset{(\mathit{\rho },\mathit{\gamma },\mathit{d})}{\mathrm{\text{min}}}c={\mathit{F}}^T\mathit{d}\hfill \\ \mathrm{s}.\mathrm{t}.\{\begin{array}{c}\mathit{K}(\mathit{\rho },\mathit{\gamma })\mathit{d}=\mathit{F},\hfill \\ {\mathrm{V}}^{\text{macro}}(\mathit{\rho })\le {\widehat{V}}^{\text{macro}},\hfill \\ {\mathrm{V}}^{\text{lattice}}(\mathit{\gamma })\le {\widehat{V}}^{\text{lattice}},\hfill \\ ϵ\mathbf{1}\le \mathit{\rho }\le \mathbf{1},\hfill \\ {\mathrm{l}}_{\mathrm{b}}\mathbf{1}\le \mathit{\gamma }\le u_b\mathbf{1},\hfill \end{array}\hfill \end{array}$$(14) ), the stiffness and volumes are given as functions of ${\rho }_e^{\text{heav}}={\rho }_e^{\text{heav}}({\rho }_e^{\text{filt}}(\mathit{\rho }))$ and ${\gamma }_e^{\text{filt}}={\gamma }_e^{\text{filt}}(\mathit{\gamma }),$ i.e., $\mathit{K}({\mathit{\rho }}^{\text{heav}},{\mathit{\gamma }}^{\text{filt}}),\hspace{0.17em}{V}^{\text{macro}}({\mathit{\rho }}^{\text{heav}})$ and $V^{\text{lattice}}({\mathit{\gamma }}^{\text{filt}}).$ The sensitivities in (Eqs. 16–20) are then calculated using the chain rule, e.g., (24) $$\begin{array}{c}\frac{\partial c}{\partial {\rho }_e}=\frac{\partial c}{\partial {\rho }_m^{\text{heav}}}\frac{\partial {\rho }_m^{\text{heav}}}{\partial {\rho }_m^{\text{filt}}}\frac{\partial {\rho }_m^{\text{filt}}}{\partial {\rho }_e},\hfill \\ \frac{\partial c}{\partial {\gamma }_e}=\frac{\partial c}{\partial {\gamma }_m^{\text{filt}}}\frac{\partial {\gamma }_m^{\text{filt}}}{\partial {\gamma }_e}.\hfill \end{array}$$(24)

4. The Gyroid lattice structure

The proposed multi-scale topology optimization framework is demonstrated in the next section by using the frame-based Gyroid lattice structure, which is formulated as an implicit-based surface, i.e., (25) $$f^{\text{lattice}}=\hspace{0.17em}\text{sin}\hspace{0.17em}(\omega x)\hspace{0.17em}\text{cos}\hspace{0.17em}(\omega y)+\hspace{0.17em}\text{sin}\hspace{0.17em}(\omega y)\hspace{0.17em}\text{cos}\hspace{0.17em}(\omega z)+\hspace{0.17em}\text{sin}\hspace{0.17em}(\omega z)\hspace{0.17em}\text{cos}\hspace{0.17em}(\omega x)+\kappa ,$$(25) where ω is controlling the size of the period and $\kappa =\kappa (\gamma )$ is the scalar level controlling the relative density of the lattice. A cube of 1x1x1 periods of the frame-based Gyroid lattice ($f^{\text{lattice}}\ge 0$) for $\gamma =0.2$ is depicted in grey in Figure 4. The corresponding shell-based Gyroid lattice structure was considered in Strömberg (2021a).

Figure 4. The material interpolation laws f_ij for the frame-based Gyroid structure obtained by numerical homogenization of RVEs for prescribed normal strains and simple shear strains. The markers represent the results from the finite element analysis of the RVEs.

The effective elastic properties of the lattice are obtained by using linear finite element analysis (FEA) using Abaqus/Standard of RVEs subjected to 3 prescribed normal strains, e.g., ${\mathit{ϵ}}^1=\{1,0,0,0,0,0\},$ and 3 prescribed simple shear strains, e.g., ${\mathit{ϵ}}^4=\{0,0,0,1,0,0\},$ in orthogonal directions. The corresponding stresses are then given by ${\mathit{\sigma }}^i={\mathit{C}}_e{\mathit{ϵ}}^i,$ e.g., ${\mathit{\sigma }}^1=\{f_{11}{c}_{11},f_{21}{c}_{21},f_{31}{c}_{31},0,0,0\}$ and ${\mathit{\sigma }}^4=\{0,0,0,f_{44}{c}_{44},0,0\}.$ The RVE used for the prescribed normal strains is a cube with 3x3x3 periods of the Gyroid cell, and the simple shear strains are treated using a rectangular cuboid with 9x3x3 periods of the Gyroid cell in order to reduce the influence of the free boundaries, see Figure 4. Approximately, 2 million linear tetrahedral elements are used for the 3x3x3 cube and 6 million elements for the 9x3x3 cuboid. This numerical homogenization approach is simple and straight-forward to implement (see e.g., Steven 1997 and Li, Dai, et al. 2019) and the computational time is only few minutes for each RVE on a PC workstation. An alternative treatment that would speed up the computational time would be to use a single periodic cell as presented in grey in Figure 4 with periodic boundary conditions (see e.g., Omairey, Dunning, and Sriramula 2019).

From the numerical homogenization, it is clear that the lattice exhibits cubic symmetric behavior, although the geometry has no cubic symmetry, see this by rotating the Gyroid cell depicted in grey in Figure 4 by 90 degrees. The cubic behavior of the Gyroid structure appears clearly when Young’s modulus is plotted in different directions in a polar plot, an example of such polar plot of Young’s modulus is found in the lower right corner of Figure 4 for γ_e=0.2.

As a result from the procedure above, see also Strömberg (2021a), the effective elastic matrix for the frame-based Gyroid structure is given by (26) $${\mathit{C}}_e=\left[\begin{array}{cccccc}{f}_{11}{c}_{11}& f_{21}{c}_{21}& f_{21}{c}_{21}& 0& 0& 0\\ f_{21}{c}_{21}& f_{11}{c}_{11}& f_{21}{c}_{21}& 0& 0& 0\\ f_{21}{c}_{21}& f_{21}{c}_{21}& f_{11}{c}_{11}& 0& 0& 0\\ 0& 0& 0& f_{44}{c}_{44}& 0& 0\\ 0& 0& 0& 0& f_{44}{c}_{44}& 0\\ 0& 0& 0& 0& 0& f_{44}{c}_{44}\end{array}\right],$$(26) where the material interpolation laws $f_{\mathit{\text{ij}}}=f_{\mathit{\text{ij}}}({\gamma }_e)$ are presented in Figure 4. These functions are given as the following convex combination of RAMP and SIMP basis functions: (27) $$\begin{array}{c}{f}_{\mathit{\text{ij}}}=\frac{{\alpha }_1^{\mathit{\text{ij}}}{\gamma }_e}{1+(1-{\gamma }_e)}+\frac{{\alpha }_2^{\mathit{\text{ij}}}{\gamma }_e}{1+2(1-{\gamma }_e)}+\frac{{\alpha }_3^{\mathit{\text{ij}}}{\gamma }_e}{1+3(1-{\gamma }_e)}+\dots \\ \frac{{\alpha }_4^{\mathit{\text{ij}}}{\gamma }_e}{1+5(1-{\gamma }_e)}+{\alpha }_1^{\mathit{\text{ij}}}{\gamma }_e+{\alpha }_6^{\mathit{\text{ij}}}{\gamma }_e^2+{\alpha }_7^{\mathit{\text{ij}}}{\gamma }_e^3,\end{array}$$(27) where ${\alpha }_k^{\mathit{\text{ij}}}\ge 0,\hspace{0.17em}\sum _k{\alpha }_k^{\mathit{\text{ij}}}=1$ are regression coefficients presented in Table 1. Concerning regression models with convex combination constraints, see e.g., Strömberg (2021b).

Table 1. The regression coefficients for the material interpolation laws f_ij in (Eq. 27(27) $$\begin{array}{c}{f}_{\mathit{\text{ij}}}=\frac{{\alpha }_1^{\mathit{\text{ij}}}{\gamma }_e}{1+(1-{\gamma }_e)}+\frac{{\alpha }_2^{\mathit{\text{ij}}}{\gamma }_e}{1+2(1-{\gamma }_e)}+\frac{{\alpha }_3^{\mathit{\text{ij}}}{\gamma }_e}{1+3(1-{\gamma }_e)}+\dots \\ \frac{{\alpha }_4^{\mathit{\text{ij}}}{\gamma }_e}{1+5(1-{\gamma }_e)}+{\alpha }_1^{\mathit{\text{ij}}}{\gamma }_e+{\alpha }_6^{\mathit{\text{ij}}}{\gamma }_e^2+{\alpha }_7^{\mathit{\text{ij}}}{\gamma }_e^3,\end{array}$$(27) ) for the frame-based Gyroid structure.

fij	${\alpha }_1^{\mathit{\text{ij}}}$	${\alpha }_2^{\mathit{\text{ij}}}$	${\alpha }_3^{\mathit{\text{ij}}}$	${\alpha }_4^{\mathit{\text{ij}}}$	${\alpha }_5^{\mathit{\text{ij}}}$	${\alpha }_6^{\mathit{\text{ij}}}$	${\alpha }_7^{\mathit{\text{ij}}}$
f11	0	0.4380	0.0055	0	0	0.1325	0.4240
f21	0	0	0.1380	0.6653	0	0	0.1967
f44	0.2267	0.0002	0.0001	0	0.0059	0.7671	0

5. Numerical examples

The proposed multi-scale topology optimization approach is implemented and demonstrated in this section for the frame-based Gyroid lattice structure presented in the previous section. The implementation is done in our in-house toolbox, TopoBox, using Matlab and Fortran. The pre- and post-processing of models and results are done in Abaqus/CAE.

A three-dimensional L-shaped benchmark similar to the one presented by Borrvall and Petersson (2001) is considered for different settings using the implementation. The design domain is presented in Figure 5, which is discretized using 157946 linear tetrahedral elements. The four corners of the back-side is fixed and the lower opposite corner is subjected to a single vertical force F_z or a weighted load case using F_x, F_y and F_z with equal magnitudes and equal weights. Young’s modulus and Poisson’s ratio for the isotropic bulk material is 2.1E5 N/mm² and 0.3, respectively. Notice that this benchmark includes both bending and torsion.

Figure 5. The L-shaped design domain subjected to a single load case F_z or a weighted load case (F_x, F_y, F_z). The compliance for the single load case is compared using no grading or optimal grading.

The L-shaped benchmark is first studied by only optimizing the macro-layout for almost constant relative density of the lattice structure and constant volume of bulk material. This is done by the formulation in (Eq. 14(14) $$\begin{array}{c}\underset{(\mathit{\rho },\mathit{\gamma },\mathit{d})}{\mathrm{\text{min}}}c={\mathit{F}}^T\mathit{d}\hfill \\ \mathrm{s}.\mathrm{t}.\{\begin{array}{c}\mathit{K}(\mathit{\rho },\mathit{\gamma })\mathit{d}=\mathit{F},\hfill \\ {\mathrm{V}}^{\text{macro}}(\mathit{\rho })\le {\widehat{V}}^{\text{macro}},\hfill \\ {\mathrm{V}}^{\text{lattice}}(\mathit{\gamma })\le {\widehat{V}}^{\text{lattice}},\hfill \\ ϵ\mathbf{1}\le \mathit{\rho }\le \mathbf{1},\hfill \\ {\mathrm{l}}_{\mathrm{b}}\mathbf{1}\le \mathit{\gamma }\le u_b\mathbf{1},\hfill \end{array}\hfill \end{array}$$(14) ) by setting ${\xi }^{\text{bulk}}=0.2,\hspace{0.17em}{\xi }_{\text{macro}}^{\text{lattice}}\in$[0.2:0.1:0.7], and $l_b={\xi }_{\text{macro}}^{\text{lattice}}-0.025$ and $u_b={\xi }_{\text{macro}}^{\text{lattice}}+0.025,$ where 0.025 is the move limit step in the SLP algorithm. By using (Eq. 13(13) $${\xi }_{\text{macro}}^{\text{lattice}}{\xi }^{\text{macro}}={\xi }^{\text{bulk}}.$$(13) ), the corresponding values on ${\xi }^{\text{macro}}$ are [1.0000,0.6667, 0.5000,0.4000,0.3333,0.2857]. The optimal compliance c as function of ${\xi }_{\text{macro}}^{\text{lattice}}$ is plotted in Figure 5. The compliance c is normalized with the maximum value obtained for ${\xi }_{\text{macro}}^{\text{lattice}}$=0.2, denoted $c_{\text{max}}.$ The density solutions ${\gamma }_e^{\text{filt}}$ for elements with ${\rho }_e^{\text{heav}}>0.6$ are also plotted in Figure 5, this is done for ${\xi }_{\text{macro}}^{\text{lattice}}$=0.3, 0.4 and 0.7. These solutions are similar to the one presented in Borrvall and Petersson (2001). Observe that the bulk material is constant for all solutions presented in Figure 5. The same is true for the results presented next in Figures 6 and 7.

Figure 6. Optimal combinations of macro-layouts and local gradings for the L-shaped benchmark subjected to a single load case.

Figure 7. The weighted compliance for the L-shaped benchmark with and without local grading. A crisp macro-layout and the corresponding local grading are plotted for different views.

Next, optimal local grading is included by modifying the above setting by letting l_b=0.1 and u_b=0.8. The obtained optimal compliance as function of ${\xi }_{\text{macro}}^{\text{lattice}}$ is also plotted in Figure 5 using the same normalization in order to compare the results with the results obtained previously using no grading. A significant improvement in compliance is observed for ${\xi }_{\text{macro}}^{\text{lattice}}<$0.6. For ${\xi }_{\text{macro}}^{\text{lattice}}$=0.2, the improvement is almost a factor of 10. For ${\xi }_{\text{macro}}^{\text{lattice}}>$0.6, both curves are converging toward each other as expected. The results obtained for the optimal grading are studied in more detail in Figure 6, where the density solutions ${\gamma }_e^{\text{filt}}$ are plotted on elements with ${\rho }_e^{\text{heav}}>0.6$ as previous. For ${\xi }_{\text{macro}}^{\text{lattice}}$=0.4, the density plot is also sliced and the corresponding histogram of ${\gamma }_e^{\text{filt}}$ is depicted. It is clear for the single load case that the majority of all values on ${\gamma }_e^{\text{filt}}$ converge toward the lower and upper limit, l_b and u_b, respectively. This is explained by the material interpolation laws f_ij that penalize intermediate values in a similar way as SIMP and RAMP, see Figure 4. One should also observe the small difference obtained in the compliance for the different volumes of macro layout when optimal grading is included. The highest and lowest compliance values only differ 12.3%.

The L-shaped benchmark is also considered for the weighted compliance problem using $(F_x,F_y,F_z)$ with equal magnitudes and equal weights in the formulation of the compliance. That is, (Eq. 14(14) $$\begin{array}{c}\underset{(\mathit{\rho },\mathit{\gamma },\mathit{d})}{\mathrm{\text{min}}}c={\mathit{F}}^T\mathit{d}\hfill \\ \mathrm{s}.\mathrm{t}.\{\begin{array}{c}\mathit{K}(\mathit{\rho },\mathit{\gamma })\mathit{d}=\mathit{F},\hfill \\ {\mathrm{V}}^{\text{macro}}(\mathit{\rho })\le {\widehat{V}}^{\text{macro}},\hfill \\ {\mathrm{V}}^{\text{lattice}}(\mathit{\gamma })\le {\widehat{V}}^{\text{lattice}},\hfill \\ ϵ\mathbf{1}\le \mathit{\rho }\le \mathbf{1},\hfill \\ {\mathrm{l}}_{\mathrm{b}}\mathbf{1}\le \mathit{\gamma }\le u_b\mathbf{1},\hfill \end{array}\hfill \end{array}$$(14) ) is solved for three equilibrium equations $\mathit{K}{\mathit{d}}_i={\mathit{F}}_i$ (i = 1,2,3) and (28) $$c=\sum _{i=1}^3{\mathit{F}}_i^T{\mathit{d}}_i.$$(28)

Now, ${\xi }^{\text{bulk}}=0.1$ and solutions are generated for ${\xi }^{\text{macro}}\in$[0.25:0.05:0.5], i.e., [0.4,0.3333,0.2857,0.25,0.222,0.2]. Again, the solutions using optimal grading are significantly better than the solutions obtained for constant lattice density, see Figure 7. For ${\xi }_{\text{macro}}^{\text{lattice}}$=0.2, the optimal graded solution is approximately 5 times better, and, for ${\xi }_{\text{macro}}^{\text{lattice}}$=0.4, more than 2 times better. The densities for the solution with lowest compliance is also plotted in Figure 7. The optimal macro layout defined by ${\rho }_e^{\text{heav}}>$0.6 is plotted in green for three different views. The optimal grading ${\gamma }_e^{\text{filt}}$ is then plotted and sliced on this macro layout. In addition, the histogram of ${\gamma }_e^{\text{filt}}$ is presented at the lower right corner of Figure 7.

As can be seen from the plots in Figure 7, a crisp macro layout plotted in green is obtained, although the weighted compliance of three load cases is solved for only 10% bulk material. This might be explained by the intermediate values appearing in the local grading solution ${\gamma }_e^{\text{filt}}.$ Still, the number of elements with intermediate values is low, which can be seen from the histogram plot. Thus, most values of ${\gamma }_e^{\text{filt}}$ still converge toward the lower and upper limit, l_b and u_b, respectively.

For the optimal solution presented in Figure 7, the convergence of the weighted compliance normalized by the initial compliance value, the volume fractions in (Eq. 29(29) $$\begin{array}{c}{\widehat{\xi }}^{\text{bulk}}=\frac{V^{\text{bulk}}(\{{\rho }_e^{\text{heav}}>{\rho }_t\},\mathit{\gamma })}{V^{\text{bulk}}(\mathbf{1},\mathbf{1})},\hfill \\ {\widehat{\xi }}^{\text{macro}}=\frac{V^{\text{bulk}}(\{{\rho }_e^{\text{heav}}>{\rho }_t\},\mathbf{1})}{V^{\text{bulk}}(\mathbf{1},\mathbf{1})},\hfill \\ {\widehat{\xi }}_{\text{macro}}^{\text{lattice}}=\frac{V^{\text{bulk}}(\{{\rho }_e^{\text{heav}}>{\rho }_t\},\mathit{\gamma })}{V^{\text{bulk}}(\{{\rho }_e^{\text{heav}}>{\rho }_t\},\mathbf{1})}\mathrm{ }\text{for}\mathrm{ }{V}^{\text{bulk}}(\{{\rho }_e^{\text{heav}}>{\rho }_t\},\mathbf{1})>0.\hfill \end{array}$$(29) ) and the measure of discreteness introduced by Sigmund (2007) are plotted in Figure 8. In particular, the convergence of the volume fractions in (Eq. 12(12) $$\begin{array}{c}{\xi }^{\text{bulk}}=\frac{V^{\text{bulk}}}{V^{\text{design}}},\hfill \\ {\xi }^{\text{macro}}=\frac{V^{\text{macro}}}{V^{\text{design}}},\hfill \\ {\xi }_{\text{macro}}^{\text{lattice}}=\frac{V_{\text{macro}}^{\text{lattice}}}{V^{\text{macro}}},\hfill \end{array}$$(12) ) are studied for three different thresholds on ${\rho }_e^{\text{heav}}>{\rho }_t$=0.95, 0.6 and 0 such that the volume fractions in (Eq. 12(12) $$\begin{array}{c}{\xi }^{\text{bulk}}=\frac{V^{\text{bulk}}}{V^{\text{design}}},\hfill \\ {\xi }^{\text{macro}}=\frac{V^{\text{macro}}}{V^{\text{design}}},\hfill \\ {\xi }_{\text{macro}}^{\text{lattice}}=\frac{V_{\text{macro}}^{\text{lattice}}}{V^{\text{macro}}},\hfill \end{array}$$(12) ) are numerically evaluated as (29) $$\begin{array}{c}{\widehat{\xi }}^{\text{bulk}}=\frac{V^{\text{bulk}}(\{{\rho }_e^{\text{heav}}>{\rho }_t\},\mathit{\gamma })}{V^{\text{bulk}}(\mathbf{1},\mathbf{1})},\hfill \\ {\widehat{\xi }}^{\text{macro}}=\frac{V^{\text{bulk}}(\{{\rho }_e^{\text{heav}}>{\rho }_t\},\mathbf{1})}{V^{\text{bulk}}(\mathbf{1},\mathbf{1})},\hfill \\ {\widehat{\xi }}_{\text{macro}}^{\text{lattice}}=\frac{V^{\text{bulk}}(\{{\rho }_e^{\text{heav}}>{\rho }_t\},\mathit{\gamma })}{V^{\text{bulk}}(\{{\rho }_e^{\text{heav}}>{\rho }_t\},\mathbf{1})}\mathrm{ }\text{for}\mathrm{ }{V}^{\text{bulk}}(\{{\rho }_e^{\text{heav}}>{\rho }_t\},\mathbf{1})>0.\hfill \end{array}$$(29)

Figure 8. Convergence in weighted normalized compliance (0.0078), volume fractions and measure of discreteness. The corresponding histogram of ${\rho }_e^{\text{heav}}$ is also plotted.

Notice, that (Eq. 12(12) $$\begin{array}{c}{\xi }^{\text{bulk}}=\frac{V^{\text{bulk}}}{V^{\text{design}}},\hfill \\ {\xi }^{\text{macro}}=\frac{V^{\text{macro}}}{V^{\text{design}}},\hfill \\ {\xi }_{\text{macro}}^{\text{lattice}}=\frac{V_{\text{macro}}^{\text{lattice}}}{V^{\text{macro}}},\hfill \end{array}$$(12) ) is recovered from (Eq. 29(29) $$\begin{array}{c}{\widehat{\xi }}^{\text{bulk}}=\frac{V^{\text{bulk}}(\{{\rho }_e^{\text{heav}}>{\rho }_t\},\mathit{\gamma })}{V^{\text{bulk}}(\mathbf{1},\mathbf{1})},\hfill \\ {\widehat{\xi }}^{\text{macro}}=\frac{V^{\text{bulk}}(\{{\rho }_e^{\text{heav}}>{\rho }_t\},\mathbf{1})}{V^{\text{bulk}}(\mathbf{1},\mathbf{1})},\hfill \\ {\widehat{\xi }}_{\text{macro}}^{\text{lattice}}=\frac{V^{\text{bulk}}(\{{\rho }_e^{\text{heav}}>{\rho }_t\},\mathit{\gamma })}{V^{\text{bulk}}(\{{\rho }_e^{\text{heav}}>{\rho }_t\},\mathbf{1})}\mathrm{ }\text{for}\mathrm{ }{V}^{\text{bulk}}(\{{\rho }_e^{\text{heav}}>{\rho }_t\},\mathbf{1})>0.\hfill \end{array}$$(29) ) when ρ_t=0. Furthermore, by following Sigmund (2007), the measure of discreteness is defined by (30) $$M_{\text{nd}}=\sum _{e=1}^{n_{\text{el}}}\frac{4{\rho }_e^{\text{heav}}(1-{\rho }_e^{\text{heav}})}{n_{\text{el}}}\times 100\%.$$(30)

In Figure 8, the convergence of the volume fractions in (Eq. 29(29) $$\begin{array}{c}{\widehat{\xi }}^{\text{bulk}}=\frac{V^{\text{bulk}}(\{{\rho }_e^{\text{heav}}>{\rho }_t\},\mathit{\gamma })}{V^{\text{bulk}}(\mathbf{1},\mathbf{1})},\hfill \\ {\widehat{\xi }}^{\text{macro}}=\frac{V^{\text{bulk}}(\{{\rho }_e^{\text{heav}}>{\rho }_t\},\mathbf{1})}{V^{\text{bulk}}(\mathbf{1},\mathbf{1})},\hfill \\ {\widehat{\xi }}_{\text{macro}}^{\text{lattice}}=\frac{V^{\text{bulk}}(\{{\rho }_e^{\text{heav}}>{\rho }_t\},\mathit{\gamma })}{V^{\text{bulk}}(\{{\rho }_e^{\text{heav}}>{\rho }_t\},\mathbf{1})}\mathrm{ }\text{for}\mathrm{ }{V}^{\text{bulk}}(\{{\rho }_e^{\text{heav}}>{\rho }_t\},\mathbf{1})>0.\hfill \end{array}$$(29) ) are plotted for ${\rho }_e^{\text{heav}}>0.95,$ where ${\widehat{\xi }}^{\text{bulk}},\hspace{0.17em}{\widehat{\xi }}^{\text{macro}}$ and ${\widehat{\xi }}_{\text{macro}}^{\text{lattice}}$ converge to 0.0992, 0.247 and 0.402, respectively. By switching the threshold from ρ_t=0.95 to ρ_t=0.6 in (Eq. 29(29) $$\begin{array}{c}{\widehat{\xi }}^{\text{bulk}}=\frac{V^{\text{bulk}}(\{{\rho }_e^{\text{heav}}>{\rho }_t\},\mathit{\gamma })}{V^{\text{bulk}}(\mathbf{1},\mathbf{1})},\hfill \\ {\widehat{\xi }}^{\text{macro}}=\frac{V^{\text{bulk}}(\{{\rho }_e^{\text{heav}}>{\rho }_t\},\mathbf{1})}{V^{\text{bulk}}(\mathbf{1},\mathbf{1})},\hfill \\ {\widehat{\xi }}_{\text{macro}}^{\text{lattice}}=\frac{V^{\text{bulk}}(\{{\rho }_e^{\text{heav}}>{\rho }_t\},\mathit{\gamma })}{V^{\text{bulk}}(\{{\rho }_e^{\text{heav}}>{\rho }_t\},\mathbf{1})}\mathrm{ }\text{for}\mathrm{ }{V}^{\text{bulk}}(\{{\rho }_e^{\text{heav}}>{\rho }_t\},\mathbf{1})>0.\hfill \end{array}$$(29) ), one obtains instead 0.0995, 0.249 and 0.400, respectively, and if ${\rho }_e^{\text{heav}}>0$ in (Eq. 29(29) $$\begin{array}{c}{\widehat{\xi }}^{\text{bulk}}=\frac{V^{\text{bulk}}(\{{\rho }_e^{\text{heav}}>{\rho }_t\},\mathit{\gamma })}{V^{\text{bulk}}(\mathbf{1},\mathbf{1})},\hfill \\ {\widehat{\xi }}^{\text{macro}}=\frac{V^{\text{bulk}}(\{{\rho }_e^{\text{heav}}>{\rho }_t\},\mathbf{1})}{V^{\text{bulk}}(\mathbf{1},\mathbf{1})},\hfill \\ {\widehat{\xi }}_{\text{macro}}^{\text{lattice}}=\frac{V^{\text{bulk}}(\{{\rho }_e^{\text{heav}}>{\rho }_t\},\mathit{\gamma })}{V^{\text{bulk}}(\{{\rho }_e^{\text{heav}}>{\rho }_t\},\mathbf{1})}\mathrm{ }\text{for}\mathrm{ }{V}^{\text{bulk}}(\{{\rho }_e^{\text{heav}}>{\rho }_t\},\mathbf{1})>0.\hfill \end{array}$$(29) ), then one gets as expected 0.1, 0.25 and 0.4, respectively, which correspond to the prerequisites of the problem. The converged values obtained for these three different thresholds on ${\rho }_e^{\text{heav}}$ are indeed very similar. The explanation for this is that the smooth Heaviside filter in (Eq. 23(23) $${\rho }_e^{\text{heav}}=\frac{\mathrm{\text{tanh}}(\beta \eta )+\mathrm{\text{tanh}}(\beta ({\rho }_e^{\text{filt}}-\eta ))}{\mathrm{\text{tanh}}(\beta \eta )+\mathrm{\text{tanh}}(\beta (1-\eta ))},$$(23) ) efficiently removes almost all intermediate densities. This can also be seen by studying $M_{\text{nd}}$ and the histogram of ${\rho }_e^{\text{heav}}$ in Figure 8. $M_{\text{nd}}$ converges to 0.68%, which is just slightly larger than the values reported by Sigmund (2007) for the 2 D MBB test example using two different Heaviside filters. By studying the convergence plots, one clearly observes the effect of the Heaviside filter on the compliance, the volume fractions and $M_{\text{nd}}$ when it is activated after 50 iterations.

The optimal combinations of macro layouts ${\rho }_e^{\text{heav}}$ and local gradings ${\gamma }_e^{\text{filt}}$ can be converted to implicit-based geometries, which in turn can be exported to stl-files using marching cubes. A procedure for this is presented in (Strömberg 2019), which involves the following steps:

1. The macro-layout solution ${\rho }_e^{\text{heav}}$ is converted to an implicit-based surface $f^{\text{macro}}$ using a support vector machine (Strömberg 2019), where $f^{\text{macro}}\ge 0$ defines the macro layout.

2. The local grading solution ${\gamma }_e^{\text{filt}}$ is mapped to κ in (Eq. 25(25) $$f^{\text{lattice}}=\hspace{0.17em}\text{sin}\hspace{0.17em}(\omega x)\hspace{0.17em}\text{cos}\hspace{0.17em}(\omega y)+\hspace{0.17em}\text{sin}\hspace{0.17em}(\omega y)\hspace{0.17em}\text{cos}\hspace{0.17em}(\omega z)+\hspace{0.17em}\text{sin}\hspace{0.17em}(\omega z)\hspace{0.17em}\text{cos}\hspace{0.17em}(\omega x)+\kappa ,$$(25) ) using Shepard interpolations (Shepard 1968). Then, $f^{\text{lattice}}\ge 0$ defines the local graded lattice structure for the whole design domain.

3. $f_{\text{macro}}^{\text{lattice}}=\mathrm{\text{min}}(f^{\text{macro}},f^{\text{lattice}})\ge$0 defines the optimal combination of macro layout and local grading. This is plotted to the left in Figure 9 for the optimal solution presented previously. By zooming in, the grading appears clearly.

4. A thin cover for the lattice designs can be generated by $f^{\text{cover}}=\mathrm{\text{min}}(f^{\text{macro}},-(f^{\text{macro}}-\chi )),$ where χ is a scalar parameter setting the thickness of the cover.

5. The final engineering design is then obtained by $\mathrm{\text{max}}(f_{\text{macro}}^{\text{lattice}},f^{\text{cover}}).$ This is illustrated in Figure 9 for the L-shaped benchmark.

Figure 9. Stl-files of the optimal solution for the L-shaped benchmark by using implicit-based geometry. The cut of the design is done in order to show the lattice inside the cover.

By using the proposed multi-scale framework, ultra-lightweight design of graded TPMS-based structures can easily be generated. This is demonstrated for the L-shaped benchmark for the single load case with $l_b=0.1,\hspace{0.17em}{u}_b=0.7,\hspace{0.17em}{\xi }^{\text{macro}}$=0.1 and ${\xi }_{\text{macro}}^{\text{lattice}}$=0.35, implying ${\xi }^{\text{bulk}}$ as low as 0.035. The multi-scale results and the corresponding design as stl-files are presented in Figure 10. The relative lattice density ${\gamma }_e^{\text{filt}}$ is plotted on elements with ${\rho }_e^{\text{heav}}>$0.6 in the upper left corner of the figure. The corresponding graded frame-based Gyroid structure is plotted to the right. Again, the grading appears clearly by zooming. A cover is then generated and added to the lattice. The corresponding final design with cover is depicted sliced in the figure. The volume fraction of the cover is 2% taken separtely, but it is only adding 1.1% extra volume when the union of the lattice is performed. Thus, the total volume of bulk material is only 4.8% of the original design domain, which is considered to be an ultra-lightweight design, see e.g., Zhang, Sato, and Yanagimoto (2021). If no grading is included for this setting of the benchmark, then the compliance increases a factor of 2 and the macro layout becomes less crisp.

Figure 10. A ultra-lightweight design of the L-shaped benchmark for a single load case.

As an illustration of a real engineering application, the design domain of the established GE-bracket is optimized for a weighted load case of three orthogonal directions of forces acting on the hole of the pin and fixed at the four bolt holes, see Figure 11. The optimal combination of macro layout and local grading for ${\xi }^{\text{macro}}$=0.25 and ${\xi }_{\text{macro}}^{\text{lattice}}$=0.4, implying 10% bulk material, are solved using (Eq. 14(14) $$\begin{array}{c}\underset{(\mathit{\rho },\mathit{\gamma },\mathit{d})}{\mathrm{\text{min}}}c={\mathit{F}}^T\mathit{d}\hfill \\ \mathrm{s}.\mathrm{t}.\{\begin{array}{c}\mathit{K}(\mathit{\rho },\mathit{\gamma })\mathit{d}=\mathit{F},\hfill \\ {\mathrm{V}}^{\text{macro}}(\mathit{\rho })\le {\widehat{V}}^{\text{macro}},\hfill \\ {\mathrm{V}}^{\text{lattice}}(\mathit{\gamma })\le {\widehat{V}}^{\text{lattice}},\hfill \\ ϵ\mathbf{1}\le \mathit{\rho }\le \mathbf{1},\hfill \\ {\mathrm{l}}_{\mathrm{b}}\mathbf{1}\le \mathit{\gamma }\le u_b\mathbf{1},\hfill \end{array}\hfill \end{array}$$(14) ). The obtained optimal combinations of density solutions ${\rho }_e^{\text{heav}}$ and ${\gamma }_e^{\text{filt}}$ are presented in Figure 11. In particular, ${\gamma }_e^{\text{filt}}$ is plotted on elements with ${\rho }_e^{\text{heav}}>$0.6 and then sliced. The macro layout defined by ${\rho }_e^{\text{heav}}>$0.6 is also plotted in green. By transferring the density solutions to implicit based geometry as discussed above, the corresponding stl-file is generated and presented to the left in Figure 11. In addition, a thin cover is added using Boolean operations in order to a obtain a robust engineering design presented to the center and the right of Figure 11.

Figure 11. Optimal combination of macro-layout and local graded frame-based Gyroid structure for the well-known GE-bracket.

6. Concluding remarks

In this paper, a new innovative multi-scale topology optimization framework for generating optimal combinations of macro-layouts and local graded TPMS-based lattice structures is proposed and implemented. The formulation is done by introducing two density variables, one for defining the macro layout and the other one defining the relative density of the lattice structure. The density for the macro layout is governed by RAMP and the relative density of the lattice structure is governed by effective elastic properties obtained by numerical homogenization of representative volume elements of the TPMS-based lattice. By solving benchmarks in 3 D using the proposed implemented formulation, it is demonstrated that optimal designs of macro-layouts with graded Gyroid structures are generated that are significantly stiffer than the corresponding optimal designs using constant lattice densities. It is also shown that the difference in compliance for designs with different volume fractions of macro layouts is small when local grading is performed. Furthermore, the optimal macro-layouts are typically crisp even for multiple load cases, and by combining low volume fractions of macro-layout and lattice structure, designs with low volume of bulk material are easily generated. In conclusion, the proposed multi-scale topology optimization formulation for optimal combinations of macro layouts and local graded TPMS-based lattice structures is most promising for generating ultra-lightweight designs for additive manufacturing. Even if the stiffness becomes slightly lower than the stiffness obtained for a corresponding optimal solid structure without lattices, the optimal lattice design might be beneficial from other aspects such as buckling loads and frequency response. To include these kind of aspects in the presented framework is a topic of future research.

Disclosure statement

The author declares that he has no conflict of interest.

Additional information

Funding

This work was financed by Vinnova in the project RoboDAM.