Level set-based heterogeneous object modeling and optimization

Accepted Manuscript Level set-based heterogeneous object modeling and optimization Jikai Liu, Qian Chen, Yufan Zheng, Rafiq Ahmad, Jinyuan Tang, Yongsheng Ma

PII: DOI: Reference:

S0010-4485(18)30414-7 https://doi.org/10.1016/j.cad.2019.01.002 JCAD 2663

To appear in:

Computer-Aided Design

Received date : 3 August 2018 Accepted date : 6 January 2019 Please cite this article as: J. Liu, Q. Chen, Y. Zheng et al. Level set-based heterogeneous object modeling and optimization. Computer-Aided Design (2019), https://doi.org/10.1016/j.cad.2019.01.002 This is a PDF file of an unedited manuscript that has been accepted for publication. As a service to our customers we are providing this early version of the manuscript. The manuscript will undergo copyediting, typesetting, and review of the resulting proof before it is published in its final form. Please note that during the production process errors may be discovered which could affect the content, and all legal disclaimers that apply to the journal pertain.

Level set-based heterogeneous object modeling and optimization Jikai Liu1,2, Qian Chen3, Yufan Zheng5, Rafiq Ahmad5, Jinyuan Tang4*, Yongsheng Ma5* 1

Center for Advanced Jet Engineering Technologies (CaJET), School of Mechanical Engineering, Shandong University, Jinan, China

2

Key Laboratory of High Efficiency and Clean Mechanical Manufacture, Shandong University, Ministry of Education, China

3

Department of Mechanical Engineering and Materials Science, University of Pittsburgh, Pittsburgh, PA 15261, United States 4

State Key Laboratory of High-Performance Complex Manufacturing, Central South University, Changsha, China 5

Department of Mechanical Engineering, University of Alberta, Edmonton, Canada

*Corresponding authors. Emails: [email protected]; [email protected]

Abstract This paper presents a level set-based heterogeneous object (HO) modeling and optimization method. This HO model employs multiple level set functions to build the geometry, utilizes zero-value level set contours as material source profiles, and realizes functionally graded material blending with a signed distance-based blending function. More importantly, this HO model supports the concurrent structure and material optimization because of the unified level set framework for both structure and material composition representation. Beyond macro HO, heterogeneous meta-material optimization will be addressed as well. This new model remedies the shortage of traditional HO models that focus more on modeling but less on optimization. About the numerical optimization, design update with the sensitivity result will be carefully discussed, since there includes infeasible terms (in domain integration format). Two strategies will be explored to address this issue: ignoring the infeasible part of the sensitivity, or transforming the sensitivity result into a purely boundary integration-based expression. A few numerical examples will be studied to prove the effectiveness of the proposed HO modeling and optimization method. Keywords: Heterogeneous object modeling, Heterogeneous object optimization; Level set; Functionally graded material; Sensitivity analysis

1. Introduction 1.1 Heterogeneous object modeling Heterogeneous object (HO) modeling has become an important research area in recent years because of the extensive use of composite materials and the advancement of heterogeneous manufacturing process, e.g. 3D printing [1]. Extensive research works have been conducted in the past two decades and several HO modeling methods have been proposed. A literature survey can be found in [2]. Among the proposed methods, voxel based model, cellular model, and control feature based model are the most popular and a brief literature survey will be presented in the following paragraphs.

1

Voxel-based model [3–6] discretizes the geometry into voxels and defines material compositions of each voxel. Therefore, voxel-based model could be directly used for finite element analysis and has the potential to capture the highly irregular heterogeneous material composition distribution. Cellular model is an alternative of the geometry discretization which belongs to a larger scale. Kumar and Dutta [7] and Kumar et al. [8] applied the r-set to modeling the geometry and each divided sub-regions was assigned a certain material class. Bhashyam et al. [9] developed a prototype system with graphical user interface, material composition function library, and property estimation rule library, which is effective and efficient in designing and analyzing HOs. Later, Shin and Dutta [10] and Shin et al. [11] presented a constructive representation of HO model by manipulating the heterogeneous sub-regions with heterogeneous Boolean operations. For cellular model, the material composition function is necessary for each individual cell. Control feature-based model developed by Siu and Tan [12] employs source profile features as references to determine the functionally graded material (FGM) distribution with the source profile-based material composition functions. This method was also applied to composite laminates, where the continuously varying material compositions, fiber aspect ratios, and fiber orientations were concurrently controlled by source profile features and material composite functions [13]. Later, heterogeneous feature tree was proposed for complex HO modeling, which included multiple levels of source profile features and the non-regular Boolean operations [14,15]. In addition, the authors integrated CAE tools with the HO modeling system to perform HO design validation and modification [16]. Liu et al. [17] realized the local composition control by employing single-control feature, multi-control feature, and Laplace equation based approaches. Biswas et al. [18] and Xu and Shaw [19] both employed the distance field for HO modeling which in nature is an extension of the control feature-based method by using the implicit source profiles. Such distance field-based methods enable HO modeling based on freeform geometry and even point cloud [20]. Recently, Gupta and Tandon [21] developed a convolution material primitive-based method which highlighted the compound material blending. B-spline-based model can be regarded as an extension of control feature-based model. Qian and Dutta [22,23] developed a diffusion-based method to determine the material composition distribution based on B-spline source profiles. This method was also applied to designing heterogeneous turbine blade [24]. Yang and Qian [25] employed B-spline finite element method to unify the design and analysis model; the heterogeneous lofting algorithms were used to determine the material composition distribution between B-spline source profiles. Samanta and Koc [26] applied B-spline for complex geometry and material feature modeling and an optimization process was employed to enhance the HO generation. Recently, Ozbolat nad Koc [27] and Samanta et al. [28] enhanced the method by realizing material blending between freeform B-spline curves via ruling line generation. Other HO models like radial basis function-based model [29] and stochastic Vironoi diagram representation [30,31] are also powerful under certain scenarios. 1.2 Heterogeneous object optimization Practically for a design task, it is expected to not only build the HO model, but also optimize it to derive the optimal solution. However, this aspect attracts much less attention. Chen and Feng [32,33] performed HO optimization guided by the Axiomatic Design method. Both geometry and material information of the HO were created through optimization. However, a clear limitation exists in this method that: it procedurally optimized the geometry and then the material, which lacks overall optimality. In a recent work, Kou et al. [34] applied particle swarm optimization to optimally distributing FGMs for a well-defined geometry.

2

Level set model [35,36] is relatively a new HO representation. To be specific, each level set function is a closed contour for material/material (or material/void) interface representation. Because of overlapping, n level set functions can represent m material phases, and the relationship between n and m depends on the specific multi-level set interpolation. This capacity makes level set-based HO model suitable for multimaterial structural representation and optimization [37–42]. Level set method has also been used for coating structure design [43,44] Additionally, level set method was applied to HO modeling and optimization with FGMs [45,6,46]. A prominent characteristic of level set method is that it provides a unified mathematical framework for HO modeling and optimization that well supports concurrent structural and material composition optimization. There are also some other concurrent shape/topology and material optimization methods, such as the FGM-SIMP [47] and iso-geometric [48], which have realized similar design effects. 1.3 Research highlight Level set-based HO model is one of the a few exceptions that could support both aspects of modeling and optimization. However, the existing level set-based HO models generally have the material information separately defined from the geometry (i.e. the level set function) [45,6,46], which can be categorized as a voxel-based model. The limitation of voxel-based model is the high memory requirement to store data [2] and the non-trivialness to manually manipulate the material distribution. For example, in [45], the material compositions were element-wisely defined and the graded material composition distribution was achieved by constraining the gradient of the material composition variable. In [6,46], the element-wisely defined material composition variables were just optimized without constraining the gradients. Therefore, the main motivation of this research is to develop a new level set-based HO model, wherein each level set function corresponds to a material source profile and the material compositions can be calculated based on the signed distance information. This apparently belongs to a control feature-based model. In this way, only knowing the individual level set contours (representing the different material source profiles) and the material composition function would be sufficient to store and build the HO model. This model is more compact and exact in data representation [2]. In addition, parameterized level set method can be used for CSG (Constructive Solid Geometry) type feature-based CAD modeling [49]. Therefore, the proposed method would also be friendly to feature-based HO modeling and optimization, for which manual manipulation of the HO model would be straightforward. On the other hand, the design space of the control feature-based HO model is restricted compared to the voxel-based model, because the material heterogeneity cannot be freely defined. Even so, the control feature-based HO model has its unique advantage to model structures that have the surface-to-surface material grading. For example, the natural bamboo structure has the fiber density increasing from the inside surface to the outside surface so that the structure can have uniform strength [50,51]; see Fig. 1. Some aerospace structures have the material compositions from pure ceramic of the outside surface gradually to the pure metal of the inner surface, so that the outside surface has lower coefficients of thermal expansion and high thermal/corrosion protection, and the inside surface provides the load carrying capability [52].

3

Fig. 1 Cross-section of the natural bamboo structure [51] Then, from the algorithm aspect, the proposed HO optimization method inherits the Ersatz material approach in handling the boundary elements but the formulated problem is more complex. Different from [45,6,46], the material constitutive model in the newly formulated problem is a function of the level set functions. Therefore, there involves several terms in the sensitivity result, including terms in domain integral format, which greatly complicates the problem solution since update of the discrete level set field is only performed at a narrow band around the boundary area. Therefore, new techniques are proposed to process the sensitivity result so that design update can be correctly performed. 2. Heterogeneous object modeling with level set method 2.1 Level set based geometric modeling For geometric modeling, B-rep (Boundary representation) is widely adopted by commercial CAD systems, which represents the CAD geometry by explicit geometric entities, e.g. faces, edges, and vertices, and their topological relationship. However, in this work, CSG (Constructive solid geometry) is adopted, because level set function represents the geometry in an implicit form [53] which is insensitive to topological changes [49]. This characteristic greatly facilitates the structural topology optimization. To be specific, the level set function, Φ(𝐗): 𝑅 𝑛 ⟼ 𝑅, defines any freeform geometry in the implicit form, as demonstrated in Eq. (1). Φ(𝐗) > 0, 𝐗 ∈ Ω/ ∂Ω {Φ(𝐗) = 0, 𝐗 ∈ ∂Ω Φ(𝐗) < 0, 𝐗 ∈ D/Ω

(1)

In implementation of Eq. (1), the level set function can be either discretely defined or parametrically defined through B-Spline [49] or Radial basis function [54,55]. Level set method also supports form feature-based modeling [49,56]. For instance, a circle feature can be represented by: Φ(𝐗) = 𝑅𝑐 − 𝑠𝑞𝑟𝑡((𝑥 − 𝑥0 )2 + (𝑦 − 𝑦0 )2 )

(2)

and a square feature by: 𝐻𝑠 2

Φ(𝐗) = 𝑚𝑖𝑛 {𝐻𝑠 2

𝐻𝑠 + [(𝑥 − 𝑥0 )𝑐𝑜𝑠𝜃 + (𝑦 − 𝑦0 )𝑠𝑖𝑛𝜃], 2 } 𝐻 𝑦0 )𝑐𝑜𝑠𝜃], 2𝑠 + [−(𝑥 − 𝑥0 )𝑠𝑖𝑛𝜃 + (𝑦 − 𝑦0 )𝑐𝑜𝑠𝜃]

− [(𝑥 − 𝑥0 )𝑐𝑜𝑠𝜃 + (𝑦 − 𝑦0 )𝑠𝑖𝑛𝜃],

− [−(𝑥 − 𝑥0 )𝑠𝑖𝑛𝜃 + (𝑦 −

(3)

where (𝑥0 , 𝑦0 ) is the center coordinate; 𝑅𝑐 is the circle radius; 𝐻𝑠 and 𝜃 are the square length and orientation, respectively. Then, with the individually represented geometry primitives, a complex geometry can be constructed through Boolean operations as: Φ1 ∪ Φ2 = 𝑚𝑎𝑥(Φ1 , Φ2 ) Φ1 ∩ Φ2 = 𝑚𝑖𝑛(Φ1 , Φ2 ) Φ1 \ Φ2 = 𝑚𝑖𝑛(Φ1 , −Φ2 )

(4)

Finally, the geometry is formulated as: 4

Φs (𝐗) = 𝐶(Φf1 , Φf2 , Φf3 , … )

(5)

where 𝐶 is the set of Boolean operations, and 𝑛 represents the number of geometry primitives. 2.2 Level set based material blending For material blending, the individual level set contours are used as material source profiles and a signed distance-based material blending function is proposed to determine the local material compositions. To be specific, an important characteristic of level set function is the signed distance information (i.e., Eq. (6) is satisfied). Then, the level set value at any point indicates its closest distance from the structural boundary (zero-value level set contour) and the sign indicates the point to be either solid (positive) or void (negative). An example of signed distance field is demonstrated in Fig. 2. (6)

|∇Φ(𝐗)| = 1

Fig. 2 Signed distance field Then, the material blending can be realized based on the signed distance information. Each level set function is treated as a material source profile. Since a signed distance field is associated with each level set contour, the signed distance-based material blending can be achieved. Specific form of the level setbased material blending function is demonstrated in Eq. (7). 𝑉𝑖 =

[1/|Φi (𝐗)|]

𝑘

𝑘 ∑𝑚 𝑗=1[1/|Φj (𝐗)|]

for any 𝐗 which satisfies 𝐻(Φ(𝐗)) > 0

(7)

Φ(𝐗) = 𝐶(Φ1 , Φ2 , … , Φm ) in which, Φi (𝐗) is the ith level set function value at point X and | | means the absolute. 𝑉𝑖 is the volume fraction of the material type 𝑖, which is associated to the material source profile Φi (𝐗). 𝑚 is the total number of level set functions used for HO modeling. It is worth noticing that an index 𝑘 is involved in the function to control the order of smoothness. For instance, Figure 3 plots the 1st order and 2nd order one dimensional material blending. Another point to clarify is that, the level set functions Φ𝑖 (𝐗) used for the 5

material blending are different from the level set functions Φ𝑓𝑖 used to construct the geometry in Eq. (5). Φ𝑓𝑖 means the a basic feature primitive, such as a hole or a rectangle. Differently, Φ𝑖 (𝐗) represents a material source profile. It can be a basic feature primitive as demonstrated in Eq. (2) and (3); more importantly, it can also be a high-level feature model that is the consequence of some Boolean operations. In this way, the material source profile can be more flexibly defined. This material blending rule is similar to the Shepard interpolation [57–60], which uses the inverse of distances as the weights to evaluate the material composition rates. Each weight falls in the range between 0 and 1, and the sum of the weights equals to 1. More importantly, the weights are differentiable on the level set function.

(a) Linear material blending

(b) Quadratic material blending

Fig. 3 Material blending in different orders Also in Eq. (7), 𝐻( ) is the Heaviside function which is defined as: 𝐻(Φ) = 0 Φ < 0 (8) { 𝐻(Φ) = 1 Φ ≥ 0 The approximate form of the Heaviside function as employed in [61] will be used in the later numerical implementations.

2.3 Examples for HO modeling This subsection presents some examples of the level set-based HO modeling. First, Figure 4 demonstrates a few HO models with two source profile features and different orders of smoothness. The rectangular frame and the interior hole are the two zero-value level set contours, i.e. the material source profiles. All results demonstrated in Fig. 4 have realized the smooth material transitions between the source profiles.

6

(a) Linear material blending

(b) Quadratic material blending

(c) Third-order material blending

(d) Fourth-order material blending

Fig. 4 Material blending with two material source profiles Then, Figure 5 presents a few HO models with three source profile features and different orders of smoothness. The rectangular frame and the two interior holes are the three material source profiles. In this case, we explored the different spatial positions of the interior holes. Especially for the results demonstrated in Fig. 5b and 5d, the two interior holes overlap and form a single heterogeneous material source profile, which indicates that the level set-based HO model can support heterogeneous source profile-based HO modeling.

7

(a) Linear blending without overlap

(b) Linear blending with overlap

(c) Quadratic blending without overlap

(d) Quadratic blending with overlap

Fig. 5 Material blending with three material source profiles The level set-based HO model also supports point and line sources, as demonstrated in Fig. 6.

(a) Linear blending with point source

(b) Linear blending with line source

Fig. 6 Material blending with point and line sources

3. Heterogeneous object optimization with level set method In this section, the level set-based HO model proposed in the last section will be optimized to extremize the structural performance. A highlight would be the concurrent structural and material optimization. 3.1 Problem definition Taking linear elastic structures for example, the material elasticity tensor is expressed by: 𝑚

𝐄𝐵 = ∑ 𝑖

[1/|Φ𝑖 (𝐗)|]

𝑘

𝑘 ∑𝑚 𝑗=1[1/|Φ𝑗 (𝐗)|]

𝐄𝑖

(9)

8

for any 𝐗 which satisfies 𝐻(Φ(𝐗)) > 0 Φ(𝐗) = 𝐶(Φ1 , Φ2 , … , Φm ) in which 𝐄𝑖 is the elasticity tensor of the ith material type. Here, it is worth noticing that, the local elastic properties are approximated through a simple linear mixture according to the material volume fractions. This approximation is widely used by FGM structure optimization implementations [34,46] for the sake of simplicity. Note that, this mixture rule does not satisfy the Hashin–Shtrikman bounds [62], and a more rigorous model of the material properties [45] should be followed in case that physical realization and validation of the FGM structure is targeted. Structural compliance is employed as the design objective and the related topology optimization problem is formulated as: min. 𝐽(𝐄B , 𝐮, Φ) = ∫ 𝐄B 𝐞(𝐮)𝐞(𝐮)𝐻(Φ)𝑑Ω 𝐷

𝑠. 𝑡. 𝑎(𝐄B , 𝐮, 𝐯, Φ) = 𝑙(𝐯, Φ),

∀𝐯∈𝑈

𝑉 = ∫ 𝐻(Φ)𝑑Ω ≤ 𝑉𝑚𝑎𝑥

(10)

𝐷

𝑎(𝐄B , 𝐮, 𝐯, Φ) = ∫ 𝐄B 𝐞(𝐮)𝐞(𝐯)𝐻(Φ)𝑑Ω 𝐷

𝑙(𝐯, Φ) = ∫ 𝐩𝐯𝐻(Φ)𝑑Ω + ∫ 𝛕𝐯𝛿(Φ)|∇Φ|𝑑Ω 𝐷

𝐷

where 𝐮 is the displacement vector, 𝐯 is the test vector, and 𝐞(𝐮) is the strain. 𝑈 = {𝐯 ∈ 𝐻1 (Ω)𝑑 |𝐯 = 0 𝑜𝑛 ΓD } is the space of kinematically admissible displacement field. 𝐩 is the body force and 𝛕 is the boundary traction force. 𝛿(Φ) is the Dirac delta function, which is defined as: 𝛿(Φ) = 0 Φ ≠ 0 { 𝛿(Φ) = +∞ Φ = 0

+∞

∫

𝛿(Φ)𝑑Φ = 1

(11)

−∞

3.2 Problem solution In order to solve the optimization problem, material derivative and the adjoint method are employed to perform the shape sensitivity analysis. The Lagrangian is formulated as, 𝐿 = 𝐽(𝐄B , 𝐮, Φ) + 𝑎(𝐄B , 𝐮, 𝐰, Φ) − 𝑙(𝐰, Φ) + 𝜆(∫ 𝐻(Φ)𝑑Ω − 𝑉𝑚𝑎𝑥 )

(12)

𝐷

in which 𝐰 is the adjoint variable. Material derivative of the Lagrangian is performed as, 𝐿′ = 𝐽′(𝐄B , 𝐮, Φ) + 𝑎′(𝐄B , 𝐮, 𝐰, Φ) − 𝑙′(𝐰, Φ) + 𝜆(∫ 𝐻(Φ)𝑑Ω)′

(13)

𝐷

9

where, 𝑚 ′(

𝐽 𝐄B , 𝐮, Φ) = ∑ ∫ 𝑖=1 𝐷

𝜕𝐄B 𝐞(𝐮)𝐞(𝐮)𝐻(Φ)𝑉𝑛 𝑖 |∇Φ𝑖 |𝑑Ω + 2 ∫ 𝐄B 𝐞(𝐮)𝐞(𝐮′)𝐻(Φ)𝑑Ω 𝜕Φ𝑖 𝐷 (14)

𝑚

+ ∑ ∫ 𝐄B 𝐞(𝐮)𝐞(𝐮)𝛿(Φ) 𝑖=1 𝐷 𝑚 ′(

𝑎 𝐄𝐵 , 𝐮, 𝐰, Φ) = ∑ ∫ 𝑖=1 𝐷

𝜕Φ 𝑖 𝑉 |∇Φ𝑖 |𝑑Ω 𝜕Φ𝑖 𝑛

𝜕𝐄B 𝐞(𝐮)𝐞(𝐰)𝐻(Φ)𝑉𝑛 𝑖 |∇Φ𝑖 |𝑑Ω + ∫ 𝐄B 𝐞(𝐮′)𝐞(𝐰)𝐻(Φ)𝑑Ω 𝜕Φ𝑖 𝐷 (15)

𝑚

𝜕Φ 𝑖 + ∫ 𝐄B 𝐞(𝐮)𝐞(𝐰′)𝐻(Φ)𝑑Ω + ∑ ∫ 𝐄B 𝐞(𝐮)𝐞(𝐰)𝛿(Φ) 𝑉 |∇Φ𝑖 |𝑑Ω 𝜕Φ𝑖 𝑛 𝐷 𝐷 𝑖=1

𝑙′(𝐰, Φ) = ∫ 𝛕𝐰′𝛿(Φ)|∇Φ|𝑑Ω

(16)

𝐷 ′

𝑚

𝜆 (∫ 𝐻(Φ)𝑑Ω) = 𝜆 ∑ ∫ 𝛿(Φ) 𝐷

𝑖=1 𝐷

𝜕Φ 𝑖 𝑉 |∇Φ𝑖 |𝑑Ω 𝜕Φ𝑖 𝑛

(17)

It is noted that, only design-independent boundary traction forces are considered and the loading area will be non-designable. Substituting Eq. (14-17) into Eq. (13). Collect all the terms including 𝐰 ′ and the sum is shown in Eq. (18) which naturally equals to zero. ∫ 𝐄B 𝐞(𝐮)𝐞(𝐰 ′ )𝐻(Φ)𝑑Ω − ∫ 𝛕𝐰′𝛿(Φ)|∇Φ|𝑑Ω = 0 𝐷

(18)

𝐷

Then, collect the terms containing 𝐮′ and make the sum equal to zero, that is: 2 ∫ 𝐄B 𝐞(𝐮)𝐞(𝐮′)𝐻(Φ)𝑑Ω + ∫ 𝐄B 𝐞(𝐮′)𝐞(𝐰)𝐻(Φ)𝑑Ω = 0 𝐷

(19)

𝐷

Through solving Eq. (19), the adjoint variable 𝐰 = −𝟐𝐮 can be derived. Finally, by collecting the remaining terms, the sensitivity analysis result is obtained as, 𝑚

𝑚

𝑖=1

𝑖=1

𝜕Φ 𝑖 𝜕𝐄B 𝐿 = ∑ ∫ 𝑅𝛿(Φ) 𝑉𝑛 |∇Φ𝑖 |𝑑Ω − ∑ ∫ 𝐞(𝐮)𝐞(𝐮)𝐻(Φ)𝑉𝑛 𝑖 |∇Φ𝑖 |𝑑Ω 𝜕Φ 𝜕Φ 𝑖 𝑖 𝐷 𝐷 ′

(20)

𝑅 = −𝐄B 𝐞(𝐮)𝐞(𝐮) + 𝜆 Here, it is worth noticing that, the second term in Eq. (20) has the velocities 𝑉𝑛 𝑖 defined inside the whole material domain. However, the non-zero level set contours cannot be freely evolved because of the signed distance feature. For this reason, it is non-trivial to determine the boundary velocities to ensure the Lagrangian to change in the steepest descent direction. To fix this issue, we neglect the second term in 10

Eq. (20) inspired by [63,64], and with this simplification, the optimization problem still converges. Beyond that, a more rigorous treatment of Eq. (20) will be explored in Section 6. Then, by following Eq. (21), 𝑉𝑛 𝑖 = −𝑅

(21)

𝐿 could be guaranteed to change in the descent direction, as shown in Eq. (22): 𝑚 ′

𝐿 = − ∑ ∫ 𝑅 2 𝛿(Φ) 𝑖=1 𝐷

𝜕Φ |∇Φ𝑖 |𝑑Ω 𝜕Φ𝑖

(22)

With the sensitivity result, design update is performed by solving the Hamilton-Jacobi equation; see Eq. (8), through the standard up-winding scheme [53]. (23)

Φ𝑡 = 𝑉𝑛 |∇Φ(𝐗)|

4. Numerical examples In this section, a few numerical examples will be studied to demonstrate the effectiveness of the proposed HO optimization. 4.1 Michell structure First, the Michell structure problem is studied to minimize the structural compliance under the maximum material volume ratio of 0.4. Given the boundary condition, a point force of magnitude 1 is loaded at the center point of the bottom edge, and the two bottom corners are fixed with zero displacement. The design domain has the size of 100-by-50. As demonstrated in Fig. 7, two schemes will be explored: (i) the interior source profile employs the Young’s modulus of 1.3 and the exterior source profile has the Young’s modulus of 0.65; and (ii) the interior source profile employs the Young’s modulus of 0.65 and the exterior source profile has the Young’s modulus of 1.3. The Poisson ratio is assumed to be constantly 0.4. Young’s modulus = 0.65

Young’s modulus = 1.3

(a) Problem setup of scheme 1

Young’s modulus = 1.3

Young’s modulus = 0.65

(b) Problem setup of scheme 2

Fig. 7 Initial problem setups of the two different schemes

11

(a) Scheme 1 with 𝑘 = 1 (compliance = 10.70)

(b) Scheme 1 with 𝑘 = 2 (compliance = 11.26)

(c) Scheme 1 with 𝑘 = 3 (compliance = 11.47)

(d) Scheme 2 with 𝑘 = 1 (compliance = 7.84)

(e) Scheme 2 with 𝑘 = 2 (compliance = 7.59)

(f) Scheme 2 with 𝑘 = 3 (compliance = 7.53)

Fig. 8 Optimization results of the Michell structure problem (the color bar represents the Young’s modulus) Optimization results are demonstrated in Fig. 8, about which the following conclusions can be drawn. (1) If the interior source profile is stronger than the exterior source profile; see Fig. 8(a-c), the linear material blending leads to the best structural performance, because the interior source profile would affect a larger area, which means more strong materials distributed inside the design domain, especially the two side supporting legs. (2) If the interior source profile is weaker than the exterior source profile; see Fig. 8(d-f), the linear material blending leads to the worst structural performance, still because the interior source profile would affect a larger area compared to the higher-order material blendings. We can see from Fig. 8(e-f) that, the two side supporting legs are almost fully occupied by the strong material when 𝑘 = 2 𝑜𝑟 3. (3) By comparing the two schemes, the two side supporting legs of scheme 1 are thicker than that of scheme 2, because the former is mainly composed of the weak material so that more materials accumulate there to bear the load. Hence, the design result of scheme 1 is more compact. 4.2 Cantilever structure 12

This subsection will investigate the cantilever problem. Again, the design problem is to minimize the structural compliance under the maximum material volume ratio of 0.4. Given the boundary condition, a point force of magnitude 1 is loaded at the center point of the right-side edge, and the left-side edge is fixed with zero displacement. The design domain has the size of 100-by-50. As demonstrated in Fig. 9, the two schemes as introduced in the last example will be re-studied with the same material properties. Young’s modulus = 0.65

Young’s modulus = 1.3

Young’s modulus = 1.3

Young’s modulus = 0.65

(a) Problem setup of scheme 1

(b) Problem setup of scheme 2

Fig. 9 Initial problem setups of the two different schemes

(a) Scheme 1 with 𝑘 = 1 (compliance = 99.15)

(b) Scheme 1 with 𝑘 = 2 (compliance = 99.93)

(c) Scheme 1 with 𝑘 = 3 (compliance = 101.59)

(d) Scheme 2 with 𝑘 = 1 (compliance = 81.23)

(e) Scheme 2 with 𝑘 = 2 (compliance = 77.99)

(f) Scheme 2 with 𝑘 = 3 (compliance = 77.71)

Fig. 10 Optimization results of the cantilever problem (the color bar represents the Young’s modulus) 13

Optimization results are demonstrated in Fig. 10 and the following discussions are presented through a detailed analysis. (1) If the interior source profile is stronger than the exterior source profile; see Fig. 10(a-c), the linear material blending (𝑘 = 1) derives the best structural performance, because the interior source profile would affect a larger area, which means more strong materials distributed inside the design domain. This matches the observations in the Michell structure example. (2) If the interior source profile is weaker than the exterior source profile; see Fig. 10(d-f), the linear material blending (𝑘 = 1) leads to the worst structural performance, still because the interior source profile would affect a larger area compared to the higher-order material blendings. (3) By comparing the two schemes, scheme 2 has realized much better structural performance because of more strong materials used and correctly distributed. In addition, the optimized shape of scheme 2 is much closer to conventional topology optimization results. 4.3 Cantilever structure with topological changes In this subsection, the cantilever problem is studied again by placing more interior voids inside the design domain to produce topological changes. Figure 11 presents the initial problem setup. Young’s modulus = 1.3

Young’s modulus = 0.65

Fig. 11 Initial problem setup of the cantilever structure problem

(a) Result with 𝑘 = 1 (compliance = 81.78)

(b) Result with 𝑘 = 2 (compliance = 85.36)

14

(c) Result with 𝑘 = 3 (compliance = 85.61)

Fig. 12 Topology optimization results of the cantilever problem (the color bar represents the Young’s modulus) From Fig. 12, we can see that the structural topologies effectively evolve to the optima. However, a different phenomenon has been observed as compared to the last cantilever example. The linear material blending this time derives the best structural performance. For the reason, the involved structural topology is more complex than the last cantilever example, where more interior ribs are involved which occupy a large area. Therefore, a lower-order material blending makes the interior ribs stronger; refer to Fig. 12a. Besides, the right half of Fig. 12a shrinks a bit in the vertical direction, which further enhances the interior supporting ribs. In addition, the results in Fig. 10(d-f) employ better structural performance than the ones in Fig. 12, even though their topology structure is less optimal. For the reason, the average material Young’s modulus employed in Fig. 10(d-f) is bigger than that of Fig.12. It is interesting to summarize that, a better structural topology does not produce a better structural performance for the FGM structures. 4.4 Michell structure with heterogeneous source profile This subsection studies the Michell structure problem with heterogeneous source profiles. The initial problem setups are demonstrated in Fig. 13a and 13c. The heterogeneous source profile in Fig. 13a employs the Young’s modulus combination of 1.3 and 0.33, while the heterogeneous source profile in Fig. 13c employs the Young’s modulus combination of 1.3 and 0.98. Only linear material blending is considered in these examples.

(a) Young’s modulus combination of 1.3 and 0.33

(b) Optimization result of (a)

15

(c) Young’s modulus combination of 1.3 and 0.98

(d) Optimization result of (c)

Fig. 13 Optimization results of the Michell structure problem with heterogeneous source profiles The optimization results are demonstrated in Fig. 13b and 13d, respectively. We can see in both results, the heterogeneous source profiles have evolved into being homogeneous with only the strong source profile remaining. This is reasonable given the design objective of structural compliance minimization. 4.5 Varying order of smoothness As summarized from the numerical results in the last section, the order of smoothness has a strong influence in the structural performance. Hence, this section intends to explore the level set-based HO optimization with a varying order of smoothness. If the 𝑘 employed in Eq. (9) is also applied as a variable, the related sensitivity result would be: 𝜕𝐿 𝜕𝑘 𝜕𝐄B 𝜕𝑘 = −∫ 𝐞(𝐮)𝐞(𝐮) 𝐻(Φ)𝑑Ω 𝜕𝑘 𝜕𝑡 𝜕𝑘 𝜕𝑡 𝐷 𝑤ℎ𝑒𝑟𝑒 𝑚

𝑘

𝑘

[1/|Φ𝑖 (𝐗)|] ∗ 𝑙𝑛[1/|Φ𝑖 (𝐗)|] ∗ ∑𝑚 𝜕𝐄B 𝑗=1[1/|Φ𝑗 (𝐗)|] =∑ 𝐄𝑖 𝑘 𝑚 2 𝜕𝑘 (𝐗) {∑ [1/| Φ |] } 𝑗 𝑗=1 𝑖 𝑚

−∑ 𝑖

𝑘

(24)

𝑘

[1/|Φ𝑖 (𝐗)|] ∗ ∑𝑚 𝑗=1{[1/|Φ𝑗 (𝐗)|] ∗ 𝑙𝑛[1/|Φ𝑗 (𝐗)|]} 𝑘

2 {∑𝑚 𝑗=1[1/|Φ𝑗 (𝐗)|] }

𝐄𝑖

Therefore, 𝜕𝑘 𝜕𝐄B (25) = 𝐞(𝐮)𝐞(𝐮) 𝜕𝑡 𝜕𝑘 In addition, the smoothness 𝑘 should be constrained by both upper and lower bounds given certain engineering requirements. Therefore, the constraints as presented in Eq. (26) are employed. 𝑘 < 𝑘̅ 𝑘>𝑘 ̅ where 𝑘 is the upper bound and 𝑘 is the lower bound.

(26)

To prove the effect of the varying 𝑘, the examples presented in Fig. 7a and 9a are re-studied. 𝑘 is bounded within the range [1,3]. The optimization results are demonstrated in Fig. 14 and the structural compliance has been reduced to 10.65 and 97.59, respectively. We have observed the phenomenon that 𝑘 tends to be larger at areas filled with stronger materials, and smaller at areas with weaker materials. And the overall effect is to have more strong materials distributed inside the optimal design. 16

(a) Young’s modulus distribution

(b) 𝑘 value distribution

(c) Young’s modulus distribution

(d) 𝑘 value distribution

Fig. 14 Optimization results with varying 𝑘 4.6 3D L-bracket problem At the end, a 3D L-bracket problem is demonstrated for compliance minimization under the maximum material volume ratio of 0.4. Referring to Fig. 15, a point force of magnitude 1 is loaded at the center of the right-side top edge, and the top surface is fixed with zero displacement. The design domain has the size of 80*80*20. Two schemes will be explored: (i) the side frame employs the Young’s modulus of 1.0 and the interior holes have the Young’s modulus of 0.5; and (ii) the side frame employs the Young’s modulus of 0.5 and the interior holes have the Young’s modulus of 1.0. The Poisson ratio is assumed to be constantly 0.4. Note that, the front, back, and top surfaces are not assigned of any material source, and the linear material blending is considered.

17

Fig. 15 Initial problem setups of the two different schemes

(a) Optimization result of scheme 1 (compliance = 19.73)

(b) Optimization result of scheme 2 (compliance = 25.10) Fig. 16 Optimization results of the 3D L-bracket problem The optimization results are demonstrated in Fig. 16. Similar to the former cases, the scheme 2 that has the weaker exterior material source contracts more to from more internal structures that are made of the strong material; and the scheme 1 derives the geometry that is more approaching the conventional topology design.

5. Heterogeneous meta-material optimization Meta-material modeling is critical for porous structure optimization, which homogenizes the properties of local representative volumes, so that the macro structural analysis can be performed with a coarse mesh. Because of the employed coarse non-conformal mesh, the computational cost can be greatly saved. Even 18

though the representative volume is already in the micro or meso scale, a heterogeneous design is still meaningful for multi-functioning considerations, e.g., surfaces connecting other representative volumes should be stiff, while surfaces forming the voids should employ certain physical or chemical properties, because of the contact to media. Definitely, the level set-based HO model proposed in this work can be trivially extended to model metal-material, and the details about heterogeneous meta-material optimization will be presented in this section. In the past, both meta-material optimization [65–68] and two-scale topology optimization [69] have been conducted under the level set framework, but heterogeneous meta-material optimization has never been addressed. Therefore, a related method will be developed in this section. Under the level set framework, the homogenized elasticity tensor of the representative volume can be calculated by Eq. (27). 1 (27) ∫ 𝐄 (𝐞0 − 𝐞(𝐮∗ ))(𝐞0 − 𝐞(𝐮∗ )) 𝐻(Φ)𝑑Ω |𝑌| 𝐷 B where, |𝑌| is the representative volume area, 𝐞0 is the applied unit strain fields, e.g. (1,0,0), (0,1,0), and (0,0,1). 𝐮∗ is the perturbed displacement field which is obtained by solving Eq. (28), which is Y-period. 𝐄B H =

∫ 𝐄B (𝐞0 − 𝐞(𝐮∗ ))𝐞(𝐯) 𝐻(Φ)𝑑Ω = 0,

∀𝐯∈𝑈

(28)

𝐷

Therefore, to design material with specified properties, the optimization problem is formulated as follows: 1 min. 𝐽 = (𝐸𝐵𝑖𝑗 𝐻 − 𝐸𝑖𝑗 )2 2 𝑠. 𝑡. 𝑎(𝑬𝐵 , 𝒖∗ , 𝒗, 𝛷) = 𝑙(𝒗, 𝛷), ∀𝒗 ∈ 𝑈 𝑎(𝑬𝐵 , 𝒖∗ , 𝒗, 𝛷) = ∫ 𝑬𝐵 𝒆(𝒖∗ )𝒆(𝒗)𝐻(𝛷)𝑑𝛺

(29)

𝐷

𝑙(𝒗, 𝛷) = ∫ 𝑬𝐵 𝒆0 𝒆(𝒗)𝐻(𝛷)𝑑𝛺 𝐷

where, 𝐸𝑖𝑗 is the targeted value of the homogenized 𝐸𝐵𝑖𝑗 𝐻 . There is no inequality constraints in Eq. (29) so that the design result would not be unique. By performing a similar process as conducted in subsection 3.2, the sensitivity result can be derived as: 𝑉𝑛 = −(𝐸𝐵𝑖𝑗 𝐻 − 𝐸𝑖𝑗 )𝐄B (𝐞0 − 𝐞(𝐮∗ ))(𝐞0 − 𝐞(𝐮∗ ))

(30)

Then, a few numerical examples are studied to prove the effectiveness of the proposed heterogeneous meta-material optimization method. The input material distribution is demonstrated in Fig. 17. Two material source profiles are employed where the strong material has the Young’s modulus of 1.3 and the weak material employs the Young’s modulus of 0.65. The Poisson ratio is constantly 0.4. Correspondingly, a few targeted values of the homogenized elasticity tensor elements are utilized as the objective, and the related optimization results are demonstrated in Fig. 18.

19

Young’s modulus = 1.3

Young’s modulus = 0.65

Fig. 17 Input material distribution of the meta-material

(a) Meta-material unit with 𝐸11 = 0.8, 𝐸22 = 0.4, and the related repetitive structure

(b) Meta-material unit with 𝐸11 = 0.7, 𝐸22 = 0.5, and the related repetitive structure

20

(c) Meta-material unit with 𝐸11 = 0.6, 𝐸22 = 0.6, and the related repetitive structure

(d) Meta-material unit with 𝐸12 = 0.2, and the related repetitive structure

(e) Meta-material unit with 𝐸12 = 0.15, and the related repetitive structure

21

(f) Meta-material unit with 𝐸12 = 0.1, and the related repetitive structure

(g) Meta-material unit with 𝐸12 = 0.05, and the related repetitive structure

Fig. 18 Heterogeneous meta-material optimization results 6. Sensitivity result transformation As discussed earlier, the sensitivity result of Eq. (20) has a domain integral term, which is nonimplementable, and this term is ignored in the numerical examples of Section 4. Even though good design effect has been achieved, we would like to explore a more rigorous approach to dealing with this nonimplementable term, i.e. performing sensitivity result transformation, to see how much improvement can be made with the extra numerical efforts. Specifically, the sensitivity result transformation presented in the authors’ earlier work [41] will be used, and the details will be briefly presented below. As the basis of this transformation, a corollary is cited from [70]: Corollary 1. For a 2D integrable function Φ(𝐗) as demonstrated in Fig. 19, ∫ Φ(𝐗)𝑑Ω = ∫ (∫ Ω

𝜕Ω

Φ(𝐙)(1 − 𝑑Ω (𝐙)𝜅(𝐘))𝑑𝐙)𝑑Γ

ray𝜕Ω (𝐘)∩Ω

(31)

Two definitions are presented below to interpret this corollary.

22

Definition 1. For any 𝐗 ∈ 𝑅 𝑛 , Π𝜕𝛺 (𝐗) ≔ {𝐘0 ∈ 𝜕Ω, |𝐗 − 𝐘0 | = 𝑖𝑛𝑓𝐘∈𝜕Ω |𝐗 − 𝐘|} is the set of projections of 𝐗 on 𝜕Ω. When Π𝜕Ω (𝐗) reduces to a single point, it is called the projection 𝑃𝜕Ω (𝐗) of 𝐗 onto 𝜕Ω. Definition 2. For any 𝐘 ∈ 𝜕Ω, 𝑟𝑎𝑦𝜕Ω (𝐘) ≔ {𝐙 ∈ 𝑅 𝑛 , 𝑑Ω 𝑖𝑠 𝑑𝑖𝑓𝑓𝑒𝑟𝑒𝑛𝑡𝑖𝑎𝑏𝑙𝑒 𝑎𝑡 𝐙 𝑎𝑛𝑑 𝑃𝜕Ω (𝐙) = 𝐘} is the ray emerging from 𝒀.

Fig. 19 Schematic plot of corollary 1 [41] In Corollary 1, 𝑑Ω (𝐙) represents the distance which equals to Φ(𝐙) in case that the signed distance regulation is satisfied. By applying Corollary 1, the sensitivity result of Eq. (20) can be changed into： 𝑚 ′

𝐿 = ∑ ∫ 𝑅𝛿(Φ) 𝑖=1 𝐷 𝑚

− ∑ ∫ (∫ 𝑖=1 𝐷

ray𝜕Ω

𝜕Φ 𝑖 𝑉 |∇Φ𝑖 |𝑑Ω 𝜕Φ𝑖 𝑛

𝜕𝐄B 𝐞(𝐮)𝐞(𝐮)(1 − 𝑑Ω (𝐙)𝜅(𝐘))𝑑𝐙) 𝛿(Φ𝑖 )𝑉𝑛 𝑖 |∇Φ𝑖 |𝑑Ω 𝜕Φ 𝑖 (𝐘)∩Ω

(32)

𝑅 = −𝐄B 𝐞(𝐮)𝐞(𝐮) + 𝜆 So far, sensitivity result transformation has been completed. The Michell structure problem of Subsection 5.1 and the cantilever problem of Sub-section 5.2 will be re-studied. Note that only the situation of 𝑘 = 1 will be re-considered. Correspondingly, the optimization results are shown in Fig. 18 and 19, respectively.

(a) Optimization result 1 of the Michell structure problem (compliance =10.78), which corresponds to Fig. 7a

23

(b) Optimization result 2 of the Michell structure problem (compliance =7.39), which corresponds to Fig. 7d

Fig. 18 Optimization results of the Michell structure problem

(a) Optimization result 1 of the cantilever problem (compliance =87.79), which corresponds to Fig. 9a

(b) Optimization result 2 of the cantilever problem (compliance =76.54), which corresponds to Fig. 9d

Fig. 19 Optimization results of the cantilever problem The following conclusions can be drawn from the modified optimization results: 24

(1) The modified results from Fig. 18 and 19 have reduced the structural compliance by -0.74%, 5.74%, 11.46%, and 5.76%, respectively, which means the full sensitivity can generally lead to better structural optimal. The only exception is Fig. 18a which very likely has been trapped at a local optimum. (2) It is interesting to find out that the internal material source profile has contracted into a point in Fig. 18(b), so that more strong materials were distributed inside the design domain.

At the end, we would like to explore one more example that one material source profile keeps its shape during optimization while its position will be optimized. Specifically, the MBB structure problem will be studied and the problem setup is shown in Fig. 20. The material volume fraction limit is 50%. We can see that there is a black circle inside the design domain, which corresponds to the fixed-shape material source profile, i.e., its shape will not vary but its position will be optimized. This fixed-shape profile employs Young’s modulus of 0.65, while all other boundaries have the Young’s modulus of 1.3. The Poisson ratio is assumed to be constantly 0.4. In addition, two initial guesses will be used to perform the optimization, since this type of moving feature problem is very sensitive to the feature’s starting position [71–73].

(a) Boundary condition

(b) Initial guess 1

(c) Initial guess 2

Fig. 20 Setup of the MBB problem: only the right half is shown in (b) and (c) because of symmetry Then, the optimization for each initial guess is performed with two different sensitivity results, i.e. the simplified sensitivity of Eq. (20) and the full sensitivity of Eq. (32). And the optimization results are shown in Fig. 21. We can conclude from the results that, optimization with the full sensitivity can lead to better structural optimal, which is at least true for this case.

25

(a) Optimization result with initial guess 1 with the simplified sensitivity of Eq. (20): compliance = 66.42

(b) Optimization result with initial guess 1 with full sensitivity of Eq. (32): compliance = 65.77

(c) Optimization result with initial guess 2 with the simplified sensitivity of Eq. (20): compliance = 67.02

(d) Optimization result with initial guess 2 with full sensitivity of Eq. (32): compliance = 65.80

Fig. 21 Optimization results

7. Conclusion This paper introduces a level set-based heterogeneous object (HO) modeling and optimization method. The individual level set contours are utilized as material source profiles, and the heterogeneous material composition distribution is calculated based on a signed distance-based material blending function. More importantly, a level set-based HO optimization method is developed, which has demonstrated unique advantages in designing functionally graded structures. A large number of numerical examples have been studied to demonstrate the effectiveness of the proposed optimization method. We have found that, the material blending smoothness has demonstrated a strong influence in the optimization result. Another interesting finding is that, placing more interior holes inside the initial design domain may not lead to better structural optima, which is very different from the conventional level set topology optimization. Functionally graded meta-material optimization has also been successfully performed, which to the best of the authors’ knowledge, is a never studied topic for topology optimization. About the numerical optimization algorithm, another main contribution of this paper is that, both simplified sensitivity result and the full sensitivity have been used for implementation. Even though the latter requires extra numerical efforts to implement, it generally would lead to better optimal results. For future work, it is intended to perform the concurrent two-scale HO optimization by bridging the macro- and micro-scale HO optimization techniques. Another important issue to address is that, the shape derivative is used to evolve the structure so that the algorithm lacks the capability of hole nucleation. This leads to the initial guess dependency issue. 26

Therefore, certain hole nucleation mechanism should be introduced into the proposed method; for instance, the material removal scheme of bi-directional evolutionary structural optimization (BESO) was recently introduced into level set framework for stable hole nucleation-type topological changes [74]. On the other hand, since each boundary profile is associated with a material source, the difficulty arises about how to assign the appropriate material source to the newly produced boundary. These issues will be focused in our forth-coming research.

Acknowledgments The authors would like to acknowledge the support from the Qilu Young Scholar award, Shandong University, and the Open Research Fund of Key Laboratory of High Performance Complex Manufacturing, Central South University [grant number Kfkt2016-07].

Reference [1]

[2] [3]

[4] [5] [6]

[7] [8] [9] [10] [11] [12] [13]

Liu J, Gaynor AT, Chen S, Kang Z, Suresh K, Takezawa A, et al. Current and future trends in topology optimization for additive manufacturing. Struct Multidiscip Optim 2018;57:2457–2483. doi:10.1007/s00158-018-1994-3. Kou XY, Tan ST. Heterogeneous object modeling: A review. Comput-Aided Des 2007;39:284–301. doi:10.1016/j.cad.2006.12.007. Jackson TR, Liu H, Patrikalakis NM, Sachs EM, Cima MJ. Modeling and designing functionally graded material components for fabrication with local composition control. Mater Des 1999;20:63– 75. doi:10.1016/S0261-3069(99)00011-4. Wu X, Liu W, Wang MY. A CAD modeling system for heterogeneous object. Adv Eng Softw 2008;39:444–53. doi:10.1016/j.advengsoft.2007.03.002. Wang S, Chen N, Chen C-S, Zhu X. Finite element-based approach to modeling heterogeneous objects. Finite Elem Anal Des 2009;45:592–6. doi:10.1016/j.finel.2009.02.005. Liu J, Duke K, Ma Y. Computer-aided design–computer-aided engineering associative featurebased heterogeneous object modeling. Adv Mech Eng 2015;7:1687814015619767. doi:10.1177/1687814015619767. Kumar V, Dutta D. An Approach to Modeling & Representation of Heterogeneous Objects. J Mech Des 1998;120:659–67. doi:10.1115/1.2829329. Kumar V, Burns D, Dutta D, Hoffmann C. A framework for object modeling. Comput-Aided Des 1999;31:541–56. doi:10.1016/S0010-4485(99)00051-2. Srinivas Bhashyam, Ki Hoon Shin, Debashish Dutta. An integrated CAD system for design of heterogeneous objects. Rapid Prototyp J 2000;6:119–35. doi:10.1108/13552540010323547. Shin K-H, Dutta D. Constructive Representation of Heterogeneous Objects. J Comput Inf Sci Eng 2001;1:205–17. doi:10.1115/1.1403448. Shin K-H, Natu H, Dutta D, Mazumder J. A method for the design and fabrication of heterogeneous objects. Mater Des 2003;24:339–53. doi:10.1016/S0261-3069(03)00060-8. Siu YK, Tan ST. ‘Source-based’ heterogeneous solid modeling. Comput-Aided Des 2002;34:41–55. doi:10.1016/S0010-4485(01)00046-X. Siu YK, Tan ST. Modeling the material grading and structures of heterogeneous objects for layered manufacturing. Comput-Aided Des 2002;34:705–16. doi:10.1016/S0010-4485(01)00200-7.

27

[14] Kou XY, Tan ST. A hierarchical representation for heterogeneous object modeling. Comput-Aided Des 2005;37:307–19. doi:10.1016/j.cad.2004.03.006. [15] Kou XY, Tan ST, Sze WS. Modeling complex heterogeneous objects with non-manifold heterogeneous cells. Comput-Aided Des 2006;38:457–74. doi:10.1016/j.cad.2005.11.009. [16] Kou XY, Tan ST. A systematic approach for Integrated Computer-Aided Design and Finite Element Analysis of Functionally-Graded-Material objects. Mater Des 2007;28:2549–65. doi:10.1016/j.matdes.2006.10.024. [17] Liu H, Maekawa T, Patrikalakis NM, Sachs EM, Cho W. Methods for feature-based design of heterogeneous solids. Comput-Aided Des 2004;36:1141–59. doi:10.1016/j.cad.2003.11.001. [18] Biswas A, Shapiro V, Tsukanov I. Heterogeneous material modeling with distance fields. Comput Aided Geom Des 2004;21:215–42. doi:10.1016/j.cagd.2003.08.002. [19] Xu A, Shaw LL. Equal distance offset approach to representing and process planning for solid freeform fabrication of functionally graded materials. Comput-Aided Des 2005;37:1308–18. doi:10.1016/j.cad.2005.01.005. [20] Zhou H, Liu Z, Lu B. Heterogeneous object modeling based on multi-color distance field. Mater Des 2009;30:939–46. doi:10.1016/j.matdes.2008.07.002. [21] Gupta V, Tandon P. Heterogeneous object modeling with material convolution surfaces. ComputAided Des 2015;62:236–47. doi:10.1016/j.cad.2014.12.005. [22] Qian X, Dutta D. Physics-Based Modeling for Heterogeneous Objects. J Mech Des 2003;125:416– 27. doi:10.1115/1.1582877. [23] Qian X, Dutta D. Feature-based design for heterogeneous objects. Comput-Aided Des 2004;36:1263–78. doi:10.1016/j.cad.2004.01.012. [24] Qian X, Dutta D. Design of heterogeneous turbine blade. Comput-Aided Des 2003;35:319–29. doi:10.1016/S0010-4485(01)00219-6. [25] Yang P, Qian X. A B-spline-based approach to heterogeneous objects design and analysis. ComputAided Des 2007;39:95–111. doi:10.1016/j.cad.2006.10.005. [26] Samanta K, Koc B. Feature-based design and material blending for free-form heterogeneous object modeling. Comput-Aided Des 2005;37:287–305. doi:10.1016/j.cad.2004.03.005. [27] Ozbolat IT, Koc B. Multi-directional blending for heterogeneous objects. Comput-Aided Des 2011;43:863–75. doi:10.1016/j.cad.2011.04.002. [28] Samanta K, Ozbolat IT, Koc B. Optimized normal and distance matching for heterogeneous object modeling. Comput Ind Eng 2014;69:1–11. doi:10.1016/j.cie.2013.12.010. [29] Yoo D-J. Heterogeneous object modeling using the radial basis functions. Int J Precis Eng Manuf 2013;14:1133–40. doi:10.1007/s12541-013-0154-3. [30] Kou XY, Tan ST. A simple and effective geometric representation for irregular porous structure modeling. Comput-Aided Des 2010;42:930–41. doi:10.1016/j.cad.2010.06.006. [31] Kou XY, Tan ST. Microstructural modelling of functionally graded materials using stochastic Voronoi diagram and B-Spline representations. Int J Comput Integr Manuf 2012;25:177–88. doi:10.1080/0951192X.2011.627948. [32] Chen K-Z, Feng X-A. Computer-aided design method for the components made of heterogeneous materials. Comput-Aided Des 2003;35:453–66. doi:10.1016/S0010-4485(02)00069-6. [33] Chen K-Z, Feng X-A. CAD modeling for the components made of multi heterogeneous materials and smart materials. Comput-Aided Des 2004;36:51–63. doi:10.1016/S0010-4485(03)00077-0. [34] Kou XY, Parks GT, Tan ST. Optimal design of functionally graded materials using a procedural model and particle swarm optimization. Comput-Aided Des 2012;44:300–10. doi:10.1016/j.cad.2011.10.007. [35] Wang MY, Wang X. “Color” level sets: a multi-phase method for structural topology optimization with multiple materials. Comput Methods Appl Mech Eng 2004;193:469–96. doi:10.1016/j.cma.2003.10.008. [36] Wang MY, Wang X. A level-set based variational method for design and optimization of heterogeneous objects. Comput-Aided Des 2005;37:321–37. doi:10.1016/j.cad.2004.03.007. 28

[37] 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. doi:10.1016/j.cma.2013.10.003. [38] Wang Y, Luo Z, Kang Z, Zhang N. A multi-material level set-based topology and shape optimization method. Comput Methods Appl Mech Eng 2015;283:1570–86. doi:10.1016/j.cma.2014.11.002. [39] Vermaak N, Michailidis G, Parry G, Estevez R, Allaire G, Bréchet Y. Material interface effects on the topology optimizationof multi-phase structures using a level set method. Struct Multidiscip Optim 2014;50:623–44. doi:10.1007/s00158-014-1074-2. [40] Cui M, Chen H, Zhou J. A level-set based multi-material topology optimization method using a reaction diffusion equation. Comput-Aided Des 2016;73:41–52. doi:10.1016/j.cad.2015.12.002. [41] Liu J, Ma Y. A new multi-material level set topology optimization method with the length scale control capability. Comput Methods Appl Mech Eng 2018;329:444–63. doi:10.1016/j.cma.2017.10.011. [42] Zhang W, Song J, Zhou J, Du Z, Zhu Y, Sun Z, et al. Topology optimization with multiple materials via moving morphable component (MMC) method. Int J Numer Methods Eng 2018;113:1653–75. doi:10.1002/nme.5714. [43] 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. doi:10.1016/j.cma.2017.09.017. [44] Fu J, Li H, Gao L, Xiao M. Design of shell-infill structures by a multiscale level set topology optimization method. Comput Struct 2019;212:162–72. doi:10.1016/j.compstruc.2018.10.006. [45] Xia Q, Wang MY. Simultaneous optimization of the material properties and the topology of functionally graded structures. Comput-Aided Des 2008;40:660–75. doi:10.1016/j.cad.2008.01.014. [46] Dunning PD, Brampton CJ, Kim HA. Simultaneous optimisation of structural topology and material grading using level set method. Mater Sci Technol 2015;31:884–94. doi:10.1179/1743284715Y.0000000022. [47] Almeida SRM, Paulino GH, Silva ECN. Layout and material gradation in topology optimization of functionally graded structures: a global–local approach. Struct Multidiscip Optim 2010;42:855–68. doi:10.1007/s00158-010-0514-x. [48] Taheri AH, Hassani B. Simultaneous isogeometrical shape and material design of functionally graded structures for optimal eigenfrequencies. Comput Methods Appl Mech Eng 2014;277:46–80. doi:10.1016/j.cma.2014.04.014. [49] Chen J, Shapiro V, Suresh K, Tsukanov I. Shape optimization with topological changes and parametric control. Int J Numer Methods Eng 2007;71:313–46. doi:10.1002/nme.1943. [50] Nogata F, Takahashi H. Intelligent functionally graded material: Bamboo. Compos Eng 1995;5:743–51. doi:10.1016/0961-9526(95)00037-N. [51] Tan T, Rahbar N, Allameh SM, Kwofie S, Dissmore D, Ghavami K, et al. Mechanical properties of functionally graded hierarchical bamboo structures. Acta Biomater 2011;7:3796–803. doi:10.1016/j.actbio.2011.06.008. [52] Gupta A, Talha M. Recent development in modeling and analysis of functionally graded materials and structures. Prog Aerosp Sci 2015;79:1–14. doi:10.1016/j.paerosci.2015.07.001. [53] Osher S, Fedkiw R. Level Set Methods and Dynamic Implicit Surfaces. vol. 153. New York, NY: Springer New York; 2003. [54] Wang S, Wang MY. Radial basis functions and level set method for structural topology optimization. Int J Numer Methods Eng 2006;65:2060–90. doi:10.1002/nme.1536. [55] Wang SY, Lim KM, Khoo BC, Wang MY. An extended level set method for shape and topology optimization. J Comput Phys 2007;221:395–421. doi:10.1016/j.jcp.2006.06.029. [56] Liu J, Ma Y-S. 3D level-set topology optimization: a machining feature-based approach. Struct Multidiscip Optim 2015;52:563–82. doi:10.1007/s00158-015-1263-7. [57] Kang Z, Wang Y. Structural topology optimization based on non-local Shepard interpolation of density field. Comput Methods Appl Mech Eng 2011;200:3515–25. doi:10.1016/j.cma.2011.09.001. 29

[58] Kang Z, Wang Y. A nodal variable method of structural topology optimization based on Shepard interpolant. Int J Numer Methods Eng 2012;90:329–42. doi:10.1002/nme.3321. [59] Xia Q, Shi T. Optimization of composite structures with continuous spatial variation of fiber angle through Shepard interpolation. Compos Struct 2017;182:273–82. doi:10.1016/j.compstruct.2017.09.052. [60] Xia Q, Shi T. A cascadic multilevel optimization algorithm for the design of composite structures with curvilinear fiber based on Shepard interpolation. Compos Struct 2018;188:209–19. doi:10.1016/j.compstruct.2018.01.013. [61] Wang MY, Wang X, Guo D. A level set method for structural topology optimization. Comput Methods Appl Mech Eng 2003;192:227–46. doi:10.1016/S0045-7825(02)00559-5. [62] Sigmund O. Design of multiphysics actuators using topology optimization – Part II: Two-material structures. Comput Methods Appl Mech Eng 2001;190:6605–27. doi:10.1016/S00457825(01)00252-3. [63] Guo X, Zhang W, Zhong W. Explicit feature control in structural topology optimization via level set method. Comput Methods Appl Mech Eng 2014;272:354–78. doi:10.1016/j.cma.2014.01.010. [64] Brampton CJ, Wu KC, Kim HA. New optimization method for steered fiber composites using the level set method. Struct Multidiscip Optim 2015;52:493–505. doi:10.1007/s00158-015-1256-6. [65] Wang Y, Luo Z, Zhang N, Kang Z. Topological shape optimization of microstructural metamaterials using a level set method. Comput Mater Sci 2014;87:178–86. doi:10.1016/j.commatsci.2014.02.006. [66] Challis VJ, Xu X, Zhang LC, Roberts AP, Grotowski JF, Sercombe TB. High specific strength and stiffness structures produced using selective laser melting. Mater Des 2014;63:783–8. doi:10.1016/j.matdes.2014.05.064. [67] Wang Y, Luo Z, Zhang N, Wu T. Topological design for mechanical metamaterials using a multiphase level set method. Struct Multidiscip Optim 2016:1–16. doi:10.1007/s00158-016-1458-6. [68] Wang X, Mei Y, Wang MY. Level-set method for design of multi-phase elastic and thermoelastic materials. Int J Mech Mater Des 2004;1:213–39. doi:10.1007/s10999-005-0221-8. [69] Wang Y, Wang MY, Chen F. Structure-material integrated design by level sets. Struct Multidiscip Optim 2016;54:1145–56. doi:10.1007/s00158-016-1430-5. [70] Allaire G, Dapogny C, Delgado G, Michailidis G. Multi-phase structural optimization via a level set method. ESAIM Control Optim Calc Var 2014;20:576–611. [71] Zhang W, Xia L, Zhu J, Zhang Q. Some Recent Advances in the Integrated Layout Design of Multicomponent Systems. J Mech Des 2011;133:104503–104503. doi:10.1115/1.4005083. [72] Liu J, Cheng L, To AC. Arbitrary void feature control in level set topology optimization. Comput Methods Appl Mech Eng 2017;324:595–618. doi:10.1016/j.cma.2017.06.021. [73] Liu J, Yu H, Ma Y. Minimum void length scale control in level set topology optimization subject to machining radii. Comput-Aided Des 2016;81:70–80. doi:10.1016/j.cad.2016.09.007. [74] Xia Q, Shi T, Xia L. Stable hole nucleation in level set based topology optimization by using the material removal scheme of BESO. Comput Methods Appl Mech Eng 2019;343:438–52. doi:10.1016/j.cma.2018.09.002.

30

Highlights    

Presents a level set-based heterogeneous object (HO) modeling and optimization method Features in topology optimization of functionally graded structures Explores the optimal designs with different material blending coefficients Develops the accurate sensitivity result rather than ignoring the non-implementable part of the



original sensitivity Investigates functionally graded meta-material optimization