A new approach to eliminating enclosed voids in topology optimization for additive manufacturing

Journal Pre-proof A new approach to eliminating enclosed voids in topology optimization for additive manufacturing Yulin Xiong (Investigation) (Writing - original draft), Song Yao (Writing - review and editing) (Validation), Zi-Long ZhaoWriting - review and editing, Validation), Yi Min Xie (Conceptualization) (Supervision)

PII:

S2214-8604(19)31551-9

DOI:

https://doi.org/10.1016/j.addma.2019.101006

Reference:

ADDMA 101006

To appear in:

Additive Manufacturing

Received Date:

9 September 2019

Revised Date:

22 November 2019

Accepted Date:

13 December 2019

Please cite this article as: Xiong Y, Yao S, Zhao Z-Long, Xie YM, A new approach to eliminating enclosed voids in topology optimization for additive manufacturing, Additive Manufacturing (2019), doi: https://doi.org/10.1016/j.addma.2019.101006

This is a PDF file of an article that has undergone enhancements after acceptance, such as the addition of a cover page and metadata, and formatting for readability, but it is not yet the definitive version of record. This version will undergo additional copyediting, typesetting and review before it is published in its final form, but we are providing this version to give early visibility of the article. 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. © 2019 Published by Elsevier.

A new approach to eliminating enclosed voids in topology optimization for additive manufacturing

Yulin Xionga,b, Song Yaoa, Zi-Long Zhaob, Yi Min Xieb,* a

Centre for Innovative Structures and Materials, School of Engineering, RMIT University, Melbourne 3001, Australia *

Corresponding author. E-mail address: [email protected] (Y.M. Xie).

ro of

b

Key Laboratory of Traffic Safety on Track, Ministry of Education, School of Traffic & Transportation Engineering, Central South University, Changsha 410075, China

-p

Highlights:

Abstract

na

lP

re

 A new approach to controlling the structural connectivity is proposed.  Graph theory is introduced in combination with the BESO technique to eliminate enclosed voids.  Highly efficient structural designs without enclosed voids can be achieved.  The new approach solves a problem of practical importance in additive manufacturing.

Jo

ur

Topology optimization is increasingly used in lightweight designs for additive manufacturing (AM). However, conventional optimization techniques do not fully consider manufacturing constraints. One important requirement of powder-based AM processes is that enclosed voids in the designs must be avoided in order to remove and reuse the unmelted powder. In this work, we propose a new approach to realizing the structural connectivity control based on the bidirectional evolutionary structural optimization technique. This approach eliminates enclosed voids by selectively generating tunnels that connect the voids with the structural boundary during the optimization process. The developed methodology is capable of producing highly efficient structural designs which have no enclosed voids. Furthermore, by changing the radius and the number of tunnels, competitive and diverse designs can be achieved. The effectiveness of the approach is demonstrated by two examples of three-dimensional structures. Prototypes of the obtained designs without enclosed voids have been fabricated using AM. 1

Jo

ur

na

lP

re

-p

ro of

Keywords: Topology optimization; BESO; Powder bed fusion; Manufacturing constraint; Structural connectivity control

2

1. Introduction Additive manufacturing (AM) has become a powerful tool in a wide range of industries in recent years. Different from traditional manufacturing, AM builds 3D components on a layerby-layer basis [1], and thus it can be used to fabricate highly complex structures [2]. Various AM techniques have been developed, including selective laser sintering (SLS) [3] and electron beam melting (EBM) [4]. In order to improve the performance of products, some studies focus on the enhancement of the properties of printed materials [5]. Peng et al. reinforced the composites significantly by using short and continuous carbon fibers [6]. Other researchers are devoted to developing design methods for AM to maximize the performance of products [7,8].

lP

re

-p

ro of

Topology optimization has received much attention as an efficient computer-aided design technique [9]. The common goal of topology optimization is to maximize the structural performance while satisfying certain constraints. It is a powerful design method to explore the great potential of AM [10]. Several optimization methods have been developed over the last three decades, e.g., the homogenization method [11,12], the solid isotropic material with penalization (SIMP) method [13,14], the evolutionary structural optimization (ESO) [15,16], the bi-directional evolutionary structural optimization (BESO) [17,18], and the level set method [19–21]. These methods are applied widely in architecture [22], automotive [23], and aerospace industries [24,25]. Some optimized structural designs have been successfully applied in engineering applications such as railway vehicles [26] and aircrafts [27]. Topology optimization is not only used in mechanical design, but also applied in multidisciplinary design, including thermodynamics [28,29], biomechanics [30], acoustics [31,32], micro-structured materials [33–36], and nano-photonic design [37–39].

ur

na

Manufacturability plays a vital role in AM technologies. Some manufacturing limitations of AM need to be considered in topology optimization. A large amount of effort has been devoted to, e.g., the feature size limitation [40,41], manufacturable inclination angle [42,43], and fabrication tolerance [44]. There are two difficulties to impose the manufacturing constraints into topology optimization. Firstly, it is difficult to establish a general method to describe all manufacturing limitations. Different approaches need to be developed to satisfy the manufacturing constraints in topology optimization. Secondly, imposing manufacturing constraints into topology optimization might reduce structural efficiency. It remains a challenge to satisfy manufacturing constraints while guaranteeing structural performance.

Jo

The structural connectivity constraint, one of the major manufacturing limitations in AM, is integrated into the BESO technology in this work. It is required that the final design should have no enclosed voids, which is necessary for powder-based AM technologies. If the optimized structures have enclosed voids, the unmelted powder would be left in the enclosed voids, which will increase the weight of the structures [45]. In order to achieve light-weight designs and reuse unmelted powder, it is of significant importance to remove the unmelted powder from enclosed voids. Several researchers have investigated this manufacturing constraint. Li et al. [46] and Liu et al. [47] proposed a virtual temperature approach to eliminating enclosed voids generated by topology optimization for structural connectivity. In their approach, the void parts are assumed to be filled with a virtual heat-conducting material, 3

so that the structural connectivity is controlled by constraining the maximum temperature. However, an added heat conduction optimization of the structure would incur the additional computational cost and could result in a significant reduction of the structural performance, as we shall demonstrate through an example in this paper. Zhou and Zhang [48] proposed a method for controlling the void parts using closed B-splines and super-ellipses. In their work, the void parts need to be initially distributed and kept outside the structure. In these void parts, the variations outside the structure will be constrained to satisfy the structural connectivity constraint, which might substantially decrease the performance of the optimized structure.

ro of

Most recently, Zhao et al. proposed an effective approach to controlling the structural connectivity in topology optimization [49]. Using this approach, the structural performance and the effect of the structural complexity control can be well balanced. This approach has been extended and successfully applied in the morphological optimization of biological organs [50]. However, it remains a challenging issue to eliminate enclosed voids in an optimized design.

re

-p

In the present work, an efficient structural connectivity control (SCC) approach is proposed based on the BESO method. This approach can evaluate the structural connectivity and find the enclosed voids during the optimization process. In order to eliminate enclosed voids, the approach would find possible paths that connect the voids with the structural boundary. The structural connectivity constraint is satisfied mainly by altering the elements around these paths during the optimization process. Thus, the SCC approach can generate a material distribution which is similar to the conventional optimized result, without significantly sacrificing the structural performance.

na

lP

The paper is organized as follows. In Section 2, the SCC approach is developed based on the BESO method. In Section 3, two examples are used to demonstrate the effectiveness of this approach. In Section 4, the main conclusions of this study are summarized.

2. Structural connectivity control in topology optimization

ur

2.1. BESO method

Jo

In this work, the SCC approach is developed based on the BESO method. The BESO method can be used in multidisciplinary structural optimization, such as stiffness optimization [51,52], natural frequency optimization [53], and the optimal design of functional gradient materials [54], biological materials [30] and concurrent optimization for structures and materials [55]. The BESO method is widely recognized owing to its high-quality topology solutions, simple to understand and implement, and excellent computational efficiency. To demonstrate the effectiveness of the approach, a static, linear and elastic structural design problem is considered. The topology optimization problem is compliance minimization subject to a volume constraint. The mathematical description of the problem is described as below [56]: 4

1 𝑇 1 𝑓 𝑢 = 𝑢𝑇 𝐾𝑢 2 2 n ∑i=1 Vi xi Subject to ∶ 𝑉𝑓 ∗ − n =0 ∑i=1 Vi Minimize

𝐶=

(1)

𝑥𝑖 = 𝑥𝑚𝑖𝑛 𝑜𝑟 1, where 𝐶, 𝑓, 𝑢, 𝐾 are the objective function, global force vector, displacement vector, and global stiffness matrix, respectively. 𝑥𝑖 is the 𝑖 th design variable with candidate values of either 1 for the solid elements or prescribed 𝑥𝑚𝑖𝑛 for the void elements. Here 𝑉𝑓 ∗ is the target volume fraction. In general, the volume fraction of material decreases from 100% to the target iteratively.

𝛼𝑖 =

ro of

As a gradient-based method, the design variable 𝑥𝑖 of element 𝑖 updates based on its sensitivity 𝛼𝑖 . The gradient of compliance with respect to the design variable is derived from the adjoint method [57]: 𝜕𝐶 1 𝜕𝐾𝑖 1 𝑝−1 = − 𝑢𝑖𝑇 𝑢𝑖 = − 𝑝𝑥𝑖 𝑢𝑖𝑇 𝐾𝑖 𝑢𝑖 , 𝜕𝑥𝑖 2 𝜕𝑥𝑖 2

(2)

re

-p

where 𝐾𝑖 and 𝑢𝑖 are the stiffness matrix and displacement of the element, respectively. The exponent 𝑝 is the penalization power from the material interpolation scheme (typically 𝑝 = 3 [13]). This raw sensitivity is usually processed in order to produce a mesh-independent solution [58,59]. The elementary sensitivity filter scheme is invoked,

lP

∑𝑛𝑗=1 𝑊𝑗 𝛼𝑗 𝛼̅𝑖 = 𝑛 . ∑𝑗=1 𝑊𝑗

(3)

na

The weight function parameter 𝑊𝑗 is defined as 𝑊𝑗 = max(0, 𝑟𝑓 − 𝑑𝑖𝑗 ),

(4)

Jo

ur

where the parameter 𝑟𝑓 is the filter radius and 𝑑𝑖𝑗 is the distance between the centers of elements 𝑖 and 𝑗. In order to achieve a convergent solution, it is suggested to further average the sensitivity with its historical data. Averaging the sensitivity of the current iteration with that of the previous iteration is used [60] 1 (5) 𝛼̃ = (𝛼̅𝑖𝑘 + 𝛼̅𝑖𝑘−1 ). 2 Then the elements update scheme is devised based on the optimality criteria for the soft-kill BESO method [57]. In an iteration, a threshold is used to evaluate whether the elements shall be changed. The threshold is determined by the relative ranking of sensitivity and target volume of next iteration 𝑉 𝑘+1 . The 𝑉 𝑘+1 is defined based on the current volume 𝑉 𝑘 and the evolutionary ratio 𝑒𝑟:

5

(6) 𝑉 𝑘+1 = 𝑉 𝑘 (1 ± 𝑒𝑟). The design variables of solid elements are switched from 1 to 𝑥𝑚𝑖𝑛 if their sensitivities are lower than the threshold, and the design variables of void elements are switched from 𝑥𝑚𝑖𝑛 to 1 if the sensitivities are higher than the threshold. The optimization algorithm is implemented by Python and linked to Abaqus. The validity and convergence of the basic BESO computational framework have been analyzed and verified [61].

2.2. Structural connectivity

Jo

ur

na

lP

re

-p

ro of

The structural connectivity is related to the distribution of enclosed voids in the structure. Structures without enclosed voids are referred to as simply connected structures, while structures with enclosed voids are multiply-connected. In order to create designs without enclosed voids, the structural connectivity is evaluated first. In the BESO method, the elementbased definition is used to characterize structures. As illustrated in Fig. 1a and 1b, the blue elements are solid, whose design variables are 1, and the others are void, whose design variables are 𝑥𝑚𝑖𝑛 . The connected-component labelling algorithm [62,63] is integrated into the element-based structural definition to evaluate structural connectivity.

6

ro of -p re lP na

ur

Fig. 1 Illustration of the BESO representation of the structural connectivity: (a) an optimized structure, (b) the element-based model and (c) its structural connectivity.

Jo

The proposed SCC approach categorizes the elements into different sets, according to the design variables of both the target element and its neighbors. As shown in Fig. 1c, the elements are determined to belong to a solid elemental set, i.e., the blue part, and void elemental sets, i.e., the other parts. Face-connectivity or edge-connectivity can be used to define the connectivity of the elements. For example, in 3D, face-connectivity means connected elements share a surface, whereas, with edge-connectivity, they share one edge or two nodes. These definitions can be not only used for pixels or voxels, but also for irregular mesh in the finite element model. The algorithm labels whether two neighboring elements belong to the same set according to their design variables. Fig. 2 shows the example and the decision tree of the algorithm used for labelling an element, where element a, b, d and e are adjacent to target element c. In the 7

na

lP

re

-p

ro of

decision tree, 1 indicates that element c and its neighboring element have the same design variable value, and 0 indicates that they have different values. The labeled and unlabeled elements are highlighted in red and blue, respectively. The algorithm first compares element c with the labelled element a. Elements c and a belong to the same set if they have the same design variable value. Otherwise, the algorithm compares element c with other elements and labels them as a new set. This algorithm can be adapted for identifying the structural elemental sets based on the binary result generated by the BESO method, without any post-processing of the results.

ur

Fig. 2 Decision tree of elemental sets labelling for 2D edge-connectivity neighboring elements.

Jo

To find the enclosed voids, the algorithm determines if the void part is connected to the boundary. If so, this void part is identified as the non-closed void part, and if not, it is an enclosed void. The result of the structural connectivity is shown in Fig. 1c, the purple parts represent the non-closed void parts in the structure, and elements in other colors (yellow, red and green) are enclosed voids.

2.3. Topology modelling with hierarchical graphs In order to find the paths that connect the enclosed voids with the structural boundary, the element-based model needs to be converted into graphs. However, if millions of elements are 8

converted into millions of vertices in the graph, the computational cost of the algorithm is unacceptably high. In order to improve the computational efficiency of the SCC approach, we use a hierarchical graph scheme, in which the vertex in a high-level graph represents a set of elements, and the vertex in the lowest level graph represents an element. Two-level graphs are used in the following examples. The proposed approach in 3D can be easily extended to 2D cases. The difference between 2D and 3D cases is that only planar tunnels can be generated in 2D, but spatial tunnels can be generated in 3D. In this work, the approach mainly focuses on 3D cases. The 2D example is used for illustration purpose only.

ro of

The hierarchical graphs are built based on structural connectivity. The void elemental sets are identified to the enclosed-voids sets 𝑆e or the boundary sets 𝑆Γ . The solid elemental set is divided into subsets of elements, which are identified as the boundary sets 𝑆Γ or the solid interior sets 𝑆Ω . An example of the division of the solid elemental set is shown in Fig. 3a, where the solid elemental set is divided by the grey lines. Obviously, the coarser the division of the solid elemental set is, the fewer vertices will be in the level-1 graph, and the more vertices will be in the level-2 graph. To balance the efficiency of each graph, it is suggested that the number of vertices 𝑛𝑣 in each level can be set as 𝑛

𝑙 𝑛𝑣 = √ 𝑛𝑒 ,

(7)

-p

where 𝑛𝑙 and 𝑛𝑒 are the numbers of the level of the hierarchical graphs and the number of the elements, respectively. Based on the elements sets, the level-1 graph can be built as follow:

re

Given an undirected graph 𝐺 (1) = (𝑉 (1) , 𝐸 (1) ) with the vertices set 𝑉 (1) and the edges set 𝐸 (1) (1)

as the level-1 graph model. Each vertex represents the elemental sets. The edge 𝑒𝑖,𝑗 ∈ 𝐸 (1) (1)

lP

represents the connection between the subsets 𝑖 and 𝑗. For each edge 𝑒𝑖,𝑗 , the weight function (1)

(1)

𝑤𝑖,𝑗 is associated with the weight of the relative vertex 𝑣𝑖 (1)

(1)

(1)

(1)

The factor 𝜙𝑖

(1)

(8)

is defined as, (1)

1 = { 0.5

if 𝑣𝑖

(1)

(1)

∈ 𝑉Γ

(1)

ur

(1)

𝜙𝑖

(1)

+ 𝜙𝑗 𝑤𝑗 .

na

𝑤𝑖,𝑗 = 𝜙𝑖 𝑤𝑖

(1)

and 𝑣𝑗 , formulated as

0

(1) 𝑉Γ and

(9)

if 𝑣𝑖 ∈ 𝑉Ω otherwise,

(1)

Jo

where 𝑉Ω are the boundary vertices and solid interior vertices in the level-1 graph. The factor of enclosed-voids vertices is set as 0 because it has no effect on the generation of (1)

tunnels. To minimal the cost of the structural performance, the weight 𝑤𝑖 sum of the sensitivity of the elements in the vertex

is donated by the

(1) 𝑣𝑖 :

𝑛𝑖 (1) 𝑤𝑖

= ∑𝛼 ̃𝑘 ,

(10)

𝑘=1

9

Jo

ur

na

lP

re

-p

ro of

where 𝑛𝑖 is the number of elements in this vertex. The level-1 graph is shown in Fig. 3b, the blue and purple vertices are the boundary vertices, and the cyan vertices are the solid interior vertices. The yellow, red and green vertices represent different enclosed-voids vertices. There is no edge between the boundary vertices because they are the end of the paths.

Fig. 3 Illustration of the topology modeling with the hierarchical graph: (a) division of the solid elemental sets, (b) level-1 graph and (c) an example of converting a level-1 vertex to level-2 vertices. 10

The elements of the level-1 solid vertex would be converted to the vertices in the level-2 graph. For example, the level-1 vertex (the red circle in Fig. 3b) is converted to level-2 vertices. In Fig. 3c, each of the vertices in the level-2 graph is an element. As all elements in the level-1 vertex are interior solid elements, they are highlighted by the same color. In the level-2 graph 𝐺 (2) = (𝑉 (2) , 𝐸 (2) ) , the vertices represent solid elements or enclosed voids. The weight (2)

(2)

(2)

function 𝑤𝑖,𝑗 of the edge 𝑒𝑖,𝑗 , is associated with the weight of the vertex 𝑣𝑖 (2)

(2)

(2)

(2)

Similarly, the factor 𝜙𝑖

(2)

𝜙𝑖

(2)

1 = { 0.5

if 𝑣𝑖

(2)

(11)

is defined as, (2)

∈ 𝑉Γ

(2)

(12)

if 𝑣𝑖 ∈ 𝑉Ω otherwise,

0 (2)

(2)

+ 𝜙𝑗 𝑤𝑗 .

ro of

(2)

𝑤𝑖,𝑗 = 𝜙𝑖 𝑤𝑖

(2)

and 𝑣𝑗 :

(2)

where 𝑉Γ and 𝑉Ω are the boundary vertices and solid interior vertices in the level-2 graph. To control the radius of the tunnels, the sensitivity of elements within a given radius of the

𝑟 ∑𝑛𝑘=1 𝛼 ̃𝑘 = , 𝑛𝑟

re

(2) 𝑤𝑖

(2)

of vertex 𝑣𝑖

-p

(2)

element 𝑖 are taken into account. The weight 𝑤𝑖

is defined as

(13)

lP

where 𝑛𝑟 is the number of elements within a given radius of the element 𝑖.

na

2.4. Pathfinding scheme for eliminating enclosed voids

Jo

ur

Fig. 4a shows a 3D structure with three enclosed voids. In the SCC approach, the graph theory and BESO technique are combined to eliminate the enclosed voids by creating tunnels that connect them to the structural boundary. The pathfinding scheme is proposed to find and evaluate the possible paths along which the tunnels will be built. Since the enclosed voids can also be eliminated by connecting to other enclosed voids, the algorithm will also find paths between enclosed voids (Fig. 4b). To minimize the influence of tunneling on the optimized material layout, only one tunnel is adopted for each enclosed void, e.g., the yellow tunnels in Fig. 4b. Fig. 5 shows the flowchart of the SCC approach. The proposed approach is able to impose the structural connectivity constraint without significantly reducing the structural performance.

11

ro of -p re lP na ur

Jo

Fig. 4 Illustration of the pathfinding scheme for the elimination of enclosed voids: designs (a) without and (b) with tunnels.

12

ro of -p re lP na ur Jo

Fig. 5 Flowchart of the SCC approach.

Dijkstra’s shortest path algorithm [64,65] is here introduced based on the hierarchical graphs. Take the graph in Fig. 3b as an example, the SCC approach finds the level-1 tunnels connecting (1)

the enclosed-voids vertices to the boundary. For each enclosed-voids vertex 𝑣𝑒𝑖 , the SCC

13

approach traverses all boundary vertices in the level-1 graph to find all the shortest paths (1)

(1)

between them. Each pathfinding from 𝑣𝑒𝑖 to a boundary vertex 𝑣𝜏∗ can be described as (1)

(1)

(1)

(1)

(1)

(1)

(1)

(1)

(1)

Find Path 𝑃(𝑣𝑒𝑖 → 𝑉𝜏 ) = (𝑣𝑒𝑖 , 𝑣1𝑠𝑡 , 𝑣2𝑛𝑑 , … , 𝑣𝑙 , 𝑣𝜏∗ ), 𝑣𝜏∗ ∈ 𝑉𝜏 𝑙−1 (1)

(14)

(1)

Minimize ∑ 𝑤𝑘,𝑘+1 + 𝑤𝑝,𝜏∗ , 𝑘=1 (1)

(1)

(1)

where 𝑣1𝑠𝑡 , 𝑣2𝑛𝑑 and 𝑣𝑙 represents the first vertex, the second vertex and the last solid interior vertex in the path, respectively. Among all these paths, the shortest one of them is the (1)

path of the level-1 tunnel from the vertex 𝑣𝑒𝑖 to the boundary. (1)

(1)

and 𝑣𝑒𝑗 can be described as (1)

(1)

(1)

(1)

(1)

(1)

(1)

ro of

Similarly, pathfinding of the shortest level-1 tunnels between the enclosed-voids vertex 𝑣𝑒𝑖

Find Path 𝑃 (𝑣𝑒𝑖 → 𝑣𝑒𝑗 ) = (𝑣𝑒𝑖 , 𝑣1𝑠𝑡 , 𝑣2𝑛𝑑 , … , 𝑣𝑙 , 𝑣𝑒𝑗 ) 𝑙−1

(15)

-p

Minimize ∑ 𝑤𝑘,𝑘+1 . 𝑘=1 (1)

(1)

represent the first, second, and last vertex in this path,

re

(1)

Here, the vertex 𝑣1𝑠𝑡 , 𝑣2𝑛𝑑 and 𝑣𝑙 respectively.

lP

The pathfinding result of the level-1 tunnels is shown in Fig. 6a. Each enclosed void has two or more paths which connect them to the structural boundary. To reduce the influence on the material distribution and the structural performance caused by tunneling, only one path is adopted for each enclosed cavity.

na

In the SCC approach, the MST algorithm is used to remove the unnecessary tunnels. To implement this algorithm, a graph 𝐺 (𝑚𝑠𝑡) is built for describing the result of the level-1 shortest pathfinding, which can be defined as follow:

ur

Given an undirected graph 𝐺 (𝑚𝑠𝑡) = (𝑉 (𝑚𝑠𝑡) , 𝐸 (𝑚𝑠𝑡) ) with the vertices set 𝑉 (𝑚𝑠𝑡) and the edges set 𝐸 (𝑚𝑠𝑡) . In graph 𝐺 (𝑚𝑠𝑡) , a vertex represents the structural boundary and other vertices (𝑚𝑠𝑡)

represent the enclosed voids. The vertices are connected to each other. The edge 𝑒𝑖,𝑗

Jo

𝐸 (𝑚𝑠𝑡) represents the connection between the vertex (𝑚𝑠𝑡)

𝑤𝑖,𝑗

(𝑚𝑠𝑡)

of edge 𝑒𝑖,𝑗 (1)

vertex 𝑣𝑗

(𝑚𝑠𝑡) 𝑣𝑖

and

(𝑚𝑠𝑡) 𝑣𝑗 .

∈

The weight function (1)

is equal to the total weight of the shortest path between vertex 𝑣𝑖

and

in the level-1 graph.

The MST graph 𝐺 (𝑚𝑠𝑡) is illustrated in Fig. 6b, the blue vertex represents the boundary of the structure, and the yellow, red and green vertices represent different enclosed voids. Based on the graph 𝐺 (𝑚𝑠𝑡) , the MST algorithm generates a subgraph in which the sum of the weights of the tunnels is minimized. The tunnels in the subgraph will be reserved, and the others will be deleted. Besides, this subgraph must satisfy two conditions as follows. First, if there are 𝑚 14

vertices in 𝐺 (𝑚𝑠𝑡) , then the subgraph has 𝑚 − 1 edges. Second, each vertex in the subgraph has at least one edge connected to it. The MST algorithm can be defined as: (𝑚𝑠𝑡)

Find Graph 𝐺𝑠𝑢𝑏

(𝑚𝑠𝑡)

= (𝑉 (𝑚𝑠𝑡) , 𝐸𝑠𝑢𝑏 )

𝑚 𝑚−𝑖 (𝑚𝑠𝑡) (𝑚𝑠𝑡) 𝑥𝑖,𝑗

Minimize ∑ ∑ 𝑤𝑖,𝑗 𝑖

𝑗

𝑚 𝑚−𝑖 (𝑚𝑠𝑡)

Subject to: ∑ ∑ 𝑥𝑖,𝑗 𝑖

=𝑚−1

(16)

𝑗

𝑚−𝑖 (𝑚𝑠𝑡)

> 0, ∀𝑣𝑖

∈ 𝑉 (𝑚𝑠𝑡)

𝑗 (𝑚𝑠𝑡)

𝑥𝑖,𝑗 (𝑚𝑠𝑡)

where 𝑥𝑖,𝑗 (𝑚𝑠𝑡)

set 𝐸𝑠𝑢𝑏

= 0 𝑜𝑟 1, (𝑚𝑠𝑡)

is the variable for determining whether 𝑒𝑖,𝑗 (𝑚𝑠𝑡)

(𝐸𝑠𝑢𝑏

ro of

(𝑚𝑠𝑡)

∑ 𝑥𝑖,𝑗

∈ 𝐸 (𝑚𝑠𝑡) belongs to the sub-edges

⊆ 𝐸 (𝑚𝑠𝑡) ), 𝑚 is the number of vertices 𝑉 (𝑚𝑠𝑡) . The proposed approach is

Jo

ur

na

lP

re

-p

able to generate more than one tunnels by searching multiple paths of level-1 tunnels, which should be retained in the results of the MST algorithm.

15

ro of -p re lP na ur Jo Fig. 6 (a) Pathfinding result of the level-1 graph, (b) the MST graph, (c) the level-2 elementbased model, (d) the level-2 graph and (e) the optimized design with tunnels. 16

The red edges in Fig. 6b represent the paths of tunnels which will be retained, and the other edges represent the paths that need to be deleted. These reserved tunnels in the level-1 graph are considered as the possible region of the level-2 tunnels. The subsets of the elements, which corresponding to these tunnels, are illustrated in Fig. 6c. The grey elements in the level-2 elements-based model will not be converted to vertices in the level-2 graph. Obviously, the number of vertices in the level-2 graph is significantly reduced. Fig. 6d shows the level-2 graph defined in the level-2 elements-based model.

(2)

(2)

(2)

(2)

(2)

(2)

(2)

ro of

Then the SCC approach finds the level-2 paths which have minimal structural performance cost. Here, the enclosed-voids vertices are the source of the paths in the MST result, and the boundary vertices are corresponding to the end of these paths. Similarly, the level-2 pathfinding problem can be defined as (2)

(2)

Find Path 𝑃(𝑣𝑒𝑖 → 𝑉𝜏 ) = (𝑣𝑒𝑖 , 𝑣1 , 𝑣2 , … , 𝑣𝑝 , 𝑣𝜏∗ ), 𝑣𝜏∗ ∈ 𝑉𝜏 𝑝−1 (2)

(17)

(2)

Minimize ∑ 𝑤𝑘,𝑘+1 + 𝑤𝑝,𝜏∗ . 𝑘=1

-p

The red edges in Fig. 6d represent the paths of the level-2 tunnels.

lP

re

It should be noted that the SCC approach would change the sensitivities of the elements in the level-2 tunnels rather than delete them. Thus, the BESO technique will remove these elements and add elements in other places to satisfy the volume constraint. The sensitivities of the elements 𝑒𝑡 of the level-2 tunnels, including the elements that are represented by the blue vertices inside the red circle and their neighboring elements within the given radius, are changed in the optimization process. (2)

(2)

where the weight 𝑤𝑖 1 0

is from the Eq. (13), and

𝑖 ∈ 𝑒𝑡 𝑖 ∉ 𝑒𝑡 .

ur

𝜂𝑖 = {

(18)

na

𝛼̃𝑖 = 𝛼̃𝑖 − 𝜂𝑖 𝑤𝑖 ,

(19)

Jo

The final design of the example is illustrated in Fig. 6e. Remarkably, the tunnels could be generated in any iteration. However, to maximize optimization objective, it is suggested to start tunneling when the targeted volume fraction has been reached. The optimization is completed if the result has no enclosed voids and the convergence criterion is achieved. The convergence criterion is defined as [66] 𝑘−𝑡+1 |∑𝑁 − 𝐶 (𝑘−𝑁−𝑡+1) )| 𝑡=1(𝐶 ≤ 𝜖, (𝑘−𝑡+1) ∑𝑁 𝑡=1 𝐶

(20)

17

where 𝑘 is the current iteration number, 𝑁 is an integral number, 𝜖 is convergent error and 𝐶 (𝑘) is the compliance of the 𝑘 th iteration.

3. Numerical examples 3.1. A platform

Jo

ur

na

lP

re

-p

ro of

In the first example, a platform subjected to uniform pressure 𝑃 = 1 × 10−6 on the top surface is considered. It is pinned on the red region of the bottom surface, as shown in Fig. 7. A quarter structure is discretized into 40 × 40 × 30 elements for finite element analysis. The element size is 2 × 2 × 2. Three layers of elements on the top of the structure are chosen as a nondesign domain. The volume constraint 𝑉𝑓 = 30%.

Fig. 7 Design domain and boundary conditions of the platform.

The optimization result without SCC is shown in Fig. 8. The sectional view shows that there are six cavities enclosed in the platform. The design without enclosed voids is further shown in Fig. 9, where the radius of the tunnels is 6. Compared with the optimized design in Fig. 8, some parts, e.g., the X-shaped part, are broken for connecting the central hole with the other 18

Jo

ur

na

lP

re

-p

ro of

holes. The middle parts in the left and right sides of the platform (Fig. 8) are broken to connect the inside of the structure with its boundary.

19

ro of -p re lP

Jo

ur

na

Fig. 8 Optimized structure without SCC: (a) finite element model and (b) smoothed model.

20

ro of -p re lP na ur Jo

Fig. 9 Optimized design of the platform without enclosed voids.

Fig. 10 compares the evolutionary histories of the structural compliance with and without SCC. The compliances are normalized by that of the optimized result without SCC. A jump point is caused by tunneling (blue curve). Sometimes multiple tunnels, which are close to each other, are generated by the SCC approach. However, they may merge together into one tunnel in the following iterations of the BESO.

21

ro of

Fig. 10 Evolution history of the compliance.

Jo

ur

na

lP

re

-p

The finite element model of the 1-path design was smoothed and fabricated by a 3D printer using nylon powder (Fig. 11). The excess material can be easily removed from the prototype through the tunnels by using a cleaning system.

Fig. 11 3D printed model of the platform design without enclosed voids. 22

Jo

ur

na

lP

re

-p

ro of

The proposed approach is able to generate more tunnels by searching multiple shortest paths. For example, the enclosed void can be connected to the structural boundary by two tunnels, which are the first and second shortest tunnels. The symmetric design without enclosed voids is obtained and shown in Fig. 12, where the enclosed voids are required to connect with others through two tunnels. The X-shaped part is shown in Fig. 8 is broken into four parts. There are two exits at the left and right side of the structure for powder removal. By controlling the number of tunnels, competitive and diverse designs can be achieved. Besides, designs with multiple tunnels are easier to remove the powder than designs with a single tunnel, especially for complex designs.

Fig. 12 Optimized design of the platform with multiple tunnels. 23

The compliances of these two designs are about 0.6% and 2% larger than the design without SCC. Both the designs are structurally highly efficient.

3.2. A cantilever beam

na

lP

re

-p

ro of

A cantilever beam is optimized in this example (Fig. 13). A torque is applied on the free end. The structure is pin on the red surface as the illustration in the figure. The whole model is discretized into 40 × 40 × 120 elements and the element size is 0.5 × 0.5 × 0.5. Three layers of elements on both right and left surfaces are set as non-designable. The target volume fraction of material is 40%.

Jo

ur

Fig. 13 Design domain and boundary conditions of the cantilever beam.

24

ro of -p re lP na ur

Jo

Fig. 14 Optimized structure without SCC: (a) finite element model and (b) smoothed model.

The finite element model and the smoothed model of the optimized structure without SCC are shown in Fig. 14. There is a large enclosed cavity in this structure. Three cases with different tunnel numbers 𝐶1 = 1, 𝐶2 = 4 and 𝐶3 = 8 are considered, i.e., the enclosed hole is required to connect with the boundary through 1, 4 and 8 tunnels, respectively. The numbers of tunnels are set as 4 and 8 in order to generate symmetric designs. The radius of the tunnel is 2. The translucent views of the results are shown in Fig. 15. The volume constraint is satisfied in all cases. The 8-tunnel structure is approximately symmetric because of the computational error. Theoretically, the tunnel in the 1-path design could be generated at any position of the 8 tunnels in the 8-path result. However, the hierarchical graphs are ordered, which means that the 25

Jo

ur

na

lP

re

-p

ro of

algorithm would generate the tunnel by a certain sequence. Similarly, the tunnels in the 4-path case are generated by the same sequence. The compliance values of the resulting designs without enclosed voids are only 0.7%, 3%, and 6% higher than that of the design without SCC respectively. In work by Liu et al. [47], two designs without enclosed voids were created and the ratios of their structural rigidities to that of the unconstrained optimized design are as low as 0.39 and 0.03, indicating very low structural efficiency. For the same example, Zhou and Zhang [48] obtained three designs without enclosed voids. However, their compliance values increased by 72.56%, 13.07%, and 12.67%, respectively, compared to that of the original optimized solution. These results clearly demonstrate the superior performance of the proposed approach compared to existing methods.

26

ro of

Fig. 15 designs without enclosed voids of the cantilever beam: (a) 1-path, (b) 4-path, and (c) 8-path.

Jo

ur

na

lP

re

-p

The smoothed model of the 8-path design is shown in Fig. 15c. It was fabricated by a 3D printer using nylon powder. The printed structures are shown in Fig. 16. It is clearly seen from the left half of the structure that powders can be removed through the tunnels.

Fig. 16 3D printed model of the 8-path cantilever beam without enclosed voids.

27

4. Conclusions In this study, we have proposed a structural connectivity control approach based on the bidirectional evolutionary structural optimization technique. In this approach, the structural connectivity is evaluated during the optimization process. The shortest path is found to connect enclosed voids with the structural boundary. The enclosed voids are eliminated by generating tunnels along these paths. A hierarchical graph scheme is used to improve the computational efficiency of the approach.

ro of

In order to reduce the effect of generating tunnels on the structural performance, the proposed approach finds the paths which have the minimum influence on the objective function. In doing so, the developed methodology is capable of producing highly efficient structural designs which have no enclosed voids.

-p

Two examples are used to demonstrate the effectiveness of the algorithms. The results show that the proposed approach is capable of creating structurally efficient designs without enclosed voids. Furthermore, this technique can generate diverse and competitive designs by changing the radius and the number of tunnels.

re

CRediT Author Statement Yulin Xiong: Investigation, Writing – Original Draft

lP

Song Yao: Writing - Review & Editing, Validation

na

Zi-Long Zhao: Writing - Review & Editing, Validation Yi Min Xie: Conceptualization, Supervision

Jo

ur

Declaration of interests The authors declare that they have no known competing financial interests or personal relationships that could have appeared to influence the work reported in this paper.

Acknowledgements This work was supported by the Australian Research Council (DP160101400), the National Key R&D Program of China (2016YFB1200602-33), the National Natural Science Foundation of China (51778283, U1334208, 51275532), and the National Science and Technology Support Program (2015BAG13B01-05).

28

Reference Berman B. 3-D printing: the new industrial revolution. Bus Horiz 2012;55:155–62. doi:10.1016/j.bushor.2011.11.003.

[2]

Attaran M. The rise of 3-D printing: the advantages of additive manufacturing over traditional manufacturing. Bus Horiz 2017;60:677–88. doi:10.1016/j.bushor.2017.05.011.

[3]

Kruth JP, Mercelis P, Van Vaerenbergh J, Froyen L, Rombouts M. Binding mechanisms in selective laser sintering and selective laser melting. Rapid Prototyp J 2005;11:26–36. doi:10.1108/13552540510573365.

[4]

Harrysson OLA, Cansizoglu O, Marcellin-Little DJ, Cormier DR, West HA. Direct metal fabrication of titanium implants with tailored materials and mechanical properties using electron beam melting technology. Mater Sci Eng C 2008;28:366–73. doi:10.1016/j.msec.2007.04.022.

[5]

Botelho EC, Figiel, Rezende MC, Lauke B. Mechanical behavior of carbon fiber reinforced polyamide composites. Compos Sci Technol 2003;63:1843–55. doi:10.1016/S0266-3538(03)00119-2.

[6]

Peng Y, Wu Y, Wang K, Gao G, Ahzi S. Synergistic reinforcement of polyamide-based composites by combination of short and continuous carbon fibers via fused filament fabrication. Compos Struct 2019;207:232–9. doi:10.1016/j.compstruct.2018.09.014.

[7]

Rosen D. Design for additive manufacturing: past, present, and future directions. J Mech Des 2014;136:090301. doi:10.1115/1.4028073.

[8]

Salonitis K. Design for additive manufacturing based on the axiomatic design method. Int J Adv Manuf Technol 2016;87:989–96. doi:10.1007/s00170-016-8540-5.

[9]

Liu J, Ma Y. A survey of manufacturing oriented topology optimization methods. Adv Eng Softw 2016;100:161–75. doi:10.1016/j.advengsoft.2016.07.017.

[10]

Chahine G, Smith P, Kovacevic R. Application of topology optimization in modern additive manufacturing. 21st Annu Int Solid Free Fabr Symp - An Addit Manuf Conf SFF 2010 2010:606–18.

[11]

Bendsøe M. Optimal shape design as a material distribution problem. Struct Optim 1989;202:193–202.

[12]

Bendsøe MP, Kikuchi N. Generating optimal topologies in structural design using a homogenization method. Comput Methods Appl Mech Eng 1988;71:197–224. doi:10.1016/0045-7825(88)90086-2.

[13]

Bendsøe MP, Sigmund O. Material interpolation schemes in topology optimization. Arch Appl Mech (Ingenieur Arch 1999;69:635–54. doi:10.1007/s004190050248.

[14]

Bendsøe MP, Sigmund O. Topology optimization: theory, methods, and applications. Berlin, Heidelberg: Springer Berlin Heidelberg; 2003. doi:10.1007/978-3-662-05086-6.

[15]

Xie YM, Steven GP. A simple evolutionary procedure for structural optimization. Comput Struct 1993;49:885–96. doi:10.1016/0045-7949(93)90035-C.

Jo

ur

na

lP

re

-p

ro of

[1]

[16]

Xie YM, Steven GP. Optimal design of multiple load case structures using an evolutionary procedure. Eng Comput 1994;11:295–302. doi:10.1108/02644409410799290.

[17]

Huang, Xie YM, Burry. Advantages of bi-directional evolutionary structural optimization (BESO) over evolutionary structural optimization (ESO). Adv Struct Eng 2007;10:727–37. doi:10.1260/136943307783571436.

[18]

Huang X, Xie YM. A new look at ESO and BESO optimization methods. Struct Multidiscip Optim 2008;35:89–92. doi:10.1007/s00158-007-0140-4.

[19]

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.

29

Allaire G, Jouve F, Toader AM. Structural optimization using sensitivity analysis and a level-set method. vol. 194. 2004. doi:10.1016/j.jcp.2003.09.032.

[21]

Luo Z, Wang MY, Wang S, Wei P. A level set-based parameterization method for structural shape and topology optimization. Int J Numer Methods Eng 2008;76:1–26. doi:10.1002/nme.2092.

[22]

Tang J, Xie YM, Felicetti P. Conceptual design of buildings subjected to wind load by using topology optimization. Wind Struct 2014;18:21–35. doi:10.12989/was.2014.18.1.021.

[23]

Cavazzuti M, Baldini A, Bertocchi E, Costi D, Torricelli E, Moruzzi P. High performance automotive chassis design: a topology optimization based approach. Struct Multidiscip Optim 2011;44:45–56. doi:10.1007/s00158-010-0578-7.

[24]

Zhu JH, Zhang WH, Xia L. Topology optimization in aircraft and aerospace structures design. Arch Comput Methods Eng 2016;23:595–622. doi:10.1007/s11831-015-9151-2.

[25]

Aage N, Andreassen E, Lazarov BS, Sigmund O. Giga-voxel computational morphogenesis for structural design. Nature 2017;550:84–6. doi:10.1038/nature23911.

[26]

Das R, Jones R, Xie YM. Optimal topology design of industrial structures using an evolutionary algorithm. Optim Eng 2011;12:681–717. doi:10.1007/s11081-010-9132-0.

[27]

Luo Z, Yang J, Chen L. A new procedure for aerodynamic missile designs using topological optimization approach of continuum structures. Aerosp Sci Technol 2006;10:364–73. doi:10.1016/j.ast.2005.12.006.

[28]

Xie YM, Li Q, Querin OM, Steven GP. Shape and topology design for heat conduction by Evolutionary Structural Optimization. Int J Heat Mass Transf 1999;42:3361–71. doi:10.1016/s0017-9310(99)000083.

[29]

Haertel JHK, Engelbrecht K, Lazarov BS, Sigmund O. Topology optimization of a pseudo 3D thermofluid heat sink model. Int J Heat Mass Transf 2018;121:1073–88. doi:10.1016/j.ijheatmasstransfer.2018.01.078.

[30]

Zhao ZL, Zhou S, Feng XQ, Xie YM. On the internal architecture of emergent plants. J Mech Phys Solids 2018;119:224–39. doi:10.1016/j.jmps.2018.06.014.

[31]

Du J, Olhoff N. Minimization of sound radiation from vibrating bi-material structures using topology optimization. Struct Multidiscip Optim 2007;33:305–21. doi:10.1007/s00158-006-0088-9.

[32]

Du J, Yang R. Vibro-acoustic design of plate using bi-material microstructural topology optimization. J Mech Sci Technol 2015;29:1413–9. doi:10.1007/s12206-015-0312-x.

[33]

Gao J, Xue H, Gao L, Luo Z. Topology optimization for auxetic metamaterials based on isogeometric analysis. Comput Methods Appl Mech Eng 2019;352:211–36. doi:10.1016/j.cma.2019.04.021.

[34]

Wu J, Luo Z, Li H, Zhang N. Level-set topology optimization for mechanical metamaterials under hybrid uncertainties. Comput Methods Appl Mech Eng 2017;319:414–41. doi:10.1016/j.cma.2017.03.002.

[35]

Sigmund O. Materials with prescribed constitutive parameters:an inverse homogenization problem. J Solids Struct 1994;31.

Jo

ur

na

lP

re

-p

ro of

[20]

[36]

Huang X, Radman A, Xie YM. Topological design of microstructures of cellular materials for maximum bulk or shear modulus. Comput Mater Sci 2011;50:1861–70. doi:10.1016/j.commatsci.2011.01.030.

[37]

Sigmund O, Hougaard K. Geometric properties of optimal photonic crystals. Phys Rev Lett 2008;100. doi:10.1103/PhysRevLett.100.153904.

[38]

Zhou S, Townsend S, Xie YM, Huang X, Shen J, Li Q. Design of fishnet metamaterials with broadband negative refractive index in the visible spectrum. Opt Lett 2014;39:2415. doi:10.1364/OL.39.002415.

[39]

Sigmund O, Chigrin DN. Topology optimization in nano-photonics. AIP Conf Proc 2009;26:26–8. doi:10.1063/1.3253908.

30

Lazarov BS, Wang F. Maximum length scale in density based topology optimization. Comput Methods Appl Mech Eng 2017;318:826–44. doi:10.1016/j.cma.2017.02.018.

[41]

Lazarov BS, Wang F, Sigmund O. Length scale and manufacturability in density-based topology optimization. Arch Appl Mech 2016;86:189–218. doi:10.1007/s00419-015-1106-4.

[42]

Leary M, Merli L, Torti F, Mazur M, Brandt M. Optimal topology for additive manufacture: a method for enabling additive manufacture of support-free optimal structures. Mater Des 2014;63:678–90. doi:10.1016/j.matdes.2014.06.015.

[43]

Wang Y, Gao J, Kang Z. Level set-based topology optimization with overhang constraint: towards support-free additive manufacturing. Comput Methods Appl Mech Eng 2018;339:591–614. doi:10.1016/j.cma.2018.04.040.

[44]

Sigmund O. Manufacturing tolerant topology optimization. Acta Mech Sin 2009;25:227–39. doi:10.1007/s10409-009-0240-z.

[45]

Meisel N, Williams C. An investigation of key design for additive manufacturing constraints in multimaterial three-dimensional printing. J Mech Des 2015;137:111406. doi:10.1115/1.4030991.

[46]

Li Q, Chen W, Liu S, Fan H. Topology optimization design of cast parts based on virtual temperature method. CAD Comput Aided Des 2018;94:28–40. doi:10.1016/j.cad.2017.08.002.

[47]

Liu S, Li Q, Chen W, Tong L, Cheng G. An identification method for enclosed voids restriction in manufacturability design for additive manufacturing structures. Front Mech Eng 2015;10:126–37. doi:10.1007/s11465-015-0340-3.

[48]

Zhou L, Zhang W. Topology optimization method with elimination of enclosed voids. Struct Multidiscip Optim 2019. doi:10.1007/s00158-019-02204-y.

[49]

Zhao ZL, Zhou S, Cai K, Xie YM. A direct approach to controlling the topology in structural optimization. Comput Struct 2019;106141. doi:10.1016/j.compstruc.2019.106141.

[50]

Zhao ZL, Zhou S, Feng XQ, Xie YM. Morphological optimization of scorpion telson. J Mech Phys Solids 2020;135:103773. doi:10.1016/j.jmps.2019.103773.

[51]

Huang X, Xie YM, Burry MC. A new algorithm for bi-directional evolutionary structural optimization. JSME Int J Ser C 2006;49:1091–9. doi:10.1299/jsmec.49.1091.

[52]

Huang X, Xie YM. Bidirectional evolutionary topology optimization for structures with geometrical and material nonlinearities. AIAA J 2007;45:308–13. doi:10.2514/1.25046.

[53]

Huang X, Zuo ZH, Xie YM. Evolutionary topological optimization of vibrating continuum structures for natural frequencies. Comput Struct 2010;88:357–64. doi:10.1016/j.compstruc.2009.11.011.

[54]

Radman A, Huang X, Xie YM. Maximizing stiffness of functionally graded materials with prescribed variation of thermal conductivity. Comput Mater Sci 2014;82:457–63. doi:10.1016/j.commatsci.2013.10.024.

[55]

Huang X, Zhou S, Xie YM, Li Q. Topology optimization of microstructures of cellular materials and composites for macrostructures. Comput Mater Sci 2013;67:397–407. doi:10.1016/j.commatsci.2012.09.018.

Jo

ur

na

lP

re

-p

ro of

[40]

[56]

Huang X, Xie YM. Evolutionary topology optimization of continuum structures: methods and applications. vol. 91. Chichester, UK: John Wiley & Sons, Ltd; 2010. doi:10.1002/9780470689486.

[57]

Huang X, Xie YM. A further review of ESO type methods for topology optimization. Struct Multidiscip Optim 2010;41:671–83. doi:10.1007/s00158-010-0487-9.

[58]

Sigmund O. Morphology-based black and white filters for topology optimization. Struct Multidiscip Optim 2007;33:401–24. doi:10.1007/s00158-006-0087-x.

[59]

Sigmund O, Petersson J. Numerical instabilities in topology optimization: a survey on procedures dealing with checkerboards, mesh-dependencies and local minima. Struct Optim 1998;16:68–75. doi:10.1007/BF01214002.

31

Huang X, Xie YM. Convergent and mesh-independent solutions for the bi-directional evolutionary structural optimization method. Finite Elem Anal Des 2007;43:1039–49. doi:10.1016/j.finel.2007.06.006.

[61]

Zuo ZH, Xie YM. A simple and compact Python code for complex 3D topology optimization. Adv Eng Softw 2015;85:1–11. doi:10.1016/j.advengsoft.2015.02.006.

[62]

Fiorio C, Gustedt J. Two linear time Union-Find strategies for image processing. Theor Comput Sci 1996;154:165–81. doi:10.1016/0304-3975(94)00262-2.

[63]

Wu K, Otoo E, Shoshani A. Optimizing connected component labeling algorithms. Univ Calif 2005:1965. doi:10.1117/12.596105.

[64]

Dijkstra EW. A note on two problems in connexion with graphs. Numer Math 1959;1:269–71. doi:10.1007/BF01386390.

[65]

Karlsson RG, Poblete P V. An O(mlog logD) algorithm for shortest paths. Discret Appl Math 1983;6:91–3. doi:10.1016/0166-218X(83)90104-X.

[66]

Huang X, Xie YM. Bi-directional evolutionary topology optimization of continuum structures with one or multiple materials. Comput Mech 2009;43:393–401. doi:10.1007/s00466-008-0312-0.

Jo

ur

na

lP

re

-p

ro of

[60]

32