## Abstract

This paper proposes a novel density-based method for structural design considering restrictions of multi-axis machining processes. A new mathematical formulation based on Heaviside function is presented to transform the design field into a geometry which can be manufactured by multi-axis machining process. The formulation is developed for 5-axis machining, which can be also applied to 2.5D milling restriction. The filter techniques are incorporated to effectively control the minimum size of void region. The proposed method is demonstrated by solving the compliance minimization problem for different machinable freeform designs. The length to diameter (L:D) ratio geometric constraint is introduced to ensure the machinable design, where deep hole or narrow chamber features are avoided using proposed method. Several two- and three-dimensional numerical examples are presented and discussed in detail.

## 1 Introduction

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.

## 2 Mathematical Formulation for Multi-Axis Machining

### 2.1 Machining Restriction.

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*_{0} is a threshold, and the operation $\u2225\u22c5\u2225$ 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.

Element number | 10 × 10 × 10 | 25 × 25 × 25 | 50 × 50 × 50 | 75 × 75 × 75 | 100 × 100 × 100 |
---|---|---|---|---|---|

Time cost | 0.013 s | 0.21 s | 8.53 s | 18.78 s | 89.52 s |

Element number | 10 × 10 × 10 | 25 × 25 × 25 | 50 × 50 × 50 | 75 × 75 × 75 | 100 × 100 × 100 |
---|---|---|---|---|---|

Time cost | 0.013 s | 0.21 s | 8.53 s | 18.78 s | 89.52 s |

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(\u22c5)$ and $H2(\u22c5)$ 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.

*ρ*

_{v}and fictitious domain

*ρ*

_{f}can be defined as follows:

*ρ*

_{v}= 1 represents the void element, while

*ρ*

_{v}= 0 means solid element. (·)

^{(i)}means the

*i*th element.

*d*

_{fi}is weight factor for filter defined as

*r*

_{min}is the filter radius. Δ(

*f*,

*i*) denotes the distance between the center of element

*i*and

*f*.

*N*

_{f}is the set of elements

*i*for which the distance Δ(

*f*,

*i*) is less than

*r*

_{min}. More details regarding filtering techniques can be found in Ref. [30]. Therefore, the physical density

*ρ*

_{p}can be explicitly expressed as

*k*= 1, ···,

*n*), the composite physical density field can be simply expressed as

*k*milling direction. A typical projection from void density field to physical density field is demonstrated in Fig. 4.

### 2.2 Sensitivity Analysis.

*f*(

*ρ*

_{v}), the sensitivity with respect to the design variables

*ρ*

_{v}can be obtained through the chain rule as follows:

*n*is the number of machining direction. The explicit form of the term $\u2202\rho p(i)/\u2202\rho v(l)$ can be expressed as

### 2.3 Machining Restriction Considering Length to Diameter (L:D) Ratio.

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:

*D*is the diameter and

*L*is the length as shown in Fig. 6. This geometrical constraint means that there should be no material inside the circular cone envelop surface (CES) measured from arbitrary point on the design external surface. This CES geometry constraint can effectively avoid feature with small (L:D) ratio, e.g., deep holes.

*θ*is CES angle, which will be determined by CNC machine and tool shape.

**is the minimum distance from centroid of the element to red line as shown in Fig. 7, and**

*d**h*is the projection length. The sensitivity analysis process is similar to Sec. 2.2, which is omitted here.

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.

### 2.4 Objective and Constraints.

*ρ*

_{v}is the design variables, and

*C*is the objective defined by the structural compliance.

*ρ*

_{p}is the physical density field, Ω the design domain, and

*V*

_{prescribe}the prescribed volume fraction.

**is the unknown displacement field and**

*u***is the strain resolved by finite element method (FEM), and**

*ɛ***is the elastic tensor matrix. The sensitivity of objective function and constraint with respect to design variable**

*D**ρ*

_{v}can be obtained through the chain rule.

## 3 Numerical Examples

### 3.1 Machining-Based Optimization for a 2D Cantilever Beam.

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\xaf$ 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.

### 3.2 Machining-Based Optimization for 3D Cantilever Beam

#### 3.2.1 Optimized Design Without Cone Envelop Surface Constraint.

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.

*r*

_{min}are selected to produce diverse designs as shown in Fig. 20. The machining tool orientations (six directions) are shown in Fig. 19. As mentioned by Refs. [32–34], the filter size is a straightforward and effective way to control the member size of optimal design. As plotted in Fig. 20, for small radius size (

*r*

_{min}= 1), small and narrow hollow chambers are found, which may not be accessible by milling tool. However, the hollow chamber size increases after increasing the filter radius

*r*

_{min}, while the compliance values (

*C*) of multiple designs are close in this case

#### 3.2.2 Optimized Design With Cone Envelop Surface Constraint.

In this section, the results with CES constraint are presented. The CES angle limitation is chosen as $\theta *=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.

### 3.3 Machining-Based Optimization for 3D Messerschmitt–Bölkow–Blohm (MBB) Beam.

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\xaf=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.

## 4 Conclusion

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.

## Acknowledgment

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

## Conflict of Interest

There are no conflicts of interest.

## Data Availability Statement

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.

## References

*CNC Machining Handbook: Building, Programming, and*

*Implementation*