Topology optimization based on reduction methods with applications to multiscale design and additive manufacturing

Abstract

Advanced manufacturing processes such as additive manufacturing offer now the capability to control material placement at unprecedented length scales and thereby dramatically open up the design space. This includes the considerations of new component topologies as well as the architecture of material within a topology offering new paths to creating lighter and more efficient structures. Topology optimization is an ideal tool for navigating this multiscale design problem and leveraging the capabilities of advanced manufacturing technologies. However, the resulting design problem is computationally challenging as very fine discretizations are needed to capture all micro-structural details. In this paper, a method based on reduction techniques is proposed to perform efficiently topology optimization at multiple scales. This method solves the design problem without length scale separation, i.e., without iterating between the two scales. Ergo, connectivity between space-varying micro-structures is naturally ensured. Several design problems for various types of micro-structural periodicity are performed to illustrate the method, including applications to infill patterns in additive manufacturing.

Key words： multiscale topology optimization; micro-structure; additive manufacturing; reduction techniques; substructuring; static condensation; super-element

Introduction

Manufacturing technologies are rapidly advancing and now offer the capability to control material placement at unprecedented length scales. This gives tremendous design freedom to designers. Whether being the lattice infill of an additively manufactured part or a standard structure, manufacturing constraints still do exist and must be accounted for during the design process [¹–⁶]. Nonetheless, these new manufacturing technologies offer the possibility to fabricate structures for which the macrostructural layout as well as the underlying micro-structure, or architecture of material have been specifically designed. That said, to fully leverage the capabilities of these new technologies, the design process must rigorously account for both scales to achieve the best performances.

Computational topology optimization provides a systematic, mathematically driven framework for navigating this new design challenge. Traditional implementations of structural optimization occur at the macroscale level, i.e., a selected material is distributed over a design domain such that an objective function, e.g., compliance or structural mass, is minimized while satisfying a set of constraints. To achieve the optimal material distribution, the homogenization method [⁷] or the SIMP interpolation model [⁸,⁹] are commonly used. Topology optimization is not restricted to the optimal distribution of a single material, but several materials can be considered [¹⁰,¹¹]. Topology optimization has also been employed to design novel materials with extreme properties [¹²,¹³] such as negative Poisson’s ratio [¹⁴,¹⁵], thermal expansion coefficient [¹⁶], fluid permeability [¹⁷,¹⁸], piezoelectric properties [¹⁹] and phononic properties [²⁰,²¹], to name a few. To tailor the effective properties of the designed material, the inverse homogenization method [²²] has been widely employed. Alternatively, in cases where homogenization may not apply, such as optimizing energy absorption considering material plasticity [²³,²⁴], one can create a finitely periodic representation of the material at significantly increased computational cost. When designing novel materials, the optimization process is usually carried out with a general idea of future potential needs rather than focusing on a specific application. Ergo, when applying the optimized material to a specific application or an existing structural design, the resulting design may not be optimal due to the uncoupled characteristic of the process. The idea of designing simultaneously the structure and the material is thus appealing since it would lead to a structure exhibiting an optimized macrostructural layout with architectured materials that account for the structure boundary conditions.

To capture micro-structural details, fine discretization is needed as it determines the resolution of the design optimization process. However, with increasing degrees of freedom (dofs), computational cost grows drastically and subsequently prevents the use of very fine resolutions. In order to consider multiscale optimization, a hierarchical approach has been developed wherein the macroscopic layout and the micro-structures are designed iteratively [²⁵]. The hierarchical approach relies on the homogenization method which determines the homogenized macroscopic response of a unit cell [²⁶,²⁷]. Homogenized responses are then used to optimize the macroscopic layout at a moderated computational cost. Moreover, due to the presence of a weak coupling between the two scales, the solution of the optimization problem at the microscale level can easily be parallelized [²⁸]. The hierarchical approach has also been employed to solve nonlinear problems [²⁹,³⁰]. More recently, the two-scale problem was linearized and a level set based topology optimization was adopted at both scales to design the structure and the material [³¹]. Again, the homogenization theory is used to bridge the scales.

For a single micro-structure repeated throughout the computational domain, it was demonstrated that an optimized periodic structure with sufficient finite geometric periodicity converges to optimized material unit cells with an infinite geometric periodicity obtained using homogenization [³²,³³]. However, it must be noted that this is true only for certain design problems, particularly for cases where localization due to boundary condition effects is not an issue and where the homogenized moduli to be optimized is dominant and known a priori. The scale effect of micro-structure upon the optimal topology solution was investigated and discussed in Ref. [³⁴].

A common issue of these two-scale optimization approaches using the homogenization theory concerns the lack of connectivity between unit cells, especially when the material design varies in space. The homogenization theory assumes the micro-structure is infinitely periodic, yet in practice optimization is used to actively vary micro-structure in space. A consequence of this inconsistency is that structural features and load paths can be disconnected, meaning regions that are predicted to offer large effective stiffness in reality offer none. Ensuring the connectivity between micro-structures is thus of utmost importance to ensure the numerically predicted performance of the optimized design. This topic has received more attention in recent years, notably in the contexts of free material optimization [³⁵] and functionally graded materials [³⁶]. More recently, a contrast-independent spectral conditioner [³⁷] based on the multiscale finite element method [³⁸] was applied to solve large structures with fully-resolved micro-structural details. In the latter approach, the authors avoid the length scale separation, i.e., the micro- and macroscale problems are treated in an integrated way, without iterations between the two scales.

In this paper, we propose an efficient, flexible and easy-to-implement method to perform topology optimization at multiple scales while ensuring the manufacturability of micro-structural details. The method relies on reduction techniques, also known as condensation or substructuring. The idea of using reduction methods was introduced in Ref. [³⁴] for the purpose of dealing with conventional designs of materials and structures in a unified way. In the proposed approach, the micro and macro-scales are fully coupled, i.e., the two-scale design problem is treated in an integrated way, without length scale separation. The use of reduction techniques has a double benefit. While they highly improve the computational efficiency, they also achieve the connection between space-varying micro-structures if appropriate boundary nodes are retained in the reduction process. It is also clear that the micro-structural features are designed with full awareness of the structural boundary conditions, as all applied loads are carried by distinct solid features. The idea for the condensation-based multiscale approach was first presented and highlighted by the authors of Ref. [³⁹]. Herein we expand on this initial work to include full derivations and algorithmic details, and provide extensions to consider additional cases of periodicity, infill design of existing components, and connections to length scale.

The layout of the paper is as follow. Section 2 briefly recalls the principle of reduction techniques and describes the proposed method to solve multiscale topology optimization problem. The optimization problem formulation and solution algorithm are discussed in Section 3. Several examples for various kinds of micro-structural periodicity are solved and the optimality of lattice infills is investigated for linear elastic material in Section 4. A discussion on the methods that decouple the micro— and macro-scales is provided in Sections 5 and 6 concludes the paper.

Multiscale topology optimization based on reduction techniques

Generalities of reduction techniques

Reduction techniques are commonly employed to reduce the size of a static or dynamic problem. The adopted reduction method originating from the work of Turner et al. [⁴⁰] is hereafter briefly summarized.

Let us consider the static equilibrium equation of a structure, expressed as

(1)\mathit{K}\mathit{d}=\mathit{F},

where \mathit{K} is the positive definite global stiffness matrix considering the free dofs, \mathit{d} is the vector of free nodal displacements, and \mathit{F} is the force vector composed of nodal applied forces and forces due to displacement boundary conditions. Resulting from a user-defined input, the static equilibrium Eq. (1) can be partitioned into two sets of dofs, named the retained (subscript r) and condensed (subscript c) dofs, leading to

(2)\left[\begin{array}{ll}{\mathit{K}}_{\mathrm{rr}}& {\mathit{K}}_{\mathrm{rc}}\\ {\mathit{K}}_{\mathrm{cr}}& {\mathit{K}}_{\mathrm{cc}}\end{array}\right]\left[\begin{array}{l}{\mathit{d}}_r\\ {\mathit{d}}_c\end{array}\right]=\left[\begin{array}{l}{\mathit{F}}_r\\ {\mathit{F}}_c\end{array}\right],

where we note {\mathit{K}}_{\mathrm{rr}} and {\mathit{K}}_{\mathrm{cc}} are both positive definite due to the positive definiteness of \mathit{K}.

The goal of the reduction method is to transform the original problem of dimension n, where n equals the number of dof, into a reduced size problem by temporarily removing all condensed dof. To achieve this reduction, the second equation of Eq. (2) is removed by expressing the displacement d_c, as follows:

(3){\mathit{d}}_c={\mathit{K}}_{\mathrm{cc}}^{-\mathrm{1}}({\mathit{F}}_c-{\mathit{K}}_{\mathrm{cr}}{\mathit{d}}_r).

Then, inserting that result into the first equation of Eq. (2), the condensed stiffness equation reads:

(4)\overline{\mathit{K}}{\mathit{d}}_r=\overline{\mathit{F}},

where \overline{\mathit{K}} and \overline{\mathit{F}} are named respectively the condensed stiffness matrix and the condensed force vector and they are defined as:

(5)\overline{\mathit{K}}={\mathit{K}}_{\mathrm{rr}}-{\mathit{K}}_{\mathrm{rc}}{\mathit{K}}_{\mathrm{cc}}^{-\mathrm{1}}{\mathit{K}}_{\mathrm{cr}},

(6)\overline{\mathit{F}}={\mathit{F}}_r-{\mathit{K}}_{\mathrm{rc}}{\mathit{K}}_{\mathrm{cc}}^{-\mathrm{1}}{\mathit{F}}_c.

From Eq. (4), the displacement vector of the retained dof {\mathit{d}}_r can be computed at a reduced cost. The displacements of the condensed dof {\mathit{d}}_c can be subsequently recovered if needed by using the Eq. (3).

Multiscale topology optimization based on reduction techniques

The simultaneous optimization of micro-structural details and structural layout requires dealing with very fine meshes since the latter determines the resolution of the optimization. However, the computational cost severally grows with increasing mesh size, leading to extreme computational demands. Improving computational efficiency and scalability has been the focus of many recent topology optimization works, including advanced preconditioning strategies [⁴¹,⁴²], approximation strategies [⁴³] and parallel processing implementations [⁴⁴], among others.

In this paper, a method based on reduction techniques is proposed to efficiently perform high resolution topology optimization, including simultaneous design of what can be classified as micro-structural unit cells and structure topology. To reduce the problem size, each unit cell of the structure is condensed in a super-element, i.e., each unit cell possesses its own condensed stiffness matrix {\overline{\mathit{K}}}_{\mathrm{se}} by applying Eq. (5). Likewise, the condensed load vector {\overline{\mathit{F}}}_{\mathrm{se}} can be applied by employing Eq. (6). Consequently, solving the macroscopic static equilibrium problem is performed at an extreme reduced cost. During the condensation process, the retained nodes are selected as the border nodes of the unit cells. This selection is paramount as it ensures the connectivity between space-varying micro-structures.

Figure 1 illustrates the proposed approach where a unit-cell is highlighted and where a two-level mesh can be observed. At the lower scale, the standard finite element mesh is defined within the super-element, while at the larger scale, the representative mesh is the coarser mesh composed of the super-elements that embed the lower scales.

Fig.1 Design domain discretization with a zoom on a unit-cell/super-element. H is the size of the square unit cell and h is the size of the finite element mesh. It must be noted that the method is not restricted to square unit cell.

Full size|PPT slide

In addition to a two-level mesh, the condensation process can be performed in chain leading to a multi-level reduction. Embedding super-elements within super-elements further enhances the computation speed. However, there is an optimal number of layers beyond which the creation time of super-elements outweighs the benefit of solving reduced systems and additional layers no longer benefit the overall computation time. Figure 2 illustrates a two-level reduction approach. When considering a multi-level reduction, the master unit cell can be defined at different level. Referring to Fig. 2, the master unit cell can either be defined at the first level (e.g., S{E}_{\mathrm{1},\mathrm{1}}) or at the second level (e.g., S{E}_{\mathrm{2},\mathrm{1}}). In this paper, the master unit cell is always defined at the highest level.

Fig.2 Tree structure of a multi-level reduction.

Full size|PPT slide

The optimization problem and solution algorithm

Regularization of the topology optimization problem

Topology optimization problems are known to be ill-posed and prone to checkerboard issues and mesh dependency. Several methods have been proposed to regularize the problem [⁴⁵,⁴⁶]. In this paper, the Heaviside projection method [⁴⁷] is adopted to stabilize the problem by imposing a minimum length scale on structural features.

Using the Heaviside projection method, the independent design variables are not the element volume fractions \mathit{\rho } but instead auxiliary variables \mathit{\phi }. These variables can be located at any point in space [⁴⁸] and are projected onto element space to determine the element volume fractions \mathit{\rho } that define topology. Herein we locate these variables at the nodes of the finite element mesh and the projection occurs over the minimum length scale radius {\mathit{r}}_{\min } via the following regularized Heaviside function:

(7){\rho }_e(\mathit{\phi })=\mathrm{1}-{\mathrm{e}}^{-\beta {\mu }_e(\mathit{\phi })}+{\displaystyle \frac{{\mu }_e(\mathit{\phi })}{{\phi }_{\max }}}{\mathrm{e}}^{-\beta {\phi }_{\max }},

where {\mu }_e is a proximity-based linear filtering of nodal design variables located within distance r_{\min } of the element centroid [⁴⁹], \beta dictates the curvature of the regularization, i.e., the projection function is linear with \beta =\mathrm{0} and approaches the Heaviside function as \beta approaches infinity, and e is the Euler’s number. The often used beta-continuation strategy is eliminated by modifying parameters of the optimizer [⁵⁰] and the upper bound {\phi }_{\max } is set to 1. The reader is referred to Ref. [⁴⁷] for further details.

Optimization problem formulation

Adopting the Heaviside projection methodology, standard topology optimization design problems can be formulated as

(8)\begin{array}{l}\underset{\mathit{\phi }}{\mathrm{minimize}C}={\overline{\mathit{L}}}^{\mathrm{T}}\overline{\mathit{d}}={\mathit{L}}^{\mathrm{T}}\mathit{d},\\ \mathrm{subject\; to}\overline{\mathit{K}}(\phi )\overline{\mathit{d}}=\overline{\mathit{F}},\\ {\displaystyle \underset{e\in \mathrm{\Omega }}{\mathrm{\Sigma }}\frac{{\rho }_e(\mathit{\phi })v_e}{V_{\mathrm{\Omega }}}}\le V_{\mathrm{f}},\\ \mathrm{0}\le {\phi }_i\le {\phi }_{\max },\forall i\in \mathrm{\Omega },\end{array}

where \Omega is the design domain, V_{\Omega } is the volume of the design domain, V_{\mathrm{f}} is the allowable volume fraction of material, and v_e is the elemental volume. The vectors \overline{\mathit{F}} and \overline{\mathit{d}} result respectively from the assembly of the condensed applied nodal load vectors and the unknown retained nodal displacement vectors for each unit cell. The matrix \overline{\mathit{K}} is the global stiffness matrix assembled from the super-element stiffness matrices {\overline{\mathit{K}}}_{\mathrm{se}}. \overline{\mathit{L}} is a problem-dependent vector and will be specified in Section 4. The absence of the symbol \overline{•} over \mathit{F}, \mathit{d}, \mathit{K} and \mathit{L} is used to refer to them as the standard elements, i.e., the assembled elements without reduction methods. Since the reduction method obeys the energy conservation principle, the objective function C can either be defined with respect to the condensed or the standard elements.

We emphasize that the design variable vector \mathit{\phi } lies at the lowest scale and that it defines the topology of the unit-cells. This means that the micro-scale topology variables are optimized while considering the condensed, macroscale response in Eq. (8).

To solve the optimization problems Eq. (8), the method of moving asymptotes [⁵¹] is employed. To mathematically motivate a 0-1 material distribution, the solid isotropic material with penalization (SIMP) method is considered [⁸,⁹]. Hence, before performing the reduction process, the element stiffness matrices are penalized as follows:

(9){\mathit{K}}_e=({\rho }_e^p(\mathit{\phi })+{\rho }_{e,\min }){\mathit{K}}_{e,\mathrm{0}},

where p\ge \mathrm{1} is the exponential penalty term, {\mathit{K}}_{e,\mathrm{0}} is the element stiffness matrix and {\rho }_{e,\min } is a small positive number preventing singularity of the global stiffness matrix, typically {\rho }_{e,\min }={\mathrm{10}}^{-\mathrm{4}}.

Gradient-based methods require a sensitivity analysis to compute the gradient of the cost and constraint functions. The sensitivities of a function f with respect to the independent design variables {\phi }_i are given by

(10){\displaystyle \frac{df}{d{\phi }_i}}={\displaystyle \underset{e\in \mathrm{\Omega }}{\mathrm{\Sigma }}\frac{\partial f}{\partial {\rho }_e}}{\displaystyle \frac{d{\rho }_e}{d{\phi }_i}}.

The term d{\rho }_e/d{\phi }_i results from the Heaviside projection method and its derivation is detailed in Ref. [⁴⁷]. The term \partial f/\partial {\rho }_e is problem-dependent and the adjoint method is employed to compute its value, yielding

(11){\displaystyle \frac{\partial f}{\partial {\rho }_e}}=-p{\rho }_e^{p-\mathrm{1}}(\mathit{\phi }){\mathit{\lambda }}_e^{\mathrm{T}}{K}_{e,\mathrm{0}}{d}_e,

where \lambda is the adjoint response vector. We note that the sensitivity analysis is carried out at the microscale level since it involves the design variables \phi that are defined at the microscale level. Therefore, the sensitivity analysis requires recovering the condensed displacements using Eq. (3) whereupon the sensitivity is evaluated in a classical manner.

The flowchart of the optimization process is given in Fig. 3, where a light grey box highlights the novelty of this paper.

Fig.3 Flowchart of the optimization process.

Full size|PPT slide

Examples

The proposed approach is demonstrated on the four design problems illustrated in Fig. 4. The first three are compliance minimization problems and since this problem is equivalent to minimizing the external work, it follows that \mathit{L}=\mathit{F} and for sensitivity analysis, that{\mathit{\lambda }}_e={\mathit{d}}_e. The fourth example considers the design of an inverter compliant mechanism wherein the goal is to maximize the negative displacement at the output point for a given force at the input point. It follows that the vector \mathit{L} is a dummy load vector composed of zeros except at the entry associated with the dofs of the output point displacement whose value is set to 1. For the sensitivity analysis, the element vector {\mathit{\lambda }}_e is recovered from the condensed adjoint response vector \overline{\mathit{\lambda }} that solves \overline{\mathit{K}}(\mathit{\phi })\overline{\mathit{\lambda }}=\overline{\mathit{L}}. The inverter problem is further described in Refs. [⁴⁶,⁵²].

Fig.4 Design problems. (a) Double clamped beam; (b) cantilever beam with a distributed load; (c) Michell beam; (d) invert compliant mechanism.

Full size|PPT slide

All the design problems consider a material with a Young’s modulus of E=\mathrm{1} and Poisson’s ratio of \nu =\mathrm{0.3}, which is initially uniformly distributed over the design domain. Four-node bilinear elements (Q4), with the plane stress assumption and a unit thickness, are employed to mesh the structure. A continuation method is applied to the exponent parameter p of the SIMP law where p is increased from 1 to 5 in unit increments. The examples typically consider a one-level reduction process, i.e., super-elements are not embedded within super-elements except when explicitly mentioned. The optimization processes are deemed converged when the relative change of the objective function and the relative constraint violation are respectively less than {\mathrm{10}}^{-\mathrm{4}} and {\mathrm{10}}^{-\mathrm{6}} between iterations. We note that for intermediate steps of the continuation process, the convergence criteria on the relative change of the objective function is relaxed to {\mathrm{10}}^{-\mathrm{2}}.

The reduction approach gives total freedom to treat different configurations of the unit cell periodicity. The designer may decide to design specific unit cells and to replicate these unit cells through the structure (Fig. 5). A first possibility is to design a single unit cell, i.e., one master cell, and to replicate the optimized micro-structure throughout the whole structure, leading to a fully periodic structure (Fig. 5(a)). The designer can also create layered structures by defining several master cells and replicating each optimized micro-structure through a pre-defined number of layers (Fig. 5(b)) [³⁷]. This configuration offers the possibility of defining graded material structures by assigning varying properties to the master unit cells. When designing the selected unit cell topology, features at the boundaries have a minimum length scale of half of the imposed r_{\min } with the adopted projection method. To circumvent this issue, the designer can account for the neighborhood of the master unit cells when performing the projection. The last configuration considers the fully free-form design case where each unit cell is a master cell. Hence, each cell can have a different micro-structure and the optimizer has a total freedom to distribute material over the design domain.

The first two configurations lead to major savings as the proposed method only needs to condense and store the stiffness matrix of the few master cells which contain the independent design variables. In the third configuration, each cell must be condensed independently and therefore computational savings result only in solving the large system of equilibrium equations as a series of smaller systems.

Fully periodic structure design

The first numerical application is carried out on the double clamped beam design problem (Fig. 4(a)) for the fully periodic case (Fig. 5(a)) with the following parameters: L=2, B=1, F=1, h=H/40, r_{\min }=2h, V_{\mathrm{f}}=\mathrm{0.5}, and \beta =\mathrm{50}. The design problem is performed for an increasing number of unit cells and the different scenarios are characterized by the parameter M=B/H representing the number of unit cells along the vertical axis.

Fig.5 Three different types of periodicity. (a) Fully periodic structure; (b) user-defined selective periodicity; (c) fully free-form periodicity (no imposed periodicity). The unit cells containing design variables are shaded and letters identify the unit cell micro-structures.

Full size|PPT slide

The optimized solutions are presented in Fig. 6 where it is observed that the optimized micro-structures converge and that the topologies are similar for M\ge \mathrm{4}. Indeed, the proposed method implicitly factors in the size of the unit cell. A key point of emphasis is that structural features are present at the location of the applied load (the center point of the macroscale beam) in all cases. This is a key difference from using hierarchical methods where a unit cell is designed using inverse homogenization, which lacks awareness of the macroscale loads and boundary conditions. Figure 7 illustrates the optimized structures for M=\{\mathrm{2},\mathrm{4},\mathrm{64}\}. The results are in good agreement with Refs. [³³,³⁷] although the optimization methods slightly differ. As observed in these two references, the compliance increases with an increase in the number of unit cells. This augmentation stems from the size of the unit cell, as reducing the size confines the design space and thus, the optimizer has less freedom to minimize the objective function.

Fig.6 Optimized micro-structures of the double clamped beam design problem. (a) M=2, C=9.24; (b) M=4, C=11.13; (c) M=8, C=12.29; (d) M=16, C=14.02; (e) M=64, C=16.75.

Full size|PPT slide

Fig.7 Optimized designs of the double clamped beam. (a) M=2; (b) M=4; (c) M=64.

Full size|PPT slide

Continuing along with a fully periodic structure, the cantilever beam design problem (Fig. 4(b)) is performed for several values of M with the following parameters: L=2, B=1, F=0.01 per unit length, h= H/40, r_{\min }=2h, V_{\mathrm{f}}=\mathrm{0.5}, and \beta =\mathrm{64}. The optimized structures with a zoom on the associated micro-structure are illustrated in Fig. 8. This example is particularly interesting as it again highlights a benefit of using an integrated method without length scale separation. It is observed that the micro-structures possess a vertical bar on the right edge to support the distributed load on the right beam edge. When decoupling the scales of the design problem, the optimization process fails to capture accurately the loading applied at the macroscale level, requiring non-design subdomains to be enforced to circumvent this issue (see for instance Ref. [³³]). In contrast, this issue is handled directly with the proposed method. As previously observed, the compliance also increases as the unit cell size decreases.

Fig.8 Optimized micro-structures and structures of the cantilever beam design problem, where M is the number of unit cells along the vertical axis. (a) M=2, C=1.098×10^–2; (b) M=4, C=1.240×10^–2; (c) M=8, C=1.298×10^–2; (d) M=16, C=1.336×10^–2.

Full size|PPT slide

Efficiency of the reduction approach

To demonstrate the efficiency of the developed method, the results of the double clamped beam design problem with M=64 are further investigated. This problem can be classified as large scale problem since the model consists of 13107200 finite elements with a total number of dof equal to 26229762. The computation has been performed on a cluster equipped with AMD Opteron Processor 6220 CPUs, 132 GB memory, running on CentOS 6.7 and Matlab 2016a. The study has been carried out with Matlab restricted and not restricted to a single computational thread.

Using the proposed method, the problem size is reduced by 95%, i.e., it only consists of 1309698 retained dofs, and a single design iteration with a one-level reduction takes around 41 min. Solving this problem without the reduction method must be possible but the computation time for a single iteration took so long (more than one day) that the design process was aborted.

Multi-level reduction can be employed to lessen even more the computation time. Table 1 summarizes the results where the mesh of the smallest entity as well as the number of super-elements at each level are specified. As the master unit cell is always defined at the highest reduction level in this paper, the four treated examples all have M=64. Using multi-thread computation, the computation time can be reduced using up to 3 reduction levels. However, an additional reduction level is not beneficial as the super-elements become so small the construction and the communication costs outweigh the gains in solving smaller systems of equations. It is possible that a more refined mesh or different computational implementation would lead to additional reduction levels being advantageous. Using a single computational thread, the average computation time is less sensitive to multi-level reduction.

Tab.1 Computation time using multi-level reduction for the double clamped beam having 26229762 dofs

Reduction level	Mesh: Smallest entity	Number of SE per level along x-axis	Number of SE per level along y-axis	Multi thread: Average time per iteration/min	Single thread: Average time per iteration/min
1	40×40	[128]	[64]	41	42
2	20×20	[2;128]	[2;64]	34	41
3	10×10	[2;2;128]	[2;2;64]	34	41
4	5×5	[2;2;2;128]	[2;2;2;64]	37	43

Note: Columns 3 and 4 indicate the number of super-elements embedded within each super-element for each reduction level. The last number corresponds to the highest level and is thus the number of super-elements composing the reduced global stiffness matrix.

For other periodicities with several master unit cells, the computation time can also be improved by parallelizing the creation of super-elements. It is worth noticing that the proposed method can be easily implemented in an existing topology optimization code. The main modification concerns the creation of a like-tree structure containing the super-element hierarchy for solution of the finite element equations.

The proposed method is an “exact” method and thus it solves micro-structural details to full accuracy. Therefore, the method may be slower compared to approximate methods, such as the method using a spectral coarse basis preconditioner [³⁷] for instance. For the double clamped beam design with M=64, the proposed method is 5 times slower, i.e., 42 min per design iteration compared to 8 min in Ref. [³⁷]. Time comparison should be taken only as an indicative basis since the code implementation and computer hardware are different. However, the proposed method seems to achieve convergence in a relatively small number of iterations. Adopting the same stopping criteria on the objective function ({\mathrm{10}}^{-\mathrm{3}}), the proposed method requires 48 design iterations to achieve convergence compared to 200 in Ref. [³⁷]. With that in mind, the method is only two times slower for a full accuracy. The proposed reduction method could be faster if one could update the condensed solution in a cheap way during the optimization process. However, this would probably introduce approximation and thus, resulting in a compromise between speed and accuracy.

Layered structure design

Following the examples in Ref. [³⁷], a structure with a layered micro-structure periodicity is now considered to solve the cantilever beam design problem (Fig. 4(b)). The parameters given in Section 4.1 for this design problem are adopted.

The optimized designs are illustrated in Fig. 9 for different numbers of unit cells as well as for different layer thicknesses. It can be observed that the connectivity between space-varying micro-structures is correctly ensured. Also, the vertical feature on the right edge of the unit cell is still well generated to support the distributed load applied on the right edge of the beam.

Fig.9 Cantilever beam optimized designs-Layered micro-structure periodicity. (a) M=4, Thickness=\left\{{\displaystyle \frac{\mathrm{1}}{\mathrm{4}}},{\displaystyle \frac{\mathrm{1}}{\mathrm{4}}},{\displaystyle \frac{\mathrm{1}}{\mathrm{4}}},{\displaystyle \frac{\mathrm{1}}{\mathrm{4}}}\right\}, C=8.672×10^–3; (b) M=16, Thickness=\left\{{\displaystyle \frac{\mathrm{4}}{\mathrm{16}}},{\displaystyle \frac{\mathrm{8}}{\mathrm{16}}},{\displaystyle \frac{\mathrm{4}}{\mathrm{16}}}\right\}, C=1.071×10^–2; (c) M=18, Thickness=\left\{{\displaystyle \frac{\mathrm{3}}{\mathrm{18}}},{\displaystyle \frac{\mathrm{12}}{\mathrm{18}}},{\displaystyle \frac{\mathrm{3}}{\mathrm{18}}}\right\}, C=1.145×10^–2; (d) M=18, Thickness=\left\{{\displaystyle \frac{\mathrm{6}}{\mathrm{18}}},{\displaystyle \frac{\mathrm{6}}{\mathrm{18}}},{\displaystyle \frac{\mathrm{6}}{\mathrm{18}}}\right\}, C=1.064×10^–2.

Full size|PPT slide

Fully free-form topology optimization

In the previous sections, the micro-structure periodicity is prescribed and consequently, micro-structural details are obtained. In this section, the fully free-form periodicity, i.e., the most general case, is investigated (Fig. 5(c)).

The double clamped beam design problem is solved with the following parameters: L=2, B=1, F=1, r_{\min }=2h, V_{\mathrm{f}}=\mathrm{0.3}, and \beta =\mathrm{64}. In Fig. 10(a) with M=16 and h=H/10, a classical result is recovered where we note the absence of distinct porous micro-structures. The same design problem is performed with a finer mesh to capture smaller micro-structural details if any (M=16 and h=H/40). The computation is first carried out with r_{\min }=8h to have the same minimum length scale as in Fig. 10(a). Figures 10(a) and 10(b) illustrate similar optimized designs while the objective function of the finer mesh is a little bit higher due to the finer discretization. Reducing the filtering radius size to r_{\min }=2h, the optimized design exhibits a few smaller features (Fig. 10(c)). However, this is attributed to the change of length scale rather than to the creation of micro-structures.

Fig.10 Optimized designs for the double clamped beam (a, b and c) and the compliant inverter mechanism (d) without imposing any periodicity on the unit cell. (a) M=16, h=H/10, C=9.02; (b) M=16, h=H/40, C=9.87, r_{\min }=\mathrm{8}h; (c) M=16, h=H/40, C=9.66; (d) M=10, h=H/12, C_{\mathrm{pen}}=-\mathrm{2.046}, C_{\mathrm{unpen}}=-\mathrm{0.236}.

Full size|PPT slide

These designs offer significantly improved structural efficiency as the compliances are slightly below the best periodic solutions (Fig. 5(a)) even though using significantly less volume of material (30% compared to 50%). We note these findings support the logic that optimal micro-structures are fully solid or fully void when optimizing elastic stiffness, and that porous micro-structures such as lattices should only arise when considering design objectives other than elastic stiffness, such as buckling or functional properties [³⁹,⁴⁴,⁵³].

To further investigate the fully free-form case, the challenging design problem of the inverter compliant mechanism is solved with the following parameters: L=120, F=1, k_{\mathrm{in}}=\mathrm{1}, k_{\mathrm{out}}=\mathrm{0.001}, r_{\min }= 2h, V_{\mathrm{f}}=\mathrm{0.25}, \beta =\mathrm{64}, M=10, and h=H/12. Again, a classical topology optimization result is recovered with an objective function of C_{\mathrm{pen}}=-\mathrm{2.046} (Fig. 10(d)). Considering unpenalized element density, i.e., p=\mathrm{1}, the objective function of the final design is C_{\mathrm{unpen}}=-\mathrm{0.236}. As previously noticed, the optimized design indicates a lack of micro-structural details or infill patterns.

To validate the non-optimality of porous infill patterns, these two examples are solved with an additional local volume constraint for each unit cell. This is similar in concept to the maximum length scale constraint in Ref. [⁵⁴], which when relaxed slightly was shown to effectively generate infill-like patterns [⁵⁵]. Here, we impose this constraint on each super-element, rather than each element, providing convenient implementation and reducing the number of local constraints, albeit resulting in a different restriction than the preceding works. These constraints prevent a super-element from being fully filled of material and thus promote infill pattern. Mathematically, the design problem Eq. (8) is supplemented with the following constraints:

(12){\displaystyle \sum _{e\in {\mathrm{\Omega }}_{\mathrm{se},k}}\frac{{\rho }_e(\mathit{\phi })v_e}{V_{\mathrm{se},k}}}\le LVF,\hspace{1em}k=\mathrm{1},\mathrm{2},\mathrm{...},n_{\mathrm{se}},

where n_{\mathrm{se}} is the number of super-elements, {\mathrm{\Omega }}_{\mathrm{se},k} is the design domain of the super-element k, V_{\mathrm{se},k}is the volume of the super-element k, and LVF is the allowable local volume fraction of material. The results for the double clamped beam and for the inverter compliant mechanism design problems are respectively illustrated in Figs. 11 and 12 for several LVF values. For both design problems, it is observed that local volume constraints promote the appearance of lattice-like features that negatively impact the stiffness-based objective function, the smaller the local volume fraction, the worse the objective function value.

Fig.11 Optimized design of the double clamped beam considering local volume constraints to promote lattice-like features (M=16, h=H=40). (a) C=9.95, LVF=0.9; (b) C=10.33, LVF=0.7; (c) C=10.91, LVF=0.5.

Full size|PPT slide

Fig.12 Optimized design of the compliant inverter considering local volume constraints (M=10, h=H/12). (a) C_{\mathrm{pen}}=-\mathrm{1.973}, C_{\mathrm{unpen}}=-\mathrm{0.537}, LVF=0.9; (b) C_{\mathrm{pen}}=-\mathrm{1.901}, C_{\mathrm{unpen}}=-\mathrm{0.355}, LVF=0.8; (c) C_{\mathrm{pen}}=-\mathrm{1.733}, C_{\mathrm{unpen}}=-\mathrm{0.215}, LVF=0.7.

Full size|PPT slide

Imposing minimum and maximum local volume constraints

In order to promote infill patterns, local volume constraints have been incorporated in the design problem to prevent the super-elements from being completely filled with material [⁵⁵].

In this section, in addition to the previous maximum local volume constraints, minimum local volume constraints are also considered to force a more uniform distribution of material within the design domain. These new constraints impose a minimum volume of material for each super-element and they are mathematically expressed as

(13){\displaystyle \sum _{e\in {\mathrm{\Omega }}_{\mathrm{se},k}}\frac{{\rho }_e(\mathit{\phi })v_e}{V_{\mathrm{se},k}}}\ge LV{F}_{\mathrm{m}},\hspace{1em}k=\mathrm{1},\mathrm{2},\mathrm{...},n_{\mathrm{se}},

where LV{F}_{\mathrm{m}} is the minimum volume fraction of material required in each unit cell.

The double clamped beam design problem is solved for several LV{F}_{\mathrm{m}} with the following parameters: L=2, B=1, F=1, r_{\min }=2h, V_{\mathrm{f}}=\mathrm{0.3}, LVF=\mathrm{0.5}, and \beta =\mathrm{64}. As the LV{F}_{\mathrm{m}} value increases, the material is better distributed within the design domain and finer features are generated while minimum length scale is still satisfied (Fig. 13). Here again, the more the design freedom is constrained, the larger the magnitude of the objective function. To satisfy the local constraints on the minimum amount of material LV{F}_{\mathrm{m}}, “light grey areas” may exist in the optimized designs (Figs. 13(a) and 13(c)). Those areas are not optimal for minimizing the compliance but they stem from the optimization problem formulation. In these cases, the algorithm has chosen to waste material to satisfy the constraints rather than spend more material resource to create a distinct but structurally inefficient feature.

The Michell beam design problem is now performed for several LVF with the following parameters: L=40, B=25, F=1, r_{\min }=2h, V_f=\mathrm{0.5}, LV{F}_{\mathrm{m}}=\mathrm{0.1} and \beta =\mathrm{64}. For large LVF values, the optimizer tends to generate a topology very similar to the classical topology optimization design of the Michell beam (Fig. 14(a)). “Grey areas” also exist in the very low strain energy regions of the optimized designs to satisfy the constraints on the minimum amount of material within each super-element. It is clear this material offers no structural benefit and thus SIMP penalization is ineffective in these regions. If material is required to exist in these regions, additional penalization or a metric related to this requirement is needed. Reducing the LVF value, the material is better distributed within the design domain (Fig. 14). With LVF=\mathrm{0.5}, the optimized design fills the macroscale space fully (Fig. 14(c)). One can observe two blurred areas at the top and bottom of the right edge, where strain energies are near zero. As already observed, the compliance increases as the LVF becomes more and more restrictive. Similar observations about the non-optimality of Michell-type truss structure are pointed out in Ref. [⁵³].

Fig.13 Optimized designs of the double clamped beam with minimum and maximum local volume constraints (M=8, h=H/40). (a) C=10.685, LVF_m=0.05; (b) C=11.431, LVF_m=0.2; (c) C=14.717, LVF_m=0.3.

Full size|PPT slide

Fig.14 Optimized designs of the Michell beam with minimum and maximum local volume constraints (M=20, h =H/10). (a) C=43.804, Max LVF=0.9; (b) C=46.785, Max LVF =0.7; (c) C=64.112, Max LVF=0.5.

Full size|PPT slide

Re-design of the inner structure of a part

During the lifetime of a system, it may occur that some components have been over-designed and that the replacement of those components with porous designs that can now be realized by additive manufacturing would be of great benefits. Optimization techniques can be used to improve the over-designed components but more constraints usually exist as the new design should fit within the existing system. This situation is here considered where we assume the doubled clamped beam design illustrated in Fig. 10(b) turns out to be over-designed. Hereafter, several lighter designs are proposed. The design domain consists of the inner part of the component and it is required to remain within the geometric envelope of the existing component, that can be observed in Fig. 15 (thin contour line).

The fully free-form approach is adopted to re-design the inner part of the component and thus infill patterns can be generated. The optimization is performed for several reduced allowable volume fractions of material with the same parameters as for the initial design. The optimized designs are illustrated on the left column of Fig. 15. The same optimization processes are performed for a smaller minimum length scale and are illustrated on the right column of Fig. 15. It can be seen that infill patterns are not promoted and that, as expected, the structure becomes less stiff as the allowable amount of material is reduced.

Fig.15 Double clamped beam re-design (M=16, h=H/40, r_min=8h for left column and r_min=2h for right column). (a) C=11.20, V_f=0.25; (b) C=11.06, V_f=0.25; (c) C=14.34, V_f=0.2; (d) C=13.92, V_f=0.2; (e) C=20.00, V_f=0.15; (f) C=18.83, V_f=0.15.

Full size|PPT slide

Discussion on two-step approaches

The SIMP interpolation law is generally adopted to force the design towards a 0-1 solution. Indeed, with p\ge \mathrm{1}, intermediate densities are “uneconomical” since the stiffness for intermediate densities is small with respect to the cost of material.

The issue of interpreting “grey” material is dodged when the optimization process actually achieves a black-and-white design. However, the issue of constructing a meta-material that mimics the properties of “grey” material is often raised when the optimized design is not totally clear of grey. Several multiscale optimization methods rely on the use of small values of p where-upon the methods take advantage of grey areas to artificially introduce micro-structures. For instance, a standard topology optimization problem is first performed for the macroscale design with p=1. Subsequently, the “grey” areas are replaced by pre-defined and convenient isotropic micro-structures.

We would like to briefly go back on the construction of a material model mimicking the SIMP interpolation model. Considering isotropic linear elastic materials, Hashin and Shtrikman [⁵⁶,⁵⁷] established the limits of possible isotropic material properties that can be achieved by constructing materials with micro-structure. The limits are usually referred to as the Hashin–Shtrikman bounds. The conditions under which the SIMP interpolation model respect the Hashin–Shtrikman bounds are clearly derived in Ref. [¹⁰] and thus the conditions that SIMP may represent a realizable physical model. For example, p\ge \mathrm{3} when using base materials with Poisson’s ratio of 1/3.

Consequently, it is physically impossible to identify micro-structures mimicking the “grey” region properties when using p=1 since this leads to micro-structures that are beyond the Hashin–Shtrikman bounds. Ergo, these two-step methods introduce a clear decoupling between the micro- and macroscopic scales and the resulting approach is not consistent in such cases.

Concluding remarks

Reduction techniques have been employed to perform topology optimization. Condensing the unit cells in super-elements reduces drastically the size of the design problem while still solving micro-structural details to full accuracy without length scale separation. The proposed method is flexible and any periodicity of the micro-structure can be easily enforced. Fully periodic and layered structures lead to major computation time saving as only a couple of unit cells must be condensed. Also, the proposed method can be easily implemented in existing topology optimization codes.

Several classical examples have been solved to illustrate the capabilities of the method. From the results, it appears that micro-structure and infill patterns do not seem to be optimal for stiffness-based design problems with linear elastic material when a finite length scale is considered, agreeing with findings in Ref. [⁵³]. Finite length scale is to be attributed since it is related to the resolution of the manufacturing process. Infill patterns appear here only when restricting the allowable volume of material within each super-element. Compared to other methods that introduce an additional constraint per finite element to promote infill pattern (see Ref. [⁵⁸] for instance), the proposed approach offers an alternative that uses less constraints and is synergistic with the reduction approach. Introducing local volume constraints was originally introduced in Ref. [⁵⁴] for imposing a maximum length scale on the design problem. These updated local volume constraint approaches therefore also influence maximum length scale, though in a less direct manner.

Adopting an integrated approach and fully resolving micro-structural details, the optimization process promotes “solid” unit cells for stiffness-based design problems with linear elastic material when imposing a finite length scale. However, for more advanced design problems, infill patterns seem to give better performances such as increasing the critical buckling load [⁵⁵] or the structural robustness [⁵⁸]. In the latter two studies, the improvements are due to positive side effects. For future work, it would be interesting to identify and formulate a design problem that generates infill pattern as the optimal result of the design process.

Acknowledgements

This research was supported in part by the US National Science Foundation (NSF) under Grant Number 1462453, and in part by the National Aeronautics and Space Administration (NASA) under Grant Number 80NSSC18K0428. Any opinions, findings, and conclusions or recommendations expressed in this article are those of the author(s) and do not necessarily reflect the views of NSF or NASA.

References

1 Gao W, Zhang Y, Ramanujan D, The status, challenges, and future of additive manufacturing in engineering. Computer Aided Design, 2015, 69: 65–89

2 Mumtaz K A, Vora P, Hopkinson N. A method to eliminate anchors/supports from directly laser melted metal powder bed processes. In: Proceedings of the Solid Freeform Fabrication Symposium. Austin, 2011, 55–64

3 Hussein A, Hao L, Yan C, Advanced lattice support structures for metal additive manufacturing. Journal of Materials Processing Technology, 2013, 213(7): 1019–1026

4 Leary M, Merli L, Torti F, Optimal topology for additive manufacture: A method for enabling additive manufacture of support-free optimal structures. Materials & Design, 2014, 63: 678–690

5 Gaynor A T, Meisel N A, Williams C B, Topology optimization for additive manufacturing: Considering maximum overhang constraint. In: Proceedings of the 15th AIAA/ISSMO Multidisciplinary Analysis and Optimization Conference. Atlanta: AIAA, 2014, 1–8

6 Gaynor A T, Guest J K. Topology optimization considering overhang constraints: Eliminating sacrificial support material in additive manufacturing through design. Structural and Multidisciplinary Optimization, 2016, 54(5): 1157–1172

7 Bendsøe M P, Kikuchi N. Generating optimal topologies in structural design using a homogenization method. Computer Methods in Applied Mechanics and Engineering, 1988, 71(2): 197–224

8 Bendsøe M P. Optimal shape design as a material distribution problem. Structural Optimization, 1989, 1(4): 193–202

9 Rozvany G I N, Zhou M, Birker T. Generalized shape optimization without homogenization. Structural Optimization, 1992, 4(3‒4): 250–252

10 Bendsøe M P, Sigmund O. Material interpolation schemes in topology optimization. Archive of Applied Mechanics, 1999, 69(9‒10): 635–654

11 Gaynor A T, Meisel N A, Williams C B, Multiple-material topology optimization of compliant mechanisms created via PolyJet three-dimensional printing. Journal of Manufacturing Science and Engineering, 2014, 136(6): 061015

12 Cadman J E, Zhou S, Chen Y, On design of multi-functional microstructural materials. Journal of Materials Science, 2013, 48(1): 51–66

13 Osanov M, Guest J K. Topology optimization for architected materials design. Annual Review of Materials Research, 2016, 46(1): 211–233

14 Sigmund O. Tailoring materials with prescribed elastic properties. Mechanics of Materials, 1995, 20(4): 351–368

15 Andreassen E, Lazarov B S, Sigmund O. Design of manufacturable 3D extremal elastic microstructure. Mechanics of Materials, 2014, 69(1): 1–10

16 Sigmund O, Torquato S. Composites with extremal thermal expansion coefficients. Applied Physics Letters, 1996, 69(21): 3203–3205

17 Guest J K, Prévost J H. Optimizing multifunctional materials: Design of microstructures for maximized stiffness and fluid permeability. International Journal of Solids and Structures, 2006, 43(22‒23): 7028–7047

18 Challis V J, Guest J K, Grotowski J F, Computationally generated cross-property bounds for stiffness and fluid permeability using topology optimization. International Journal of Solids and Structures, 2012, 49(23‒24): 3397–3408

19 Silva N E C, Fonseca O J S, Kikuchi N. Optimal design of piezoelectric microstructures. Computational Mechanics, 1997, 19(5): 397–410

20 Sigmund O, Søndergaard Jensen J. Systematic design of phononic band-gap materials and structures by topology optimization. Philosophical Transactions. Series A, Mathematical, Physical, and Engineering Sciences, 2003, 361(1806): 1001–1019

21 Rupp C J, Evgrafov A, Maute K, Design of phononic materials/structures for surface wave devices using topology optimization. Structural and Multidisciplinary Optimization, 2007, 34(2): 111–121

22 Sigmund O. Materials with prescribed constitutive parameters: An inverse homogenization problem. International Journal of Solids and Structures, 1994, 31(17): 2313–2329

23 Lotfi R, Ha S, Carstensen J V, Topology optimization for cellular material design. MRS Proceedings, 2014, 1662: mrsf13-1662-vv03-08

24 Carstensen J V, Guest J K. Topology optimization of nonlinear cellular materials. In: Proceedings of the 17th AIAA/ISSMO Multidisciplinary Analysis and Optimization Conference. Washington, D.C.: AIAA, 2016

25 Rodrigues H C, Guedes J M, Bendsøe M P. Hierarchical optimization of material and structure. Structural and Multidisciplinary Optimization, 2002, 24(1): 1–10

26 Sanchez-Palencia E. Non-Homogeneous Media and Vibration Theory. Berlin: Springer, 1980

27 Guedes J M, Kikuchi N. Preprocessing and postprocessing for materials based on the homogenization method with adaptive finite element methods. Computer Methods in Applied Mechanics and Engineering, 1990, 83(2): 143–198

28 Coelho P G, Cardoso J B, Fernandes P R, Parallel computing techniques applied to the simultaneous design of structure and material. Advances in Engineering Software, 2011, 42(5): 219–227

29 Nakshatrala P B, Tortorelli D A, Nakshatrala K B. Nonlinear structural design using multiscale topology optimization. Part I: Static formulation. Computer Methods in Applied Mechanics and Engineering, 2013, 261 ‒ 262: 167–176

30 Xia L, Breitkopf P. Concurrent topology optimization design of material and structure within FE² nonlinear multiscale analysis framework. Computer Methods in Applied Mechanics and Engineering, 2014, 278: 524–542

31 Sivapuram R, Dunning P D, Kim H A. Simultaneous material and structural optimization by multiscale topology optimization. Structural and Multidisciplinary Optimization, 2016, 54(5): 1267–1281

32 Xie Y M, Zuo Z H, Huang X, Convergence of topological patterns of optimal periodic structures under multiple scales. Structural and Multidisciplinary Optimization, 2012, 46(1): 41–50

33 Zuo Z H, Huang X, Yang X, Comparing optimal material microstructures with optimal periodic structures. Computational Materials Science, 2013, 69: 137–147

34 Zhang W, Sun S. Scale-related topology optimization of cellular materials and structures. International Journal for Numerical Methods in Engineering, 2006, 68(9): 993–1011

35 Schury F, Stingl M, Wein F. Efficient two-scale optimization of manufacturable graded structures. SIAM Journal on Scientific Computing, 2012, 34(6): B711–B733

36 Radman A, Huang X, Xie Y M. Topology optimization of functionally graded cellular materials. Journal of Materials Science, 2013, 48(4): 1503–1510

37 Alexandersen J, Lazarov B S. Topology optimization of manufacturable microstructural details without length scale separation using a spectral coarse basis preconditioner. Computer Methods in Applied Mechanics and Engineering, 2015, 290: 156–182

38 Lazarov B S. Topology optimization using multiscale finite element method for high-contrast media. In: Lirkov I, Margenov S, WaŚniewski J, eds. Large-Scale Scientific Computing. LSSC 2013. Lecture Notes in Computer Science, vol 8353. Berlin: Springer, 2014, 339–346

39 Osanov M, Carstensen J V, Tromme E, Topology optimization for additive manufacturing: New projection-based design algorithms. In: Proceedings of the 17th AIAA/ISSMO Multidisciplinary Analysis and Optimization Conference. Washington, D.C.: AIAA, 2016

40 Turner M J, Clough R W, Martin H C, Stiffness and deflection analysis of complex structures. Journal of the Aeronautical Sciences, 1956, 23(9): 805–823

41 Wang S, de Sturler E, Paulino G H. Large-scale topology optimization using preconditioned Krylov subspace methods with recycling. International Journal for Numerical Methods in Engineering, 2007, 69(12): 2441–2468

42 Amir O. Revisiting approximate reanalysis in topology optimization: On the advantages of recycled preconditioning in a minimum weight procedure. Structural and Multidisciplinary Optimization, 2015, 51(1): 41–57

43 Amir O, Stolpe M, Sigmund O. Efficient use of iterative solvers in nested topology optimization. Structural and Multidisciplinary Optimization, 2010, 42(1): 55–72

44 Aage N, Andreassen E, Lazarov B S. Topology optimization using PETSc: An easy-to-use, fully parallel, open source topology optimization framework. Structural and Multidisciplinary Optimization, 2015, 51(3): 565–572

45 Sigmund O, Petersson J. Numerical instabilities in topology optimization: A survey on procedures dealing with checkerboards, mesh-dependencies and local minima. Structural Optimization, 1998, 16(1): 68–75

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

47 Guest J K, Prévost J H, Belytschko T. Achieving minimum length scale in topology optimization using nodal design variables and projection functions. International Journal for Numerical Methods in Engineering, 2004, 61(2): 238–254

48 Guest J K, Smith Genut L C. Reducing dimensionality in topology optimization using adaptive design variable fields. International Journal for Numerical Methods in Engineering, 2010, 81(8): 1019–1045

49 Bruns T E, Tortorelli D A. Topology optimization of non-linear elastic structures and compliant mechanisms. Computer Methods in Applied Mechanics and Engineering, 2001, 190(26‒27): 3443–3459

50 Guest J K, Asadpoure A, Ha S H. Eliminating beta-continuation from heaviside projection and density filter algorithms. Structural and Multidisciplinary Optimization, 2011, 44(4): 443–453

51 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

52 Sigmund O. On the design of compliant mechanisms using topology optimization. Mechanics of Structures and Machines, 1997, 25(4): 493–524

53 Sigmund O, Aage N, Andreassen E. On the non-optimality of Michell structures. Structural and Multidisciplinary Optimization, 2016, 54(2): 361–373

54 Guest J K. Topology optimization with multiple phase projection. Computer Methods in Applied Mechanics and Engineering, 2009, 199(1‒4): 123–135

55 Clausen A, Aage N, Sigmund O. Exploiting additive manufacturing infill in topology optimization for improved buckling load. Engineering, 2016, 2(2): 250–257

56 Hashin Z. The elastic moduli of heterogeneous materials. Journal of Applied Mechanics, 1962, 29(1): 143–150

57 Hashin Z, Shtrikman S. A variational approach to the theory of the elastic behavior of multiphase materials. Journal of the Mechanics and Physics of Solids, 1963, 11(2): 127–140

58 Wu J, Aage N, Westermann R, Infill optimization for additive manufacturing—Approaching bone-like porous structures. IEEE Transactions on Visualization and Computer Graphics, 2018, 24(2): 1127–1140