Topology optimization subject to additive manufacturing constraints

In topology optimization the goal is to find the ideal material distribution in a domain subject to external forces. The structure is optimal if it has the highest possible stiffness. A volume constraint ensures filigree structures, which are regulated via a Ginzburg–Landau term. During 3D printing overhangs lead to instabilities. As a remedy an additive manufacturing constraint is added to the cost functional. First order optimality conditions are derived using a formal Lagrangian approach. With an Allen-Cahn interface propagation the optimization problem is solved iteratively. At a low computational cost the additive manufacturing constraint brings about support structures, which can be fine tuned according to demands and increase stability during the printing process.

focuses on eliminating instabilities during the printing process which are caused by overhangs. These are problematic, because printing layer by layer can cause sagging of material. Avoiding these instabilities is mandatory for a reliable and reproducible printing process, which is necessary for a safe application of additive manufacturing on an industrial level.
Commonly, external support structures are added to alleviate that problem. On the downside, this method uses more material and takes longer to print. Removing support structures after the printing process is both time and labor intensive.
Topology optimization is concerned with the automatic optimal design of structures for a given physical setting, see [9]. From a mathematical point of view, topology optimization belongs to the area of optimal control of partial differential equations, [39]. Physical settings could for example include fluid dynamics, heat conduction or, as in the present setting, mechanics. Realizing the resulting designs can be challenging with traditional manufacturing methods, but is possible via additive manfucturing. For a broad overview of trends in topology optimization for additive manufacturing see [32].
Topology optimization is conducted by solving the arising partial differential equations and iteratively updating the design to improve a given objective function. Three common approaches are the SIMP method [8], the level set method [42] and the phase field method.
The problem of overhangs has been approached in different ways. As long as overhangs do not extend past an allowable self-supporting angle, they are printable. This is used in the work of [25] by incorporating projection operations to limit overhang angles. Another approach is to determine wether enough supporting material is present under a certain point. This method was developed in [30] and requires regular meshes. The problem of instabilities is here alleviated by incorporating an additive manufacturing constraint (AMC) in the modeling stage. Rigidness of the structure during the building process is taken into consideration.
When executing the topology optimization, the AMC is accounted for via a penalty term. While being printed, the resulting structure shall be stable without relying on manually added supports. This approach is introduced in the work of [2] and results are presented in [1]. In contrast to the phasefield approach used in the present paper, the authors use a level-set method, where the topology is perturbed by solving a Hamilton-Jacobi equation with a velocity term derived from the shape gradient.
As a novel approach the AMC is incorporated into the phase field method for topology optimization. Optimality conditions are derived. The approach is validated numerically by showing that the AMC increases stability during the manufacturing process.
This paper is structured in the following way: After explaining linear elasticity, the theoretical foundation of topology optimization is laid out in Sect. 2.3. The AMC is introduced in Sect. 2.4 and incorporated into the optimal control problem. First order optimality conditions are derived via a formal Lagrangian approach in Sect. 3. The discretization is explained in Sect. 4. To examine the influence of the AMC, numerical experiments are done with and without it in Sect. 5.

Problem setting
The aim of this chapter is to define the optimal control problem. After modeling the mechanics, topology optimization is explained and the AMC is formalized. First of all some notations are introduced.

Notation
Let ⊂ R d , d = 2, 3 be a bounded, regular domain and with a piecewise Lipschitz boundary . Assume there exists a Dirichlet boundary D ⊂ ∂ with surface area | D | > 0. The space of square integrable functions L 2 ( ) is used to define the Sobolev spaces For more details see [39]. The Frobenius inner product for the matrices M, N is defined by the pairwise sum of element-products With the fourth order tensor C, the product CM is defined via

Linear elasticity problem
Using the displacement u : → R d , the linearized strain tensor is defined. The distribution of material in is described by a phase field ϕ ∈ L ∞ ( ) with almost everywhere in . Here, ϕ = 0 describes void and ϕ = 1 represents areas containing material. In a physically accurate setting each point in space either does or does not contain material, i.e. ϕ ∈ {0, 1}, leading to a sharp transition. However, in the realm of optimization a smooth transition between material and void is desired in order to calculate derivatives. This is achieved by explicitly allowing impure phases, i.e. states with 0 < ϕ < 1.
The transition is seen as the interface between material and void. Consider the symmetric fourth order stiffness tensor C(ϕ) with continuously differentiable entries. It is assumed that the derivative of the stiffness tensor C (ϕ) is globally Lipschitz continuous. A transition function with cubic behavior for ϕ ∈ [0, 1] is employed. This is not essential in the phase field approach, but empirically leads to a faster convergence.
For ζ > 0 the transition function is defined via where s r is a monotone C 1,1 -function such that s is in C 1,1 . The elasticity tensor for the whole domain is defined in the following way where for a quadratic matrix E. The constants λ mat ≥ 0 and μ mat > 0 are called Lamé parameters. The constant > 0 is a small parameter, which will be related later to the interface width when defining the Ginzburg-Landau term in Equation (8).
A load f is acting on a part of the boundary labeled f and a volume force g is present where material ϕ is located. The above definition warrants 0 = C void C mat , which ensures low stiffness in void, but avoids degeneracy.
With the outer normal vector n, the strong form of the mechanical system is defined via The weak formulation is written as The superscript m is used throughout this paper in connection with the mechanical system. A unique u m exists, as proven for the phase field method in [11]. A more general and thorough discussion on existence theory for elasticity problems can be found in [20].

Topology optimization
This section is based on [14]. The aim of topology optimization is to distribute a limited amount of material in an optimal sense. Without limiting the material usage, the solution of the optimal control problem defined in Sect. 2.5 would be trivial: ϕ ≡ 1 on , i.e. covering the whole domain with material.
An additional constraint makes the problem more interesting. Let m be the mass parameter with 0 < m < 1. The amount of material is restricted via a volume constraint where | | denotes the Lebesque measure of the domain . The admissible set is defined as The objective is to find a material distribution ϕ ∈ G m and a corresponding solution of the elasticity problem u : → R d such that the mean compliance is minimized. If C(ϕ) would be allowed to become zero for ϕ = 0, boundary forces could be avoided trivially by not placing material where these forces act. On the other hand, soft surrogate material close to f leads to large displacements, increasing the compliance. This ensures that optimal solutions always have material placed where boundary forces act. However, this minimization problem is not well-posed as explained in [3]. The regularity of the solution is not ensured.
In computational examples this can lead to a checkerboard solution. Checkerboarding is the frequent occurrence of jumps between material and void, which is not desirable, see [37]. The algorithm might produce so-called porous solutions, which can be thought of as sponge-like microstructures. In many industrial cases this type of material is hard to realize and therefore undesirable. The ill-posedness can be alleviated by adding a perimeter regularization which was proven in [7]. The paper [38] explains that the perimeter can be approximated by the Ginzburg-Landau term for > 0 and ϕ ∈ H 1 ( , R) ∩ L ∞ ( ). For convergence properties as the interface parameter approaches zero see [40]. The first term penalizes transitions between material and void through the gradient of the material distribution. The second term contains a potential ψ ∈ C 1,1 (R, R) with ψ ≥ 0. A commonly used potential is the double well potential The aim of the potential is to penalize impure phases, which are states where the phase field ϕ is not equal to 0 or 1.

Lemma 1
There exist positive constants¯ ,¯ and such that for all symmetrical M, N ∈ R d×d \ {0} and all ϕ ∈ R the following relationships hold: Proof Follows from construction of C(ϕ).

Additive manufacturing constraint
In additive manufacturing a structure is created in a layer-by-layer approach. Following With a layer thickness l > 0, the layer from height hl to height h is defined via The Dirichlet boundary for the AMC problems differs from the Dirichlet boundary of the mechanical problem. In the mechanical problem it describes areas where the structure is supported when forces are applied. In the constraint problems the process of additive manufacturing is simulated. Homogeneous Dirichlet boundary conditions are assumed for all parts of the geometry touching the building plate 0 . Notice that 0 does not change for different intermediate shapes. During the manufacturing process no surface forces are applied to the structure. Each intermediate shape h is only subjected to gravity g. It is assumed that gravity is the primary opponent to successful 3D printing, which is the case for example for fused deposition modeling. Other issues, like thermal problems in powder bed fusion [29], are not considered here. As opposed to [2], where gravity acts on the whole intermediate structure, it is here restricted to the newest, upmost layer. Otherwise the compliances of the lower layers of the structure would be taken into account multiple times, whereas the top layer is only considered once. This would lead to designs with a heavy bias towards stability during the beginning of the printing process. Support structures would mainly arise in the lower parts of the design. Restricting gravity to the upmost layer can also be motivated physically: Since the lower layers of the structure are already cooled down and hardened, they will not be displaced much under the influence of gravity. Only applying gravity to the newly printed layer can be thought of as implicity considering the phase change of molten to hardened plastic filament. For a given height h > 0 the AMC system is defined as Note that one cannot just define the system on the layer L h , since the Dirichlet boundary is not part of that domain. The weak formulation is written as is a weak solution of the AMC system (AMCS), if it fulfills the weak formulation (11) on the corresponding space h .
The mechanical system (MS) describes the load case where a force is acting on the structure. On the other hand the AMC system (AMCS) models the printing process, where only gravity is considered. To help distinguish between both systems the superscript c is used in conjunction with the constraint problems.
In order to judge the stiffness during the manufacturing process, the compliance of the layer L h is defined For large overhang angles gravity will cause large displacements and therefore increase the compliance. Since compliance is the inverse of stiffness, low compliance values are desirable, because they coincide with high stiffness during the printing process. A penalty function is defined with H > 0 via The idea is to slice the domain into N layers and define successively growing domains A test geometry is created to show the effect of the weighting term 1 h . In Fig. 1(a) the material is represented by the red areas. Two vertical beams are connected via eight horizontal beams, which have large overhangs. The bottom beam touches the building plate and is therefore not critical. Each of the other horizontal beams would lead to problems during additive manfucturing. However, Fig. 1(b) shows that without weighting, the AMC penalty term P d increases for higher layers. This is caused by a linearly growing amount of material beneath the current layer, which is also subject to displacement. The algorithm would focus on decreasing overhangs towards the top of the domain. With weighting the AMC penalty term is more evenly distributed independent of layer height. As desired, the algorithm would decrease overhangs everywhere in the domain. Up to this point, the height h has been a continuous parameter. To avoid double subscripts it is now used as an index of the domain. For example, u c h is the displacement of the structure up to the h-th layer. The layer L h corresponds to h \ h-1 . The following notation is used: The outer integral of P(u c , ϕ) is approximated via a Riemann sum. For the discretized penalty function the constant factor H N is dropped, yielding Not only the final structure shall bend as little as possible, but also all intermediate structures. Towards calculating the penalty function an additional N elasticity problems have to be solved.

Formal derivation of optimality conditions
The aim of this chapter is to derive optimality conditions. A rigorous mathematical analysis was conducted in [11], where the AMC could be incorporated. Here, only the formal Lagrange approach is covered, see [39].
With the weak formulation of the mechanical system (5) and the weak forms of the AMC systems (11) the Lagrange functional can be defined via The Gâteaux derivatives of the Lagrange function with respect to the states u m , u c h , h = 1, . . . , N in the directions q m , q h , h = 1, . . . , N respectively are calculated. This results in the adjoint equations for h = 1, . . . , N . Note that these adjoint equations match the state equations (5), (11) for The derivative of the Lagrange function with respect to ϕ in direction ω is Since G m is convex, any optimal control satisfies the corresponding variational inequality, see [39, Thm. 1.2.]. The results are summarized in the following theorem.

Discretization
Observations and theoretical derivations of linear elasticity take place in the threedimensional space. However, two-dimensional models are computationally less demanding and will therefore be studied here. An Allen-Cahn interface propagation is explained, which allows stationary problems to be solved iteratively. The Primal-Dual Active-Set Method is presented to solve the optimal control problem.
The goal is to construct a thin, 3D-printable structure. The plane stress model is feasible since it assumes that the z-dimension is very small in comparison to the other dimensions.
The Lamé parameters have to be adjusted, as explained in [22].
The goal is to iteratively solve the stationary problem. The idea of a phase field interface propagations stems from the work of [19] and [4]. After an initial material distribution ϕ 0 is set, pseudo timestepping is used to evolve the solution. The negative derivative of the reduced cost functional is chosen as the direction of interface movement, thus where τ is the length of a pseudo time-step. Thus Towards better stability the equation is discretized in time using a semi-implicit approach Note that the time-step τ is just a numerical construct and does not represent a physically relevant time, hence the word pseudo. When formulating the optimality system the box constraint 0 ≤ ϕ ≤ 1 was taken into account via the space G m . In accordance with [39] Karush-Kuhn-Tucker multipliers μ + , μ -∈ L p ( ) with 1 < p 2 are introduced such that 0 ≤ ϕ ≤ 1 a.e. in , (μ + , ϕ -1) L p = (μ -, -ϕ) L p = 0 in .
In an optimum the following holds:
If 0 < ϕ < 1, then μ -= μ + = 0. Otherwise, if The regularity of μ -, μ + follows from standard arguments by making use of the regularity of ∂L ∂ϕ , see [26]. This incorporation of the box constraint via Karush-Kuhn-Tuckermultipliers leads to the Primal-Dual Active-Set Method, see [14]. In each step the active sets are updated using the result from the previous iteration. If there are no further changes in these sets, the algorithm has converged and terminates. The pseudocode is displayed in Algorithm 1.
Intuitively one can think of material flowing to the areas where it has the largest impact on decreasing the Lagrange functional. Notice that ϕ is not constrained per se and thus might end up being significantly smaller than 0 or larger than 1 in parts of the domain. The idea behind this algorithm is to enforce the box constraints. Start with a ϕ 0 that does fulfill all constraints and μ 0 + , μ 0 -≡ 0, leading to A 1 + and A 1 -being empty and thus μ 1 + , μ 1 -≡ 0. Eventually ϕ k will exceed 1 or fall below 0 in parts of the domain. Then the corresponding active set A k+1 + or A k+1 -will be nonempty and in turn ϕ set to 1 or 0, respectively. This ensures that the box constraints are met.
The potential ψ introduced in Sect. 2.3 has a local maximum in ϕ = 0.5. As an initialization ϕ 0 ≡ 0.5 is used throughout our numerical experiments. In [14] it is explained why numerical oscillations can occurr if the active set parameter c is chosen too small with respect to the mesh resolution. In our numerical experiments we did not observe oscillations when setting c = 2 h 2 min , with h min being the minimum mesh size. The elasticity equations and the gradient equation are solved via the Finite Element Method. According to the work of [27] the Primal-Dual Active-Set-Method can be interpreted as a semismooth Newton method. This yields local superlinear convergence.

Numerical examples
The goal of this section is to show the effect of the AMC on the optimal topologies. Standard topology optimization examples of an MBB beam and a cantilever were chosen. The code is implemented using P1 finite elements in FEniCS, see [33] and [5]. The main computational challenge arises from a high number of layers leading to many intermediate elasticity problems that need to be solved in each iteration. If done in serial, this is quite time consuming. Here, multiprocessing is used to solve the different elasticity problems concurrently on different processors. This leads to a considerable speedup.
Here, only linearly elastic, homogeneous, isotropic materials are considered. For a formal definition see [17]. Without explicitly mentioning them, SI-units are used throughout. Young's modulus E is set to 3.5×10 6 and the Poisson ratio ν is assumed to be 0.3. According to [35], this corresponds to the widely used 3D printing filament PLA. Computations are done for a surface force f = (0, -1) T and the maximum number of iteration steps is set to 100. The results of the topology optimization depend on the following three parameters: the Ginzburg-Landau parameter , the Ginzburg-Landau term prefactor γ and AMC prefactor β.

Parameter study for and γ
The findings of this section do not take the AMC into account. This is achieved by setting β = 0. The parameters , γ can be varied to influence the topologies resulting from the optimization. The parameter γ is the prefactor of the Ginzburg-Landau term and Algorithm 1 Primal-Dual Active-Set-Method 1: Initialize ϕ 0 and set μ 0 + , μ 0 -≡ 0, k = 0.

3: Set
s.t. ϕ k+1 dx = m for (ϕ k+1 , λ k+1 ). thus controls the influence of the perimeter regularization, which can be interpreted as a surface tension. Larger values of γ cause the material to lump together to simple, nonfiligree structures. A sensibly small γ allows for the creation of more filigree structures while still avoiding checkerboarding. The parameter regulates the influence of the two summands within the Ginzburg-Landau term. When is increased, the gradient term receives a higher weight. Penalizing the derivative causes less jumps to occur. However, impure phases with 0 < ϕ < 1 are more likely to appear. Often, is referred to as the interface thickness parameter. If it is set too large, filigree structures cannot emerge. When is decreased, a transition from material to void is penalized less. A smaller interface thickness brings about more filigree structures. If is very small, the potential term dominates. Fine, pure-phased, porous structures arise. It helps to decrease the mesh size, which allows for thinner beams to be displayed. These effects of the parameters and γ are summarized in Fig. 2.
Finding a suitable pairing can be very time consuming. In order to alleviate this problem, the parameters and γ are chosen adaptively as suggested in [23]. The following turned out to be a good choice: with c GL set to 0.1.

Cantilever beam
The setting of a cantilever beam is shown in Fig. 3. The structure is fixed at the left hand side and a downward force is applied in the middle of the right hand side. The domain is discretized using 26,667 degrees of freedom. The material usage is fixed at 40%. Ten layers In order to test the influence of AMC, calculations are done for different factors β. The results can be seen in Fig. 4. Blue areas represent void and red areas correspond to material. The white parts are the interface, where 0 < ϕ < 1. Figure 4(a) depicts the structure without influence of the AMC, i.e. β = 0. The material is concentrated in major beams. Not a lot of support structures are present to hold the major beams in place. Generally, the severity of an overhang can be measured by its angle to the normal axis. As a rule of thumb, an angle of over 45 • is critical. Here diagonal structures can be seen that have an angle to the vertical axis of well above 45 • , making it difficult to print the structure.
The penalty term is introduced by successively increasing the factor β. A change in topology can be seen when β reaches 0.01, see Fig. 4(b). Supporting structures arise, which hold up the major beams while they are being printed. This leads to less deformations during additive manufacturing. The small beams are mostly vertically oriented and thus can be manufactured easily.
As the AMC parameter β is increased, more and more support structures arise, which can be observed in Figs. 4(e) and 4(f ). More material is placed near the bottom of the domain. Beams are oriented more vertically than before and intersect each other. Smaller beams deform less, which makes the structure more stable.
Notice that opposed to the geometric approach in [25], overhangs are tackled indirectly and thus are not fully avoided.
The resulting values for the compliance, Ginzburg-Landau term and AMC penalty term can be found in Table 1. As one would expect the compliance increases as the influence

Figure 5
Setting of an MBB beam of the AMC gets larger. However, for the largest β value it only increases by around 10%, whereas the AMC penalty term is decreased by 80%, making the printing process more stable.

MBB beam
The setting of an MBB beam can be seen in Fig. 5. The structure is fixed in both directions at the bottom left hand side, but allowed to move in the horizontal direction at the bottom right hand side. A vertical force is applied in the middle of the top of the domain. Making use of the symmetry axis in the middle of the MBB beam, only half of the domain is discretized. The 400 × 133 crossed P1 finite elements lead to a total of 106,934 degrees of freedom. Different from the cantilever example, the number of layers is increased from 10 to 100, making multiprocessing even more valueable. In [2] the computation for the same number of layers using 300 × 100 Q1 elements, corresponding to 30,401 degrees of freedom, is reported to take 237 hours. With more than three times the degrees of freedom, computations here take less than 3 hours. When 30,401 degrees of freedom are used, this time reduces to 24 minutes. The amount of material is set to 20%. As before, β is increased to see the influence of the AMC. The results can be seen in Fig. 6. Without the AMC large overhangs arise. Especially the long horizontal beam at the top of the structure is difficult to print. It is not supported from below, likely causing molten material to sag during the manufacturing processing. As β is increased, a lot of smaller supporting structures are created. In Fig. 6(f ) the area below the long horizontal beam is evenly filled with supports. One also notices that the thick 45 • beam near the force boundary is becoming thinner and more vertical as β is increased. More support structures and vertically aligned beams make it easier to manufacture the design. In Table 2 the results are gathered. As β is increased, the compliance is increased by less than 10%, but the AMC penalty term decreases by almost 90%. While the structure with β = 64 deforms slightly more during its original use case, it is considerablly more stable during the printing process.

Conclusion
Instabilities during additive manufacturing motivated research on this topic. The AMC was incorporated into the optimal control problem. Optimality conditions arose from a formal Lagrangian approach. The numerical approach was discussed. Experiments showed that the AMC decreases overhangs and brings about support structures.
To the author's best knowledge this is the first time the AMC is applied to the phase field method for topology optimization. First order optimality conditions were derived. An Allen-Cahn interface propagation iteratively solved the stationary optimality problem. After a pseudo time discretization the Primal Dual Active Set Algorithm was executed. It was discussed how the topology can be influenced by setting the parameters and γ . After applying the AMC noticeable changes occured. More support structures emerged, leading to smaller intermediate compliances and in turn to a higher stiffness during the manufacturing process. The influence of the penalty term is adjustable, which allows the structures to be fine-tuned according to demands. Less labor intensive post processing is necessary, giving this approach an advantage over manually added support structures. However, opposed to the geometric approach, large overhang angles are tackled indirectly and thus are not fully avoided. It would be beneficial to print different structures to find out how this influences the manufacturing process. Further research on this topic could include the incorporation of a micro field to represent lattice structures. This helps create printable light-weight constructions. Here, a phase flow approach was used. An implementation of the projected gradient method, see [15], could speed up the convergence. Additionally, one could implement an adaptive mesh, which is fine in the interface region and coarser otherwise. This would reduce computational costs and make three dimensional calculations less demanding. An integration of stress constraints is desirable for industrial applications. On a larger scope there are many more problems in additive manufacturing which can be tackled with optimization methods.