Space-time topology optimization for additive manufacturing

The design of optimal structures and the planning of (additive manufacturing) fabrication sequences have been considered typically as two separate tasks that are performed consecutively. In the light of recent advances in robot-assisted (wire-arc) additive manufacturing which enable addition of material along curved surfaces, we present a novel topology optimization formulation which concurrently optimizes the structure and the fabrication sequence. For this, two sets of design variables, i.e., a density field for defining the structural layout, and a time field which determines the fabrication process order, are simultaneously optimized. These two fields allow to generate a sequence of intermediate structures, upon which manufacturing constraints (e.g., fabrication continuity and speed) are imposed. The proposed space-time formulation is general, and is demonstrated on three fabrication settings, considering self-weight of the intermediate structures, process-dependent critical loads, and time-dependent material properties.


Introduction
Recent advances in additive manufacturing (AM, also known as 3D printing) enable the fabrication of structures with unprecedented geometric complexity. The benefits of this manufacturing flexibility are probably best exploited in combination with the design of structures by topology optimization (TO). TO aims at finding the optimal distribution of material under a given set of constraints (Bendsøe and Sigmund 2004). The optimized structures are often very complex from a geometric perspective. Without applying additional constraints to reduce complexity, the optimized structures are difficult or impossible to produce by conventional manufacturing technologies. Such extra constraints nevertheless compromise the structural optimality. In the past years, impressive progress has been made on the integration of topology optimization for additive manufacturing. For an overview of research on this topic, we refer to a recent survey article by Liu et al. (2018a). In particular, the developments have been focusing on AM constraints and/or characteristics such as the overhang angle (e.g., Wu et al. 2016b;Gaynor and Guest 2016;Mirzendehdel and Suresh 2016;Qian 2017;Langelaar 2017Langelaar , 2018van de Ven et al. 2018;Garaigordobil et al. 2018;Allaire et al. 2017a;Wang et al. 2019), infill structures (e.g., Wu et al. 2017Wu et al. , 2018Groen et al. 2019;Fu et al. 2018;Clausen et al. 2016;Garner et al. 2019), thermal residual stresses (e.g., Allaire and Jakabčin 2018), and material anisotropy, i.e., due to the deposition path (e.g., Liu et al. 2018b;Dapogny et al. 2019).
In additive manufacturing, structures are fabricated progressively, i.e., by adding material in an incremental manner. The fabrication sequence is typically planned after the structure has been designed or optimized. In commonly used AM processes such as fused deposition modeling and selective laser sintering, given an optimized structure with a specific orientation, the structure is sliced into a set of planar layers. The planar layers are added to the structure one upon another, by extruding small flattened strings of molten material or by melting and fusing powder material.
The AM platforms often have three degrees of freedom, allowing three dimensional translation of the printer head or the structure under construction. The introduction of rotational degrees of freedom into AM platforms has further increased the fabrication flexibility. For instance, using a robotic arm to continuously rotate the structure during construction, it allows to deposit material along curved layers (Dai et al. 2018). The increased flexibility in production further enlarges the design space with the planning of the fabrication sequence.
As mentioned, the optimization of structures and the planning of the fabrication sequence are typically performed separately. In topology optimization, it mostly concerns the mechanical performance of the final structure as a whole, and does not evaluate the mechanical properties of the unfinished structure during the fabrication process. Consider the fabrication of a large scale structure using wire and arc additive manufacturing (Williams et al. 2016). The mechanical properties of the structure at all intermediate stages shall also comply with certain criteria. In general, a number of aspects, including self-weight, material curing and solidification, thermal dissipation and distortion, are influenced by the fabrication sequence. These aspects in turn affect the (mechanical) performance of the structure at both the intermediate and final stages.
In this paper, we make the first step towards the concurrent optimization of structural layout and the corresponding fabrication sequence, which we shall call space-time topology optimization. The space-time topology optimization uses two sets of design variables. The first set represents the structural layout by a density field which is standard, as in traditional density-based approaches (Sigmund 2001). The second set encodes the fabrication sequence by a time field, with the ascending order indicating the incremental addition of structural material. We present a general formulation where the objective function could take into account the structural properties of both intermediate structures as well as the complete structure. To this end, a sequence of intermediate structures is defined by the density and the time field. We impose general constraints on intermediate structures, regarding fabrication continuity and process speed. This general formulation is demonstrated by integrating a few simplified yet meaningful aspects that are associated with the fabrication sequence, including self-weight of the intermediate structures, process-dependent loads, and timedependent mechanical properties (e.g., in a curing process).
The present work is related to a few recent papers which dealt with prescribed fabrication sequence in topology optimization. Allaite et al. (2017a, b) proposed a novel mechanical constraint functional, which mimics the layer by layer construction process featured by additive manufacturing technologies. This constraint aggregates the self-weights of all the intermediate structures. Amir and Mass (2018) proposed a formulation which integrates the loading and support conditions during construction. This formulation effectively reduces temporary supports or scaffolds for fabricating the optimized layouts. Bruggi et al. (2018) developed a formulation for optimizing support structures in problems involving a time-dependent construction process. Allaire and Jakabčin (2018) introduced a model for shape and topology optimization, taking into account the effects of the thermal stresses on intermediate structures during the additive construction process. In the approaches described above, the fabrication sequence is prescribed, and in particular, a sequence of planar layers. In contrast to these approaches, in the present work, the fabrication sequence is optimized together with the structure. The proposed method forms a perfect match with recent advances in additive manufacturing which enable flexible fabrication beyond consecutive planar layers.
We note that the term space-time topology optimization was used by Jensen (2009) in a different context, i.e., to optimize the point-wise, time-dependent material properties for transient problems. It was outlined for one-dimensional wave propagation in an elastic rod, taking time-dependent Young's modulus as design variables. In this paper, the temporal domain is used to encode the fabrication sequence. The structural analysis in our examples concerns a series of static equilibria, evaluated at specific timepoints during the fabrication process.
The remainder of this paper is organized as follows. In Section 2, we present the formulation including the general objective function, and constraints on intermediate structures regarding fabrication continuity and speed. This general formulation is followed by an example to explain the consequences of the constraints (Section 3). We then demonstrate the space-time optimization concept on three scenarios, considering self-weight of the intermediate structures (Section 4), process-dependent loads (Section 5), and time-dependent material properties (Section 6). After a discussion on parameters and alternative formulations in Section 7, we present conclusions in Section 8.

Space-time topology optimization
The structural layout and the fabrication sequence are described by two fields. In addition to a density field (ρ) known from traditional topology optimization, a time field (t) is introduced to encode the order of material deposition. The objective function (J ) is abstractly defined as a function of these two fields, by where the first term (J complete ) measures the structural property (e.g., compliance) of the entire structure, while the second term (J process ) measures properties of intermediate structures during the manufacturing process.
In this section, we first present the generation of intermediate structures from the density and the time field. We then present example constraints reflecting fabrication requirements, i.e., volume constraints and continuity constraints on intermediate structures.

Intermediate structures
Using a finite element discretization of the design space, each element is associated with a (pseudo) density value ρ e ∈ [0, 1] and a (pseudo) time value t e ∈ [0, 1]. The density value indicates whether the element is empty (ρ e = 0) or solid (ρ e = 1) in the final (complete) structure. The time value indicates the time at which the material associated with the element is added to the structure. Thus, a relatively larger time value indicates that this element is fabricated later. As in traditional density-based approaches, the density value eventually converges to either 0 or 1. However, it shall be noted that the time value shall be continuous.
In space-time topology optimization, structural properties may be evaluated not only for the final (complete) structure, but also for the incomplete structure, i.e., at intermediate stages of the manufacturing process. In the proposed space-time topology optimization, at a specific time T , the elements with a time value t e ≤ T have been added to the structure. The intermediate structure at time T is thus determined by This is illustrated in Fig. 1. The discrete density field is illustrated in (a), and the continuous time field in (b), with a few iso-contours. Figure 1(c) and (d) depict the incomplete structure, i.e., intermediate stages of the fabrication process, at T = 0.2, and 0.4, respectively.
The intermediate structure at time T can be extracted by (2). To avoid the use of conditional equations which are not differentiable, we make use of a filtering technique to generate the intermediate structure from the density and the time field. This process is visualized in Fig. 2.
Two sets of design variables, φ for density and τ for time, are used in optimization. First, in order to avoid checkerboard patterns, convolution operators are applied to smooth both fields. This results inφ and t =τ , with the tilde indicating smoothed continuous fields. It is worth noting that both fields need to be smoothed. A checkerboard pattern in one of the fields leads to a checkerboard pattern in intermediate structures, since intermediate structures are specified by the multiplication of the two fields, as will be introduced shortly. We use the convolution operator as in classical density-based approaches for smoothing. This yields and where v i is the area or volume of an element, and the weighting function is defined as where r is the filter radius, x e and x i are the positions of the centroid of element e and its neighbor element i ∈ S e = {i | w(x i , r) > 0}, respectively. We also note that the filter radii, r t for time and r d for density, can take different values. Besides avoiding checkerboard patterns, r d also regulates the thickness of resulting structures. Following the smoothing operator, a (smoothed) Heaviside projection is applied to obtain discrete values ρ =φ andt, with the bar indicating projected discrete fields. For the density field, the projection converts a density value smaller (or larger) than a given threshold to close to 0 (or 1), by where β d is a positive number to control the sharpness of the step function, and η = 0.5 is the density threshold value.  This projection has been discussed, for instance, by Wang et al. (2011). For the time field, a projection is used to convert a time value smaller (or larger) than a given threshold, T in time, to close to 1 (or 0). This is achieved bȳ where β t , similar to β d , controls the projection sharpness, and T is the threshold time value.
The intermediate structure at time T is thus defined by This function can be interpreted as a segmentation of the density field ρ by the iso-contour of t = T , as visualized in Fig. 1.

Volume constraints on intermediate structures
In additive manufacturing, the structure is fabricated incrementally. A parameter involved in this process is the fabrication speed, i.e., the amount of material which can be deposited per unit of time. To incorporate the fabrication speed in the space-time optimization, we divide the time range [0, 1] by a finite number (N + 1) of uniformly distributed timepoints, denoted by Consequently, there are N intervals (also called stages in the following), each with a duration of 1 N . The number of stages (N) is prescribed, and thus the specific time T i when an intermediate structure is evaluated is determined. For simplicity, we assume a constant fabrication speed; the maximum volume of the complete structure (V 0 ) is equally added during each of the uniform time intervals. In other words, the increment in volume during each time interval is bounded by V 0 N , i.e., where V [T i ] and V [T i−1 ] denote the total material volume processed up to timepoints T i and T i−1 , respectively. The initial volume, V [T 0 ] , is prescribed as 0. For compliance minimization as studied in this paper, since the optimization always uses the maximum amount of available material volume, this is equivalent to where v e is the area or volume of an element. Since in this paper a uniform finite element discretization is used, v e is constant for all elements (v e = v 0 ).

Continuity constraints on intermediate structures
During the incremental additive fabrication process, it is for most manufacturing processes essential that the material is deposited on material which has been deposited previously. If the material is to be deposited in isolation, auxiliary structures are temporarily needed to hold such an isolated structural fragment in place. To avoid the processing costs associated with additional supports, we thus formulate a continuity constraint to prevent isolated material patches. An isolated material patch during the fabrication process can be associated with a local minimum in the time field; all its adjacent elements have a larger time value and thus will be fabricated later. Therefore, isolated material patches can be prevented by requiring where N e denotes the set of elements adjacent to element e. M is the set of active elements in the design domain, i.e., all elements except those which are prescribed as the starting point/region for the fabrication process (i.e., with t e = 0).

Relaxation
The continuity constraint (13) is not differentiable, and it applies to a large number of elements. To facilitate numerical optimization, we present an aggregated formulation. This formulation involves two steps.
First, the non-differentiable function (13) is approximated by a continuous function. To this end, we make use of the fact that t ≤ 1 and, consequently, min i∈N e (t i ) can be rewritten as: where max i∈N e (1 − t i ) can be approximated with a pnorm : As a result, g(t e ) is approximated by Second, the per-element constraints (13) are aggregated. A straightforward way to aggregate (16) is by making use of the max function, which can be approximated with a p-norm. However, applying a p-norm on top of another p-norm (i.e (15)), both with p as large as 50, leads to a highly non-linear response function. Our initial numerical tests showed that the optimization convergence using this function is far from ideal. Therefore, we rather aggregate (16) by computing the average of a function defined on g(t e ) where # denotes the number of elements in a set, is a small constant, and the function H is defined by The left-hand side in (18) becomes 1 #(M) if there exists a local minimum in the set of active elements (M), i.e x = g(t e ) > 0. Therefore, by assigning a value that is smaller than 1 #(M) ( = 10 −9 in this work), (18) would effectively avoid local minima.
The step function H is approximated by where β m controls the sharpness of projection. We note that a Heaviside projection-based aggregation has recently been used to control overhang angle (Qian 2017;Wang et al. 2019) and local stresses (Wang and Qian 2018). A detailed comparison between the Heaviside projection-based aggregation and the p-norm is provided in Wang and Qian (2018).

Sensitivity analysis
This subsection presents derivatives of the constraints which we proposed in the previous two subsections.

Sensitivity analysis of volume constraints
For the volume constraint defined in (11), its derivative regarding φ e at time T i is given as: where ∂ρ k ∂φ k follows as and ∂φ k ∂φ e is calculated based on the definition ofφ e in (3). Similarly, the derivative of constraint (11) regarding τ e at time T i is given as: where and ∂t k ∂τ e is calculated with the definition of t e in (4).

Sensitivity analysis of continuity constraints
To simplify the derivation, the constant 1 #(M) in (18) is temporally dropped out. Thus, the derivative of (18) regarding the time variable τ e is where with and

Demonstration of manufacturing constraints
In this section, a simplified space-time optimization problem is used to demonstrate the consequences of the fabrication constraints. The test problem is compliance minimization under the assumption of linear elasticity. The design domain and boundary conditions are illustrated in Fig. 3. In contrast to a single constraint on the total material volume as in traditional compliance minimization, here a series of volume constraints is imposed. Each of them corresponds to the maximum material volume processed up to a specific time, as discussed in Section 2.2. Furthermore, the manufacturing continuity constraint as discussed in Section 2.3 is included. Let N denote the number of prescribed time intervals, the problem is described by where T denotes the transpose operator, U is the displacement vector, K the stiffness matrix, and F the force vector. The stiffness matrix K is assembled from element stiffness matrices defined by k e = E e (ρ e )k 0 , where k 0 is the stiffness matrix of a solid element with unit Young's modulus and E e (ρ e ) is the Young's modulus corresponding to element e, interpolated via the modified solid isotropic material with penalization (SIMP), given by where E 0 is the Young's modulus of a solid element, E min a small term assigned to prevent the global stiffness matrix from becoming singular, and q a penalization Fig. 3 The design domain and boundary conditions. The design domain is discretized into 120×40 square bilinear finite elements factor (typically q = 3). In this test formulation, the objective is to reduce the compliance of the entire structure. Structural properties (e.g., compliance) of intermediate structures are not included here, and will be discussed in following sections where process-dependent loads and material properties are introduced. The optimization problem is solved using the method of moving asymptotes (Svanberg 1987). The derivative of (29) and (31) regarding to the design variables τ e and φ e is standard. The derivatives of (34) and (35) have been given in Section 2.4. Figure 4 compares optimized density and time fields for three configurations where the prescribed start of manufacturing (indicated by blue quads) is different. The first row shows the initial time fields. They are initialized by computing the distance of each element to the start region. The distance field is then normalized by the maximum distance value among all elements. While the time field is initialized by a monotonic field, the density field is initialized uniformly, using the target volume fraction. The second row shows the optimized continuous time fields, where the curves indicate the boundary of different manufacturing stages (a total of 8 in this example). The time fields are used to color the optimized structures, shown in the last row. Although the optimized structures are different, their compliance values are very close. In fact, the time field plays no role in the objective function. It affects only the segmentation of the structures into stages. Figure 5 shows a sequence of intermediate stages, illustrating stages of the additive manufacturing process. This demonstration verifies that the optimized fields satisfy continuity and volume constraints. Figure 6 compares the optimization using two different initial time fields, the top row with a uniform time field (a) and the bottom row with a monotonic time field initialized by a normalized distance field corresponding to the start of the fabrication process (d). The optimized time field in (e) is monotonic from left to right, while the optimized time field in (b) exhibits some variations. Two local maxima in the top and bottom middle of the time field in (b) can be observed. Local maxima do not violate the continuity constraint, but lead to a more complex manufacturing sequence. A monotonic initial time field is used in this paper, if not explicitly stated otherwise. Figure 7 shows optimized results with 10, 20, and 30 time intervals (i.e., manufacturing stages). The fabrication starts from the top left in the design domain. The initial time field has been shown in the last column of Fig. 4. The fabrication granularity increases as the number of time intervals is increased, allowing a finer planning of the fabrication process. It is also verified that the optimized fields satisfy the continuity and volume constraints on intermediate structures.  . 6 Optimization with different initial time fields. The initial time field in (a) is uniform, while the initial time field in (d) is the normalized distance field corresponding to the correct start (blue quad). The corresponding optimized time fields and structures are shown in the same row. The optimized time field (b) has several local maxima. In contrast, using a monotonic field as the initial time field, no extra local maxima are found in (e), even with a smaller number of iterations Fig. 7 Optimized time fields and structures with different numbers of manufacturing stages: 10, 20, and 30 from left to right, respectively. The startpoint is the top left element and the initial time field is shown in Fig. 4 (top right subfigure). The compliances are listed at the bottom

Self-weight of intermediate structures
In additive manufacturing, the deformation of intermediate structures due to gravity may become significant, e.g., fabrication using flexible materials such as TPU (thermoplastic polyurethane). To prevent such deformation, self-weight of intermediate structures has been considered to design selfsupporting structures and supporting scaffolds (Allaire et al. 2017a, b;Bruggi et al. 2018;Amir and Mass 2018). In our formulation, we concurrently optimize the structural layout and its fabrication sequence to minimize the compliance due to self-weight of intermediate structures. The problem is described by where K(ρ [T i

Sensitivity analysis
The sensitivity of the objective function regarding the density design variable is given as follows: where k j (ρ) and U j are the element stiffness matrix and displacement vector of finite element j for the complete structure, k j (ρ [T i ] ) and U [T i ] j are the element stiffness matrix and displacement vector of finite element j at the i th manufacturing stage, i.e., considering the structure deposited until time T i . According to the definition of k j (ρ) in Section 3 and by using the chain rule, ∂k j (ρ) ∂φ e is defined as: At time T i , where the calculation of ∂ρ e ∂φ e and ∂φ e ∂φ e are described in Section 2.4, (22).
For the time field, the sensitivity of the objective function c with respect to design variable τ e is given by: At time T i , , we should first define G(ρ [T i ] ) which is directly related to ρ [T i ] . Let L denote the connectivity matrix between finite elements and their nodes. L is a sparse matrix with dimension of 2n v × n e , where n v and n e are the number of nodes and finite elements, respectively. The non-zero entries of L are where v is a node in the finite element grid, V e is the set of the (four) nodes of element e, g e v is the magnitude of gravity for node v assigned by element e and it is one quarter of the gravity of e. The index 2v indicates the y-component of gravity of node v. Since we assume the direction of gravity is downwards, the x-component of gravity is zero, i.e., L(2v − 1, e) = 0.
The gravity of the intermediate structure at time T i is given by Since , with respect to φ e and τ e are given by and where L(j, k) is the entry of L in the j th row and k th column.

Numerical results
To demonstrate the space-time optimization considering self-weight of intermediate structures, we setup an experiment with the same design domain and boundary conditions as shown in Fig. 3. The gravity of intermediate structures is included. The gravity of a solid element is assigned a value such that the gravity of the final structure is 1. Note that the magnitude of the external force load (F ) is also 1. Fabrication starts from the left boundary of the design domain. The number of manufacturing stages is N = 8. In Fig. 8, we present the structures and time fields optimized with four different weighting factors, i.e., α i = 0.001, 0.1, 0.4, and 0.6. From the optimized structures (second row) it can be observed that as the influence of self-weight increases, the solutions are characterized by an increased number of solid elements in the vicinity of the fixing location (left edge). To better visualize the distribution of solid elements, we vertically divide the design domain into 12 equal subdomains, and calculate the number of the solid elements within each subdomain. The histogram shown in the last row clearly confirms that more solid elements accumulate to the left.
The compliance values from these tests are summarized in Table 1. The second row corresponds to α i = 0, i.e., the objective is independent of the gravity load. The compliance of intermediate structures due to gravity is reported for each stage. As the weighting factor α i increases, take stage 4 for example, the compliance due to gravity decreases from 3.17 to 1.09. The compliance of the final structure due to the external load is reported in the second last column. As α i increases, this compliance value also increases. This increase is mild; with α i = 0.6, the compliance is increased by 4.10% (last column), from 157.17 to 163.62. This is accompanied by a significant drop in compliance due to gravity, e.g., the compliance of stage 8 decreases from 28.39 to 16.75.
Further increasing the relative weighting factor leads to convergence issues which are typical for design-dependent loads (Bruyneel and Duysinx 2005). In the limit of an infinite weighting factor, the objective is only measured for compliance due to gravity. In this case, the least compliance is obtained by not depositing any material. We observed

Process-dependent critical loads
When manufacturing a structure with a robot platform moving along the structure (see Fig. 9), the weight of the robot platform can be substantial. The location of the robot platform, and thus the load, is coupled with the construction process. This is sketched in Fig. 10 (left). The robot starts from the top left corner (p 0 ), and consecutively moves a fixed step rightwards. At each location, it can put material within the range it can reach, depicted by a circular sector for the initial location. It is assumed that the amount of material deposited by the robot from each location is the same, i.e., the fabrication speed is constant. Since each point in the design space can be reached by the robot from a few locations, its fabrication time is bounded by a lower and upper bound. The bounds, visualized in Fig. 10 (middle), are computed based on the operation radius of the robotic arm (r). The time interval assigned to respectively. An example is illustrated with a green quad in Fig. 10 (left). The element x e is reachable by the robot arm from p 5 to p 7 . Therefore the lower bound is the starting time at manufacturing location p 5 which is T 5 , and its upper bound is the ending time at manufacturing location p 7 which is T 8 .
When the robot is located at p i , the intermediate structure it fabricates shall be able to support the robot at its next location, p i+1 . Thus the compliance due to the weight of the robot at p i+1 is included in the objective function. Located at the last fabrication location p 7 , the robot will finish the complete structure. For the complete structure, the compliance is measured for an external force (F ) applied at the bottom right. This is formulated as Eq. 31, Eq. 32, Eq. 34, Eq. 35.
This formulation is largely similar to the formulation considering the self-weight presented in the previous section. The first difference lies in (52), as the robotic weight W p i+1 r is independent of the design variables. The superscript p i+1 indicates the location of the weight. The second difference is the lower and upper bounds (53). Since the load does not depend on the density, the sensitivity analysis is a simplified version from the previous section, and is omitted here.
The optimized time field and structure are shown in Fig. 10 (right) where the compliance of the final structure due to the external load is also reported. Figure 11 shows optimized results for two different robotic motions, i.e., the robot moves up and down, and two robots located at top and bottom move rightwards at the same pace. The intermediate In the objective function, α i is used to balance between the compliance of the entire structure due to the external load (F = 1) and the compliances of the intermediate structures due to the robotic weight W r = 0.5. In the above examples, α i is 0.5. A set of 8 different α i values is used to demonstrate its influence on the optimization results. The compliances are summarized in Table 2. As α i increases, i.e., the weight of the robot plays a more significant role in the objective function, the compliance values of intermediate structures associated with the robot weight naturally decrease. For instance, at Stage 7, the compliance drops from 109.00 (α i = 0.001) to 41.59 (α i = 0.01), and 28.85 (α i = 0.5). It is observed that beyond α i = 0.5 the change in compliance is relatively small. When α i increases from 0.001 to 0.5, in contrast to the rapid change in the compliance of intermediate structures due to the robot, the compliance of the entire structures due to the external load changes mildly, as can be seen from the compliance listed in the second last column, and the relative change in the last column.
The above examples are generated with 8 manufacturing stages. To demonstrate the scalability of our framework for a larger number of manufacturing stages, it is tested with 10, 12, and 16 manufacturing stages, as shown in Fig. 14. The robot motion is depicted in Fig. 10 (left) and α i = 0.5. The compliances of the final structures are listed at the bottom. It is observed that the compliance of the final structure increases along with an increasing in the number of manufacturing stages. This is due to the fact that an   The last column (diff.) indicates the difference in compliance of the final structure, compared to α i = 0.001 increasingly larger number of process-dependent loads are included in the objective, and thus effectively reduce the significance of the external load. In Fig. 15, we show convergence plots of the compliances for the entire structure and intermediate structures. This plot corresponds to the optimization with two simutaneous robots as shown at the bottom of Fig. 11. All the compliance values reduce rapidly in the first 100 iterations. The changes in the compliances become small after about 200 iterations. The bumps in these curves are caused by the increasing of the projection parameter β d in (6) for the density field. This figure demonstrates that the optimization process converges well.

Time-dependent material properties
Until this section, we have assumed that the Young's modulus of a finite element depends solely on the density (36). If the fabrication process is taken into account, the material properties may also be influenced by the time at which the material is deposited. For instance, material curing or solidification (and thus the stiffness) may be timedependent. Figure 16 illustrates a monotonic function of Young's modulus regarding time. Assume the fabrication of the entire structure is finished at timepoint T = 1. At this timepoint, the Young's modulus of an element that has been filled with material at t e ∈ [0, 1] is calculated by With the time-dependent material properties, the compliance minimization problem is updated, Eqs. 31 − 35, where the stiffness matrix K(ρ, t) is constructed with the Young's modulus interpolated using both the density and the time field (see (55)). We test this formulation on four scenarios, including two monotonically increasing functions, and two quadratic functions which open downwards. These functions serve the purpose of demonstrating the influence of time-dependent material properties on structural design and performance. In Fig. 17 each row shows optimized results corresponding to a different time interpolation function (left). Next to the function, from left to right, the optimized time fields, optimized structures colored by the time field, and optimized structures colored by the time-dependent Young's modulus. In these examples, the time range [0, 1] is divided into five equal intervals. In each interval the same amount of material volume is allowed.
The optimization with time-dependent material properties is to some extent similar to the optimization of multiple materials (e.g., Hvejsel and Lund 2011;Zuo and Saitou 2017). The difference is that here the different materials are ordered by a time variable, i.e., the moment they are produced during the manufacturing.

Parameters
The numerical optimization process involves some parameters. Table 3 summarizes the parameters which take constant values. A continuation approach is applied to projection parameters. The density projection parameter, β d , starts from 1 and is increased every 20 iterations, by an increment of 2 for the first 200 iterations, and after that by an increment of 4, until it reaches 50. The time projection parameter, β t , starts from 10 and is increased by 5 every 30 iterations, until it reaches 50.

Continuity constraints
In Section 2.3.1, the continuity constraint is relaxed by approximating the maximum function using a p-norm, followed by a smoothed Heaviside projection. Due to the approximation error, while this constraint is numerically satisfied, it may still lead to local minima in the resultant time field. We note that these local minima are not visible from the visualization, since the minimum value is very close to its neighbors, with a difference of 10 −3 . To completely eliminate local minima, we introduce an alternative formulation for the continuity constraint, with Here, # denotes the number of elements in a set. M is the set of active elements, i.e., all elements except those which are prescribed as the starting point/region for the fabrication process. N e is the set of neighboring elements. γ is a small constant which is set to 10 −9 .
As γ approaches 0, this constraint effectively restricts t e towards the mean value of its neighbors ( i∈Ne t i #(N e ) ). This alternative constraint is sufficient but not necessary, while the constraint (18) is sufficient and necessary. Figure 18  Results generated with two different continuity constraints, i.e., using (18) (left) and (59) (right) compares the optimization results using (18) on the left and (59) on the right. The time field on the right is smoother. This difference is attributed to the fact that the alternative formulation is more restrictive. This new formulation involves a quadratic term, as opposed to highly nonlinear p-norm and Heaviside projection as in (18). This is a useful alternative if smoothness is desired.
In Fig. 18 (left) we observe local maxima in the optimized time field. While these features comply with the continuity constraint, they pose some challenges for manufacturing, since it essentially requires later stages to fill some enclosed voids that have been created from previous stages. In 3D such enclosed voids are not accessible. The enclosed voids can be better detected in Figs. 12 and 13 where the full sequence is shown. As can be seen from Fig. 18 (right), the new continuity constraint effectively avoids both local minima and local maxima (i.e., enclosed voids) in optimized time fields, and thus improves manufacturability in this regard.
Both continuity constraints are defined exclusively on the time field, i.e., without considering the density field which defines the structural layout. A further investigation of Fig. 4 (left and right) reveals a potential manufacturing problem resulting from this. In Fig. 4 (second row, left), the optimized time field is monotonic, with its global minimum being located at the bottom left corner (indicated by a small blue quad) which is prescribed as the starting point of fabrication. Mapping this time field to the optimized structure, visualized in the bottom row, left, shows that the top left patch (dark red, fabricated in the first stage) is not connected to the starting point (the small blue quad) in the first construction stage. We envision a solution to this problem can be devised by defining the continuity constraint on a modified time field. Specifically, for elements which have a density value of (close to) zero, modifying their time value to 1. This modification could be realized by a (series of) differentiable projection. Figure 19 shows the convergence plot of the continuity constraint (18) for the test in Fig. 4 (right). The red curve Fig. 19 A convergence plot of the continuity constraint (18) corresponds to the aggregated value, H(t) in (18), while the blue dash line indicates the threshold , which is 10 −9 in this test. After a few oscillations at the beginning of the optimization process, the constraint is satisfied, i.e., H(t) < .

Volume fraction and design domain
We have performed tests to demonstrate that the proposed method works well with different problem settings. These tests considered self-weight and external loads (see Section 4). A weighting factor α i = 0.5 has been used to balance the compliance due to the external load and due to gravity. The continuity constraint, (59), is used. Figure 20 shows optimized designs with three different volume fractions: from left to right, V 0 = 0.3, 0.4, and 0.5. The top row shows the optimized time fields, while the bottom row shows the optimized structural layouts, colored by the corresponding time fields.
A staircase-shaped design domain is shown in Fig. 21. The design domain and boundary conditions are illustrated in Fig. 21 (left). The volume fraction is 0.6. The dimension of the design domain is 90 × 50. An external force F is applied at the top-right corner, and the gravity of intermediate structure is considered. The optimized time field and the structural layout colored by the time field are shown in the middle and right, respectively. The black polygon in the middle and right indicates the boundary of the design domain.

Extension to 3D
The proposed formulation is directly applicable to 3D. Figure 22 shows a 3D test domain, and the resulting sequence of intermediate structures. The problem formulation follows (29)-(35) in Section 3. The left face of the design domain is fixed, and fabrication is supposed to start from there. A distributed load is applied to the bottom on Space-time topology optimization on a 3D design domain. The isosurface (blue) of the optimized structure is extracted with a density threshold of 0.5. The surface that coincides with the domain boundary is indicated in red the right face. The domain is discretized into 96 × 48 × 48 cubic trilinear finite elements. The finite element analysis is performed using a geometric multigrid solver (Amir et al. 2014;Wu et al. 2016a). The time field is initialized by a distance field corresponding to the starting face. The new continuity constraint (59) is applied. β t starts from 10 and is increased by 10 every 10 iterations until it reaches 80. The volume fraction is 0.12. The filter radius, for both the time and density field, is √ 3. The other parameters take the same value as listed in Table 3.

Conclusions
In this paper, we have presented a general formulation for simultaneous design of the structural layout and the manufacturing sequence, referred to as space-time topology optimization. In addition to a density field for capturing the structural layout, a time field is introduced to encode the manufacturing process. The intermediate structures which correspond to stages of the manufacturing, are generated from these two fields. Constraints for fabrication continuity and process speed are imposed. The potential of the proposed space-time optimization is demonstrated with three fabrication considerations -self-weight of the intermediate structure, process-dependent loads due to a moving manufacturing platform, and process time dependent material properties. Clearly, these examples are by no means exhaustive with respect to the potential of the formulation. The convergence and influence of some key parameters are evaluated by an extensive parameter study.
The proposed formulation opens up a new direction in the integration of topology optimization and advanced manufacturing techniques. Extending this formulation from 2D to 3D is straightforward. As future work, we are particularly interested in considering manufacturing introduced distortion which highly depends on the manufacturing sequence.

Replication of results
Important details for replication of results have been described in the manuscript. The Matlab code is made open source, and available upon request.

Conflict of interests
The authors declare that they have no conflict of interest.
Open Access This article is distributed under the terms of the Creative Commons Attribution 4.0 International License (http:// creativecommons.org/licenses/by/4.0/), which permits unrestricted use, distribution, and reproduction in any medium, provided you give appropriate credit to the original author(s) and the source, provide a link to the Creative Commons license, and indicate if changes were made.