# Doubly 1-Bit Quantized Massive MIMO

Italo Atzeni,1 Antti Tölli,1 Duy H. N. Nguyen,2 and A. Lee Swindlehurst3
This work was supported by the Research Council of Finland (336449 Profi6,
346208 6G Flagship, 348396 HIGH-6G, and 357504 EETCAMD), by the European
Commission (101095759 Hexa-X-II), and by the National Science Foundation
(ECCS-2146436, CCF-2225575, and CCF-2225576). 1 Centre for Wireless
Communications, University of Oulu, Finland
2 Department of Electrical and Computer Engineering, San Diego State
University, USA
3 Department of Electrical Engineering and Computer Science, University of
California, Irvine, USA
Emails: {italo.atzeni<EMAIL_ADDRESS><EMAIL_ADDRESS><EMAIL_ADDRESS>
###### Abstract

Enabling communications in the (sub-)THz band will call for massive multiple-
input multiple-output (MIMO) arrays at either the transmit- or receive-side,
or at both. To scale down the complexity and power consumption when operating
across massive frequency and antenna dimensions, a sacrifice in the resolution
of the digital-to-analog/analog-to-digital converters (DACs/ADCs) will be
inevitable. In this paper, we analyze the extreme scenario where both the
transmit- and receive-side are equipped with fully digital massive MIMO arrays
and 1-bit DACs/ADCs, which leads to a system with minimum radio-frequency
complexity, cost, and power consumption. Building upon the Bussgang
decomposition, we derive a tractable approximation of the mean squared error
(MSE) between the transmitted data symbols and their soft estimates. Numerical
results show that, despite its simplicity, a doubly 1-bit quantized massive
MIMO system with very large antenna arrays can deliver an impressive
performance in terms of MSE and symbol error rate.

###### Index Terms:

1-bit ADCs, 1-bit DACs, doubly massive MIMO, (sub-)THz communications.

##  I Introduction

Beyond-5G wireless systems are envisioned to exploit the large amount of
bandwidth in the sub-THz and THz frequency ranges (0.1–0.3 THz and 0.3–3 THz,
respectively) [1]. To overcome the severe pathloss and penetration loss,
massive multiple-input multiple-output (MIMO) arrays will be needed at either
the transmit- or receive-side, or at both [2]. Fully digital massive MIMO
architectures provide highly flexible beamforming and large-scale spatial
multiplexing while avoiding the complex beam-management schemes of their
hybrid analog-digital counterparts. However, the power consumption of the
digital-to-analog/analog-to-digital converters (DACs/ADCs), which scales
linearly with the bandwidth and exponentially with the number of resolution
bits [3], poses significant practical challenges. In this context, low-
resolution DACs/ADCs have been considered as a promising enabler of truly
massive, fully digital antenna arrays [4, 5, 6]. Simple 1-bit DACs/ADCs can
also alleviate the overall complexity and power consumption of the radio-
frequency (RF) chains. For instance, 1-bit DACs at the transmitter allow the
use of low-cost power amplifiers that are not constrained to operate with
backoff [7], whereas 1-bit ADCs at the receiver relax the requirements on the
automatic gain control. Hence, the performance loss arising from the 1-bit
quantization can be compensated for by adding very low-cost antennas and RF
chains.

Massive MIMO systems with low-resolution data converters have been generally
studied assuming coarse quantization at the base station (either in the DACs
[4, 8] or in the ADCs [5, 9]) and full-resolution user equipment. In this
paper, we analyze the extreme scenario where both the transmit- and receive-
side are equipped with fully digital massive MIMO arrays [2] and 1-bit
DACs/ADCs, which is referred to in the following as doubly 1-bit quantized
massive MIMO. Indeed, combining 1-bit DACs and ADCs leads to a fully digital
system with minimum RF complexity, cost, and power consumption. Furthermore,
the short wavelengths characterizing the (sub-)THz band and the use of 1-bit
DACs/ADCs can pave the way for compact and energy-efficient massive MIMO
arrays that can also be implemented at the user equipment (e.g., on cars and
unmanned aerial vehicles). Only a handful of works in the literature have
considered MIMO systems with both 1-bit DACs and 1-bit ADCs. For instance,
[10] studied the 1-bit multiple-input single-output capacity with perfect
channel state information (CSI) at both the transmitter and receiver.
Moreover, [11] described the 1-bit MIMO capacity with imperfect CSI assuming
that the number of transmit or receive antennas tends to infinity. Lastly,
[12] showed that a proper combination of transmit beamforming and equiprobable
signaling allows the system to operate close to the 1-bit MIMO capacity.

$\mathbf{s}\longrightarrow\underbrace{\mathbf{x}=\mathbf{W}\mathbf{s}}_{\textnormal{precoding}}\longrightarrow\underbrace{\mathbf{t}=Q_{\textnormal{\tiny{TX}}}(\mathbf{x})}_{\textnormal{1-bit
DACs}}\longrightarrow\sqrt{\rho}\mathbf{t}$$\mathbf{y}=\sqrt{\rho}\mathbf{H}\mathbf{t}+\mathbf{z}\longrightarrow\underbrace{\mathbf{r}=Q_{\textnormal{\tiny{RX}}}(\mathbf{y})}_{\textnormal{1-bit
ADCs}}\longrightarrow\underbrace{\hat{\mathbf{s}}=\mathbf{V}^{\mathrm{H}}\mathbf{r}}_{\textnormal{combining}}$TransmitterReceiver$\mathbf{H}$
Figure 1: Diagram of the considered doubly 1-bit quantized massive MIMO
system.

The goal of this paper is to provide a tractable analytical framework that
lends itself to the performance analysis and optimization of doubly 1-bit
quantized massive MIMO systems. Considering a point-to-point system, we build
upon the Bussgang decomposition [13] to unfold the relation between the
transmitted data symbols (at the input of the transmitter’s 1-bit DACs) and
the soft-estimated symbols acquired via linear combining of the quantized
received signal (at the output of the receiver’s 1-bit ADCs). Assuming perfect
channel estimation, we derive a tractable approximation of the mean squared
error (MSE) between the transmitted data symbols and their soft estimates as
well as the combining strategy that minimizes it. This approximation is
accurate in the regime of a large number of transmit antennas, which makes the
signal at the input of the receiver’s 1-bit ADCs approximately Gaussian.
Numerical results show that, despite its simplicity, a doubly 1-bit quantized
massive MIMO system with very large antenna arrays can deliver an impressive
performance in terms of MSE and symbol error rate (SER), which is not far from
that of a massive MIMO system with full-resolution DACs and 1-bit ADCs.

##  II System Model

We consider a point-to-point doubly 1-bit quantized massive MIMO system, where
a transmitter equipped with $N$ antennas and 1-bit DACs transmits $K$ data
streams to a receiver with $M$ antennas and 1-bit ADCs, with $K\leq\min(N,M)$.
Such a point-to-point system may represent, e.g., a wireless backhaul
scenario; however, the following discussion can be easily extended to a multi-
user uplink or downlink scenario. To model the 1-bit DACs and ADCs, we
introduce the 1-bit quantization function
$Q_{\textnormal{\tiny{A}}}(\cdot):\mbox{$\mathbb{C}$}^{B\times
1}\to\mathcal{Q}_{\textnormal{\tiny{A}}}$, with
$\mathcal{Q}_{\textnormal{\tiny{A}}}\triangleq\sqrt{\frac{\eta_{\textnormal{\tiny{A}}}}{2}}\\{\pm
1\pm j\\}^{B\times 1}$ and

$\displaystyle
Q_{\textnormal{\tiny{A}}}(\mathbf{b})\triangleq\sqrt{\frac{\eta_{\textnormal{\tiny{A}}}}{2}}\Big{(}\mathrm{sgn}\big{(}\mathrm{Re}[\mathbf{b}]\big{)}+j\,\mathrm{sgn}\big{(}\mathrm{Im}[\mathbf{b}]\big{)}\Big{)}.$
(1)

Note that the output of the 1-bit quantization function is a vector of scaled
quadrature phase-shift keying (QPSK) symbols. In the following, we use the
subscripts $\textnormal{A}=\textnormal{TX}$ and
$\textnormal{A}=\textnormal{RX}$ to indicate “transmitter” and “receiver”,
respectively.

Let $\mathbf{H}\in\mbox{$\mathbb{C}$}^{M\times N}$ denote the channel matrix
between the transmitter and the receiver. In this paper, we assume that
$\mathbf{H}$ is perfectly known at both the transmitter and receiver, and we
leave the analysis with imperfect CSI for future work. The transmitter aims at
conveying the data symbol vector $\mathbf{s}\in\mbox{$\mathbb{C}$}^{K\times
1}$ to the receiver. As in [4, 8], we consider a quantized linear precoding
strategy, whereby the precoding matrix is designed based on $\mathbf{H}$ and
independently of $\mathbf{s}$ and the subsequent quantization step.111An
alternative approach, which goes beyond the scope of this paper, is symbol-
level precoding, whereby the analog signal at the output of the 1-bit DACs is
designed based on $\mathbf{H}$ and $\mathbf{s}$ [7]. Symbol-level precoding
outperforms its quantized linear counterpart at the cost of higher complexity.
In this setting, $\mathbf{s}$ is first precoded as

$\displaystyle\mathbf{x}\triangleq\mathbf{W}\mathbf{s}\in\mbox{$\mathbb{C}$}^{N\times
1}$ (2)

where $\mathbf{W}\in\mbox{$\mathbb{C}$}^{N\times K}$ is the precoding matrix,
and then quantized at the 1-bit DACs as

$\displaystyle\mathbf{t}\triangleq Q_{\textnormal{\tiny{TX}}}(\mathbf{x}).$
(3)

Here, the scaling factor $\eta_{\textnormal{\tiny{TX}}}$ of the 1-bit
quantization function in (1) is fixed as
$\eta_{\textnormal{\tiny{TX}}}=\frac{1}{N}$ to satisfy the power constraint
$\|\mathbf{t}\|^{2}=1$.

Subsequently, the analog signal $\mathbf{t}$ is transmitted over the channel
with transmit power $\rho$ and the signal arriving at the receiver is given by

$\displaystyle\mathbf{y}\triangleq\sqrt{\rho}\mathbf{H}\mathbf{t}+\mathbf{z}\in\mbox{$\mathbb{C}$}^{M\times
1}$ (4)

where $\mathbf{z}\sim\mathcal{C}\mathcal{N}(\mathbf{0},\mathbf{I}_{M})$ is a
vector of additive white Gaussian noise (AWGN). Since the AWGN has unit
variance, $\rho$ can be interpreted as the transmit signal-to-noise ratio
(SNR). Then, $\mathbf{y}$ is quantized at the 1-bit ADCs as

$\displaystyle\mathbf{r}\triangleq
Q_{\textnormal{\tiny{RX}}}(\mathbf{y})=Q_{\textnormal{\tiny{RX}}}\big{(}\sqrt{\rho}\mathbf{H}Q_{\textnormal{\tiny{TX}}}(\mathbf{W}\mathbf{s})+\mathbf{z}\big{)}.$
(5)

Here, the scaling factor $\eta_{\textnormal{\tiny{RX}}}$ of the 1-bit
quantization function in (1) is fixed as
$\eta_{\textnormal{\tiny{RX}}}=\rho+1$ so that the element-wise variance of
the output coincides with that of the input when the channel elements have
unit variance and $N\to\infty$ (see (15) in the following). Note that the
doubly 1-bit quantized signal in (5) is what is observed at the receiver.
Finally, the receiver acquires a soft estimate of $\mathbf{s}$ via linear
combining of the digital signal $\mathbf{r}$ as

$\displaystyle\hat{\mathbf{s}}\triangleq\mathbf{V}^{\mathrm{H}}\mathbf{r}\in\mbox{$\mathbb{C}$}^{K\times
1}$ (6)

where $\mathbf{V}\in\mbox{$\mathbb{C}$}^{M\times K}$ is the combining matrix.

In this paper, we consider the MSE between $\mathbf{s}$ and its soft estimate
$\hat{\mathbf{s}}$ in (6) as a performance metric, which is given by

$\displaystyle\varepsilon\triangleq\frac{1}{K}\mathbb{E}_{\mathbf{s},\mathbf{z}}[\|\mathbf{s}-\hat{\mathbf{s}}\|^{2}].$
(7)

We observe that deriving the above MSE requires obtaining a tractable
expression for $\mathbf{r}$ in (5) as a function of $\mathbf{s}$.

##  III Linearization Via the Bussgang Decomposition

In this section, we consider Gaussian data symbols, i.e.,
$\mathbf{s}\sim\mathcal{C}\mathcal{N}(\mathbf{0},\mathbf{I}_{K})$. Then, we
express the doubly 1-bit quantized signal observed at the receiver as a linear
function using the Bussgang decomposition [13], which allows one to write the
output of a nonlinear system as a scaled version of the input plus an
uncorrelated distortion. Finally, we derive a tractable approximation of the
MSE in (7).

###  III-A Linearization at the Transmitter

$\displaystyle\mathbf{C}_{\mathbf{t}}$
$\displaystyle=\frac{2}{\pi}\eta_{\textnormal{\tiny{TX}}}\Big{(}\arcsin\big{(}\mathrm{Diag}(\mathbf{C}_{\mathbf{x}})^{-\frac{1}{2}}\mathrm{Re}[\mathbf{C}_{\mathbf{x}}]\mathrm{Diag}(\mathbf{C}_{\mathbf{x}})^{-\frac{1}{2}}\big{)}+j\,\arcsin\big{(}\mathrm{Diag}(\mathbf{C}_{\mathbf{x}})^{-\frac{1}{2}}\mathrm{Im}[\mathbf{C}_{\mathbf{x}}]\mathrm{Diag}(\mathbf{C}_{\mathbf{x}})^{-\frac{1}{2}}\big{)}\Big{)},$
(14) $\displaystyle\tilde{\mathbf{C}}_{\mathbf{r}}$
$\displaystyle\triangleq\frac{2}{\pi}\eta_{\textnormal{\tiny{RX}}}\Big{(}\arcsin\big{(}\mathrm{Diag}(\mathbf{C}_{\mathbf{y}})^{-\frac{1}{2}}\mathrm{Re}[\mathbf{C}_{\mathbf{y}}]\mathrm{Diag}(\mathbf{C}_{\mathbf{y}})^{-\frac{1}{2}}\big{)}+j\,\arcsin\big{(}\mathrm{Diag}(\mathbf{C}_{\mathbf{y}})^{-\frac{1}{2}}\mathrm{Im}[\mathbf{C}_{\mathbf{y}}]\mathrm{Diag}(\mathbf{C}_{\mathbf{y}})^{-\frac{1}{2}}\big{)}\Big{)}$
(21)

Let us define

$\displaystyle\mathbf{C}_{\mathbf{x}}$
$\displaystyle\triangleq\mathbb{E}_{\mathbf{s}}[\mathbf{x}\mathbf{x}^{\mathrm{H}}]=\mathbf{W}\mathbf{W}^{\mathrm{H}}\in\mbox{$\mathbb{C}$}^{N\times
N},$ (8) $\displaystyle\mathbf{C}_{\mathbf{t}}$
$\displaystyle\triangleq\mathbb{E}_{\mathbf{s}}[\mathbf{t}\mathbf{t}^{\mathrm{H}}]\in\mbox{$\mathbb{C}$}^{N\times
N}.$ (9)

Then, we use the Bussgang decomposition to linearize $\mathbf{t}$ in (3) with
respect to $\mathbf{x}$ in (2) (and, thus, with respect to $\mathbf{s}$) as

$\displaystyle\mathbf{t}=\mathbf{G}_{\textnormal{\tiny{TX}}}\mathbf{x}+\mathbf{d}_{\textnormal{\tiny{TX}}}$
(10)

where $\mathbf{d}_{\textnormal{\tiny{TX}}}\in\mbox{$\mathbb{C}$}^{N\times 1}$
is a zero-mean, non-Gaussian distortion vector that is uncorrelated with
$\mathbf{x}$ (and, obviously, with $\mathbf{s}$) and

$\displaystyle\mathbf{G}_{\textnormal{\tiny{TX}}}\triangleq\mathbb{E}_{\mathbf{s}}[\mathbf{t}\mathbf{x}^{\mathrm{H}}]\mathbf{C}_{\mathbf{x}}^{-1}\in\mbox{$\mathbb{C}$}^{N\times
N}$ (11)

is the Bussgang gain matrix. Indeed,
$\mathbf{G}_{\textnormal{\tiny{TX}}}\mathbf{x}$ and
$\mathbf{d}_{\textnormal{\tiny{TX}}}$ are the minimum MSE estimate of
$\mathbf{t}$ given $\mathbf{x}$ and the corresponding estimation error,
respectively. Moreover, according to (10), $\mathbf{C}_{\mathbf{t}}$ in (9)
can be obtained as

$\displaystyle\mathbf{C}_{\mathbf{t}}=\mathbf{G}_{\textnormal{\tiny{TX}}}\mathbf{C}_{\mathbf{x}}\mathbf{G}_{\textnormal{\tiny{TX}}}+\mathbf{C}_{\mathbf{d}_{\textnormal{\tiny{TX}}}}$
(12)

with
$\mathbf{C}_{\mathbf{d}_{\textnormal{\tiny{TX}}}}\triangleq\mathbb{E}[\mathbf{d}_{\textnormal{\tiny{TX}}}\mathbf{d}_{\textnormal{\tiny{TX}}}^{\mathrm{H}}]\in\mbox{$\mathbb{C}$}^{N\times
N}$.

Since $\mathbf{s}$ is Gaussian, we have that
$\mathbf{G}_{\textnormal{\tiny{TX}}}$ is a diagonal matrix that can be
computed in closed form following well-known steps (see, e.g., [5]) as

$\displaystyle\mathbf{G}_{\textnormal{\tiny{TX}}}=\sqrt{\frac{2}{\pi}\eta_{\textnormal{\tiny{TX}}}}\mathrm{Diag}(\mathbf{C}_{\mathbf{x}})^{-\frac{1}{2}}.$
(13)

Note that, when $K<N$ (i.e., when the number of data streams is strictly
smaller than the number of transmit antennas, which is generally the case in
downlink massive MIMO systems), $\mathbf{C}_{\mathbf{x}}$ in (8) is rank-
deficient and $\mathbf{G}_{\textnormal{\tiny{TX}}}$ cannot be computed as in
(11);222When $\mathbf{C}_{\mathbf{x}}$ is rank-deficient, its inverse in (11)
can be replaced by its pseudoinverse, as described in [14]. In this case, the
resulting $\mathbf{G}_{\textnormal{\tiny{TX}}}$ may not be diagonal, although
$\mathbf{G}_{\textnormal{\tiny{TX}}}\mathbf{x}$ will be the same whether
$\mathbf{G}_{\textnormal{\tiny{TX}}}$ is obtained via (11) or (13).
nonetheless, it can still be obtained via (13), which does not involve the
inversion of $\mathbf{C}_{\mathbf{x}}$. Furthermore, $\mathbf{C}_{\mathbf{t}}$
in (9) can be computed in closed form as in (14) at the top of the page [5].

###  III-B Linearization at the Receiver

Let us define

$\displaystyle\mathbf{C}_{\mathbf{y}}$
$\displaystyle\triangleq\mathbb{E}_{\mathbf{s},\mathbf{z}}[\mathbf{y}\mathbf{y}^{\mathrm{H}}]=\rho\mathbf{H}\mathbf{C}_{\mathbf{t}}\mathbf{H}^{\mathrm{H}}+\mathbf{I}_{M}\in\mbox{$\mathbb{C}$}^{M\times
M},$ (15) $\displaystyle\mathbf{C}_{\mathbf{r}}$
$\displaystyle\triangleq\mathbb{E}_{\mathbf{s},\mathbf{z}}[\mathbf{r}\mathbf{r}^{\mathrm{H}}]\in\mbox{$\mathbb{C}$}^{M\times
M}$ (16)

where we point out that the above expectations are taken also over the AWGN.
Then, we use the Bussgang decomposition to linearize $\mathbf{r}$ in (5) with
respect to $\mathbf{y}$ in (4) as

$\displaystyle\mathbf{r}=\mathbf{G}_{\textnormal{\tiny{RX}}}\mathbf{y}+\mathbf{d}_{\textnormal{\tiny{RX}}}=\mathbf{G}_{\textnormal{\tiny{RX}}}\big{(}\sqrt{\rho}\mathbf{H}(\mathbf{G}_{\textnormal{\tiny{TX}}}\mathbf{W}\mathbf{s}+\mathbf{d}_{\textnormal{\tiny{TX}}})+\mathbf{z}\big{)}+\mathbf{d}_{\textnormal{\tiny{RX}}}$
(17)

where $\mathbf{d}_{\textnormal{\tiny{RX}}}\in\mbox{$\mathbb{C}$}^{M\times 1}$
is a zero-mean, non-Gaussian distortion vector that is uncorrelated with
$\mathbf{y}$ and

$\displaystyle\mathbf{G}_{\textnormal{\tiny{RX}}}\triangleq\mathbb{E}_{\mathbf{y}}[\mathbf{r}\mathbf{y}^{\mathrm{H}}]\mathbf{C}_{\mathbf{y}}^{-1}\in\mbox{$\mathbb{C}$}^{M\times
M}$ (18)

is the Bussgang gain matrix. Moreover, according to (17),
$\mathbf{C}_{\mathbf{r}}$ in (16) can be obtained as

$\displaystyle\mathbf{C}_{\mathbf{r}}\triangleq\mathbf{G}_{\textnormal{\tiny{RX}}}\mathbf{C}_{\mathbf{y}}\mathbf{G}_{\textnormal{\tiny{RX}}}+\mathbf{C}_{\mathbf{d}_{\textnormal{\tiny{RX}}}}$
(19)

with
$\mathbf{C}_{\mathbf{d}_{\textnormal{\tiny{RX}}}}\triangleq\mathbb{E}[\mathbf{d}_{\textnormal{\tiny{RX}}}\mathbf{d}_{\textnormal{\tiny{RX}}}^{\mathrm{H}}]\in\mbox{$\mathbb{C}$}^{M\times
M}$.

At this stage, we have achieved our goal to express the doubly 1-bit quantized
signal observed at the receiver as a linear function of $\mathbf{s}$. However,
since $\mathbf{y}$ is not Gaussian due to the 1-bit DACs at the
transmitter,333This can be also observed from (17), which includes the non-
Gaussian distortion vector $\mathbf{d}_{\textnormal{\tiny{TX}}}$. we have that
$\mathbf{G}_{\textnormal{\tiny{RX}}}$ is generally not diagonal. Consequently,
$\mathbf{C}_{\mathbf{r}}$ and $\mathbf{G}_{\textnormal{\tiny{RX}}}$ are not
available in closed form, which makes the expression in (17) not tractable.
Despite that, we make the following observation: each element of
$\mathbf{u}\triangleq\mathbf{H}\mathbf{t}\in\mbox{$\mathbb{C}$}^{M\times 1}$
is a weighted sum of $N$ scaled QPSK symbols. Therefore, when the number of
transmit antennas $N$ is large and the channel elements are weakly correlated,
we have that $\mathbf{u}$ is approximately Gaussian, i.e.,
$\mathbf{u}\mathrel{\dot{\sim}}\mathcal{C}\mathcal{N}(\mathbf{0},\mathbf{H}\mathbf{C}_{\mathbf{t}}\mathbf{H}^{\mathrm{H}})$.
This approximation becomes asymptotically exact as $N\to\infty$; nonetheless,
with i.i.d. channel elements, we observe through numerical simulations that it
is already quite accurate for $N\geq 8$. In this context, we have that
$\mathbf{y}$ is approximately Gaussian, i.e.,
$\mathbf{y}\mathrel{\dot{\sim}}\mathcal{C}\mathcal{N}(\mathbf{0},\mathbf{C}_{\mathbf{y}})$.
Exploiting this argument, $\mathbf{G}_{\textnormal{\tiny{RX}}}$ in (18) can be
approximated as
$\mathbf{G}_{\textnormal{\tiny{RX}}}\approx\tilde{\mathbf{G}}_{\textnormal{\tiny{RX}}}$,
where

$\displaystyle\tilde{\mathbf{G}}_{\textnormal{\tiny{RX}}}\triangleq\sqrt{\frac{2}{\pi}\eta_{\textnormal{\tiny{RX}}}}\mathrm{Diag}(\mathbf{C}_{\mathbf{y}})^{-\frac{1}{2}}$
(20)

is obtained following similar steps as for (13). Furthermore,
$\mathbf{C}_{\mathbf{r}}$ in (16) can be approximated as
$\mathbf{C}_{\mathbf{r}}\approx\tilde{\mathbf{C}}_{\mathbf{r}}$, where
$\tilde{\mathbf{C}}_{\mathbf{r}}$ is given in (21) at the top of the page and
results from following similar steps as for (14).

###  III-C Tractable Approximation of the MSE

Building upon the linearizations at the transmitter and at the receiver
described in Sections III-A and III-B, respectively, we approximate the MSE in
(7) as follows.

######  Proposition 1.

The MSE in (7) can be approximated as $\varepsilon\approx\tilde{\varepsilon}$,
with

$\displaystyle\tilde{\varepsilon}\triangleq
1+\frac{1}{K}\mathrm{tr}(\mathbf{V}^{\mathrm{H}}\tilde{\mathbf{C}}_{\mathbf{r}}\mathbf{V})-\frac{2}{K}\sqrt{\rho}\mathrm{tr}\big{(}\mathrm{Re}[\mathbf{V}^{\mathrm{H}}\tilde{\mathbf{G}}_{\textnormal{\tiny{RX}}}\mathbf{H}\mathbf{G}_{\textnormal{\tiny{TX}}}\mathbf{W}]\big{)}.$
(22)

###### Proof:

See Appendix A ∎

Our numerical simulations in Section IV show that the approximate MSE
$\tilde{\varepsilon}$ in (22) behaves as an upper bound on the true MSE
$\varepsilon$. Moreover, for a fixed precoding matrix $\mathbf{W}$,444Note
that (22) depends on the specific choice of the precoding matrix $\mathbf{W}$
and also on $\mathbf{C}_{\mathbf{x}}$ in (8) through
$\mathbf{G}_{\textnormal{\tiny{TX}}}$ in (13),
$\tilde{\mathbf{G}}_{\textnormal{\tiny{RX}}}$ in (20), and
$\tilde{\mathbf{C}}_{\mathbf{r}}$ in (21). we have that $\tilde{\varepsilon}$
is a convex quadratic function of the combining matrix $\mathbf{V}$. Hence, we
can minimize the approximate MSE with respect to $\mathbf{V}$, which yields

$\displaystyle\mathbf{V}^{\star}\triangleq\underset{\mathbf{V}}{\operatornamewithlimits{argmin}}\,\tilde{\varepsilon}=\sqrt{\rho}\tilde{\mathbf{C}}_{\mathbf{r}}^{-1}\tilde{\mathbf{G}}_{\textnormal{\tiny{RX}}}\mathbf{H}\mathbf{G}_{\textnormal{\tiny{TX}}}\mathbf{W}.$
(23)

##  IV Numerical Results

In this section, we evaluate the performance of the considered doubly 1-bit
quantized massive MIMO system in terms of MSE between the transmitted data
symbols and their soft estimates. We consider far-field propagation and
generate the channel matrix $\mathbf{H}$ using the discrete physical channel
model [15]. The transmitter and receiver are equipped with square uniform
planar arrays of $\sqrt{N}\times\sqrt{N}$ (resp. $\sqrt{M}\times\sqrt{M}$)
half-wavelength spaced antennas and are placed with their broadsides facing
each other. Between them lies a cluster of $10^{2}$ scatterers, which give
rise to as many independent propagation paths. For both the transmitter and
receiver, the scatterers are confined within an angle spread of
$\frac{\pi}{6}$ around the broadside direction in both the azimuth and
elevation. The channels are normalized such that their elements have unit
variance and the pathloss is incorporated into the transmit SNR $\rho$. The
following numerical results are obtained by averaging over $10^{3}$
independent channel realizations. For each realization of $\mathbf{H}$, the
precoding matrix $\mathbf{W}$ comprises the $K$ principal right eigenvectors
of $\mathbf{H}$, whereas the combining matrix $\mathbf{V}$ is computed as in
(23), which minimizes the approximate MSE in (22). The values of the
simulation parameters $N$, $M$, $K$, and $\rho$ are reported above each
figure.

40057678410241296160000.050.10.150.2$N=M$MSE$K=16$, $\rho=10$
dB$\tilde{\varepsilon}$ (1-bit DACs/ADCs)$\varepsilon$ (1-bit
DACs/ADCs)$\varepsilon$ (1-bit ADCs) Figure 2: MSE for 1-bit DACs/ADCs and
full-resolution DACs/1-bit ADCs versus number of transmit/receive antennas.
The true MSE $\varepsilon$ is obtained via Monte Carlo simulations.
4005767841024129616000.050.10.150.2$N$$\tilde{\varepsilon}$$K=16$, $\rho=10$
dB$M=400$$M=1024$$M=1600$ Figure 3: Approximate MSE for 1-bit DACs/ADCs versus
number of transmit antennas for different numbers of receive antennas. The
dotted curves are obtained by switching $N$ and $M$.

Figure 2 plots the MSE versus the number of transmit/receive antennas, with
$N=M$. The approximate MSE $\tilde{\varepsilon}$ in (22) behaves as a tight
upper bound on the true MSE $\varepsilon$ in (7), where the latter is obtained
via Monte Carlo simulations with $10^{3}$ independent realizations of the data
symbol vector $\mathbf{s}$ for each realization of $\mathbf{H}$. Furthermore,
we observe that truly massive antenna arrays at both the transmitter and
receiver are necessary to achieve impressive values of the MSE, e.g., about
$5\times 10^{-2}$ for $N=M=1600$. Nonetheless, the performance of the
considered doubly 1-bit quantized massive MIMO system is not far from that of
a massive MIMO system with full-resolution DACs and 1-bit ADCs. Specifically,
replacing the 1-bit DACs with full-resolution ones reduces the MSE not even by
a factor of two at the cost of much higher RF complexity and power consumption
at the transmitter.

248163264$0$$0.1$$0.2$$0.3$$0.4$$0.5$$K$$\tilde{\varepsilon}$$\rho=10$
dB$N=M=400$$N=M=1024$$N=M=1600$ Figure 4: Approximate MSE for 1-bit DACs/ADCs
versus number of data streams for different numbers of transmit/receive
antennas.

$-1$$0$$1$$-1$$0$$1$IQ$N=M=400$, $K=8$, $\rho=10$ dB$\hat{\mathbf{s}}$ (1-bit
DACs/ADCs): $\textrm{SER}\\!=\\!4.9\\!\times\\!10^{-2}$$\hat{\mathbf{s}}$
(1-bit ADCs): $\textrm{SER}\\!=\\!5.8\\!\times\\!10^{-3}$

$-1$$0$$1$$-1$$0$$1$IQ$N=M=1024$, $K=8$, $\rho=10$ dB$\hat{\mathbf{s}}$ (1-bit
DACs/ADCs): $\textrm{SER}\\!=\\!2.1\\!\times\\!10^{-3}$$\hat{\mathbf{s}}$
(1-bit ADCs): $\textrm{SER}\\!=\\!4.7\\!\times\\!10^{-6}$

$-1$$0$$1$$-1$$0$$1$IQ$N=M=1600$, $K=8$, $\rho=10$ dB$\hat{\mathbf{s}}$ (1-bit
DACs/ADCs): $\textrm{SER}\\!=\\!1.4\\!\times\\!10^{-4}$$\hat{\mathbf{s}}$
(1-bit ADCs): $\textrm{SER}\\!=\\!1.1\\!\times\\!10^{-8}$

Figure 5: Soft-estimated symbols for 1-bit DACs/ADCs and full-resolution
DACs/1-bit ADCs with 16-PSK data symbols for different numbers of
transmit/receive antennas. The magenta solid marks represent the transmitted
data symbols.

Figure 3 illustrates the approximate MSE versus the number of transmit
antennas for different numbers of receive antennas (solid lines) and versus
the number of receive antennas for different numbers of transmit antennas
(dotted lines). We observe that increasing the number of either transmit or
receive antennas produces roughly the same effect. However, since the
combining matrix is optimized for a given channel and precoding matrix, the
second option provides slightly better results in this case. Figure 4 depicts
the approximate MSE versus the number of data streams for different numbers of
transmit/receive antennas, with $N=M$. For each configuration, there is an
optimal number of data streams: on the one hand, judiciously increasing the
number of data streams generates a useful scrambling of the 1-bit quantized
signals at the $M$ receive antennas [9]; on the other hand, the inter-stream
interference becomes dominant for large values of $K$.

##  V Data Detection

In this section, we briefly evaluate the data detection performance with non-
Gaussian data symbols in terms of SER. In this respect, we point out that
$\mathbf{x}$ in (2) may be approximately Gaussian even with non-Gaussian data
symbols when $K$ is large. Figure 5 plots the soft-estimated symbols with
16-PSK (phase-shift keying) data symbols for different numbers of
transmit/receive antennas, with $N=M$. As the number of antennas increases,
the dispersion of the soft-estimated symbols around the transmitted data
symbols reduces noticeably, which translates into an improved SER performance.
For $N=M=1600$, a remarkable SER in the order of $10^{-4}$ is obtained.
Nonetheless, an acceptable SER (considering the absence of coding) in the
order of $10^{-2}$ is achieved already for $N=M=400$.

##  VI Conclusions

Among fully digital systems, doubly 1-bit quantized massive MIMO systems are
endowed with minimum RF complexity, cost, and power consumption. In this
setting, we derived a tractable approximation of the MSE between the
transmitted data symbols and their soft estimates as well as the combining
strategy that minimizes it. We showed that, despite its simplicity, a doubly
1-bit quantized massive MIMO system with very large antenna arrays can deliver
an impressive performance in terms of MSE and SER, which is not far from that
of a massive MIMO system with full-resolution DACs and 1-bit ADCs. Future work
will analyze the impact of imperfect CSI along with the overall energy
efficiency.

##  Appendix A Proof of Proposition 1

We begin by writing (7) as

$\displaystyle\varepsilon$
$\displaystyle=\frac{1}{K}\mathbb{E}_{\mathbf{s},\mathbf{z}}\big{[}\|\mathbf{s}-\mathbf{V}^{\mathrm{H}}\mathbf{r}\|^{2}\big{]}$
(24)
$\displaystyle=\frac{1}{K}\mathbb{E}_{\mathbf{s},\mathbf{z}}\big{[}\mathrm{tr}(\mathbf{s}\mathbf{s}^{\mathrm{H}})+\mathrm{tr}(\mathbf{V}^{\mathrm{H}}\mathbf{r}\mathbf{r}^{\mathrm{H}}\mathbf{V})-2\mathrm{tr}\big{(}\mathrm{Re}[\mathbf{V}^{\mathrm{H}}\mathbf{r}\mathbf{s}^{\mathrm{H}}]\big{)}\big{]}$
(25)
$\displaystyle=1+\frac{1}{K}\mathrm{tr}(\mathbf{V}^{\mathrm{H}}\mathbf{C}_{\mathbf{r}}\mathbf{V})-\frac{2}{K}\mathrm{tr}\big{(}\mathrm{Re}\big{[}\mathbf{V}^{\mathrm{H}}\mathbb{E}[\mathbf{r}\mathbf{s}^{\mathrm{H}}]\big{]}\big{)}.$
(26)

Then, by plugging (17) into the remaining expectation term in (26), we obtain

$\displaystyle\mathbb{E}[\mathbf{r}\mathbf{s}^{\mathrm{H}}]$
$\displaystyle=\mathbb{E}_{\mathbf{s},\mathbf{z}}\Big{[}\big{(}\mathbf{G}_{\textnormal{\tiny{RX}}}\big{(}\sqrt{\rho}\mathbf{H}(\mathbf{G}_{\textnormal{\tiny{TX}}}\mathbf{W}\mathbf{s}+\mathbf{d}_{\textnormal{\tiny{TX}}})+\mathbf{z}\big{)}+\mathbf{d}_{\textnormal{\tiny{RX}}}\big{)}\mathbf{s}^{\mathrm{H}}\Big{]}$
(27)
$\displaystyle=\sqrt{\rho}\mathbf{G}_{\textnormal{\tiny{RX}}}\mathbf{H}\mathbf{G}_{\textnormal{\tiny{TX}}}\mathbf{W}\mathbb{E}_{\mathbf{s}}[\mathbf{s}\mathbf{s}^{\mathrm{H}}]$
(28)

where (28) follows from the fact that $\mathbf{d}_{\textnormal{\tiny{TX}}}$,
$\mathbf{z}$, and $\mathbf{d}_{\textnormal{\tiny{RX}}}$ are all uncorrelated
with $\mathbf{s}$. Hence, plugging (28) into (26) yields

$\displaystyle\varepsilon=1+\frac{1}{K}\mathrm{tr}(\mathbf{V}^{\mathrm{H}}\mathbf{C}_{\mathbf{r}}\mathbf{V})-\frac{2}{K}\sqrt{\rho}\mathrm{tr}\big{(}\mathrm{Re}[\mathbf{V}^{\mathrm{H}}\mathbf{G}_{\textnormal{\tiny{RX}}}\mathbf{H}\mathbf{G}_{\textnormal{\tiny{TX}}}\mathbf{W}]\big{)}.$
(29)

Finally, the approximate MSE in (22) is obtained by replacing
$\mathbf{G}_{\textnormal{\tiny{RX}}}$ with
$\tilde{\mathbf{G}}_{\textnormal{\tiny{RX}}}$ in (20) and
$\mathbf{C}_{\mathbf{r}}$ with $\tilde{\mathbf{C}}_{\mathbf{r}}$ in (21). ∎

## References

  * [1] N. Rajatheva, I. Atzeni, E. Björnson _et al._ , “White paper on broadband connectivity in 6G,” 2020. [Online]. Available: http://jultika.oulu.fi/files/isbn9789526226798.pdf
  * [2] S. Buzzi and C. D’Andrea, “Energy efficiency and asymptotic performance evaluation of beamforming structures in doubly massive MIMO mmWave systems,” _IEEE Trans. Green Commun. and Netw._ , vol. 2, no. 2, pp. 385–396, 2018.
  * [3] I. Atzeni, A. Tölli, and G. Durisi, “Low-resolution massive MIMO under hardware power consumption constraints,” in _Proc. Asilomar Conf. Signals, Syst., and Comput. (ASILOMAR)_ , 2021.
  * [4] A. K. Saxena, I. Fijalkow, and A. L. Swindlehurst, “Analysis of one-bit quantized precoding for the multiuser massive MIMO downlink,” _IEEE Trans. Signal Process._ , vol. 65, no. 17, pp. 4624–4634, 2017.
  * [5] Y. Li, C. Tao, G. Seco-Granados, A. Mezghani, A. L. Swindlehurst, and L. Liu, “Channel estimation and performance analysis of one-bit massive MIMO systems,” _IEEE Trans. Signal Process._ , vol. 65, no. 15, pp. 4075–4089, 2017.
  * [6] S. Jacobsson, G. Durisi, M. Coldrey, U. Gustavsson, and C. Studer, “Throughput analysis of massive MIMO uplink with low-resolution ADCs,” _IEEE Trans. Wireless Commun._ , vol. 16, no. 6, pp. 1304–1309, 2017.
  * [7] A. Li, C. Masouros, A. L. Swindlehurst, and W. Yu, “1-bit massive MIMO transmission: Embracing interference with symbol-level precoding,” _IEEE Commun. Mag._ , vol. 59, no. 5, pp. 121–127, 2021.
  * [8] S. Jacobsson, G. Durisi, M. Coldrey, T. Goldstein, and C. Studer, “Quantized precoding for massive MU-MIMO,” _IEEE Trans. Commun._ , vol. 65, no. 11, pp. 4670–4684, 2017.
  * [9] I. Atzeni and A. Tölli, “Channel estimation and data detection analysis of massive MIMO with 1-bit ADCs,” _IEEE Trans. Wireless Commun._ , vol. 21, no. 6, pp. 3850–3867, 2022.
  * [10] Y. Nam, H. Do, Y. Jeon, and N. Lee, “On the capacity of MISO channels with one-bit ADCs and DACs,” _IEEE J. Sel. Areas Commun._ , vol. 37, no. 9, pp. 2132–2145, 2019.
  * [11] A. Bazrafkan and N. Zlatanov, “Asymptotic capacity of massive MIMO with 1-bit ADCs and 1-bit DACs at the receiver and at the transmitter,” _IEEE Access_ , vol. 8, pp. 152 837–152 850, 2020.
  * [12] A. Lozano, “1-bit MIMO for terahertz channels,” in _Proc. Int. ITG Workshop Smart Antennas (WSA)_ , 2021.
  * [13] Ö. T. Demir and E. Björnson, “The Bussgang decomposition of nonlinear systems: Basic theory and MIMO extensions [lecture notes],” _IEEE Signal Process. Mag._ , vol. 38, no. 1, pp. 131–136, 2020.
  * [14] E. Björnson, L. Sanguinetti, and J. Hoydis, “Hardware distortion correlation has negligible impact on UL massive MIMO spectral efficiency,” _IEEE Trans. Commun._ , vol. 67, no. 2, pp. 1085–1098, 2019\.
  * [15] A. M. Sayeed, “Deconstructing multiantenna fading channels,” _IEEE Trans. Signal Process._ , vol. 50, no. 10, pp. 2563–2579, 2002.


# Topology Optimization for Manufacturing with Accessible Support Structures

Amir M. Mirzendehdel, Morad Behandish, and Saigopal Nelaturi  Palo Alto
Research Center (PARC), 3333 Coyote Hill Road, Palo Alto, California 94304

###### Abstract

Metal additive manufacturing (AM) processes often fabricate a near-net shape
that includes the as-designed part as well as the sacrificial support
structures that need to be machined away by subtractive manufacturing (SM),
for instance multi-axis machining. Thus, although AM is capable of generating
highly complex parts, the limitations of SM due to possible collision between
the milling tool and the workpiece can render an optimized part non-
manufacturable. We present a systematic approach to topology optimization (TO)
of parts for AM followed by SM to ensure removability of support structures,
while optimizing the part’s performance. A central idea is to express the
producibility of the part from the near-net shape in terms of accessibility of
every support structure point using a given set of cutting tool assemblies and
fixturing orientations. Our approach does not impose any artificial
constraints on geometric complexity of the part, support structures, machining
tools, and fixturing devices. We extend the notion of inaccessibility measure
field (IMF) to support structures to identify the inaccessible points and
capture their contributions to non-manufacturability by a continuous spatial
field. IMF is then augmented to the sensitivity field to guide the TO towards
a manufacturable design. The approach enables efficient and effective design
space exploration by finding nontrivial complex designs whose near-net shape
can be 3D printed and post-processed for support removal by machining with a
custom set of tools and fixtures. We demonstrate the efficacy of our approach
on nontrivial examples in 2D and 3D.

###### keywords:

Design for Manufacturing , Topology Optimization , Support Structures,
Accessibility Analysis , Multi-Axis Machining , Configuration Space , Hybrid
Manufacturing

††journal: Computer-Aided Design

##  1 Introduction

Many additive manufacturing (AM) technologies require sacrificial support
structures for various reasons such as carrying unsupported weight before
solidification and prevent “drooling” (e.g., fused deposition modeling),
restraining the part to a substrate (e.g., stereo-lithography), or disposing
the excess heat generated by intensive laser power to prevent residual
stresses, thermal distortion, burning, and warping in metal AM (e.g., powder-
bed fusion) [1]. Support structures can increase the AM cost directly and
indirectly, especially in metal-based processes where material, fabrication
time and energy, and post-fabrication cleanup costs are significant. For
instance, recent studies have shown that depending on the heat treatment and
porosity of support structures, support removal by machining can cause
significant tool wear [2]. Nonetheless, support removal via milling remains
popular in post-processing of metal AM parts [3, 4].

Despite many efforts in eliminating the support structures by optimizing the
build orientation or the part, in many practical cases it is not feasible to
fabricate AM parts without support structures. In metal AM, the sacrificial
supports need to be removed by subtractive manufacturing (SM) such as multi-
axis computer numerically-controlled (CNC) machining. This common “hybrid”
(AM-then-SM) workflow is illustrated in Fig. 1.

Topology optimization (TO) is widely used as an automated design tool to
search the expanded design space of AM parts for high-performance and light-
weight designs. TO commonly generates parts of complex geometry with intricate
features that can pose manufacturing challenges. The novel concept of
inaccessibility measure field (IMF) was introduced in [5, 6] to locally
quantify one’s ability to access different points of an intermediate design
using a set of cutting tool assemblies and fixturing orientations. The
continuous IMF was used to penalize the TO sensitivity field, steering the
updates to improve not only physics-based performance criteria (e.g.,
stiffness), but also manufacturability (e.g., machining). In this paper, we
extend this TO framework from a pure-SM to an AM-then-SM workflow in which the
optimization process is constrained to generate designs whose corresponding
near-net shape is accessible at every support structure point for a given set
of cutting tool assemblies and fixturing orientations. The near-net shape is
defined as the super-set of the final as-manufactured shape and its
corresponding support structures. Our approach applies to parts, support
structures, tool assemblies, fixture, and build platforms of arbitrary shape,
as well as high-axis machining degrees of freedom.

Figure 1: Metal AM is typically a hybrid AM-then-SM process. The figure shows
(a) an as-designed model; (b) an intermediate step in the AM process (e.g.,
powder-bed fusion [7]) at a sub-optimal build orientation (selected
arbitrarily for illustration); (c) the near-net including both the target part
and support structured; (d) an intermediate step in the SM process (e.g.,
multi-axis machining with a vise fixture holding the build platform) to remove
the support structures; and (e) the final as-manufactured part.

###  1.1 Related work

In this section, we review recent advances in design for manufacturing (DfM),
particularly design for AM (DfAM), and hybrid (additive followed by
subtractive) manufacturing (DfHM) with a focus on TO as an automated design
tool.

Over the past few years, TO has become increasingly popular in industrial
engineering applications to recommend viable topologies in early stages of
part design [8, 9]. The main challenge in adoption of TO for detailed/final
design is to fill the gap between the ‘as-designed’ and ‘as-manufactured’
models that results in a costly trial-and-error cycle.

Considering both traditional and modern manufacturing processes, numerous DfM
and DfAM approaches have been developed, where typically the manufacturing
constraints are simplified under reasonable assumptions, reduced to straight-
forward design rules, and are imposed through filtering/projection methods
[10]. Traditional manufacturing constraints include ensuring minimum feature
size [11, 12, 13], casting [14, 15], laser-cutting/extrusion [16, 17], and
milling [18, 19, 20]. For instance, proposed approaches for incorporating
multi-axis machining constraint in TO typically assume simple geometries for
machine tools, ignore surrounding fixturing devices, and abstract tool motion
by tracing straight lines. Recently, a general approach was presented [5, 6]
to incorporate accessibility constraints for high-axis machining via
continuous measures, using well-established configuration space
(${\mathsf{C}}-$space) concepts widely used in spatial planning and robotics
[21], with no restrictive assumptions on geometric complexity (of part and
tool), tool motion, or machining setup. On the other hand, DfAM methods have
mainly focused on reducing support material by optimizing the part geometry
[22, 23, 24, 25, 26], optimizing support structures [27, 28, 29], optimizing
build orientation [30, 31], TO for anisotropic materials [32, 33], minimizing
residual stresses [34, 35, 36], and geometric post-processing [37, 38].

DfHM approaches are still fairly limited and focus on AM-then-SM processes.
For example, Liu et al. [39] proposed a TO method by imposing a casting
constraint on free-form AM design domain and feature fitting on shape-
preserved segments for machined surfaces. Lou et al. [40] presented a TO
formulation based on a nonlinear virtual temperature method to ensure that
inaccessible enclosed voids are self-supporting, while allowing the free
boundary to have overhang surfaces. No accessibility analysis with respect to
tool geometry or machining setup was considered. More detailed reviews of
manufacturing-oriented TO can be found in [41, 8, 42].

###  1.2 Contributions & Outline

This paper proposes a TO method for hybrid AM-then-SM manufacturing, where
given a prescribed build orientation and multi-axis machining setup, the near-
net shape corresponding to the optimized design is guaranteed to have
removable supports. The multi-axis machining setup comprises a collection of
machining tool assemblies, fixturing devices, and motions including
translations and rotations. There is no artificial and restricting assumptions
on geometric complexity of tools and fixtures, or motion. Our proposed
manufacturability-constrained optimization incorporates a sophisticated
manufacturing constraint at early stages of design to drastically reduce
downstream trial and error costs. More specifically, the contributions of this
paper are:

  1. 1.

Extending the ‘inaccessibility measure field’ to AM-then-SM hybrid
manufacturing to capture accessibility of support structures.

  2. 2.

Formulating a density-based TO framework that incorporates accessibility
constraint for support removal given a realistic multi-axis machining setup
with multiple geometrically complex cutting tool assemblies, fixturing
devices, and arbitrary motions.

  3. 3.

Proposing a support removal planning algorithm to gradually remove supports
until the part is detached from the platform.

  4. 4.

Demonstrating the effectiveness of our method by solving multiple benchmark
and realistic examples in 2D and 3D.

##  2 Proposed Method

In this section, we will extend the TO formulation, specifically density-based
TO, to incorporate accessibility of support structures given a multi-axis
machining setup.

###  2.1 Topology Optimization Formulation

The design problem can be formulated as the following constrained optimization
problem:

$\displaystyle\mathop{\text{Minimize}}\limits_{\Omega\subseteq\Omega_{0}}\quad$
$\displaystyle\varphi(\Omega),$ (1a) such that
$\displaystyle[\textbf{K}_{\Omega}][\textbf{u}_{\Omega}]=[{\mathbf{f}}],$ (1b)
$\displaystyle V_{\Omega}\leq V_{\text{target}},$ (1c) $\displaystyle
V_{\Gamma}/V_{S}\leq\epsilon,$ (1d)

where the following definitions and notations are used:

  * $\Omega_{0}\subset{\mathds{R}}^{3}$: 3D design domain;

  * $\Omega\subseteq\Omega_{0}$: 3D pointset representing the (intermediate or final) design, to be updated iteratively;

  * $\rho_{\Omega}:\Omega_{0}\to[0,1]$: The density field that implicitly represents the design as a superlevel-set;

  * $\phi(\Omega)\in{\mathds{R}}$: The objective function value (e.g., compliance) for a given design;

  * $[{\mathbf{K}}_{\Omega}]\in{\mathds{R}}^{n\times n}$: Finite element stiffness matrix for a voxelization of the domain into $n$ elements;

  * $[{\mathbf{u}}_{\Omega}]\in{\mathds{R}}^{n}$: Displacement vector for the voxelization;

  * $[{\mathbf{f}}]\in{\mathds{R}}^{n}$: External force vector for the voxelization;

  * $V_{\text{target}}\in{\mathds{R}}$: Target volume upperbound;

  * $V_{\Omega}\in{\mathds{R}}$: The current design’s volume;

  * $O\subset{\mathds{R}}^{3}$: Obstacle geometry defined as the union of part, platform, and fixturing devices ($O=\Omega~{}\cup~{}P~{}\cup~{}F$), where there must be no collision between cutting tools and obstacle $O$;

  * $V_{S}\in{\mathds{R}}$: Support structure volume;

  * $V_{\Gamma}\in{\mathds{R}}$: Secluded (inaccessible) support volume; and

  * $\epsilon>0$: A small tolerance for non-manufacturability (e.g., to enable convergence in the presence of numerical inaccuracies).

In Section 2.2, we explain the accessibility analysis required to find the
secluded supports $\Gamma$ in terms of the cutting tool assemblies $T$,
obstacle $O$, and support structures $S$.

Without considering the support accessibility constraint in (1d), the
constrained TO problem defined in (1a) through (1c) can be expressed as the
following Lagrangian minimization problem:

$\mathcal{L}_{\Omega}:=\varphi(\Omega)+~{}\lambda_{1}(\dfrac{V_{\Omega}}{V_{\text{target}}}-1)+[\lambda_{2}]^{\mathrm{T}}\Big{(}[\textbf{K}_{\Omega}][{\mathbf{u}}_{\Omega}]-[{\mathbf{f}}]\Big{)}.$
(2)

Taking the derivative w.r.t. design variables denoted by the prime symbol
$(\cdot)^{\prime}$ and using the chain rule we have:

$\displaystyle\mathcal{L}_{\Omega}^{\prime}$
$\displaystyle=\varphi^{\prime}(\Omega)+\lambda_{1}\dfrac{{V}_{\Omega}^{\prime}}{V_{\text{target}}}+[\lambda_{2}]^{\mathrm{T}}\Big{(}[\textbf{K}_{\Omega}][{\mathbf{u}}_{\Omega}]\Big{)}^{\prime},$
(3)
$\displaystyle=\Big{(}[\dfrac{\partial\varphi}{\partial{\mathbf{u}}}]+[\lambda_{2}]^{\mathrm{T}}[\textbf{K}_{\Omega}]\Big{)}[{\mathbf{u}}_{\Omega}^{\prime}]$
$\displaystyle+\lambda_{1}\dfrac{{V}_{\Omega}^{\prime}}{V_{\text{target}}}+[\lambda_{2}]^{\mathrm{T}}[\textbf{K}_{\Omega}^{\prime}][{\mathbf{u}}_{\Omega}].$
(4)

(4) reduces to (5) via the adjoint method [43]:

$\displaystyle\mathcal{L}_{\Omega}^{\prime}=\lambda_{1}\dfrac{{V}_{\Omega}^{\prime}}{V_{\text{target}}}+[\lambda_{2}]^{\mathrm{T}}[\textbf{K}_{\Omega}^{\prime}][{\mathbf{u}}_{\Omega}],$
(5)
$\displaystyle\text{if}~{}[\lambda_{2}]:=-[\textbf{K}_{\Omega}]^{-1}[\dfrac{\partial\varphi}{\partial{\mathbf{u}}}].$

Here the objective function is compliance, i.e.,
$\varphi(\Omega):=[{\mathbf{f}}]^{\mathrm{T}}[{\mathbf{u}}_{\Omega}]$, thus
$[\lambda_{2}]=-[{\mathbf{u}}_{\Omega}]$.

###  2.2 Accessibility Analysis for Support Structures

In this section, we briefly discuss the IMF as a means to quantify collision
at each point ${\mathbf{x}}\in\Omega_{0}$ given a design $\Omega$, a set of
tool assemblies $T\subset{\mathds{R}}^{3}$, sampled cutting points
$K\subset{\mathds{R}}^{3}$ and rotations $\Theta\subset{\mathrm{SO}(3)}$,
platform $P\subset{\mathds{R}}^{3}$, and fixturing devices
$F\subset{\mathds{R}}^{3}$. In particular, we review the convolution-based
formulation of IMF, more details can be found in [6].

We define the collision between obstacle $O=\Omega\cup P\cup F$ and cutting
tool $T$ under rotation $R\in\Theta$ and translation ${\mathbf{t}}$ as:

$\mathsf{vol}\big{[}O\cap(R,{\mathbf{t}})T\big{]}=\langle{\mathbf{1}}_{O},{\mathbf{1}}_{(R,{\mathbf{t}})T}\rangle=({\mathbf{1}}_{O}\ast\tilde{{\mathbf{1}}}_{RT})({\mathbf{t}}),$
(6)

where $\mathsf{vol}[\cdot]$ denotes the volume of a 3D pointset which can be
expressed as the inner product of the indicator functions of the obstacle
${\mathbf{1}}_{O}$ and the transformed (rotated and translated) tool
${\mathbf{1}}_{(R,{\mathbf{t}})T}$. It has been shown in [44] that the inner
product can be written as a convolution for objects in relative translation,
namely, the obstacle $O$ and the rotated tool $RT$. The convolution operator
$\ast$ is defined as:

$({\mathbf{1}}_{O}\ast\tilde{{\mathbf{1}}}_{RT})({\mathbf{t}})=\int_{O}{\mathbf{1}}_{O}({\mathbf{x}}){\mathbf{1}}_{-RT}({\mathbf{t}}-{\mathbf{x}})dv[{\mathbf{x}}],$
(7)

where
${\mathbf{1}}_{-RT}({\mathbf{t}})={\mathbf{1}}_{RT}(-{\mathbf{t}})=\tilde{{\mathbf{1}}}_{RT}({\mathbf{t}})$
is the reflection of the rotated tool about the origin, hence
${\mathbf{1}}_{-RT}({\mathbf{t}}-{\mathbf{x}})={\mathbf{1}}_{RT}({\mathbf{x}}-{\mathbf{t}})=\tilde{{\mathbf{1}}}_{RT}({\mathbf{t}}-{\mathbf{x}})$
is the indicator function of the moving object (i.e., rotated tool assembly)
translated to the query points ${\mathbf{x}}\in O$. The convolution of (7) can
be efficiently evaluated via Fast Fourier Transforms (FFT).

As illustrated in Fig. 2, we define the IMF as a mapping
$f_{\text{IMF}}:{\mathds{R}}^{3}\to{\mathds{R}}$ that assigns a scalar value
to each point ${\mathbf{x}}\in{\mathds{R}}^{3}$, where for a given tool
assembly $T=(H\cup K)$ at orientation $R$ we take a pointwise minimum of
shifted convolutions for different choices of sharp points in the cutter
${\mathbf{k}}\in K$:

$f_{\text{IMF}}({\mathbf{x}};O,R,T,K):=\min_{{\mathbf{k}}\in
K}~{}{\mathsf{vol}\big{[}O\cap(R,{\mathbf{x}})(T-{\mathbf{k}})\big{]}}.$ (8)

Figure 2: Consider a 2D part and a 2D tool assembly. At a given query point
${\mathbf{x}}\in\Omega_{0}$, the IMF is computed by looking at different
rotations, e.g., $R_{1}\in{\mathrm{SO}(2)}$ of the tool assembly. For the
oriented tool $T\to R_{1}T$, the origin is shifted to different sharp points
${\mathbf{k}}_{1},{\mathbf{k}}_{2},{\mathbf{k}}_{3},\ldots\in K$; e.g., for
the (rotated) sharp point $R_{1}{\mathbf{k}}_{j}\in R_{1}K$, the tool is
translated as $R_{1}T\to(R_{1}T-R_{1}{\mathbf{k}}_{j})$. The new origin is
brought in contact with the query point, hence
$(R_{1}T-R_{1}{\mathbf{k}}_{j})\to(R_{1}T-R_{1}{\mathbf{k}}_{j})-{\mathbf{x}}$.
This is repeated for all candidate sharp points and the IMF is computed as the
minimum over all rotations and sharp points.

Since, the maximum collision volume happens when the entire tool assembly $T$
collides with the obstacle $O$, $f_{\text{IMF}}$ defined in (8) is normalized
so that its value is within [0,1] everywhere. The normalized IMF value denoted
as $\bar{f}_{\text{IMF}}$ can be written as:

$\bar{f}_{\text{IMF}}({\mathbf{x}};O,R,T,K):=\min_{{\mathbf{k}}\in
K}~{}\dfrac{\mathsf{vol}\big{[}O\cap(R,{\mathbf{x}})(T-{\mathbf{k}})\big{]}}{\mathsf{vol}\big{[}T\big{]}}.$
(9)

To use (9) in a density-based TO loop, we define the function
$\rho_{O}:O\to[0,1]$ as:
$\rho_{O}({\mathbf{x}}):=\rho_{\Omega}({\mathbf{x}})+{\mathbf{1}}_{P}({\mathbf{x}})+{\mathbf{1}}_{F}({\mathbf{x}})$,
in which $\rho_{\Omega}({\mathbf{x}})$ is obtained directly from TO. 111Also,
note that the implicit assumption here is that the fixture does not collide
with neither the platform nor the part – otherwise the function will not be
within [0,1]. 222The + sign indicates disjoint union in the limit when
$\rho_{\Omega}=1_{\Omega},\text{ and }1_{\Omega}+1_{P}+1_{F}=1_{(\Omega\cup
P\cup F)}$.

The IMF in terms of intermediate densities is defined as:

$\bar{f}_{\text{IMF}}({\mathbf{x}};\rho_{O},T,K):=\min_{{\mathbf{k}}\in
K}~{}\dfrac{(\rho_{O}\ast\tilde{{\mathbf{1}}}_{RT})({\mathbf{x}}-R{\mathbf{k}})}{\mathsf{vol}\big{[}T\big{]}}.$
(10)

Since $\rho_{O}\in[0,1]$,
$0\leq\bar{f}_{\text{IMF}}({\mathbf{x}},\rho_{O})\leq\bar{f}_{\text{IMF}}({\mathbf{x}};O)\leq
1$. It is also worth noting that the (6) relates “volume” of collision to a
convolution and is only valid for indicator functions. Thus, the convolution
of (10) is not exactly a volumetric measure, however it is bounded by (9) and
as the TO approaches a solid-void design, the values of the two equations
converge. The procedure for computing IMF for a rotated tool is described in
Algorithm 1.

Algorithm 1 Compute $[\bar{f}_{\text{IMF}}]$ for a rotated tool.

procedure IMF($[\rho_{O}],R,[{\mathbf{1}}_{T}],[{\mathbf{1}}_{K}]$)

$[{\mathbf{1}}_{RT}]\leftarrow\textsc{Rotate}([{\mathbf{1}}_{T}],R)$
$\triangleright$ Re-sampling

$[\tilde{{\mathbf{1}}}_{RT}]\leftarrow\textsc{Reflect}([{\mathbf{1}}_{RT}])$

$[g]\leftarrow\textsc{Convolve}([\rho_{O}],[\tilde{{\mathbf{1}}}_{RT}])$
$\triangleright$ FFT-based

Initialize $[{f}_{\text{IMF}}]\leftarrow[0]$

for all ${\mathbf{k}}\in K$ do

$[h]\leftarrow\textsc{Translate}([g],-R{\mathbf{k}})$

$[{f}_{\text{IMF}}]\leftarrow\min([f_{\text{IMF}}],[h])$ $\triangleright$ Over
sharp points

end for

$[\bar{f}_{\text{IMF}}]\leftarrow[{f}_{\text{IMF}}]/\mathsf{vol}[T]$

return ($[\bar{f}_{\text{IMF}}]$)

end procedure

Next, given available orientations $\Theta\subset{\mathrm{SO}(3)}$ for each
tool assembly, $f_{\text{IMF}}$ for each tool can be written as:

$\bar{f}_{\text{IMF}}({\mathbf{x}};\rho_{O}):=\min_{R\in\Theta}\bar{f}_{\text{IMF}}({\mathbf{x}};\rho_{O},T,K)$
(11)

Subsequently, the combined IMF for all $n_{T}$ tool assemblies is computed as:

$\bar{f}_{\text{IMF}}({\mathbf{x}};\rho_{O}):=\min_{1\leq i\leq
n_{T}}\bar{f}_{\text{IMF}}({\mathbf{x}};\rho_{O},T_{i},K_{i})$ (12)

IMF of (12) can be computed as described in Algorithm 2.

Algorithm 2 Compute overall $[\bar{f}_{\text{IMF}}]$ of (12).

procedure
IMF($[\rho_{\Omega}],[{\mathbf{1}}_{P}],[{\mathbf{1}}_{F}],[{\mathbf{1}}_{H_{i}}],[{\mathbf{1}}_{K_{i}}],\\{\Theta_{i}\\};n_{T}$)

Define
$[\rho_{O}]\leftarrow[\rho_{\Omega}]+[{\mathbf{1}}_{P}]+[{\mathbf{1}}_{F}]$
$\triangleright$ Disjoint union

Initialize $[\bar{f}_{\text{IMF}}]\leftarrow[0]$ $\triangleright$ IMF for all
the tools

for $i\leftarrow 1$ to $n_{T}$ do

Define
$[{\mathbf{1}}_{T_{i}}]\leftarrow[{\mathbf{1}}_{H_{i}}]+[{\mathbf{1}}_{K_{i}}]$
$\triangleright$ Implicit union

for all $R\in\Theta_{i}$ do

$\gamma\leftarrow\textsc{IMF}([\rho_{O}],R,[{\mathbf{1}}_{T_{i}}],[{\mathbf{1}}_{K_{i}}])$
$\triangleright$ Alg. 1

$[\bar{f}_{\text{IMF}}]\leftarrow\min([\bar{f}_{\text{IMF}}],\gamma)$
$\triangleright$ (12)

end for

end for

return ($[\bar{f}_{\text{IMF}}]$)

end procedure

###  2.3 Topology Optimization for Accessible Supports

To incorporate the support accessibility constraint for multi-axis machining,
we modify the sensitivity field $\mathcal{S}_{\Omega}$ as follows:

$\mathcal{S}_{\Omega}:=(1-w_{\text{acc}})~{}\bar{\mathcal{S}}_{\varphi}+w_{\text{acc}}~{}\bar{\mathcal{S}}_{\text{IMF}},$
(13)

where $0\leq w_{\text{acc}}<1$ is the filtering weight for accessibility. In
this paper, $w_{acc}$ gradually increases from 0 to 0.5. This allows the TO to
start from a close-to-optimal design w.r.t. performance and gradually modify
it w.r.t. support accessibility. $\bar{\mathcal{S}}_{\varphi}$ is the
normalized sensitivity field with respect to the objective function, noting
that the volume constraint is satisfied with the optimality criteria iteration
[45]. $\bar{\mathcal{S}}_{\text{IMF}}$ is the normalized accessibility filter
defined in terms of the normalized IMF over the design obtained with the
physical density $\tilde{\rho}_{\Omega}$. The physical density
$\tilde{\rho}_{\Omega}$ is obtained from the Heaviside projection [46]:

$\tilde{\rho}_{\Omega}({\mathbf{x}})=1-e^{-\beta\rho_{\Omega}({\mathbf{x}})}+\rho{}_{\Omega}({\mathbf{x}})e^{-\beta},$
(14)

Our initial TO experiments showed that a layer-wise penalization of the
support inaccessibility measure field helps the optimizer to converge to a
better solution. The heuristic essentially gives preference in retaining
material close to platform if performance sensitivity and collision values are
comparable. Further, since the constraint is imposed as a soft constraint, in
some pathological cases a small amount of inaccessible supports might remain
towards the end of optimization, where there may be low-density regions with a
small non-zero IMF value, for which the sensitivity filter is not enough to
make the barely-secluded regions fully accessible. In this case, we continue
optimization while also imposing a filter directly on density values by
gradually increasing the density at the secluded regions. In our experiments,
this helps TO to converge to a valid solution within a few iterations.

$\bar{\mathcal{S}}_{\text{IMF}}({\mathbf{x}}):=\begin{dcases}-\bar{L}^{q}f_{\text{IMF}}({\mathbf{x}};\tilde{\rho}_{O})&\text{if}~{}{\mathbf{x}}\in\Omega\cup
S,\\\ 0&\text{otherwise}.\end{dcases}$ (15)

where for the $k^{th}$ layer of $n_{z}$ layers,
$\bar{L}_{k}^{q}=\left(\dfrac{n_{z}-k+1}{n_{z}}\right)^{q}$. For all examples,
we used a layer-wise penalization coefficient of $q=4$. We use $\beta:=1$ for
2D and 3D examples of Section 3, which in our implementation seems to provide
a good physical density with sufficient smoothness for computing IMF.

Algorithm 3 describes the TO with support accessibility constraint, where b
and $\alpha$ denote the build direction and the overhang angle for support
generation in AM.

Algorithm 3 TO with multi-axis machining constraint for support removal.

procedure TO($\Omega_{0},V_{target},\textbf{b},\alpha,q$)

Initialize
$[\rho_{\Omega}]\leftarrow[V_{\text{target}}/\textsc{Integral}([{\mathbf{1}}_{\Omega_{0}}])]$

Initialize $\Delta\leftarrow\infty$

Initialize $iter\leftarrow 0$

Initialize $w_{acc}\leftarrow 0$

${S_{0}}\leftarrow\textsc{\text{SuppGen}}(\Omega_{0},\textbf{b},\alpha)$
$\triangleright$ Support gen.

$[{\mathbf{1}}_{N_{0}}]\leftarrow[{\mathbf{1}}_{\Omega_{0}}][{\mathbf{1}}_{S_{0}}]$
$\triangleright$ Initial near-net shape

$\bar{L}^{q}\leftarrow\textsc{LayerCoefs}([{\mathbf{1}}_{N_{0}}])$
$\triangleright$ $\bar{L}_{k}=\dfrac{n_{z}-k+1}{n_{z}}$

while $\Delta>\delta$ and $iter<I$ do

$[\tilde{\rho}_{\Omega}]\leftarrow\textsc{Heaviside}([\rho_{\Omega}],\beta)$
$\triangleright$ Projection

${S}\leftarrow\textsc{\text{SuppGen}}([\tilde{\rho}_{\Omega}],b,\alpha)$
$\triangleright$ Support gen.

$[{\mathbf{1}}_{N}]\leftarrow[{\mathbf{1}}_{\Omega}][{\mathbf{1}}_{S}]$
$\triangleright$ Near-net shape

$[{\mathbf{u}}]\leftarrow\textsc{\text{FEA}}([\tilde{\rho}_{\Omega}],[{\mathbf{f}}])$
$\triangleright$ Solve FEA

$\varphi(\Omega)\leftarrow\textsc{\text{Evaluate}}([\rho_{\Omega}],[{\mathbf{u}}])$
$\triangleright$ Obj. func.

$[\bar{\mathcal{S}}_{\varphi}]\leftarrow\textsc{\text{Gradient}}([\rho_{\Omega}],[{\mathbf{u}}])$$\triangleright$
Sensitivity

$[{f}_{\text{IMF}}]\leftarrow\textsc{IMF}([\tilde{\rho}_{\Omega}],[{\mathbf{1}}_{F}],\ldots;n_{T})$

$\triangleright$ Call Algorithm 2 with obvious arguments

$[{\mathbf{1}}_{\Gamma}]\leftarrow[{\mathbf{1}}_{S}][{f}_{\text{IMF}}>\lambda]$

$[\bar{\mathcal{S}}_{\text{IMF}}]\leftarrow[{\mathbf{1}}_{N}][-\bar{L}^{q}{f}_{\text{IMF}}]$

$[\mathcal{S}]\leftarrow(1-w_{\text{acc}})[\bar{\mathcal{S}}_{\varphi}]+w_{\text{acc}}[\bar{\mathcal{S}}_{\text{IMF}}]$
$\triangleright$ (13)

$[\rho_{\Omega}^{\text{new}}]\leftarrow\textsc{\text{Update}}([\rho_{\Omega}],[\mathcal{S}])$
$\triangleright$ Update density

$\Delta\leftarrow\textsc{Integrate}([\rho_{\Omega}^{\text{new}}]-[\rho_{\Omega}])$
$\triangleright$ Vol. diff.

$iter\leftarrow iter+1$ $\triangleright$ Iter. counter

$[\rho_{\Omega}]\leftarrow[\rho_{\Omega}^{\text{new}}]$ $\triangleright$ For
next iteration

if $iter>I_{acc}$ then $\triangleright$ e.g., $I_{acc}=20)$

$w_{acc}\leftarrow w_{acc}+0.01$ $\triangleright$ Update weight

end if

if $iter>I_{\rho}$ then $\triangleright$ e.g., $I_{\rho}=150$

$\rho_{\Gamma}=\min(\rho_{\Gamma}+0.5,1.0)$$\triangleright$ Penalize density
at secluded regions

$\bar{\mathcal{S}}_{\Gamma}=\min(-\bar{L}^{q},\bar{\mathcal{S}}_{\Gamma})$
$\triangleright$ Penalize sensitivity at secluded regions

end if

end while

return $[\rho_{\Omega}]$

end procedure

##  3 Results

In this section, we will present benchmark and realistic examples in 2D and
3D. All results are generated using a density-based implementation, where the
optimality criteria method was used to update the density field.
All examples are run on a desktop machine with Intel Corei7-7820X CPU with 8
processors running at 4.5 GHz, 32 GB of host memory, and an NVIDIA GeForce GTX
1080 GPU with 2,560 CUDA cores and 8 GB of device memory.

###  3.1 Benchmark Example: Cantilever Beam

First, we consider a simple cantilever beam example in 2D and 3D. The loading
conditions are shown in Fig. 3. We use material properties of Stainless Steel
with Young’s modulus of $E=270~{}GPa$ and Poisson ratio of $\nu=0.3$. In 2D,
we set the volume fraction to 0.5 and resolution at 100$\times$50\. In 3D, the
volume fraction is 0.3 and the resolution at 25$\times$51$\times$25.

Figure 3: Cantilever beam boundary conditions in 2D and 3D.

In each example, we solve both accessibility-constrained and unconstrained TO
and report the nonzero secluded support volume for the latter
$V_{\Gamma_{unc}}>0$ whose prevention comes at the cost of an increase in the
resulting compliance values ($\phi_{con}>\phi_{unc}$). An example of TO with
no accessibility constraint on supports is shown in Fig. 4a where the build
direction is along +Z, and the milling tool can only approach the part from
left ($\Theta=\\{0\\}$). Figures 4b and 4c illustrate the overhang points and
support structures of the optimized design with $90^{\circ}$ overhang angle,
meaning each point needs to be supported by some material directly underneath
it. Figure 4d shows the normalized IMF, which is used to identify accessible
and secluded support structures depicted in Fig. 4e. Observe that the platform
is also considered as an obstacle. Figure 4f shows the filtered IMF
$\bar{\mathcal{S}}$ defined in (15).

Figure 4: TO for cantilever beam in 2D without accessibility constraint on
support structures for $90^{\circ}$ overhang angle: (a) optimized cantilever
beam at volume fraction of 0.5 and the oriented end-mill tool, (b) overhang
points, (c) support structures (d) normalized IMF, (e) accessible and secluded
supports, and (f) filtered IMF.

The 2D results for the cantilever beam example under different build and tool
orientations are illustrated in Fig. 5. Performance and support accessibility
results are summarized in Table 2. As the number of tool orientations
increase, the performance typically improves. However, for design of Fig. 5g
the compliance is slightly higher than Fig. 5a and 5e, since in 2D the
optimizer cannot exploit the higher degrees of freedom to generate a
drastically different design and it converges to a slightly worse locally
optimum design as the problem is not convex. Density and compliance
convergence plots are shown in Figures 6 and 7, respectively. Observe that the
support accessibility constraint is imposed after a few (e.g., 20) iterations
and the graphs for with and without accessibility constraint coincide up to
that point. As $w_{acc}$ gradually increases similar to continuation method
widely used for density penalization, the constrained design deviates from the
nominal unconstrained one to find a better performance (typically higher
compliance or ${\varphi_{\text{con}}}/{\varphi_{\text{unc}}}\geq 1$) with no
secluded supports. Figure 8 compares the values of total support volume and
secluded support volume throughout the optimization. As expected, the support
structure volume can drastically change from one design to another as overhang
surfaces are created or removed.

Figure 5: TO for cantilever beam in 2D with accessibility constraint on
support structures for $90^{\circ}$ overhang angle: (a) $0^{\circ}$ tool angle
with $+Z$ build orientation, (b) $180^{\circ}$ tool angle with $+Z$ build
orientation, (c) $0^{\circ}$ tool angle with $-Z$ build orientation, (d)
$180^{\circ}$ tool angle with $-Z$ build orientation, (e)
$\left\\{0^{\circ},180^{\circ}\right\\}$ tool angle with $+Z$ build
orientation, (f) $\pm 45^{\circ}$ tool angle with $+Z$ build orientation, (g)
$\left\\{0^{\circ},180^{\circ},\pm 45^{\circ},\pm 135^{\circ},\pm
90^{\circ}\right\\}$ tool angle with $+Z$ build orientation. Figure 6: Density
convergence plots for cantilever beam TO of Fig. 5a without and with support
accessibility constraint for $90^{\circ}$ overhang angle. Figure 7: Compliance
convergence plots for cantilever beam TO of Fig. 5a without and with support
accessibility constraint.  Figure 8: Volume of accessible and inaccessible
support structure for cantilever beam TO of Fig. 5a.

The overhang angle in previous examples were conservatively chosen to be
$90^{\circ}$. Let us consider the square cantilever beam example of Fig. 9a,
where the build direction is along +Z and
$[\Theta]=\\{R(\theta)~{}|~{}\theta=2\pi
i/8~{}\text{for}~{}i=0,1,\ldots,7\\}$. The unconstrained optimized design is
shown in Fig. 9b. Observe that there are some isolated regions that are
classified as accessible where there is no connected path for the cutting tool
from outside of the part to that location. As was previously mentioned, the
convolution-based IMF only checks for necessary condition of accessibility;
that is whether the tool would collide with surrounding obstacles if cutter is
placed at the point of interest. A sufficient condition would be to check for
existence of a connected path for the tool to reach the point of interest. In
the optimization, this usually is resolved in subsequent iterations as the
surrounding secluded regions are penalized and become solid. In other words,
as the boundary of the design evolves, these isolated accessible regions
become secluded and correctly penalized until no such region exists, as
illustrated in Fig. 9c with $\alpha=90^{\circ}$ and Fig. 9d with
$\alpha=45^{\circ}$. Although we have not explicitly imposed any self-
supporting constraint, the TO has automatically created a self-supporting
enclosed void.

Figure 9: Difference of optimized designs for different overhang angles. (a)
problem setup, (b) near-net optimized design without accessibility constraint,
(c) near-net optimized design with accessibility constraint at $90^{\circ}$
overhang angle, and (d) near-net optimized design with accessibility
constraint at $45^{\circ}$ overhang angle.

Figure 10 illustrates the voxel representation of the optimized cantilever
beam example in 3D with no support accessibility constraint.

(a)

(b)

Figure 10: Voxel representation of the optimized cantilever beam with no
support accessibility constraint at 0.3 volume fraction.

Figure 11a illustrates the optimized cantilever beam example in 3D with no
support accessibility constraint and the corresponding IMF. The cutting tool
is approaching the part from the left, build direction is along +Z, and the
platform is also considered as an obstacle. Figure 11b shows the optimized
design with accessible and secluded support structures where
$\dfrac{V_{\Gamma_{unc}}}{V_{\mathsf{S}_{unc}}}=0.43$. The 3D results for the
cantilever beam example under different build and tool orientations are
illustrated in Fig. 12. Performance and support accessibility results are
summarized in Table 2. Observe that in 3D, as the number of available tool
orientations increase, the performance of the optimized design improves. In
the case of Fig. 12f, there are 14 tool orientations and a very small amount
of support in the unconstrained design is secluded
($\dfrac{V_{\Gamma_{unc}}}{V_{\mathsf{S}_{unc}}}$ is almost 0). Nonetheless,
imposing the additional accessibility constraint slightly changes the design
and compliance is slightly increased
($\dfrac{\varphi_{\text{con}}}{\varphi_{\text{unc}}}=1.12$).

Table 3 provides a benchmark study on computational time for FEA, support
generation, and IMF computation at different voxel resolutions. IMF is based
on FFT algorithm that is implemented on GPU and is evaluated very efficiently
(Fractions of a second).

Figure 11: TO for cantilever beam in 3D at $0.3$ volume fraction without
accessibility constraint on support structures for $45^{\circ}$ overhang
angle: (a) IMF by $0^{\circ}$ tool angle with $+Z$ build orientation, (b)
accessible and secluded supports. Figure 12: TO for cantilever beam in 3D with
accessibility constraint on support structures for $45^{\circ}$ overhang
angle: (a) $R_{X}(0)$ tool angle with $+Z$ build orientation, (b) $R_{X}(\pi)$
tool angle with $+Z$ build orientation, (c) $R_{X}(0)$ tool angle with $-Z$
build orientation, (d) $R_{X}(\pi)$ tool angle with $-Z$ build orientation,
(e) $\left\\{R_{X}(\pm\pi/4),R_{Z}(\pm\pi/4),\right\\}$ tool angle with $+Z$
build orientation, (f)
$\\{R_{X}(0),R_{X}(\pi),R_{X}(\pm\pi/2),R_{X}(\pm\pi/4),R_{X}(\pm
3\pi/4),R_{Z}(\pm\pi/2),$
$R_{Z}(\pm\pi/4),R_{Z}(\pm 3\pi/4)\\}$ tool angle with $+Z$ build orientation.
Table 1: Summary of cantilever beam results in 2D.

$b$ | $\Theta$ | $\dfrac{\varphi_{\text{con}}}{\varphi_{\text{unc}}}$ | $\dfrac{V_{\Gamma_{unc}}}{V_{\mathsf{S}_{unc}}}$
---|---|---|---
$(0,+1)$ | $\\{R(0)\\}$ | 1.53 | 0.62
$(0,+1)$ | $\\{R(\pi)\\}$ | 2.66 | 1.00
$(0,-1)$ | $\\{R(0)\\}$ | 1.93 | 0.72
$(0,-1)$ | $\\{R(\pi)\\}$ | 3.33 | 0.70
$(0,+1)$ | $\\{R(0),R(\pi)\\}$ | 1.37 | 0.62
$(0,+1)$ | $\\{R(\pm\pi/4)\\}$ | 1.89 | 0.74
$(0,+1)$ | $\\!\begin{aligned} \\{&R(0),R(\pi),R(\pm\pi/2),\\\ &R(\pm\pi/4),R(\pm 3\pi/4)\end{aligned}$ | 1.62 | 0.62

Table 2: Summary of cantilever beam results in 3D.

$b$ | $\Theta$ | $\dfrac{\varphi_{\text{con}}}{\varphi_{\text{unc}}}$ | $\dfrac{V_{\Gamma_{unc}}}{V_{\mathsf{S}_{unc}}}$
---|---|---|---
$(0,0,+1)$ | $\\{R_{X}(0)\\}$ | 2.23 | 0.43
$(0,0,+1)$ | $\\{R_{X}(\pi)\\}$ | 5.13 | 0.91
$(0,0,-1)$ | $\\{R_{X}(0)\\}$ | 3.42 | 0.64
$(0,0,-1)$ | $\\{R_{X}(\pi)\\}$ | 10.00 | 0.79
$(0,0,+1)$ | $\\{R_{X}(\pm\pi/4),R_{Z}(\pm\pi/4)\\}$ | 1.18 | 0.029
$(0,0,+1)$ | $\\!\begin{aligned} \\{&R_{X}(0),R_{X}(\pi),\\\ &,R_{X}(\pm\pi/2),R_{X}(\pm\pi/4),\\\ &R_{X}(\pm 3\pi/4),R_{Z}(\pm\pi/2),\\\ &R_{Z}(\pm\pi/4),R_{Z}(\pm 3\pi/4)\\}\end{aligned}$ | 1.12 | 0.00

Table 3: Computational time of FEA, support generation, and IMF computation. Resolution | FEA (s) | Supp. Gen. (s) | IMF (s)
---|---|---|---
25$\times$51$\times$25 | 1.74 | 0.0006 | 0.009
50$\times$100$\times$50 | 12.13 | 0.002 | 0.05
75$\times$150$\times$75 | 39.9 | 0.007 | 0.1

###  3.2 GE Bracket

Next, let us consider the GE-bracket problem of Fig. 13. The material is
similar to the previous example. The target volume fraction is 0.3 and the TO
resolution is about 80,000 hexahedral finite elements. In all examples, we
consider the build platform as an additional obstacle in IMF calculation.
Figure 14a depicts the optimized design without support accessibility
constraint built in +Z direction with about 9% of supports being secluded. The
optimized design with accessibility constraint is shown in Fig. 14b with no
secluded supports, while compliance increases about 11% and the overall
support ratio $\dfrac{V_{\mathsf{S}_{con}}}{V_{\mathsf{S}_{unc}}}=0.65$.
Figure 15 illustrates the optimized bracket without and with accessibility
constraint with build direction along -Y and 3 tool orientations. The design
with no secluded support has 10% higher compliance value than the
unconstrained case and
$\dfrac{V_{\mathsf{S}_{con}}}{V_{\mathsf{S}_{unc}}}=0.74$. Figure 16
illustrated the optimized design build along +X direction with two tool
orientations. Imposing accessibility constraint results in 10% increase in
compliance while the overall support volume reduces by about 7.5%. Results for
GE bracket examples are summarized in Table 4.

Figure 13: GE bracket geometry and boundary conditions.

(a) Without accessibility constraint.

(b) With accessibility constraint.

Figure 14: Optimized GE bracket at 0.3 volume fraction and $+Z$ build
direction.

(a) Without accessibility constraint.

(b) With accessibility constraint.

Figure 15: Optimized GE bracket at 0.3 volume fraction and $-Y$ build
direction.

(a) Without accessibility constraint.

(b) With accessibility constraint.

Figure 16: Optimized GE bracket at 0.3 volume fraction and $+X$ build direction. Table 4: Summary of GE bracket results. $b$ | Tool Directions | $\dfrac{\varphi_{\text{con}}}{\varphi_{\text{unc}}}$ | $\dfrac{V_{\Gamma_{unc}}}{V_{\mathsf{S}_{unc}}}$
---|---|---|---
$(0,0,+1)$ | $\begin{dcases}(0,+1,0)\\\ (+1,0,0)\\\ (-1,0,0)\end{dcases}$ | 1.11 | 0.09
$(0,-1,0)$ | $\begin{dcases}(0,0,-1)\\\ (+1,0,0)\\\ (-1,0,0)\end{dcases}$ | 1.10 | 0.06
$(+1,0,0)$ | $\begin{dcases}(0,+1,0)\\\ (-1,0,0)\end{dcases}$ | 1.10 | 0.05

###  3.3 Quad-Copter

Figure 17 illustrates the geometry and simplified loading condition on a quad-
copter considering lift and drag forces while retaining some battery-housing
case and upper frame of the copter. The underlying material is assumed to be
Aluminum 7075 with Young’s modulus $E=70GPa$ and $\nu=0.33$. The target volume
fraction is 0.35 and the design domain is discretized into about 200,000
hexahedral elements.

The optimized result without support accessibility constraint is shown in Fig.
18, where the build direction $\textbf{b}=\left\\{0.3,0.3,1.0\right\\}$ and
the set of 14 tool orientations:

$\displaystyle\Theta=\\{$ $\displaystyle R_{X}(0),R_{X}(\pi),$
$\displaystyle,R_{X}(\pm\pi/2),R_{X}(\pm\pi/4),$ $\displaystyle R_{X}(\pm
3\pi/4),R_{Z}(\pm\pi/2),$ $\displaystyle R_{Z}(\pm\pi/4),R_{Z}(\pm 3\pi/4)\\}$

Despite a large volume of seemingly unnecessary supports, this setup ensures
that the initial design domain is manufacturable, i.e., all supports are
accessible. Figure 18a depicts the optimized design, accessible supports, and
secluded supports in the absence of support constraint where the volume of
secluded supports is about 35 $cm^{3}$. Figure 18b illustrates the optimized
design with accessibility constraint, where there is no secluded support
region. The compliance values between the two scenarios are quite comparable,
$\dfrac{\varphi_{\text{con}}}{\varphi_{\text{unc}}}=1.02$. Figure 19
illustrated different views of the final designs with secluded support region.

Figure 17: Quad-copter geometry and boundary conditions.

(a) Without accessibility constraint.

(b) With accessibility constraint.

Figure 18: Optimized quad-copter at 0.35 volume fraction and multi-axis
milling tool. (a) without accessibility constraint, and (b) with accessibility
constraint.

(a) Without accessibility constraint.

(b) With accessibility constraint.

Figure 19: Optimized quad-copter bracket at 0.35 volume fraction.

###  3.4 Upright

Finally, let us consider the upright example of Fig. 20 with relative pressure
values applied on different surfaces. The material properties are those of
Ti6Al4V, with $E=113.8$ GPa and $\nu=0.34$. The design domain is discretized
into about 100,000 hexahedral finite elements.

Figure 21 illustrates the additive and multi-axis machining setup where the
build directions is along $\textbf{b}=\left\\{0.25,0.25,0.85\right\\}$, a VISE
fixture holding the platform, and 3 different cutting tools approaching the
part from different angles. Specifically, there is a total of 26 tool
orientations with:

$\displaystyle\Theta_{1}=\\{$ $\displaystyle R_{X}(0),R_{X}(\pi),$
$\displaystyle,R_{X}(\pm\pi/2),R_{X}(\pm\pi/4),$ $\displaystyle R_{X}(\pm
3\pi/4),R_{Z}(\pm\pi/2),$ $\displaystyle R_{Z}(\pm\pi/4),R_{Z}(\pm 3\pi/4)\\}$
$\displaystyle\Theta_{2}=\\{$ $\displaystyle R_{X}(\pm\pi/2),R_{X}(-\pi/4),$
$\displaystyle R_{Z}(\pm 3\pi/4),-0.88\pi\\}$ $\displaystyle\Theta_{3}=\\{$
$\displaystyle R_{X}(0),R_{X}(\pi),$ $\displaystyle
R_{X}(\pm\pi/2),R_{Z}(\pm\pi/2),$

The tools are initially oriented along +Y. Figure 22 illustrates different
views of the optimized upright without and with accessibility constraint. The
unconstrained TO results in 7.14 $cm^{3}$, while the accessibility-constrained
design exhibits an 18% increase in compliance.

We have established that IMF ensures that a support removal plan exists for
the optimized design shown in Fig. 22b for the machining setup shown in Figure
20. Although, the focus of this paper is not on generating a support removal
plan, we propose a simple strategy to gradually remove support structures from
the near-net shape assuming the part is only fixturable through the platform:

  1. 1.

Initialize the near-net shape and the machining setup

  2. 2.

Starting from the top layer along build direction, find the tool and
orientation that removes the maximum volume from that layer based on:

    1. (a)

Compute IMF for each rotated tool via Algorithm 1

    2. (b)

Take $\tau$ sub-levelset of IMF to find accessible regions

    3. (c)

Select orientated tool if accessible volume is the largest

  3. 3.

Update the near-net shape

  4. 4.

If there is no supports in this layer, go to next layer

Note that the obstacle used to compute IMF for TO does not include support
structure, since we are assuming that as long the tools do not collide with
the part, platform, or fixture, there is a sequence by which the surrounding
support structures can be removed to make a particular support region
accessible. In other words, in TO formulation we are only interested in
whether the support can be removed not how. On the other hand, to generate
valid plans, obstacle for support removal planning is the union of the
supports, part, platform, and fixturing device. Thus, performing task 2 for a
single layer (along build direction) and small $\tau$ results in a very
conservative plan that may be inefficient and cumbersome. The method can be
relaxed by allowing removal from multiple layers and slightly larger $\tau$ to
tolerate some level of collision with the support structures to allow
volumetric removal in each step. To penalize collision with part or fixture,
we assign higher values for these regions when we compute the convolution of
(7). Here, we allow for removal of material from 10% of total number of layers
at each step with $\tau=0.005$ and penalizing the non-support regions by a
factor of 1000.

Figure 23 illustrates 8 sampled intermediate steps is support removal
planning, where different tools are used to machine as much support structures
as possible while ensuring that the part remains connected to the platform
until all supports that are in contact with the part are removed. In this
case, overall 63% of the supports need to be machined.

Figure 20: Upright geometry and boundary conditions.

(a) Without accessibility constraint.

(b) With accessibility constraint.

Figure 21: Optimized upright bracket at 0.35 volume fraction, with build
direction of (0.25, 0.25,0.85), VISE fixture holding the build platform, and
three milling tools and 26 tool orientations (tool 1: 14, tool 2: 6, and tool
3: 6 orientations).

(a) Without accessibility constraint.

(b) With accessibility constraint.

Figure 22: Optimized upright bracket at 0.35 volume fraction.

(a) Tool 2 (0.14%)

(b) Tool 1 (0.79%)

(c) Tool 3 (2.97%)

(d) Tool 3 (5.82%)

(e) Tool 3 (1.05%)

(f) Tool 1 (1.76%)

(g) Tool 2 (0.70%)

(h) Tool 2 (1.23%)

Figure 23: Sampled steps of support removal plan where supports are gradually
removed (green) from top to bottom of the current near-net shape (gray) until
the part is detached from the platform. In this case, overall 63% of the
supports need to be machined.

##  4 Conclusion

We presented a general approach to ensure removability of support structures
upfront at the early stages of design for a hybrid AM-then-SM process. In
particular, we extend the concept of ‘inaccessibility measure field’ (IMF)
introduced in [6] to near-net shapes produced by many metal additive
processes. IMF is a continuous field defined using well-established
mathematical concepts in spatial planning in terms of convolutions in
configuration space.

We extend the density-based TO formulation to incorporate accessibility of
supports given a multi-axis machining setup comprised of a collection of
realistic tool assemblies and fixturing devices. We define a normalized
support accessibility filter that is augmented to the performance sensitivity
field by an adaptive weight factor $w_{acc}$. To improve convergence, for the
first few iterations (say 20), the accessibility constraint is ignored
($w_{acc}=0$). Once an initial approximation of the design is obtained,
$w_{acc}$ is gradually increased up to an upper-bound (say 0.5). Our
preliminary experiments showed that a layer-wise penalization of the support
inaccessibility measure field helps the optimizer to converge to a better
design. Further, since the constraint is imposed as a soft constraint, in some
pathological cases a small amount of inaccessible supports remain towards the
end of optimization. In this case we continue optimization while also imposing
a filter directly on density values. In our experiments, this helps TO to
converge to a valid solution within a few iterations. The effectiveness of our
proposed method is demonstrated through benchmark and realistic examples in 2D
and 3D.

In this work, the part remains fixturable throughout the optimization process
since we assumed that it is sufficient for the fixturing device to only the
platform. This may be limiting in support removal planning, where in many
cases we need to hold the workpiece itself or design a soft-jaw to ensure
support removal without damaging the part. Future work will focus on
incorporating a general fixturability constraint in TO.

## References

  * [1] J. Jiang, X. Xu, J. Stringer, Support structures for additive manufacturing: a review, Journal of Manufacturing and Materials Processing 2 (4) (2018) 64.
  * [2] V. Tripathi, A. Armstrong, X. Gong, G. Manogharan, T. Simpson, E. De Meter, Milling of inconel 718 block supports fabricated using laser powder bed fusion, Journal of Manufacturing Processes 34 (2018) 740–749.
  * [3] C. Höller, T. Hinterbuchner, P. Schwemberger, P. Zopf, R. Pichler, F. Haas, Direct machining of selective laser melted components with optimized support structures, Procedia Cirp 81 (2019) 375–380.
  * [4] W. Hintze, R. von Wenserski, S. Junghans, C. Möller, Finish machining of ti6al4v slm components under consideration of thin walls and support structure removal, Procedia Manufacturing 48 (2020) 485–491.
  * [5] A. M. Mirzendehdel, M. Behandish, S. Nelaturi, Exploring feasible design spaces for heterogeneous constraints, Computer-Aided Design 115 (2019) 323–347.
  * [6] A. M. Mirzendehdel, M. Behandish, S. Nelaturi, Topology optimization with accessibility constraint for multi-axis machining, Computer-Aided Design 122 (2020) 102825.
  * [7] V. Bhavar, P. Kattire, V. Patil, S. Khot, K. Gujar, R. Singh, A review on powder bed fusion technology of metal additive manufacturing, Additive manufacturing handbook (2017) 251–253.
  * [8] J. Liu, A. T. Gaynor, S. Chen, Z. Kang, K. Suresh, A. Takezawa, L. Li, J. Kato, J. Tang, C. C. Wang, et al., Current and future trends in topology optimization for additive manufacturing, Structural and Multidisciplinary Optimization 57 (6) (2018) 2457–2483.
  * [9] T. Dbouk, A review about the engineering design of optimal heat transfer systems using topology optimization, Applied Thermal Engineering 112 (2017) 841–854.
  * [10] J. Norato, B. Bell, D. A. Tortorelli, A geometry projection method for continuum-based topology optimization with discrete elements, Computer Methods in Applied Mechanics and Engineering 293 (2015) 306–327.
  * [11] O. Sigmund, Morphology-based black and white filters for topology optimization, Structural and Multidisciplinary Optimization 33 (4-5) (2007) 401–424.
  * [12] J. K. Guest, Imposing maximum length scale in topology optimization, Structural and Multidisciplinary Optimization 37 (5) (2009) 463–473.
  * [13] M. Zhou, B. S. Lazarov, F. Wang, O. Sigmund, Minimum length scale in topology optimization by geometric constraints, Computer Methods in Applied Mechanics and Engineering 293 (2015) 266–282.
  * [14] L. Harzheim, G. Graf, A review of optimization of cast parts using topology optimization, Structural and multidisciplinary optimization 31 (5) (2006) 388–399.
  * [15] Q. Li, W. Chen, S. Liu, H. Fan, Topology optimization design of cast parts based on virtual temperature method, Computer-Aided Design 94 (2018) 28–40.
  * [16] K. Ishii, S. Aomura, Topology optimization for the extruded three dimensional structure with constant cross section, JSME International Journal Series A Solid Mechanics and Material Engineering 47 (2) (2004) 198–206.
  * [17] V.-N. Hoang, H. Nguyen-Xuan, Extruded-geometric-component-based 3d topology optimization, Computer Methods in Applied Mechanics and Engineering 371 (2020) 113293.
  * [18] M. Langelaar, Topology optimization for multi-axis machining, Computer Methods in Applied Mechanics and Engineering (2019).
  * [19] A. T. Gaynor, T. E. Johnson, Eliminating occluded voids in additive manufacturing design via a projection-based topology optimization scheme, Additive Manufacturing (2020) 101149.
  * [20] N. Morris, A. Butscher, F. Iorio, A subtractive manufacturing constraint for level set topology optimization, Structural and Multidisciplinary Optimization (2020) 1–16.
  * [21] T. Lozano-Perez, Spatial planning: A configuration space approach, in: Autonomous robot vehicles, Springer, 1990, pp. 259–271.
  * [22] A. M. Mirzendehdel, K. Suresh, Support structure constrained topology optimization for additive manufacturing, Computer-Aided Design 81 (2016) 1–13.
  * [23] M. Langelaar, Topology optimization of 3d self-supporting structures for additive manufacturing, Additive Manufacturing 12 (2016) 60–70.
  * [24] J. Wu, C. C. Wang, X. Zhang, R. Westermann, Self-supporting rhombic infill structures for additive manufacturing, Computer-Aided Design 80 (2016) 32–42.
  * [25] X. Qian, Undercut and overhang angle control in topology optimization: a density gradient based integral approach, International Journal for Numerical Methods in Engineering 111 (3) (2017) 247–272.
  * [26] G. Allaire, C. Dapogny, R. Estevez, A. Faure, G. Michailidis, Structural optimization under overhang constraints imposed by additive manufacturing technologies, Journal of Computational Physics 351 (2017) 295–328.
  * [27] E. Malekipour, A. Tovar, H. El-Mounayri, Heat conduction and geometry topology optimization of support structure in laser-based additive manufacturing, in: Mechanics of Additive and Advanced Manufacturing, Volume 9, Springer, 2018, pp. 17–27.
  * [28] M. Zhou, Y. Liu, C. Wei, Topology optimization of easy-removal support structures for additive manufacturing, Structural and multidisciplinary optimization (2020).
  * [29] C. Wang, X. Qian, Optimizing support for heat dissipation in additive manufacturing, in: ASME 2020 International Design Engineering Technical Conferences and Computers and Information in Engineering Conference, American Society of Mechanical Engineers Digital Collection, 2020.
  * [30] A. Chandrasekhar, T. Kumar, K. Suresh, Build optimization of fiber-reinforced additively manufactured components, Structural and Multidisciplinary Optimization 61 (1) (2020) 77–90.
  * [31] C. Wang, X. Qian, Simultaneous optimization of build orientation and topology for additive manufacturing, Additive Manufacturing (2020) 101246.
  * [32] A. M. Mirzendehdel, B. Rankouhi, K. Suresh, Strength-based topology optimization for anisotropic parts, Additive Manufacturing 19 (2018) 104–113.
  * [33] P. Zhang, J. Liu, A. C. To, Role of anisotropic properties on topology optimization of additive manufactured load bearing structures, Scripta Materialia 135 (2017) 148–152.
  * [34] L. Cheng, X. Liang, J. Bai, Q. Chen, J. Lemon, A. To, On utilizing topology optimization to design support structure to prevent residual stress induced build failure in laser powder bed metal additive manufacturing, Additive Manufacturing 27 (2019) 290–304.
  * [35] G. Allaire, L. Jakabčin, Taking into account thermal residual stresses in topology optimization of structures built by additive manufacturing, Mathematical Models and Methods in Applied Sciences 28 (12) (2018) 2313–2366.
  * [36] T. Miki, T. Yamada, Topology optimization for considering distortion in additive manufacturing, arXiv preprint arXiv:2008.06153 (2020).
  * [37] M. K. Swierstra, D. K. Gupta, M. Langelaar, Automated and accurate geometry extraction and shape optimization of 3d topology optimization results, arXiv preprint arXiv:2004.05448 (2020).
  * [38] S. C. Subedi, C. S. Verma, K. Suresh, A review of methods for the geometric post-processing of topology optimized models, Journal of Computing and Information Science in Engineering (2020) 1–17.
  * [39] J. Liu, A. C. To, Topology optimization for hybrid additive-subtractive manufacturing, Structural and Multidisciplinary Optimization 55 (4) (2017) 1281–1299.
  * [40] Y. Luo, O. Sigmund, Q. Li, S. Liu, Additive manufacturing oriented topology optimization of structures with self-supported enclosed voids, Computer Methods in Applied Mechanics and Engineering 372 (2020) 113385.
  * [41] J. Liu, Y. Ma, A survey of manufacturing oriented topology optimization methods, Advances in Engineering Software 100 (2016) 161–175.
  * [42] J. Plocher, A. Panesar, Review on design and structural optimisation in additive manufacturing: Towards next-generation lightweight structures, Materials & Design (2019) 108164.
  * [43] M. P. Bendsøe, M. Philip, O. Sigmund, Topology Optimization–Theory, Methods, and Applications, 2nd Edition, Springer-Verlag Berlin Heidelberg, 2009\.  doi:10.1007/978-3-662-05086-6.
  * [44] M. Lysenko, S. Nelaturi, V. Shapiro, Group morphology with convolution algebras, in: Proceedings of the 14th ACM symposium on solid and physical modeling, ACM, 2010, pp. 11–22.
  * [45] O. Sigmund, A 99 line topology optimization code written in matlab, Structural and multidisciplinary optimization 21 (2) (2001) 120–127.
  * [46] E. Andreassen, A. Clausen, M. Schevenels, B. S. Lazarov, O. Sigmund, Efficient topology optimization in Matlab using 88 lines of code, Structural and Multidisciplinary Optimization 43 (1) (2011) 1–16.  doi:10.1007/s00158-010-0594-7.