Multiscale concurrent topology optimization for cellular structures with multiple microstructures based on ordered SIMP interpolation

Computational Materials Science 155 (2018) 74–91

Contents lists available at ScienceDirect

Computational Materials Science journal homepage: www.elsevier.com/locate/commatsci

Multiscale concurrent topology optimization for cellular structures with multiple microstructures based on ordered SIMP interpolation

T

⁎

Yan Zhang, Mi Xiao, Hao Li, Liang Gao , Sheng Chu State Key Laboratory of Digital Manufacturing Equipment and Technology, Huazhong University of Science and Technology, 1037 Luoyu Road, Wuhan, Hubei 430074, China

A R T I C LE I N FO

A B S T R A C T

Keywords: Topology optimization Cellular structures Multiscale concurrent design Ordered SIMP interpolation Parametric level set method

This paper presents a novel multiscale concurrent topology optimization method for cellular structures with various types of microstructures to obtain a superior structural performance at an affordable computation cost. At macroscale, different microstructures are regarded as different materials in a multi-material topology optimization. An ordered solid isotropic material with penalization (SIMP) interpolation model, which is efficient for solving multi-material topology optimization problems without introducing any new variables, is employed to generate an overall structure layout with predetermined multiple materials. At microscale, each macroscale element is viewed as an individual microstructure. All macroscale elements with the same material are represented by a unique microstructure. A parametric level set method (PLSM) is applied to evolve the shape and topology of the representative microstructures, whose effective properties are evaluated by a numerical homogenization method. The connectivity of adjacent microstructures is guaranteed by a kinematical connective constraint. Using the proposed method, the macrostructural topology as well as the locations and configurations of microstructures can be simultaneously optimized. Numerical examples are provided to test the effectiveness of the proposed method.

1. Introduction Cellular structures made of porous microstructures are artificially engineered to have unconventional effective properties, such as being ultralight but with excellent structural performance, high specific stiffness/strength, and multifunctionalities in a broad range of missioncritical applications [1,2]. Over the past few decades, cellular structures have received increasing attention owing to the flexibility in tailoring specific effective properties by optimizing the layout and geometry of the microstructures rather than the constituent materials [3–5]. Particularly, cellular structures are commonly composed of a number of microstructures with ordinary material constituents, e.g., metals or plastics [3,4]. Thus, the geometry of the microstructures and their layouts within the macrostructure provide great potential for designing novel cellular structures with either superior performance or multifunctionalities by some computational design technologies [5,6], such as the popular topology optimization method [7–9]. Topology optimization is a highly developed tool for a diversity of both structural design and material design. It has been widely used in mechanical, automotive, and aerospace industries since the works of Bendsøe and Kichuchi [7] and Rozvany et al. [10]. Essentially, topology

⁎

optimization is a numerical iterative procedure that distributes the materials inside a fixed design reference domain to seek the best material layout, such that the objective function is optimized for a given set of boundary conditions. Up to the present, some topology optimization methods have been developed, such as the homogenization method [7,11], solid isotropic material with penalization (SIMP) [12,13], evolutionary structural optimization (ESO) [14,15], and level set method (LSM) [16,17]. Because the superior mechanical properties of cellular structures can be achieved by tailoring the microscopic configurations of material unit cells, the term material here refers to the effective properties of macrostructures rather than those of the constituent material in the microstructures. In recent years, a powerful systematic design approach has been developed by incorporating a numerical homogenization method into topology optimization to optimize the geometrical configuration of the material unit cell in order to achieve either a desired effective property or an extreme effective property, such as negative Poisson’s ratios [3,18], extreme shear or bulk moduli [19,20], extreme thermal expansion coefficients [21,22], multiphysics and multifunctional properties [23,24], functionally graded properties [5,25,26] and many others [27,28]. A further step has been made in [29,30],

Corresponding author. E-mail address: [email protected] (L. Gao).

https://doi.org/10.1016/j.commatsci.2018.08.030 Received 26 April 2018; Received in revised form 8 July 2018; Accepted 14 August 2018 0927-0256/ © 2018 Elsevier B.V. All rights reserved.

Computational Materials Science 155 (2018) 74–91

Y. Zhang et al.

macrostructure is composed of various types of microstructures, have been conducted. For instance, Zhang and Sun [46] proposed a design element method and Li et al. [4] developed a hierarchical multiscale topology optimization method for both structure and materials, where various microstructures of each layer are designed for a layered structure. Xu and Cheng [47] investigated a two-scale concurrent design method with multiple micro heterogeneous materials based on a criterion for principal stress orientation. Li et al. [48] proposed a concurrent design for cellular structures with multi-patch microstructures using a discrete element density optimization method, aiming to improve the structural performance. Based on a discrete element density optimization method, Deng and Chen [49] developed a robust concurrent topology optimization approach to design a multiscale structure consisting of multiple porous materials under a random loading uncertainty. From the viewpoint of the above studies, it is reasonable and efficient to achieve a multiscale concurrent design for a macrostructure with various types of microstructures. Essentially, different types of microstructures, which are used to assemble macrostructures in multiscale concurrent topology optimization, can be regarded as different types of materials in a multi-material topology optimization. Therefore, this study proposes a novel multiscale concurrent topology optimization method for cellular structures composed of multiple microstructures, in which both the macrostructural topology and locations of microstructures are determined by a multi-material topology optimization algorithm. So far, several popular multi-material models have been successfully used in topology optimization multi-material structures, including the conventional multi-material density model [50,51], the “color” level set model [52], and the multi-material level set model [53]. Recently, Zhang et al. [54] presented a bi-material model based on independent point-wise density interpolation (iPDI) [55] for bimaterial microstructural topology optimization with chiral auxetic metamaterials. As an alternative multi-material model, the ordered SIMP interpolation model used in this study is efficient for solving multi-material topology optimization problems without introducing any new variables. Moreover, it can significantly reduce the risk of obtaining a local optimum when a considerable number of microstructures are used in the macro design domain [56]. In this framework, all macroscale elements with the same material (i.e., the same relative element density) are represented by a unique material microstructure. The relative element density of each macroscale element is regarded as the volume constraint of the microstructural design. At the microscale, only a given number of representative microstructures are required to be topologically optimized, where the PLSM integrated with a numerical homogenization method [57] is used. In this way, the macrostructural topology, as well as the configurations of microstructures and their distribution within the macrostructure can be simultaneously optimized.

which is to optimize the material microstructures for the macrostructure with known boundary conditions. It is noted that the majority of the above works mainly focus on the design of the material microstructure itself and essentially belong to a single-scale design approach. Actually, for pursuing higher-performance structures, it is desirable to simultaneously consider the material microstructures and the boundary conditions of the macrostructure to fully explore a vast design space of both scales. Hence, to effectively implement the concurrent design of the macrostructure and its materials for seeking the optimal distribution of the best material, the demand for an efficient multiscale computational method has become increasingly intense. From the perspective of a multiscale concurrent design, the optimal designs should be simultaneously implemented at both the macroscale and microscale, such that the most compatible configurations of the macrostructure and microstructure can be achieved in response to the boundary conditions. Currently, some interesting studies have been conducted, where the concurrent design optimization is performed for both the macrostructural topology and the element-wise varying microstructures. For example, Rodrigues and his co-workers [31,32] proposed a hierarchical computational procedure for simultaneously optimizing both the material distribution and material mechanical property for each microstructural element; then, the work had been extended to 3D [33]. Based on a free material optimization, Schury et al. [34] proposed a two-scale optimization method to optimize an inhomogeneous macrostructure consisting of many structurally graded materials. Xia and Breitkopf [35,36] investigated a concurrent design of a structure and material under the FE2 nonlinear multiscale analysis framework, in which the microstructures vary from point to point. Sivapuram et al. [37] introduced a new multiscale optimization approach to simultaneously optimize a structure and material, where the positions of all microstructures are predefined and fixed during the multiscale concurrent optimization. Wang et al. [38] developed a multiscale concurrent design method to achieve an optimal structure based on connectable graded microstructures, which are spatially varying. In most of the above studies, the microstructures vary from element to element throughout the entire macro design domain, and the design freedom of the concurrent design for both structure and material can be fully explored. However, the considerable number of material microstructures to be topologically optimized results in an intensive computational cost. To address such issue, several studies assume that the material microstructures are identical throughout the entire macrostructure. For instance, Fujii et al. [39] investigated the compliance minimization problem of a macrostructure by optimizing a uniform material microstructure. Liu et al. [40] presented a concurrent topology optimization method to simultaneously achieve the optimal structure and material microstructure for minimum system compliance with uniformly distributed microstructures. Wang et al. [41] proposed a structure–material integrated design method to simultaneously generate optimized material distribution patterns and optimized configuration of periodically arranged microstructures over a material region of the structure. Chen et al. [42] presented a concurrent design approach of structure and material microstructure, in which the material microstructure is assumed to be uniform in the macrostructure. Yan et al. [43] proposed a two-scale topology optimization approach to concurrently design a structure under the assumption that the macrostructure is uniformly composed of a porous thermal insulation material. Long et al. [44,45] proposed a concurrent design method for optimizing composite structures consisting of a unique microstructure, considering the effect of constituent phases with distinct Poisson’s ratios. Obviously, compared with a macrostructure consisting of elementwise microstructures, a macrostructure with a uniform material microstructure can significantly reduce the computational cost. However, this design approach may significantly restrict the design space at both scales, leading to an inferior structural performance. Recently, for pursuing higher-performance structures with an affordable computational cost, several interesting works, where the

2. Multiscale concurrent topology optimization Fig. 1 schematically illustrates the proposed multiscale concurrent optimization of a cantilever composed of multiple types of microstructures, aiming to maximize the macrostructural performance subject to the global volume constraint. Here, the macrostructural topology as well as the locations and configurations of the microstructures can be concurrently optimized. The proposed method optimizes a representative microstructure for a group of microstructures having the same relative density (local volume fraction). It can be inferred that the proposed multiscale concurrent optimization can achieve a design with superior macrostructural performance than the method with a uniform microstructure distribution in the entire macro design domain. Furthermore, it can achieve a significantly lower computational cost than the method with microstructures varying from element to element in the entire macrostructure. In the following, to facilitate the description and discussion, the superscripts “MA” and “MI” are applied to denote 75

Computational Materials Science 155 (2018) 74–91

Y. Zhang et al.

Elastic modulus

Ee

Material 3

Material 2

Material 1

Void

Intermediate density

ρe

Fig. 3. Ordered multi-material SIMP interpolation of elastic modulus.

F

weight by using the classical SIMP interpolation while encouraging the convergence to a black-and-white design. Based on the attractive performance of the classical SIMP interpolation, the ordered interpolation of the elastic modulus of multiple materials can be constructed. As depicted in Fig. 3, the multiple materials are sorted in ascending order of the normalized density variable ρTi : Fig. 1. Schematic of the macrostructure with six representative microstructures.

ρi = ρTi / ρmax , (i = 1, 2, 3, …, w )

(2)

where ρmax is the maximum among all the candidate densities, and w is the number of materials. Then, the ordered SIMP interpolation of the elastic modulus of multiple materials is formulated as

the physical quantities at macroscale and microscale, respectively. 2.1. Ordered multi-material SIMP interpolation at macroscale

Ee (ρe ) = F (ρe ) E0 In the classical SIMP interpolation, the relationship between the element density ρe and the elastic modulus Ee in the equilibrium analysis is written as

Ee (ρe ) = ρeP E0, ρe ∈ [ρmin , 1], P > 1 and e = 1, 2, …, NE

+ BE is the extended power function. The scaling where F (ρe ) = coefficient AE and translation coefficient BE for ρe ∈ [ρi , ρi + 1] are given as

(1)

AE =

where NE denotes the total number of finite elements (FEs), E0 the modulus of solid element, ρmin a small positive value (e.g., 0.001) to avoid singularity in the numerical process, and P is a penalty factor. As shown in Fig. 2, for the interpolation curve in Eq. (1), as the element density ρe in the neighborhood of void (0,0) approaches to 0, the rate of decrease in elastic modulus Ee is small. In contrast, as the element density ρe in the neighborhood of solid (1,1) approaches to 1, the rate of increase in elastic modulus Ee becomes large. This means that the optimized structure tends to have a higher stiffness and lighter

Increase

Largely Increase

Decrease

Ee

Ei−Ei + 1 ρiP −ρiP+ 1

and BE = Ei−AE ρiP

(4)

where Ei is the elastic modulus of the ordered material i. According to the variable thickness sheet (VTS) method [8,12,58], the elastic modulus Ei can be calculated with the relative density ρi , i.e., Ei = ρi E0 . From the above analysis, it is noted that the topologically optimized multi-material structure consisting of any type of material can be obtained by means of the ordered multi-material SIMP interpolation model without introducing any new variables. In this study, an efficient multi-material topology optimization method is employed to optimize the material distribution at macroscale, where the macro element with the same material is represented by a unique microstructure. According to the ordered multi-material SIMP interpolation of the MA MA elastic modulus in Eq. (3), the elasticity tensor Dijkl (ρe ) can be expressed as

Largely Decrease

Material 1(1,1)

Elastic modulus

(3)

AE ρeP

∼0 MA MA Dijkl (ρe ) = F (ρeMA )(De )ijkl , ρeMA ∈ [ρmin , 1]

dEe d ρe dEe d ρe

, so ρe Ee

, so ρe

denotes the relative density of the e macro element, and where ∼0 (De )ijkl is a variable used to identify the equivalent base material property for the eth macro element. During the concurrent optimization process, the elasticity tensor MA MA Dijkl (ρe ) of the eth macro element is assumed to be equivalent to the H (uMI , Φ MI ) of one microstructure with the same elasticity tensor Dijkl density, i.e.,

Ee

MA MA H Dijkl (ρe ) = Dijkl (uMI , Φ MI )

Void(0,0) Intermediate density

(5) th

ρeMA

ρe

H Dijkl (uMI ,

Φ MI )

(6)

can be evaluated via the numerical homogenizawhere MA MA tion method. Compared to the isotropic material properties Dijkl (ρe ) of the ordered SIMP interpolation model, it is noted that the

Fig. 2. Classical power function interpolation of elastic modulus. 76

Computational Materials Science 155 (2018) 74–91

Y. Zhang et al. H (uMI , Φ MI ) may be orthotropic to homogenized elasticity tensor Dijkl flexibly supply different directional stiffness in concurrent optimization. ∼0 According to Eqs. (5) and (6), (De )ijkl can be calculated as follows:

∼0 H (De )ijkl = Dijkl (uMI , Φ MI )/ F (ρeMA )

Actually, only the normal velocity field υ determines the movement of the structure boundary, and the first-order H–J PDE can be rewritten as [16,17]

∂Φ MI (x, t ) −υ·|∇Φ MI (x, t )| = 0, ∂t

(7)

υ = V·n =

At the macroscale, the density-based topology optimization problem to minimize the compliance objective subject to structural weight constraints, can be stated as

)=

=

∑ e=1

∑

MA ρeMA V0−Vmax ⩽ 0,

e=1

F = KU, ρmin ⩽ ρeMA ⩽ 1.

(13)

= [φ1MI (x), φ2MI (x), …, φNMI (x)] [α1MI (t ), α 2MI (t ), …, αNMI (t )]T is

(i = 1, 2, …, N )

(14)

where ri (x) = ‖x−x i ‖/ rd , in which rd is the radius of the support domain of the knot x i , and ‖x−x i ‖ describes the distance from the current sample knot x to the knot x i . As indicated in Eq. (14), only the knots within the support domain of the current knot x contribute to its function value. It can be observed from Eq. (13) that the time and spatial terms associated with the dynamic level set function are decoupled via the CSRBF interpolation. Substituting Eq. (13) into Eq. (11), the H–J PDE is transformed into a new form of an ordinary differential equation:

where C is the structural compliance and G denotes the global volume constraint. The design variable ρeMA denotes the relative density of the eth macro element, and NE denotes the total number of FEs over the entire macrostructure. F, U, and K are the external load vector, displacement vector, and global stiffness matrix at macroscale, respectively. F (•) is the ordered multi-material SIMP interpolation function. u and k 0 are respectively the element displacement vector and the element stiffness matrix corresponding to a solid material elasticity MA property. V0 is the volume of element, and Vmax denotes the maximum volume of the macrostructure.

φ MI (x)

2.3. Parametric level set method at microscale

dα MI (t ) −υ ·|∇φ MI (x) α MI (t )| = 0 dt

(15)

Thus, the normal velocity field can be rewritten as

In the present work, PLSM is employed to evolve the shape and topology of the material microstructure at microscale. The first step of PLSM is to implicitly represent the moving boundary of the structure as the zero level set of a higher-dimensional level set function Φ MI , which is Lipschitz-continuous [16,17]. Defining a fixed Eulerian reference domain D MI ⊂ Rd (d = 2, 3) , which contains all admissible shapes ΩMI , each region can be described by using the level set function Φ MI as follows:

Φ MI (x) > 0, ∀ x ∈ ΩMI ⧹∂ΩMI (Solid) ⎧ ⎪ MI Φ (x) = 0, ∀ x ∈ ∂ΩMI (Boundary) ⎨ MI ⎪ Φ (x) < 0, ∀ x ∈ D MI ⧹ΩMI (Void) ⎩

υ=

φ MI (x) α̇ MI |∇φ MI (x) α MI (t )|

, where, α̇ MI (t ) =

dα MI (t ) dt

(16)

It is obvious that the normal velocity field is naturally extended to the entire design domain, as φ MI (x) , α MI (t ) , and α̇ MI (t ) in the above equation are all evaluated over the total set of knots in this design domain rather than only the knots on the boundary of the structure. Therefore, no additional numerical scheme is required to extend the velocity field from the boundary to the entire design domain. Furthermore, the expansion coefficients α MI are the only unknown terms to be recognized as the design variables. Thus, the original shape and topology optimization has been converted into a general size optimization with α MI regarded as a set of generalized size parameters [60,61]. In this way, many well-established efficient gradient-based algorithms [62,63] can be directly applied to solve the PLSM-based optimization model.

(9)

∂ΩMI

where x denotes the point coordinates in and denotes the boundary of the structure. When introducing the pseudo-time t to enable the dynamic motion of the shape deformations, actually the moving of the structure boundary toward its optimum is equivalent to solving the following first-order Hamilton–Jacobi partial differential equation (H–J PDE):

∂Φ MI (x, t ) + V ·∇Φ MI (x, t ) = 0 ∂t

φiMI (x) αiMI (t )

φiMI (x) = max{0, (1−ri (x))}4 (4ri (x) + 1),

(8)

D MI ,

(12)

is the vector of the CSRBFs where the vector of generalized and α MI (t ) = expansion coefficients. N is the total number of fixed knots over the design domain. The CSRBFs mentioned above are referred to as a type of widely applied radial basis functions with C2 continuity [59]:

F (ρeMA ) uT k 0u

NE

S. t . G (ρeMA ) =

∑ i=1

φ MI (x)

=

⎟

N

NE

Min

⎜

Φ MI (x, t ) = φ MI (x) α MI (t ) =

Find ρeMA (e = 1, 2, …, NE ) UT KU

d x (t ) d x (t ) ⎞ ⎛ ∇Φ MI (x, t ) ⎞ ·n = ⎛ − MI dt ⎝ dt ⎠ ⎝ |∇Φ (x, t )| ⎠

Here, n is the unit outward normal vector pointing from the interior to the exterior on the boundary of the structure, owing to the sign convention of the level set function in this study [5]. To eliminate the numerical difficulties of directly solving the firstorder H–J PDE [17], the completely supported radial basis function (CSRBF) [59] is introduced to approximate the level set function Φ MI (x, t ) , as follows:

2.2. Macrostructural topology optimization

FT U

(11)

where the normal velocity field υ can be denoted as follows:

∼0 In this study, it is noted that (De )ijkl is no longer the natural property of the solid material element e, as it practically acts as a temporal value varying in the concurrent optimization. Before macrostructure optimi∼0 zation during the concurrent optimization, the terms (De )ijkl of each microstructure need to be calculated by using Eq. (7), where H Dijkl (uMI , Φ MI ) and F (ρeMA ) can be obtained from the previous iteration. ∼0 Then, the new (De )ijkl of each microstructure will be used to update ρeMA .

C (ρeMA

Φ MI (x, 0) = Φ0MI (x)

2.4. Microstructural topology optimization From the multi-material distribution obtained by Eq. (8), the multimaterials (such as M materials) are recognized as multiple macro elements with different relative densities, and the macro element with the same relative density is represented by a type of representative microstructure. It can be observed that there are M-2 types of different

(10)

where V = d x/ dt consisting of the normal velocity field υ and the tangential velocity field γ is the velocity field of the structure boundary. 77

Computational Materials Science 155 (2018) 74–91

Y. Zhang et al.

representative microstructures (excluding void and solid) within the entire design domain to be topologically optimized at microscale. The macroscale relative element densities are regarded as the volume fraction constraint for the representative microstructures. Thus, the microscale topology optimization of the microstructures can be stated as

a (uMA , v MA, uMI , Φ MI ) =

l (v MA) =

MI −ρm V0MA ⩽ 0, S. t . Gm (α MI ) = ∫ΩMI H (Φ MI (α MI )) dΩm m

m, n

max

(17)

where M is the number of total types of representative microstructures and N is the total number of knots in the micro level set grid for a representative microstructure. αmMI, n is the actual design variable at microscale, which indicates the expansion coefficient of the CSRBF interpolation of the nth knot in the micro level set grid of the mth representative microstructure. εij and εkl represent the strain fields, where i, j, k, l = 1, 2, …, d and d is the spatial dimension. uMA is the macroscale displacement field, which is calculated by substituting the effecH of the material microstructure into the corretive properties Dijkl H sponding state equation F = K (Dijkl ) U at the macroscale, as given in MA Eq. (8). v is the macroscale virtual displacement field belonging to the kinematically admissible displacement field U (uMA) . uMI is the MI of the microscale displacement field defined at the design domain Ωm mth representative microstructure. V0MA denotes the volume of a macroscale FE with solid material. ρm indicates the relative density of the mth macroscale FE. Gm is the local volume fraction constraint for the mth representative microstructure, subjected to the volume fraction of the ∼ MI = 0.001 and α ∼ MI = 1 are the macroscale element, i.e., ρm V0MA . α min max ∼ MI respectively. α ∼ MI are regularized design lower and upper bounds of α m, n m, n variables, which will be used in the optimization algorithm to facilitate the numerical implementation.The effective properties of the microstructures evaluated via the numerical homogenization method may not be accurate because of the non-periodicity between two different neighboring microstructures, where the macrostructure is composed of various microstructures. However, as shown in [25,26], when the change in property gradient between two different neighboring microstructures is sufficiently small, the numerical homogenization method can practically be applied to approximate the effective properties of the material microstructure. In this study, based on the small parameter perturbation of the displacement, the effective elasticity tensor for each material microstructure occupying domain ΩMI is given by H Dijkl (uMI , Φ MI (α MI ))

=

1 |ΩMI |

∫Ω

∀ v MI (kl) ∈ U (ΩMI )

f v MAdΩMA +

∫Ω

MA

τ v MAdΓMA

(21)

As a gradient-based method, the optimality criteria (OC) algorithm [13,62] is employed to update the design variables at both scales in this study. It is necessary to calculate the first-order derivatives of the objective function and the constraints with respect to the design variables for guiding the search direction of the OC algorithm, where the design variables are the element densities for the macrostructure optimization problem and the expansion coefficients of the CSRBF interpolation for the microstructure optimization. 3.1. Sensitivity analysis at macroscale The sensitivity of the structure compliance objective function with respect to the macroscale design variable ρeMA is expressed by [8]

∂C (ρeMA ) ∂ρeMA

=−

dF (ρeMA ) dρeMA

uT k 0u

(22)

F (ρeMA )

Based on Eq. (3), the derivative of macroscale design variable ρeMA can be given as

dF (ρeMA ) dρeMA

∗ 0(ij ) (εpq −εpq (uMI (ij) )) Dpqrs (εrs0(kl)

∂C (ρeMA ) ∂ρeMA

(18)

= AE P (ρeMA ) P − 1

with respect to the

(23)

= −AE P (ρeMA ) P − 1uT k 0u

(24)

The sensitivity of the volume constraint with respect to the macroscale design variable ρeMA is obtained by

∂G (ρeMA ) ∂ρeMA

= V0

(25)

As shown in Fig. 3 and Eqs. (3) and (4), the interpolation function of the elastic modulus is non-smooth, leading to discontinuity of their first-order derivatives in Eqs. (22)–(24) at the interpolation points. Although the discontinuity of derivatives could reduce mathematical rigor, and in general, cause numerical instabilities to the optimization process, it is observed that this is not an issue in the numerical example that we attempted. This may be because design variables hitting exactly an interpolation point is an extremely rare event with the 18 digits, double precision numbers used to implement the density design variables ρeMA . Even if it occurs, the dynamic move limit in Eq. (26) (see

0(ij ) ∗ (εpq −εpq (uMI (ij) )) Dpqrs εrs∗ (v MI (kl) ) H (Φ MI ) dΩMI

= 0,

MA

3. Sensitivity analysis and optimization algorithm

where ΩMI is the design domain of the microstructure, and |ΩMI | is the area or the volume of the microstructure Dpqrs corresponding to the elasticity tensor of the solid material. H (Φ MI ) is the Heaviside function 0(ij ) is the to indicate different parts of the design domain [16,17,60]. εpq macroscopic test unit strains containing three unit vectors, namely (1 0 0)T , (0 1 0)T , and (0 0 1)T , for two dimensions. The displace∗ ment uMI (ij) of locally varying strain fields εpq (uMI (ij) ) can be calculated by solving the following equations at microscale with the given mac0(ij ) : roscopic strain εpq MI

∫Ω

Thus, Eq. (22) can be rewritten as MI

−εrs∗ (uMI (kl) )) H (Φ MI ) dΩMI

∫Ω

H εij (uMA) Dijkl (uMI , Φ MI ) εkl (v MA) dΩMA

where f is the body force in the macro domain ΩMA , and τ is the traction along the boundary ΓMA . As mentioned above, the topology optimization for cellular structure with multiple microstructures has been formulated through a multiscale concurrent manner. The macroscale and microscale optimizations are actually bridged by the relative density of macroscale element (the volume fraction of the microstructure) and the effective property of the microstructures. The relative density of macroscale element can be obtained by the multi-material topology optimization algorithm integrated with the ordered SIMP interpolation model, and the effective property of the microstructures is evaluated via the numerical homogenization method. As PLSM is employed to topologically optimize the microstructures in the entire macrostructure design domain, the multiscale design result with the distinct interfaces and smoothed boundaries can facilitate the fabrication without additional extensive post-processing.

H Min C (α MI ) = ∫ΩMA εij (uMA) Dijkl (uMI , Φ MI (α MI )) εkl (uMA) dΩMA

min

MA

(20)

Find αmMI, n (m = 1, 2, …, M ; n = 1, 2, …, N )

a (uMA , v MA, uMI , Φ MI ) = l (v MA), ∀ v MA ∈ U (uMA), ∼ MI ⩽ α ∼ MI ⩽ α ∼ MI . α

∫Ω

(19)

The bilinear energy form and the linear load form in Eq. (17) are respectively expressed as 78

Computational Materials Science 155 (2018) 74–91

Y. Zhang et al.

Substituting Eqs. (29)–(31) into Eq. (32) yields

Section 3.2) can also suppress the numerical instabilities [56].

∫Ω

3.2. Updating the design variables at macroscale

MA MA (k ) ρe

MA ⎧ max[ρmin , ρeMA (k ) −move (k ) ] < [BeMA (k ) ]ς ρeMA (k ) MA ⎨[B MA (k ) ]ς ρ MA (k ) < min[ρ MA (k ) + move (k ) 1] e e ⎩ e MA MA (k ) ρe

MA

H ∂Dijkl (uMI , Φ MI (α MI ))

εij (uMA)

∂C (α MI ) =− ∂t

if

if [BeMA (k ) ]ς

∫Ω

∂t

∫Ω

MA

εij (uMA)

H MI MI ∂Dijkl (um , Φm (α MI ))

⩽ max[ρmin , ρeMA (k )

∂t

−move (k ) ]

∂ρeMA (k )

/max(ϑ, ΛMA (k )

∂G (ρeMA ) ∂ρeMA (k )

=

where

MI δ (Φm )

=

H MI MI ∂Dijkl (um , Φm (α MI ))

)

∂t

(27)

+ εij ∂a ∂t

Φ MI (α MI ))

∂t

∂l = ∂t

∂t

∂t

N

=

MAdΩMA + MA f v̇

∫Ω

MA

H MI MI ∂Dijkl (um , Φm (α MI ))

∂αmMI, n εkl (uMA)] dΩMA

H εij (uMA) Dijkl (uMI , Φ MI (α MI )) εkl (v̇ MA) dΩMA = MA

τ v̇ MAdΓMA

⎛⎜ 1 MI ⎝ |Ωm |

∫Ω

MI m

∗ 0(ij ) MI (ij ) −εpq (εpq (um )) Dpqrs (εrs0(kl)

∑

H MI MI ∂Dijkl (um , Φm (α MI )) ∂αmMI, n · MI ∂t ∂αm, n

(37)

=

1 MI |Ωm |

∫Ω

MI m

∗ 0(ij ) MI (ij ) −εpq (εpq (um )) Dpqrs (εrs0(kl)

MI (kl) MI MI −εrs∗ (um )) φmMI, n (x)·δ (Φm ) dΩm

M

∂αmMI, n

(29)

(38)

∫Ω

MA

N

∫Ω

MA

εij (uMA)[ ∑

H MI MI ∂Dijkl (um , Φm (α MI ))

n=1

∂αmMI, n

·

] εkl (uMA) dΩMA

(39)

On the other hand, the derivative of the objective function C (α MI ) with respect to t can be expressed by the chain rule as

(30)

∂C (α MI ) = ∂t

f v̇ MAdΩMA

M

N

∑ ∑ m=1 n=1

MI ∂C (α MI ) ∂αm, n · MI ∂t ∂αm, n

(40)

Again, comparing the corresponding terms in Eqs. (39) and (40), the derivative of the objective function C (α MI ) with respect to the design variables αmMI, n can be obtained as

(31)

Differentiating the macroscale equilibrium in Eq. (17) equation with respect to t yields

∂a (uMA , v MA, uMI , Φ MI ) ∂l (v MA) = ∂t ∂t

N

=

(28)

∂C (α MI ) =−∑ ∂t m=1

εkl (v MA) dΩMA

∫Ω

(35)

Recalling Eq. (34) and substituting Eq. (37) into Eq. (34) yields

MAd ΓMA MA τ v̇

+

∗ 0(ij ) MI (ij ) −εpq (εpq (um )) Dpqrs (εrs0(kl)

Comparing the corresponding terms in Eqs. (36) and (37), it can be H found that the derivative of the effective elasticity tensor Dijkl with respect to the design variables αmMI, n can be expressed as

The following conjugate equation can be established considering the terms v̇ MA contained in Eqs. (29) and (30):

∫Ω

∑

n=1

∂t

∫Ω

MI m

is the derivative of the Heaviside function

n=1

H MI MI ∂Dijkl (um , Φm (α MI ))

H + εij (uMA) Dijkl (uMI , Φ MI (α MI )) εkl (v̇ MA)] dΩMA

+ ∫ΩMA εij (uMA)

∫Ω

(36)

H = ∫ΩMA [εij (u̇ MA) Dijkl (uMI , Φ MI (α MI )) εkl (v MA)

H (uMI , Φ MI (α MI )) ∂Dijkl

1 MI |Ωm |

On the other hand, the derivative of the effective elasticity tensor H Dijkl with respect to t can be expressed by the chain rule as

H [2εij (u̇ MA) Dijkl (uMI , Φ MI (α MI )) εkl (uMA)

H ∂Dijkl (uMI , (uMA)

(34)

∂αmMI, n MI (kl) MI MI ⎞ −εrs∗ (um )) φmMI, n (x)·δ (Φm ) dΩm ⎟· ⎠ ∂t

In this section, the shape derivative is introduced to conduct the sensitivity of boundary perturbations with respect to the microscale time variable t by following a method similar to that by Allaire et al. [17]. The shape derivatives of C (α MI ) , a (uMA , v MA, uMI , Φ MI ) , and l (v MA) can be respectively calculated by MA

εkl (uMA) dΩMA

MI H (Φm ) , and ξ is chosen as 2–4 times the mesh size from numerical experience [48,64]. υmMI denotes the normal velocity within the design domain of the mth representative microstructure. By substituting the normal velocity υmMI defined in Eq. (16) into Eq. (35), we obtain

3.3. Sensitivity analysis at microscale

∫Ω

∂t

ξ 1 · MI 2 2 π (Φm ) +ξ

where ϑ is a small positive constant to avoid zero terms, and ΛMA is a Lagrangian multiplier that can be found by using a bi-sectioning algorithm.

∂C (α MI ) = ∂t

H ∂Dijkl (uMI , Φ MI (α MI ))

MI (kl) MI MI MI −εrs∗ (um )) υmMI ·|∇Φm |δ (Φm ) dΩm

where k denotes the iteration step number starting from 1; move (k ) is a positive and iteration-dependent dynamic moving limit, i.e., move (k ) = max(β (k ) move (0), m min ) , where move (0) = 0.15 is the initial move limit, β = 0.96 is the scale coefficient, and m min = 0.001 is the minimal move limit to improve the accuracy of the OC algorithm in this study. The numerical damping parameter ς MI = 1/2, (0 < ς MI < 1) is employed to stabilize the iterative process, and BeMA (k ) can be defined by

∂C (ρeMA )

(33)

Based on [3,4,48], the derivative of the effective elasticity tensor H Dijkl of the material microstructure with respect to t can be expressed as

(26)

BeMI (k ) = −

εkl (v MA) dΩMA

Considering that the compliance optimization problem is self-adjoint, the shape derivative of the objective function can be rewritten by substituting Eq. (33) into Eq. (28) as follows:

ρeMA (k + 1) if min[ρeMA (k ) + move (k ) 1] ⩽ [BeMA (k ) ]ς

H (uMI , Φ MI (α MI )) εkl (v MA) dΩMA εij (u̇ MA) Dijkl

=−

Based on the obtained sensitivities of the structure compliance objective and the volume constraint, the following well-established OC algorithm is applied to update the macroscale design variable ρeMA : ⎧ min[ρeMA (k ) + move (k ) 1], ⎪ ⎪ ⎪[BeMA (k ) ]ς MA ρ MA (k ) , e ⎪ ⎪ = ⎨ ⎪ ⎪ ⎪ max[ρ , ρ MA (k ) −move (k ) ], min e ⎪ ⎪ ⎩

MA

∂C (α MI ) =− ∂αmMI, n

(32)

∫Ω

MA

εij (uMA)

H MI MI ∂Dijkl (um , Φm (α MI ))

∂αmMI, n

εkl (uMA) dΩMA

(41)

Similarly, the derivative of the local volume constraints Gm (α MI )

79

Computational Materials Science 155 (2018) 74–91

Y. Zhang et al.

with respect to the design variables αmMI, n can be given by

∂Gm (α MI ) = ∂αmMI, n

∫Ω

MI m

∼ MI (k + 1) α m, n

MI MI φmMI, n (x)·δ (Φm ) dΩm

∼ MI (k ) ⎧ min[(η MI + 1) α m, n ⎪ ⎪ ⎪ ⎪ ∼ MI (k ), ⎪[BmMI, n(k ) ]ς MI α m, n ⎪ = ⎨ ⎪ ⎪ ⎪ ∼ MI (k ) ⎪ max[(1−η MI ) α m, n ⎪ ⎪ ⎩

(42)

3.4. Updating the design variables at microscale As indicated in Eq. (17), it is not very easy to specify the fixed lower and upper bounds of the actual design variables αmMI, n during the optimization. Hence, in order to facilitate the numerical implementation at the microscale, the well-established OC algorithm is also applied to ∼ MI with the fixed lower update the regularized design variables α m, n MI ∼ MI ∼ bounds α min = 0.001 and upper bounds α max = 1. Then, the actual design ∼ MI . variables αmMI, n are renewed by using the updated α m, n ∼ MI (k ) can be calculated as follows: The regularized design variables α m, n MI (k ) MI (k ) ∼ MI (k ) = αm, n −αm,min α m, n (k ) MI (k ) αm,max−αmMI,min

where k denotes the iteration step number starting from 1. (k ) αmMI,max are defined by (k ) αmMI,min = 2 × min(αmMI, n(k ) ),

(k ) αmMI,max = 2 × max(αmMI, n(k ) )

[BmMI, n(k ) ]ς

MI ∼ MI (k ) α m, n

if ∼ MI (k ) α ∼ MI ] < [B MI (k ) ]ς MI α ∼ MI (k ) ⎧ max[(1−η MI ) α min m, n m, n m, n ∼ MI ], α min

∼ MI (k ) < min[(η MI + 1) α ∼ MI (k ) α ∼ MI ] ⎨[B MI (k ) ]ς MI α max m, n m, n ⎩ m, n MI (k ) ς MI ∼ MI (k ) if [B ] α ⩽ m, n

m, n

∼ MI (k ) α ∼ MI ] max[(1−η MI ) α min m, n

(45) where the moving limit η MI (0 < η MI < 1) and the damping parameter ς MI are employed to stabilize the iterative process, and BmMI, n(k ) can be defined by

BmMI, n(k ) = −

(43) (k ) αmMI,min

∼ MI ], if min[(η MI + 1) α ∼ MI (k ) α ∼ MI ] ⩽ α max max m, n

MI ∂C (α MI ) MI (k ) ∂Gm (α ) /max(ϑ, Λm ) MI MI ∂αm, n ∂αm, n

(46) MI Λm

and

is a where ϑ is a small positive constant to avoid zero terms, and Lagrangian multiplier that can be updated by using a bi-sectioning algorithm. Finally, when the convergent criterion of OC algorithm is satisfied, ∼ MI (k + 1) the actual design variables αmMI, n(k + 1) can be obtained by α m, n

(44)

Based on the Kuhn–Tucker conditions [8,62], the regularized design ∼ MI can be iteratively updated by using a heuristic scheme as variables α m, n follows:

∼ MI (k + 1) × (α MI (k ) −α MI (k ) ) + α MI (k ) αmMI, n(k + 1) = α m, n m,max m,min m,min

Fig. 4. Flowchart of numerical implementation of the proposed method. 80

(47)

Computational Materials Science 155 (2018) 74–91

Y. Zhang et al.

L

Design Domain

Microstructure 1

H

F

Microstructure 2 Fig. 6. Sketch of MBB beam. Table 1 Types of materials and corresponding relative density.

Kinematical Connectors Fig. 5. Sketch of predefined connectors between two adjacent microstructures.

4. Numerical implementation The flowchart of the proposed method, which mainly contains two stages, i.e., material distribution optimization at macroscale and concurrent optimization at both scales, is shown in Fig. 4. At macroscale, MA based on the elastic tensor Dijkl calculated via Eq. (5), which is the ordered multi-material interpolation model, the macro FE analysis is implemented to calculate the displacement field U and the macrostructure compliance C. Then, the sensitivities of the macrostructure compliance and the global volume constraint with respect to the macro design variable ρMA can be calculated via Eqs. (24) and (25). The OC algorithm given in Eq. (26) is applied to update the macro design variable ρMA. The macro material distribution composed of multiple predefined materials (i.e., multiple microstructures to be topologically optimized) cannot be obtained until the convergent condition is satisfied. Then, the concurrent optimization at both scales will be conducted. In the microscale optimization, the effective elasticity tensors DH of multiple microstructures are calculated by using the homogenization method (Eq. (18)) with the induced displacement field uMI obtained by the micro FE analysis (Eq. (19)). Based on the macro displacement field U and the macro material distribution ρMA from the above macroscale optimization, the sensitivity information is calculated by using Eqs. (41) and (42). The OC algorithm given in Eqs. (45) and (47) is employed to update the micro design variable αMI. In the macroscale optimization, the macro FE analysis is implemented to obtain the macro displacement fields U and the macrostructural compliance C, with the new DH obtained by the microscale optimization in the last iteration. The sensitivities for the macroscale optimization can be calculated via Eqs. (7), (24), and (25). The OC algorithm given in Eq. (26) is applied to update the macro design variable ρMA. The concurrent optimization at microscale and macroscale is repeated until the overall stiffness of the macrostructure is maximized. As discussed in [25,65,66], it is necessary to tackle the connectivity issue between the adjacent microstructure to produce a manufacturable design when the macrostructure is featured with various microstructures. In this study, a kinematically connective constraint approach [25], as depicted in Fig. 5, is employed to guarantee the connectivity of the adjacent microstructures. The predefined connectors can be achieved by setting the non-design regions at the same positions within the design domain of all the microstructures. Note that these predefined connectors will more or less limit the design space in the concurrent optimization, and thus, they may slightly compromise the performance of the optimized design. However, in engineering, it is acceptable to reasonably sacrifice structural performance in order to achieve a manufacturable design.

Different types of materials

Corresponding material density ρ MA

3 6 9 14

0.20, 0.50, 0.80 0.14, 0.28, 0.42, 0.56, 0.7, 0.84 0.10, 0.20, 0.30, 0.40, 0.50, 0.60, 0.70, 0.80, 0.90 0.07, 0.13, 0.20, 0.26, 0.33, 0.40, 0.46, 0.54, 0.60, 0.67, 0.73, 0.80, 0.87, 0.93 0.05, 0.10, 0.15, 0.2, 0.25, 0.3, 0.35, 0.4, 0.45, 0.5, 0.55, 0.6, 0.65, 0.7, 0.75, 0.8, 0.85, 0.9, 0.95

19

concurrent optimization problems. For all the examples, the materials are subject to plane stress conditions. The Young’s modulus of the base material is E0 = 10 , and the Poisson’s ratio is μ = 0.3. The microstructures are discretized into 50 × 50 = 2500 four-node elements. The optimization will terminate when the difference of the objective function values between two successive iterations is lower than 10−4, or a maximum of 200 iteration steps of the concurrent optimization is reached. To provide a fair comparison, the macrostructure compliance of the multiscale concurrent design is directly calculated by using the effective properties of the optimized representative microstructures.

5.1. MBB beam with a concentrated load The MBB beam with length L = 150 and height H = 30 is shown in Fig. 6. It is loaded with a concentrated vertical force of F = 5 at the center of the top edge and supported on rollers at the bottom right corner and on fixed supports at the bottom left corner. A mesh with 150 × 30 = 4500 four-node elements is applied to discretize the macro design domain. The objective function is to minimize the macrostructure compliance under a global volume constraint of 40%. Here, five cases are discussed to determine the available number of representative microstructures employed to concurrently design the macrostructure, where the macrostructure is assumed to be composed of 3, 6, 9, 14, and 19 representative microstructures, respectively. As mentioned above, each representative microstructure is regarded as a type of material. The types of materials and the corresponding relative density of each case is listed in Table 1. As shown in Fig. 7, the macrostructure compliance of the multiscale concurrent design with 3 representative microstructures is 214.3908. The optimized macrostructural topology is depicted in Fig. 7(a), and the color1 bar in Fig. 7(b) shows the relative density of each material within the macrostructural topology (i.e., volume fraction of each representative microstructure), and the representative microstructures corresponding to each material are placed above the color bar in ascending order of volume fraction from left to right in Fig. 7(b). According to the microstructure distribution from the final macrostructural topology in Fig. 7(a), the result of the multiscale concurrent design shown in Fig. 7(c) is assembled with three types of representative microstructures in Fig. 7(b). Similar to Fig. 7, Figs. 8–11 show the result of the multiscale concurrent design with macrostructure

5. Numerical examples In this section, three numerical examples are presented to demonstrate the effectiveness of the proposed method for multiscale

1 For interpretation of color in Fig. 7, the reader is referred to the web version of this article.

81

Computational Materials Science 155 (2018) 74–91

Y. Zhang et al.

(b) Three representative microstructures with different volume fractions

(a) Optimized macrostructural topology

(c) Multiscale concurrent design of MBB Fig. 7. Three representative microstructures with compliance of 214.3908.

(b) Six representative microstructures with different volume fractions

(a) Optimized macrostructural topology

(c) Multiscale concurrent design of MBB Fig. 8. Six representative microstructures with compliance of 208.8187.

1

(b) Nine representative microstructures with different volume fractions

(a) Optimized macrostructural topology

(c) Multiscale concurrent design of MBB Fig. 9. Nine representative microstructures with compliance of 203.8208.

multiscale concurrent optimization to obtain a superior performance of the macrostructure at an affordable computational cost. Thus, in the subsequent example, unless otherwise specified, 9 representative microstructures are used in the multiscale concurrent optimization of the macrostructure. To fully demonstrate the proposed method, the result of MBB with global volume constraints of 50% is also presented. The final macrostructural topology is shown in Fig. 12(a), and the result of the multiscale concurrent optimization assembled with 9 representative microstructures is depicted in Fig. 12(b). The iterative histories of the objective function and global volume constraint are shown in Fig. 13. The iterative histories of the local volume constraints of each

compliance values of 208.8187, 203.8208, 203.0075, and 202.4900, which are composed of 6, 9, 14, and 19 representative microstructures, respectively. Obviously, a more excellent macrostructural performance can be obtained with more representative microstructures. It is easily understood that more representative microstructures can enable a broader design space. However, too many representative microstructures will lead to more microstructures being topologically optimized, significantly increasing the computational cost. From the above result of the multiscale concurrent design, it can also be observed that the macrostructural compliance decreases slightly when the number of representative microstructures is greater than 9. Thus, it is advisable that 9 representative microstructures are applied to implement the 82

Computational Materials Science 155 (2018) 74–91

Y. Zhang et al.

(b) Fourteen representative microstructures with different volume fractions

(a) Optimized macrostructural topology

(c) Multiscale concurrent design of MBB Fig. 10. Fourteen representative microstructures with compliance of 203.0075.

(b) Nineteen representative microstructures with different volume fractions

(a) Optimized macrostructural topology

(c) Multiscale concurrent design of MBB Fig. 11. Nineteen representative microstructures with compliance of 202.4900.

(a) Optimized macrostructural topology

(b) Multiscale concurrent design with 9 representative microstructures Fig. 12. Compliance is 160.8174 for MBB with 50% global volume fraction.

properties as indicated in Table 2, which may provide better performances than the microstructure with isotropic properties with regard to the compliance optimization problem. The multiscale concurrent optimization can flexibly supply different directional stiffness with different configurations of representative microstructures. All the representative microstructures within the entire macrostructure are well connected by the setting of predefined connectors in different microstructures as shown in Fig. 3. Fig. 13 shows the iterative process of the multiscale concurrent optimization. At the first 123 iterative steps, the distinct

representative microstructure are plotted in Fig. 14. The configurations, properties, and percentages of these representative microstructures are provided in Table 2 and Fig. 15, respectively. It can be observed from Fig. 12 and Figs. 7–11 that the external frame regions of the multiscale concurrent design or the main path of force transfer is occupied by a number of completely solid or highdensity microstructure to effectively bear loading and resist deformation, which accords with the well-accepted design results for the MBB beam. All the representative microstructures show orthotropic 83

Computational Materials Science 155 (2018) 74–91

Y. Zhang et al.

Fig. 13. Iterative histories of the objective function and global volume constraint.

Fig. 14. Iterative histories of local volume constraint for each representative microstructure.

results using SIMP, under the corresponding global volume constraints of 50% and 40%. The monoscale macrostructural optimization is depicted in Fig. 16. It is obvious that the multiscale concurrent optimization has a lower macrostructural compliance (superior macrostructural performance) than the monoscale macrostructural optimization. The reason is that the proposed multiscale concurrent optimization method can greatly extend the design space by simultaneously optimizing the macrostructural topology, as well as the microstructural configurations and their corresponding distributions within the macrostructure.

material distribution at macroscale is achieved by using the multi-material topology optimization algorithm integrated with the ordered SIMP interpolation model. The subsequent 97 iterative steps show the iterative process of the simultaneous optimization of the macrostructural topology as well as the configurations and distributions of the microstructures. It can be observed from Figs. 13 and 14 that the global volume constraint is conserved all the time, whereas the local volume constraints of all representative microstructures are satisfied gradually during the optimization process. The violation of local volume constraints accounts for the increase in the objective function at the first stages of the simultaneous optimization. When the local volume constraints of the representative microstructures are satisfied, the macrostructural compliance starts to reduce gradually. Normally, it is well known that the monoscale macrostructural optimization can perform better than both the monoscale microstructural optimization and the multiscale concurrent optimization with uniform microstructures in most compliance optimization problems [41,48]. To show the benefits of the proposed multiscale concurrent optimization method, the optimized results composed of 9 representative microstructures are directly compared with the monoscale macrostructural

5.2. Long cantilever beam with multiple loads As depicted in Fig. 17, the long cantilever beam with length L = 100 and height H = 40 is evenly loaded with five vertical forces of F = 2.5 at the top edge, and all the degrees of freedom at the left edge are constrained. The design domain of macrostructure is discretized into 100 × 40 = 4000 four-node elements. The objective function is to minimize the macrostructural compliance under the global volume fraction constraints. In this example, four different types of global 84

Computational Materials Science 155 (2018) 74–91

Y. Zhang et al.

Table 2 Configurations and properties of each representative microstructure. Volume fraction

Microstructure

Level set surface

Effective elasticity tensor DH

0.1

0 ⎤ ⎡ 0.2757 0.2670 0 ⎥ ⎢ 0.2670 0.2889 0 0.2533 ⎦ ⎣ 0

0.2

0 ⎤ ⎡ 0.5992 0.5689 0 ⎥ ⎢ 0.5689 0.5993 0 0.5251⎦ ⎣ 0

0.3

0 ⎤ ⎡ 0.9785 0.8845 0 ⎥ ⎢ 0.8845 0.9842 0 0.8003 ⎦ ⎣ 0

0.4

0 ⎤ ⎡1.4730 1.2298 0 ⎥ ⎢1.2298 1.4909 0 1.0896 ⎦ ⎣ 0

0.5

0 ⎤ ⎡1.8247 1.1799 0 ⎥ ⎢1.1799 3.0593 0 1.1793 ⎥ ⎢ 0 ⎣ ⎦

0.6

0 ⎤ ⎡ 4.0987 1.2397 0 ⎥ ⎢1.2397 3.1469 0 1.2811⎦ ⎣ 0

0.7

0 ⎤ ⎡5.6221 1.4457 0 ⎥ ⎢1.4457 4.2221 0 1.5664 ⎦ ⎣ 0

0.8

0 ⎤ ⎡ 7.0202 1.8705 0 ⎥ ⎢1.8705 5.8313 0 0 2.0481 ⎣ ⎦

0.9

0 ⎤ ⎡ 9.1309 2.2224 0 ⎥ ⎢ 2.2224 6.8301 0 2.6618 ⎦ ⎣ 0

As depicted in Figs. 18–21, the macrostructure compliance values of the multiscale concurrent design under the global volume constraints of 60%, 40%, 20%, and 10%, are 392.0606, 608.6006, 1290.1458, and 2944.8087, respectively. As shown in Fig. 22, the macrostructural compliance values of the monoscale macrostructural optimization are 402.4065, 690.0025, 1852.4094, and 8452.2462 under the

volume fraction constraints are chosen, namely 60%, 40%, 20%, and 10%, and 9 representative microstructures are used. In order to demonstrate the performance of the proposed multiscale concurrent optimization method, the results of monoscale macrostructural optimization under the same global volume fraction constraints are provided for comparison.

85

Computational Materials Science 155 (2018) 74–91

Y. Zhang et al.

Fig. 15. Percentage of each representative microstructure.

More importantly, it can be found from the above results that the multiscale concurrent optimization is much better than the monoscale macrostructural optimization when the global volume fraction constraint is limited below 20%. This is because the monoscale macrostructural optimization cannot easily obtain a clear “black-white” solution under the constraints of a relatively small volume fraction, as shown in Fig. 22(c) and (d). By contrast, the multiscale concurrent optimization can achieve an optimized result with clear and distinct structural topologies by fully utilizing the information of microstructures at microscale and incorporating the interaction at both macrostructure and microstructures. It can be observed from Figs. 18–21 that the solid materials or the higher-density microstructures are mainly distributed at the main path of force transfer or the upper and lower edges to more effectively transfer loads and resist deformation for a long cantilever loaded with multiple forces. Hence, the proposed multiscale concurrent optimization method can significantly improve the macrostructural performance by incorporating the interaction at both scales, especially when the global volume fraction constraint (i.e., the overall material usage) is at a low level.

(a) Compliance is 216.5226 with global volume fraction of 40%

(b) Compliance is 169.2639 with global volume fraction of 50% Fig. 16. Monoscale macrostructure designs of MBB beam using SIMP.

F

F

Design Domain

F

F

H

F

5.3. Short cantilever beam with multiple cases

L

A short cantilever beam with length L = 80 and height H = 60 with multiple cases is shown in Fig. 23. The vertical force F1 = 5 is loaded at the bottom right corner, and the vertical force F2 = −5 is loaded at the top right corner. The left edge of the short cantilever beam is fixed. A mesh with 80 × 60 = 4800 four-node elements is applied to discretize the design domain. The objection function is to minimize the macrostructural compliance under various global volume fraction constraints, i.e., 40%, 20%, and 10%, and 9 representative microstructures are also used in this example. For comparison, the SIMP-based monoscale macrostructural optimization with the same volume constraints is also given as follows. Figs. 24–26 show the multiscale concurrent optimization under the different global volume fraction constraints of 40%, 20%, and 10%. Their corresponding macrostructural compliance values are 253.8597, 528.4621, and 1454.0011, respectively. All the representative

Fig. 17. Sketch of cantilever beam with multiple loads.

corresponding global volume fraction constraints. All the representative microstructures to be topologically optimized with the proposed method are well connected due to the predefined connectors as shown in Fig. 3. It is obvious that the compliance values of the multiscale concurrent optimization in all numerical examples are lower than those using the monoscale macrostructural optimization. It demonstrates again that the proposed multiscale concurrent optimization method can fully make use of different microstructures to broaden the design space and sufficiently incorporate the interaction at both scales to achieve a superior macrostructural performance.

86

Computational Materials Science 155 (2018) 74–91

Y. Zhang et al.

(b) Nine representative microstructures with different volume fractions

(a) Optimized macrostructural topology

(c) Multiscale concurrent design of cantilever with multiple loads Fig. 18. Compliance is 392.0606 with 60% global volume fraction.

(b) Nine representative microstructures with different volume fractions

(a) Optimized macrostructural topology

(c) Multiscale concurrent design of cantilever with multiple loads Fig. 19. Compliance is 608.6006 with 40% global volume fraction.

concurrent optimization can produce much stiffer structures than the monoscale macrostructural optimization, especially when the volume fraction constraints are relatively small, such as under 20%. As shown in Figs. 24–26, the higher-density elements or solid material elements are mainly distributed at the upper and lower edges of the macro design domain, and the lower-density elements are placed at the middle region. It is suitable for the short cantilever shown in Fig. 23, which can more effectively resist the deformation caused by the external

microstructures within the entire macro design domain maintain a good connectivity in geometry by the setting of the predefined connector in each microstructure. They also reveal the remarkable orthotropic properties to maximize the macrostructural performance. Fig. 27 shows the monoscale macrostructural optimization under the same volume fraction constraints, with the corresponding macrostructural compliance of 265.2747, 650.6535, and 2767.2622. In accordance with the result of the long cantilever beam with multiple loads, the multiscale

87

Computational Materials Science 155 (2018) 74–91

Y. Zhang et al.

(b) Nine representative microstructures with different volume fractions

(a) Optimized macrostructural topology

(c) Multiscale concurrent design of cantilever with multiple loads Fig. 20. Compliance is 1290.1458 with 20% global volume fraction.

Fig. 21. Compliance is 2944.8087 with 10% global volume fraction.

6. Conclusions

loads.More interestingly, under the relatively small volume fraction constraints (i.e., 20% or 10%), the macrostructure of the proposed multiscale concurrent optimization is more stable than that provided by the monoscale macrostructural optimization, when the external load F1 or F2 slightly slopes on the vertical direction. This also reveals the excellent capacity of the proposed multiscale concurrent optimization for multiple cases, especially when the global volume fraction constraint is at a relatively low level.

This paper has presented a novel multiscale concurrent topology optimization method for cellular structures with multiple types of microstructures. In this method, the macrostructural topology, as well as the locations and configurations of microstructures are simultaneously optimized using a unified optimization procedure. At the macroscale, the multi-material topology optimization algorithm integrated with the

88

Computational Materials Science 155 (2018) 74–91

Y. Zhang et al.

(a) Compliance is 402.4065 with 60% global volume fraction

(b) Compliance is 690.0025 with 40% global volume fraction

(c) Compliance is 1852.4094 with 20% global volume fraction

(d) Compliance is 8452.2462 with 10% global volume fraction

Fig. 22. Monoscale macrostructure designs of cantilever with multiple loads.

F2

L

H

Design Domain

the macroscale element is regarded as the volume constraint of the microstructural design. The number of representative microstructures is same as the types of predefined multi-materials. All the representative microstructures are topologically optimized by integrating the PLSM with the numerical homogenization method. Each representative microstructure is designed for all microstructures with the same volume fraction. Thus, the total number of microstructures in the macrostructure is significantly reduced, leading to a considerable reduction in the computational cost of the multiscale concurrent optimization. The numerical examples reveal the excellent capacity of the proposed multiscale concurrent optimization in improving the macrostructural performance compared to the monoscale macrostructural optimization, especially when the overall material usage is at a low level. The future work is to extend the proposed multiscale concurrent optimization method to other application fields, such as dynamics, acoustics, and functionally graded and multifunctional materials.

F1

Fig. 23. Sketch of short cantilever beam with multiple cases.

ordered SIMP interpolation model is used to achieve the optimized macrostructural topology and optimal locations of microstructures. At the microscale, all macroscale elements with the same material are represented by a unique microstructure. The relative element density of

(a) Optimized macrostructural topology

(b) Nine representative microstructures with different volume fractions

(c) Multiscale concurrent design of cantilever with multiple cases

Fig. 24. Compliance is 253.8597 with 40% global volume fraction.

89

Computational Materials Science 155 (2018) 74–91

Y. Zhang et al.

(a) Optimized macrostructural topology

(b) Nine representative microstructures with different volume fractions

(c) Multiscale concurrent design of cantilever with multiple cases

Fig. 25. Compliance is 528.4621 with 20% global volume fraction.

(a) Optimized macrostructural topology

(b) Nine representative microstructures with different volume fractions

(c) Multiscale concurrent design of cantilever with multiple cases

Fig. 26. Compliance is 1454.0011 with 10% global volume fraction.

(a) Compliance is 265.2747 with 40% global volume fraction

(b) Compliance is 650.6535 with 20% global volume fraction

(c) Compliance is 2767.2622 with 10% global volume fraction

Fig. 27. Monoscale macrostructure design of cantilever with multiple cases.

7. Authors contributions

Acknowledgments

The first author Yan Zhang conceived the central idea, analysed most of the data, and wrote the initial draft of the paper. The corresponding author Liang Gao contributed to refining the ideas and revising the manuscript. All the remaining authors discussed the data and were involved in writing and revising the manuscript.

This work was supported by the National Basic Scientific Research Program of China [grant number JCKY2016110C012], and the National Natural Science Foundation of China [grant numbers 51705166, 51675196 and 51421062], and the Program for HUST Academic Frontier Youth Team.

90

Computational Materials Science 155 (2018) 74–91

Y. Zhang et al.

References [1] [2] [3] [4] [5] [6] [7] [8] [9] [10] [11] [12] [13] [14] [15] [16] [17] [18] [19] [20] [21] [22] [23] [24] [25] [26] [27] [28] [29] [30] [31] [32] [33]

[34] [35] [36] [37]

R.M. Christensen, Int. J. Solids Struct. 37 (2000) 93–104. S.C. Han, J.W. Lee, K. Kang, Adv. Mater. 27 (2015) 5506–5511. Y. Wang, Z. Luo, N. Zhang, Z. Kang, Comput. Mater. Sci. 87 (2014) 178–186. H. Li, Z. Luo, N. Zhang, L. Gao, T. Brown, Comput. Methods Appl. Mech. Eng. 309 (2016) 453–475. H. Li, Z. Luo, L. Gao, P. Walker, Comput. Methods Appl. Mech. Eng. 328 (2018) 340–364. G. Rozvany, Struct. Multidiscip. Optim. 21 (2001) 90–108. M.P. Bendsøe, N. Kikuchi, Comput. Methods Appl. Mech. Eng. 71 (1988) 197–224. M.P. Bendsoe, O. Sigmund, Topology Optimization: Theory, Methods, and Applications, Springer, Berlin, Heidelberg, 2003. O. Sigmund, K. Maute, Struct. Multidiscip. Optim. 48 (2013) 1031–1055. G.I. Rozvany, M. Zhou, T. Birker, Struct. Optim. 4 (1992) 250–252. G. Allaire, Shape Optimization by the Homogenization Method, Springer Science & Business Media, 2012. M.P. Bendsøe, O. Sigmund, Arch. Appl. Mech. 69 (1999) 635–654. O. Sigmund, Struct. Multidiscip. Optim. 21 (2001) 120–127. Y.M. Xie, G.P. Steven, Comput. Struct. 49 (1993) 885–896. X. Huang, Y. Xie, Comput. Mech. 43 (2009) 393. M.Y. Wang, X. Wang, D. Guo, Comput. Methods Appl. Mech. Eng. 192 (2003) 227–246. G. Allaire, F. Jouve, A.-M. Toader, J. Comput. Phys. 194 (2004) 363–393. E. Andreassen, B.S. Lazarov, O. Sigmund, Mech. Mater. 69 (2014) 1–10. X. Huang, A. Radman, Y. Xie, Comput. Mater. Sci. 50 (2011) 1861–1870. J. Gao, H. Li, L. Gao, M. Xiao, Adv. Eng. Softw. 116 (2018) 89–102. O. Sigmund, S. Torquato, J. Mech. Phys. Solids 45 (1997) 1037–1067. O. Sigmund, S. Torquato, Appl. Phys. Lett. 69 (1996) 3203–3205. J.K. Guest, J.H. Prévost, Int. J. Solids Struct. 43 (2006) 7028–7047. N. de Kruijf, S. Zhou, Q. Li, Y.-W. Mai, Int. J. Solids Struct. 44 (2007) 7092–7109. S. Zhou, Q. Li, J. Mater. Sci. 43 (2008) 5157–5167. A. Radman, X. Huang, Y.M. Xie, J. Mater. Sci. 48 (2012) 1503–1510. X. Huang, S. Zhou, G. Sun, G. Li, Y.M. Xie, Comput. Methods Appl. Mech. Eng. 283 (2015) 503–516. E. Andreassen, J.S. Jensen, Struct. Multidiscip. Optim. 49 (2014) 695–705. X. Huang, S. Zhou, Y. Xie, Q. Li, Comput. Mater. Sci. 67 (2013) 397–407. W. Su, S. Liu, Struct. Multidiscip. Optim. 42 (2010) 243–254. H. Rodrigues, J.M. Guedes, M.P. Bendsoe, Struct. Multidiscip. Optim. 24 (2002) 1–10. P. Coelho, H. Rodrigues, Struct. Multidiscip. Optim. 52 (2015) 91–104. P. Coelho, P. Fernandes, J. Guedes, H. Rodrigues, Struct. Multidiscip. Optim. 35

[38] [39] [40] [41] [42] [43] [44] [45] [46] [47] [48] [49] [50] [51]

[52] [53] [54] [55] [56] [57] [58] [59] [60] [61] [62] [63] [64] [65] [66]

91

(2008) 107–115. F. Schury, M. Stingl, F. Wein, SIAM J. Scient. Comput. 34 (2012) B711–B733. L. Xia, P. Breitkopf, Comput. Methods Appl. Mech. Eng. 278 (2014) 524–542. L. Xia, P. Breitkopf, Comput. Methods Appl. Mech. Eng. 286 (2015) 147–167. R. Sivapuram, P.D. Dunning, H.A. Kim, Struct. Multidiscip. Optim. 54 (2016) 1267–1281. Y. Wang, F. Chen, M.Y. Wang, Comput. Methods Appl. Mech. Eng. 317 (2017) 84–101. D. Fujii, B. Chen, N. Kikuchi, Int. J. Numer. Meth. Eng. 50 (2001) 2031–2051. L. Liu, J. Yan, G. Cheng, Comput. Struct. 86 (2008) 1417–1425. Y. Wang, M.Y. Wang, F. Chen, Struct. Multidiscip. Optim. 54 (2016) 1145–1156. W. Chen, L. Tong, S. Liu, Comput. Struct. 178 (2017) 119–128. X. Yan, X. Huang, G. Sun, Y.M. Xie, Compos. Struct. 120 (2015) 358–365. K. Long, P.F. Yuan, S. Xu, Y.M. Xie, Eng. Optim. 50 (2017) 599–614. K. Long, D. Han, X. Gu, Comput. Mater. Sci. 129 (2017) 194–201. W. Zhang, S. Sun, Int. J. Numer. Meth. Eng. 68 (2006) 993–1011. L. Xu, G. Cheng, Struct. Multidiscip. Optim. 57 (2018) 2093–2107. H. Li, Z. Luo, L. Gao, Q. Qin, Comput. Methods Appl. Mech. Eng. 331 (2018) 536–561. J. Deng, W. Chen, Struct. Multidiscip. Optim. 56 (2017) 1–19. O. Sigmund, Comput. Methods Appl. Mech. Eng. 190 (2001) 6605–6627. Y. Luo, Z. Kang, Layout Design of Reinforced Concrete Structures Using TwoMaterial Topology Optimization with Drucker-Prager Yield Constraints, SpringerVerlag, New York Inc, 2013. M.Y. Wang, S. Zhou, Comput. Methods Appl. Mech. Eng. 193 (2004) 469–496. Y. Wang, Z. Luo, Z. Kang, N. Zhang, Comput. Methods Appl. Mech. Eng. 283 (2015) 1570–1586. H. Zhang, Y. Luo, Z. Kang, Compos. Struct. (2018). Z. Kang, Y. Wang, Comput. Methods Appl. Mech. Eng. 200 (2011) 3515–3525. W. Zuo, K. Saitou, Struct. Multidiscip. Optim. 55 (2016) 477–491. E. Andreassen, C.S. Andreasen, Comput. Mater. Sci. 83 (2014) 488–495. O. Sigmund, N. Aage, E. Andreassen, Struct. Multidiscip. Optim. 54 (2016) 361–373. H. Wendland, Adv. Comput. Math. 4 (1995) 389–396. Z. Luo, M.Y. Wang, S. Wang, P. Wei, Int. J. Numer. Meth. Eng. 76 (2008) 1–26. Z. Luo, L. Tong, M.Y. Wang, S. Wang, J. Comput. Phys. 227 (2007) 680–705. M. Zhou, G. Rozvany, Comput. Methods Appl. Mech. Eng. 89 (1991) 309–336. K. Svanberg, Int. J. Numer. Meth. Eng. 24 (1987) 359–373. Z. Luo, L. Tong, Z. Kang, Comput. Struct. 87 (2009) 425–434. J.E. Cadman, S. Zhou, Y. Chen, Q. Li, J. Mater. Sci. 48 (2013) 51–66. J. Alexandersen, B.S. Lazarov, Comput. Methods Appl. Mech. Eng. 290 (2015) 156–182.