An Advection-Diffusion based Filter for Machinable Designs in Topology Optimization
AAn Advection-Di ff usion based Filter for Machinable Designs in TopologyOptimization Lukas C. Høghøj ∗ , Erik A. Tr¨a ff Department of Mechanical Engineering, Section for Solid Mechanics, Technical University of Denmark, Kgs. Lyngby, Denmark
Abstract
This paper introduces a simple formulation for topology optimization problems ensuring manufacturability by ma-chining. The method distinguishes itself from existing methods by using the advection-di ff usion equation with Robinboundary conditions to perform a filtering of the design variables. The proposed approach is less computationallyexpensive than the traditional methods used. Furthermore, the approach is easy to implement on unstructured meshesand in a distributed memory setting. Finally, the proposed approach can be performed with few to no continuationsteps in any system parameters. Applications are demonstrated with topology optimization on unstructured mesheswith up to 64 million elements and up to 29 milling tool directions. Keywords:
Manufacturability, Machining, Milling, Filtering, Topology Optimization
1. Introduction
Topology optimization (TO) is a widely adopted method for structural optimization [1, 2]. The main advantage ofTO is its ability to generate structures of arbitrary topology regardless of the prescribed initial conditions, i.e. little tono prior knowledge of the optimal structure is required. Due to this property, topology optimization is often used inthe initial conceptualization of new designs. If a manufacturing method, such as milling, is imposed on the designerwhen conceptualizing a part, the main advantage of topology optimization can become a disadvantage. The freedomfor TO to generate arbitrary topologies can result in structures with features, which are impossible to manufacturewith traditional machining techniques. This can be characterized as a mismatch between the constraints placed on thedesigner and the constraints used in the TO formulation.Extending the topology optimization approach to ensure manufacturability is an ongoing field of research. Awide variety of methods for topology optimization which ensures manufacturability by additive manufacturing havebeen proposed during recent years [3, 4, 5, 6, 7]. The additive manufacturing approaches usually constrain theallowed overhang of the structure in order to ensure self-support during the manufacturing process. Other types ofmanufacturing are also considered in the literature, such as Wang et al. [8], which introduces a method for ensuringtopologies that can be extruded, by utilizing the so-called PDE-filter [9] with high anisotropy.Several approaches for manufacturable TO using milling have been suggested. Gersborg and Andreasen [10]proposes a method to consider an explicitly castable or millable design by using one design variable for each row orcolumn of a structured grid. This approach is contrasted by Guest and Zhu [11], which proposes a similar method, butretains all design variables and uses cumulative summation along rows or columns as a filter to ensure that a design canbe cast or milled. The approach by cumulative summation is extended to arbitrary directions by Langelaar [12], whichmaps the densities to a structured grid aligned with the milling direction in order to perform the summation. Finally,Hur et al. [13] recently proposed a milling filter, which uses a modified variant of the advection-di ff usion equation toperform the cumulative summation for the level-set formulation, not quite unlike the proposed method. In Hur et al.[13] Dirichlet boundary conditions are used for the flow problem, which introduce the need for domain padding inorder to correctly capture non-zero physical densities on the boundary of the design domain. ∗ Corresponding author: [email protected]
Preprint submitted to Elsevier January 2021 a r X i v : . [ c s . C E ] F e b he approaches which use a cumulative sum provide a high degree of control allowing tool shapes to be taken intoaccount [12], but the cumulative summation is also computationally very expensive to evaluate. If the sum operationis precomputed, the resulting matrix will contain many non-zero values and result in very high memory usage, whichhinders application on large-scale systems. Furthermore the assembly of such a matrix is a non-trivial task, especiallyon distributed memory systems. This paper presents an approach to perform the milling tool emulation, by solvingthe advection-di ff usion equation with Robin boundary conditions. The proposed approach is conceptually simple andscales well, facilitating the solution of high resolution 3D problems. Like the cumulative summation approach, thismethod can be employed without any continuation scheme, and with little to no tuning for problem specific parameters.For now, however, the approach is limited to milling considerations without explicit tool shapes.This paper is organized as follows; Section 2 introduces the formulation necessary for the advection-di ff usionfiltering step; Section 3 introduces the optimization problem used for numerical examples; Section 4 discusses thecomputational e ffi ciency of the presented methodology; Section 5 presents two and three dimensional machinableresults, with machining in one or multiple directions and Section 6 provides a discussion of the presented methodologyand a conclusion.
2. Formulation
The milling filter is based on the approach presented by Langelaar [12], with the notable di ff erence that thecumulative sums have been replaced by solving the advection-di ff usion equation, with a dominating advective term.The filter relates the design variables to a ’physical’ density field, which guarantees that the resulting void regions arereachable by a tool direction from outside the design domain. This relation is performed through a series of filters, ofwhich most are already common in many topology optimization approaches. A brief overview of the steps used tocompose the complete milling filter is presented in Figure 1. The initial design variable is filtered in step 1, yielding theregularized design field, ˜ ρ , as discussed in Section 2.1. In step 2, the shadowing operations are performed, resulting inthe shadowed intermediate design fields, ˇ˜ ρ s , as discussed in Section 2.2. Note that a shadowing step is performed foreach prescribed tool direction. The shadowed intermediate design fields are then agglomerated in step 3, resulting inthe intermediate agglomerated variable ˇ˜ ρ , as discussed in Section 2.3. Finally, in step 4, the agglomerated variable isprojected, as discussed in Section 2.4. The sensitivity analysis of the filter is included in Section 3.2 for completeness. A density filter is initially applied to the design field. This can be performed either by a classical convolutionapproach [14], or by employing the so-called PDE-filter which solves a modified Poisson equation to perform densityfiltering [9, 15]. The reader is referred to the cited articles for details on the respective filters and their sensitivityanalysis.The density filter is introduced to regularize the design field. This is especially useful when evaluating a black andwhite field, where the gradients will be zero in most of the void region (due to the SIMP interpolation) and in partsof the solid region (due to the agglomeration of shadowing steps). By including the density filter, non-zero gradientvalues will be smoothed into the domain with zero sensitivities.In this work, the convolution approach is used for the 2D examples in Section 5.1, due to the widespread use ofthis approach for 2D codes, such as [16, 17]. The PDE-filter approach [9, 15] is used in the 3D examples, due to theincreased performance and reduced memory footprint compared to the convolution approach. ff usion The advection-di ff usion equation describes the transport of a scalar field, it is known from e.g. heat transferproblems. The steady state advection-di ff usion equation is given as: Pe u i ∂ T ∂ x i − ∂ T ∂ x i = q (1)where T is the transported field, u i the prescribed advective field, q the normalized source term and Pe the Pecletnumber which is the ratio between the advective and di ff usive transport.2: Design variable ρ
1: Density filter ˜ ρ
2: Shadowing steps ˇ˜ ρ , ˇ˜ ρ , ..3: Agglomorating fields ˇ˜ ρ
4: Heaviside projection ¯ˇ˜ ρ Figure 1: Flowchart and visualization of the composition of filters to realize the milling filter. Two tool directions are used for the milling filter,resulting in two distinct shadowed fields.
In the presented shadowing methodology, a regularized design field ˜ ρ is used as a source term, and the obtainedtransported field ˇ˜ ρ s is shadowed in the direction of the constant advection term u si of unity norm || u si || =
1, correspondingto the machining direction, where s refers to the shadowing angle index, since multiple shadowing angles may beapplied. For numerical reasons, the equation is normalized by the Peclet number. Thus, the equation solved using thefilter notation becomes u si ∂ ˇ˜ ρ s ∂ x i − Pe ∂ ˇ˜ ρ s ∂ x i = s ˜ ρ in Ω (2)where s is a factor scaling the source term on an element basis. The factor s is chosen such that the solution can gofrom 0 (void) to 1 (solid) over a single element and is further discussed in Section 2.2.2. The scaling is applied on thesource rather than on the solution, as unstructured meshes are considered.Solutions to Equation (1) are known to become numerically unstable in cases with a high Peclet number. However,while paying a price in solution accuracy, the Finite Volume formulation using an upwind di ff erence scheme on theadvection term is known to be numerically stable for high Peclet numbers. The derivations of the used finite volumeschemes are presented in Appendix A. The advection-di ff usion equation is only solved in the optimization domain. Hence, boundary conditions areimposed on free boundaries, as well as on the boundaries between the optimization- and passive solid domains - asillustrated in Figure 2.The employed boundary conditions on the free domain boundaries, Γ R , are of Robin type on boundary surfaces:ˇ˜ ρ s + n sPe ∇ ˇ˜ ρ s = Γ R (3)3 igure 2: Overview of the employed boundary conditions on the optimization domain Ω with shadowing directions θ and θ . The passive soliddomain Ω P has Dirichlet-type boundary conditions along its boundary Γ D to the optimization domain. Free boundaries of the optimization domain, Γ R , have Robin-type boundary conditions imposed. where n is the outward pointing surface normal at the boundary.On the domain boundaries adjacent to solid passive domains, Γ D , Dirichlet boundary conditions are introduced toobtain the shadow of the passive domain in the corresponding milling direction. This ensures manufacturability of thefull design, including the passive domains. ˇ˜ ρ s = Γ D (4)The implementation presented in the present work is based on unstructured meshes, hence passive void elements arenot required. If one wants to use the presented methodology with the use of passive void elements, the authors suggestto ensure that the embedded design domain is machinable with the used tool directions.Further details on the boundary conditions and their implementations are given in Appendix A. The s parameter, that is used for scaling the source term in the advection-di ff usion equation, Equation (2), is setsuch that a ˜ ρ = | u i | =
1, this parameter should be set as: s = h e (5)where h e is the average element side length.The choice of Peclet number is critical to the e ff ect of the advection-di ff usion filtering step on the obtained results.As straight walls are machinable, the conic (di ff usive) e ff ects of the filter are to be minimized. The desired e ff ect of theshadowing is hence to obtain features parallel to the specified direction, u i , with as little di ff usion as possible. This canbe achieved by using a high Peclet number. When setting Pe (cid:29)
1, the problem becomes advection dominated. Thisis also seen in Figure 3, where the field from Figure 3(a) is shadowed in a diagonal direction at three di ff erent Pecletnumbers. With Pe = .
1, as seen in Figure 3(b), the di ff usive e ff ect is clearly dominating. When raising the Pecletnumber to Pe = ff usion still has asignificant e ff ect. In the case with Pe = , seen in Figure 3(d), information is seen to almost exclusively propagatedownstream. However, inside the shadowed region, some di ff usive e ff ects are still observed. As the values inside theseregions are ˇ˜ ρ s >>
1, this is not an issue as they are thresholded using the smooth Heaviside projection, introduced inSection 2.4.It should be noted that the Peclet number is defined based on a characteristic length scale of the considered domain.This implies that when larger domains are considered, a lower Peclet number will su ffi ce to achieve similar results.For reference, Figure 3 is computed on a domain of unit side-lengths. A rule of thumb to choose the Peclet number4ased on a characteristic domain length l c is Pe thumb = l c . (6)The constant 10 is by no means a fixed rule, as acceptable results have also been found using constant factors of10 , 10 , and 10 . We do however advise ensuring that the used Peclet number satisfies Pe l c > . (a) ˜ ρ (b) Pe = . Pe =
100 (d) Pe = Figure 3: Shadowing step performed on the field (a) with di ff erent Peclet numbers. Note the di ff erent color-bars for the respective cases. The resulting fields of the shadowing for each tool direction ˇ˜ ρ s are agglomerated to a single field using the p -mean,which provides a di ff erentiable approximation to the min operator.ˇ˜ ρ e = n s n s (cid:88) s = ( ˇ˜ ρ es ) p p ≈ min s ˇ˜ ρ es (7)Here p < p -mean becomes a more accurate approximation of the minoperator, but also more non-linear, as p → −∞ . In practice p = − ff ect of a low accuracy approximation of the min operator. Whena large number of milling directions are used it becomes necessary to use a lower value of p , in order to ensure asu ffi ciently good approximation of the min operator. Other field agglomeration approaches are possible as discussedin detail in [18]. For instance, Langelaar [12] suggested using the KS functional for the agglomeration. For largeproblems, however, where ˇ˜ ρ s can take very large values, due to the high number of elements along one axis, thisresulted in numerical issues with the KS function, and the p -mean was found to be more robust in these cases. The aggregated fields are projected using the smoothened Heaviside filter [19] in order to bound the final densityfield between 0 and 1: ¯ˇ˜ ρ e = tanh( βη ) + tanh( β ( ˇ˜ ρ e − η ))tanh( βη ) + tanh( β (1 − η )) (8)where β is the projection sharpness and η the projection threshold. A desirable side e ff ect is that a black-and-whitedesign is also obtained. Throughout this article the value η = . β is chosen to either 8 or 10 dependingon the required projection sharpness.
3. Optimization formulation
The general optimization problem is formulated on the design variables ρ . ¯ˇ˜ ρ is used to denote the ’physical’ densityfield, and is obtained by performing all milling filter steps described in Section 2 on the design field ρ . The posed5ptimization problem is minimization of the compliance with a global volume constraint:min ρ ∈ R N c ( ρ ) = u (cid:62) K ( ¯ˇ˜ ρ ) u s . t . g ( ρ ) = V ( ¯ˇ˜ ρ ) V ∗ − ≤ K ( ¯ˇ˜ ρ ) u − p = ≤ ρ i ≤ i = . . . N (9)where K ( ¯ˇ˜ ρ ) is the sti ff ness matrix, p the load vector and u the resulting displacement vector. The volume constraint isenforced based on the volume of the current design, V ( ¯ˇ˜ ρ ), and the maximum allowable volume, V ∗ .The Solid Isotropic Material Penalization (SIMP) interpolation scheme is used to map the element projected designvariable ¯ˇ˜ ρ e to the corresponding Young’s modulus [14]: E ( ¯ˇ˜ ρ e ) = E min + ( ¯ˇ˜ ρ e ) p ( E max − E min ) (10)Where E min and E max are the lower and upper values of the Young’s modulus, respectively and p is the SIMP penaliza-tion parameter. The optimization problem is solved using the Method of Moving Asymptotes [20] implemented by [21] with nonstandard parameters. These parameters are necessary as the milling filter projection does not work for low projectionsharpnesses, as this results in ¯ˇ˜ ρ not being bounded correctly. This requires the usage of a constant and high β value.When optimizing with a constant projection sharpness, the problem becomes highly sensitive and the asymptotes inthe MMA algorithm hence need to be tightened, as discussed by Guest et al. [22].Therefore, to eliminate oscillatory behavior the initial asymptotes are tightened by setting the initial asympotespacing, asyinit , s = . β + . (11)Furthermore the parameter to widen the asymptotes, asyincr , is set to 1.05 and the parameter to tighten them, asydecr , to 0 .
65 instead of 1.2 and 0.7, respectively, in the standard distribution of MMA. The outer mover limit wasset to 0 . . The sensitivities of the objective- and constraint functions are obtained with respect to the final projected elementvariable, ¯ˇ˜ ρ e . The sensitivities of a function f with respect to the projected variable, ∂ f ∂ ¯ˇ˜ ρ are projected back to the designvariable, ρ , by use of the chainrule: d fd ρ = d fd ¯ˇ˜ ρ ∂ ¯ˇ˜ ρ∂ ˇ˜ ρ n s (cid:88) s = (cid:34) ∂ ˇ˜ ρ∂ ˇ˜ ρ s ∂ ˇ˜ ρ s ∂ ˜ ρ (cid:35) ∂ ˜ ρ∂ρ . (12)where the therm ∂ ˇ˜ ρ s ∂ ˜ ρ , represents the chainrule term of the shadowing step. To correct the filter for the advection-di ff usionequations the discretized milling filtering step is considered in tensor notation: A si , j ˇ˜ ρ s , j = ˜ ρ i (13)6i ff erentiating the expression with respect to the regularized field ˜ ρ i and multiplying with the sensitivities ∂ f ∂ ˇ˜ ρ s , j withrespect to the shadow yields A si , j ∂ f ∂ ˜ ρ i (cid:12)(cid:12)(cid:12)(cid:12)(cid:12) s = ∂ f ∂ ˇ˜ ρ s , j (14)which in matrix notation yields to solving the transposed filtering equation: A s (cid:124) ∂ f ∂ ˜ ρ (cid:12)(cid:12)(cid:12)(cid:12)(cid:12) s = ∂ f ∂ ˇ˜ ρ s (15)where the partial derivatives ∂ f ∂ ˜ ρ (cid:12)(cid:12)(cid:12)(cid:12) s need to be summed for all shadowing steps to obtain the full sensitivity with respectto the regularized field, as seen in Equation (12). More details on the sensitivity analysis and chain-rule terms areoutlined in Appendix B.
4. Computational E ffi ciency The computational cost of introducing the proposed milling filter depends directly on both the desired number ofmilling directions and the number of constraints used in the optimization formulation. For every milling direction anadvection-di ff usion equation needs to be solved at every design iteration. Likewise, an adjoint system of equationsneeds to be solved for every tool direction for the objective function and for every constraint every design iteration.When also accounting for the used density filter this results in a total of ( n constraints + n tools +
1) linear system solvesfor filtering, where n constraints denotes the number of constraints and n tools denotes the number of tool directions.The worst case presented in this article uses one constraint and 29 tool directions, resulting in 90 auxiliary linearsystems solutions every design iteration. While many auxiliary systems might need to be solved, the time required tosolve the finite volume problems is significantly lower than the time required to solve the linear elasticity state equation,since they are scalar problems and are hence much cheaper than the vectorial state problem. Furthermore, the systemmatrix of the advection-di ff usion problem does not change between design iterations, amortizing the computationalcost of constructing preconditioners, as these can be stored throughout the optimization process. The cost of solvingthe advection-di ff usion equation depends on a large set of parameters, where some of the most significant are; the useddiscretization of the equations, the number of elements in the mesh, and the used method to solve the resulting systemof equations.The used method for solving the resulting system of equations is also of great importance for the e ffi ciency. Ifthe system is su ffi ciently small, e.g. less than one million elements, a direct solution technique is the most viablestrategy. As the operator does not change with the design iterations, only one LU factorization is required for every tooldirection for the entirety of the optimization problem. This is the approach used in the 2D examples from Section 5.1.If the number of elements grows too large, iterative solvers are required. For the presented 3D examples theFlexible Generalized Minimal RESidual [23] method was used, with additive Schwartz preconditioning, which inturn uses incomplete LU to approximate the local solutions. All of these solvers are implemented by the Portable,Extensible Toolkit for Scientific Computation [24]. More advanced preconditioners, such as the multigrid methodemployed in solving the linear elastic state equation, can also be implemented in order to improve the computationale ffi ciency of solving the advection-di ff usion equation. This was not deemed necessary during our numerical examples,as the solution time of the auxiliary problems never exceeded reasonable limits.The discretization technique for the advection-di ff usion problem is a finite volume scheme with an upwind dif-ference scheme, which is known to be numerically stable for large Peclet numbers. Unfortunately, this means thatthe adjoint problem corresponds to a downwind di ff erence scheme, which is not so stable. It is observed that solvingthe adjoint problem requires up to 6 times the iterations before the convergence criteria is reached, compared to theforward problem for the advection-di ff usion equation. This could be mitigated by explicitly constructing the adjointoperator using a upwind scheme, although this might introduce a small error due to the di ff erence in the discretizationschemes.Using a finite volume scheme for solving the advection-di ff usion PDE means that the design field representationcorresponds to one degree of freedom per design element. Hence no mapping is required between elements and nodes.7 . Numerical examples The introduced methodology is first demonstrated on a two dimensional cantilever beam example, which was alsoused by [12]. The same cantilever beam is then optimized in three dimensions, showcasing how the methodologyworks on large scale. Finally the GE engine bracket [25] is optimized to illustrate how the methodology can be usedin an industrial example. The three dimensional examples are implemented in a preexisting in-house unstructuredtopology optimization code [26]. The optimization parameters for all three cases are given in Table 1. Note thatthe homogeneous initial value of the design vector is also stated as an optimization parameter, as the choice of thisparameter is important as the resulting initial design should neither be pure solid nor pure void. This is due to thesensitivities vanishing in SIMP at pure void and in the Heaviside projection if the agglomerated field ins ˇ˜ ρ >> Table 1: Summary of used optimization parameters
2D cantilever 3D cantilever GE BracketProjection sharpness β η . . . p
3, 5 3 3Filter radius r min .
03 0 . . Pe p -mean penalization − − − − E max . E min / E max − − → − − → − Poisson’s ratio ν . . . ρ init . .
02 0 .
002 0 . V ∗ . .
15 0 . In the following, the optimization results on a two dimensional cantilever beam are discussed. The design domainof dimensions 2 × Figure 4: Problem setup of the 2D cantilever beam example. Material can be placed inside the design domain of dimensions 2 ×
1, depicted in gray.
The domain is discretized by 200 ×
100 rectangular elements. In addition to the parameters seen in Table 1, themilling directions are selected to a variety of combinations. The SIMP penalization is set to p = p = ρ init = .
005 for cases with a single milling direction and ρ init = .
02 if multiple directions areconsidered.A reference example is optimized using the robust formulation [19]. For the reference, the β value is continuated,and the SIMP penalization power is set to p =
1. The projection thresholds are set to η dilated = . η nominal = . η eroded = .
8, for the dilated, nominal and eroded fields, respectively. This choice of threshold values should ensure aminimum length scale of 0 . r min = . (a) No milling, C nominal = .
98 (b) θ milling = C = .
68 (c) θ milling = − C = . θ milling = C = .
04 (e) θ milling = C = .
49 (f) θ milling = C = . Figure 5: Reference design and single tool direction designs for the 2D case. The final projected variable, ¯ˇ˜ ρ is shown. Results optimized with a single milling direction are shown in Figure 5. It is observed that all of the resultingdesigns are distinctly shaped by their corresponding tool direction. From Table 2 it can be seen that all of theobtained structures have a higher compliance than the reference design. The best performing design Figure 5(e) hasa compliance, which is 25% higher than the reference design. From this, it is seen that constraining the design to befiltered from a single tool direction severely limits the design freedom of the optimization algorithm.Designs optimized with multiple milling directions are shown in Figure 6. Having multiple milling directionsimpacts the obtained design, as material can be removed from multiple directions. This is most notably seen whencomparing Figures 5(b), 6(a) and 6(b). In the first case, Figure 5(b), material can only be removed from the right handside, resulting in material being placed in the lower left triangle part of the design domain. In the second case, seenin Figure 6(a), material can also be removed from below, resulting in material being placed near the upper designdomain boundary. In Table 2, it is seen that the first case, with a single direction performs better than the latter one,which indicates that a local minimum with an inferior performance has been obtained. The same local minimum wasalso observed by Langelaar [12]. The solutions find di ff erent local minima, due to the tool directions heavily a ff ectingthe initial density layout which is found using a homogeneous design variable distribution. This could potentially beavoided with an other starting guess.With a third milling direction, as seen in Figure 6(b), the obtained design can become attached to both extremitiesof the supported side. In the performance comparison, Table 2, it is also seen that this design performs better than theones presented in Figures 5(b) and 6(a) and only 24% worse than the reference design from Figure 5(a).The design obtained when optimizing with the four diagonal tool directions is seen in Figure 6(c). The design hassome features similar to the one from Figure 6(b). However, the groove going into the structure from the support planecannot be deeper with the selected tool orientations. The selected tool directions are also responsible for the triangularshape, seen inside the groove. This triangular shape does not bear any load, however it cannot be removed with theselected tool directions. The design performs 65% worse than the one seen in Figure 6(b), which is most likely due tothe high amount of non-load bearing material. 9 a) θ milling = {− , } , C = .
46 (b) θ milling = {− , , } , C = .
92 (c) θ milling = { , , , } , C = . Figure 6: Resulting designs when using multiple tool directions for the 2D case. The final projected variable, ¯ˇ˜ ρ is shown. It is noted that the performance of the designs with θ milling =
160 and θ milling = {− , , } have very similarcompliances and are also similar in a qualitative manner. It is noted that the design with θ milling = θ milling = {− , , } , seen in Figure 6(b). This might explain the slightly lower compliance of the design seen in Figure 6(b). Table 2: Comparison of compliances of the two dimensional results.
Figure θ milling C C / C ref .
98 1 . .
68 2 . .
23 3 . .
04 3 . .
49 1 . .
93 4 . {− , } .
46 3 . {− , , } .
92 1 . { , , , } .
14 1 .
65A large qualitative similarity is observed between the designs obtained by Langelaar [12] and the ones presented inFigures 5 and 6 - with the exception of the reference design (Figure 5(a)). The di ff erences in the reference designs areprobably due to di ff erent formulations being used for obtaining the reference design by [12] and the presented work. The three dimensional cantilever beam has a similar design domain as the two dimensional one from Figure 4,where the out of plane direction is modeled with the thickness 1. The applied load is a line load at the corresponding3 dimensional position. The domain is discretized with a mesh containing 31 .
25 million hexahedral elements. Theapplied filter radius corresponds to approximately 1 . β -value for the Heaviside projection sharpness is ramped up to a final value of 59 .
4, and the optimizationis process is run for 700 iterations. The threshold values are set to η = . η = . η = . . ffi cult and slow to solve due to the ill conditioning. As a remedy, a continuation scheme is appliedon the Youngs module contrast, which is reduced by a factor of 10 after 20, 40 and 60 iterations. The optimization ishence started with E min = × − E max and reaches E min = × − E max at iteration 60.10 igure 7: Reference cantilever beam obtained with no milling and a robust formulation (nominal design is shown). C nominal = . × . Several milling direction cases are considered. In the cases with a single milling direction, the p norm power isset to p = p = − p = − p norm power is changed with the varying number of tooldirections due to the changing characteristics of the p norm. A lower value of p results in a less accurate approximationof the min operator, but also results in a less non-linear operator, which is desired during the optimization.The cantilever beam is optimized using milling directions normal to the supported plane (in both directions). Theobtained design is seen in Figure 8. The obtained structure consists of two vertical plates that curve down towards theloaded line. A plate at the bottom of the domain connects the two vertical ones.In Figure 9, the design obtained with a single milling direction from top is shown. The design consists of twovertical plates, however they are not connected at the bottom of the domain. A third vertical plate connects the twoother ones along the loaded line. This third plate is higher than the plates to which it is connected. This could be toadd some bending sti ff ness along the loaded line. Figure 8: Cantilever beam with a single milling direction from the front and back. C = . × .Figure 9: Cantilever beam with a single milling direction from the top. C = . × . A cantilever beam obtained with milling only from one side, orthogonal to both the load direction and the normalof the clamped surface, is shown in Figure 10. The design cross section resembles the two dimensional cantileverbeam design with no milling constraints seen in Figure 5(a). However, this three dimensional design is not a pureextruded version of the two dimensional design, as less structure is present in the side where the milling tool comes11rom. This probably reflects that a design, where the load can be transferred to a reduced part of the domain performsbetter.A beam optimized with a similar direction, which has been rotated 45 degrees towards the supporting plane, isseen in Figure 11. The resulting cross section of the structure is seen to resemble the one from Figure 10. However,the skewed milling direction a ff ects the design drastically, both qualitatively and in terms of compliance value. It canbe observed that the angle forces material to be placed at locations that are unfavorable, as none of the other designsmake use of them (most notably the upper left side over the loaded line). Figure 10: Cantilever beam with a single milling direction from the side. C = . × .Figure 11: Cantilever beam with a single milling direction with a 45 degree angle in the plane. C = . × . The three previously shown cases with directions normal to the domain boundaries are combined to a single case,with six milling directions, where the milling directions follow the cartesian coordinate system, in positive and negativedirections. The design obtained with six milling directions is shown in Figure 12. The obtained design resembles theone obtained with no milling constraint, Figure 7. However, the hollow interior of the side walls from the reference isinfeasible with the milling filter, and has hence been replaced by holes in the structure connecting the top and bottom.
Figure 12: Cantilever beam with six milling directions, normal to the bounding box surfaces. C = . × . A cantilever beam is also optimized using a set of 29 milling directions, all coming from the northern hemisphere,as described by [12]. The design obtained with this combination of milling directions is seen in Figure 13. The designis very similar to the one using 6 directions seen in Figure 12. However, most of the holes through the structure in12he middle have been filled, and the thickness of that structure is varied instead. Some of the overhang near the onlyremaining hole next to the loaded line are only realizable due to some of the milling directions with the skewed angles.
Figure 13: Cantilever beam with 29 milling directions, coming from direction from the northern hemisphere. C = . × . The performance of the obtained designs is compared in Table 3. As expected, all designs obtained with themachinability constraints perform worse than the benchmark using the robust formulation. Furthermore, it is alsoobserved that adding more milling directions improves the performance of the design, as this also increases the designfreedom of the optimization process.
Table 3: Comparison of performance of the 3D cantliver beam designs, obtained with di ff erent milling directions and comparison with referencedesign obtained without milling. Figure Num. tools
C C / C ref . × . × . . × . . × . . × . . × . . × . ≈
11% worse, even though the obtaineddesigns are somewhat similar, qualitatively.
To demonstrate the real-world possibilities of the milling filter, the GE jet engine bracket [27] is used as anindustrial example, which is shown in Figure 14. The original bracket geometry, material properties, and load casesare adapted from [27, 25], with slight simplifications. The pinned boundary interfaces are modeled by clamping theinner surface of the bolt interface. The rigid pin is included in the model using 10 times the Young’s modulus ofthe used material. All six interfaces in the bracket model are assigned a passive solid ring, in order to ensure that aninterface exists. The compliance of the structure is minimized subject to a constraint of volume fraction V ∗ = .
137 inthe design domain, corresponding to a total bracket weight of 300 g for the design and the passive solid rings.The minimum Young’s modulus E min is updated using a simple continuation scheme, as also described in Sec-tion 5.2. The optimization begins with E min = × − E max which is decreased by an order of magnitude every 15iterations until reaching E min = × − E max at iteration 45.The large size of the computational domain (160 × × igure 14: Geometry of GE Jet Engine Bracket example. The pin with increased sti ff ness is shown in blue, while the passive rings which are keptsolid are shown in red. It should be noted that the bracket design domain is non-convex due to the two flanges which support the bolt.This non-convexity has some important implications for the solution of the advection-di ff usion equation on the bracketdomain. The boundaries will not transport information across the gaps in the design domain, treating every domainsurface as an outer surface from where a tool might remove material. This will have an e ff ect, as placing material inone flange, will not force the formulation to place material in the other flange due to the milling filter. As the flangesrepresent a small part of the domain, this error is limited, although present.As discussed in Section 2.2, the surfaces which separate the passive rings from the design domain are modeledusing a Dirichlet boundary condition in the advection-di ff usion equation in order to ensure that the passive rings alsointroduce material downstream into the design domain. Futhermore, it should be noted that while the sti ff ened pin,shown in blue in Figure 14, is included in the model, it is omitted from the visualization.The bracket examples are solved using a mesh of 64 million hexahedral elements. The results are computed on 20nodes at the DTU Sophia cluster. The computations take between 17 and 31 hours, where the milling filters and theiradjoint problems account for up to 40% of the runtime during a design iteration in the worst case with five millingdirections. While this percentage seems like a high cost of filtering, it should be noted that the total time of a designiteration at this point is 6 minutes for a design iteration, including solving four linear elastic problems with 193M dof,five advection-di ff usion equations, and ten adjoint advection-di ff usion problems. Figure 15: Reference design of the bracket example with Robust formulation and consistent boundary condition for PDE-filter. C = . × J.Figure 16: Design of the bracket example with one milling direction in the z-axis. C = . × J.
14 reference example is included for the bracket in order to compare the resulting structures and compliance values.The reference example is computed using the robust formulation [19] with β -continuation, using a PDE-based filterwith consistent boundary conditions [28], which emulate a padded domain. The robust formulation is used due to theincluded heavy-side filter, which results in discrete 0-1 designs, like the proposed milling filter. This allows for a moreaccurate comparison of compliance values, as neither method is forced to include intermediate densities in the finaldesign. Though, it should be noted that the robust formulation imposes a length-scale on the final design, which is notpresent in the milling filter. The resulting reference structure is shown in Figure 15. Figure 17: Design with five milling directions. C = . × J.Figure 18: Design with five milling directions. C = . × J.Figure 19: Design with four milling directions at an 45 degree angle. C = . × J. The first example, which is performed using only one milling direction from the top is shown in Figure 16. It canbe seen that the volume underneath the passive rings, which contain the rigid pin, is set to solid, due to the Dirichletboundary condition on the shadowing step. This prescribes the use of an already small global volume fraction, leavinglittle material for the remaining structure. All four supporting bolts are connected by thin legs, which connect in a crossshaped structure. From the resulting compliance values, shown in Table 4, it is clear that the single direction performspoorly, resulting in approximately twice the compliance value of the reference case. This is somewhat expected, giventhe restriction of available material and severely restricted design freedom.Two additional bracket examples each using five mutually orthogonal milling directions, which are visualized bythe red bars along with the resulting structures in Figures 17 and 18. The structures both di ff er from the reference bybeing more compact almost truss-like, as opposed to the shell-like reference structure. The compliance value of bothexamples with five directions are very close to the reference value, considering the design restrictions, deviating withonly 16% and 17%. From this, and similar cantilever examples with many tool directions, it can be seen that structures15 igure 20: Design with four milling directions at an 14.4 degree angle. C = . × J.Figure 21: Zoom in on the wedge-like features from the design obtained with the four milling directions at a 14.4 degree angle, seen in Figure 20. generated with milling constraints using many tool directions in general perform very well compared to referencedesigns, when taking into account the severe limitations in design freedom.Additional numerical examples are also given with more constrained tool directions, to evaluate the millingapproach itself, as structures under very limiting manufacturing constraints are generated. Figure 19 shows the bracketwith four non cartesian tool directions. Again, an agreeable structure with only 22% increase in compliance comparedto the reference is obtained. If closely examined, the structure shows two ’spikes’ in the center of the structure. Thesespikes do not carry any load, but cannot be removed by any of the given tool directions without also removing somevital load-carrying structural member. Therefore they can be compared to the phenomenon seen in Figure 6(c), wherenon-load carrying material is also present.An alternate version of the four directions is shown in Figure 20, where all four directions have a 14.4 degree angleto the vertical axis and are aligned with one of the other axes. In this case, the resulting structure is forced to use asignificant amount of material under the supporting rings in form of spikes, as can also be seen in the zoom into thecorresponding area in Figure 21. These spikes appear as there is no milling direction which allows the removal of thismaterial. Similarly, any structural member which has material removed from under it needs to have a ’wedge’-likeshape as the tools need to reach the void under the structure.Finally, a last example with a single tool direction is shown again, but this time the tool direction is at a 14.4degree angle with the vertical axis, like in the previous example, but only the direction coming from the side of thebolt flanges is considered. Again it can be seen that a thin-walled structure appears, but this time all of the wallshave the expected angle. One curious artifact of this example is that the optimizer is able to circumvent the millingfilter, and generate a hole in the bottom part of the central wall. This is understandable in a structural sense as the
Figure 22: Design with one milling directions at an 14.4 degree angle. C = . × J. able 4: Comparison of compliance of the GE jet engine bracket examples. Figure Num. tools
C C / C ref
15 0 4 . × J 116 1 9 . × J 2 . . × J 1 . . × J 1 . . × J 1 . . × J 1 . . × J 2 . ff ect occursdue to an interesting interplay between the advection-di ff usion filter and the Heaviside projection. As can be seen inFigure 3(d), the advection-di ff usion-filtered value can slowly decrease along the advection direction. This is partly dueto the numerical di ff usion due to the used upwind scheme, which is necessary for numerical stability. If the value getsbelow the Heaviside parameter η , which is 0 . η value, although the example isincluded here for completeness.
6. Conclusion
The article presents a formulation for performing topology optimization of manufacturable structures throughmilling by using a PDE-based alternative to a cumulative summation. The proposed method has been shown togenerate manufacturability by machining in a number of numerical examples in two and three dimensions, resulting intopologies where all void regions are reachable through a tool direction from outside the domain. The proposed methodhas also shown itself to scale well to larger topology optimization problems, as shown by the numerical examplesconsisting of 64 million elements. From the large-scale examples it can also be seen that the increased number ofelements allows thin shell members in the structures to be resolved, notably in the bracket. This would not be possibleon lower resolutions.The proposed method for performing the cumulative summation by solving the advection-di ff usion equation, canalso be considered a significant simplification of the mapping process [12], since it is not necessary to keep track ofmultiple meshes, and their relative orientation.Some disadvantages of the PDE-based formulation can also be noted from the numerical examples. The advection-di ff usion equation is unable to transport material through regions which have not been meshed, such as between thetwo flanges of the bracket example. When a milling direction across flanges is chosen, no coupling between the flangesoccurs.The transport problem could in principle be solved by developing special boundary conditions, which couple thefield values across the gap in the mesh. This would require a quite non-trivial e ff ort in terms of mathematical andsoftware development. An easier alternative, would be to use an auxiliary mesh of the bounding box of the domain,and solve the advection-di ff usion problems on this mesh, somewhat similar to the approach used in [12]. Solvingthe advection-di ff usion equation in the bounding box is still expected to scale better with the number of elements,compared to performing a cumulative summation.Another disadvantage is the lack of control in the tool-shape, which is possible when explicitly computing thecumulative summations [12]. A further challenge is the di ff usivity term of the advection-di ff usion equation, which isneeded for numerical stability. It can be seen from the presented results that the di ff usive terms can be kept very small.While it is di ffi cult to guarantee a tool shape based on the solution of the advection-di ff usion equation, it couldbe possible to ensure a minimum member thickness of the features in the resulting design. This could be done byproviding eroded, nominal, and dilated variations of the design field during Heaviside projection phase. These fieldscould then be used to provide a robust formulation based on [19]. This extension of the milling formulation is left tofurther work, as it considered somewhat a separate issue from the advection-di ff usion based approach.17t is found that the chosen milling directions have a large impact on the compliance of the resulting structures,especially when few milling directions are used. This is to be expected, as the milling filter poses a large restriction onthe possible structures. Acknowledgments
The authors would like to thank Assoc. Prof. Dr. Schousboe Andreasen, Assoc. Prof. Dr. Aage, and Prof. Dr.Sigmund for technical discussions, proof-reading, and supervising our PhDs.The authors acknowledge the support of the Villum Foundation for funding the aforementioned PhDs through theVillum Investigator Project InnoTop.
Appendix A. Finite Volume implementation with upwind scheme
The advection-di ff usion equation from Equation (2) is integrated over a control volume. After applying the Gaussdivergence theorem, the finite volume form of the equation over a control volume V with boundaries S is obtained[29]: (cid:90) S n i u si ˇ˜ ρ s dS − (cid:90) S Pe n i ∇ ˇ˜ ρ s dS = (cid:90) V ˜ ρ dV (A.1)Finite di ff erence schemes are used to describe the fluxes through the faces of the cell. The stencil for the respectiveface consists of the two cells adjacent to the face. The di ff usive flux through the face is described using a CenterDi ff erence Scheme (CDS): n i ∇ ˇ˜ ρ s (cid:12)(cid:12)(cid:12) f = (cid:12)(cid:12)(cid:12)(cid:12)(cid:12) dd · A | A | (cid:12)(cid:12)(cid:12)(cid:12)(cid:12) ˇ˜ ρ s (cid:12)(cid:12)(cid:12) N − ˇ˜ ρ s (cid:12)(cid:12)(cid:12) P | d | (A.2)where A is the normal area vector and d the distance between the cell center and the neighbor cell center [30]. Thesubscript f refers to the face value, P to the cell and N to the neighbor adjacent to the face.In order to guarantee numerical stability, an upwind di ff erence scheme is used for the advection coe ffi cient in thehence, the interpolated value at the face is the one from the upstream cell:if : ( n i u si ) f > ρ s (cid:12)(cid:12)(cid:12) f = ˇ˜ ρ s (cid:12)(cid:12)(cid:12) N if : ( n i u si ) f < ρ s (cid:12)(cid:12)(cid:12) f = ˇ˜ ρ s (cid:12)(cid:12)(cid:12) P (A.3)where n i is the outward pointing surface normal.The Robin boundary conditions are introduced in Equation (3), using a linear interpolation of the ˇ˜ ρ s value on thecell boundary and CDS Equation (A.2) on the corresponding gradient, the following boundary condition contributionis found for the boundary face: ˇ˜ ρ s (cid:12)(cid:12)(cid:12) N = − ˇ˜ ρ s (cid:12)(cid:12)(cid:12) P − n i sPe A d | d || A | + n i sPe A d | d || A | (A.4)Dirichlet boundary conditions are implemented adjacent to passive domains. Here, the surface contribution isgiven as: ˇ˜ ρ s = x ∈ ∂ Ω (A.5)the right hand side in these cells is corrected with: − a f (A.6)where a f is the surface contribution of the boundary face.18 ppendix B. Sensitivity analysis The sensitivity of the full milling filter can be found by applying the chain rule on all operations, as seen inEquation (12).The derivative of the compliance w.r.t. the used density field can be found for each element as [14] dcd ¯ˇ˜ ρ e = − u (cid:62) d K d ¯ˇ˜ ρ e u = − u e (cid:62) d K e d ¯ˇ˜ ρ e u e (B.1)where the e superscripts denote the element density, deformations, and local sti ff ness matrix.The partial derivative of the Heaviside projection can be found analytically for every entry of the density field as ∂ ¯ˇ˜ ρ e ∂ ˇ˜ ρ e = β (cid:16) − tanh ( β ( ˇ˜ ρ e − η )) (cid:17) tanh( βη ) + tanh( β (1 − η ))) (B.2)The partial derivative of the p-norm agglomeration for a given field i can be found analytically for every element e as ∂ ˇ˜ ρ e ∂ ˇ˜ ρ es = n n s (cid:88) s = ( ˇ˜ ρ es ) p p − n ( ˇ˜ ρ es ) p − (B.3)The final partial derivative ∂ ˜ ρ∂ρ depends on the chosen density filtering strategy. In both cases applying the partialderivative corresponds to performing the density filtering on the accumulated sensitivity dcd ˜ ρ . References [1] M. P. Bendsøe, O. O. Sigmund, Topology optimization : theory, methods, and applications, Springer, 2003.[2] O. Sigmund, K. Maute, Topology optimization approaches, Structural and Multidisciplinary Optimization 48 (2013) 1031–1055. doi: .[3] A. T. Gaynor, J. K. Guest, Topology optimization considering overhang constraints: Eliminating sacrificial support material in additivemanufacturing through design, Structural and Multidisciplinary Optimization 54 (2016) 1157–1172. doi: .[4] M. Langelaar, Topology optimization of 3D self-supporting structures for additive manufacturing, Additive Manufacturing (2016) 11.[5] Y. Mass, O. Amir, Topology optimization for additive manufacturing: Accounting for overhang limitations using a virtual skeleton, AdditiveManufacturing 18 (2017) 58–73. doi: .[6] X. Qian, Undercut and overhang angle control in topology optimization: A density gradient based integral approach: Undercut and overhangangle control in topology optimization, International Journal for Numerical Methods in Engineering 111 (2017) 247–272. doi: .[7] K. Zhang, G. Cheng, L. Xu, Topology optimization considering overhang constraint in additive manufacturing, Computers & Structures 212(2019) 86–100. doi: .[8] B. Wang, Y. Zhou, K. Tian, G. Wang, Novel implementation of extrusion constraint in topology optimization by Helmholtz-type anisotropicfilter, Structural and Multidisciplinary Optimization 62 (2020) 2091–2100. doi: .[9] B. S. Lazarov, O. Sigmund, Filters in topology optimization based on helmholtz-type di ff erential equations, International Journal forNumerical Methods in Engineering 86 (2011) 765–781. doi: .[10] A. R. Gersborg, C. S. Andreasen, An explicit parameterization for casting constraints in gradient driven topology optimization, Structural andMultidisciplinary Optimization 44 (2011) 875–881. doi: .[11] J. K. Guest, M. Zhu, Casting and milling restrictions in topology optimization via projection-based algorithms, in: Volume 3: 38th DesignAutomation Conference, Parts A and B, American Society of Mechanical Engineers, 2012, p. 913–920. doi: .[12] M. Langelaar, Topology optimization for multi-axis machining, Computer Methods in Applied Mechanics and Engineering (2019) 27.[13] D. Y. Hur, Y. Sato, T. Yamada, K. Izui, S. Nishiwaki, Level-set based topology optimization considering milling directions via fictitiousphysical model, Mechanical Engineering Journal (2020). doi: .[14] M. P. Bendsøe, O. Sigmund, Topology Optimization: Theory, Methods, and Applications, Springer, 2003.[15] B. S. Lazarov, F. Wang, O. Sigmund, Length scale and manufacturability in density-based topology optimization, Archive of AppliedMechanics 86 (2016) 189–218. doi: .[16] A. Clausen, M. Schevenels, B. S. Lazarov, O. Sigmund, E ffi cient topology optimization in matlab using 88 lines of code, Structural andMultidisciplinary Optimization 43 (2011) 1–16. doi: .[17] F. Ferrari, O. Sigmund, A new generation 99 line Matlab code for compliance topology optimization and its extension to 3D, Structural andMultidisciplinary Optimization 62 (2020) 2211–2228. doi: .[18] G. J. Kennedy, J. E. Hicken, Improved constraint-aggregation methods, Computer Methods in Applied Mechanics and Engineering 289(2015) 332–354. doi: .
19] F. Wang, B. S. Lazarov, O. Sigmund, On projection methods, convergence and robust formulations in topology optimization, Structural andMultidisciplinary Optimization 43 (2011) 767–784. doi: .[20] K. Svanberg, The method of moving asymptotes—a new method for structural optimization, International Journal for Numerical Methods inEngineering 24 (1987) 359–373. doi: .[21] N. Aage, B. S. Lazarov, Parallel framework for topology optimization using the method of moving asymptotes, Structural and MultidisciplinaryOptimization 47 (2013) 493–505. doi: .[22] J. K. Guest, A. Asadpoure, S.-H. Ha, Eliminating beta-continuation from heaviside projection and density filter algorithms, Structural andMultidisciplinary Optimization 44 (2011) 443–453. doi: .[23] Y. Saad, A Flexible Inner-Outer Preconditioned GMRES Algorithm, SIAM Journal on Scientific Computing 14 (1993) 461–469. doi: .[24] S. Balay, S. Abhyankar, M. Adams, J. Brown, P. Brune, K. Buschelman, L. Dalcin, A. Dener, V. Eijkhout, W. Gropp, D. Karpeyev, D. Kaushik,M. Knepley, D. May, L. C. McInnes, R. T. Mills, T. Munson, K. Rupp, P. Sanan, B. Smith, S. Zampini, H. Zhang, H. Zhang, PETSc Web page,2019. URL: .[25] W. Carter, D. Erno, D. Abbott, C. Bruck, G. Wilson, J. Wolfe, G. Finkhousen, A. Tepper, R. Stevens, The ge aircraft engine bracket challenge:an experiment in crowdsourcing for mechanical design concepts, in: Solid Freeform Fabrication Symposium, 2014, pp. 1402–1411.[26] E. A. Tr¨a ff , O. Sigmund, N. Aage, Topology optimization of ultra high resolution shell structures, Thin-Walled Structures 160 (2021)107349. URL: . doi: https://doi.org/10.1016/j.tws.2020.107349 .[27] Ge jet engine bracket challenge webpage, https://grabcad.com/challenges/ge-jet-engine-bracket-challenge , 2014. Ac-cessed: 2020-12-04.[28] M. Wallin, N. Ivarsson, O. Amir, D. Tortorelli, Consistent boundary conditions for pde filter regularization in topology optimization, Structuraland Multidisciplinary Optimization (2020). doi: .[29] J. H. Ferziger, M. Peri´c, Computational Methods for Fluid Dynamics, Springer Berlin Heidelberg, 2002. doi: .[30] C. B. Dilgen, S. B. Dilgen, D. R. Fuhrman, O. Sigmund, B. S. Lazarov, Topology optimization of turbulent flows, Computer Methods inApplied Mechanics and Engineering 331 (2018) 363–393. doi: ..