Selfsupporting Topology Optimization for Additive Manufacturing
Abstract
The paper presents a topology optimization approach that designs an optimal structure, called a selfsupporting structure, which is ready to be fabricated via additive manufacturing without the usage of additional support structures. Such supports in general have to be created during the fabricating process so that the primary object can be manufactured layer by layer without collapse, which is very timeconsuming and waste of material.
The proposed approach resolves this problem by formulating the selfsupporting requirements as a novel explicit quadratic continuous constraint in the topology optimization problem, or specifically, requiring the number of unsupported elements (in terms of the sum of squares of their densities) to be zero. Benefiting form such novel formulations, computing sensitivity of the selfsupporting constraint with respect to the design density is straightforward, which otherwise would require lots of research efforts in general topology optimization studies. The derived sensitivity for each element is only linearly dependent on its sole density, which, different from previous layerbased sensitivities, consequently allows for a parallel implementation and possible higher convergence rate. In addition, a discrete convolution operator is also designed to detect the unsupported elements as involved in each step of optimization iteration, and improves the detection process 100 times as compared with simply enumerating these elements. The approach works for cases of general overhang angle, or general domain, and produces an optimized structures, and their associated optimal compliance, very close to that of the reference structure obtained without considering the selfsupporting constraint, as demonstrated by extensive 2D and 3D benchmark examples.
keywords:
selfsupporting, topology optimization, explicit quadratic constraints, additive manufacturing, discrete convolution1 Introduction
Topology optimization aims to generate an optimal material distribution within a design domain under certain geometric or physical constraints. Since its introduction in late the 1980s Bendsoe88 , this problem has attracted wide industrial and academic interests due to its large potentiality in engineering applications and its intrinsic mathematical challenges. Topology optimization has developed in many different forms, such as: homogenization Bendsoe88 , density (SIMP) Bendsoe89 , evolutionary approaches (BESO) xie1993simple ; Huang2010 , level set Wang2003 ; VanDijk2013 , or more recently IGA (isogeometric analysis) Qian2013 , to name a few. See Sigmund2013 for a recent and comprehensive review on this topic.
The complex geometric designs produced by topology optimization show the approach’s superiority in balancing the geometric distribution and the target physical performance. Such designs are however very difficult to be manufactured directly via traditional subtractive or formative manufacturing techniques Sigmund2006On ; Serphos2014Incorporating ; Liu2016 . On the other hand, rapidly developing additive manufacturing technologies have the promise to overcome the barrier between the potentiality that the topology optimization approaches can provide and the limitations that traditional manufacturing technologies can fabricate. In reality, additive manufacturing is a natural counterpart to topology optimization in that they have very versatile capability to quickly generate and realize new components not existing before Lang20161 ; Lang20162 .
Despite the enhanced geometric freedom associated with additive manufacturing, specific design rules must still be satisfied in order to ensure manufacturability. The fabrication overhang angle is an example of a rule which is of paramount importance to ensure that the part will not collapse when fabricating the designed structure layer by layer. A structure satisfying such an overhang angle constraint is called selfsupporting. For example, Thomas Thomas2010 identified 45 degree as the typical maximum overhang angle with a large number of experiments. For a non selfsupporting structure, its geometry has to be modified or additional support structures need to be generated. Modifying the geometry will ultimately reduce the structure’s physical performance, while additional support raises the issue of automatic and minimum volume support design Majhi1999 ; Vanek2014Clever ; Wang2013 ; Dumas2014Bridging ; Wu2016 , and further postprocessing activities to remove the unwanted supports. In the case that the support is made of the same material as the main component, such as the selective laser melting (SLM) process using metals, it is extremely difficult to remove out the support structure. Particularly, when the generated supports are embedded within a closed volume of the model, it is impossible to remove them.
The best strategy to resolve the issue of topology optimization for additive manufacturing is perhaps to design a completely selfsupporting structure, via topology optimization, that can be fabricated directly without the usage of support materials. Brackett et al first suggested including the overhang angle constraints into the topology optimization process Brackett2011topology , but does not produce a complete selfsupporting structure. The first selfsupporting structure built from topology optimization is due to the pioneering work of Gaynor and Guest in 2014 Gaynor2014 (and a very recent journal version Gaynor2016 ), which is achieved via introducing a wedgeshaped filter during the topology optimization process. Also very recently, an excellent work was also conducted by Langelaar in 2D Lang20161 and in 3D Lang20162 via introducing a novel selfsupporting filter into the topology optimization process, which is achieved via building smooth approximation to the minimum and maximum functions. Impressive 2D and 3D examples were also shown in these studies Gaynor2016 ; Lang20161 ; Lang20162 .
In this paper, an alternative novel selfsupporting topology optimization approach is proposed to generate a structure of optimal physical performance that does not need any additional support materials. It is achieved via carefully formulating the selfsupporting constraint as an explicit quadratic function with respect to the design density, specifically, requiring the number of unsupported elements (in terms of the sum of their densities) to be zero. Benefiting from the novel quadratic formulation, the selfsupporting sensitivity for each element is straightforward to compute, and is only linearly dependent on density of the element itself; notice that designing a proper filter and computing the associated sensitivity usually requires lots of research efforts in general topology optimization framework Sigmund2013 ; Gaynor2016 ; Lang20161 ; Lang20162 . The derived sensitivity does not involve density information of any other elements, and thus allows for a parallel implementation, which is particularly important for 3D problem of high DOFs. Previous approaches Gaynor2016 ; Lang20161 ; Lang20162 have a nonlinear layerbased sensitivity expressions, and may thus inhibit parallel implementations, as also explained by the authors Gaynor2016 ; Lang20161 ; Lang20162 . In addition, a discrete convolution operator is also designed to detect the unsupported elements as involved in each step of optimization iteration, and improves the detection process 100 times as compared with simply enumerating these elements.
Comparisons between the proposed approach and previous studies Gaynor2014 ; Gaynor2016 ; Lang20161 ; Lang20162 are also summarized in Table 1.
Method 



Sensitivity  








MinMax 






Not needed 

The remainder of the paper is organized as follows. The novel formulation of selfsupporting topology optimization, together with its overall numerical procedure, is presented in Section 2. Several numerical techniques behind the proposed approach are detailed in Section 3. Extensive 2D and 3D examples are demonstrated in Section 4. pFinally the paper is concluded in Section 5.
2 Problem statement and approach overview
In the section, the selfsupporting constraint is formulated as a quadratic continuous function in terms of the element density, and is integrated within the classical SIMP framework Bendsoe89 for selfsupporting topology optimization. Following on from this, the proposed numerical approach to resolve the problem is outlined.
2.1 Supported and unsupported elements
The supported elements generally stand for the structural elements that can be fabricated via an additive manufacturing technology without collapse with respect to the fabrication process. They are defined here using the concepts of a maximum printable supporting angle, or overhang angle, which is first assumed to be 45 degree following previous study Brackett2011topology . Extensions of the approach to general overhang angles are also explained later. We also assume that the printing direction is following the positivey axis direction in both 2D and 3D for ease of explanation.
First consider a 2D discrete structured mesh model consisting of square elements , that is,
where are the indices increasing along the and axes respectively. Without confusion, we also use to represent a square element without explicitly mentioning its indices . In addition, a density matrix of size is also associated to , where an entry value or respectively represents a solid or void element of .
As illustrated in Fig. 1(a), given a solid element in (in orange), it is supported, or called a supported element, if one of the three blue elements below it is solid. We formulate the selfsupporting condition in a continuous form as follows: an element is supported if
Correspondingly, the supporting set of model is the set of all supported elements within , that is,
(1) 
Similarly, given a 3D structured mesh model consisting of cubic elements,
where are the indices increasing along the axes respectively, the supporting set of is similarly defined (see also Fig. 1(b)):
(2) 
Note also here that only five elements are included here as other element do not form an appropriate overhang angle with the orange element.
Correspondingly, the set of unsupported elements of a model is
(3) 
2.2 Formulation of selfsupporting topology optimization
The selfsupporting topology optimization problem aims to find the optimal material distribution within a design domain under certain boundary conditions. As widely studied before, the problem of minimum compliance or equivalently maximum stiffness is examined here. Following the classical SIMP framework Bendsoe89 , the problem of selfsupporting topology optimization is formulated here as an optimization problem with an additional explicit selfsupporting constraint. The constraint is carefully reformulated using a simple quadratic function with respect to the density, specifically, requiring the number of unsupported elements (in terms of the sum of square of their densities) is zero. Details are explained below.
The problem of selfsupporting topology optimization is stated as: find the optimal density distribution ,
(4) 
where is the vector of design variables (element densities) to be computed, is the vector of global displacements and is the global stiffness matrix. The objective function is the structure’s compliance, defined as
is the nodal force vector, and are the material volume and design domain volume, is the prescribed volume fraction, is the index set of unsupported elements as defined in (3), and is a small parameter closed to . A penalty parameter , usually set as , is applied here for the 0,1 convergence of , or specifically,
where is the element stiffness matrix associated with an element in the model and the associated element density.
The only difference of the above conventions in (4) with previous SIMPbased formulations is that it has an additional constraint to meet the selfsupporting requirement. This condition is based on the observation that when the sum of the element densities of unsupported elements tends to , all the elements are selfsupported; the square is used here so that its derivative is not constant. The simple quadratic expression allows for a straightforward sensitivity derivation of the selfsupporting constraint, and ultimately results in a linear expression.
2.3 Approach overview
The selfsupporting topology optimization problem (4) is readytosolve using the MMA approach noticing that the sensitivity is straightforward to compute (as can be further seen in Section 3.1). On the other hand, in order to further improve the approach’s convergence and computational efficiency, the overall optimization process is carefully designed, as plotted in Fig. 2 and detailed below.
Firstly, whether an element is selfsupporting is dependent on the printing direction. Different printing directions produce different optimization structures. Thus, if the printing direction is chosen arbitrarily, it may produce a structure totally different from the benchmark supportneeded structure obtained without considering the selfsupporting constraint, with a possible worse physical performance. In some very special case, the optimization approach may not converge. Thus, an appropriate printing direction is first set via generating a coarse structure via topology optimization without considering the selfsupporting constraints. It chooses the coordinate axis direction with the least number of unsupported elements as the print direction.
The criteria to generate this initial coarse structure, or to stop the initial optimization process stops, is based on the measure of nondiscreteness , proposed in Sigmund2007Morphology . It stops when
(5) 
where
and is the number of elements of the domain, is the density an element . The valve parameter is set , which corresponds to a structure of average density of or .
The generated coarse structure is then set as an initial structure for selfsupporting topology optimization, together with the chosen printing direction. During the optimization process, in order to prevent the selfsupporting constraint hindering the formation of loadcarrying structures, the constraint is imposed in a soft way using a similar strategy as previously performed in work Guo2014 . Specifically, the value of the tolerance parameter in (4) is decreased continuously from a relatively large value to the predefined tolerance during the course of optimization.
During the process of selfsupporting topology optimization, a blackwhite filter is additionally applied in the last few steps to produce a totally 01 density distribution so that the discrete convolution finds exactly the selfsupporting elements without the influence of gray elements. This is achieved using the Heaviside projection filter as originally designed by Guest et al Guest2004Filter ; Andreassen2011Efficient . Otherwise, gray elements may not be strong enough to support a black element above it, and special care has to be considered Lang20161 .
Lastly, in order to have a completely selfsupporting structure (without any unsupported elements), the selfsupporting constraint is added in a strict way in the last few iteration steps (when the number of unsupported elements is no longer decreasing) via removing nonself supporting elements (whose number is at most 5 is all tests given in this paper). This strategy has an ignorable influence on the final structure’s compliance, noticing the nature of optimization approach, i.e. the approach may fluctuate between selfsupporting constraints and target optimization and the volume faction constraint. Consequently, all these remaining elements are not essential in determining the structure’s physical compliance. This is very different from removing unsupported elements from the supportneeded structure in a postprocessing step.
3 Numerical aspects
3.1 Sensitivity analysis
The sensitivity of the selfsupport constraint is straightforward to compute from (4), and given as:
(6) 
where is the set of unsupported elements.
Note here that the sensitivity for each element is only dependent on the density of the element itself, without information of density of any other elements, and can be easily implemented in parallel. Previous approaches Gaynor2016 ; Lang20161 ; Lang20162 have shown skilled techniques in designing selfsupporting filters, and the selfsupporting sensitivity of each element was expressed in terms of the densities information of the layers below it. The approach has its great freedom in computing selfsupporting structure, but also requires a further improvement to overcome it inhabitation of parallel processing, as also explained by the authors Gaynor2016 ; Lang20161 ; Lang20162 .
Derivations of the sensitivities of the objective function or of the volume constraint , involved in (4), are totally the same as those done in previous studies Sigmund2001A ; Sigmund2013 . Integrating these sensitivities with thickness control can be achieved using the Heaviside filter Guest2004Filter , as will be demonstrated in Section 4. Details are not further explained here.
3.2 Discrete convolution for efficient unsupported element detection
Computing the sensitivity (6) of the selfsupporting constraints requires detecting the set of all the unsupported elements. They can be easily detected via enumerating all the discrete elements of not satisfying the property in (1) or (2). Such process is however very timeconsuming. In order to further accelerate this process, a novel convolution operator is designed for such detections, which acceleration the process of detecting the unsupported element with a speedup of 100 times as compared with the approach of directly enumerating them elementwise.
pGiven a discrete structure of size in 2D, we can see from (1) that an element is supported if the summation of the densities of its supporting elements is larger than zero, or specifically,
The newly introduced selfsupporting convolution operator is designed based on this observation. Specifically, suppose the overhang angle is . The associated 2D selfsupporting kernel matrix of size is defined in Fig. 3. A new matrix is then computed via performing the convolution between the density matrix and , or
(7) 
where is the convolution between matrices and , whose element is defined as
and the sign function
(8) 
Correspondingly, we have the set of supported elements of the discrete structure ,
(9) 
for defined in (7).
The above expression assumes a 01 density distribution of , while in the SIMP approach as studied here, the density matrix usually has entry value ranging from 0 to 1. A specific value usually needs to be set to replace 0 in the sign function (8) for practical applications.
The basic procedure of the convolution computation is further shown in Fig. 5 and explained below. For each element under consideration, a matrix centering at is selected. This is then followed by its Hadamard product, i.e. the element by element product between the matrices,
(10) 
where is the rotation of matrix at a degree of . The convolution value of element is the summation of all the values in the derived matrix .
The above procedure works for every element . For the boundary elements, an additional loop of void elements are added. Note also that the bottom elements are always taken as supported considering the fact that they are always supported by the baseboard of the fabrication device.
Once the set of supported elements of is determined from (9), the set of unsupported element is derived consequently,
(11) 
The convolution procedure in 3D is similar to that in 2D, and the corresponding kernel matrix is shown in Fig. 4. Extension of the approach to general overhang angle will be laterp explained in 3.3.
We also compare in Table 2 the computational time in detecting the supported elements using the proposed selfsupporting convolution operator and using direct enumeration element by element. Almost a 2 order of speedup is observed from the results. The time of detecting supported elements is not ignorable compared with FE computations. For example, one step of FE for the size of just takes 0.4163 seconds, and the enumeration time takes 0.5214 seconds. This is very important for the practical usage of the proposed approach, particularly on 3D complex structures with millions of elements, considering that FE computations can be implemented in parallel and accelerating the detection of supported elements then becomes important.
Domain size  Enumeration (s)  Convolution (s)  Speedup 

0.0094  0.0001  94  
0.0380  0.00035  106.16  
0.5214  0.0034  152.11  
0.1419  0.0010  139.13  
4.7974  0.0324  147.82  
114.7641  0.7605  150.88 
Note that various convolution operators have been designed and used as filters in topology optimization for design control, for example for removing checkerboard patterns or thickness control D1995Checkerboard . They usually aim to compute an element’s density or sensitivity via averaging those of the elements around it. In contrast to these researches, the convolution is used here for detecting unsupported elements. It does not change the element density.
3.3 Extension to general overhang angle
The above described procedure is mainly described for an overhang angle of for ease of explanations. Its extension to a general overhang angle is further explained below. The only difference is the construction of the kernel matrix involved in (10).
As shown in Fig. 6(a), given an overhang angle , a straight line passing through the center of an element and with a slope angle is drawn. Then the first element in each column whose centers are below line are taken as the supporting elements with respect to element . Their densities are set and others’ , which together determines a matrix . Rotating with respect to the center at a degree of gives the convolution kernel matrix , as shown in Fig. 6(b). The above procedure works for building convolution Kernel using multiple layers; the more layers taken the more accurate of the built kernel matrix for detecting the supported elements.
4 Examples
Extensive 2D and 3D examples are performed to test performance of the proposed approach. For illustration purposes, the material, load and geometry data are chosen to be dimensionless. The Young’s modulus and Poisson’s ratio of the solid material are set as and for all examples. The penalization factor is set to a value of . The minimal thickness is set to be 2, and the overhang angle is set to if not explicitly specified. The printing direction is selected during each optimization process and marked in the example figure. The 2D examples were implemented in Matlab, and the 3D examples were implemented in C++ and GPU for parallel computations on a computer of 3.2G CPU, 8.0G RAM and GeForce GTX 970 Graphic card.
The examples include the classical Cantilever beam, MBB, a 2D square example and three 3D examples. The Cantilever beam is used to illustrate various aspects of the approach: basic performance on a rectangular domain or a general domain, iteration process, thickness control, different volume fractions, different overhang angles, different types of external forces. The MBB shows the approach’s ability in handling constraints of multiple print directions. The 2D square example demonstrates the approach’s performance in case of complex topological structure obtained at distributed external forces. The 3D examples are further used to demonstrate the approach’s ability in handling complex 3D models of millions of DOFs via parallel implementation.
Following previous studies Gaynor2014 ; Gaynor2016 ; Lang20161 ; Lang20162 , we measure the ability of a selfsupporting topology optimization approach in maintaining the structure’s physical performance using the compliance ratio
(12) 
where is respectively the compliance of the structure computed with or without considering selfsupporting constraint.
The computational results for 2D examples are first summarized in Table 3; cases of 3D examples are explained later.
Beam  9000  24  92.7  92.8  100.11% 
Beam (hole)  7755  25  115.4  117.4  101.73% 
Beam ()  9000  24  92.7  92.8  100.11% 
Beam ()  9000  18  92.4  92.8  101.43% 
Beam ()  9000  12  92.6  92.9  100.32% 
Beam (concentrated)  14400  21  322.9  323.0  100.31% 
Beam (distributed)  14400  12  255.4  256.1  100.27% 
Beam (mixed)  14400  60  31971.1  32514.3  101.70% 
Beam (vf=0.6)  9000  24  92.7  92.8  100.11% 
Beam (vf=0.5)  9000  97  105.8  106.8  100.95% 
Beam (vf=0.4)  9000  37  127.4  128.3  100.71% 
Beam (vf=0.25)  9000  55  196.8  202.9  103.10% 
Beam (Angle=30)  9000  12  105.8  113.6  107.37% 
Beam (Angle=45)  9000  97  105.8  106.8  100.95% 
Beam (Angle=60)  9000  425  105.8  141.3  133.55% 
MBB  38400  1874  185.7  191.3  103.02% 
Square  22500  474  1312.8  1544.7  117.66% 
.
4.1 Cantilever beam example
The Cantilever beam, as shown in Fig. 7, is first tested. The model on the left has a rectangular domain, and has a target volume fraction of . The model on the right has general domain made via cutting a circular hole within the left one, and has a target volume fraction of . Both models are fixed on the left edge with an external force exerted on the middle point of its right edge. The print direction is determined from left to right.
Without considering the selfsupporting structure, the structure in Fig. 8(a),(c) are obtained where the elements in red are those that cannot be successfully printed out. The proposed selfsupporting topology optimization approach results in the structures in (b) and (d), both of which do not contain any unsupported elements. We can see from the results that the range containing unsupported elements in (a),(c) moves upward in (b),(d) to adapt the requirement of selfsupporting. In addition, it is also very interesting to notice that various parts of (b) or (d) are different from those of (a) or (c), for the structures’ maximal physical performance, although simultaneously maintaining their overall structures. The structures computed with or without selfsupporting constraints have a very close compliance, of a compliance ratio respectively of and for the left and right examples.
In handling the right model of a general design domain, we work on the rectangular domain following the procedure below. In each step of the optimization iteration, the density of each element within the circular domain is set , and then the convolution operation (detailed in Section 3.2) is performed in the whole rectangular domain to detect the unsupported elements. The above two steps are repeated until convergence.
4.1.1 Iteration performance
The iteration process of the example given in Fig. 7(a) is further explored by examining the variations of the structure’s topology, the number of unsupported elements and the compliance and volume fraction of the derived structure, as shown in Figs. 9 and 10.
The iteration process is divided into the following main steps (see also Fig. 9). Firstly, a topology optimization step without considering the selfsupporting constraint is performed, and results in the ”gray” structure in (a). After this, the relaxed selfsupporting constraint is added in the optimization iteration step, producing a structure in (b). The derived ”gray” structure is then transformed into a blackwhite structure using the Heaviside project filter, as given in (c). After this, the selfsupporting topology optimization process is iterated to reduce the number of unsupported elements while simultaneously optimizing its physical performance and maintaining its volume fraction, producing the structures in (d),(e),(f), and ultimately the final structure in (g). The unsupported elements are marked red in Figs. 9(c)(g), and illustrated in the caption. Their number is gradually decreased during the optimization iteration process.
Figure 10 shows the variation of the structure’s number of unsupported elements, target compliance and volume fraction in (a),(b) and (c) respectively. As can be seen, as the iteration step increases, the number of unsupported elements decrease until finally reaching zero. However, fluctuations in the number of unsupported elements may happen during the iteration process. The structures’ compliance and volume fraction decrease and finally reaches to a stable state.
4.1.2 Thickness control
The proposed approach is also able to control the structure’s thickness to meet different device requirements, with addition direct usage of a density filter for thickness control. Fig. 11 shows the obtained selfsupporting structures of thicknesses , and . The overall structures of these three selfsupporting structures are similar, and as the minimum thickness increases, the slender beams are removed gradually. Consequently, the smaller the thickness required the more details preserved in the final structures. The associated compliances of the three structures are very close to each other, as summarized in Table 3, of a ratio to their reference structures .
4.1.3 Different volume fractions
Performance of the approach is also tested at constraints of different volume fractions, and the computed structures are shown in . Such selfsupporting structures become harder to obtain for small value of volume fractions. As can be observed from the results, the selfsupporting constraints can still be satisfied although the number of elements decreases as the volume fraction becomes smaller.
4.1.4 Different overhang angles
As have been explained previously in Section 3.3, the proposed approach can also work for overhang angle different from via using different convolution Kernel matrices . We demonstrate its performance still using the Cantilever beam example in Fig. 22(a), for three different overhang angles: . The associated convolution kernel matrices for angles of are also shown in Fig. 13(d),(e). As can be observed from the examples, as the overhang angle becomes bigger, the boundary edges moves upward to meet the selfsupporting constraints. The structure’s minimal compliance is still well kept at these different angles, with respectively of , and . It is also noticed that the larger overhang angle deteriorates the structures’ physical performance. Similar phenomenon were also observed in previous studies Gaynor2016 ; Lang20161 ; Lang20162 .
4.2 2D Cantilever beam example with different types of forces
In order to test the ability of the proposed approach in selecting the print direction and its performance in finding the optimal structure at different external loadings, the 2D Cantilever beam example in Fig. 14 is tested under different types of forces, respectively of concentrated force, distributed force and mixed forces. In this example, the design domain is discretized into square FE elements. The volume fraction is and the minimum thickness is . The concentrated force is exerted on the middle point of the right edge and points downward. The distributed force is exerted evenly on the bottom, top and right edges of the model, while the case of mixed forces takes them both into account.
For each of the three cases, different print directions were chosen and shown in Fig. 14 by the proposed approach. The corresponding optimal structures are also shown in Figs. 15 respectively. As can be observed from the results and the summary in Table. 3, different boundary conditions may require different print directions and produce different optimized structures, which all can be handled successfully via the proposed approach. Compliance of the reference structure are maintained at a compliance ratio respectively of , and .
4.3 2D MBB example constrained by more than one print direction
The proposed approach is also able to simultaneously take into account more than one print direction constraints, provided they do not conflict with each other. This is illustrated using the classical MBB problem in Fig. 17. Due to the symmetry of the model’s structure and boundary conditions, only half of the computational domain is used here which consists of square FE mesh elements. The volume fraction is and the minimum thickness is .
The aim is to produce a selfsupporting structure maintaining the mirror symmetry of the original model. Thus the selfsupporting requirement has to be added in both directions: from right to left and from left to right for the half sized structure in Fig. 17(b). As a result, a selfsupporting in both directions is obtained in Fig. 17, as compared with the supportneeded structure in Fig. 17(a), where the unsupported elements are plotted in red. The compliance of the supportneeded structure and the selfsupporting structure are respectively 185.7 and 191.3, at a ratio of .
Note that the 2D MBB problem of the same domain size was also tested by Gaynor and Guest in Gaynor2016 , where the print direction was manually set from bottom to top. Such setting thus does not require constraints of multiple constraints.
4.4 Complex internal structure at distributed force
The proposed approach is also able to produce selfsupporting structure for complex internal structure, as demonstrated using the square example at distributed forces in Fig. 18(a). The computed supportneeded structure and selfsupporting structure are respectively given in Fig. 18(b) and (c). It can be seen from the results that the original supportneeded structure has many small flat edges which prevent the structure to be fabricated without a large number of additional supports. Such unsupported elements have disappeared in the optimized selfsupporting structure of Fig. 18(b), despite the structure’s high complexity. In addition, the selfsupporting structure also has a close compliance to that of the original supportneeded structure, respectively of and at a compliance ratio of . A selfsupporting structure is very necessary for such highly complex shape, as computing the supports or its removal would be extremely troublesome if not impossible.
4.5 3D examples
The proposed approach allows for a parallel implementation, which is particularly important for complex 3D examples of high DOFs. We have implemented the approach in GPU for parallel computations, and it works efficiently for examples of millions of elements within almost an hour. The test examples includes the classical benchmark examples: a 3D wheel, a 3D Cantilever, and a newly devised examples of a 3D desk. The domain sizes and their associated computational time is summarized in Table 4.
Example  Wheel in Fig. 19  Cantilever(a) in Fig. 19  Cantilever(b) in Fig. 19  Desk in Fig. 23p 

Size  
Time  35m  11m  17m  66m 
4.5.1 3D Wheel
The 3D wheel example in Fig. 19 consists of cubic mesh elements. The four corners of the bottom face are fixed and a concentrated force is exerted on the middle point of the bottom face. The target volume fraction is and the minimum thickness is . The print direction is chosen as top to bottom.
The computed selfsupporting structure using the proposed approach is shown in Fig. 20(b), as compared with its counterpart of supportneeded structure in Fig. 20(a). The corresponding structure slices at of both the selfsupporting and supportneeded structures are also shown and compared in Fig. 20(f)(g)(h). As can be seen from the results, the originally flat regions of the supportneeded structure, which cannot be fabricated without supports, have been optimized to meet the selfsupporting requirement. The resulting structure is totally selfsupporting for direct fabrication purpose and its compliance is , very close to that of the original supportneeded one of , of a ratio
4.5.2 3D Cantilever
Two different 3D Cantilever examples are tested here as illustrated in Fig. 22: one of size exerted by point loadings, and another one of size exerted by edge loadings; a same example to the latter was also studied in Lang20162 . The target volume fraction is and the minimum thickness is . The print direction is chosen as top to bottom. The computed selfsupporting structures are shown in Fig. 22, which is directly to be fabricated without any additional support materials.
4.5.3 3D Desk
A more complex 3D desk problem, as shown in Fig. 23, is designed to further test performance of the proposed approach. The example is of size , and consist of millions of elements. In this example, the four bottom corners are fixed and the top face is exerted by a uniform distributed force pointing downward. The target volume fraction is , the minimum thickness is and the print direction is chosen from right to left. The final generated selfsupporting structure using the proposed approach is shown in Fig. 23, taking minutes. It is also interesting to note that the four legs of the desk are not totally solid but take a porous bonelike structure to balance the constraint of the object weight and the target compliance.
4.6 Summary
As can be observed from these examples, using the proposed selfsupporting topology optimization approach, a completely selfsupporting structure is achieved that has a very similar geometric shape to that of the supportneeded one without considering the selfsupporting constraint. The edges or faces of the original supportneeded structure, are alighted toward the print direction so that all elements can be successfully fabricated. In addition, the difference in compliance between the selfsupporting structure and the supportneeded structure is maintained within a very small or negligible range. All these enhancements demonstrate the strength of the proposed approach in designing the selfsupporting structures and simultaneously maintaining their optimized structure and physical performance.
5 Conclusion
A novel selfsupporting topology optimization approach is provided in this paper applicable to additive manufacturing. The usage of convolution operator and the associated numerical techniques enables the selfsupporting structure to be reliably generated with high efficiency and robustness. Extensive 2D and 3D examples are provided to test the approach’s performance. The final derived structures are completely selfsupporting and have a compliance very close to the optimal supportneeded structure, proving its high effectiveness.
The proposed approach is at present implemented using regular square or cubic elements, which are dominated in researches of topology optimization. On the other hand, the overall framework also works for general domains consisting of irregular quad or hex elements, but is also limited by the fact that the convolution operator presented in Section 3.2 is no longer applicable as they become different for different elements. Thus the element by element enumeration has to be taken, and will reduce the computational efficiency.
The proposed approach chooses an appropriate printing direction before taking into account the selfsupporting constraint during the optimization process. The strategy aims to provide a good initial value to improve the convergence and to result in a possible smaller target compliance. On other hand, if ignoring the step of choosing the print angle and setting an arbitrary print direction, the approach may fail to produce a converged result. For example, an optimized selfsupporting structure is hard to obtain for the MBB example at a different direction. Such phenomenon may not prevent the approach in generating a selfsupporting structure suitable for additive manufacturing, but may hinder its usage for specific applications. It deserves further researches efforts to improve this.
The proposed approach can also be extended to porous interior designs of 3D freeform structures that do not need any additional supports within its interior. Such supports would otherwise be very difficult to remove. This novelty will be very useful in some additive manufacturing technologies and is currently an active area of research. Furthermore, besides the selfsupporting requirements, other constraints such as hanging bridge also need to be considered in the optimization process so that a designed structure can be robustly fabricated.
Acknowledgements
The work described in this paper is partially supported by the NSF of China (No. 61472356, 61210007).
References
References
 (1) M. Bendsoe, N. Kikuchi, Generating optimal topologies in structural design using a homogenization method, Computer Methods in Applied Mechanics and Engineering 71 (2) (1988) 197–224.
 (2) M. P. Bendsøe, Optimal shape design as a material distribution problem, Structural optimization 1 (4) (1989) 193–202.
 (3) Y. M. Xie, G. P. Steven, A simple evolutionary procedure for structural optimization, Computers & structures 49 (5) (1993) 885–896.
 (4) X. Huang, Y. M. Xie, A further review of ESO type methods for topology optimization, Structural and Multidisciplinary Optimization 41 (2010) 671–683.
 (5) M. Y. Wang, X. Wang, D. Guo, A level set method for structural topology optimization, Computer Methods in Applied Mechanics and Engineering 192 (2003) 227–246.
 (6) N. P. van Dijk, K. Maute, M. Langelaar, F. van Keulen, Levelset methods for structural topology optimization: a review, Structural and Multidisciplinary Optimization 48 (3) (2013) 437–472.
 (7) X. Qian, Topology optimization in Bspline space, Computer Methods in Applied Mechanics and Engineering 265 (2013) 15–35.
 (8) O. Sigmund, K. Maute, Topology optimization approaches: a comparative review, Structural and Multidisciplinary Optimization 48 (6) (2013) 1031–1055.
 (9) O. Sigmund, On topology Optimization with Manufacturing Constraints, Springer Netherlands, 2006.
 (10) M. R. Serphos, Incorporating AMspecific manufacturing constraints into topology optimization.
 (11) J. Liu, Y. Ma, A survey of manufacturing oriented topology optimization methods, Advances in Engineering Software 100 (2016) 161–175.
 (12) M. Langelaar, An additive manufacturing filter for topology optimization of printready designs, Structural and Multidisciplinary Optimization (2016) 1–13.
 (13) M. Langelaar, Topology optimization of 3D selfsupporting structures for additive manufacturing, Additive Manufacturing 12 (2016) 60–70.
 (14) T. Daniel, The development of design rules for selective laser melting, University of Wales.
 (15) J. Majhi, R. Janardan, J. Schwerdt, M. Smid, P. Gupta, Minimizing support structures and trapped area in twodimensional layered manufacturing, Computational Geometry 12 (34) (1999) 241–267.
 (16) J. Vanek, J. Galicia, B. Benes, Clever support: Efficient support structure generation for digital fabrication, Computer Graphics Forum 33 (5) (2014) 117–125.
 (17) W. Wang, T. Y. Wang, Z. Yang, L. Liu, X. Tong, W. Tong, J. Deng, F. Chen, X. Liu, Costeffective printing of 3D objects with skinframe structures, ACM Transactions on Graphics 32 (6) (2013) 1–10.
 (18) J. Dumas, J. Hergel, S. Lefebvre, Bridging the gap: Automated steady scaffoldings for 3D printing, Acm Transactions on Graphics 33 (4) (2014) 1–10.
 (19) J. Wu, C. L. Wang, X. Zhang, R. Westermann, Selfsupporting rhombic infill structures for additive manufacturing, ComputerAided Design (2016) 1–11.
 (20) D. Brackett, I. Ashcroft, R. Hague, Topology optimization for additive manufacturing, in: Proceedings of the Solid Freeform Fabrication Symposium, Austin, TX, 2011, pp. 348–362.
 (21) A. T. Gaynor, J. K. Guest, Topology optimization for additive manufacturing: Considering maximum overhang constraint, AIAA AVIATION 2014  15th AIAA/ISSMO Multidisciplinary Analysis and Optimization Conference (June) (2014) 1–9.
 (22) A. T. Gaynor, J. K. Guest, Topology optimization considering overhang constraints: Eliminating sacrificial support material in additive manufacturing through design, Structural & Multidisciplinary Optimization 54 (5) (2016) 1157–1172.
 (23) O. Sigmund, Morphologybased black and white filters for topology optimization, Structural & Multidisciplinary Optimization 33 (45) (2007) 401–424.
 (24) X. Guo, W. Zhang, W. Zhong, Explicit feature control in structural topology optimization via level set method, Computer Methods in Applied Mechanics and Engineering 272 (2014) 354–378.
 (25) J. Guest, J. Prevost, T. Belytschko, Achieving minimum length scale in topology optimization using nodal design variables and projection functions, International journal for numerical methods in engineering 61 (2) (2004) 238–254.
 (26) E. Andreassen, A. Clausen, M. Schevenels, B. S. Lazarov, O. Sigmund, Efficient topology optimization in matlab using 88 lines of code, Structural & Multidisciplinary Optimization 43 (1) (2011) 1–16.
 (27) O. Sigmund, A 99 line topology optimization code written in matlab, Structural & Multidisciplinary Optimization 21 (2) (2001) 120–127.
 (28) A. D az, O. Sigmund, Checkerboard patterns in layout optimization, Structural & Multidisciplinary Optimization 10 (1) (1995) 40–45.