Design of shell-infill structures by a multiscale level set topology optimization method

Computers and Structures 212 (2019) 162–172

Contents lists available at ScienceDirect

Computers and Structures journal homepage: www.elsevier.com/locate/compstruc

Design of shell-infill structures by a multiscale level set topology optimization method Junjian Fu, Hao Li, Liang Gao ⇑, Mi Xiao State Key Lab of Digital Manufacturing Equipment and Technology, Huazhong University of Science and Technology, Wuhan 430074, China

a r t i c l e

i n f o

Article history: Received 11 April 2018 Accepted 10 October 2018

Keywords: Shell-infill structures Level set method Numerical homogenization method Multiscale topology optimization

a b s t r a c t Shell-infill structures are widely used in additive manufacturing (AM) to retain the external appearance and reduce the printing costs. However, shell-infill structures are generally designed by a mono-scale topology optimization method. In this case, these shell-infill structures are not optimal designs due to their incompatible shell and infill layouts. In order to obtain the optimized shell and infill simultaneously, this paper presents a multiscale level set topology optimization method for designing shell-infill structures. In macroscale optimization, two distinct level sets of a single level set function (LSF) are used to represent the interface of the shell and the infill, respectively. The thickness of the shell is assumed to be uniform, which is guaranteed by a level set reinitialization scheme. In microscale optimization, the pattern of the microstructure is optimized to achieve the optimal macro structural performance. The numerical homogenization method is applied to evaluate the effective elasticity matrix of the microstructural infill. The shell-coated macro structure and the micro infill are optimized concurrently to achieve the optimal shell-infill design with prescribed volume fractions. Both 2D and 3D examples are investigated to demonstrate the effectiveness of the presented method. Ó 2018 Elsevier Ltd. All rights reserved.

1. Introduction A shell-infill structure is composed of a solid outer shell and a porous infill. The solid shell is the external wall to mainly preserve the outline shape of the structure. While the internal infill is generally featured as honeycomb-like cellular meso-structures to maintain the strength of the structure. The design of the shellinfill structure is inspired by nature, such as plant stems, bird beaks and human bones, where cellular meso-structures coexist with solid shells. These bio-structures are evolved to withstand external loads and keep their functions in a dynamic balance for millions of years. Shell-infill structures have special applications in coated implants, architectures and additive manufacturing (AM). For example, in AM technologies, e.g. Fused Deposition Modeling (FDM), printing solid parts requires large amounts of material and long printing time, leading to high costs. Therefore, a great number of additive manufactured or 3D printed parts are not printed as solid. Printing parts with a solid shell and a porous infill is a compromise between the cost and the strength [1]. However, designing of shell-infill structures by experience and repeated simulation is inefficiency. And the optimal structural performance is

⇑ Corresponding author. E-mail address: [email protected] (L. Gao). https://doi.org/10.1016/j.compstruc.2018.10.006 0045-7949/Ó 2018 Elsevier Ltd. All rights reserved.

generally not guaranteed. In this situation, it is still a challenging issue in designing shell-infill structures. Topology optimization [2] is considered as a powerful design method for weight reduction, optimal structural design and multi-functional design in AM. A plenty of topology optimization methods, including the homogenization method [2,3], the density method [4,5], the level set method (LSM) [6,7], the evolutionary structural optimization method (ESO) [8,9] and the moving morphable components (MMC) method [10,11], have developed their own ways for AM. Topology optimization of overhang-constraint structures [12–15], multi-material structures [16,17], microstructures [18,19], multiscale structures [20] and shell-infill structures [21,22] are tailored specifically for AM. Some recent work in regards to the design of shell-infill structures based on topology optimization have been proposed by Clausen originally [21,22]. In the method, the shell-infill or coating structures are designed based on Solid Isotropic Material with Penalization (SIMP) scheme [4]. By introducing multiple smoothing and projection strategies, the material interpolation identifies the interface or the shell of structures by density spatial gradients. The thickness of the interface can be highly uniform. Experiments show that the shell-infill structure exhibits increased performance under buckling loads [23]. This design method is extended to design AM oriented hierarchical shell-infill structures [24]. Apart

J. Fu et al. / Computers and Structures 212 (2019) 162–172

from the density method, level set based topology optimization method can also be used to design shell-infill structures. As a viable alternative, the LSM handles the boundary problems easily by its implicit nature. A level set topology optimization method considering interface effects has been presented earlier [25]. Wang and Kang [26] realize that the LSM is a good choice for shell-infill or coated structures. Only one level set function (LSF) is employed in the description of both the shell and infill, which facilitates the structural representation and sensitivity derivation. And the uniform thickness of the shell is naturally guaranteed by a level set reinitialization scheme. Nevertheless, the majority of current studies are focused on the single scale design in which both the shell and the infill are treated as separated materials [21,26]. Although the infill can be designed as predefined microstructures, the single scale topology optimization of the shell-infill structure limits the structural performance, keeping the design away from its real optimum. It should be noted that even though Wu’s method [24] generates bone-like infill structures with a solid shell, it is strictly a macroscale design with multiple sub-area constraints. A shell-infill structure is a typical structure that should be designed by a multiscale method. Because the upscale shell with a thickness has direct influence on the downscale microstructural infill and vice versa. In order to fully investigate the characteristics and the performance of shell-infill structures, it is desirable to extend the mono-scale design method to a multiscale topology optimization method. The multiscale topology optimization is a concurrent design method that optimizes the macrostructure and material microstructures simultaneously with a prescribed objective function. The numerical homogenization method is applied to bridge the microstructures and macroscale material properties [27]. The first attempt for the multiscale topology optimization of space varying microstructures is conducted by Rodrigues [28] in 2002. After that, Liu et al. [29] propose a multiscale topology optimization method by assuming the uniformity of material microstructures in macroscale. Multiscale topology optimization has gained its prosperity with the development of AM technologies. Many multiscale topology optimization methods and strategies have been proposed by researchers in recent years. As classified by the method, there are SIMP based [28–33], LSM based [34–38] and ESO based [39–41] multiscale topology optimization method. But different methods generally result in similar optimized designs. As for the optimization strategies, there are multiscale optimization strategies with uniform [29,34,39], element-wise [28,40,41], layer-wise [30,36], graded [38,42] and multi-patch microstructures [35,37]. Generally speaking, these strategies are proposed as compromises between the computational cost and the structural performance. As for the applications, there are compliance design [33,37], frequency design [43,44] and multi-materials design [45]. The multiscale structures optimized by these methods are composed of complex macrostructures and force-adapted microstructures. This provides us an insight into the design of shell-infill structures by a multiscale method. However, the above multiscale optimized results do not consider the shell effects, i.e. no shell is incorporated into the macrostructure. Consequently, the multiscale design method cannot be applied directly to the shell-infill structural design. It is necessary to extend the current multiscale topology optimization method to a new one that considers the shell and infill simultaneously. This paper presents a multiscale topology optimization method for shell-infill structures. The shell-infill structural representation is based on the LSF. By specifying parameters like shell thickness and infill density (volume fraction), optimal shell layout and infill pattern are optimized concurrently. Both the shell and the infill are optimized with the same macro objective function. The infill is assumed to be uniform microstructures distributed over the

163

whole design domain. Manufacturing constraints, like overhang constraints [11–14], can be implemented independently on these uniform microstructures. The remaining of this paper is organized as follows: Section 2 deals with the parametric level set method (PLSM). Section 3 introduces the geometric and physical representation of the shell-infill. Section 4 develops the mathematical model and performs the sensitivity analysis for the multiscale shell-infill topology optimization. Section 5 investigates some 2D and 3D examples to demonstrate the effectiveness of the presented method. Section 6 concludes the paper. 2. Level set method and its parametrization LSM is initially proposed in 1988 to track a moving boundary [6]. Rather than following the boundary itself, LSM builds the boundary into one higher dimensional surface as the zero contour. Instead of using the conventional LSM based topology optimization [46,47], this paper focuses on a further developed PLSM due to its simplicity of implementation [48,49]. 2.1. Level set method In the framework of LSM, the structure X is described by the following definition

8 > < Uðx; tÞ > 0 x 2 X Uðx; tÞ ¼ 0 x 2 @ X > : Uðx; tÞ < 0 x 2 D=ðX [ @ XÞ

ð1Þ

where X describes the structure. @ X is the interface of the structure represented implicitly as the zero contour of LSF. D is the design domain that contains the structure completely. The pseudo time t adds dynamics to the implicit interface @ X. Taking the derivative of Uðx; tÞ ¼ 0 on both sides with respect to t and giving the unit outward normal n ¼ rU=jrUj on the interface, we obtain the level set equation as a Hamilton-Jacobi equation

@ Uðx; tÞ  V n jrUj ¼ 0 @t

ð2Þ

where V n is the normal velocity field. The signed distance property of the LSF is preserved by solving the reinitialization equation [46,47].

@U ¼ SðUÞð1  jrUjÞ @t

ð3Þ

with

U

SðUÞ ¼ qffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffi U2 þ n2 jrUj2

ð4Þ

where SðUÞ is the sign function. n is a small number specified as the mesh size. 2.2. Parametrization of LSM In this study, the authors prefer a PLSM that can be solved by gradient based mathematical optimizers. The scalar LSF is parametrized by a CSRBF based interpolation.

Uðx; tÞ ¼

Ne X

ui ðxÞai ðtÞ ¼ uT ðxÞaðtÞ

ð5Þ

i¼1

where N e is the number of knot. The vector of spatial shape functions uðxÞ and the vector of time dependent expansion coefficients aðtÞ are given as

164

J. Fu et al. / Computers and Structures 212 (2019) 162–172

uðxÞ ¼ ½ u1 ðxÞ u2 ðxÞ    uN ðxÞ T

ð6Þ

aðtÞ ¼ ½ a1 ðtÞ a2 ðtÞ    aN ðtÞ T

ð7Þ

where ui ðxÞ is the Wendland’s CSRBF with C2 continuity [50]

ui ðxÞ ¼ ðmaxð0; 1  rÞÞ4  ð4r þ 1Þ rðkx  xi kÞ ¼

qffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffi ðx  xi Þ2 þ ðy  yi Þ2 dmI

ð8Þ

ð9Þ

x is the coordinate matrix of the sample knot. dmI denotes the radius of influence domain for a specific knot ðxi ; yi Þ. Substituting Eq. (5) into Eq. (2), the parametrized level set equation is written by

u ðxÞa_ ðtÞ  V n jruðxÞ aðtÞj ¼ 0 T

T

ð10Þ

With this parametrization, the original time dependent Hamilton-Jacobi equation is decoupled to a system of ordinary differential equations. Gradient based mathematical optimizers, such as optimality criteria method (OC) [4] and the method of moving asymptotes (MMA) [51], can be used to update the design variables. It should be noted that the reinitialization scheme is not a necessity for the PLSM based topology optimization method [48]. The introducing of reinitialization into this paper is especially for the design of shell-infill structures to preserve a uniform thickness of the shell. Global reinitialization of the LSF may prevent emerging of new holes during the evolution, the authors simply initialize the design domain with reasonably sufficient holes to include the admissible topology configurations. Although incorporating topological derivatives into LSM is a solution for the hole nucleation [52], this is beyond the research scope of this paper.

jrUj ¼ 1. The 0-k level set area in Fig. 1(a) is contoured straightforwardly with uniform thickness k. Contrarily, the un-reinitialized LSF in Fig. 1(b) has a too steep slope. For one thing, the value of the thickness of the shell is unknown and out of control. For another thing, the thickness of the shell is not uniform enough. According to the definition of shell-infill in (11), the material property C inside a point of the macro structural domain is defined as follows

C ¼ CS ðHðUMA Þ  HðUMA  kÞÞ þ CI ðuMI ; UMI ÞHðUMA  kÞ

where CS is the elasticity property of the shell, generally prescribed as the property of an isotropic solid material. CI is the elasticity property of the infill, which is dependent on the microstructural LSF UMI and fluctuation displacement uMI . HðÞ is the Heaviside function. k denotes the thickness of the shell. 3.2. The homogenization method It is assumed that the infill domain is composed of periodically distributed microstructures that have a much smaller size than that of the macro structure. The homogenization method [3] is applied to predict the effective elasticity matrix CI of periodic microstructures in the infill domain. Under the framework of PLSM, the effective elasticity property is calculated as MI MI C ijkl I ðu ; U Þ ¼

Z

1 jDMI j

T

DMI

0ðijÞ 0ðklÞ ðepq  epq ðuMIðijÞ ÞÞ Cpqrs ðers

 ers ðuMIðklÞ ÞÞHðUMI ÞdXMI

ð13Þ

where i; j; k; l ¼ 1; 2: Cpqrs is the elasticity property of solid material. jDMI j is the volume of the domain of the unit cell. epq denotes the applied test strain field, including unit strain in the horizontal direc0ðijÞ

tion ½1; 0; 0T , unit strain in the vertical direction ½0; 1; 0T and unit

3. Geometric and physical representation of the shell-infill

shear strain ½0; 0; 1T (for 2D problems) [27]. The periodic fluctuating displacement uMI can be obtained by solving the following micro equilibrium equation

3.1. Representation of the shell-infill Two distinct level sets of a single LSF are utilized to describe both the structural domain of the shell and the infill as shown in Fig. 1. In this research, a shell with uniform thickness is considered. The external interface of the shell is referred as CS . This interface is easily captured by UMA ðx; tÞ ¼ 0, known as the zero level set of the LSF. The uniform shell can be obtained by offsetting the external interface CS to an internal interface CI . The internal interface CI is defined by U ðx; tÞ ¼ k, the k level set of the LSF. The distance k between CS and CI is equivalent to the thickness of the shell. The area surrounded by the two interfaces is the shell part. The area inside the interface CI is the infill part. In the framework of a LSF representation, a shell-infill structure with a uniform shell thickness is defined as follows.

aMI ðuMI ; v MI ; UMI Þ ¼ l ðv MI ; UMI Þ 8v MI 2 U MI MI

with

aMI ðuMI ; v MI ; UMI Þ ¼

x 2 XI ðinfillÞ x 2 CI ðinterface between infill and shellÞ x 2 XS ðshellÞ x 2 CS ðinterface between shell and voidÞ x 2 XV ðvoidÞ ð11Þ

The LSF must be maintained as a signed distance function to achieve a uniform shell with the thickness under control. Fig. 1 are two LSFs with and without reinitialization, respectively. The reinitialized LSF has uniform gradients on the whole domain, i.e.

Z DMI

ð14Þ

ðepq ðuMIðijÞ ÞÞ Cpqrs ers ðv MIðklÞ ÞHðUMI ÞdXMI T

ð15Þ

MA

8 MA > U ðx; tÞ > k > > > > MA > > < U ðx; tÞ ¼ k 0 < UMA ðx; tÞ < k > > > > UMA ðx; tÞ ¼ 0 > > > : MA U ðx; tÞ < 0

ð12Þ

l ðv MI ; UMI Þ ¼ MI

Z DMI

0ðijÞ ðepq Þ Cpqrs ers ðv MIðklÞ ÞHðUMI ÞdXMI T

ð16Þ

where v MI is the virtual displacement field in microscale, belonging to the space U MI which includes all the kinematically admissible displacements. 4. Multiscale topology optimization of shell-infill structures The shell-infill structure is optimized in an integrated manner. The compliance of the macro structure is dependent not only on the macroscopic layout but also on the microscopic pattern of the microstructure. 4.1. Multiscale shell-infill design model The compliance minimization of a shell-infill structure is defined as following

165

J. Fu et al. / Computers and Structures 212 (2019) 162–172

a) The LSF intersected by Φ MA = 0 and Φ MA = k with reinitialization and the 0-k contour

b) The LSF intersected by Φ MA = 0 and Φ MA = k without reinitialization and the 0-k contour Fig. 1. Illustration of the shell-infill structure represented by a LSF.

Z 1 Min:JðuMA ;uMI ;UMA ;UMI Þ¼ ðeT ðuMA ÞC ðUMA ;UMI ;uMI ÞeðuMA ÞÞdXMA 2 DMA 8 MA MA MI G ðU ;U Þ¼ > > R   > > > R HðUMI ÞdXMI MA > MA MA MA DMI > ðHð U ÞHð U kÞÞþ Hð U kÞ dXMA f V MA 60 MA > MI > D V > > > > MI MI R > a ðu ;v ;U Þ¼l ðv MA ;UMA Þ 8v MA 2U MA > > > > >aMI ðuMI ;v MI ;UMI Þ¼lMI ðv MI ;UMI Þ 8v MI 2U MI > > > > >aMA 6aMA 6aMA > > max > : min MI MI aMI min 6 a 6 amax ð17Þ

where the superscript ‘MA’ and ‘MI’ indicate quantities belong to macroscale and microscale, respectively. uMA is the macro displacement filed. uMI is the micro displacement filed. UMA and UMI are the LSF of the macro and micro structure respectively. e is the strain field of macro structure. V MA and V MI are the macro and micro total MA

MI

volume of their design domain. f and f are prescribed macro and micro volume fraction, respectively. aMA and aMI are the expansion coefficients for CSRBF-based interpolation in macroscale and microscale, respectively. The bilinear and linear terms of the macroscale equilibrium equation are given as

aMA ðuMA ; v MA ; UMA Þ ¼

Z DMA

ðeT ðuMA ÞC ðUMA ; UMI ; uMI Þeðv MA ÞÞdXMA ð18Þ

l

MA

ðv MA ; UMA Þ ¼

Z

MA

fv DMA

Z

þ where

DMA

HðUMA ÞdXMA

pvMA dðUMA ÞjrUMA jdXMA

ð19Þ

v MA is the virtual macro displacement field belonging to the

space U MA spanned by the kinematically admissible set of macro displacements. f is the body force density of macro structure and

p is the traction force on macro boundary. dðÞ is the partial derivative of the Heaviside function. The smoothed Heaviside function HðUÞ is defined as

HðUÞ ¼

8 > < f;

U < D

3ð1fÞ U ð > 4 D

:

3

 3UD3 Þ þ 1þf ; D 6 U 6 D 2

ð20Þ

U>D

1:

where f ¼ 0:001 is a small positive number. D is half of the bandwidth. The derivative of the Heaviside function is known as Dirac function, which is derived as

( dðUÞ ¼

3ð1fÞ ð1 4D

0:

2

 UD2 Þ; jUj 6 D jUj > D

ð21Þ

UMA and UMI are two LSFs in Eq. (17) serving as the macro structure and the microstructure, respectively. The layout of the shell is defined on the macro LSF UMA . The microstructural configuration and effective elasticity property are determined by the micro LSF

UMI . The design variables of the mathematical model can be either the macro coefficients aMA or the micro coefficients aMI . When it comes to aMA , it is the upscale optimization of the macro structure with a shell. While design variables aMI are for the downscale optimization of the microstructures of the infill. These two optimization processes are realized in a concurrent way. The main difference between the multiscale model for the shellinfill structure and a general structure is that a solid shell is introduced to the macroscale structure. The variable k in (17) denotes the level set value that controls the thickness of the shell. The definition of geometry and material properties for the shell-infill are embedded in the objective function. Accordingly, the macro volume constraint is newly constructed to count both the shell and the infill. The micro volume fraction directly influences the macro structural design, which will be investigated in detail in Section 5.

166

J. Fu et al. / Computers and Structures 212 (2019) 162–172

4.2. Sensitivity analysis in macroscale As the objective function is related to both scales, the derivative with respect to macro and micro variables should be derived separately. The material derivative and the adjoint variable method are used to derive the shape sensitivity of the objective and the constraint. The material derivative of the macro objective is given by _ MA ; UMA Þ ¼ Jðu

Z

MA 

e ðu Þ C ðU ; U ;u Þeðu ÞdX T

DMA

Z



DMA

Z 



MI

MI

MA

MA

eT ðruMA  V MA ÞC ðUMA ; UMI ; uMI ÞeðuMA ÞdXMA þ

1 2

MA

MI

D

@ UMA 1  MA dXMA þ 2 @t

Z D

MI

MA

MA

eT ðuMA ÞCS eðuMA ÞdðUMA Þ MA

ð22Þ

0

where ðÞ or ðÞ denotes the material derivative. ðÞ ¼ @ðÞ=@t

where

dðUMA  kÞ dðU

MA

Þ

Þ þ CI ðuMI ; UMI Þ

dðUMA  kÞ dðUMA Þ ð24Þ

The shell-infill topology optimization problem considered in this paper is a special case of the multi-materials problem. The stress/strain tensor along the interface between the shell and infill would be discontinuous due to the distinct elasticity properties. It  should be noted that ðuMA Þ belongs to the admissible set and will be handled by the adjoint method. The rigorous derivation for the level set based topology optimization involving multi-materials is given by Guo et al. [53]. Comparing with the single material problem, there is a jump term in the shape derivative of multi-materials problem. When level set description with extended finite element model is used, the jump term must be calculated exactly. However, this term can be smeared out since the ersatz material is adopted for finite element analysis (FEA) in this study [54]. Considering that the compliance minimization problem is self-adjoint, we obtain the shape derivative by following a similar way to [53]

_ MA ; UMA Þ ¼  Jðu

Z

e ðu ÞCU ðU ; U ; u Þeðu ÞdðU T

DMA

MA



MA

MI

MI

MA

MA

@ UMA Þ MA dX @t ð25Þ

Parametrizing the macro LSF by Eq. (10) and substituting into Eq. (25) N Z X _ MA ; UMA Þ ¼  Jðu i

DMA

MA eT ðuMA ÞCU ðUMA ; UMI ;uMI ÞeðuMA ÞdðUMA ÞuMA Þ i ðx

daMA i dt

MA

dXMA ð26Þ

where N is the number of knot in macroscale. Applying the chain rule for the macro objective in (17)

@Jðu ; U @tMA MA

MA

N Þ X @JðuMA ; UMA Þ daMA i ¼ MA @ aMA dt i i

DMA

MA ðeT ðuMA ÞCU ðUMA ; UMI ; uMI ÞeðuMA ÞÞdðUMA ÞuMA ÞdXMA i ðx

ð28Þ

Similarly, based on the material derivative and chain rule, the shape derivative of macro volume constraint with respect to macro CSRBF expansion coefficient aMA is given by i

Z DMA

MA LðUMA ; UMI ÞdðUMA ÞuMA ÞdXMA i ðx

dðUMA  kÞ

LðUMA ; UMI Þ ¼ 1 

dðU

MA

MA

represents the macro spatial derivative. V MA is the macro velocity filed. The above formula is derived under the assumption that the thickness k of the shell remains constant during the optimization process. Formula (22) is further simplified as Z  _ MA ;UMA Þ¼ eT ðuMA Þ C ðUMA ;UMI ;uMI ÞeðuMA ÞdXMA Jðu DMA Z 1  eT ðruMA V MA ÞC ðUMA ;UMI ;uMI ÞeðuMA ÞdXMA þ 2 DMA Z MA @ U  eT ðuMA ÞCU ðUMA ;UMI ;uMI ÞeðuMA ÞdðUMA Þ MA dXMA ð23Þ @t DMA

CU ðUMA ; UMI ; uMI Þ ¼ CS ð1 

Z

ð29Þ

where

@ UMA MA dX @tMA



@JðuMA ; UMA Þ ¼ @ aMA i

@GMA ðUMA ; UMI Þ ¼ @ aMA i

e ðu ÞðCI ðu ; U Þ  CS Þeðu ÞdðU  kÞ MA T



MA

Comparing equation (26) with (27), we obtain the shape derivative of the macro objective with respect to the macro CSRBF expansion coefficient aMA i

ð27Þ

Þ

R þ

DMI

HðUMI ÞdXMI dðUMA  kÞ V MI

dðUMA Þ

ð30Þ

4.3. Sensitivity analysis in microscale Noting that the macro LSF UMA and the shell thickness k are not dependent on the micro variables, the micro shape derivative of the objective is expressed by following a similar derivation to [37]

_ MI ; UMI Þ ¼  Jðu

Z MA

D

ðeT ðuMA Þ

@CI ðuMI ; UMI Þ HðUMA @tMI

 kÞeðuMA ÞÞdXMA

ð31Þ

The derivative of effective elasticity matrix CI ðuMI ; UMI Þ with respect to t MI is given by

@C I ðuMI ; UMI Þ 1 ¼  MI @t MI jD j

Z MI

bðuMI ÞdðUMI Þ

D

@ UMI dXMI @tMI

ð32Þ

where

bðuMI Þ ¼



T



 MIðijÞ 0ðklÞ e0ðijÞ Þ Cpqrs ers  ers ðuMIðklÞ Þ pq  epq ðu



ð33Þ

By parametrizing the micro LSF, Eq. (32) is re-written as

Z M X daMI @C I ðuMI ; UMI Þ 1 j MI ¼ bðuMI ÞdðUMI ÞuMI dXMI j ðx Þ MI MI MI MI @t j jD D dt j ð34Þ where M is the number of knot in microscale. Applying the chain rule for @C I ðuMI ; UMI Þ=@t MI MI M @C I ðuMI ; UMI Þ X @C I ðuMI ; UMI Þ daj ¼ MI MI MI @ aj @t dt j

ð35Þ

Comparing (34) with (35), we obtain the shape derivative of the effective elasticity matrix with respect to the micro CSRBF expansion coefficient aMI j

@C I ðuMI ; UMI Þ 1 ¼  MI @ aMI jD j j

Z DMI

MI MI bðuMI ÞdðUMI ÞuMI j ðx ÞdX

ð36Þ

Applying the chain rule for @JðuMI ; UMI Þ=@tMI and substituting (36) into (31), we obtain the shape derivative of the objective with respect to micro CSRBF expansion coefficient aMI j

@JðuMI ; UMI Þ ¼ @ aMI j

Z MA

D

ðeT ðuMA Þ

@CI ðuMI ; UMI Þ HðUMA  kÞeðuMA ÞÞdXMA @ aMI j ð37Þ

167

J. Fu et al. / Computers and Structures 212 (2019) 162–172

The micro volume constraint with respect to micro CSRBF is similarly derived as expansion coefficient aMI j

@GMI ðUMI Þ ¼ @ aMI j

Z DMI

MI MI dðUMI ÞuMI j ðx ÞdX

ð38Þ

5. Numerical examples In the multiscale design for shell-infill structures, the material property of the shell is not dependent on the optimized microstructure. It is the property of the micro infill that is determined by the optimized microstructure. Details of the monoscale design of 2D and 3D microstructures [55,56] are not provided in this paper. Essentially, the material property of the infill could be different from element to element. The element-wise optimization strategy may lead to a great computational burden and difficulties in manufacturing. As a result, in order to facilitate the manufacturing and optimization process, uniform microstructures are assumed on the whole infill domain. The following numerical examples are investigated to demonstrate the effectiveness of the proposed method. This section provides multiscale shell-infill designs for some popular structures. The authors firstly compare the shell-infill designs with non-shell designs. The influences of the infill volume fraction on the shellinfill structures are also discussed. The method is finally extended to some other 2D and 3D examples. As defined previously, all the examples are limited to compliance minimization. The mathematical optimizer of the multiscale optimization is OC. The optimization terminates when the absolute difference of two successive objectives is less than 1e5 or the maximum loop number is reached. It is also assumed that the shell is made of a solid material with elasticity modulus Es = 1 and Poisson’s ratio vs = 0.3, respectively. While the infill is composed of porous microstructures whose base material property is the same as the solid material used in shell. Unless otherwise stated, the thickness of the shell is prescribed as one elemental size, i.e. k = 1. Because the shell is designed on the macro structure, the level set reinitialization scheme is imposed only on macro LSF to guarantee a shell with uniform thickness. 5.1. The influence of the shell on the infill structures The first example is presented to investigate the influence of the shell on the multiscale infill design. The macro base structure is a beam with aspect ratio L:H = 1:2 (L for X direction, H for Y direction). The left boundary of the design domain is fixed. Two load cases indicated in Fig. 2 are applied in this example. For case A, a concentrated load F = 1 is applied on the middle point of the right boundary. For case B, the concentrated load is applied on the lower corner of the right side. The macro design domain of the cantilever beam is discretized by 20  40 plane stress quadrilateral elements. The micro design domain is a square discretized by 40  40 plane stress quadrilateral elements. The macro and micro volume fraction in Eq. (17) are prescribed as fMA = 1 and fMI = 0.5 respectively. For comparison of the shell effects, only microstructures of the infill are optimized. The macro structures are not involved in the multiscale optimization. But they are simply coated without/with a shell. The thickness of the shell is set to one unit elemental size. This is a frequently encountered problem in AM, in which the print model is generally pre-determined by engineers, except for the shell and infill parameters are left to be specified. Table 1 and Fig. 3 show that both the topology and the effective elasticity matrix of the optimized microstructures have changed after the incorporation of a shell. The number of microstructures

within the macrostructure is interpolated from 20  40 to 8  16 for better visualization. For case A, both of the optimized infill microstructures exhibit orthogonal isotropy. A uniform coating of the macro structure brings relatively more strong materials in the Y directional boundaries and less in the X directional boundaries since the X-Y ratio of the macro structure is 1:2. After the incorporating of shell, there is a significant decrease of the Y directional bulk modulus CI ð2; 2Þ, e.g. decreased from 0.2804 (without a shell) to 0.2060 (with a shell) in case A. This phenomenon is observed because there are too much strong shell materials have already been distributed in Y direction, the optimized microstructures have paid less strength on CI ð2; 2Þ. For case B, the optimized infill microstructures exhibit anisotropy since the concentrated load is applied on the lower corner of the design domain. The load transmission in case B is more complex. Still, a relative significant decrease of CI ð2; 2Þ is observed in Fig. 3. The authors only present two cases in this example to demonstrate the effects by adding a shell to the conventional multiscale design. The shell dose have some effects on the microstructural design. However, it is difficult to find a regular rule for the variations of the microstructural topology and effective elasticity matrix after considering the shell. Thus, generating a shell by a post processing method based on the traditional multiscale design without a shell cannot guarantee the optimal microstructures. It is desirable to incorporate the shell into the mathematical model of the multiscale topology optimization such that the optimal macro shell and the micro infill are obtained simultaneously. 5.2. The influence of the shell on macro structures Not only does the shell affects the optimized microstructures of the infill, it also has some influences on the macro structure. The macro structure is a cantilever beam with aspect ratio L:H = 2:1. The left boundary of the design domain is fixed. As illustrated in Fig. 4, a concentrated load F = 1 is applied on the lower right corner. The macro design domain of the cantilever beam is discretized by 80  40 plane stress quadrilateral elements. The micro design domain of the microstructure is a square discretized by 40  40 plane stress quadrilateral elements. The macro and micro volume constraint are prescribed as fMA = 0.4 and fMI = 0.7, respectively. The number of microstructures within the macrostructure is interpolated from 80  40 to 16  8 for better visualization. Table 2 shows that both the macro structure and microstructure are changed compared with the counterparts without a shell. In the

a) Load case A

b) Load case B

Fig. 2. The design domain of two load cases.

168

J. Fu et al. / Computers and Structures 212 (2019) 162–172

Table 1 The optimized shell-infill structure and effective elasticity matrix. Load cases

Macrostructure

Microstructure

Effective elasticity matrix 2 3 0:2201 0:1140 0:0000 4 0:1140 0:2804 0:0000 5 0:0000 0:0000 0:1018

Case A

2

0:2361 4 0:1481 0:0000

2

Case B

0:2206 4 0:0864 0:0222

2

0:2117 4 0:1273 0:0393

a) Load case A

0:1481 0:2060 0:0000

3 0:0000 0:0000 5 0:1353

0:0864 0:3171 0:0269

3 0:0222 0:0269 5 0:0684

0:1273 0:2067 0:0521

3 0:0393 0:0521 5 0:1207

b) Load case B

Fig. 3. Comparison of the effective elasticity matrices in bar graphs.

macroscopic design, the structure with a shell is thinner than the structure without a shell. Because giving a total volume, the shell consumes more materials on a unit area than the infill does. The shell-coated structure indicates materials are distributed more efficiently on the force transmission path, leading to a stiffer structure (Fig. 5). In microscopic design, one may find that the changes on the topology and effective elasticity matrix are not significant enough compared with the example in Section 5.1, because the shell is coated along the boundary without a directional (X or Y) pattern. It should be noted that the multiscale optimization with uniform microstructures in this study is actually a multi-objective multiscale optimization for all elements with identical weights, which means the changes on the topology as well as the effective elasticity matrix are averaged by all elements.

5.3. The influence of the micro volume fraction on shell-infill structures This numerical example is designed to investigate the influences of micro volume fraction of the infill on the shell-infill structures. The design domain, as illustrated in Fig. 6, is discretized by 90  30 quadrilateral elements with aspect ratio L:H = 3:1. The X and Y displacements of four corners are fixed. A concentrated load F = 1 is applied on the center of the structure. The macro structural volume fraction is given as fMA = 0.4 in this example. Table 3 shows the shell-infill designs with fixed macro volume fraction but varying micro volume fractions. The number of microstructures within the macrostructure is interpolated from 90  30 to 18  6 for better visualization. Volume fractions of the microstructures of the infill are given as 0.4, 0.5, 0.6, 0.7, 0.8, 0.9 and 1.0, respectively. The solid shells with thickness of one

169

J. Fu et al. / Computers and Structures 212 (2019) 162–172

Fig. 4. Design domain of a cantilever beam.

elemental size are well preserved in all the macro designs. The area of the macro structure decreases as the increase of the microstructural volume fraction, because the density of the infill is becoming larger when the total macro volume fraction is a fixed value. A small micro volume fraction fMI signifies a highly porous infill, leading to a large-area macro structure. On the contrary, a large micro volume fraction fMI contributes to a small-area macro structure. It will be seen that the design approaches a pure macro topology optimization with fMI = 1. Fig. 7 shows the compliance of the shell-infill design with varying microstructural volume fractions. It demonstrates that the larger the micro volume fraction, the smaller the structural compliance. One may prefer a stiff solid design by prescribing a larger micro volume fraction in microscale. Although the shell-infill structure is not as stiff as the solid design with the same volume fraction, it is reported to have better stability with buckling loadings by slightly sacrificing its stiffness [23]. Based on the volume definition of the multiscale optimization formulation in (17), the macro volume fraction fMA has a maximum value once the shell thickness k and the micro volume fraction fMI are specified. When the shell thickness equals to zero, i.e. k = 0, the maximum macro volume fraction fMA equals to the micro volume fraction fMI. In this case, the microstructures are periodically distributed over the whole macro design domain. When the shell thickness is larger than zero, i.e. k > 0, the maximum macro volume fraction fMA can be slightly larger than the micro volume fraction fMI. Because the shell is assumed to be composed of solid materials, and thus the shell is always denser than the infill. Compared with the case k = 0, more materials are needed to fill the design domain with microstructures and an outer shell.

5.4. Shell-infill design of a MBB beam This example investigates the shell-infill design for the popular MBB beam. The design domain with aspect ratio L:H = 6:1 is

Fig. 5. Iteration history of multiscale shell-infill structures.

Fig. 6. Design domain of four-node fixed structure.

discretized by 180  30 quadrilateral elements. As shown in Fig. 8, the left lower corner is fixed and the right lower corner is simply supported. A concentrated load F = 1 is applied at the center of the upper boundary. The macro and micro volume fraction are given as fMA = 0.4 and fMI = 0.6, respectively. The MBB beam is designed in its full length not a half-length. The shell thickness is highly uniform since a reinitialization scheme is imposed on the macro LSF. Both of the optimized macrostructure with a shell and the microstructure are given in Fig. 9. The number of microstructures within the macrostructure is interpolated from 180  30 to 36  6 for better visualization. As indicated in Fig. 10, the micro volume constraint is satisfied after several iterations. However, the topology of the microstructure may not converge, which will have directly on the convergence of the macro objective and constraint. The compliance of the shell-infill MBB converges to 146.47 after 92 iterations. And the whole optimization process is very stable.

Table 2 The multiscale shell-infill design. Thickness of shell k=1

k=0

Macrostructure

Microstructure

Effective elasticity 2 0:5944 0:1549 4 0:1549 0:2681 0:0217 0:0197

2

0:6291 4 0:1277 0:0239

0:1277 0:2932 0:0237

matrix 3 0:0217 0:0197 5 0:1704

3 0:0239 0:0237 5 0:1455

170

J. Fu et al. / Computers and Structures 212 (2019) 162–172

Table 3 Shell-infill design with varying fMI. fMA

fMI

0.4

0.4

Macro design with the shell

0.5

0.6

0.7

0.8

0.9

1.0

Micro infill design

Effective elasticity matrix 2 3 0:1239 0:0964 0:0000 4 0:0964 0:1972 0:0000 5 0:0000 0:0000 0:0825 2

0:2476 4 0:1201 0:0000 2

0:3545 4 0:1772 0:0000 2

0:4971 4 0:1997 0:0000 2

0:6559 4 0:2337 0:0000 2

0:8574 4 0:2694 0:0000 2

1:0989 4 0:3297 0:0000

0:1201 0:2453 0:0000

3 0:0000 0:0000 5 0:1094

0:1772 0:2778 0:0000

3 0:0000 0:0000 5 0:1657

0:1997 0:3806 0:0000

3 0:0000 0:0000 5 0:1998

0:2337 0:5413 0:0000

3 0:0000 0:0000 5 0:2456

0:2694 0:7927 0:0000

3 0:0000 0:0000 5 0:3037

0:3297 1:0989 0:0000

3 0:0000 0:0000 5 0:3846

a) The shell-infill structure

0.5450 0.0672 0.0000 CI = 0.0672 0.1222 0.0000 0.0000 0.0000 0.0913 b) The microstructure and elasticity matrix Fig. 9. Shell-infill design of the MBB beam.

Fig. 7. The compliance of the shell-infill design with varying fMI.

Fig. 8. The design domain of the MBB beam.

Fig. 10. Iteration history of the MBB design.

J. Fu et al. / Computers and Structures 212 (2019) 162–172

171

Fig. 13. Optimized shell-infill design.

Fig. 11. The macro design domain of a 3D cantilever beam.

for some AM techniques. A closed shell can be achieved by adding passive variables constraint on the boundary nodes of the LSF. The microstructure infill is analyzed and optimized based on the homogenization theory. There is a micro FEA in each optimization iteration. For the discrete equilibrium equation of microstructural FEA, it is a linear system with multiple right-hand sides. Solving this system of linear equations for 3D problems are computationally challenging with the refinement of mesh, because there are six right-hand sides for 3D problems. To accelerate the solving, a block conjugate gradient method [57,58] is suggested for the 3D microstructural FEA. 6. Conclusions

Fig. 12. The shell (red) and infill (blue). (For interpretation of the references to color in this figure legend, the reader is referred to the web version of this article.)

5.5. Extensions to a 3D problem A 3D design case is utilized to further demonstrate the generality of the proposed method. As shown in Fig. 11, the macro design domain is a hexahedron with length-width-height ratio 2:1:2. All the degree of freedoms on the four corners of the left face are fixed. A concentrated load F = 1 is applied at the middle point of right face of the structure. The macro hexahedral domain and the micro cubic domain are discretized by 20  10  20 and 14  14  14 eight-node hexahedral elements, respectively. The macro structural volume fraction and the microstructural volume fraction are set as fMA = 0.3 and fMI = 0.8, respectively. The optimization model and the shell-infill representation stated previously can be extended to 3D cases straightforwardly. The main difference is that the LSF is a 4D hyper surface in 3D shell-infill structural design. As indicated in Fig. 12, the red part of the structure is the volume intersected by 0 < U < 1, which is the shell of this structure. The blue part of the structure is the MA

volume of the infill represented by UMA > 1. The design is cut to exhibit the uniform microstructures inside the shell. The whole optimized shell-infill design is given in Fig. 13. The reader may find that some boundaries of the design domain is designed without shells. The authors leave the shell open here. On the one hand, it is straightforward for readers to observe the shell and infill of the structure. On the other hand, this open shell is extremely important in cleaning up the un-melted metal powder

This paper presents a multiscale level set topology optimization method for designing shell-infill structures. By taking advantage of the LSM, the shell and the infill are simply represented by two distinct level sets of a single LSF. The level set reinitialization scheme is utilized to preserve a uniform thickness of the shell. The homogenization method is applied to evaluate the effective elasticity matrix of the microstructural infill. The optimized shell-coated macro structure and micro infill are obtained simultaneously. Numerical examples show that the optimal design is different from the traditional multiscale optimization design when a shell is considered. Both 2D and 3D examples demonstrate the effectiveness of the proposed method. The presented method is not limited to the design of shell-infill structures for AM. It can also be extended to the design of bio-inspired structures with shell and infill. Acknowledgments This research is partially supported by National Basic Scientific Research Program of China (JCKY2016110C012), China Equipment Pre-research Program (41423010102), National Natural Science Foundation of China (51705166, 51675196 and 51721092), China Postdoctoral Science Foundation (2017M612446) and Program for HUST Academic Frontier Youth Team. References [1] Chai X, Chai H, Wang X, et al. Fused deposition modeling (FDM) 3D printed tablets for intragastric floating delivery of domperidone. Sci Rep 2017;7 (1):2829. [2] Bendsøe MP, Kikuchi N. Generating optimal topologies in structural design using a homogenization method. Comput Methods Appl Mech Eng 1988;71 (2):197–224. [3] Hassani B, Hinton E. A review of homogenization and topology optimization I— homogenization theory for media with periodic structure. Comput Struct 1998;69(6):707–17. [4] Bendsoe MP, Sigmund O. Topology optimization: theory, methods, and applications. Springer Science & Business Media; 2013. [5] Sigmund O, Maute K. Topology optimization approaches. Struct Multidiscip Optim 2013;48(6):1031–55.

172

J. Fu et al. / Computers and Structures 212 (2019) 162–172

[6] Osher S, Sethian JA. Fronts propagating with curvature-dependent speed: algorithms based on Hamilton-Jacobi formulations. J Comput Phys 1988;79 (1):12–49. [7] Van Dijk NP, Maute K, Langelaar M, et al. Level-set methods for structural topology optimization: a review. Struct Multidiscip Optim 2013;48(3):437–72. [8] Xie YM, Steven GP. A simple evolutionary procedure for structural optimization. Comput Struct 1993;49(5):885–96. [9] Xia L, Xia Q, Huang X, et al. Bi-directional evolutionary structural optimization on advanced structures and materials: a comprehensive review. Arch Comput Methods Eng 2016:1–42. [10] Guo X, Zhang W, Zhong W. Doing topology optimization explicitly and geometrically-a new moving morphable components based framework. J Appl Mech 2014;81(8):081009. [11] Zhang W, Yuan J, Zhang J, et al. A new topology optimization approach based on Moving Morphable Components (MMC) and the ersatz material model. Struct Multidiscip Optim 2016;53(6):1243–60. [12] Gaynor AT, Guest JK. Topology optimization considering overhang constraints: eliminating sacrificial support material in additive manufacturing through design. Struct Multidiscip Optim 2016;54(5):1157–72. [13] Langelaar M. An additive manufacturing filter for topology optimization of print-ready designs. Struct Multidiscip Optim 2017;55(3):871–83. [14] Guo X, Zhou J, Zhang W, et al. Self-supporting structure design in additive manufacturing through explicit topology optimization. Comput Methods Appl Mech Eng 2017;323:27–63. [15] Zhang W, Zhou L. Topology optimization of self-supporting structures with polygon features for additive manufacturing. Comput Methods Appl Mech Eng 2018;334:56–78. [16] Wang Y, Luo Z, Kang Z, et al. A multi-material level set-based topology and shape optimization method. Comput Methods Appl Mech Eng 2015;283:1570–86. [17] Vogiatzis P, Chen S, Wang X, et al. Topology optimization of multi-material negative Poisson’s ratio metamaterials using a reconciled level set method. Comput Aided Des 2017;83:15–32. [18] Cadman JE, Zhou S, Chen Y, et al. On design of multi-functional microstructural materials. J Mater Sci 2013;48(1):51–66. [19] Clausen A, Wang F, Jensen JS, et al. Topology optimized architectures with programmable Poisson’s ratio over large deformations. Adv Mater 2015;27 (37):5523–7. [20] Xia L, Breitkopf P. Recent advances on topology optimization of multiscale nonlinear structures. Arch Comput Methods Eng 2017;24(2):227–49. [21] Clausen A, Aage N, Sigmund O. Topology optimization of coated structures and material interface problems. Comput Methods Appl Mech Eng 2015;290:524–41. [22] Clausen A, Andreassen E, Sigmund O. Topology optimization of 3D shell structures with porous infill. Acta Mech Sin 2017;33(4):778–91. [23] Clausen A, Aage N, Sigmund O. Exploiting additive manufacturing infill in topology optimization for improved buckling load. Engineering 2016;2 (2):250–7. [24] Wu J, Clausen A, Sigmund O. Minimum compliance topology optimization of shell–infill composites for additive manufacturing. Comput Methods Appl Mech Eng 2017;326:358–75. [25] Allaire G, Dapogny C, Delgado G, et al. Multi-phase structural optimization via a level set method. ESAIM: Control, Optim Calculus Variat 2014;20 (2):576–611. [26] Wang Y, Kang Z. A level set method for shape and topology optimization of coated structures. Comput Methods Appl Mech Eng 2018;329:553–74. [27] Andreassen E, Andreasen CS. How to determine composite material properties using numerical homogenization. Comput Mater Sci 2014;83:488–95. [28] Rodrigues H, Guedes JM, Bendsoe MP. Hierarchical optimization of material and structure. Struct Multidiscip Optim 2002;24(1):1–10. [29] Liu L, Yan J, Cheng G. Optimum structure with homogeneous optimum trusslike material. Comput Struct 2008;86(13–14):1417–25. [30] Zhang W, Sun S. Scale-related topology optimization of cellular materials and structures. Int J Numer Meth Eng 2006;68(9):993–1011. [31] Alexandersen J, Lazarov BS. Topology optimisation of manufacturable microstructural details without length scale separation using a spectral coarse basis preconditioner. Comput Methods Appl Mech Eng 2015;290:156–82. [32] Deng J, Yan J, Cheng G. Multi-objective concurrent topology optimization of thermoelastic structures composed of homogeneous porous material. Struct Multidiscip Optim 2013;47(4):583–97.

[33] Groen JP, Sigmund O. Homogenization-based topology optimization for highresolution manufacturable microstructures. Int J Numer Meth Eng 2018;113 (8):1148–63. [34] Wang Y, Wang MY, Chen F. Structure-material integrated design by level sets. Struct Multidiscip Optim 2016;54(5):1145–56. [35] Sivapuram R, Dunning PD, Kim HA. Simultaneous material and structural optimization by multiscale topology optimization. Struct Multidiscip Optim 2016;54(5):1267–81. [36] Li H, Luo Z, Zhang N, et al. Integrated design of cellular composites using a level-set topology optimization method. Comput Methods Appl Mech Eng 2016;309:453–75. [37] Li H, Luo Z, Gao L, et al. Topology optimization for concurrent design of structures with multi-patch microstructures by level sets. Comput Methods Appl Mech Eng 2018;331:536–61. [38] Li H, Luo Z, Gao L, et al. Topology optimization for functionally graded cellular composites with metamaterials by level sets. Comput Methods Appl Mech Eng 2018;328:340–64. [39] Yan X, Huang X, Zha Y, et al. Concurrent topology optimization of structures and their composite microstructures. Comput Struct 2014;133:103–10. [40] Xia L, Breitkopf P. Concurrent topology optimization design of material and structure within FE2 nonlinear multiscale analysis framework. Comput Methods Appl Mech Eng 2014;278:524–42. [41] Xia L, Breitkopf P. Multiscale structural topology optimization with an approximate constitutive model for local material microstructure. Comput Methods Appl Mech Eng 2015;286:147–67. [42] Wang Y, Chen F, Wang MY. Concurrent design with connectable graded microstructures. Comput Methods Appl Mech Eng 2017;317:84–101. [43] Liu Q, Chan R, Huang X. Concurrent topology optimization of macrostructures and material microstructures for natural frequency. Mater Des 2016;106:380–90. [44] Vicente WM, Zuo ZH, Pavanello R, et al. Concurrent topology optimization for minimizing frequency responses of two-level hierarchical structures. Comput Methods Appl Mech Eng 2016;301:116–36. [45] Da DC, Cui XY, Long K, et al. Concurrent topological design of composite structures and the underlying multi-phase materials. Comput Struct 2017;179:1–14. [46] Wang MY, Wang X, Guo D. A level set method for structural topology optimization. Comput Methods Appl Mech Eng 2003;192(1–2):227–46. [47] Allaire G, Jouve F, Toader AM. Structural optimization using sensitivity analysis and a level-set method. J Comput Phys 2004;194(1):363–93. [48] Luo Z, Tong L, Kang Z. A level set method for structural shape and topology optimization using radial basis functions. Comput Struct 2009;87(7– 8):425–34. [49] Luo Z, Tong L, Ma H. Shape and topology optimization for electrothermomechanical microactuators using level set methods. J Comput Phys 2009;228(9):3173–81. [50] Wendland H. Scattered data approximation. Cambridge University Press; 2004. [51] Svanberg K. The method of moving asymptotes—a new method for structural optimization. Int J Numer Meth Eng 1987;24(2):359–73. [52] Burger M, Hackl B, Ring W. Incorporating topological derivatives into level set methods. J Comput Phys 2004;194(1):344–62. [53] Guo X, Zhang W, Zhong W. Stress-related topology optimization of continuum structures involving multi-phase materials. Comput Methods Appl Mech Eng 2014;268:632–55. [54] Zhang W, Song J, Zhou J, et al. Topology optimization with multiple materials via moving morphable component (MMC) method. Int J Numer Meth Eng 2018;113(11):1653–75. [55] Du Y, Li H, Luo Z, et al. Topological design optimization of lattice structures to maximize shear stiffness. Adv Eng Softw 2017;112:211–21. [56] Gao J, Li H, Gao L, et al. Topological shape optimization of 3D micro-structured materials using energy-based homogenization method. Adv Eng Softw 2018;116:89–102. [57] Feng YT, Owen DRJ, Peric´ D. A block conjugate gradient method applied to linear systems with multiple right-hand sides. Comput Methods Appl Mech Eng 1995;127(1–4):203–15. [58] Amir O, Stolpe M, Sigmund O. Efficient use of iterative solvers in nested topology optimization. Struct Multidiscip Optim 2010;42(1):55–72.