Self-supporting Topology Optimization for Additive Manufacturing
Dengyang Zhaoa , Ming Lia,∗, Yusheng Liua
a State

Key Laboratory of CAD&CG, Zhejiang University, Hangzhou, China

arXiv:1708.07364v1 [] 24 Aug 2017

Abstract
The paper presents a topology optimization approach that designs an optimal structure, called a self-supporting 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 time-consuming and waste of material.
The proposed approach resolves this problem by formulating the self-supporting 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
self-supporting 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 layer-based 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 self-supporting
constraint, as demonstrated by extensive 2D and 3D benchmark examples.
Keywords: self-supporting, topology optimization, explicit quadratic constraints, additive manufacturing, discrete convolution

1. 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 [1], 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 [1], density (SIMP) [2], evolutionary approaches (BESO) [3, 4], level set [5, 6], or more recently IGA
(iso-geometric analysis) [7], to name a few. See [8] 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 [9, 10, 11]. 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 [12, 13].
∗ Corresponding

author: liming@cad.zju.edu.cn

Preprint submitted to Elsevier

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 self-supporting. For example, Thomas [14] identified 45 degree as the typical maximum
overhang angle with a large number of experiments. For a non
self-supporting 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 [15, 16, 17, 18, 19], and
further post-processing 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 self-supporting 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 [20],
August 25, 2017

but does not produce a complete self-supporting structure. The
first self-supporting structure built from topology optimization
is due to the pioneering work of Gaynor and Guest in 2014 [21]
(and a very recent journal version [22]), which is achieved
via introducing a wedge-shaped filter during the topology optimization process. Also very recently, an excellent work was
also conducted by Langelaar in 2D [12] and in 3D [13] via introducing a novel self-supporting 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 [22, 12, 13].
In this paper, an alternative novel self-supporting 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
self-supporting 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 self-supporting 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 [8, 22, 12, 13]. 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 [22, 12, 13] have
a nonlinear layer-based sensitivity expressions, and may thus
inhibit parallel implementations, as also explained by the authors [22, 12, 13]. 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 [21, 22, 12, 13] are also summarized in Table 1.
The remainder of the paper is organized as follows. The
novel formulation of self-supporting 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.

y
45°

y
0

0

x

x

z
(a) 2D case

(b) 3D case

Figure 1: A supported element (in orange) is supported by one of the supporting
elements in blue in 2D and 3D.

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 [20]. Extensions of the approach to general overhang angles are also
explained later. We also assume that the printing direction is
following the positive-y axis direction in both 2D and 3D for
ease of explanation.
First consider a 2D discrete structured mesh model M consisting of square elements e(n, m), that is,
M = {e(n, m)| 1 ≤ n ≤ N, 1 ≤ m ≤ M},
where n, m are the indices increasing along the x and y axes
respectively. Without confusion, we also use e to represent a
square element without explicitly mentioning its indices n, m.
In addition, a density matrix ρ of size N × M is also associated to M, where an entry value ρ(n, m) = 1 or 0 respectively
represents a solid or void element e(n, m) of M.
As illustrated in Fig. 1(a), given a solid element e(n, m) in
M (in orange), it is supported, or called a supported element, if
one of the three blue elements below it is solid. We formulate
the self-supporting condition in a continuous form as follows:
an element e(n, m) ∈ M is supported if
X
ρ(r, m − 1) > 0.
n−1≤r≤n+1

Correspondingly, the supporting set MS of model M is the
set of all supported elements within M, that is,
X
MS = {e(n, m) ∈ M|m = 1 or
ρ(r, m − 1) > 0}. (1)
n−1≤r≤n+1

Similarly, given a 3D structured mesh model M consisting
of cubic elements,
M = {e(n, m, l)| 1 ≤ n ≤ N, 1 ≤ m ≤ M, 1 ≤ l ≤ L},

2. Problem statement and approach overview

where n, m, l are the indices increasing along the x, y, z axes respectively, the supporting set of M is similarly defined (see also
Fig. 1(b)):
X
MS = {e(n, m, l) ∈ M|m = 1 or
ρ(r, m − 1, s) > 0}.

In the section, the self-supporting constraint is formulated as
a quadratic continuous function in terms of the element density, and is integrated within the classical SIMP framework [2]
for self-supporting topology optimization. Following on from
this, the proposed numerical approach to resolve the problem is
outlined.

n−1≤r≤n+1,
l−1≤s≤l+1

(2)
2

Table 1: Difference between self-supporting topology optimization approaches

Method
Gaynor and
Guest [21, 22]
Lang [12, 13]
The proposed

Supported Element
Detection
Layerwise
Enumeration
Layerwise
Enumeration

Self-supporting
Constraint Prescription
Area occupation

Convolution

Min-Max
Quadratic
function

2.2. Formulation of self-supporting topology optimization
The self-supporting 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 [2], the problem of self-supporting topology optimization
is formulated here as an optimization problem with an additional explicit self-supporting 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 self-supporting topology optimization is
stated as: find the optimal density distribution ρ,

















Not needed

Sensitivity
Nonlinear
Layer dependent
Nonlinear
Layer dependent
Linear
Not dependent

where K0e is the element stiffness matrix associated with an element e in the model M and ρe the associated element density.
The only difference of the above conventions in (4) with previous SIMP-based formulations is that it has an additional conP
straint U(ρ) = e∈U ρ2e ≤  to meet the self-supporting requirement. This condition is based on the observation that when the
sum of the element densities of unsupported elements tends to
0, all the elements are self-supported; the square is used here
so that its derivative is not constant. The simple quadratic expression allows for a straightforward sensitivity derivation of
the self-supporting constraint, and ultimately results in a linear
expression.

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
M is
MU = M \ M S .
(3)

min c(u, ρ), s.t.
ρ∈RN×M

Smooth
Approximation
Heaviside
function
Norm
approximation

2.3. Approach overview
The self-supporting topology optimization problem (4) is
ready-to-solve 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 self-supporting 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 support-needed structure obtained without considering the self-supporting 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 self-supporting
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 non-discreteness Mnd , proposed in [23]. It stops when

(4)

K(ρ)u = f(ρ),
V(ρ)/V0 ≤ f,
P
U(ρ) = e∈MU ρ2e ≤ ,
0 < ρe ≤ 1, e = 1, . . . , s,

where ρ is the vector of design variables (element densities) to
be computed, u is the vector of global displacements and K is
the global stiffness matrix. The objective function c(u, ρ) is the
structure’s compliance, defined as
c(u, ρ) = uT Ku.

0
Mnd < Mnd
,

f(ρ) is the nodal force vector, V(ρ) and V0 are the material volume and design domain volume, f is the prescribed volume
fraction, MU is the index set of unsupported elements as defined in (3), and  > 0 is a small parameter closed to 0. A
penalty parameter p, usually set as p = 3, is applied here for
the 0,1 convergence of ρ, or specifically,

where

(5)

4ρe (1 − ρe )
× 100%,
n
and n is the number of elements of the domain, ρe is the den0
sity an element e. The valve parameter Mnd
is set 0.36, which
corresponds to a structure of average density of 0.1 or 0.9.
PN

Mnd =

Ke = ρep K0e ,
3

e=1

N

Topology optimization
without self-support constraint

Initial
structure

Y

Determine the
print direction

Topology optimization
with self-support constraint

Self-support
structure

Figure 2: Main procedure of the overall approach.

The generated coarse structure is then set as an initial structure for self-supporting topology optimization, together with
the chosen printing direction. During the optimization process,
in order to prevent the self-supporting constraint hindering the
formation of load-carrying structures, the constraint is imposed
in a soft way using a similar strategy as previously performed
in work [24]. Specifically, the value of the tolerance parameter
 in (4) is decreased continuously from a relatively large value
to the predefined tolerance 0 during the course of optimization.
During the process of self-supporting topology optimization,
a black-white filter is additionally applied in the last few steps
to produce a totally 0-1 density distribution so that the discrete
convolution finds exactly the self-supporting elements without
the influence of gray elements. This is achieved using the
Heaviside projection filter as originally designed by Guest et
al [25, 26]. Otherwise, gray elements may not be strong enough
to support a black element above it, and special care has to be
considered [12].
Lastly, in order to have a completely self-supporting structure (without any unsupported elements), the self-supporting
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 non-self 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 self-supporting 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 support-needed structure in a post-processing step.

1

1

1

0

0

0

0

0

0

Figure 3: The 2D self-supporting convolution kernel matrix H to detect supported elements for an overhang angle of 45◦ .

have shown skilled techniques in designing self-supporting filters, and the self-supporting 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 [22, 12, 13].
Derivations of the sensitivities of the objective function (¸ρ)
or of the volume constraint V(ρ), involved in (4), are totally
the same as those done in previous studies [27, 8]. Integrating
these sensitivities with thickness control can be achieved using
the Heaviside filter [25], 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 self-supporting constraints requires detecting the set MU of all the unsupported elements. They can be easily detected via enumerating all the discrete elements e of M not satisfying the property in (1) or (2).
Such process is however very time-consuming. 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
element-wise.
pGiven a discrete structure M of size N × M in 2D, we can
see from (1) that an element e(n, m) ∈ M is supported if the
summation of the densities of its supporting elements is larger
than zero, or specifically,
X
ρ(r, m − 1) > 0.

3. Numerical aspects
3.1. Sensitivity analysis
The sensitivity of the self-support constraint U(ρ) is straightforward to compute from (4), and given as:


∂U 
2ρe if e ∈ MU
=
(6)
0
∂ρe 
if e < MU ,

n−1≤r≤n+1

where MU is the set of unsupported elements.
Note here that the sensitivity for each element e is only dependent on the density ρe of the element itself, without information of density of any other elements, and can be easily implemented in parallel. Previous approaches [22, 12, 13]

The newly introduced self-supporting convolution operator
is designed based on this observation. Specifically, suppose the
overhang angle is 45◦ . The associated 2D self-supporting kernel matrix H of size 3 × 3 is defined in Fig. 3. A new matrix ρS
4

0 0 0
0 0 0
0 0 0

The convolution procedure in 3D is similar to that in 2D, and
the corresponding kernel matrix H 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 self-supporting
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 600 × 400 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.

0 1 0
1 1 1
0 1 0
Figure 4: The 3D self-supporting convolution kernel matrix H to detect supported elements for an overhang angle of 45◦ .

is then computed via performing the convolution between the
density matrix ρ and H, or
ρ̃S = sign(ρ ∗ H),

(7)

where ρ ∗ H is the convolution between matrices ρ and H,
whose (n, m) element is defined as
(ρ ∗ H)(n, m) =

1 X
1
X

Table 2: Time-cost comparison between direct enumeration and using convolution for detecting supported elements in 2D and 3D cases

ρ(n − i, m − j) · H(i + 1, j + 1),

Domain size

i=−1 j=−1

80 × 40
320 × 80
600 × 400
40 × 40 × 20
100 × 100 × 100
500 × 200 × 100

and the sign function



0 if x = 0,
sign(x) = 

1 if x > 0.

(8)

Correspondingly, we have the set of supported elements MS of
the discrete structure M,
MS = {e(n, m) | ρ̃S (n, m) = 1},

ρe ,

Convolution (s)

Speedup

0.0094
0.0380
0.5214
0.1419
4.7974
114.7641

0.0001
0.00035
0.0034
0.0010
0.0324
0.7605

94
106.16
152.11
139.13
147.82
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 [28]. 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.

(9)

for ρ̃S defined in (7).
The above expression assumes a 0-1 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 e
under consideration, a 3 × 3 matrix ρe centering at e is selected.
This is then followed by its Hadamard product, i.e. the element
by element product between the matrices,
Ce = H̃

Enumeration (s)

3.3. Extension to general overhang angle
The above described procedure is mainly described for an
overhang angle of 45◦ 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 H involved
in (10).
As shown in Fig. 6(a), given an overhang angle θ, a straight
line L passing through the center of an element e and with a
slope angle θ is drawn. Then the first element in each column
whose centers are below line L are taken as the supporting elements with respect to element e. Their densities are set 1 and
others’ 0, which together determines a matrix H̄. Rotating H̄
with respect to the center e at a degree of 180◦ gives the convolution kernel matrix H, 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.

(10)

where H̃ is the rotation of matrix H at a degree of 180◦ . The
convolution value of element e is the summation of all the values in the derived matrix Ce .
The above procedure works for every element e. 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 MS of supported elements of M is determined
from (9), the set of unsupported element is derived consequently,
MU = M \ M S .
(11)
5

ρe

ρ
1

1

1

1

1

0

1

1

0

1

1

1

0

0

0

0

0

0

Rotate

1

1

1

1

1

0

1

1

0

0

0

0

0

0

0

1

1

1

H

0

Ce
0

0

0

0

0

1

1

0

Σ

2

2

H

Figure 5: The discrete convolution procedure in 2D for detecting supported elements using the designed self-supporting Kernel matrix H in Fig. 3.

θ

(a) Supporting elements

1

0

0

0

0

0

1

0

1

1

1

1

1

0

0

0

0

0

0

0

0

0

0

0

0

0

0

0

0

0

0

0

0

0

0

compliance ratio
C
,
(12)
Cre f
where Cre f , C is respectively the compliance of the structure
computed with or without considering self-supporting constraint.
The computational results for 2D examples are first summarized in Table 3; cases of 3D examples are explained later.

(b) Convolution Kernel matrix H

4.1. Cantilever beam example
The Cantilever beam, as shown in Fig. 7, is first tested. The
model on the left has a 150 × 60 rectangular domain, and has a
target volume fraction of 0.6. 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 0.5. 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 self-supporting 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
self-supporting 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 self-supporting. 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 self-supporting constraints have a very close
compliance, of a compliance ratio respectively of 100.11% and
101.73% 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 0, 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.

Figure 6: Procedure for building a self-supporting convolution Kernel matrix
for a general overhang angle θ.

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 E = 1 and υ = 0.3 for all examples. The penalization factor is set to a value of 3. The minimal thickness is
set to be 2, and the overhang angle is set to 45◦ 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 [21, 22, 12, 13], we measure
the ability of a self-supporting topology optimization approach
in maintaining the structure’s physical performance using the

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
6

Beam
Beam (hole)
Beam (r = 1.5)
Beam (r = 2)
Beam (r = 3)
Beam (concentrated)
Beam (distributed)
Beam (mixed)
Beam (vf=0.6)
Beam (vf=0.5)
Beam (vf=0.4)
Beam (vf=0.25)
Beam (Angle=30)
Beam (Angle=45)
Beam (Angle=60)
MBB
Square

#M

#MU

Cre f

C

C/Cre f

9000
7755
9000
9000
9000
14400
14400
14400
9000
9000
9000
9000
9000
9000
9000
38400
22500

24
25
24
18
12
21
12
60
24
97
37
55
12
97
425
1874
474

92.7
115.4
92.7
92.4
92.6
322.9
255.4
31971.1
92.7
105.8
127.4
196.8
105.8
105.8
105.8
185.7
1312.8

92.8
117.4
92.8
92.8
92.9
323.0
256.1
32514.3
92.8
106.8
128.3
202.9
113.6
106.8
141.3
191.3
1544.7

100.11%
101.73%
100.11%
101.43%
100.32%
100.31%
100.27%
101.70%
100.11%
100.95%
100.71%
103.10%
107.37%
100.95%
133.55%
103.02%
117.66%

Table 3: Summary of the numerical results for various tested 2D examples: #M, element number of the design domain; #MU , Cre f , number of unsupported
elements of structured computed from topology optimization without considering self-supporting constraints and its associated compliance; C, compliance of the
self-supporting structure computed using the proposed approach. The (Cantilever) beam examples are described in Figs. 7 and 14, the MBB example in Fig. 16, the
2D squares in Fig. 18, the 3D examples in Figs. 19,22 and 23

.

150

30

60

Print direction
(a) Rectangular domain

150

60

(a)
Supportneeded

30

30
20

(b)
Selfsupporting

(c)
Supportneeded

(d)
Selfsupporting

Figure 8: Computational results for the Cantilever beam example in Fig. 7. See
also Table 3 for more details.

75
Print direction

side project filter, as given in (c). After this, the self-supporting
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

(b) General domain
Figure 7: The Cantilever beam examples with a minimum thickness of 1.5.

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 self-supporting 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 black-white structure using the Heavi7

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
45◦ via using different convolution Kernel matrices H. We
demonstrate its performance still using the Cantilever beam
example in Fig. 22(a), for three different overhang angles:
30◦ , 45◦ , 60◦ . The associated convolution kernel matrices H for
angles of 30◦ , 60◦ 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 C/Cre f respectively
of 107.37%,100.95% and 133.55%. It is also noticed that the
larger overhang angle deteriorates the structures’ physical performance. Similar phenomenon were also observed in previous
studies [22, 12, 13].

1
0.9
0.8
0.7
0.6
0.5
0.4
0.3
0.2
0.1
0

(a) step 20

(b) step 44

(c) step 45 (22)

4.2. 2D Cantilever beam example with different types of forces

(d) step 54 (12)

(e) step 82 (6)

(f) step 101 (3)

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
240 × 60 square FE elements. The volume fraction is 0.6 and
the minimum thickness is 2. 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 100.31%, 100.27% and
101.70%.

(g) step 166 (0)

Figure 9: Key frame figures at different iteration steps, where the number in
brackets stands for the number of unsupported elements.

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 self-supporting structures of thicknesses 1.5, 2 and 3. 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
C/Cre f = 100.11%, 101.43%, 100.32%.

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. 16. Due to the symmetry of
the model’s structure and boundary conditions, only half of the
computational domain is used here which consists of 160 × 30
square FE mesh elements. The volume fraction is 0.5 and the
minimum thickness is 1.5.
The aim is to produce a self-supporting 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. 16(b). As a result, a self-supporting in both directions is

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 self-supporting structures become harder to
obtain for small value of volume fractions. As can be observed
from the results, the self-supporting constraints can still be satisfied although the number of elements decreases as the volume
fraction becomes smaller.
8

35

Unsupported elements

30

200

25

180

20

160

15

140

Compliance

Volume fraction

0.68
0.66
0.64

10

120

Iteration steps
5
0
20

(20,96.5) (45,94.8)

100

40

60

80

100

120

140

160

80
100

(a) Unsupported elements

Iteration 0.62
steps

Iteration steps

0.6

101

102

103

0

(b) Compliance

50

100

150

(c) Volume fraction

Figure 10: Variations of the structure’s number of unsupported elements, compliance and volume fraction during the optimization iteration steps.

(a) rmin = 1.5

(b) rmin = 2

(c) rmin = 3

(a) 0.6

(b) 0.5

(c) 0.4

(d) 0.25

Figure 12: Self-supporting structures computed using the proposed approach at
different volume fractions. See also Table 3 for more computational details.

Figure 11: Numerical results for the Cantilever beam example in Fig. 7(a) at
different thicknesses of 1.5, 2 and 3. See also Table 3 for more computational
details.

for such highly complex shape, as computing the supports or
its removal would be extremely troublesome if not impossible.

obtained in Fig. 17, as compared with the support-needed structure in Fig. 17(a), where the unsupported elements are plotted
in red. The compliance of the support-needed structure and the
self-supporting structure are respectively 185.7 and 191.3, at a
ratio of 103.02%.
Note that the 2D MBB problem of the same domain size was
also tested by Gaynor and Guest in [22], where the print direction was manually set from bottom to top. Such setting thus
does not require constraints of multiple constraints.

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.

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 support-needed structure and selfsupporting structure are respectively given in Fig. 18(b) and (c).
It can be seen from the results that the original support-needed
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
self-supporting structure of Fig. 18(b), despite the structure’s
high complexity. In addition, the self-supporting structure also
has a close compliance to that of the original support-needed
structure, respectively of 1544.7 and 1312.8 at a compliance ratio of 117.66%. A self-supporting structure is very necessary

4.5.1. 3D Wheel
The 3D wheel example in Fig. 19 consists of 100 × 100 × 100
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 0.25 and the
minimum thickness is 2. The print direction is chosen as top to
bottom.
The computed self-supporting structure using the proposed
approach is shown in Fig. 20(b), as compared with its counterpart of support-needed structure in Fig. 20(a). The corresponding structure slices at x = 25, 35, 45 of both the self-supporting
and support-needed structures are also shown and compared in
9

Table 4: Time comparison for 3D examples

Example

Wheel in Fig. 19

Cantilever(a) in Fig. 19

Cantilever(b) in Fig. 19

Desk in Fig. 23p

Size
Time

100 × 100 × 100
35m

100 × 50 × 50
11m

150 × 50 × 50
17m

120 × 120 × 80
66m

240
30
60

Distributed force
print direction

Concentrated force print direction

Mixed force print direction
Figure 14: The Cantilever beam example with distributed external forces.
(a) 30◦

(b) 45◦

(c) 60◦

1

1

0

0

0

1

1

1

0

1

0

0

1

1

1

0

0

0

1

0

0

0

0

0

0

0

0

0

0

0

0

0

0

0

0

0

0

0

0

0

0

0

0

0

0

0

0

0

0

0

(d) Kernel matrix H at 30◦

The example is of size 120 × 120 × 80, and consist of 1.152
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 0.3, the
minimum thickness is 2 and the print direction is chosen from
right to left. The final generated self-supporting structure using
the proposed approach is shown in Fig. 23, taking 66 minutes.
It is also interesting to note that the four legs of the desk are not
totally solid but take a porous bone-like structure to balance the
constraint of the object weight and the target compliance.

(e) Kernel matrix H at 60◦

4.6. Summary
As can be observed from these examples, using the proposed
self-supporting topology optimization approach, a completely
self-supporting structure is achieved that has a very similar geometric shape to that of the support-needed one without considering the self-supporting constraint. The edges or faces of
the original support-needed structure, are alighted toward the
print direction so that all elements can be successfully fabricated. In addition, the difference in compliance between the
self-supporting structure and the support-needed structure is
maintained within a very small or negligible range. All these
enhancements demonstrate the strength of the proposed approach in designing the self-supporting structures and simultaneously maintaining their optimized structure and physical performance.

Figure 13: Self-supporting structures obtained using the proposed approach
at different overhang angles, and their associated convolution Kernel matrices.
See also Table 3 for more computational details.

Fig. 20(f)(g)(h). As can be seen from the results, the originally
flat regions of the support-needed structure, which cannot be
fabricated without supports, have been optimized to meet the
self-supporting requirement. The resulting structure is totally
self-supporting for direct fabrication purpose and its compliance is 87.32, very close to that of the original support-needed
one of 86.16, of a ratio 1.01%
4.5.2. 3D Cantilever
Two different 3D Cantilever examples are tested here as illustrated in Fig. 22: one of size 100 × 50 × 50 exerted by point
loadings, and another one of size 150 × 50 × 50 exerted by edge
loadings; a same example to the latter was also studied in [13].
The target volume fraction is 0.3 and the minimum thickness
is 2. The print direction is chosen as top to bottom. The computed self-supporting structures are shown in Fig. 22, which is
directly to be fabricated without any additional support materials.

5. Conclusion
A novel self-supporting topology optimization approach is
provided in this paper applicable to additive manufacturing.
The usage of convolution operator and the associated numerical techniques enables the self-supporting 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
support-needed structure, proving its high effectiveness.

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.
10

240
80

(a) Half domain
Figure 16: The 2D MBB example: due to its symmetry only one half of the
model needs to be studied in the topology optimization process.

(a)
force

Point

(b)
force

Point

(c)
Distr.
forces

(d)
Distr.
forces

(a) Support-needed

(b) Self-supporting
(e) Point and distributed forces

Figure 17: The topology optimization result of the MBB model without (a)
or with (b) self-supporting structure. The self-supporting constraint is needed
simultaneously in two different directions: from left to right and from right to
left, so that the mirror symmetry of the original structure can be maintained.
See also Table 3 for more computational details.

a self-supporting 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 free-form 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
self-supporting 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.

(f) Point and distributed forces
Figure 15: The derived optimal structures of the Cantilever beam example in
Fig. 14 under different types of external forces. In all cases, structure similarity and compliance closeness is observed between the structures with or without considering the self-supporting constraints. Different print directions were
chosen for these different cases for ease of convergence and optimal physical
performance. See also Table 3 for more computational details.

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 self-supporting 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 self-supporting structure is
hard to obtain for the MBB example at a different direction.
Such phenomenon may not prevent the approach in generating

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.

11

Print direction

150

150

(a) Problem

(b) Support-needed

(c) Self-supporting

100
100

100

Print direction

Figure 18: The problem of a square under a distributed force, and its computational results. See also Table 3 for more computational details.

(a) Support-needed

(b) Self-supporting

Figure 19: The wheel example.
(c) Support-needed, x = (d) Support-needed, x = (e) Support-needed, x =
25
35
45

[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, Level-set methods
for structural topology optimization: a review, Structural and Multidisciplinary Optimization 48 (3) (2013) 437–472.
[7] X. Qian, Topology optimization in B-spline 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 AM-specific 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 print-ready designs, Structural and Multidisciplinary Optimization
(2016) 1–13.
[13] M. Langelaar, Topology optimization of 3D self-supporting 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 two-dimensional layered manufacturing, Computational Geometry 12 (3-4) (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, Cost-effective printing of 3D objects with skin-frame

(f) Self-supporting, x = (g) Self-supporting, x = (h) Self-supporting, x =
25
35
45
Figure 20: The self-supporting structure generated using the proposed approach
and its comparison with that of the support-needed one generated without the
consideration of self-supporting constraint. Both the 3D structure and their
slices at x = 25, 35, 45 are shown.

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, Self-supporting rhombic infill structures for additive manufacturing, Computer-Aided 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

12

50
150(100)
50

Print direction

F2
F1

Figure 21: A 3D Cantilever example of size 100 × 50 × 50 at a point loading
F1 , and of size 150 × 50 × 50 at an edge loading F2 .

[24]

[25]

[26]

[27]
[28]

(a) Self-supporting

(b) Slice in x = 25

(c) Self-supporting (d) Self-supporting (e) Slice in x = (f) Slice in x =
15
20
Figure 22: The generated self-supporting structures, together with its slice, for
two difference cases: (a),(b) for structure of size 100×50×50 at a point loading
F1 , and (c),(d),(e),(f) of size 150×50×50 at an edge loading F2 for the example
in Fig. 21.

100
100

80

Figure 23: A 3D desk example.

13

Print direction

[23]

constraints: Eliminating sacrificial support material in additive manufacturing through design, Structural & Multidisciplinary Optimization 54 (5)
(2016) 1157–1172.
O. Sigmund, Morphology-based black and white filters for topology optimization, Structural & Multidisciplinary Optimization 33 (4-5) (2007)
401–424.
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.
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.
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.
O. Sigmund, A 99 line topology optimization code written in matlab,
Structural & Multidisciplinary Optimization 21 (2) (2001) 120–127.
A. Dłaz, O. Sigmund, Checkerboard patterns in layout optimization,
Structural & Multidisciplinary Optimization 10 (1) (1995) 40–45.

(a) Self-supporting structure

(b) Slice in x = 25

(c) Slice in x = 35

Figure 24: The generated self-supporting structure and the associated slices at
x = 23 and x = 35, for the 3D desk example in Fig. 23.

14

