A Novel Mathematical Formulation for Density-Based Topology Optimization Method Considering Multi-Axis Machining Constraint

Topology optimization has become an important tool to determine the optimal shape for maximum performance subject to given design constraints. These performance objective functions and design constraints affect the applicability of the design and typically include the mass, compliance, stress, or displacement due to loading conditions. Research for modern topology optimization first started with the work of Bendsøe and Kikuchi [1], and new methods have since been implemented such as the solid isotropic material with penalization (SIMP) method [2] and the level-set method [3], both of which can be found in commercial software. For the SIMP method, each element is assigned a density variable to control the material distribution and remove material from gradient-based optimization. For the level-set method, the boundary of the shape geometry is the zero level set of an implicit function. The parameters defining the function are optimized to produce the optimal shape.

These methods can produce complex geometries not typically considered by a human engineer or designer. Consequently, it can be difficult to manufacture these designs and may require extensive post-processing [4]. To control the geometric complexity, other new methods have been developed such as moving morphable components (MMC) [5,6], moving morphable voids (MMV) [7], and geometry projection [8–15]. MMC and MMV use controlled geometries to add or remove material from the design domain and determine the optimal shape. Geometry projection methods [8,10] use geometries to describe the density distribution on a background mesh. Additionally, integrating the topology optimization process with computer-aided design (CAD) has been addressed as well to further reduce post-processing [16,17].

This paper mainly focuses on the multi-axis machining constraints, which is a widely used technique in subtractive manufacturing for metal component production. For multi-axis machining, the relative position and orientation of cutting tools and workpiece can be manipulated in 4 or 5 degrees of freedom. The unnecessary material is removed by machining tool until the desired shape achieves. Compared with 2.5D milling, the multi-axis machining allows more design freedom. Besides machining, metal additive manufacturing is also a well-known technique to print freeform design with high design freedom, while the material strength and fatigue properties of parts made of AM still have large gap compared with metal machining. For high-strength aerospace or naval structures, multi-axis machining is an ideal option for manufacturing parts. Thus, topology optimization approaches considering machining constraint is a necessary and valuable research to achieve a trade-off between the available manufacturing technologies. In recent years, a few effective approaches for multi-axis machining-based topology optimization are proposed. For 2.5D milling, Liu and Ma [18] propose an explicit feature-based level-set method, where the feature fitting algorithm is incorporated into the boundary evolvement process. Furthermore, Liu and To [19] proposed a novel CAD-based topology optimization system for milling constraint, where feature and dynamic modeling history is incorporated into the optimization process. For multi-axis machining, Mirzendehdel et al. [20] presented a topology optimization framework using convolutions in configuration space to enable manufactured design using multi-axis machining, where an inaccessibility measure field (IMF) of design domain is introduced to identify nonmanufacturable features. Morris et al. [21] proposed a level-set-based topology optimization method for multi-axis machining, where the advection velocity at every iteration is modified to ensure the manufacturability conditions. This level-set-based method simulates the subtractive process by cutting materials accessible to the machine in every iteration until optimization converges. Langelaar [22] proposed a density-based approach incorporating multi-axis restrictions, where a filter technique based on KS aggregation function [23] is introduced to transform an input design field into a manufacturable geometry. Mirzendehdel et al. [24] proposed framework to optimize the build orientation with respect to removability of support structures, where a Pareto-optimality criterion is implemented to achieve a trade-off between build orientation and machining setup. Meanwhile, Mirzendehdel et al. [25] further extended the notion of IMF to support structures of additive manufacturing to identify the inaccessible points, and the IMF is integrated into the sensitivity field to guide the TO results. This approach enables optimal design can be additive manufactured and post-processed for support removal using machining techniques. Besides the multi-axis machining, several other related topology optimization studies regarding casting constraints can also be found in Refs. [26–29].

The aim of this paper is to propose a novel density-based approach to optimize parts for multi-axis machining restrictions, where a concise aggregation-free density projection formulation is demonstrated. The proposed method has concise mathematical expression without aggregation function, which is easy for sensitivity analysis and numerical implementation. The length to diameter (L:D) ratio constraint is taken into consideration in the proposed scheme. The current method is in the framework of density-based method, which is able to inherit the existing sensitivity from standard density method with chain rule for different physical problems. The subtractive manufacturing constraints can be satisfied after density projection. The remainder of this paper is organized as follows. Section 2 presents the mathematical formulation of density projection algorithm. Section 3 demonstrates three typical topology optimization problems with machining constraint from 2D to 3D to verify the effectiveness of the proposed method, followed by conclusions in Sec. 4.

The milling head is generally described by a cylinder capped with a hemisphere oriented in the milling direction. The end of the head cylinder is assumed to extend infinitely far away from the cutting surface. For 3-axis machining, the workpiece is still while the cutting tool moves along the three mutually perpendicular axes to mill the part, which is the most widely used technique to create a metal part. While for 5-axis machining process, the milling tool can be oriented to approach the surface from an arbitrary direction. The machinable part is determined by whether all its surface points are accessible by a tool bit without any intersection with the interior of the part. As described by Ref. [22], machinability means that the density field must be monotonically increasing in the insertion direction without any internal holes. To enforce the milling constraints, a Heaviside function-based projection method is proposed here to construct the monotonically increasing density field in the milling direction.

As shown in Fig. 2, d₀ is a threshold, and the operation $∥⋅∥$ denotes the Euclidean norm. Note that the computational cost will increase with the increase of the element number. Table 1 demonstrates the computational costs for rectangle domain based on matlab code [20] for one milling direction. The element set $M$ is just computed once and store the information in a matrix, which will be reused during optimization.

As demonstrated in Fig. 1, the gray element in physical field ρ_p is computed through summation of blue element $(M)$ in fictitious density field ρ_f according to Eq. (1). $H1(⋅)$ and $H2(⋅)$ denote the approximate Heaviside functions defined as

The shape of Heaviside functions is plotted in Fig. 3.

The reason for introducing two Heaviside functions with different parameters is due to numerical stability during optimization process after numerical experiments. In general, low coefficient values are preferred, which will make the objective function smooth enough. However, a low coefficient value will weaken the filter effect, which is critical to enable the machinable design. Meanwhile, a higher coefficient value will generate less intermediate density in final optimal designs. Thus, a reasonable coefficient for Heaviside function needs to be chosen to achieve a trade-off between smoothness and filter effects.

Besides the density field must be monotonically increasing in the insertion direction. The other crucial machinable constraint is the length to diameter (L:D) ratio [31]. The (L:D) ratio is an important concept for machinable design, which has not been considered in the previous topology optimization research. As shown in Fig. 5, the hole is too deep for the machine tool to reach. Although this case satisfies the density nondecreasing requirement along the insertion direction, the part still cannot be machinable due to the inaccessible deep hole. This situation can be resolved if the cone angle geometrical constraints are introduced. As shown in Fig. 6, a cylinder cone angle θ is defined to further constrain design as follows:

A simple schematic is shown here to demonstrate the difference between two projection schemes. As shown in Fig. 8, the projection design mapped from fictitious field generates small deep hole if no CES constraint is applied, while no deep hole appears once the CES constraint (θ = 30 deg) is implemented.

The first 2D numerical example for machining-constrained topology optimization is demonstrated in Fig. 9. The design domain is uniformly meshed by 100 × 100 quad elements with unit length. The loading F = 1 is applied on the right bottom corner, and left side is fully fixed. The volume fraction constraint $V¯$ is chosen as 0.2. The elastic constants are chosen as follows: Elastic modulus E = 1 and Poisson’s ratio μ = 0.3. The filter radius for design variable is chosen as r_min = 4. The initial design is plotted in Fig. 11. To make a comparison, different milling direction constraints are applied to produce an optimized design. The reference solution without any manufacturing constraints is demonstrated in Fig. 10. Obviously, the reference solution cannot be manufactured if machining operations are limited in the x, y-plane. The initial and optimized designs for single milling orientation are demonstrated in Fig. 11, and convergence history is plotted in Fig. 12. Compared with the reference solution, the optimized part for unidirectional milling restriction shows strong limitation of design freedom, which also has great impact on structural compliance. Note that the designs with single milling orientation constraint are considerably more compliant with respect to reference design. For multiple tool orientations, the optimized designs are shown in Fig. 13, and convergence history is plotted in Fig. 14. The structural performance with multiple directions is better than the single orientation, while the compliance value is still higher than the reference design. This simple 2D numerical example proves the effectiveness of proposed method to force the optimization process toward a different solution considering multi-axis machining constraints.

In this subsection, we focus on verifying the effectiveness of proposed method for 3D designs. In the first test example, a three-dimensional cantilever beam example is presented for compliance minimization (Fig. 15). The unit force is applied at the midpoint of right bottom edge, and the left end is fully fixed. The design domain is discretized by a 144 × 48 × 48 hexahedral mesh with unit element length. The elastic constants are chosen as follows: Elastic modulus E = 1 and Poisson’s ratio μ = 0.3. The volume fraction constraint is set to be 0.3. The filter radius for multi-axis machining optimization is chosen as r_min = 4. A reference solution is demonstrated in Fig. 17, where no machining constraints are applied. As shown in Fig. 17, the reference optimized design consists of hollow chamber, which is not accessible by milling tools. It is worth to mention that the initialization of fictitious field for optimization is shown in Fig. 16. To ensure manufacturability of the optimized design, the results of two orientation machining optimizations for the cantilever beam are plotted in Fig. 18. The orientation of machine tool is described by the vector (a, b, c) in the local coordinate as shown in Fig. 15. The compliance value of designs obtained from different orientations (Fig. 18) is close and slightly higher than the reference design. For multi-axis machining constraints, the optimized result is demonstrated in Fig. 19 with 6 different milling orientations. However, compliances differ only slightly among all designs, which denotes the number of different directions has little effect on the current issue.

In this section, the results with CES constraint are presented. The CES angle limitation is chosen as $θ*=30deg$⁠. The filter radius is r_min = 1.5. In this example, we use the optimized design from standard density-based result to initialize the fictitous field. The initialization of fictitous domain and corresponding void field are plotted in Fig. 21. Based on numerical experiments, level-set initialization method or initialization based on standard density result are better than uniform field. The problem for the uniform field initialization is that the convergence rate is very slow. Thus, we do not recommend applying uniform density distribution to initialize the design domain. Different milling directions are applied to generate multiple machinable solutions. The initial designs and optimized designs are demonstrated in Fig. 22. Note that the initial design is generated through density projection based on the void field (Fig. 21(a)). Obviously, compared with results in Fig. 20, no deep hole or narrow chamber is found in the final design. The optimal compliance values with CES constraint are generally higher than optimal designs without CES constraint.

In this section, an MBB design example is presented for compliance minimization. The unit force is applied at the center of top face. Left and right bottom edges are fully fixed as shown in Fig. 23. Due to the symmetry, only half of the MBB beam is chosen for optimization and is discretized by a 144 × 48 × 48 hexahedral mesh with unit element length. The material properties are the same as the previous example. The filter size for design variables is selected as r_min = 4. The volume constraint is chosen as $V¯=0.2$⁠. The reference result without machining constraints obtained from standard density-based method is shown in Fig. 25. This reference design is not manufacturable through machining as it contains several inaccessible internal surfaces. Similar to the previous example, we start by considering milling operation in two opposite directions. The initialization of fictitious field for optimization is demonstrated in Fig. 24. Several machinable designs are generated through the proposed method as demonstrated in Fig. 26. The machinable MBB designs have very different optimized configurations compared with the reference case, while the compliance values of machinable designs are remarkably competitive with the reference one.

For a multi-axis milling scenario, the set of milling orientations is extended as shown in Fig. 27. As plotted in Fig. 27, the machinable design is obtained by six mutually orthogonal milling directions. The compliance value (C = 10.05) of optimized design is slightly outperformed than the reference design (C = 10.72). Clearly, the reference design is only a local optimum as gradient-based method is only able to converge to a local minimum.

In recent years, developing advanced topology optimization methods for conventional subtractive manufacturing becomes a new trend in this field. In this paper, a novel mathematical formulation to impose multi-axis machining restrictions in the framework of density-based method is proposed. A simple density mapping method based on Heaviside function for multi-axis restrictions is demonstrated in detail, where no aggregation functions (e.g., KS-function [23]) are involved. Several 2D and 3D numerical examples are demonstrated to validate the effectiveness of proposed method. In general, some small features generated by topology optimization are difficult to reach like a thin and deep pocket. For multi-axis machining, deep and narrow areas require specialized tooling and are time-consuming to produce. These narrow areas may require additional equipment setup and increase the cost of a component. Increasing the filter radius as described in the proposed method can effectively reduce these narrow regions. Besides, the CES constraint we introduced in this paper can effectively remove these narrow or deep void features. We need to mention that the current method is assuming that the multi-axis tool’s motion is only to be translational from outside the part, which is a limitation of current design method. Compared with Amir’s method [20], the proposed method does not consider the fixturing devices for the machining setup, which will be further investigated and extended to realistic geometry in the future. Besides, optimizing the machining direction is another research direction which will be explored in the future.

The financial support for this work from National Science Foundation (CMMI-1634261) is gratefully acknowledged.

There are no conflicts of interest.

The datasets generated and supporting the findings of this article are obtainable from the corresponding author upon reasonable request. The authors attest that all data for this study are included in the paper. No data, models, or code were generated or used for this paper.