# Regularization in Space-Time Topology Optimization for
Multi-Axis Additive Manufacturing

Weiming Wang<EMAIL_ADDRESS>Kai Wu<EMAIL_ADDRESS>Fred van Keulen
<EMAIL_ADDRESS>Jun Wu<EMAIL_ADDRESS>Department of Mechanical,
Aerospace and Civil Engineering, The University of Manchester, Manchester,
United Kingdom  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

In additive manufacturing, the fabrication sequence has a large influence on
the quality of manufactured components. While planning of the fabrication
sequence is typically performed after the component has been designed, recent
developments have demonstrated the possibility and benefits of simultaneous
optimization of both the structural layout and the corresponding fabrication
sequence. This is particularly relevant in multi-axis additive manufacturing,
where rotational motion offers enhanced flexibility compared to planar
fabrication. The simultaneous optimization approach, called space-time
topology optimization, introduces a pseudo-time field to encode the
manufacturing process order, alongside a pseudo-density field representing the
structural layout. To comply with manufacturing principles, the pseudo-time
field needs to be monotonic, i.e., free of local minima. However, explicitly
formulated constraints proposed in prior work are not always effective,
particularly for complex structural layouts that commonly result from topology
optimization.

In this paper, we introduce a novel method to regularize the pseudo-time field
in space-time topology optimization. We conceptualize the monotonic additive
manufacturing process as a virtual heat conduction process starting from the
surface upon which a component is constructed layer by layer. The virtual
temperature field, which shall not be confused with the actual temperature
field during manufacturing, serves as an analogy for encoding the fabrication
sequence. In this new formulation, we use local virtual heat conductivity
coefficients as optimization variables to steer the temperature field and,
consequently, the fabrication sequence. The virtual temperature field is
inherently free of local minima due to the physics it resembles. We
numerically validate the effectiveness of this regularization in space-time
topology optimization under process-dependent loads, including gravity and
thermomechanical loads.

###### keywords:

Space-time topology optimization , fabrication sequence optimization , multi-
axis additive manufacturing , heat equations

##  1 Introduction

The past few years have seen tremendous progress in computational design for
additive manufacturing [1] and topology optimization approaches in particular
[2, 3]. In topology optimization, structural design is formulated as an
optimization problem of distributing a limited amount of material in a design
domain to maximize its performance, for example, load-bearing capacity. The
highly optimized designs often have complex geometry, which makes them
difficult to produce. While the flexibility of additive manufacturing makes it
possible to produce geometrically complex components, it is not without
challenges. In the past decade, a focus in topology optimization has been on
the incorporation of additive manufacturability, to avoid optimized designs
that cannot be 3D printed directly or need costly post-processing. The
additive manufacturability is commonly characterized by geometric features
such as the overhang angle (e.g., [4, 5, 6, 7]) and enclosed voids (e.g. [8,
9]). Let us refer to the review article by Liu et al. [10] for developments up
to 2018.

A recent development in topology optimization for additive manufacturing has
been on the integration of the physics involved in the layer-by-layer additive
process into design optimization [11]. This means simulating the physics, such
as self-weight and thermomechanical loads, acting on intermediate structures
during the fabrication process. Allaire et al. [12] proposed one of the first
approaches in this emerging direction, addressing overhang constraints, which
had typically been treated as a geometric criterion [4, 5, 6, 7]. By
restricting the deformation of intermediate structures due to their self-
weight, the final optimized structure exhibits fewer overhanging features,
reducing the need for auxiliary supports during fabrication. The self-weight
of intermediate structures was considered similarly by Amir and Mass [13],
Bruggi et al. [14], and Haveroth et al. [15]. Starting from self-weight, more
sophisticated and computationally demanding physics in the fabrication process
has been incorporated in topology optimization, particularly for mitigating
undesired effects due to thermomechanical loads [16, 17, 18, 19]. These
process-informed topology optimization approaches, by incorporating a model of
the layer-by-layer manufacturing process in design optimization, provide a
more detailed level of quality control compared to focusing on geometric (and
topological) features of optimized designs.

In the aforementioned topology optimization approaches, the structural layout
is optimized to generate high-performance components. In addition to
structural design, the fabrication sequence, according to which the component
is produced line-by-line and layer-by-layer, significantly influences
manufacturability and the quality of fabricated components. In wire and arc
additive manufacturing (WAAM), for example, substantial distortion may occur
during fabrication as well as after cutting off the build plate [20]. The
distortion is largely influenced by the fabrication sequence. Typically, the
fabrication sequence is planned after the structural layout has been designed.
This sequential workflow results in a gap between the quality of digital
designs and their physical counterparts. Earlier work tackled this problem by
assuming a prescribed fabrication sequence, e.g., considering factors like
material anisotropy due to a prescribed toolpath strategy in the optimization
of the structural layout [21, 22]. More recent efforts have started optimizing
the toolpath for powder bed fusion [23, 24] and layers for multi-axis additive
manufacturing [25, 26]. Boissier et al. [23, 24] represented the toolpath as a
continuous curve, formulating it as a shape optimization problem.

Our work is motivated by emerging robot-assisted multi-axis additive
manufacturing (e.g., WAAM), where the rotation of robotic arms offers more
flexibility than traditional 3D printing, allowing not only translational but
also rotational motion. With rotation, the fabrication sequence is not limited
to planar layers. This brings in benefits such as turning an overhanging
feature into a self-supporting one, by appropriately rotating the build plate
[27]. In our previous work [25, 26], the fabrication sequence was encoded by a
pseudo-time field, where a larger pseudo-time indicates that the corresponding
location is to be materialized later in the fabrication process. Isocontours
of the pseudo-time field segment the structural layout into (curved) layers.
This parameterization is akin to density-based topology optimization in the
sense that each location in the design domain carries an optimization
variable, opening up the solution space of fabrication sequences. This
parameterization also enables the simultaneous optimization of the structural
layout and the fabrication sequence, i.e., by optimizing the density and time
field concurrently, referred to as space-time topology optimization [25].

Parameterizing the fabrication process by a pseudo-time field opens up a huge
solution space. However, the sequence represented by the isocontours of the
pseudo-time field is not guaranteed to be manufacturable. It is thus necessary
to regularize the pseudo-time field to ensure that the optimized fabrication
sequence meets manufacturing requirements. A critical requirement that we
focus on in this paper is the continuity of fabrication: material can only be
deposited on the already produced portion of the part, besides directly on the
build plate. In the pseudo-time field, the violation of this requirement can
be characterized by local minima, i.e., the surrounding points are associated
with larger time values, and thus have not been materialized yet. A local
minimum may present inside the domain or on the boundary of the domain. Both
are undesirable, except it is adjacent to the build plate. An illustration of
these violations is illustrated in Fig. 1.

Figure 1: (a) Illustration of a component and a pseudo-time field, with a
smaller value indicating that the corresponding point is materialized earlier
in fabrication. (b) The isolines of the pseudo-time field segment the part
into a set of consecutive layers. (c) A series of intermediate structural
layouts during the fabrication process. Note that, in stages 1 to 3, a patch
of material, corresponding to a local minimum in the pseudo-time field, is in
isolation from the rest of the component. This renders the fabrication
sequence infeasible.

In our prior work [25, 26], we addressed these violations by introducing
explicitly formulated constraints. While these constraints function well in
most cases, their effectiveness depends on the initialization of the pseudo-
time field. In this paper, we present a novel regularization method to
implicitly handle these constraints. Our starting point is to initiate the
pseudo-time field using geodesic distances. Unlike Euclidean distances,
geodesic distances are dependent on the structural layout and free of local
minima. This initialization strategy works for fixed structural layouts, but
is not directly applicable to space-time optimization, where the structural
layout and the fabrication sequence are optimized concurrently. Because the
structural layout evolves over iterations, and is not well defined during
optimization, i.e., grey elements as encountered in density-based topology
optimization. To solve these problems, we reformulate this initialization
strategy using a partial differential equation, a heat equation in particular.
This equation enables us to encode intermediate densities and evolving
structural layouts in its parameters.

The remainder of this paper is organized as follows. In Section 2, we present
a simplified problem setup of space-time optimization where the structural
layout is fixed. This problem is solved by initialiating the fabrication
sequence based on geodesic distances. In Section 3, this solution to the
simplified problem setup is extended to address problems in space-time
optimization where the structural layout and the fabrication sequence are
optimized concurrently. The effectiveness of the proposed approach is
evaluated in Section 4 for design optimization in 2D and 3D under process-
dependent self-weight and thermomechanical loads. In Section 5, we conclude
the paper with some ideas for future work.

##  2 Fabrication sequence optimization

In this section, we consider a simplified version of space-time topology
optimization. Here, the structural layout is prescribed, and we optimize the
fabrication sequence to produce it. This reduced problem is intended to
illustrate the problem of local minima, demonstrate the effects of the initial
sequence in optimization, and highlight the efficacy of a novel initialization
strategy. These insights inspire our approach to solving the problem in the
full space-time topology optimization.

###  2.1 Sequence parameterization

We parameterize the fabrication sequence to produce a component using a
pseudo-time field defined over the domain of the component. The component is
discretized using a Cartesian grid, resulting in quadrilateral elements in 2D
and hexahedral elements in 3D. To mimic the additive manufacturing process,
each element $e$ is associated with a scalar value $t_{e}$, indicating the
pseudo-time at which the material associated with the element is added to the
structure. Thus, an element with a larger pseudo-time indicates that this
element is materialized later than an element with a smaller one. For
simplicity, we also refer to the pseudo-time field as the time field.

Figure 2: Schematic of curved layers for fabricating an $L$-shaped component.
Two successive intermediate structures during fabrication are illustrated on
the left and right.

The time field $\bm{t}$ serves as design variables in sequence optimization.
It is bounded to be between $0$ and $1$. The isocontours by
$\eta_{j}=\frac{j}{N}$, with $j=0,...,N$, thus partition the component into
$N$ layers. The partition is illustrated in Fig. 2. Each layer consists of
elements whose time value belongs to the corresponding time interval. For
instance, the $j$-th layer is denoted as
$L_{j}=\\{e|t_{e}\in[\eta_{j-1},\eta_{j})\\},j=1,2,...,N$. We note that the
time field is a continuous scalar field, mimicking the progressive additive
manufacturing process. The isocontours of the time field, and the resulting
layers between two consecutive isocontours, are not restricted to be planar
but can be curved.

We describe the intermediate structures during fabrication by a series of
pseudo-density fields (or density fields, for simplicity),
$\bm{\rho}^{\\{j\\}}$, with $j=0,1,...,N$. Two successive intermediate
structures are illustrated in Fig. 2. The superscript $\\{j\\}$ indicates the
$j$-th stage during fabrication, i.e., layers from the first up to the $j$-th
have been fabricated. $j=0$ refers to the start of fabrication,
$\rho^{\\{0\\}}_{e}=0$, $\forall e$, while $j=N$ refers to the completion of
fabrication, $\rho^{\\{N\\}}_{e}=1$, $\forall e$.

The evolving intermediate structure is derived from the time field. Let
$\bm{1}$ denote the density field of a prescribed, solid component. Assuming
that the density field and time field are discretized using the same Cartesian
grid, the intermediate structure at the $j$-th stage is described by an entry-
wise multiplication between the density and the time field,

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

Here $\bar{\bm{t}}^{\\{\eta\\}}$ is a function of the time field $\bm{t}$. It
converts time values less than $\eta$ to $1$ and time values greater than
$\eta$ to $0$, in a differentiable manner. It is implemented using a smoothed
Heaviside function, a commonly employed approach in density-based topology
optimization [28, 29],

$\bar{{t}}^{\\{\eta\\}}=1-\frac{\tanh(\beta_{t}\eta)+\tanh(\beta_{t}({t}-\eta))}{\tanh(\beta_{t}\eta)+\tanh(\beta_{t}({1}-\eta))},$
(2)

where $\beta_{t}$ is a positive number to control the projection sharpness. A
continuation scheme is applied to increase this sharpness parameter after a
certain number of iterations.

###  2.2 Sequence optimization

The sequence optimization is generically formulated as

$\displaystyle\underset{\bm{t}}{\text{minimize}:}$ $\displaystyle
c(\bm{t})=f(\bm{\rho}^{\\{1\\}}(\bm{t}),...,\bm{\rho}^{\\{j\\}}(\bm{t}),...,\bm{\rho}^{\\{N\\}}(\bm{t})),$
(3) $\displaystyle\text{subject to}:$ $\displaystyle
t_{e}\geq\min\\{t_{i}|i\in\mathcal{N}_{e}\\},\forall
e\in\mathcal{M}\backslash\mathcal{S}_{0},$ (6) $\displaystyle
t_{e}\leq\max\\{t_{i}|i\in\mathcal{N}_{e}\\},\forall
e\in\mathcal{M}\backslash\mathcal{B},$ $\displaystyle t_{e}\in[0,1],\forall
e.$

The objective concerns the evolving structure during fabrication
$\bm{\rho}^{\\{j\\}}$, which is a function of the time field (cf. Eq. 1). In
Eq. 6 and Eq. 6, $\mathcal{N}_{e}$ denotes the set of elements adjacent to
element $e$. Eq. 6 requests the time value of each element to be larger than
the minimum among its adjacent elements. This requirement applies to all
elements in the component ($\mathcal{M}$), except those designated as the
start region for fabrication ($\mathcal{S}_{0}=\\{e|t_{e}=0\\}$), i.e., on the
build plate. An element where this condition is violated indicates that the
element needs to be materialized in midair, i.e., before any of its adjacent
elements have been produced. Similarly, Eq. 6 requests the time value of each
element to be smaller than the maximum among its adjacent elements. This
requirement applies to all elements in the component ($\mathcal{M}$), except
those on the boundary of the component ($\mathcal{B}$). An element where this
condition is violated indicates that the element needs to be materialized
within an enclosed region, i.e., all of its adjacent elements have been
produced. Both Eq. 6 and Eq. 6 are essential for the validity of a fabrication
sequence.

The challenges to solving this formulation are twofold. Firstly, the
constraints (Eq. 6 and Eq. 6) involve $\min$ and $\max$ functions, which are
non-differential. Secondly, the number of constraints is large, since both Eq.
6 and Eq. 6 apply to all elements in the component. In previous work [25, 26],
these constraints were reformulated using approximation and aggregation
schemes. For instance,

$\frac{1}{{n}(\mathcal{M})}\sum_{e\in\mathcal{M}}||t_{e}-\textrm{mean}_{i\in\mathcal{N}_{e}}(t_{i})||^{2}\leq\epsilon,$
(7)

where $n$ denotes the number of elements in a set. The mean function
calculates the average of time values of adjacent elements, i.e.,
$\textrm{mean}_{i\in\mathcal{N}_{e}}(t_{i})=\frac{1}{{n}(\mathcal{N}_{e})}\sum_{i\in\mathcal{N}_{e}}t_{i}$.
$\epsilon$ is a small constant for numerical stability (e.g.,
$\epsilon=10^{-6}$). While this approximation works well in most cases, its
efficacy depends on the initialization of the fabrication sequence.

Figure 3: Comparison of the Euclidean and the geodesic distance field as the
initialization of the time field for fabrication sequence optimization. (a) A
2D component. The Euclidean distance and the geodesic distance from location
$p$ to the build plate at the bottom are shown by the green line and the black
polyline, respectively. (b) The Euclidean distance field. (c) The geodesic
distance field. (d) The optimized time field using (b) as initialization. (e)
The optimized time field using (c) as initialization. The sequences
corresponding to the time fields in (d) and (e) are shown in (f) and (g),
respectively.

###  2.3 Initialization of fabrication sequence

The initialization of the fabrication sequence greatly influences the validity
of optimized sequences, especially when the structural layout is complex.
Distance fields serve as a good option for initialization as they can mimic
the progression of additive manufacturing. A distance field is a scalar field,
representing the minimal distance from each point in the domain to a certain
shape, i.e., the build plate. Depending on the distance metric, the distance
field has quite different features. Fig. 3 compares the effects of two
distance metrics as initialization for sequence optimization. The component is
illustrated in Fig. 3 (a), where a horizontal build plate is located at the
bottom111In multi-axis additive manufacturing, the build plate is not
necessarily flat nor located at the bottom.. The two distance metrics under
comparison are the Euclidean distance and geodesic distance. The Euclidean
distance measures the shortest path (i.e., a straight line) between a point
and the build plate. The Euclidean distance is irrespective of the shape of
the component. In contrast, the geodesic distance measures the shortest path,
within the component, from a point to the build plate. The shortest Euclidean
path and geodesic path for a single point are illustrated in Fig. 3 (a). The
corresponding distance fields, together with isocontours, are visualized in
Fig. 3 (b) and (c). Both distance fields are normalized by dividing the
distance values over the corresponding largest distance. In Fig. 3 (b), the
isocontours of the Euclidean distance field are planar. On the right-hand side
of the domain, there is an overhanging feature, and there the isocontours are
separated by the void region. In contrast, in Fig. 3 (c), the geodesic
distance field is free of local minima. The isocontours of the geodesic
distance field in the top-right of the component are curved. The optimized
fabrication sequences starting from these two initial sequences are shown in
Fig. 3 (d) and (e), while the corresponding layered fabrication processes are
shown in (f) and (g), respectively. It can be seen in Fig. 3 (f) that during
fabrication, an isolated layer appears at stage 6. This shows that the
optimization is incapable of eliminating the local minimum inherited from the
initial field. While this violates the constraint in Eq. 6, it numerically
satisfies the approximated version (Eq. 7). Setting a smaller threshold
($\epsilon$) proves to be ineffective. Even though the approximated constraint
is not precisely satisfied, the optimization is incapable of steering the
location of the local minimum such that it disappears. This experiment shows
the importance of a valid initial sequence that can mimic the progression of
additive manufacturing and is free of local minima. The geodesic distance
field fulfills these requirements.

##  3 Regularization of space-time topology optimization with heat equations

As discussed in the previous section, using the geodesic distance field to
initialize the time field is promising for optimizing the fabrication sequence
when the structural layout is fixed. However, it is not directly applicable to
space-time topology optimization, where the structural layout, including its
topology, evolves during the optimization process.

In this section, we present a novel approach to regularize the time field in
order to avoid the occurrence of local minima in the context of space-time
topology optimization. Specifically, we relate the time field to the solution
of a heat equation, i.e., a temperature field. The temperature field is
monotonic and free of local minima. Furthermore, the heat equation is defined
on the structural layout, and the resulting temperature field is sensible to
changes in the structural layout. The monotonicity and the sensibility to
shape changes are essential features of the geodesic distance field.
Furthermore, by adapting the parameters in the heat equation we can
accommodate intermediate density values as encountered during optimization.
Our adoption of the heat equation is inspired by advancements in computer
graphics, where research has demonstrated that the heat equation offers a
computationally efficient means for calculating an approximate geodesic
distance field [30]. In topology optimization for additive manufacturing, heat
equations were previously used for identifying enclosed voids and consequently
avoiding them [31, 32].

###  3.1 Optimization workflow

The workflow of space-time topology optimization using heat equations for
regularizing the time field is illustrated in Fig. 4. We employ a density
field ($\bm{\rho}$, $\rho_{e}\in[0,1]$) to represent the complete structure
and a time field ($\bm{t}$, $t_{e}\in[0,1]$) to encode the fabrication
sequence. While the density field eventually converges to either 0 or 1, the
time field is continuous. The product of $\bm{\rho}$ and $\bm{t}$ is projected
into a series of intermediate structures, using increasing threshold values
$\eta_{j}=\frac{j}{N}$, with $j=0,...,N$ and $N$ being the prescribed number
of layers. Rather than taking $\bm{\rho}$ and $\bm{t}$ as design variables, we
introduce two auxiliary fields, $\bm{\phi}$ and $\bm{\mu}$. The density field
$\bm{\rho}$ is computed from design field $\bm{\phi}$, following a density
filter ($\tilde{\bm{\phi}}$) and a smoothed Heaviside projection
($\bar{\tilde{\bm{\phi}}}$). This process
($\bm{\phi}\to\tilde{\bm{\phi}}\to\bm{\rho}=\bar{\tilde{\bm{\phi}}}$) is
commonly used in density-based topology optimization, and thus not included in
the workflow illustration.

The novel element in the workflow is the transformation from $\bm{\mu}$ and
$\bm{\rho}$, through ${\bm{\tau}}$, to $\bm{t}$. We use a heat equation to
obtain a monotonic time field defined on an evolving structural layout
$\bm{\rho}$ with heterogeneous thermal diffusivity $\bm{\mu}$. The solution of
the heat equation is a virtual temperature field $\bm{\tau}$. It is important
to note that this virtual temperature field shall not be confused with the
temperature field coming from the thermal process involved in metal additive
manufacturing. Here the heat equation is meant to create a monotonic field,
rather than predicting the actual temperature distribution during fabrication.
Conceptually, it works as follows. An artificial heat source is placed on the
build plate (at the bottom of the domain), from where the heat diffuses
through the structure. This creates a continuous temperature distribution
within the structure. The heat source has a constant temperature of $1$. At a
location in the structure further away from the heat source, its temperature
approaches $0$. The time field, ascending from the build plate, is defined as
$\bm{t}=\bm{1}-\bm{\tau}$.

Figure 4: The workflow of deriving a series of intermediate structures
$\bm{\rho}^{\\{j\\}}$ from the thermal diffusivity field $\bm{\mu}$ and
density field $\bm{\rho}$.

###  3.2 Heat equation

To obtain a monotonic temperature field, we construct the following heat
equation,

$\nabla\cdot(\kappa\nabla\tau)-\alpha_{\texttt{T}}\tau=0.$ (8)

Here $\nabla$ is the vector differential operator. This partial differential
equation represents steady-state heat conduction within a heterogeneous
medium, featuring spatially-varying thermal diffusivity $\kappa(x,y)$, and a
drain term $\alpha_{\texttt{T}}\tau$, where $\alpha_{\texttt{T}}$ denotes a
constant drain rate. The introduction of the drain term serves a crucial
purpose. In the absence of this term, achieving a steady-state temperature
distribution with spatial variation would necessitate employing an
inhomogeneous Dirichlet boundary condition. However, when dealing with a
structure that evolves during space-time topology optimization, there is no
clear rationale for determining the ideal location of a heat sink. An
alternative approach to achieve a spatially-varying temperature distribution
is to resort to transient heat conduction, and calculate the temperature
distribution at a specific time [30, 33].

The heat equation is subject to the following boundary conditions: At the
boundary adjacent to the build plate $\Gamma_{0}$, the temperature is
constant, ${\tau\big{|}}_{\Gamma_{0}}=1$. Other boundary segments of the
structure are thermally insulated,
${(\nabla\tau\cdot\overrightarrow{n})}\big{|}_{\Gamma\setminus\Gamma_{0}}=0$.
Here $\nabla\tau$ represents the spatial gradient of the temperature field,
and $\overrightarrow{n}$ denotes the normal direction on the boundary of the
structure.

The solution of the heat equation, and the corresponding time field together
with fabrication sequences, are illustrated in Fig. 5. In this example the
thermal diffusivity field ${\kappa}{(x,y)}$ is homogeneous. The temperature
field and, consequently, the fabrication sequence are determined by the
thermal diffusivity field. To allow the optimization of the fabrication
sequence, the thermal diffusivity field is related to design variables
$\bm{\mu}$, augmented by the unitless pseudo-density $\bm{\rho}$,

$\bm{\kappa}=\bm{\rho}\bm{\mu}.$ (9)

The multiplication by the pseudo-density field constrains heat conduction
within the structural layout, excluding void regions from the fabrication
sequence planning. Additionally, this allows pseudo-densities to take
intermediate values during the optimization process, thereby facilitating
gradient-based numerical optimization.

Figure 5: Heat conduction on a given component (a). The start points are the
nodes on the bottom boundary. The generated virtual temperature field by
solving Equation (10) is shown in (c). The light yellow indicates the largest
temperature, and the dark blue indicates the smallest virtual temperature. The
corresponding time field is shown in (c). The sequence generated with the time
field shown in (d) which is composed of 10 layers.

We note interesting similarities and large differences to the filters based on
Helmholtz‐type partial differential equations (PDEs) [34]. In our formulation,
the optimization variables (thermal diffusivity) in the PDE and the solution
of the PDE (virtual temperature) have different physical interpretations. The
desirable characteristic of the PDE is to attain a monotonic scalar field. In
contrast, in the PDE filter [34], the optimization variables and the solution
both represent density fields. The intention behind it is to get a smoothed
scalar field, avoiding checkerboard patterns.

####  3.2.1 Numerical discretization

We discretize Eq. 8 using the finite element method, resulting in a linear
form,

$\bm{K}_{\texttt{T}}\bm{T}=\bm{b},$ (10)

where $\bm{T}$ is the nodal temperature vector. $\bm{b}$ is the thermal load
vector corresponding to the boundary condition. $\bm{K}_{\texttt{T}}$ is the
global system matrix, assembled from element matrices:

$\bm{K}_{\texttt{T}}={\sum}_{e}{\bm{K}^{e}_{\texttt{T}}}+{\sum}_{e}{\bm{K}^{e}_{\texttt{C}}}.$
(11)

Here the symbol ${\sum}$ stands for matrix assembly. The first element matrix
corresponds to the thermal diffusivity. It is derived from the shape functions
$N_{i}$ and $N_{j}$ by using the variational principle:

$\bm{K}^{e}_{\texttt{T}}=\int_{\Omega_{e}}{\kappa_{e}\left(\frac{\partial
N_{i}}{\partial x}\frac{\partial N_{j}}{\partial x}+\frac{\partial
N_{i}}{\partial y}\frac{\partial N_{j}}{\partial y}\right)d\Omega}.$ (12)

The second element matrix, corresponding to the drain term, can be derived
similarly,

$\bm{K}^{e}_{\texttt{C}}=\int_{\Omega_{e}}{\alpha_{\texttt{T}}N_{i}N_{j}d\Omega}.$
(13)

In our implementation, instead of calculating a consistent matrix
$\bm{K}^{e}_{\texttt{C}}$ from the shape functions using Eq. 13, we adopt a
lumped matrix for simplicity, assuming that the drainage is lumped at the
nodes.

The nodal temperature distribution from solving Eq. 10 is transformed into an
elemental temperature field using bi-linear interpolation. In the matrix form,
this is written as $\bm{\tau}=\bm{G}\bm{T}$. $\bm{G}$ is a sparse
transformation matrix for converting nodal values to elemental values.
Afterwards, the time field is computed by $\bm{t}=\bm{1}-\bm{\tau}$.

####  3.2.2 Parameter analysis

The constant drain rate $\alpha_{\texttt{T}}$ is a key parameter in the heat
equation. Alongside the thermal diffusivity $\kappa$, this parameter plays a
significant role in determining the temperature distribution. To be able to
generate layers of comparable size, the temperature distribution shall be
smoothly varying between 0 and 1. Fig. 6 shows that a too small or large drain
rate results in layers with significant variations in size. The domain is
homogeneous in terms of thermal diffusivity $\kappa$, taking a value of $0.1$.
The heat source, with a constant temperature of $1$, is placed on the bottom
of the domain. If the drain rate is large, the temperature rapidly decreases
to 0, as shown in Fig. 6 (d). After the transformation
$\bm{t}=\bm{1}-\bm{\tau}$, and with uniformly sampled thresholds
$\eta_{j}=\frac{j}{N}$, the layers are squeezed to the bottom. Reversely, a
small drain rate also leads to uneven layers (Fig. 6 b). The segmentation in
Fig. 6 (c) is more uniform, providing a good starting point for the
optimization of the thermal diffusivity field towards the desired segmentation
of the domain into layers of equal size. To ensure numerical convergence, the
thermal diffusivity variable $\mu$ is constrained within the range of $0$ to
$1$, matching the range of the density variable $\phi$.

The temperature distribution is also affected by the size of the design
domain. We introduce a size-independent parameter to help choose an
appropriate drain rate. This parameter, denoted as $\beta$, relates the drain
rate ($\alpha_{\texttt{T}}$) to the characteristic length of the design domain
($l_{c}$),

$\alpha_{\texttt{T}}=\beta/l^{2}_{c}.$ (14)

The unit of the drain rate is $\texttt{m}^{-2}$, while the unit of the length
is m. Thus $\beta$ is non-dimensional. We choose the longest path of heat
diffusion in the design domain as its characteristic length.

In Fig. 6 (e) we compare the temperature distribution under five different
$\beta$ values. When $\beta=1$, the resulting temperature smoothly decreases
from 1 on the left, to close to 0 on the right. From numerical tests we
observed good convergence in the optimization of $\mu$ and $\phi$ with
$\beta\in[0.1,1]$. Slight derivation from this range does not lead to dramatic
convergence issues. It is robust for different domain sizes and volume
fractions. In the results section, $\beta=0.1$ is consistently used for all
examples.

Figure 6: The influence of the drain rate on the temperature distribution. (a)
A rectangular component with a homogeneous thermal diffusivity. The virtual
heat source with a fixed temperature of 1 is located at the bottom of the
domain. (b, c, d) Virtual temperature distributions with $\beta=0.1$, $1$, and
$10$, respectively. The black lines represent the isocontours of the
temperature fields at temperature values of 0.2, 0.4, 0.6, and 0.8, which are
intended to segment the component into 5 layers of comparable size. (e)
Temperature distributions along a vertical line in the middle of the
component, with five different $\beta$ values.

###  3.3 Formulation of space-time topology optimization

To summarize, the space-time topology optimization using the heat equation is
generically formulated as follows.

$\displaystyle\underset{{\bm{\phi},\;{\bm{\mu}}}}{\text{minimize}:}$
$\displaystyle
f=f_{0}(\bm{\rho})+\sum_{j=1}^{N}\alpha_{j}f_{j}(\bm{\rho}^{\\{j\\}}),$ (15)
$\displaystyle\text{subject to}:$
$\displaystyle\bm{K}(\bm{\rho})\bm{U}=\bm{F},$ (22)
$\displaystyle\bm{K}^{\\{j\\}}(\bm{\rho}^{\\{j\\}})\bm{U}^{\\{j\\}}=\bm{F}^{\\{j\\}},\quad
j=1,...,N,$
$\displaystyle\bm{K}_{\texttt{T}}(\bm{\phi},\bm{\mu})\bm{T}=\bm{b},$
$\displaystyle g_{0}(\bm{\rho})=\sum\rho_{e}v_{e}-V_{0}\leq 0,$ $\displaystyle
g_{j}(\bm{\rho}^{\\{j\\}})=\sum\rho_{e}^{\\{j\\}}v_{e}-\frac{j}{N}V_{0}\leq
0,\quad j=1,...,N,$ $\displaystyle 0\leq\phi_{e}\leq 1,\forall
e\in\mathcal{M},$ $\displaystyle 0\leq\mu_{e}\leq 1,\forall e\in\mathcal{M}.$

The objective function $f$ is composed of two parts. $f_{0}$ accounts for the
performance of the final overall structure, while $f_{j}$ is related to each
intermediate structure during the fabrication process. The different
objectives are balanced by the weighting factors $\alpha_{j}$. The overall
structure and intermediate structures all have their governing mechanical
equation, i.e., Eq. 22 and Eq. 22, respectively. These mechanical equations
are independent of each other, and thus can be solved in parallel. $\bm{K}$
and $\bm{U}$ are the global stiffness matrix and displacement vector
corresponding to the overall structure, respectively, and $\bm{F}$ is the
external force. The global stiffness matrix $\bm{K}^{\\{j\\}}$ of intermediate
structures is determined by both the density variables $\phi$ and the
diffusivity variables $\mu$. The load vector $\bm{F}^{\\{j\\}}$ may also be
dependent on these variables, e.g., process-dependent thermomechanical loads.
Eq. 22 is the heat equation for obtaining a monotonic temperature field, as
introduced in the previous subsection. The overall material usage is
constrained in Eq. 22, while an even amount of material for each intermediate
structure is described by Eq. 22 to ensure a constant fabrication speed.

Here, we present the general formulation of space-time topology optimization.
For a more in-depth exploration of objective and constraint functions in
specific application scenarios, such as structural compliance and thermal-
induced distortion during fabrication, readers can refer to [25, 26]. When
considering the entire density field as fixed, the optimization model
described above simplifies into the formulation outlined in Section 2.2,
focusing solely on sequence optimization. We solve the optimization problem
using the method of moving asymptotes (MMA) [35].

###  3.4 Sensitivity analysis

The general procedure for sensitivity calculation is consistent with those
outlined in [25, 26]. Here we focus on the primary distinction that arises
from the fact that the time field is no longer a direct design variable;
rather, it functions as an implicit variable dependent on both density $\phi$
and diffusivity $\mu$.

##### Sensitivity w.r.t design variable $\bm{\mu}$

Take the objective function $f_{j}$ as an example, the sensitivity is
calculated using the adjoint method by adding an augmented Lagrangian term:

$f_{j}=f_{j}+\bm{\lambda}^{\mathsf{T}}\left(\bm{K}_{\texttt{T}}\bm{T}-\bm{b}\right),$
(23)

where $\bm{\lambda}$ is the Lagrangian multiplier. Without loss of generality,
if $f_{j}$ has no explicit component of the diffusivity design variable
$\mu_{e}$, the derivative of $f_{j}$ with respect to $\mu_{e}$ is given by:

$\displaystyle\frac{\partial f_{j}}{\partial\mu_{e}}$
$\displaystyle=\sum_{i\in\mathcal{M}}\frac{\partial f_{j}}{\partial
t_{i}}\frac{\partial
t_{i}}{\partial\mu_{e}}+\bm{\lambda}^{\mathsf{T}}\left(\frac{\partial\bm{K}_{\texttt{T}}}{\partial\mu_{e}}\bm{T}+\bm{K}_{\texttt{T}}\frac{\partial\bm{T}}{\partial\mu_{e}}\right).$
(24)

Rewrite the above equation in vector form while using the inverse relation
between the time field and nodal temperature field
$\bm{t}=\bm{1}-\bm{G}{\bm{T}}$:

$\displaystyle\frac{\partial f_{j}}{\partial\mu_{e}}$
$\displaystyle=-\left(\frac{\partial
f_{j}}{\partial\bm{t}}\right)^{\mathsf{T}}\bm{G}\frac{\partial\bm{T}}{\partial\mu_{e}}+\bm{\lambda}^{\mathsf{T}}\left(\frac{\partial\bm{K}_{\texttt{T}}}{\partial\mu_{e}}\bm{T}+\bm{K}_{\texttt{T}}\frac{\partial\bm{T}}{\partial\mu_{e}}\right)$
(25) $\displaystyle=\left(-\left(\frac{\partial
f_{j}}{\partial\bm{t}}\right)^{\mathsf{T}}\bm{G}+\bm{\lambda}^{\mathsf{T}}\bm{K}_{\texttt{T}}\right)\frac{\partial\bm{T}}{\partial\mu_{e}}+\bm{\lambda}^{\mathsf{T}}\frac{\partial\bm{K}_{\texttt{T}}}{\partial\mu_{e}}\bm{T}.$

To directly calculate the derivative of the temperature field with respect to
the thermal diffusivity, the Lagrangian multiplier $\bm{\lambda}$ is chosen to
eliminate the first term:

$-\left(\frac{\partial
f_{j}}{\partial\bm{t}}\right)^{\mathsf{T}}\bm{G}+\bm{\lambda}^{\mathsf{T}}\bm{K}_{\texttt{T}}=0.$
(26)

The derivative of $f_{j}$ regarding to $\mu_{e}$ can thus be rewritten as

$\frac{\partial
f_{j}}{\partial\mu_{e}}=\bm{\lambda}^{\mathsf{T}}\frac{\partial\bm{K}_{\texttt{T}}}{\partial\mu_{e}}\bm{T}.$
(27)

##### Lagrangian multiplier $\bm{\lambda}$

In order to solve for the Lagrangian multiplier in Eq. 26, the derivative of
$f_{j}$ with respect to the time field $\bm{t}$ needs to be calculated. Since
$f_{j}$ is the function of $\bm{\rho}^{\\{j\\}}$, from the chain rule:

$\frac{\partial f_{j}}{\partial t_{e}}=\frac{\partial
f_{j}}{\partial\bar{t}_{e}}\frac{\partial\bar{t}_{e}}{\partial t_{e}}.$ (28)

where $\partialderivative*{\bar{t}_{e}}{t_{e}}$ is calculated from the
Heaviside projection of Eq. 2. In most of the cases, $f_{j}$ contains terms
related to $\bm{U}^{\\{j\\}}$, thus another adjoint analysis is needed:

$\displaystyle\frac{\partial f_{j}}{\partial\bar{t}_{e}}$
$\displaystyle=\frac{\partial
f_{j}}{\partial\bar{t}_{e}}+\bm{\gamma}^{\mathsf{T}}\partialderivative{\left(\bm{K}^{\\{j\\}}\bm{U}^{\\{j\\}}-\bm{F}^{\\{j\\}}\right)}{\bar{t}_{e}}$
(29) $\displaystyle=\frac{\partial
f_{j}}{\partial\bm{U}^{\\{j\\}}}\frac{\partial\bm{U}^{\\{j\\}}}{\partial\bar{t}_{e}}+\bm{\gamma}^{\mathsf{T}}\left(\frac{\partial\bm{K}^{\\{j\\}}}{\partial\bar{t}_{e}}\bm{U}^{\\{j\\}}+\bm{K}^{\\{j\\}}\frac{\partial\bm{U}^{\\{j\\}}}{\partial\bar{t}_{e}}\right)$
$\displaystyle=\left(\frac{\partial
f_{j}}{\partial\bm{U}^{\\{j\\}}}+\bm{\gamma}^{\mathsf{T}}\bm{K}^{\\{j\\}}\right)\frac{\partial\bm{U}^{\\{j\\}}}{\partial\bar{t}_{e}}+\bm{\gamma}^{\mathsf{T}}\frac{\partial\bm{K}^{\\{j\\}}}{\partial\bar{t}_{e}}\bm{U}^{\\{j\\}}.$

Here $\bm{\gamma}$ is the Lagrangian multiplier related to the mechanical
equilibrium equation of intermediate structure, which is solved from:

$\frac{\partial
f_{j}}{\partial\bm{U}^{\\{j\\}}}+\bm{\gamma}^{\mathsf{T}}\bm{K}^{\\{j\\}}=0.$
(30)

The equation above simply assume that there is no explicit term about
$\bar{t}_{e}$ in $f_{j}$, and the external loads $F^{\\{j\\}}$ is independent
on the design $\bm{\rho}^{\\{j\\}}$. These are not necessary true, for
example, design-dependent load like self-weight can also be considered by
making slight modifications on the sensitivity analysis.

##### Sensitivity w.r.t design variable $\bm{\phi}$

The derivative of $f_{j}$ with respect to the density design variable
$\phi_{e}$ can be obtained using the chain rule, which aligns with the
filtering and projection procedure:

$\frac{\partial
f_{j}}{\partial\phi_{e}}=\sum_{i\in\mathcal{N}_{e}}\derivative{f_{j}}{\rho_{i}}\frac{\partial\rho_{i}}{\partial\tilde{\phi}_{i}}\frac{\partial\tilde{\phi}_{i}}{\partial\phi_{e}}.$
(31)

The derivative of $f_{j}$ regarding the physical density $\rho_{e}$ is also
calculated using the adjoint method. However, it should be noted that since
the objective function has multiple components that contains $\rho_{e}$
implicitly. As indicated in Eq. 9, the time field obtained from the heat
equation is also dependent on the density field. We thus divide the
sensitivity of $f_{j}$ into two parts. The first part is independent of the
time field $\bm{t}$, while the second part is related to the time field
$\bm{t}$.

$\displaystyle\derivative{f_{j}}{\rho_{e}}$ $\displaystyle=\frac{\partial
f_{j}}{\partial\rho_{e}}+\sum_{i\in\mathcal{M}}\frac{\partial f_{j}}{\partial
t_{i}}\frac{\partial t_{i}}{\partial\rho_{e}}.$ (32)

Similar to Eq. 25, the above equation is written into an augmented form:

$\displaystyle\derivative{f_{j}}{\rho_{e}}$ $\displaystyle=\frac{\partial
f_{j}}{\partial\rho_{e}}-\left(\frac{\partial
f_{j}}{\partial\bm{t}}\right)^{\mathsf{T}}\bm{G}\frac{\partial\bm{T}}{\partial\rho_{e}}+\bm{\lambda}^{\mathsf{T}}\left(\frac{\partial\bm{K}_{\texttt{T}}}{\partial\rho_{e}}\bm{T}+\bm{K}_{\texttt{T}}\frac{\partial\bm{T}}{\partial\rho_{e}}\right)$
(33) $\displaystyle=\frac{\partial
f_{j}}{\partial\rho_{e}}+\left(-\left(\frac{\partial
f_{j}}{\partial\bm{t}}\right)^{\mathsf{T}}\bm{G}+\bm{\lambda}^{\mathsf{T}}\bm{K}_{\texttt{T}}\right)\frac{\partial\bm{T}}{\partial\rho_{e}}+\bm{\lambda}^{\mathsf{T}}\frac{\partial\bm{K}_{\texttt{T}}}{\partial\rho_{e}}\bm{T}.$

By using the same Lagrangian multiplier as solved from Eq. 26, the above
equation is simplified into

$\displaystyle\derivative{f_{j}}{\rho_{e}}$ $\displaystyle=\frac{\partial
f_{j}}{\partial\rho_{e}}+\bm{\lambda}^{\mathsf{T}}\frac{\partial\bm{K}_{\texttt{T}}}{\partial\rho_{e}}\bm{T}.$
(34)

The first part is directly related to $\phi_{e}$, and can also be calculated
by using adjoint analysis:

$\displaystyle\frac{\partial f_{j}}{\partial\rho_{e}}$
$\displaystyle=\frac{\partial
f_{j}}{\partial\rho_{e}}+\bm{\gamma}^{\mathsf{T}}\partialderivative{\left(\bm{K}^{\\{j\\}}\bm{U}^{\\{j\\}}-\bm{F}^{\\{j\\}}\right)}{\rho_{e}}$
(35) $\displaystyle=\left(\frac{\partial
f_{j}}{\partial\bm{U}^{\\{j\\}}}+\bm{\gamma}^{\mathsf{T}}\bm{K}^{\\{j\\}}\right)\frac{\partial\bm{U}^{\\{j\\}}}{\partial\rho_{e}}+\bm{\gamma}^{\mathsf{T}}\frac{\partial\bm{K}^{\\{j\\}}}{\partial\rho_{e}}\bm{U}^{\\{j\\}}$
$\displaystyle=\bm{\gamma}^{\mathsf{T}}\frac{\partial\bm{K}^{\\{j\\}}}{\partial\rho_{e}}\bm{U}^{\\{j\\}}.$

Here $\bm{\gamma}$ is same the Lagrangian multiplier as in Eq. 30.

The sensitivity analysis of these volume constraint functions $g_{0}$ and
$g_{j}$ is a simpler version of the analysis above, as they only have explicit
parts related to their corresponding $\bm{\rho}^{\\{j\\}}$.

##  4 Results and discussion

In this section, we first demonstrate that the effectiveness of the proposed
method is independent of the initialization of the optimization variables, in
contrast to prior formulations with explicit constraints (Section 4.1). We
then evaluate our method for different application scenarios, including the
minimization of gravity-induced distortion of intermediate structures during
construction (Section 4.2), and minimization of thermal-induced distortion of
fixed structural layout (Section 4.3) and by concurrent optimization of the
structural layout and fabrication sequence (Section 4.4). Lastly, we test the
proposed method for space-time topology optimization of 3D structures (Section
4.5). In all experiments, the Young’s modulus was fixed at 1.0, and the
Poisson’s ratio was fixed at 0.3. Unless stated otherwise, in all 2D tests, we
set 20 layers and a volume fraction of 0.6. The method was implemented in
Matlab.

###  4.1 Validation of avoiding local minima

The first example is intended to validate that our new formulation avoids the
appearance of local minima in the time field. The test is performed on a
cantilever beam, as illustrated in Fig. 7. The design domain is discretized
into a regular quadrilateral mesh with the resolution of $210\times 140$. The
objective is to reduce the compliance of the entire component, under a load on
the top-right corner. Process-dependent loads on intermediate structures are
not considered in this test, by setting $\alpha_{j}$ in Eq. 15 to zero. To
test the robustness of the new formulation, we run space-time topology
optimization with different initializations of the optimization variables. We
study in particular different initializations of the thermal diffusivity field
$\bm{\mu}$. For the density-related field $\bm{\phi}$, we initialize it with a
uniform density of the target volume fraction.

Figure 7: Illustration of the problem setup. Left: the build plate (green) is
placed at the bottom of the design domain, meaning the structure is to be
constructed from the bottom upwards. Right: the boundary condition of the
structure in its intended use scenario.

Figure 8 compares the results obtained using the original formulation [25] and
the proposed new formulation, under different initializations. All results
were generated with 500 optimization iterations. The three initializations are
shown in the first column, from top to bottom, random, monotonic, and uniform.
The second and third columns visualize the optimized time field and
corresponding fabrication sequence, respectively, obtained using the original
formulation where the optimization variables are the time field [25]. Under
random as well as uniform initialization of the time field, the optimized time
fields exhibit multiple local minima and maximum. Only the monotonic
initialization results in a clear fabrication sequence. The last two columns
visualize the optimized results with the proposed regularization. Here,
regardless of the initialization of the thermal diffusivity field $\bm{\mu}$,
the optimized time fields are monotonic. The time values increase, within the
shape of the optimized structural designs, from the bottom which is prescribed
as the build plate, to the up boundary of the domain. The areas outside of the
optimized structural designs take a time value of $1$, i.e., the upper bound
of time values. The monotonic increase of the time field is also clearly
visible from the optimized sequence shown on the right-hand side. In all
tests, the compliance values are comparable. This is because the time field in
this example is only used to segment the evolving structural layout into
layers of equal areas, and no sequence-dependent loads are applied.

Fig. 9 plots the convergence curves of the objective function (left) and
volume constraint per layer (right), corresponding to the optimization with a
uniform thermal diffusivity field. The optimization stabilizes after about 200
iterations. A continuation is applied to the projection parameter $\beta_{t}$
in Eq. 2. From these curves, we can see, that the objective function converges
well and all volume constraints are satisfied.

Figure 8: Comparison of space-time topology optimization with the time field
as optimization variables [25] and using the proposed regularization. The
three different initializations are, from top to bottom, random, monotonic,
and uniform. When using the time field as optimization variables, the time
field is initialized, while in our proposed method, the initialization is
applied to the thermal diffusivity field $\bm{\mu}$. Using the proposed
regularization, the optimized time field is monotonic, regardless of the
initial values. Figure 9: The convergence curves of space-time optimization
using the proposed regularization. Left: the objective function. Right: the
per-layer volume constraints.

###  4.2 Space-time topology optimization under gravity loads

Our next example is the concurrent optimization of the structural design and
fabrication sequence under gravity loads, a problem that has been studied in
Ref. [25]. During fabrication, the self-weight of the intermediate structure
may result in large deflection. The weight of the intermediate structure
increases as the fabrication progresses, and it is dependent on the
fabrication sequence. In this example, the objective function consists of two
parts. The first part is the compliance of the final structure, under the
boundary condition as illustrated in Fig. 7(right). The second part is the sum
of the compliance values of intermediate structures under self-weight. These
two parts are integrated with a weighting factor $\alpha$ in Eq. 15. We note
that the boundary condition for simulating the effects of gravity on the
intermediate structure is different than that for simulating the load-bearing
capacity of the complete structure.

Fig. 10 shows the optimization results with five different weighting factors.
From left to right, $0.0$, $0.1$, $2.0$, $5.0$ and $10.0$. The first row shows
the optimized time fields. The second row shows the optimized structure as
well as the fabrication sequence according to the time field. In all five
cases, the fabrication sequence starts from the bottom of the domain where the
build plate is located, and gradually moves upwards. Comparing the five
results from left to right, as the weighting factor for intermediate
structures increases, the optimized structure exhibits a progressively larger
contact region with the build plate and concentrates more material in its
lower part.

Figure 10: Space-time topology optimization with gravity loads on intermediate
structures. From left to right, the weighting factor for intermediate
structures increases from $0.0$ to $10.0$. Top: the optimized time fields.
Bottom: optimized structural layouts with the corresponding fabrication
sequence.

###  4.3 Fabrication sequence optimization for minimizing thermal-induced
distortion

In the third set of experiments, we validate our method in fabrication
sequence optimization for minimizing the distortion induced by process-
dependent thermomechanical loads. Distortion is a major issue in metal
additive manufacturing. It is caused by the high energy input and large
temperature gradients acting on the materials during fabrication. In a prior
work [26], the time field was directly optimized to obtain the fabrication
sequence. As shown in Fig. 3, the feasibility of the optimized sequence
depends on the initial values of the time field.

We test the proposed regularization on sequence optimization for three
components, shown in Fig. 11. The green line at the bottom of each component
indicates the build plate. The objective is to minimize thermal-induced
distortion. In (a) the distortion is represented by the flatness of the blue
edge, measured by the relative vertical displacements of nodes on the edge. In
(b) the distortion is represented by the flatness of a horizontal and a
vertical edge, both indicated by blue lines. In (c) the distortion is
represented by the relative vertical displacements of two nodes on the top
edge, indicated by blue dots. The distortion due to process-dependent
thermomechanical loads is simulated using the inherent strain method. We refer
readers to [26] for a detailed explanation of the simulation of the layered
manufacturing process.

The optimized results are visualized in Fig. 12. From left to right, the
optimized time field, corresponding fabrication sequence, and resulting
distortion. The optimized time field of each component shows a continuous
increase from the bottom of the domain to an edge that is furthest. This
continuous change can also be observed in the optimized fabrication sequence.
In all three examples, the measured distortion values are very small. From the
enlarged figures on the right, it can be seen that the intended flatness of
the edges is effectively achieved.

Figure 11: Illustration of the three components used in fabrication sequence
optimization. The build plate (green) is located at the bottom of each domain.
The objective is to minimize the thermal-induced distortion, measured on the
geometric features marked in blue. Figure 12: Results of fabrication sequence
optimization for distortion minimization. From left to right: the optimized
time field, fabrication sequence, and simulated displacement field.

###  4.4 Concurrent structure and sequence optimization with constraints on
thermal-induced distortion

To further demonstrate the effectiveness of the proposed method, we
concurrently optimize the structural layout and fabrication sequence. The test
was performed on the design domains shown in Fig. 12(a) and (c). The objective
is to minimize the compliance of the final structure, under the boundary
condition illustrated in Fig. 13(left). The optimization is subject to a bound
on the thermal-induced distortion. The measurement of distortion follows the
previous subsection. The target volume fractions of the two components are 0.5
and 0.6, respectively.

Fig. 13 visualizes the optimized time fields, sequences, and displacements.
While the optimized structural layouts are very complex, the time fields are
monotonic within the solid regions, demonstrating the effectiveness of the
proposed method.

Figure 13: Results of space-time topology optimization under constraints on
thermal-induced distortion. From left to right: the optimized time field,
fabrication sequence, and simulated displacement field.

###  4.5 Space-time topology optimization of 3D structures

The proposed regularization method is also tested on space-time topology
optimization of 3D structures. The design domain and boundary condition are
illustrated in Fig. 14 (a). The left side of the domain is fixed, while a
distributed load is applied on the top-right edge. The objective is to
minimize compliance under this boundary condition. This is subject to a
constraint on thermal-induced distortion, which is represented by the flatness
of the edge where the external load is applied, i.e., the relative
displacement of nodes on this edge along the $x$\- and $z$-axis due to
process-dependent thermomechanical loads. The build plate is situated at the
bottom of the domain, meaning the construction of the structure will proceed
from the bottom up. The design domain is discretized into hexahedral finite
elements with the resolution of $60\times 20\times 30$. The volume fraction is
0.2, and the number of layers is 10.

Figure 14: Space-time topology optimization of a 3D structure. (a): the design
domain and boundary condition. (b): the optimized time field. (c): simulation
of the distortion in the optimized structural layout with the corresponding
fabrication sequence. (d): Simulation of intermediate structures during the
fabrication process.

The optimized time field is shown in Fig. 14 (b), while the optimized
structure under thermal-induced distortion is visualized in (c). The flatness
of the target edge is clearly visible. Also shown in (c) are the segmented
layers, indicated by the varying colors. A simulation of the fabrication
sequence is shown in (d).

##  5 Conclusions

In this paper, we present a novel method for regularizing the time field in
space-time topology optimization. Instead of directly treating the time field
as optimization variables, we employ a static heat equation with a drain term
to derive a monotonic pseudo-temperature field. We then consider spatially
varying thermal diffusivity in the heat equation as optimization variables.
This approach effectively prevents local minima from appearing in the
optimized time field, regardless of the initial values used in optimization.
This approach is easily extendable to deal with 3D components. We demonstrate
its effectiveness in additive manufacturing scenarios, which involve process-
dependent loads such as gravity and thermomechanical loads on intermediate
structures.

A monotonically progressing time field is fundamental for ensuring the
viability of the interpreted fabrication sequence. This work thus lays a solid
foundation for further advancements in space-time optimization for multi-axis
additive manufacturing. As part of our future work, we are exploring methods
to regularize the thickness of curved layers. While multi-axis additive
manufacturing facilitates curved layers, achieving uniform thickness, or near
uniformity, remains a desirable goal.

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). Kai Wu is supported by the China Scholarship Council
(CSC). This work is also partly supported by the Dutch Research Council (NWO)
with project number 20382, “Space-time optimization for additive
manufacturing”.

## References

  * Graziosi et al. [2023] S. Graziosi, J. Wu, Z. Doubrovski, N. Meisel, M. Schütz, T. Stankovic,  Advancing design for additive manufacturing education: A focus on computational skills and competencies,  ETH Learning and Teaching Journal 4 (2023) 5–21. doi:10.16906/lt-eth.v4i1.220.
  * Brackett et al. [2011] D. Brackett, I. Ashcroft, R. Hague,  Topology optimization for additive manufacturing,  in: 2011 International Solid Freeform Fabrication Symposium, University of Texas at Austin, 2011.
  * Zegard and Paulino [2016] T. Zegard, G. H. Paulino,  Bridging topology optimization and additive manufacturing,  Structural and Multidisciplinary Optimization 53 (2016) 175–192. doi:10.1007/s00158-015-1274-4.
  * Langelaar [2016] M. Langelaar,  Topology optimization of 3d self-supporting structures for additive manufacturing,  Additive Manufacturing 12 (2016) 60–70. doi:10.1016/j.addma.2016.06.01.
  * Gaynor and Guest [2016] A. T. Gaynor, J. K. Guest,  Topology optimization considering overhang constraints: Eliminating sacrificial support material in additive manufacturing through design,  Structural and Multidisciplinary Optimization 54 (2016) 1157–1172. doi:10.1007/s00158-016-1551-x.
  * Qian [2017] X. Qian,  Undercut and overhang angle control in topology optimization: A density gradient based integral approach,  International Journal for Numerical Methods in Engineering 111 (2017) 247–272. doi:10.1002/nme.5461.
  * van de Ven et al. [2018] E. van de Ven, R. Maas, C. Ayas, M. Langelaar, F. van Keulen,  Continuous front propagation-based overhang control for topology optimization with additive manufacturing,  Structural and Multidisciplinary Optimization 57 (2018) 2075–2091. doi:10.1007/s00158-017-1880-4.
  * Liu et al. [2015] S. Liu, Q. Li, W. Chen, L. Tong, G. Cheng,  An identification method for enclosed voids restriction in manufacturability design for additive manufacturing structures,  Frontiers of Mechanical Engineering 10 (2015) 126–137. doi:10.1007/s11465-015-0340-3.
  * Donoso et al. [2022] A. Donoso, E. Aranda, D. Ruiz,  A new approach based on spectral graph theory to avoiding enclosed holes in topology optimization,  Computer Methods in Applied Mechanics and Engineering 393 (2022) 114769. doi:10.1016/j.cma.2022.114769.
  * Liu et al. [2018] 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 (2018) 2457–2483. doi:10.1007/s00158-018-1994-3.
  * Bayat et al. [2023] M. Bayat, O. Zinovieva, F. Ferrari, C. Ayas, M. Langelaar, J. Spangenberg, R. Salajeghe, K. Poulios, S. Mohanty, O. Sigmund, J. Hattel,  Holistic computational design within additive manufacturing through topology optimization combined with multiphysics multi-scale materials and process modelling,  Progress in Materials Science 138 (2023) 101129. doi:10.1016/j.pmatsci.2023.101129.
  * Allaire et al. [2017] 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. doi:10.1016/j.jcp.2017.09.041.
  * Amir and Mass [2018] O. Amir, Y. Mass,  Topology optimization for staged construction,  Structural and Multidisciplinary Optimization 57 (2018) 1679–1694. doi:10.1007/s00158-017-1837-7.
  * Bruggi et al. [2018] M. Bruggi, N. Parolini, F. Regazzoni, M. Verani,  Topology optimization with a time-integral cost functional,  Finite Elements in Analysis and Design 140 (2018) 11–22. doi:10.1016/j.finel.2017.10.011.
  * Haveroth et al. [2022] G. Haveroth, C.-J. Thore, M. Correa, R. Ausas, S. Jakobsson, J. Cuminato, A. Klarbring,  Topology optimization including a model of the layer-by-layer additive manufacturing process,  Computer Methods in Applied Mechanics and Engineering 398 (2022) 115203. doi:10.1016/j.cma.2022.115203.
  * 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.
  * 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.
  * Liu and To [2017] J. Liu, A. C. To,  Deposition path planning-integrated structural topology optimization for 3d additive manufacturing subject to self-support constraint,  Computer-Aided Design 91 (2017) 27–45. doi:10.1016/j.cad.2017.05.003.
  * Dapogny et al. [2019] C. Dapogny, R. Estevez, A. Faure, G. Michailidis,  Shape and topology optimization considering anisotropic features induced by additive manufacturing processes,  Computer Methods in Applied Mechanics and Engineering 344 (2019) 626–665. doi:10.1016/j.cma.2018.09.036.
  * Boissier et al. [2020] M. Boissier, G. Allaire, C. Tournier,  Additive manufacturing scanning paths optimization using shape optimization tools,  Structural and Multidisciplinary Optimization 61 (2020) 2437–2466. doi:10.1007/s00158-020-02614-3.
  * Boissier et al. [2023] M. Boissier, G. Allaire, C. Tournier,  Concurrent shape optimization of the part and scanning path for powder bed fusion additive manufacturing,  SIAM Journal on Control and Optimization 61 (2023) 697–722. doi:10.1137/21M1461976.
  * 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.
  * Wang et al. [2023] W. Wang, F. van Keulen, J. Wu,  Fabrication sequence optimization for minimizing distortion in multi-axis additive manufacturing,  Computer Methods in Applied Mechanics and Engineering 406 (2023) 115899. doi:10.1016/j.cma.2023.115899.
  * Dai et al. [2018] C. Dai, C. C. L. Wang, C. Wu, S. Lefebvre, G. Fang, Y.-J. Liu,  Support-free volume printing by multi-axis motion,  ACM Trans. Graph. 37 (2018) 134:1–134:14. doi:10.1145/3197517.3201342.
  * Wang et al. [2011] F. Wang, B. S. Lazarov, O. Sigmund,  On projection methods, convergence and robust formulations in topology optimization,  Structural and Multidisciplinary Optimization 43 (2011) 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.
  * Crane et al. [2013] K. Crane, C. Weischedel, M. Wardetzky,  Geodesics in heat: A new approach to computing distance based on heat flow,  Acm Transactions on Graphics 32 (2013) 13–15. doi:10.1145/2516971.2516977.
  * Liu et al. [2015] S. Liu, Q. Li, W. Chen, L. Tong, G. Cheng,  An identification method for enclosed voids restriction in manufacturability design for additive manufacturing structures,  Frontiers of Mechanical Engineering 10 (2015) 126–137. URL: https://dx.doi.org/10.1007/s11465-015-0340-3. doi:10.1007/s11465-015-0340-3.
  * Luo et al. [2020] 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. URL: https://dx.doi.org/10.1016/j.cma.2020.113385. doi:10.1016/j.cma.2020.113385.
  * Zhai et al. [2024] X. Zhai, W. Wang, F. Chen, J. Wu,  Topology optimization of differentiable microstructures,  Computer Methods in Applied Mechanics and Engineering 418 (2024) 116530. doi:10.1016/j.cma.2023.116530.
  * Lazarov and Sigmund [2011] B. S. Lazarov, O. Sigmund,  Filters in topology optimization based on helmholtz-type differential equations,  International Journal for Numerical Methods in Engineering 86 (2011) 765–781. doi:10.1002/nme.3072.
  * 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.


Genus expanded cut-and-join operators and generalized Hurwtiz numbers

Quan Zheng111E-mail<EMAIL_ADDRESS>Partially supported by NSFC
11171258 and NSFC 11571242.

Mathematics College, Sichuan University

610064, Chengdu, Sichuan, PRC

November 1, 2015

(Dedicated to Prof. An-Min Li for his 70th birthday)
Abstract

To distinguish the contributions to the generalized Hurwitz number of the
source Riemann surface with different genus, by observing carefully the
symplectic surgery and the gluing formulas of the relative GW-invariants, we
define the genus expanded cut-and-join operators. Moreover all normalized the
genus expanded cut-and-join operators with same degree form a differential
algebra, which is isomorphic to the central subalgebra of the symmetric group
algebra. As an application, we get some differential equations for the
generating functions of the generalized Hurwitz numbers for the source Riemann
surface with different genus, thus we can express the generating functions in
terms of the genus expanded cut-and-join operators.

Keyword: genus expanded cut-and-join operator; differential algebra; Hurwitz
number; generating function
Subject Classification: 14N10, 14N35, 05E10

###### Contents

  1. 1 Introduction
  2. 2 Generalized Hurwitz numbers

##  1 Introduction

Hurwitz Enumeration Problem has been studied for more than one century, which
has been extensively applied to mathematics and physics. It is an interesting
project for mathematic physics, geometry, representation theory, integrable
systems, Hodge integral, etc. Hurwitz numbers have many different expressions
in the different fields, for example, [[AMMN]], [[DLM]], [[ELSV]], [[GJ1]],
[[GJ2]],[[GJ3]], [[H]], [[MMN]], [[MMN1]], [[MMN2]]. One of the important
geometric tools to deal with Hurwitz numbers is the so-called symplectic
surgery: cutting and gluing [[IP]], [[L]],[[LR]], in the views of algebra and
differential equations, which is equivalent to the so-called cut-and-join
operators [[GJ1]], [[GJ2]], [[GJ3]],[[LZZ]], [[OP]]. The standard cut-and-join
operators can be used to deal with the almost simple Hurwitz numbers and the
almost simple double Hurwitz numbers [[GJ1]], [[GJ2]],[[GJ3]], [[LZZ]],
[[OP]]. A. Mironov, A. Morozov, and S. Natanzon, etc, have defined the
generalized cut-and-join operators in terms of the matrix Miwa variable
[[AMMN]], [[MMN]], [[MMN1]]. However, the generating functions obtained by the
generalized cut-and-join operators defined by A. Mironov, A. Morozov, and S.
Natanzon cannot distinguish the contributions of the source Riemann surface
with the different genus, which is obviously interesting and important. In our
paper, by observing carefully the symplectic surgery and the gluing formulas
of the relative GW-invariants developed by A.M. Li and Y.B. Ruan [[LR]], and
by E. Ionel and T. Parker [[IP]], we define the genus expanded cut-and-join
operators to distinguish the contributions of the source Riemann surfaces with
different genus. A little more precisely, we introduce one parameter $z$ to
“mark” the genus, then we associate to every partition $\Delta$ a genus
expanded cut-and-join operator $W(\Delta,z)$ as follows:

$W(\Delta,z)=\sum_{\Gamma,\Gamma^{\prime}\vdash
d}\frac{d!}{|C_{\Gamma^{\prime}}|}z^{d+l(\Gamma^{\prime})-l(\Delta)-l(\Gamma)}\mu_{0}^{h_{+},d}(\Gamma^{\prime},\Delta,\Gamma)p_{\Gamma}\frac{\partial}{\partial
p_{\Gamma^{\prime}}}.$ (1)

For the precise meaning of the notations in the above formula, see the section
2 and the section 3. Moreover, the power of $z$ is exactly the “lost” genus
after we have executed the symplectic cutting on the target Riemann surface,
referring to Remark 4.3.

One of our main theorems is the followings, referring the theorem 3.5:

Theorem A For a given $d$, as operators on the functions of the time-variables
$p=(p_{1},p_{2},\cdots)$, we have

${W(\Delta_{1},z)W(\Delta_{2},z)=\sum_{\Delta_{3}}z^{d-l(\Delta_{1})-l(\Delta_{2})+l(\Delta_{3})}C_{\Delta_{1}\Delta_{2}}^{\Delta_{3}}W(\Delta_{3},z)},$
(2)

where

$C_{\Delta_{1}\Delta_{2}}^{\Delta_{3}}=\frac{d!}{|C_{\Delta_{3}}|}\mu_{0}^{h,d}(\Delta_{1},\Delta_{2},\Delta_{3})$
(3)

is the structure constants of the center subalgebra of the group algebra
$\mathbb{C}[S_{d}]$.

A more interested result in this paper is that all normalized the genus
expanded cut-and-join operators with same degree form a differential algebra,
which is isomorphic to the central subalgebra of the symmetric group algebra,
referring Corollary (LABEL:cor_3.6).

As an application of the genus expanded cut-and-join operators, we define the
generating function as (referring the section 4)

$\displaystyle\Phi_{g}\\{z|(u_{1},\Delta_{1}),\cdots,(u_{n},\Delta_{n})|p^{(1)},\cdots,p^{(k)},p\\}$
(4) $\displaystyle=$ $\displaystyle\sum_{l_{1},\cdots,l_{n}\geq
0}\sum_{\Gamma,\Gamma_{1},\ldots,\Gamma_{k}}z^{2h-2}\mu_{g}^{h,d}(\underbrace{\Delta_{1},\cdots,\Delta_{1}}_{l_{1}},\cdots,\underbrace{\Delta_{n},\cdots,\Delta_{n}}_{l_{n}},\Gamma_{1},\cdots,\Gamma_{k},\Gamma)[\prod_{j=1}^{n}\frac{u_{j}^{l_{j}}}{l_{j}!}][\prod_{i=1}^{k}p^{(i)}_{\Gamma_{i}}]p_{\Gamma}$
$\displaystyle=$ $\displaystyle\sum_{l_{1},\cdots,l_{n}\geq
0}\sum_{\Gamma,\Gamma_{1},\ldots,\Gamma_{k}}z^{2h-2}({\frac{dim{\lambda}}{d!}})^{2-2g}[\prod_{\makebox{j=1}}^{n}(\phi_{\lambda}(\Delta_{j}))^{l_{j}}\frac{(u_{j})^{l_{j}}}{l_{j}!}]\phi_{\lambda}(\Gamma)[\prod_{i=1}^{k}\phi_{\lambda}(\Gamma_{i})p^{(i)}_{\Gamma_{i}}]p_{\Gamma}.$

Then we have our another main theorem (referring theorem 4.2):
Theorem B For any $i$, we have

$\displaystyle\frac{\partial\Phi_{g}\\{z|(u_{1},\Delta_{1}),\cdots,(u_{n},\Delta_{n})|p^{(1)},\cdots,p^{(k)},p\\}}{\partial
u_{i}}$ (5) $\displaystyle=$ $\displaystyle
W(\Delta_{i},z)\Phi_{g}\\{z|(u_{1},\Delta_{1}),\cdots,(u_{n},\Delta_{n})|p^{(1)},\cdots,p^{(k)},p\\}.$

Thus we can express the generating functions in terms of the genus expanded
cut-and-join operators once we know the initial values. Moreover, the
generating functions can be directly expressed by the expansion in terms of
the genus of the source Riemann surface, referring to Example 4.5.

Our paper is deeply inspired by A. Mironov, A. Morozov, and S. Natanzon’s
works [[AMMN]], [[MMN]], [[MMN1]], [[MMN2]].

##  2 Generalized Hurwitz numbers

Let $\Sigma^{h}$ be a compact (possibly disconnected) Riemann surface of genus
$h,$ and $\Sigma^{g}$ a compact connected Riemann surface of genus $g$. For a
given point set $\\{q_{1},\cdots,q_{k}\\}\in\Sigma^{g}$, which is called the
set of branch points, we call a holomorphic map $f:\Sigma^{h}\to\Sigma^{g}$ a
ramified covering of $\Sigma^{g}$ of degree $d\geq 0$ by $\Sigma^{h}$ with a
ramification type $(\Delta_{1},\cdots,\Delta_{k})$, if the preimages of
$f^{-1}(q_{i})=\\{p_{i}^{1},\cdots,p_{i}^{m_{i}}\\}$ with orders
$\Delta_{i}=(\delta^{i}_{1},\cdots,\delta^{i}_{m_{i}})$ for $i=1,\cdots,k,$
respectively. Two ramified coverings $f_{1}$ and $f_{2}$ with type
$(\Delta_{1},\cdots,\Delta_{k})$ are said to be equivalent if there is a
homeomorphism $\pi:\Sigma^{h}\to\Sigma^{h}$ such that $f_{1}=f_{2}\circ\pi$
and $\pi$ preserves the preimages and the ramification type of $f_{1}$ and
$f_{2}$ at each point $q_{i}\in\Sigma^{g}$. Let
$\mu_{g}^{h,d}(\Delta_{1},\cdots,\Delta_{k})$ be the number of equivalent
covering of $\Sigma^{g}$ by $\Sigma^{h}$ with ramification type
$(\Delta_{1},\cdots,\Delta_{k})$, which is called the generalized Hurwitz
number. Note that $\mu_{g}^{h,d}(\Delta_{1},\cdots,\Delta_{k})$ is nonzero
only if the Hurwitz formula

$(2-2g)d-(2-2h)=\sum_{i=1}^{k}(d-m_{i})$ (6)

holds, where $m_{i}$ is the numbers of the preimage of $q_{i}$. How to
determine $\mu_{g}^{h,d}(\Delta_{1},\cdots,\Delta_{k})$ is known as the
generalized Hurwitz Enumeration Problem [[H]]. In general, $\Delta_{i}$ is not
required to be a partition of $d$, for simplicity, throughout this paper, we
will assume it. Denote by $\Delta_{i}\vdash d$ if $\Delta_{i}$ is a partition
of $d.$

Let $S_{d}$ be the symmetric group of $d$ letters, $\mathbb{C}[S_{d}]$ be the
group algebra of $S_{d}$ and $Z\mathbb{C}[S_{d}]$ be the central subalgebra of
$\mathbb{C}[S_{d}]$, then in the terms of algebra, the story begins from a
simple fact as the following formula [[H]],

$\displaystyle\mu_{g}^{h,d}(\Delta_{1},\cdots,\Delta_{k})=\frac{1}{d!}[1]\prod_{j=1}^{g}\prod_{a_{j},b_{j}\in
S_{d}}[a_{j},b_{j}]C_{\Delta_{1}}\cdots C_{\Delta_{k}},$ (7)

where $[1]$ means that we take the coefficient of the identity of the product
of $g$-tuple commutators $\prod_{a,b\in S_{d}}[a,b]$ and $k$-tuple central
elements $C_{\Delta_{1}},\cdots,C_{\Delta_{k}}\in Z\mathbb{C}[S_{d}]$
corresponding to the partitions $\Delta_{1},\cdots,\Delta_{k}.$ Thus,
according to the theory of representations of the symmetric group $S_{d}$, we
arrive at the formula [[OP]],

$\displaystyle{}\mu_{g}^{h,d}(\Delta_{1},\cdots,\Delta_{k})=\sum_{\lambda\vdash
d}({\frac{dim{\lambda}}{d!}})^{2-2g}\phi_{\lambda}(\Delta_{1})\cdots\phi_{\lambda}(\Delta_{k}),$
(8)

which expresses them through the properly normalized symmetric group
characters

$\phi_{\lambda}(\Delta)=\frac{1}{dim\lambda}|C_{\Delta}|\chi_{\lambda}(\Delta),$
(9)

where $\lambda$ is a Young diagram of degree $d$, and $dim\lambda$ is the
dimension of the irreducible representation of the symmetric group $S_{d}$
corresponding to $\lambda$, $\chi_{\lambda}(\Delta)$ is the character of a
permutation $\sigma\in C_{\Delta}$ under the irreducible representation
$\lambda$, moreover, $|C_{\Delta}|$ is the number of the permutations of
$S_{d}$ with cyclic type $\Delta$.

Let $p=(p_{1},p_{2},p_{3},\cdots,)$ be indeterminantes, which are called the
time-variables, and we assume that $\Delta=(\delta_{1},\cdots,\delta_{n})$ is
a partition. We denote

$l(\Delta):=n,$ $m_{r}(\Delta):=\\#\\{i|\delta_{i}=r,\delta_{i}\in\Delta\\},$
$\Delta!:=\prod_{r\geq 1}m_{r}(\Delta)!,$ $p_{\Delta}:=p_{\delta_{1}}\cdots
p_{\delta_{n}},$

and

$\frac{\partial}{\partial
p_{\Delta}}:=\frac{1}{\Delta!}\frac{\partial}{\partial
p_{\delta_{1}}}\cdots\frac{\partial}{\partial p_{\delta_{n}}}.$

For any partition $\lambda\vdash d,$ it is well known that the symmetric group
characters $\phi_{\lambda}(\Delta)$ are related to the Schur functions
$S_{\lambda}\\{p\\}$ as following [[MMN]],

$S_{\lambda}\\{p\\}=\sum_{\Gamma^{\prime}}\frac{dim{\lambda}}{d!}\phi_{\lambda}(\Gamma^{\prime})p_{\Gamma^{\prime}}.$
(10)

Our first important lemma is about the associativity of the generalized
Hurwitz number:

######  Lemma 2.1.

[[MMN2]] For any positive integer $k>l\geq 1$, we have

$\mu_{g}^{h,d}(\Delta_{1},\cdots,\Delta_{k})=\sum_{g_{1}+g_{2}=g}\sum_{\Delta^{\prime}}\mu_{g_{1}}^{h_{1},d}(\Delta_{1},\cdots,\Delta_{l},\Delta^{\prime})\frac{d!}{|C_{\Delta^{\prime}}|}\mu_{g_{2}}^{h_{2},d}(\Delta^{\prime},\Delta_{l+1},\cdots,\Delta_{k}).$
(11)