# Fabrication Sequence Optimization for Minimizing Distortion in Multi-Axis
Additive Manufacturing

Weiming Wang<EMAIL_ADDRESS>Fred van Keulen<EMAIL_ADDRESS>Jun
Wu<EMAIL_ADDRESS>School of Mathematical Sciences, Dalian University of
Technology, Dalian, China  Department of Precision and Microsystems
Engineering, Delft University of Technology, Delft, The Netherlands
Department of Sustainable Design Engineering, Delft University of Technology,
Delft, The Netherlands

###### Abstract

Additive manufacturing of metal parts involves phase transformations and high
temperature gradients which lead to uneven thermal expansion and contraction,
and, consequently, distortion of the fabricated components. The distortion has
a great influence on the structural performance and dimensional accuracy,
e.g., for assembly. It is therefore of critical importance to model, predict
and, ultimately, reduce distortion. In this paper, we present a computational
framework for fabrication sequence optimization to minimize distortion in
multi-axis additive manufacturing (e.g., robotic wire arc additive
manufacturing), in which the fabrication sequence is not limited to planar
layers only. We encode the fabrication sequence by a continuous pseudo-time
field, and optimize it using gradient-based numerical optimization. To
demonstrate this framework, we adopt a computationally tractable yet
reasonably accurate model to mimic the material shrinkage in metal additive
manufacturing and thus to predict the distortion of the fabricated components.
Numerical studies show that optimized curved layers can reduce distortion by
orders of magnitude as compared to their planar counterparts.

###### keywords:

Fabrication sequence , multi-axis additive manufacturing , wire arc additive
manufacturing , thermal distortion , process planning , topology optimization

fn1fn1footnotetext: This work was performed when the first author worked at
the Department of Sustainable Design Engineering, Delft University of
Technology.

##  1 Introduction

Multi-axis additive manufacturing using robotic systems can reorient the print
head or the component under construction, in addition to provide translational
motion [1, 2]. This manufacturing flexibility creates new opportunities and
challenges in process planning. In conventional additive manufacturing, which
only allows for translation, the fabrication sequence is limited to the use of
planar layers. Given a fixed orientation of the component, the stacking of
successive planar layers is uniquely determined. In multi-axis additive
manufacturing, however, depositing material along curved layers becomes
possible [3, 4]. The curved layers of a component can be interpreted as an
ordered decomposition of the component according to a fabrication sequence.
Obviously, the number of possible decompositions is huge. Optimizing the
fabrication sequence for multi-axis additive manufacturing offers a new angle
to address some of the demanding issues in achieving the highest quality by
additive manufacturing. This optimization, however, is a challenging task. It
first needs a representation of the fabrication process as optimization
variables. To be suitable for gradient-based numerical optimization, this
representation is required to be differentiable — the gradient of a functional
with respect to the representation can be calculated. In a recent work [5], we
proposed to encode the continuous material deposition process by a pseudo-time
field. Each point in the component is associated with a pseudo-time. During
the fabrication process, a location with a larger pseudo-time is materialized
later than a location with a smaller pseudo-time. The effectiveness of this
encoding has been validated on compliance minimization under basic
fabrication-process-dependent loads, e.g., gravity [5]. The encoding is
expected to be applicable to more complex fabrication-process-dependent
physics, such as thermomechanical loads and microstructural developments. In
the current work, we make an important step towards this end by demonstrating
fabrication sequence optimization for minimizing distortion of the fabricated
components.

An application scenario that drives our current work is the thermomechanical
distortion in robotic wire arc additive manufacturing (WAAM). In WAAM, a metal
wire is used as feedstock and an electric arc as the energy source [1, 6]. As
a technology for producing metal components, it has a high deposition rate,
low material and equipment costs, and a large build volume. It has great
potential for fabricating large metallic parts that are commonly seen in
aerospace [2, 7], automotive [8], and maritime industries [9]. For instance,
Oak Ridge National Laboratory in the US demonstrated the fabrication of two
fully functional excavator arms that are over 2 meters tall [10, 11]. Dutch
company MX3D printed and installed a 12-meter-long bridge, using 4.5 tonnes of
stainless steel 111https://mx3d.com/. This directed energy deposition process
involves phase transformations and high temperature gradients as the metal
wire is melted and fused to form a component. This results in non-uniform
expansion and contraction of the material, and consequently non-negligible
residual stresses as well as distortion of the fabricated components [1, 2].
These issues have large impacts on the structural performance and dimensional
accuracy, e.g., for assembly. In the related field of welding, the sequence of
joining is found to significantly influence the development of residual
stresses and distortion [12, 13]. For instance, depending on the welding
sequence, the distortion in pipes can vary by more than seven times [14].
Specific to WAAM, investigations of different infill patterns and their
printing orders also confirm their influences on the residual stresses and
distortion [15, 16, 17, 18]. Mitigating residual stresses and minimizing
distortion are important topics of current interest for the deployment of WAAM
technologies [1, 11].

An important module in our fabrication sequence optimization framework is to
predict the distortion of the fabricated component. The modelling and
simulation of additive manufacturing processes is itself an active research
field [19, 20, 21, 22]. A multitude of process models have been proposed, with
different levels of accuracy and computational complexity. WAAM involves time-
dependent thermal and mechanical states along with the successive addition of
melted material. This could be modelled by a fully coupled transient
thermomechanical finite analysis. Such modelling allows gaining insights into
the development of residual stresses and distortion due to the fabrication
sequence as well as deposition rates and heat sink characteristics [16, 23,
24, 25]. Similar to multi-pass welding simulation, WAAM simulation involves
two major steps: a thermal analysis where the temperature distribution caused
by a moving heat source is determined, and a structural analysis which
accounts for a series of sequential load steps calculated in a transient
thermal analysis [26]. This coupled thermomechanical analysis however is
computationally intensive. It is prohibitive for numerical optimization of the
fabrication sequence, since numerical optimization requires simulating the
fabrication process at each optimization iteration. An appealing alternative
is the inherent strain method [27, 28], in which the residual stress and
distortion are obtained by linear elastic analysis. Due to the computational
advantages, inherent strain methods have been increasingly applied for
modelling additive manufacturing processes [29, 30, 31, 32, 33]. Inherent
strain methods have been recently integrated into structural topology
optimization for minimizing distortion. For instance, Allaire and Jakabčin
[34] proposed a level-set based topology optimization method that considers
thermal residual stresses and thermal deformations during the layer-by-layer
fabrication process. Pellens et al. [35] incorporated an inherent strain
method in density-based topology optimization to optimize support structures
for limiting the vertical displacement of layers during fabrication. Also
building upon density-based topology optimization, Misiun et al. [36]
integrated an inherent strain method to optimize the structural layout under
distortion constraints. The inherent strain method has also been used in
level-set based topology optimization for controlling distortion [37].

Our main purpose in this paper is to propose the fabrication sequence as
optimization variables for minimizing process-order-dependent structural
distortion. To this end, we adopt the inherent strain method for its
computational efficiency and adapt it to the context of multi-axis additive
manufacturing. It is chosen as an example to validate the effectiveness of our
fabrication sequence optimization framework. More sophisticated process models
are expected to work as well with our framework but are computationally
expensive. Our current validation considers some essential manufacturability
constraints but is not meant to be comprehensive. Further manufacturability
constraints, e.g., regarding the layer thickness, and interference between the
print head and already deposited part, constitute future work in the direction
opened up by the proposed framework.

The remainder of this paper is organized as follows. In Section 2 we present
the problem formulation and computational framework. We discuss 2D and 3D
numerical examples in Section 3, and conclude the paper in Section 4 with
ideas for future work.

##  2 Fabrication Sequence Optimization for Distortion Minimization

In this section, we start with an overview of the proposed framework, and
proceed to the encoding of the fabrication sequence (Section 2.1). This is
followed by a mechanical simulation model (Section 2.2) and essential
manufacturability constraints (Section 2.3). The formulation of the
optimization problem, together with sensitivity analysis, is given in Section
2.4.

Figure 1 gives an overview of the proposed computational framework. We encode
the fabrication sequence using a pseudo-time field (Fig. 1b). The pseudo-time
field is a continuous one, from which a series of intermediate structures is
calculated (Fig. 1c). In this illustration, the component is decomposed by the
pseudo-time field into 12 equally-sized layers. The distortion (Fig. 1d) is
predicted by simulating the addition of one layer at a time. In an iterative
optimization process, the distortion is minimized by updating the design
variables, i.e., the pseudo-time at each point in the domain of the component.

Figure 1: The computational framework of fabrication sequence optimization for
minimizing distortion. The sequence for an input component (a) is represented
by a continuous pseudo-time field (b), from which intermediate structures (c)
are calculated by projection operators. The simulated distortion (d) is
reduced by updating the pseudo-time field in an iterative optimization
process, leading to minimized distortion (e).

###  2.1 Fabrication sequence

We encode the fabrication sequence for additively manufacturing a component by
a pseudo-time field, $\bm{t}$. Each point in the domain of the component
(specifically, each element in a discretization of the domain) is assigned a
pseudo-time variable, $t_{e}$. It shall be emphasized that this pseudo-time
field mainly indicates the order of materialization: a larger value indicates
an element is to be fabricated later in the process. The pseudo-time field is
a normalized one, i.e., $t_{e}\in[0,1]$. For simplicity, we also refer to it
as the time field. From the time field, layers can be extracted by dividing
the range of time variables (i.e., $[0,1]$) into uniform intervals,
$[T_{j-1},T_{j}]$, $j=1,...,N$, with $T_{j}=\frac{j}{N}$ and $N$ being the
number of layers. By evaluating which interval a time value $t_{e}$ belongs
to, it can be decided which layer the corresponding element is part of. The
time field thus effectively partitions the component into $N$ layers.

Consider a component $\Omega$ in 2- or 3-dimensional space. A time field
partitions the component into $N$ non-overlapping but successive layers,
$\Omega_{i}$, $i=1,...,N$ (see Fig. 2 for illustration on a 2D $L$-shaped
component). The layered fabrication process can also be described by a series
of pseudo-density fields,
$\\{\bm{\rho}^{\\{0\\}},\bm{\rho}^{\\{1\\}},...,\bm{\rho}^{\\{j-1\\}},\bm{\rho}^{\\{j\\}},\bm{\rho}^{\\{j+1\\}},...,\bm{\rho}^{\\{N\\}}\\}$,
where the superscript $\\{j\\}$ indicates the $j$-th stage during
construction, i.e., layers from the first up to the $j$-th have been
fabricated. $\bm{\rho}^{\\{0\\}}$ refers to the beginning of fabrication,
$\rho^{\\{0\\}}_{e}=0$, $\forall e$, while $\bm{\rho}^{\\{N\\}}$ refers to the
end of fabrication, $\rho^{\\{N\\}}_{e}=1$, $\forall e$. At the $j$-th stage,
the density of an element depends on the location of its centroid ($x_{e}$),

$\rho^{\\{j\\}}_{e}=\begin{cases}1,&x_{e}\in\cup^{j}_{i=1}\Omega_{i},\\\
0,&\text{otherwise.}\end{cases}$ (1)

Figure 2: Schematic of curved layers for fabricating an $L$-shaped component.
Left: Layers of the component. Middle and right: The density field at two
consecutive stages.

To facilitate gradient-based numerical optimization, we present a
differentiable operator to partition the component based on the time field.
Let $\bm{1}$ denote the density field of a fully solid component, and assume
the density field and time field are discretized with the same mesh. The
density field consisting of the layers up to the $j$-th one is calculated by
an entry-wise multiplication,

$\bm{\rho}^{\\{j\\}}(\bm{t})=\bm{1}\circ\bar{\bm{t}}^{\\{T_{j}\\}},\quad
j=1,...,N.$ (2)

Here $\bar{\bm{t}}^{\\{T\\}}$ is projected from $\bm{t}$. This projection
converts, in a differentiable manner, time values that are smaller than $T$ to
$1$, and time values that are larger than $T$ to $0$. This continuous
projection is realized by

$\bar{\bm{t}}^{\\{T\\}}=1-\frac{\tanh(\beta_{t}T)+\tanh(\beta_{t}(\bm{t}-T))}{\tanh(\beta_{t}T)+\tanh(\beta_{t}(1-T))},$
(3)

where $\beta_{t}$ is a positive number to control the projection sharpness.
This sharpness parameter, along with a continuation scheme, is commonly used
in density-based topology optimization [38, 39]. Fig. 3 (left) and (middle)
illustrate the projection function at $T_{2}=\frac{2}{5}$ and
$T_{3}=\frac{3}{5}$, respectively, where $5$ is the total number of layers.
The increment corresponding to the $3$-rd layer is decided by
$\Delta\bm{\rho}^{\\{3\\}}=\bar{\bm{t}}^{\\{T_{3}\\}}-\bar{\bm{t}}^{\\{T_{2}\\}}$
(right).

Figure 3: Visualization of the projection function, at $T_{2}=\frac{2}{5}$
(left) and $T_{3}=\frac{3}{5}$ (middle), and the density of the $3$-rd layer,
i.e.,
$\Delta\bm{\rho}^{\\{3\\}}=\bar{\bm{t}}^{\\{T_{3}\\}}-\bar{\bm{t}}^{\\{T_{2}\\}}$
(right).

###  2.2 Mechanical simulation

Our mechanical simulation is based on an inherent strain method. Imagine a
melted material, in isolation from the build plate or previous layers. As it
solidifies, its volume shrinks (illustrated in Fig. 4). Using Voigt notation,
the thermal strain tensor in 2D is expressed by
$\underline{\bm{\varepsilon}}^{*}=[\epsilon_{x},\epsilon_{y},\epsilon_{xy}]^{\mathsf{T}}$.
The quantities of inherent strains can be determined empirically [40] or using
a detailed process simulation [41]. The shrinkage of an element due to the
thermal strain is equivalently caused by nodal forces ${\bm{f}}$ that
correspond to the thermal strain,

${\bm{f}}(\underline{\bm{\varepsilon}}^{*})=[\bm{D}]^{\mathsf{T}}[\bm{C}]\underline{\bm{\varepsilon}}^{*},$
(4)

where $[\bm{C}]$ denotes the elasticity matrix and $[\bm{D}]$ the strain-
displacement relation.

Figure 4: Illustration of the shrinkage of material due to an isotropic
thermal strain

.

Extending from an element to the addition of a new layer, the displacement of
the simulation grid $\Delta\bm{u}^{\\{j\\}}$, induced by the thermal strain of
elements in the $j$-th layer, $\Omega_{j}$, can be computed by solving

$\bm{K}^{\\{j\\}}\Delta\bm{u}^{\\{j\\}}=\bm{f}^{\\{j\\}},\quad j=1,...,N.$ (5)

Here $\bm{K}^{\\{j\\}}$ is the stiffness matrix for the entire domain
$\Omega$, corresponding to the density field $\bm{\rho}^{\\{j\\}}$. The force
vector $\bm{f}^{\\{j\\}}$ consolidates nodal forces from the thermal strain,
$\Delta{\bm{\varepsilon}^{*}}^{\\{j\\}}$. The entries of the thermal strain
vector $\Delta{\bm{\varepsilon}^{*}}^{\\{j\\}}$ are mostly zero, except for
elements that are part of $\Omega_{j}$, i.e.,

$\Delta{\bm{\varepsilon}^{*}_{e}}^{\\{j\\}}=\begin{cases}\underline{\bm{\varepsilon}}^{*},&x_{e}\in\Omega_{j},\\\
[0,0,0]^{\mathsf{T}},&\text{otherwise.}\end{cases}$ (6)

We note that the equilibrium equation Eq. (5) applies to the entire domain,
including both solid elements and the empty grid that has not been filled with
material. This is beneficial from two aspects. Firstly, the empty grid is
deformed along with the solid elements. We exploit this feature to mimic
building a new layer on top of the deformation of previous layers, as in
reality in multi-axis additive manufacturing. Secondly, it facilitates
implementation as the dimension of equilibrium equations during fabrication is
kept constant. The empty elements were referred to as quiet elements in the
literature [42]. To avoid singularity of the stiffness matrix due to the
inclusion of empty elements, a small yet non-zero density, or equivalently, a
minimum stiffness, can be assigned to them.

With the equilibrium equation (Eq. (5)), we can compute the incremental
displacement of the simulation grid due to each addition of a new layer. Under
the assumption of small displacements, the displacement field at the $j$-th
stage, $\bm{u}^{\\{j\\}}$, is computed by the summation of incremental
displacements until the current stage,

$\bm{u}^{\\{j\\}}=\textstyle{\sum^{j}_{i=1}}\Delta\bm{u}^{\\{i\\}},\quad
j=1,...,N.$ (7)

We remark that Eq. (5) can be solved in parallel for the addition of each
layer, since the incremental displacement ($\Delta\bm{u}^{\\{j\\}}$) induced
by the thermal strain of the $j$-th layer
($\Delta{\bm{\varepsilon}^{*}}^{\\{j\\}}$) is independent of the thermal
strain in prior and latter layers. This follows the analysis of a linearity
assumption that was initially proposed for simplifying the simulation of
powder-bed processes [30] with only planar layers.

####  2.2.1 Mechanical simulation with differentiable layers

The mechanical simulation model described above assumes a binary description
of the layers, i.e., each element belongs exclusively to a single layer. In
our framework for fabrication sequence optimization, however, the
decomposition into layers is not binary, but is constructed to be
differentiable with respect to the time field. In other words, an element may
partly belong to multiple layers, with different ratios. The equilibrium
equation (Eq. (5)) in mechanical simulation is thus adapted,

$\bm{K}(\bm{\rho}^{\\{j\\}}(\bm{t}))\Delta\bm{u}^{\\{j\\}}=\bm{f}(\Delta{\bm{\varepsilon}^{*}}^{\\{j\\}}(\bm{t})),\quad
j=1,...,N.$ (8)

The stiffness matrix of an element with a density $\rho_{e}\in[0,1]$ is
interpolated from that of a solid element, $\bm{K}_{0}$, by

$\bm{K}_{e}(\rho_{e})=(E_{min}+(1-E_{min})\rho^{p}_{e})\bm{K}_{0},$ (9)

where $p$ is a penalization parameter ($p=3$). $E_{min}$ is a small value
($E_{\text{min}}=10^{-9}$), introduced to prevent the stiffness matrix from
becoming singular.

For the equivalent forces of the thermal strains, on the right-hand side of
the equilibrium equation, the conditional function in Eq. (6) is replaced by

$\Delta{\bm{\varepsilon}^{*}}^{\\{j\\}}(\bm{t})|_{e}=\left(\Delta\bm{\rho}^{\\{j\\}}(\bm{t})|_{e}\right)^{q}\,\underline{\bm{\varepsilon}}^{*}.$
(10)

Here $\underline{\bm{\varepsilon}}^{*}$ is the thermal strain for an element,
as introduced previously. The increment in the density field corresponding to
the $j$-th layer is given by

$\Delta\bm{\rho}^{\\{j\\}}(\bm{t})=\bm{\rho}^{\\{j\\}}(\bm{t})-\bm{\rho}^{\\{j-1\\}}(\bm{t}),\quad
j=1,...,N,$ (11)

where, as defined previously, $\bm{\rho}^{\\{0\\}}=\bm{0}$,
$\bm{\rho}^{\\{N\\}}=\bm{1}$.

In Eq. (10) the density increment is raised to the power of $q$, similar to
the penalization $p$ in the interpolation of stiffness (Eq. (9)). If the
thermal-strain penalization ($q$) is smaller than the stiffness penalization
($p$), for instance, $p=3$ while $q=1$, the thermal strain creates a large
mechanical load for a disproportionally weak material, leading to an excessive
displacement. Based on numerical experiments, we suggest having $q\geq p$. In
the examples in this paper, we used $q=p$, if not explicitly stated otherwise.

###  2.3 Manufacturability constraints

Manufacturability constraints shall be integrated into the optimization to
avoid sequences that are infeasible from a manufacturing perspective. Here we
present the incorporation of two constraints that we consider essential. Yet
we note that these constraints are by no means exhaustive. These constraints
were formulated in our prior work of space-time topology optimization [5]. In
this current work, these constraints are integrated, with some adaptations. We
briefly reiterate the key formulation with an emphasis on adaptation.

####  2.3.1 Continuity constraints

A local minimum in the time field indicates an isolated material patch during
the fabrication process; all its adjacent elements have a larger time value
and thus have not yet been fabricated. The presence of local minima would
necessitate auxiliary structures during fabrication. A local maximum in the
time field, if not located on the contour of the component, indicates an
enclosed void, which is inaccessible; all its adjacent elements have a smaller
time value and thus have already been fabricated.

Isolated material patches and enclosed voids can be prevented by requiring

$\min_{i\in\mathcal{N}_{e}}(t_{i})\leq
t_{e}\leq\max_{i\in\mathcal{N}_{e}}(t_{i}),\;\forall e\in\mathcal{M},$ (12)

where $\mathcal{N}_{e}$ denotes the set of elements adjacent to element $e$.
The set $\mathcal{M}$ includes all elements in the domain except those that
are prescribed as the starting point/region for the fabrication process (i.e.,
with $t_{e}=0$, typically on the building plate).

The $\min$ and $\max$ functions are not differentiable. A possible approach is
to approximate them by the $p$-norm [5]. In this work, we use an alternative
approach that simplifies implementation. Specifically, we consolidate these
two constraints by an approximation,

$t_{e}=\textrm{mean}_{i\in\mathcal{N}_{e}}(t_{i}),\;\forall e\in\mathcal{M},$
(13)

with

$\textrm{mean}_{i\in\mathcal{N}_{e}}(t_{i})=\frac{\sum_{i\in\mathcal{N}_{e}}t_{i}}{{n}(\mathcal{N}_{e})},$
(14)

here, ${n}$ denotes the number of elements in a set. The mean of numbers in a
set is bounded by its minimum and maximum number. Therefore, this
approximation is conservative, in the sense that it is a sufficient but not
necessary condition of Eq. (12). The mean function ensures the smoothness of
the time field, and also effectively avoids checkerboard patterns that are
typically suppressed by convolution filters in density-based topology
optimization.

The mean equality constraint applies to all elements in the domain except the
starting element. This leads to a large number of constraints. We consolidate
them by

$g_{0}(\bm{t})=\frac{1}{{n}(\mathcal{M})}\sum_{e\in\mathcal{M}}||t_{e}-\textrm{mean}_{i\in\mathcal{N}_{e}}(t_{i})||^{2}\leq\gamma_{c},$
(15)

where $\gamma_{c}$ is a small constant which is set to $10^{-3}$ in our
examples. As $\gamma_{c}$ approaches 0, this constraint effectively restricts
$t_{e}$ towards the mean value of its neighbors. This constraint function is
quadratic, for which sensitivity analysis is straightforward.

####  2.3.2 Volume constraints on layers

For simplicity, we assume a constant fabrication speed, i.e., each layer has
the same material volume. Denote the total volume of the component by $V^{*}$.
Since the number of layers is prescribed as $N$, the volume up until the
$j$-th layer shall be $\frac{j}{N}V^{*}$,

$V^{\\{j\\}}(\bm{t})={\sum_{e}}\,\rho_{e}^{\\{j\\}}(t_{e})\cdot
v_{e}=\frac{j}{N}V^{*},\quad j=1,...,N,$ (16)

where $v_{e}$ is the volume of an element. A uniform finite element
discretization is used in our implementation, and thus $v_{e}$ is a constant
for all elements.

The above equality constraint is replaced by inequality constraints,

$-\gamma_{v}\leq g_{j}(\bm{t})=\frac{V^{\\{j\\}}}{V^{*}}-\frac{j}{N}\leq
0,\quad j=1,...,N,$ (17)

where $\gamma_{v}$ is a small constant for numerical stability ($10^{-3}$ in
our examples). The lower bound is needed in minimizing thermal strain-induced
distortion, since otherwise, a distortion of $0$ can be achieved by a trivial
solution, $\bm{t}=\bm{0}$.

###  2.4 Optimization problem formulation

The optimization aims to reduce the distortion of the fabricated component.
Depending on application scenarios, different distortion measurements can be
applied, for instance, the absolute displacement of a node (or a few nodes),
or the relative displacement between a few selected nodes, e.g., to maintain a
flat edge/surface or to maintain a right angle between edges/surfaces. These
various distortions are generically encoded by
$\bm{u}^{\mathsf{T}}\bm{Q}\bm{u}$, where
$\bm{u}=[\bm{u}^{{\\{1\\}}^{\mathsf{T}}},...,\bm{u}^{{\\{N\\}}^{\mathsf{T}}}]^{\mathsf{T}}$
encompasses the displacement vector during fabrication, and $\bm{Q}$ is a
sparse, symmetric matrix to arrive at different measurements. In the examples
shown in this paper, the objective concerns the displacement vector
$\bm{u}^{{\\{N\\}}}$, i.e., measured when the entire component has been
fabricated.

Our optimization problem is formulated as

$\displaystyle\underset{{\bm{t}}}{\min}$ $\displaystyle
d(\bm{t})=\bm{u}^{\mathsf{T}}\bm{Q}\bm{u}$ (18) $\displaystyle s.t.$
$\displaystyle\bm{K}^{\\{i\\}}(\bm{t})\Delta\bm{u}^{\\{i\\}}=\bm{f}^{\\{i\\}}(\bm{t}),\;i=1,...,N,$
(22) $\displaystyle-\gamma_{v}\leq g_{i}(\bm{t})\leq 0,\;i=1,2,...,N,$
$\displaystyle g_{0}(\bm{t})\leq\gamma_{c},$ $\displaystyle 0\leq t_{e}\leq
1,\;\forall e.$

This optimization problem is solved by gradient-based numerical optimization.
In particular, we use the method of moving asymptotes (MMA) [43].

####  2.4.1 Sensitivity analysis

Introducing Lagrange multipliers for the equilibrium equations,
$\bm{\lambda}_{i}$, $i=1,2,...,N$, the Lagrange function for the objection
function $d(\bm{t})$ is written as:

$\displaystyle
L(\bm{t})=\bm{u}^{\mathsf{T}}\bm{Q}\bm{u}+\sum_{i=1}^{i=N}\bm{\lambda}_{i}^{\mathsf{T}}\left(\bm{K}^{\\{i\\}}(\bm{t})\Delta\bm{u}^{\\{i\\}}-\bm{f}^{\\{i\\}}(\bm{t})\right).$
(23)

Then, the sensitivity of the objective function regarding the design variable
$t_{e}$ is given by

$\displaystyle\frac{\partial L(\bm{t})}{\partial t_{e}}=$ $\displaystyle
2\bm{u}^{\mathsf{T}}\bm{Q}\frac{\partial\bm{u}}{\partial
t_{e}}+\sum_{i=1}^{i=N}\bm{\lambda}_{i}^{\mathsf{T}}\left(\frac{\partial\bm{K}^{\\{i\\}}(\bm{t})}{\partial
t_{e}}\Delta\bm{u}^{\\{i\\}}(\bm{t})+\bm{K}^{\\{i\\}}(\bm{t})\frac{\partial\Delta\bm{u}^{\\{i\\}}(\bm{t})}{\partial
t_{e}}-\frac{\partial\bm{f}^{\\{i\\}}(\bm{t})}{\partial t_{e}}\right),$ (24)
$\displaystyle=$ $\displaystyle
2\bm{u}^{\mathsf{T}}\bm{Q}\sum_{i=1}^{i=N}\frac{\partial\Delta\bm{u}^{\\{i\\}}(\bm{t})}{\partial
t_{e}}+\sum_{i=1}^{i=N}\bm{\lambda}_{i}^{\mathsf{T}}\left(\frac{\partial\bm{K}^{\\{i\\}}(\bm{t})}{\partial
t_{e}}\Delta\bm{u}^{\\{i\\}}(\bm{t})+\bm{K}^{\\{i\\}}(\bm{t})\frac{\partial\Delta\bm{u}^{\\{i\\}}(\bm{t})}{\partial
t_{e}}-\frac{\partial\bm{f}^{\\{i\\}}(\bm{t})}{\partial t_{e}}\right).$

We then reorder the expression on the right-hand side,

$\displaystyle\frac{\partial L(\bm{t})}{\partial t_{e}}=$
$\displaystyle\sum_{i=1}^{i=N}\left(2\bm{u}^{\mathsf{T}}\bm{Q}+\bm{\lambda}_{i}^{\mathsf{T}}\bm{K}^{\\{i\\}}(\bm{t})\right)\frac{\partial\Delta\bm{u}^{\\{i\\}}(\bm{t})}{\partial
t_{e}}+\sum_{i=1}^{i=N}\bm{\lambda}_{i}^{\mathsf{T}}\left(\frac{\partial\bm{K}^{\\{i\\}}(\bm{t})}{\partial
t_{e}}\Delta\bm{u}^{\\{i\\}}(\bm{t})-\frac{\partial\bm{f}^{\\{i\\}}(\bm{t})}{\partial
t_{e}}\right).$ (25)

We choose $\bm{\lambda}_{i}$, satisfying

$\bm{K}^{\\{i\\}}(\bm{t})\bm{\lambda}_{i}+2\bm{Q}^{\mathsf{T}}\bm{u}=0,\;i=1,2,...,N.$
(26)

This simplifies Eq. (25) to

$\displaystyle\frac{\partial L(\bm{t})}{\partial t_{e}}=$
$\displaystyle\sum_{i=1}^{i=N}\bm{\lambda}_{i}^{\mathsf{T}}\left(\frac{\partial\bm{K}^{\\{i\\}}(\bm{t})}{\partial
t_{e}}\Delta\bm{u}^{\\{i\\}}(\bm{t})-\frac{\partial\bm{f}^{\\{i\\}}(\bm{t})}{\partial
t_{e}}\right).$ (27)

The stiffness matrix of element $e$ depends solely on its time value. Using
the chain rule, the derivative becomes

$\frac{\partial\bm{K}^{\\{i\\}}_{e}}{\partial
t_{e}}=\frac{\partial\bm{K}_{e}(\rho^{\\{i\\}}_{e})}{\partial\rho^{\\{i\\}}_{e}}\frac{\partial\rho^{\\{i\\}}_{e}}{\partial\bar{t}^{\\{i\\}}_{e}}\frac{\partial\bar{t}^{\\{i\\}}_{e}}{\partial
t_{e}},$ (28)

where
$\frac{\partial\bm{K}_{e}(\rho^{\\{i\\}}_{e})}{\partial\rho^{\\{i\\}}_{e}}$,
$\frac{\partial\rho^{\\{i\\}}_{e}}{\partial\bar{t}^{\\{i\\}}_{e}}$, and
$\frac{\partial\bar{t}^{\\{i\\}}_{e}}{\partial t_{e}}$ can be derived from Eq.
(9), Eq. (2), and Eq. (3), respectively. Then,
$\frac{\partial\bm{K}^{\\{i\\}}(\bm{t})}{\partial t_{e}}$ in Eq. (27) is
assembled from $\frac{\partial\bm{K}^{\\{i\\}}_{e}}{\partial t_{e}}$ of each
element.

The isotropic thermal strain of element $e$ contributes equivalent loads to
its four nodes in 2D or eight nodes in 3D. From Eqs. (4) and (10), we get

$\frac{\partial\bm{f}^{\\{i\\}}_{e}}{\partial
t_{e}}=q[\bm{D}]^{\mathsf{T}}[\bm{C}]\underline{\bm{\varepsilon}}^{*}\left(\Delta\rho^{\\{i\\}}_{e}\right)^{q-1}\frac{\partial\Delta\rho^{\\{i\\}}_{e}}{\partial
t_{e}},$ (29)

where $\frac{\partial\Delta\rho^{\\{i\\}}_{e}}{\partial t_{e}}$ can be derived
from Eqs. (2) and (11). Then,
$\frac{\partial\bm{f}^{\\{i\\}}(\bm{t})}{\partial t_{e}}$ in Eq. (27) is
assembled from $\frac{\partial\bm{f}^{\\{i\\}}_{e}}{\partial t_{e}}$ of each
element.

The continuity constraints are quadratic functions of $\bm{t}$, while the
volume constraints are linear functions. The sensitivities can be derived
accordingly, and are omitted here.

##  3 Numerical Results

The proposed methods have been implemented using Matlab. In this section, we
report and discuss the numerical results of our tests on 2D and 3D components.
We use bilinear quadrilateral elements for 2D, and trilinear hexahedral
elements for 3D. For demonstration purposes, the Young’s modulus of the
material is 1, and the Poisson’s ratio is 0.3. We use an isotropic thermal
strain in most examples, except in Section 3.4 where we demonstrate an
extension of our framework to take in an anisotropy thermal strain. The
isotropic thermal strain in 2D is set as $\epsilon_{x}=\epsilon_{y}=-0.01$ and
$\epsilon_{xy}=0.0$. Similarly, in 3D, $\epsilon_{i}=-0.01,i\in\\{x,y,z\\}$
and the shear component is 0.

###  3.1 2D $L$-shaped component

We first test our framework on a 2D $L$-shaped component (Fig. 5a). The domain
is discretized by a finite element grid of $120\times 80$. The nodes on the
bottom are fixed, mimicking depositing on a horizontal build plate. We assume
the fabrication starts from the element on the bottom left. A distance field
corresponding to this element is computed to initialize the design field
($\bm{t}$) (Fig. 5c). The prescribed number of layers is $N=8$.

Figure 5: (a) Design domain of an $L$-shaped component. (b) The objective is
to reduce the distortion, measured by the displacement of a node (Test-1), the
flatness of a horizontal edge (Test-2), and the right angle formed by two
edges (Test-3). (c) The fabrication starts from the bottom-left element. The
time field is initialized by the normalized distance field corresponding to
this startpoint.

####  3.1.1 Distortion measures

We demonstrate distortion minimization with different distortion measures, as
illustrated in Fig.5b. This includes minimizing the displacement of a single
node (Test-1), maintaining the flatness of an edge (Test-2), and maintaining
the perpendicularity of two edges (Test-3).

##### Test-1

The displacement of a node is computed by

$d_{1}=(u^{x}_{i})^{2}+(u^{y}_{i})^{2},$ (30)

where subscript $i$ refers to the index of the node of interest (the top-right
node in this test), and the superscripts $x$ and $y$ indicate the $x$ and $y$
components of the displacement vector, respectively.

##### Test-2

The flatness of a horizontal edge is computed by

$d_{2}=\frac{1}{{n}(\mathcal{I})}\sum_{i\in\mathcal{I}}\left(u^{y}_{i}-\textstyle{\frac{1}{{n}(\mathcal{I})}}\textstyle{\sum_{i\in\mathcal{I}}}u^{y}_{i}\right)^{2},$
(31)

where $\mathcal{I}$ denotes the set of nodes that are sampled from the edge of
interest. The number of nodes ${n}(\mathcal{I})$ is equal to or larger than
$2$. The flatness is tested with different numbers of sample nodes.

##### Test-3

The perpendicularity of a horizontal and a vertical edge is measured by the
flatness of these two edges,

$d_{3}=\frac{1}{{n}(\mathcal{I})}\sum_{i\in\mathcal{I}}\left(u^{y}_{i}-\textstyle{\frac{1}{{n}(\mathcal{I})}}\textstyle{\sum_{i\in\mathcal{I}}}u^{y}_{i}\right)^{2}+\frac{1}{{n}(\mathcal{J})}\sum_{j\in\mathcal{J}}\left(u^{x}_{j}-\textstyle{\frac{1}{{n}(\mathcal{J})}}\textstyle{\sum_{j\in\mathcal{J}}}u^{x}_{j}\right)^{2},$
(32)

where $\mathcal{I}$ and $\mathcal{J}$ denote the set of nodes from the
horizontal and vertical edge, respectively.

####  3.1.2 Results

Figure 6 (top) shows the simulated distortion of the fabricated component
using different layers (bottom). The background Cartesian grid is superimposed
for interpreting the distortion. Comparing with the distortion resulting from
horizontal layers (Fig. 6a), it can be observed that the optimized curved
layers lead to a much smaller distortion (b-d). Checking the distortion of
interest, one can find out that in Fig. 6b the top-right node locates close to
its reference position in the Cartesian grid, in (c) the two nodes on the top
form a straight horizontal line, and in (d) the top edge and right edge, each
represented by two nodes, maintain perpendicular to each other. These results
confirm the effectiveness of the proposed method.

Figure 6: (a) Planar layers (bottom) and the corresponding simulated
distortion (top). (b-d) The fabrication sequence is optimized to reduce the
distortion, with three different distortion measures, i.e., the displacement
of the top-right node (Test-1), the flatness of the top edge (Test-2), and the
right angle between the top and right edges (Test-3). Figure 7: Process
simulation with the optimized fabrication sequence (Test-2 in Fig. 6). Figure
8: Optimized sequences (bottom) and the corresponding simulated distortion
(top) for maintaining the flatness of the top edge, with the flatness being
defined by the vertical mismatch of, from left to right, 2, 3, 9, and all
nodes on the top edge.

Figure 7 shows process simulation using the optimized curved layers,
corresponding to Test-2 in Fig. 6. The entire domain of the component,
including both solid elements (i.e., already fabricated part) and the empty
grid, is simulated in each step. The empty grid passively deforms along with
the already fabricated part. While during the first few steps, the top edge in
the passive part deviates from a horizontal line, it gradually aligns
horizontally in the last steps.

Figure 8 shows distortion with the flatness of the top edge defined by
different numbers of sample nodes. The sample nodes are equally spaced on the
top edge. Compared to Test-2 in Fig. 6c, as more nodes are introduced to
define flatness, the bend in the middle of the edge is gradually flattened.

Figure 9: The simulated distortion (left) and optimized fabrication sequence
(middle), tested on a higher resolution and with 20 layers. The history of the
objective value over iterations is plotted on the right.

We further validate our method with a higher resolution and a larger number of
layers. This is shown in Fig. 9 where the resolution is $180\times 120$, and
the number of layers is set to $20$. The results are comparable to the one
shown in Fig. 6c. The plot shown in Fig. 9 (right) confirms the convergence of
the gradient-based numerical optimization. There are some oscillations on the
convergence curve due to a continuation of $\beta_{t}$ in Eq. (3). It starts
from 30 and is increased by 10 every 30 iterations, until it reaches 100. The
same setting is used in all 2D tests.

####  3.1.3 Parameter study

In Fig. 10 we study the influences of the thermal-strain penalty, $q$ in Eq.
(10). The stiffness penalty, $p$ in Eq. (9), is set to $3$. When $q<p$, the
internal force due to the thermal strain is disproportionally larger than the
stiffness, causing large distortion in finite element simulation. This
distortion can be observed in Fig. 10 (left). It also leads to convergence
problems. While the objective was reduced to a small value, one of the volume
constraints in this test was not satisfied. In contrast, when $q\geq p$,
stable optimization is observed. As $q$ changes from 5 to 10, the optimized
results show only subtle differences.

Figure 10: Distortion minimization using different thermal-strain penalties
($q$ in Eq. (10)), and the same stiffness penalty ($p=3$). From left to right:
$q<p$, $q=p$, and $q>p$. Figure 11: Left: Illustration of the smoothed
Heaviside projection (top) used in optimization and a binary 0-1 projection
(bottom) for defining layers. Right: The simulated distortion using
differentiable layers (top) and binary layers (bottom).

Differential layers (i.e., with densities very close, but not equal to 0 or 1)
have been used for simulating distortion in the optimization process. Fig. 11
illustrates the smoothed Heaviside projection (top left) and a strict 0-1
projection (bottom left). The function corresponding to the 4-th layer is
shaded for comparison. To investigate the error due to this approximation, we
perform distortion simulation with both differential layers and binary layers.
The simulation results shown on the right-hand side suggest that the overall
distortion is visually comparable. From the objective value, it is found that
the measured distortion is larger when using the binary layers in a post-
validation. This discrepancy could be reduced by increasing the sharpness of
the smoothed Heaviside projection. Reflecting on the distortion measure (Eq.
(31)), it can be calculated that the vertical mismatch is $0.18$, i.e., less
than a fifth of the unit length of a finite element or $0.15\%$ of the length
of the component. We consider this a rather small distortion.

###  3.2 2D bracket

The second example is a 2D bracket, shown in Fig. 12a. The domain is
discretized by a finite element grid of $144\times 96$. The fabrication is
supposed to start from the corner on the bottom left of the component. The
nodes on the bottom of the component are fixed in the process simulation. The
number of layers is set to 12. The distortion is measured by the average of
the displacements of nodes on the inner circle (red),

$d=\frac{1}{{n}(\mathcal{W})}{\sum_{j\in\mathcal{W}}}\left(\left(u^{x}_{j}\right)^{2}+\left(u^{y}_{j}\right)^{2}\right),$
(33)

where $\mathcal{W}$ denotes the set of nodes along the circle.

Figure 12: The fabrication sequence of a 2D bracket (a) is optimized to
minimize the average displacement of nodes on the inner circle (red). The
fabrication starts from the bottom-left corner. The objective is reduced to
$0.768$ (c) from $10.814$ (d).

Figure 12b shows the optimized fabrication sequence. The resulting distortion
is depicted in (c), while the distortion from regular horizontal layers is
shown in (d) for comparison. The distortion in (d) makes the inner circle
elliptic, while in (c) a circular hole is maintained at its intended location.

We further test the influence of the prescribed starting region of fabrication
on the optimized distortion. In Fig. 13a, the fabrication starts from the
bottom edge of the 2D shape. Fig. 13b depicts the optimized distortion. The
predicted distortion is 0.207. Comparing Fig. 13 with Fig. 12, it can be found
that while the algorithm works effectively with different starting regions,
the starting region has a large influence on the final distortion. Finding an
optimal starting region thus can be an interesting topic for further research.

Figure 13: (a) The optimized fabrication sequence where the fabrication is
supposed to start from the bottom edge of the component. The objective is
reduced to $0.207$ (b) from $10.814$ (c).

###  3.3 3D components

Our framework is directly applicable to 3D scenarios. Fig. 14a illustrates the
domain of a 3D component. It is discretized by a finite element grid of
$60\times 20\times 40$. This component is optimized with 10 layers. The
distortion is measured by the flatness and orientation of the surfaces on the
right and on the top, by calculating the mismatch of nodes along the $x$-axis
and $z$-axis, respectively,

$d=\frac{1}{{n}(\mathcal{I}_{r})}\sum_{i\in\mathcal{I}_{r}}\left(u^{x}_{i}-\textstyle{\frac{1}{{n}(\mathcal{I}_{r})}}\textstyle{\sum_{i\in\mathcal{I}_{r}}}u^{x}_{i}\right)^{2}+\frac{1}{{n}(\mathcal{I}_{t})}\sum_{i\in\mathcal{I}_{t}}\left(u^{z}_{i}-\textstyle{\frac{1}{{n}(\mathcal{I}_{t})}}\textstyle{\sum_{i\in\mathcal{I}_{t}}}u^{z}_{i}\right)^{2},$
(34)

where $\mathcal{I}_{r}$ and $\mathcal{I}_{t}$ are the sets of nodes sampled on
the right and top surfaces, respectively. A comparison of distortion from
optimized curved layers and regular planar layers can be seen in Fig. 14 (b)
and (c). The quantitative measure of distortion is $0.1173$ (b) vs. $1.7075$
(c). Fig. 15 illustrates the simulated sequence of fabrication using the
optimized layers.

Figure 14: Sequence optimization is performed to maintain the flatness and
orientation of the surfaces on the right and on the top of the given component
(a). In process simulation, nodes on the bottom surface (blue) are assumed to
be fixed. The fabrication is assumed to start from a corner on the bottom. The
optimized layers and final distortion can be seen in (b) while the distortion
of planar layers is shown in (c). Figure 15: The mechanical simulation of the
fabrication process using optimized curved layers shown in Fig. 14b.

The last example is a 3D structure designed by topology optimization (Fig.
16a). The domain of the structure is discretized by a finite element grid of
$60\times 40\times 40$. The optimized structure has a volume fraction of
$0.2$. It serves as input to our sequence optimization. For process
simulation, the bottom of the domain is fixed, and one of its corners is
prescribed as the startpoint of the fabrication. The fabrication sequence for
this component is optimized with 10 layers. The objective function is the
flatness of the surfaces on the top and on the left (i.e., similar to Eq.
(34)). The simulated distortion with optimized layers is shown in Fig. 16b. As
a comparison, the distortion from regular horizontal layers is shown in Fig.
16c. The quantitative measure of distortion is $0.0267$ (b) vs. $0.9888$ (c).
The simulation sequence with optimized curved layers is visualized in Fig. 17.

Figure 16: The given component (a) is designed by topology optimization. Its
fabrication sequence is then optimized to maintain the flatness and
orientation of the surfaces on the left and on the top. The distortion is
reduced to $0.0267$ (b) from $0.9888$ (c). Figure 17: The mechanical
simulation of the fabrication process using optimized curved layers shown in
Fig. 16b.

###  3.4 Anisotropic inherent strain

To demonstrate the extendability of our framework, we replace the isotropic
inherent strain with an anisotropic one. In 2D, the anisotropic strain
$\underline{\bm{\varepsilon}}^{*}=[\epsilon_{x},\epsilon_{y},\epsilon_{xy}]^{\mathsf{T}}$
is set as $\epsilon_{x}=-0.01$, $\epsilon_{y}=0$, and $\epsilon_{xy}=0$. Here,
the subscripts $x$ and $y$ refer to a local coordinate system associated with
each element. The $x$-axis represents the material deposition direction at the
element. The orthogonal $y$-axis is defined by the gradient of the time field.
To facilitate the computation of the gradient, we define the time field on the
nodes of the finite element grid. The gradient at each element is then
calculated by the bilinear shape function. We use the normalized gradient of
the time field to transform the anisotropic inherent strain to the global
coordinate system for finite element analysis. This transformation is
incorporated into the calculation of the nodal forces that are equivalent to
the inherent strain (Eq. 4). Accordingly, the sensitivity analysis needs to
consider the dependence of the nodal forces on the time field, Eq. 29 in
particular.

Figure 18 shows the results of fabrication sequence optimization of a 2D
component using an isotropic inherent strain (top row) and an anisotropic one
(bottom row). The 2D domain is discretized using a finite element grid of
$100\times 100$. The number of layers is $20$. The objective is to keep the
top edge of the component horizontal (Eq. 31), measured on three vertices,
i.e., at each end of the edge, and at its middle. When planar layers (a) are
used, the flatness of the top edge is the same for the isotropic and
anisotropic inherent strain (b and e, respectively). However, the average of
the vertical displacements is different, $1.49$ (b) and $2.49$ (e). The
average of vertical displacements is reduced to $-0.80$ (d) and $0.41$ (h),
respectively. In both cases, the distortion objective is significantly
reduced, from $2.57$ to $8.85\times 10^{-3}$ (d) and $2.42\times 10^{-3}$ (h).
This validates the effectiveness of the proposed fabrication sequence
optimization framework for different process models.

Figure 18: Fabrication sequence optimization of a 2D component using an
isotropic inherent strain (top row) and an anisotropic one (bottom row). The
objective is to keep the top edge of the component horizontal. In both cases,
the distortion is significantly reduced. In the isotropic case, the distortion
is reduced from $2.57$ of the planar layers (b) to $8.85\times 10^{-3}$ of the
optimized curved layers (d). In the anisotropic case, an axis of the
anisotropic strain is aligned with the gradient of the time field (g). The
distortion is reduced from $2.57$ of the planar layers (e) to $2.42\times
10^{-3}$ of the optimized curved layers (h).

###  3.5 Computational complexity

The computational complexity of the proposed framework is linear with respect
to the (prescribed) number of layers. Compared to the classic compliance
minimization problem which requires solving a single state equation, here we
need to solve $N$ state equations, with $N$ being the number of layers. Each
state equation corresponds to the deposition of a new layer (cf. Eq. 5).
Additionally, for each state equation, the sensitivity analysis in distortion
minimization involves solving the adjoint state equation (cf. Eq. 26). Thus,
roughly speaking, the complexity is about $2N$ times that of the commonly
solved compliance minimization problem. Under the assumption of small
deformation, the state equations are independent [30] and can thus be solved
in parallel, though in our implementation we did not explore this. As an
indication of the computational time, sequence optimization of the 2D
$L$-shaped component with a finite element grid of $120\times 80$ and 8 layers
in Figs. 6, 8 and 10 took 15.5 minutes on average. The higher resolution 2D
$L$-shaped component with a finite element grid of $180\times 120$ and 20
layers took 84.8 minutes. The running time of the 3D examples shown in Figs.
14 and 16 was slightly less than 15 hours. The tests were based on an
implementation using Matlab, on a PC equipped with Intel(R) Core(TM) i7-6700K
CPU @ 4.00GHz and 64GB memory. A single CPU core was used. All results shown
in this work were generated with 500 optimization iterations.

The number of layers is a prescribed parameter in our framework. It can be
selected by

$N=\frac{V^{*}}{\bar{h}\bar{l}},$ (35)

where $V^{*}$ is the total volume of the component. $\bar{h}$ is the desired
nominal thickness of layers. $\bar{l}$ is the average length of layers, which
can be estimated based on the dimension of the component.

##  4 Conclusions and Future Work

In this paper, we have presented a computational framework for fabrication
sequence optimization to reduce distortion in multi-axis additive
manufacturing (e.g., robotic wire arc additive manufacturing). The
optimization takes a pseudo-time field as design variables to mimic the
fabrication sequence. The time field segments the given component, in a
differentiable manner, into (curved) layers. Numerical results show that the
fabrication sequence has a huge influence on the distortion. By optimizing the
sequence, the distortion is reduced by more than an order of magnitude in
comparison to commonly-used planar layers. The improved dimensional accuracy
is especially relevant for large metal components. The optimization framework
is not restricted to a certain distortion model, while in this paper it is
demonstrated using inherent strains (both isotropic and anisotropic) to
simulate distortion.

As future work, we intend to carry out an experimental investigation, first to
calibrate the thermal strain, and then to validate the effectiveness of the
optimized sequence on distortion minimization. The reduction in distortion is
likely accompanied by an increase in residual stresses. In parallel to
physical tests, we plan to integrate a restriction on thermal stresses during
the additive process into our sequence optimization framework. Stress is a
local property, and the stress field evolves during fabrication. This means a
very large number of constraints has to be integrated into the optimization.

The proposed computational framework opens up a new dimension in exploring the
full potential of additive manufacturing. The current consideration of
manufacturability is by no means comprehensive. Firstly, as can be observed in
e.g. Fig. 14, the layers have varying thicknesses. By varying the material
deposition rate and movement speed of the print head [44], the layer thickness
can be adjusted during fabrication, but only to a limited extent. An idea to
reduce the variation is by imposing minimum and maximum length scales on each
individual layer. Here the robust formulation [38] may serve as a good
starting point. Secondly, while the rotation of the build plate offers a
possibility to greatly reduce or even eliminate support structures [3], it
does not guarantee support-free fabrication for complex parts with complex
curved layers. The geometric condition for support-free fabrication needs to
be further investigated, and to be incorporated into fabrication sequence
optimization. Thirdly, the rotation of the build plate potentially gives rise
to interference between the print head and the already deposited part. This
issue is coupled with support-free process planning. A potential solution to
this is to restrict the orientation of curved layers in the optimization.

Acknowledgements The authors gratefully acknowledge the support from the
LEaDing Fellows Programme at the Delft University of Technology, which has
received funding from the European Union’s Horizon 2020 research and
innovation programme under the Marie Skłodowska-Curie grant agreement No.
707404. Weiming Wang wishes to thank the National Natural Science Foundation
of China (No. 62172073), and the Natural Science Foundation of Liaoning
Province (No. 2021-MS-110).

## References

  * Williams et al. [2016] S. W. Williams, F. Martina, A. C. Addison, J. Ding, G. Pardal, P. Colegrove,  Wire+ arc additive manufacturing,  Materials Science and Technology 32 (2016) 641–647. doi:10.1179/1743284715Y.0000000073.
  * Ding et al. [2011] J. Ding, P. Colegrove, J. Mehnen, S. Ganguly, P. S. Almeida, F. Wang, S. Williams,  Thermo-mechanical analysis of wire and arc additive layer manufacturing process on large multi-layer parts,  Computational Materials Science 50 (2011) 3315–3322. doi:10.1016/j.commatsci.2011.06.023.
  * Dai et al. [2018] C. Dai, C. C. Wang, C. Wu, S. Lefebvre, G. Fang, Y.-J. Liu,  Support-free volume printing by multi-axis motion,  ACM Transactions on Graphics 37 (2018) 1–14. doi:10.1145/3197517.3201342.
  * Fang et al. [2020] G. Fang, T. Zhang, S. Zhong, X. Chen, Z. Zhong, C. C. L. Wang,  Reinforced fdm: Multi-axis filament alignment with controlled anisotropic strength,  ACM Trans. Graph. 39 (2020). doi:10.1145/3414685.3417834.
  * Wang et al. [2020] W. Wang, D. Munro, C. C. Wang, F. van Keulen, J. Wu,  Space-time topology optimization for additive manufacturing,  Structural and Multidisciplinary Optimization 61 (2020) 1–18. doi:10.1007/s00158-019-02420-6.
  * Jafari et al. [2021] D. Jafari, T. H. Vaneker, I. Gibson,  Wire and arc additive manufacturing: Opportunities and challenges to control the quality and accuracy of manufactured parts,  Materials & Design (2021) 109471\. doi:10.1016/j.matdes.2021.109471.
  * Zalameda et al. [2013] J. N. Zalameda, E. R. Burke, R. A. Hafley, K. M. Taminger, C. S. Domack, A. Brewer, R. E. Martin,  Thermal imaging for assessment of electron-beam freeform fabrication (ebf3) additive manufacturing deposits,  in: Thermosense: Thermal Infrared Applications XXXV, volume 8705, International Society for Optics and Photonics, 2013, p. 87050M. doi:10.1117/12.2018233.
  * Mueller and Mueller [2000] D. Mueller, H. Mueller,  Experiences using rapid prototyping techniques to manufacture sheet metal forming tools,  in: Proc. ISATA Conference, Dublin, volume 9, 2000. doi:10.1.1.122.2370.
  * Ya and Hamilton [2017] W. Ya, K. Hamilton,  On-demand spare parts for the marine industry with directed energy deposition: propeller use case,  in: International Conference on Additive Manufacturing in Products and Applications, Springer, 2017, pp. 70–81. doi:10.1007/978-3-319-66866-6_7.
  * Nycz et al. [2017] A. Nycz, M. W. Noakes, B. Richardson, A. Messing, B. Post, J. Paul, J. Flamm, L. Love,  Challenges in making complex metal large-scale parts for additive manufacturing: A case study based on the additive manufacturing excavator,  in: 2017 International Solid Freeform Fabrication Symposium, University of Texas at Austin, 2017\. doi:10.26153/16923.
  * Nycz et al. [2021] A. Nycz, Y. Lee, M. Noakes, D. Ankit, C. Masuo, S. Simunovic, J. Bunn, L. Love, V. Oancea, A. Payzant, C. M. Fancher,  Effective residual stress prediction validated with neutron diffraction method for metal large-scale additive manufacturing,  Materials & Design 205 (2021) 109751. doi:10.1016/j.matdes.2021.109751.
  * Romero-Hdz et al. [2016] J. Romero-Hdz, B. Saha, G. Toledo-Ramirez, D. Beltran-Bqz,  Welding sequence optimization using artificial intelligence techniques: an overview,  Int J Comput Sci Eng 3 (2016) 90–95. doi:10.14445/23488387/IJCSE-V3I11P115.
  * Beik et al. [2019] V. Beik, H. Marzbani, R. Jazar,  Welding sequence optimisation in the automotive industry: A review,  Proceedings of the Institution of Mechanical Engineers, Part C: Journal of Mechanical Engineering Science 233 (2019) 5945–5952. doi:10.1177/0954406219859909.
  * Sattari-Far and Javadi [2008] I. Sattari-Far, Y. Javadi,  Influence of welding sequence on welding distortions in pipes,  International Journal of Pressure Vessels and Piping 85 (2008) 265–274. doi:10.1016/j.ijpvp.2007.07.003.
  * Ding et al. [2015] D. Ding, Z. Pan, D. Cuiuri, H. Li,  Wire-feed additive manufacturing of metal components: technologies, developments and future interests,  The International Journal of Advanced Manufacturing Technology 81 (2015) 465–481. doi:10.1007/s00170-015-7077-3.
  * Mughal et al. [2005] M. Mughal, H. Fawad, R. Mufti, M. Siddique,  Deformation modelling in layered manufacturing of metallic parts using gas metal arc welding: effect of process parameters,  Modelling and Simulation in Materials Science and Engineering 13 (2005) 1187\. doi:10.1088/0965-0393/13/7/013.
  * Mughal et al. [2007] M. Mughal, R. Mufti, H. Fawad,  The mechanical effects of deposition patterns in welding-based layered manufacturing,  Proceedings of the Institution of Mechanical Engineers, Part B: Journal of Engineering Manufacture 221 (2007) 1499–1509. doi:10.1243/09544054JEM783.
  * Graf et al. [2018] M. Graf, A. Hälsig, K. Höfer, B. Awiszus, P. Mayr,  Thermo-mechanical modelling of wire-arc additive manufacturing (WAAM) of semi-finished products,  Metals 8 (2018) 1009\. doi:10.3390/met8121009.
  * Yang et al. [2016] Q. Yang, P. Zhang, L. Cheng, Z. Min, M. Chyu, A. C. To,  Finite element modeling and validation of thermomechanical behavior of ti-6al-4v in directed energy deposition additive manufacturing,  Additive Manufacturing 12 (2016) 169–177. doi:10.1016/j.addma.2016.06.012.
  * Körner et al. [2020] C. Körner, M. Markl, J. A. Koepf,  Modeling and simulation of microstructure evolution for additive manufacturing of metals: a critical review,  Metallurgical and Materials Transactions A 51 (2020) 4970–4983. doi:10.1007/s11661-020-05946-3.
  * Wei et al. [2021] H. Wei, T. Mukherjee, W. Zhang, J. Zuback, G. Knapp, A. De, T. DebRoy,  Mechanistic models for additive manufacturing of metallic components,  Progress in Materials Science 116 (2021) 100703. doi:10.1016/j.pmatsci.2020.100703.
  * Singh et al. [2022] P. Singh, R. M.Singari, R. Mishra,  A review of study on modeling and simulation of additive manufacturing processes,  Materials Today: Proceedings 56 (2022) 3594–3603. doi:10.1016/j.matpr.2021.12.057, first International Conference on Design and Materials.
  * Denlinger and Michaleris [2017] E. R. Denlinger, P. Michaleris,  Mitigation of distortion in large additive manufacturing parts,  Proceedings of the Institution of Mechanical Engineers, Part B: Journal of Engineering Manufacture 231 (2017) 983–993. doi:10.1177/0954405415578580.
  * Denlinger and Michaleris [2016] E. R. Denlinger, P. Michaleris,  Effect of stress relaxation on distortion in additive manufacturing process modeling,  Additive Manufacturing 12 (2016) 51–59. doi:10.1016/j.addma.2016.06.011.
  * Salonitis et al. [2016] K. Salonitis, L. D’Alvise, B. Schoinochoritis, D. Chantzis,  Additive manufacturing and post-processing simulation: laser cladding followed by high speed machining,  The International Journal of Advanced Manufacturing Technology 85 (2016) 2401–2411. doi:10.1007/s00170-015-7989-y.
  * Chiumenti et al. [2010] M. Chiumenti, M. Cervera, A. Salmi, C. A. De Saracibar, N. Dialami, K. Matsui,  Finite element modeling of multi-pass welding and shaped metal deposition processes,  Computer methods in applied mechanics and engineering 199 (2010) 2343–2359. doi:10.1016/j.cma.2010.02.018.
  * Ueda et al. [1975] Y. Ueda, K. Fukuda, K. Nakacho, S. Endo,  A new measuring method of residual stresses with the aid of finite element method and reliability of estimated values,  Journal of the Society of Naval Architects of Japan 1975 (1975) 499–507. doi:10.2534/jjasnaoe1968.1975.138_499.
  * Murakawa et al. [1996] H. Murakawa, Y. Luo, Y. Ueda,  Prediction of welding deformation and residual stress by elastic fem based on inherent strain,  Journal of the society of Naval Architects of Japan 1996 (1996) 739–751. doi:10.2534/jjasnaoe1968.1996.180_739.
  * Ma et al. [2016] N. Ma, K. Nakacho, T. Ohta, N. Ogawa, A. Maekawa, H. Huang, H. Murakawa,  Inherent strain method for residual stress measurement and welding distortion prediction,  in: International Conference on Offshore Mechanics and Arctic Engineering, volume 50008, American Society of Mechanical Engineers, 2016, p. V009T13A001. doi:10.1115/OMAE2016-54184.
  * Munro et al. [2019] D. Munro, C. Ayas, M. Langelaar, F. van Keulen,  On process-step parallel computability and linear superposition of mechanical responses in additive manufacturing process simulation,  Additive Manufacturing 28 (2019) 738–749. doi:10.1016/j.addma.2019.06.023.
  * Liang et al. [2019] X. Liang, Q. Chen, L. Cheng, D. Hayduke, A. C. To,  Modified inherent strain method for efficient prediction of residual deformation in direct metal laser sintered components,  Computational Mechanics 64 (2019) 1719–1733. doi:10.1007/s00466-019-01748-6.
  * Chen et al. [2019] Q. Chen, X. Liang, D. Hayduke, J. Liu, L. Cheng, J. Oskin, R. Whitmore, A. C. To,  An inherent strain based multiscale modeling framework for simulating part-scale residual deformation for direct metal laser sintering,  Additive Manufacturing 28 (2019) 406–418. doi:10.1016/j.addma.2019.05.021.
  * Prabhune and Suresh [2020] B. C. Prabhune, K. Suresh,  A fast matrix-free elasto-plastic solver for predicting residual stresses in additive manufacturing,  Computer-Aided Design 123 (2020) 102829. doi:10.1016/j.cad.2020.102829.
  * Allaire and Jakabčin [2018] 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 (2018) 2313–2366. doi:10.1142/S0218202518500501.
  * Pellens et al. [2020] J. Pellens, G. Lombaert, M. Michiels, T. Craeghs, M. Schevenels,  Topology optimization of support structure layout in metal-based additive manufacturing accounting for thermal deformations,  Structural and Multidisciplinary Optimization 61 (2020) 2291–2303. doi:10.1007/s00158-020-02512-8.
  * Misiun et al. [2021] G. Misiun, E. van de Ven, M. Langelaar, H. Geijselaers, F. van Keulen, T. van den Boogaard, C. Ayas,  Topology optimization for additive manufacturing with distortion constraints,  Computer Methods in Applied Mechanics and Engineering 386 (2021) 114095\. doi:10.1016/j.cma.2021.114095.
  * Miki and Yamada [2021] T. Miki, T. Yamada,  Topology optimization considering the distortion in additive manufacturing,  Finite Elements in Analysis and Design 193 (2021) 103558. doi:10.1016/j.finel.2021.103558.
  * Wang et al. [2010] F. Wang, B. S. Lazarov, O. Sigmund,  On projection methods, convergence and robust formulations in topology optimization,  Struct. Multidiscip. Optim. 43 (2010) 767–784. doi:10.1007/s00158-010-0602-y.
  * Wu et al. [2018] J. Wu, N. Aage, R. Westermann, O. Sigmund,  Infill optimization for additive manufacturing – approaching bone-like porous structures,  IEEE Transactions on Visualization and Computer Graphics 24 (2018) 1127–1140. doi:10.1109/TVCG.2017.2655523.
  * Setien et al. [2019] I. Setien, M. Chiumenti, S. van der Veen, M. San Sebastian, F. Garciandía, A. Echeverría,  Empirical methodology to determine inherent strains in additive manufacturing,  Computers & Mathematics with Applications 78 (2019) 2282–2295. doi:10.1016/j.camwa.2018.05.015.
  * Liang et al. [2018] X. Liang, L. Cheng, Q. Chen, Q. Yang, A. C. To,  A modified method for estimating inherent strains from detailed process simulation for fast residual distortion prediction of single-walled structures fabricated by directed energy deposition,  Additive Manufacturing 23 (2018) 471–486. doi:10.1016/j.addma.2018.08.029.
  * Michaleris [2014] P. Michaleris,  Modeling metal deposition in heat transfer analyses of additive manufacturing processes,  Finite Elements in Analysis and Design 86 (2014) 51–60. doi:10.1016/j.finel.2014.04.003.
  * Svanberg [1987] K. Svanberg,  The method of moving asymptotes—a new method for structural optimization,  International journal for numerical methods in engineering 24 (1987) 359–373. doi:10.1002/nme.1620240207.
  * Kuipers et al. [2020] T. Kuipers, E. L. Doubrovski, J. Wu, C. C. Wang,  A framework for adaptive width control of dense contour-parallel toolpaths in fused deposition modeling,  Computer-Aided Design 128 (2020) 102907. doi:10.1016/j.cad.2020.102907.