Manufacturing cost constrained topology optimization for additive manufacturing

Abstract

This paper presents a manufacturing cost constrained topology optimization algorithm considering the laser powder bed additive manufacturing process. Topology optimization for additive manufacturing was recently extensively studied, and many related topics have been addressed. However, metal additive manufacturing is an expensive process, and the high manufacturing cost severely hinders the widespread use of this technology. Therefore, the proposed algorithm in this research would provide an opportunity to balance the manufacturing cost while pursuing the superior structural performance through topology optimization. Technically, the additive manufacturing cost model for laser powder bed-based process is established in this paper and real data is collected to support this model. Then, this cost model is transformed into a level set function-based expression, which is integrated into the level set topology optimization problem as a constraint. Therefore, by properly developing the sensitivity result, the metallic additive manufacturing part can be optimized with strictly constrained manufacturing cost. Effectiveness of the proposed algorithm is proved by numerical design examples.

Key words： topology optimization; manufacturing cost; additive manufacturing; powder bed

Introduction

Topology optimization is a mathematical method to optimally distribute materials within a given design domain subject to the physical boundary conditions [¹]. Its key advantage is in exploring a larger design space (compared to size and shape optimization) with faster convergence (compared to global optimization). The past few years have witnessed the rapid expansion of topology optimization in covering increasingly more multidisciplinary structural design problems [²,³]. On the other hand, even though designers have greatly enjoyed the creative design power of topology optimization, manufacturing engineers often take issues with the organic shapes produced by topology optimization being difficult or even impossible to manufacture, even though many efforts have been put into this aspect [⁴]. These days, additive manufacturing (AM) technology emerges and has partially relieved many of the manufacturability issues since the layer-by-layer material deposition process of AM makes manufacturing of any complex geometry much cheaper and faster [⁵]. Hence, topology optimization for AM has recently become very popular and has attracted attention from both academia and industry. Many tough topology optimization for AM issues have been addressed including self-support design to eliminate need of the costly support structure [⁶–¹¹], material anisotropy issue to address the tool path or build direction-induced anisotropic material properties [¹²–¹⁶], multi-scale design to explore the extreme design space enabled by AM [¹⁷–²¹], and many others [²²]. These achievements lead to a better marriage between topology optimization and AM.

On the other hand, AM so far is still an expensive technique compared with conventional subtractive machining or casting for most parts. Based on our experience collaborating with the industry, the authors have received many requests to take manufacturing cost into consideration when performing topology optimization. However, no such research has been done even though this topic is very useful to the industry. Hence, the main motivation of this research is to develop a topology optimization algorithm that has the manufacturing cost as a constraint with regard to the powder bed-based metal AM process.

So far, we have seen some works on AM cost estimation. For instance, Ruffo et al. [²³] developed the cost model for low-to-medium volume laser sintering and a saw tooth shaped curve is derived of the cost per part on production volume. Then, Ruffo and Hague [²⁴] polished their model by involving the situation of mixed part production. Baumers et al. [²⁵] developed the estimation models of build time, cost, and energy consumption of the EOS M270 Direct Melting Laser Sintering (DMLS) system, wherein the full platform fill is emphasized to save average cost. Beyond only estimating the cost, some other works have utilized the cost models to instruct design for AM. For example, Huang et al. [²⁶] developed a general cost model for metal AM and innovatively used it to configure the cost minimization topology optimization subject to structural compliance constraint. Barclift et al. [²⁷] realized CAD-based cost estimation for powder-bed-based metal AM through API programming and used this tool to optimize the part build orientation.

Even though AM cost-constrained topology optimization is not totally new, there lacks a thorough and in-depth study. For instance, the part height, support volume, and build batch size all significantly affect the AM cost. They have been included in different cost models but have not been carefully addressed for topology algorithm. The reason may lie in the difficulties in deriving the relevant sensitivity results; e.g., it is non-trivial to calculate sensitivity on the part height variable. Therefore, the main contribution of this paper is to develop a cost estimation model for EOS M290 DMLS system, integrate it as a constraint of the topology optimization problem, and perform optimization to explore features of the cost-constrained topological design.

AM cost model

The cost model for powder bed-based metal AM is composed of three parts: Material cost, argon gas cost, and operation cost (including labor and utility):

(1)C_{\mathrm{AM}}=C_{\mathrm{material}}+C_{\mathrm{argon}}+C_{\mathrm{operation}}.

Among the three terms, material cost is calculated by summarizing material consumption of both the part and the support, as:

(2)C_{\mathrm{material}}=\left(V_{\mathrm{part}}\rho +V_{\mathrm{support}}\overline{\rho }\right)C_{\mathrm{material}}^{\mathrm{unit}},

where, \rho is material density and \bar{\rho } is material density of the support, which is usually around 40% of \rho because of the porous infill, C_{\mathrm{material}}^{\mathrm{unit}} indicates unit price of the metal powder.

The argon gas cost is calculated by

(3)C_{\mathrm{argon}}=\left({\displaystyle \frac{\mathrm{0.6}{T}_{\mathrm{R}}h}{L_{\mathrm{t}}}}+\mathrm{3}\right)C_{\mathrm{argon}}^{\mathrm{unit}},

where 0.6 is the average argon consumption per hour in m³, T_{\mathrm{R}} is the recoater cycle time, h is the build height, and L_{\mathrm{t}} is the layer thickness, and C_{\mathrm{argon}}^{\mathrm{unit}} indicates unit price of the argon gas. A value of 3 m³ is the volume required for the initial machine purge of argon.

Then, the operation cost is calculated using Eq. (4), which is calculated by multiplying the operation time with the labor and utility cost:

(4)\{\begin{array}{l}{C}_{\mathrm{operation}}=T_{\mathrm{operation}}(C_{\mathrm{labor}}^{\mathrm{unit}}+C_{\mathrm{utility}}^{\mathrm{unit}})\\ T_{\mathrm{operation}}=T_{\mathrm{setup}}+{\displaystyle \frac{T_{\mathrm{R}}h}{L_{\mathrm{t}}}}+\left(V_{\mathrm{part}}+{\displaystyle \frac{V_{\mathrm{support}}\overline{\rho }}{\rho }}\right)/S_{\mathrm{rate}}\end{array},

where

T_{\mathrm{setup}}

means the average machine setup time and S_{\mathrm{rate}} is the volume scanning rate.

By analyzing Eq. (2)–(4), the terms V_{\mathrm{part}}, V_{\mathrm{support}}, and h are identified as design-dependent variables, and if the design problem maintains a consistent material volume fraction, the designable variables will be further reduced to only V_{\mathrm{support}} and h. Therefore, the cost model can be re-formulated into

(5)C_{\mathrm{AM}}=a_{\mathrm{0}}+a_{\mathrm{1}}{V}_{\mathrm{support}}+a_{\mathrm{2}}h,

where

\{\begin{array}{l}{a}_{\mathrm{0}}=V_{\mathrm{part}}\rho C_{\mathrm{material}}^{\mathrm{unit}}+(T_{\mathrm{setup}}+V_{\mathrm{part}}/S_{\mathrm{rate}})(C_{\mathrm{labor}}^{\mathrm{unit}}+C_{\mathrm{utility}}^{\mathrm{unit}})+\mathrm{3}{C}_{\mathrm{argon}}^{\mathrm{unit}}\\ a_{\mathrm{1}}=\overline{\rho }{C}_{\mathrm{material}}^{\mathrm{unit}}+(C_{\mathrm{labor}}^{\mathrm{unit}}+C_{\mathrm{utility}}^{\mathrm{unit}})\overline{\rho }/(\rho S_{\mathrm{rate}})\\ a_{\mathrm{2}}=(\mathrm{0.6}{T}_R/L_t)C_{\mathrm{argon}}^{\mathrm{unit}}+(C_{\mathrm{labor}}^{\mathrm{unit}}+C_{\mathrm{utility}}^{\mathrm{unit}})T_R/L_t\end{array},

and a_{\mathrm{0}}, a_{\mathrm{1}}, and a_{\mathrm{2}} are the cost coefficients.

Level set method and problem formulation

Basic introduction to level set method

Briefly speaking, a level set function, \mathrm{\Phi }(X):R^n\to R, represents any structure in the implicit form, as

(6)

where \mathrm{\Omega } represents the material domain, D indicates the entire design domain, and thus D/\mathrm{\Omega } represents the void, X indicates the coordinate vector.

Then, the Heaviside and Dirac delta functions are adopted to formulate the numerical region and boundary integrations, respectively; see Eqs. (7) and (8).

(7)

(8)\delta (\mathrm{\Phi })={\displaystyle \frac{\partial H(\mathrm{\Phi })}{\partial \mathrm{\Phi }}},{\displaystyle {\int }_{-\infty }^{\infty }\delta (\mathrm{\Phi })\mathrm{d}\mathrm{\Phi }}=\mathrm{1.}

In addition, level set method provides an effective technique to track and update the boundary profile, i.e., with the properly calculated boundary velocities, the Hamilton-Jacobi equation can be solved to effectively and efficiently update the boundary profile. For structural optimization problems, the boundary velocities are calculated through the sensitivity analysis.

Another important characteristic of the level set function is the signed distance regulation, or in other words, Eq. (9) is satisfied (at the most area) wherein the absolute level set value at any point represents its shortest distance to the structural boundary. The signed distance characteristic enhances numerical stability and also provides a powerful tool for flexible geometric control.

(9)\left|\nabla \mathrm{\Phi }(X)\right|=\mathrm{1.}

Note that, the level set function discussed earlier is based on the discrete definition which is adopted in this work. Parametric level set functions are also available to support structural topology optimization, which however, demonstrate very different features in numerical details [²⁸–³¹].

Level set representation of the cost models

Given that the structural optimization will be performed with the level set method, the AM cost model have to be re-formulated by transforming the design-dependent variables into functions on the level set function.

Specifically, Eq. (5) is re-formulated into Eq. (10) by assuming the build direction aligning the z-axis.

(10)C_{\mathrm{AM}}=a_{\mathrm{0}}+a_{\mathrm{1}}{V}_{\mathrm{support}}(\mathrm{\Phi })+a_{\mathrm{2}}h(\mathrm{\Phi }).

Here, we would mention that it is trivial to numerically calculate the support volume if the part is represented by the level set function. The following algorithm can be used:

1) Have the signed-distance level set field \mathrm{\Phi }(X) as the input;

2) Identify the support-needed overhang areas with the following equations: Find the solid element with center coordinate (x,y,z), which satisfies:

\mathit{n}\cdot \mathit{B}\le -\cos \left(\overline{\theta }\right),

\mathit{n}=-{\displaystyle \frac{\nabla \mathrm{\Phi }(x,y,z)}{\left|\nabla \mathrm{\Phi }(x,y,z)\right|}},

where, n is the unit normal vector at the structural boundary and \mathit{B} is a unit vector in the build direction, and \bar{\theta } is the threshold overhang angle and the overhang areas with the inclination angle smaller than that would need support structures.

3) For each support-needed overhang element, the related support elements can be simply counted in the -\mathit{B} direction till encountering the next solid element.

Optimization problem formulation

Inserting the above AM cost model into the conventional compliance-minimization topology optimization problem, we formulate a new multi-objective problem, as shown below:

(11)\begin{array}{l}\min J=w_{\mathrm{1}}{\displaystyle \frac{{\displaystyle {\int }_D\mathit{D}\mathit{e}(\mathit{u})\mathit{e}(\mathit{u})H(\mathrm{\Phi })\mathrm{d}\Omega }}{Com{p}^{*}}}+(\mathrm{1}-w_{\mathrm{1}}){\displaystyle \frac{C_{\mathrm{AM}}}{C_{\mathrm{AM}}^{*}}}\\ s\mathrm{.t}.a(\mathit{u},\mathit{v},\mathrm{\Phi })=l(\mathit{v}),\forall \mathit{v}\in U_{ad}\\ {\displaystyle {\int }_{D}H(\mathrm{\Phi })\mathrm{d}\mathrm{\Omega }}=V_{\max }\\ a(\mathit{u},\mathit{v},\mathrm{\Phi })={\displaystyle {\int }_D\mathit{D}\mathit{e}(\mathit{u})\mathit{e}(\mathit{v})H(\mathrm{\Phi })\mathrm{d}\mathrm{\Omega }}\\ l(\mathit{v})={\displaystyle {\int }_{\mathrm{\Gamma }}\tau \cdot \mathit{v}\mathrm{dS}}\end{array}

in which a(\cdot ) is the energy bilinear form, l(\cdot ) is the load linear form, \mathit{u} is the deformation vector, \mathit{v} is the test vector, and \mathit{e}(\mathit{u}) is the strain, U_{\mathrm{ad}}=\left\{v\in H^{\mathrm{1}}{(\mathrm{\Omega })}^{\mathrm{d}}|\mathit{v}=\mathrm{0}\begin{array}{l}\end{array}\mathrm{on}\begin{array}{l}{\mathrm{\Gamma }}_{\mathrm{D}}\end{array}\right\} is the space of kinematically admissible displacement field, \mathit{\tau } is the boundary traction force, and the force loading boundary area is considered non-designable, V_{\max } is the upper bound of the material volume fraction, and Comp* and C_{\mathrm{AM}}^{*} are the initially calculated structural compliance and AM cost for normalization purpose.

Note that Eq. (11) can be trivially re-configured into an AM cost constraint problem, which will also be addressed in this paper.

Solution of the optimization problem

Sensitivity result of material volume-constrained compliance problem is well known as shown in Eq. (12). The sensitivity analysis process will not be specified and interested readers can refer to Refs. [³²,³³] for details.

(12)\{\begin{array}{l}{L}^{\prime }={\displaystyle {\int }_{D}R\delta (\mathrm{\Phi }){\mathrm{\Phi }}^{\prime }\mathrm{d}\mathrm{\Omega }}\\ R=-\mathit{D}\mathit{e}(\mathit{u})\mathit{e}(\mathit{u})+\lambda \end{array},

where \lambda is the Lagrange multiplier, and L is the Lagrangian constructed to solve this optimization problem.

Then, the derivative information of the AM cost model will be presented here:

(13){\displaystyle \frac{\partial C_{\mathrm{AM}}}{\partial t}}=a_{\mathrm{1}}{\displaystyle \frac{\partial V_{\mathrm{support}}(\mathrm{\Phi })}{\partial t}}+a_{\mathrm{2}}{\displaystyle \frac{\partial h(\mathrm{\Phi })}{\partial t}}.

Regarding the first term in Eq. (13), volume of the support structure is calculated based on the current geometry; however, it is non-trivial to derive an analytical expression of V_{\mathrm{support}}(\mathrm{\Phi }). Then, we quantify this derivative by considering a worst case scenario [³⁴], where an approximated sensitivity result is derived which evaluates the effect of boundary evolution on the support volume change, as

(14){\displaystyle \frac{\partial V_{\mathrm{support}}(\mathrm{\Phi })}{\partial t}}=-{\displaystyle {\int }_{D^{+}}\delta (\mathrm{\Phi })}{\mathrm{\Phi }}^{\prime }\mathrm{d}\mathrm{\Omega },

where, D^{+} indicates the overhang areas and the boundary areas at bottom of the support. Different from Ref. [³⁴], the support volume is not developed as a function of the overhang angle since homogeneous porous support structure is often adopted by metal AM.

The second term in Eq. (13) is also not derivable. Here, we assume that the normal vector of the top surface coincide with the build direction, which is generally true. Therefore, the design update of the top surface with the velocity V_N would change the height of the part by V_N\mathrm{\Delta }t. Hence, we can approximate the sensitivity of the part height by

(15){\displaystyle \frac{\partial h(\mathrm{\Phi })}{\partial t}}={\displaystyle \frac{{\displaystyle {\int }_{D^{~}}\delta (\mathrm{\Phi }){\mathrm{\Phi }}^{\prime }\mathrm{d}\mathrm{\Omega }}}{{\displaystyle {\int }_{D^{~}}\delta (\mathrm{\Phi })\left|\nabla \mathrm{\Phi }\right|\mathrm{d}\mathrm{\Omega }}}}={\displaystyle \frac{{\displaystyle {\int }_{D^{~}}\delta (\mathrm{\Phi })V_N\left|\nabla \mathrm{\Phi }\right|\mathrm{d}\mathrm{\Omega }}}{{\displaystyle {\int }_{D^{~}}\delta (\mathrm{\Phi })\left|\nabla \mathrm{\Phi }\right|\mathrm{d}\mathrm{\Omega }}}},

where D^{~} means the top surface area of the part. In Eq. (15), V_N represents the change rate of the height of the top surface. However, the numerator performs the integration along the top surface which actually reflects the area that the top surface moves through. Therefore, the surface area is calculated and employed as the denominator to make the right side of Eq. (15) a proper calculation of the height change rate.

Therefore, the complete sensitivity result of Eq. (11) is summarized below:

(16)L^{\prime }={\displaystyle {\int }_D\left(-{\displaystyle \frac{w_{\mathrm{1}}}{Comp*}}\mathit{D}\mathit{e}(\mathit{u})\mathit{e}(\mathit{u})+\lambda \right)\delta (\mathrm{\Phi }){\mathrm{\Phi }}^{\prime }\mathrm{d}\mathrm{\Omega }}+{\displaystyle \frac{(\mathrm{1}-w_{\mathrm{1}})}{C_{AM}^{*}}}\left(-a_{\mathrm{1}}{\displaystyle {\int }_{D^{+}}\delta (\mathrm{\Phi })}{\mathrm{\Phi }}^{\prime }\mathrm{d}\mathrm{\Omega }+a_{\mathrm{2}}{\displaystyle \frac{{\displaystyle {\int }_{D^{~}}\delta (\mathrm{\Phi })}{\mathrm{\Phi }}^{\prime }\mathrm{d}\mathrm{\Omega }}{{\displaystyle {\int }_{D^{~}}\delta (\mathrm{\Phi })}|\nabla \mathrm{\Phi }|\mathrm{d}\mathrm{\Omega }}}\right).

Numerical examples

For all the numerical examples, the finite element analysis is performed based on fixed quadrilateral or hexahedral meshes and artificial weak material is employed for voids in order to avoid the singularity of the stiffness matrix, which is

(17){\mathit{D}}_{\mathrm{V}}={\mathrm{10}}^{-\mathrm{3}}\mathit{D},

in which {\mathit{D}}_{\mathrm{V}} is the elasticity tensor of the void.

The volume constraint is enforced by the Augmented Lagrange method which employs the Lagrange multiplier as

(18)

in which {\mu }_k is the penalization factor and \beta is its adjustment parameter.

The data used for cost modeling is listed in Table 1, and then Eq. (10) can be specified as shown in Eq. (19).

(19)C_{\mathrm{AM}}=\mathrm{65.36}+\mathrm{0.0288}{V}_{\mathrm{part}}(\mathrm{\Phi })+\mathrm{0.0116}{V}_{\mathrm{support}}(\mathrm{\Phi })+\mathrm{5.1283}h(\mathrm{\Phi }).

Tab.1 Data used for cost modeling

Parameter	Value
\rho	4.42 g/cm3
C_{\mathrm{material}}^{\mathrm{unit}}	4500 USD/kg
C_{\mathrm{argon}}^{\mathrm{unit}}	5.13 USD/m3
T_R	9 s
L_t	0.06mm
\bar{\rho }/\rho	0.4
C_{\mathrm{labor}}^{\mathrm{unit}}+C_{\mathrm{utility}}^{\mathrm{unit}}	120 USD/h
T_{\mathrm{setup}}	0.5 h
S_{\mathrm{rate}}	3.75 mm3/s

Cantilever problem

In this example, the cantilever structure is optimized. The design domain and boundary conditions are demonstrated in Fig. 1, where a unit force is loaded at the middle of the right edge, while the left edge is fully clamped. We assume the material has a Young’s modulus of 1.0 and a Poisson’s ratio of 0.3. The maximum material volume fraction of this optimization problem is 40%. Correspondingly, the optimization results with different weight factors are shown in Fig. 2. The related data is shown in Table 2.

Fig.1 Design domain and boundary condition of the cantilever problem (size: 100×50)

Full size|PPT slide

Fig.2 Optimization results with different weight factors. (a) w_1=1.0; (b) w_1=0.9; (c) w_1=0.7; (d) w_{\mathrm{1}}=\mathrm{0.5}; (e) w_{\mathrm{1}}=\mathrm{0.3}

Full size|PPT slide

Tab.2 Data of the cantilever optimization results (AM cost in USD if the dimension is measured in mm)

Weight factor w1	Strain energy	Support volume	Part height	AM cost
1.0	39.58	2126	50	404.03
0.9	41.48	1878	50	401.15
0.7	42.35	1695	50	399.03
0.5	51.12	1477	44	365.73
0.3	59.24	1329	40	343.51

We can see from Table 2 that, the weight factor determines the balance between the strain energy and the AM cost. A larger weight factor leads to a smaller strain energy value, but causes more AM cost. Then, in order to exactly constrain the manufacturing cost, another implementation is performed by constraining the AM cost with the limit of 380. The optimization result is shown in Fig. 3, and the convergence curve is plotted in Fig. 4. We can see that the AM cost has been successfully restricted at 380, and the finally derived strain energy and AM cost are 48.02 and 377.06, respectively.

Fig.3 Cantilever design with constrain AM cost. (a) Level set representation of the design; (b) voxel view of the design (the yellow color shows the support volume)

Full size|PPT slide

Fig.4 Convergence history

Full size|PPT slide

3D L-bracket problem

Next, the 3D L-bracket problem is explored. The problem setup is shown in Fig. 5, wherein the top surface is fixed and a line load of magnitude of unity per distance is imposed at the right-side top edge. We assume the material has a Young’s modulus of 1 and a Poisson’s ratio of 0.3. The maximum material volume fraction of the optimization problem is 40%. Correspondingly, the optimization results without and with cost constraint are demonstrated in Figs. 6 and 7, respectively. The AM cost limit is set to 815. The AM cost has the unit of USD if the dimension is measured in mm.

Fig.5 Design domain and boundary condition of the L-bracket problem (size: 60×60×15)

Full size|PPT slide

Fig.6 Optimization result without cost constraint

Full size|PPT slide

Fig.7 Optimization result with cost constraint

Full size|PPT slide

In this case, the part is built from bottom to top which is straightforward as demonstrated in Fig. 5. However, the top surface of the part is applied of the zero-displacement boundary condition, which cannot be evolved to reduce the part height. Therefore, we count the part height from top to bottom because raising the bottom surface is equivalent to reducing the part height. It is assumed that, the lowest point of the bottom surface contacts the build platen. Then, the part height is defined to be the projection of the distance between the highest point and the lowest point of the solid part onto the build direction. The part height control can be applied to the top surface (cantilever case), the bottom surface (L-bracket case), or simultaneously both.

We can see that very different results are produced after imposing the cost constraint. The constrained design (Fig. 7) has the AM cost of 812.76, which is cheaper than 852.12 of the unconstrained design in Fig. 6. Given the specifications, the part height of the constrained design is 54 and the support volume is 338, both of which are smaller than 60 and 1078 for the unconstrained design. Therefore, we can see that, the large overhangs in Fig. 6 does not show in Fig. 7, and in Fig. 7, there only exists a narrow internal void which is almost self-supporting by referring to the side cross-section view. Of course, the trade-off is that the strain energy of the constrained design is 828, which is higher than the 804 for the unconstrained design.

Another issue worth mentioning is that, both optimizations have generated internal voids, which is a difficult-to-manufacture feature for powder bed-based AM process. For instance, the internal support structure can hardly be removed, and even if they are self-supporting, powders would be trapped inside and are tough to remove. Therefore, casting level set method [³⁵–³⁷] may be a better option for our optimization problem, since it produces undercut- and internal void-free topological design, which makes it trivial to post-processing the AM part.

Conclusions

AM cost-constrained level set topology optimization was performed in this work. With the AM cost for the EOS M290 DMLS system, we can see that the support volume and part height play the key role in determining the manufacturing cost, especially for the latter which significantly increases the layer number and thus the powder recoating time. Therefore, reducing the support volume and part height would be effective in decreasing the manufacturing cost. In this paper, the proposed algorithm has demonstrated the effectiveness of minimizing the AM cost by reducing the part height and support volume simultaneously. Very importantly, the sensitivities on both the support volume and part height cannot be analytically calculated, and to address this issue, approximate sensitivities have been developed. Regarding this approximation, the related technical foundation has been discussed, and the effectiveness has been proved by the numerical examples.

The current algorithm is developed under the level set framework. However, with the density-based topology optimization method, many researchers have used the density gradient information to approximate the boundary profile and related normal vectors [³⁸]. Therefore, it is trivial to numerically calculate the part height and the overhang boundaries based on the density field information. Also, performing shape optimization under the density-based framework is also possible [³⁹], which can achieve a similar design update effect as compared with level set method. Therefore, we can conclude that a similar algorithm can be developed using density-based topology optimization method, even though extensive calibrations are required.

For future work, topology optimization method without generating internal voids would need to be addressed. Unlike the current casting level set method, undercut features may be allowed but those areas should be accessible by 5-axis CNC machining, so that support structures can be removed and those surface areas can be polished for good surface quality. Another issue to address is to include post-machining cost in the AM cost model, so that we could concurrently constrain the post-machining cost, which currently can be as high as the AM cost.

Acknowledgement

The authors would like to acknowledge the support from the National Science Foundation (CMMI-1634261).

References

1 Bendsøe M P, Sigmund O. Topology Optimization. Berlin: Springer, 2004

2 Sigmund O, Maute K. Topology optimization approaches. Structural and Multidisciplinary Optimization, 2013, 48(6): 1031–1055

3 van Dijk N P, Maute K, Langelaar M, Level-set methods for structural topology optimization: A review. Structural and Multidisciplinary Optimization, 2013, 48(3): 437–472

4 Liu J, Ma Y. A survey of manufacturing oriented topology optimization methods. Advances in Engineering Software, 2016, 100: 161–175

5 Liu J, Gaynor A T, Chen S, Current and future trends in topology optimization for additive manufacturing. Structural and Multidisciplinary Optimization, 2018, 57(6): 2457–2483

6 Gaynor A T, Guest J K. Topology optimization considering overhang constraints: Eliminating sacrificial support material in additive manufacturing through design. Structural and Multidisciplinary Optimization, 2016, 54(5): 1157–1172

7 Langelaar M. An additive manufacturing filter for topology optimization of print-ready designs. Structural and Multidisciplinary Optimization, 2017, 55(3): 871–883

8 Allaire G, Dapogny C, Estevez R, Structural optimization under overhang constraints imposed by additive manufacturing technologies. Journal of Computational Physics, 2017, 351: 295–328

9 Guo X, Zhou J, Zhang W, Self-supporting structure design in additive manufacturing through explicit topology optimization. Computer Methods in Applied Mechanics and Engineering, 2017, 323: 27–63

10 Liu J, To A C. Deposition path planning-integrated structural topology optimization for 3D additive manufacturing subject to self-support constraint. Computer Aided Design, 2017, 91: 27–45

11 Zhang W, Zhou L. Topology optimization of self-supporting structures with polygon features for additive manufacturing. Computer Methods in Applied Mechanics and Engineering, 2018, 334: 56–78

12 Zhang P, Liu J, To A C. Role of anisotropic properties on topology optimization of additive manufactured load bearing structures. Scripta Materialia, 2017, 135: 148–152

13 Liu J, Yu H. Concurrent deposition path planning and structural topology optimization for additive manufacturing. Rapid Prototyping Journal, 2017, 23(5): 930–942

14 Dapogny C, Estevez R, Faure A, Shape and topology optimization considering anisotropic features induced by additive manufacturing processes. Computer Methods in Applied Mechanics and Engineering, 2018, 344: 626–665

15 Liu J, Ma Y, Qureshi A J, Light-weight shape and topology optimization with hybrid deposition path planning for FDM parts. International Journal of Advanced Manufacturing Technology, 2018, 97(1‒4): 1123–1135 doi:10.1007/s00170-018-1955-4

16 Mirzendehdel A M, Rankouhi B, Suresh K. Strength-based topology optimization for anisotropic parts. Additive Manufacturing, 2018, 19: 104–113

17 Wang Y, Chen F, Wang M Y. Concurrent design with connectable graded microstructures. Computer Methods in Applied Mechanics and Engineering, 2017, 317: 84–101

18 Wang Y, Zhang L, Daynes S, Design of graded lattice structure with optimized mesostructures for additive manufacturing. Materials & Design, 2018, 142: 114–123

19 Cheng L, Liu J, Liang X, Coupling lattice structure topology optimization with design-dependent feature evolution for additive manufactured heat conduction design. Computer Methods in Applied Mechanics and Engineering, 2018, 332: 408–439

20 Cheng L, Liu J, To A C. Concurrent lattice infill with feature evolution optimization for additive manufactured heat conduction design. Structural and Multidisciplinary Optimization, 2018, 58(2): 511–535

21 Vogiatzis P, Ma M, Chen S, Computational design and additive manufacturing of periodic conformal metasurfaces by synthesizing topology optimization with conformal mapping. Computer Methods in Applied Mechanics and Engineering, 2018, 328: 477–497

22 Liu J, To A C. Topology optimization for hybrid additive-subtractive manufacturing. Structural and Multidisciplinary Optimization, 2017, 55(4): 1281–1299

23 Ruffo M, Tuck C, Hague R. Cost estimation for rapid manufacturing—Laser sintering production for low to medium volumes. Proceedings of the Institution of Mechanical Engineers. Part B, Journal of Engineering Manufacture, 2006, 220(9): 1417–1427

24 Ruffo M, Hague R. Cost estimation for rapid manufacturing: Simultaneous production of mixed components using laser sintering. Proceedings of the Institution of Mechanical Engineers. Part B, Journal of Engineering Manufacture, 2007, 221(11): 1585–1591

25 Baumers M, Tuck C, Wildman R, Combined build-time, energy consumption and cost estimation for direct metal laser sintering. 2012. Retrieved from http://sffsymposium.engr.utexas.edu/Manuscripts/2012/2012-71-Baumers.pdf

26 Huang R, Ulu E, Kara L B, Cost minimization in metal additive manufacturing using concurrent structure and process optimization. In: Proceedings of ASME 2017 International Design Engineering Technical Conferences and Computers and Information in Engineering Conference. Cleveland, 2017, V02AT03A030

27 Barclift M, Armstrong A, Simpson T W, CAD-integrated cost estimation and build orientation optimization to support design for metal additive manufacturing. In: Proceedings of ASME 2017 International Design Engineering Technical Conferences and Computers and Information in Engineering Conference. Cleveland, 2017, V02AT03A035

28 Wang S, Wang M Y. Radial basis functions and level set method for structural topology optimization. International Journal for Numerical Methods in Engineering, 2006, 65(12): 2060–2090

29 Jiang L, Chen S. Parametric structural shape & topology optimization with a variational distance-regularized level set method. Computer Methods in Applied Mechanics and Engineering, 2017, 321: 316–336

30 Jiang L, Chen S, Jiao X. Parametric shape and topology optimization: A new level set approach based on cardinal basis functions. International Journal for Numerical Methods in Engineering, 2018, 114(1): 66–87

31 Liu J, Yu H, To A C. Porous structure design through Blinn transformation-based level set method. Structural and Multidisciplinary Optimization, 2018, 57(2): 849–864

32 Wang M Y, Wang X, Guo D. A level set method for structural topology optimization. Computer Methods in Applied Mechanics and Engineering, 2003, 192(1–2): 227–246

33 Allaire G, Jouve F, Toader A M. Structural optimization using sensitivity analysis and a level-set method. Journal of Computational Physics, 2004, 194(1): 363–393

34 Mirzendehdel A M, Suresh K. Support structure constrained topology optimization for additive manufacturing. Computer Aided Design, 2016, 81: 1–13

35 Xia Q, Shi T, Wang M Y, A level set based method for the optimization of cast part. Structural and Multidisciplinary Optimization, 2010, 41(5): 735–747

36 Allaire G, Jouve F, Michailidis G. Casting constraints in structural optimization via a level-set method. In: Proceedings of the 10th World Congress on Structural and Multidisciplinary Optimization. Orlando: HAL, 2013, hal-01088775

37 Wang Y, Kang Z. Structural shape and topology optimization of cast parts using level set method. International Journal for Numerical Methods in Engineering, 2017, 111(13): 1252–1273

38 Qian X P. Undercut and overhang angle control in topology optimization: A density gradient based integral approach. International Journal for Numerical Methods in Engineering, 2017, 111(3): 247–272

39 Gersborg A R, Andreasen C S. An explicit parameterization for casting constraints in gradient driven topology optimization. Structural and Multidisciplinary Optimization, 2011, 44(6): 875–881