High-Accuracy and Efficient Simulation of Numerical Control Machining Using Tri-Level Grid and Envelope Theory

Abstract

Virtual simulation of high-resolution multi-axis machining processes nowadays plays an important role in the production of complex parts in various industries. In order to improve the surface quality and productivity, process parameters, such as spindle speed, feedrate, and depth of cut, need to be optimized by using an accurate process model of milling, which requires both the fast virtual prototyping of machined part geometry for tool path verification and accurate determination of cutter–workpiece engagement for cutting force predictions. Under these circumstances, this paper presents an effective volumetric method that can accurately provide the required geometric information with high and stable computational efficiency under the condition of high grid resolution. The proposed method is built on a tri-level grid, which applies two levels of adaptive refinement in space decomposition to abolish the adverse effect of a large fine-level branching factor on its efficiency. Since hierarchical space decomposition is used, this multi-level representation enables the batch processing of affected voxels and minimal intersection calculations, achieving fast and accurate modeling results. To calculate the instantaneous engagement region, the immersion angles are obtained by fusing the intersection points between the bottom-level voxel edges and the cutter surface, which are then trimmed by feasible contact arcs determined using envelope theory. In a series of test cases, the proposed method shows higher efficiency than the tri-dexel model and stronger applicability in high-precision machining than the two-level grid.

1. Introduction

Multi-axis milling is a versatile machining technology used for producing delicate parts with freeform surfaces [1,2]. To verify and optimize the milling process before the parts are manufactured on real machine tools, virtual machining (VM) technology is proposed to verify NC codes, reduce cycle times and surface errors, and predict instantaneous cutting states. However, these goals cannot be achieved without accurate geometrical and physical modeling of the milling process. The former aims to compute the cutter-swept volume (CSV), machined part geometry, and cutter–workpiece engagement (CWE), while the latter conducts dynamics analysis using the obtained information.

Cutter–workpiece engagement refers to the instantaneous contact region between the cutting tool and the in-process workpiece (IPW). Based on this geometry, the in-cut edge segments can be determined for cutting force predictions, which subsequently initiate the physical modeling process. In freeform surface machining, CWE varies along a tool path, and in general, unless very simple workpiece geometry is machined, it is challenging to find an exact solution, which has attracted extensive investigation [3,4,5,6,7,8].

During the milling process, whether the cutting forces are computed or not, it is necessary to geometrically represent the machined part according to the cutting tools and their trajectories, that is, according to tool paths. This calculation process is referred to as a workpiece update in this research, and the obtained geometric information is very useful in general NC milling simulations. First, this information indicates whether the desired part can be manufactured by designed tool paths within the specified tolerance. Second, unexpected collisions between the milling tool and in-process workpieces, such as local gouging, rear gouging, and global interference, can also be detected [9,10,11]. Such collisions will lead to severely compromised surface qualities or even safety accidents.

The basic idea for NC simulations is to remove every cutter-swept volume from the blank until the final machined part is obtained [3,4,12,13]. Mathematically, this volume is defined as the totality of all 3D points that belong to the trace of a cutter moving along an arbitrary path. It can be modeled as the volume enclosed by an ingress surface at the initial location, an egress surface at the final location, and an envelope surface in between, which is the collection of the grazing points on each tool instance. According to the envelope theory, the product of the velocity and the surface normal vector on a grazing point is always equal to zero. Another way is to remove the workpiece material inside the tool instance at each cutter location (CL) point sampled along a tool path, which has also been widely used. Meanwhile, the instantaneous engagement region between the cutter and in-process workpiece is extracted and represented as a CWE map, where immersion angles are computed by a further discretization of the contact area.

In modern manufacturing, an efficient geometric modeling method will accelerate the optimization process of various machining parameters significantly. Meanwhile, with the rapid development of automotive, aerospace, and biomedical industries, the demand for high-quality machine parts is increasing sharply. However, few of the existing virtual machining technologies can trade off high efficiency against high quality satisfactorily, leading to the necessity of the present study.

The rest of this paper is structured as follows: Section 2 summarizes the geometric modeling methods utilized in virtual machining and points out their benefits and drawbacks to help identify the preferred modeling frame. Section 3 illustrates the tri-level representation of the workpiece, and Section 4 expounds on the implicit function of a generalized cutter. The simulation procedures of the proposed method are then elaborated in Section 5. Section 6 provides a quantitative comparison to demonstrate the computational improvement the proposed method brings over the tri-dexel model and two-level grid, followed by conclusions and future work in Section 7.

2. Existing Methods

A variety of geometric modeling methods have been proposed and applied in NC milling simulations, aiming to achieve both computational accuracy and efficiency in modeling the machining geometry while being flexible enough to handle different types of cutters and tool movements. These existing methods can be roughly classified into three groups: solid modeling, vector modeling, and voxel modeling. Solid modeling is able to provide the highest level of precision by utilizing the exact representation of the cutter and the workpiece geometry, while the other two groups are discrete modeling approaches, which strive to achieve higher computational speed at the cost of accuracy.

In solid modeling [14,15,16], the workpiece material removal process is simulated by making use of the geometric and topological algorithms built into the solid modeler kernel. Specifically, constructive solid geometry (CSG) and boundary representation (B-rep) are the most widely used data structures in the field of solid modeling. CSG builds a solid object by combining primitive objects with Boolean operators, while B-rep describes an object with algebraic and parametric surfaces in the forms of Bezier, B-spline, and NURBS (Non-Uniform Rational B-Splines) patches. Compared to CSG, B-rep offers a richer set of operations in terms of modifying workpiece surfaces, hence gaining more attention in modern NC milling simulations. Because the geometry of both the in-process workpiece and cutter-swept volume keeps changing, large amounts of intersections between surfaces need to be solved numerically in order to determine the intersection curves for building new patches. A typical cutter location data file usually contains a huge number of cutter movements, so solving these complex intersections has become one of the most time-consuming processes in virtual machining. Meanwhile, the computation load needed to maintain and process the tree of numerous tool swept volumes in order to generate the machined part geometry grows significantly as the program runs. In addition, the high complexity of the algorithms required to manipulate intersections in multi-axis NC simulation is more likely to result in computational breakdowns as a result of round-off errors.

The second modeling method utilizes a collection of spatial geometric elements to define the cutter and workpiece for simplified Boolean operations. The typical example of such a discrete modeling method is called Z-map representation [17,18,19,20], which represents the workpiece by a 2D array of vectors with one end on a fixed coordinate plane and the other end on the top surface of the workpiece. To simulate the cutting operations, these vectors are updated as the cutter sweeps through different regions of the grid to capture their new heights, i.e., the in-process workpiece boundary. During the milling process, the vectors can be updated easily if the cutter geometry is fairly simple. For a more complicated cutter like a bull-end mill, these vectors are usually updated by solving high-order nonlinear equations. In order to clip Z-map vectors with less computation time, Chung et al. [17] represented the cutter-swept surface generated by translational motion by an implicit function of $z=f\left(x,y\right)$ in three-axis milling and determined quadruple roots using an analytic formula; however, this approach may introduce intolerable numerical errors. Tunç et al. [18] conducted cutting force and stability simulations through the extended Z-mapping approach, where multiple dexels are connected at one vector location. Li et al. [19] developed an improved Z-map algorithm that combines servo rectangular encirclement and the angle summation method to simulate the part surface topography, while Xiao et al. [20] performed the same task based on the sequential quadratic programming algorithm. Eventually, to resolve the issue of poor sensitivity along the vector direction of any single set of parallel vectors, the most advanced tri-dexel representation, also called triple-ray representation, was proposed, which is essentially three orthogonal dexel models combined together, and it has been used by various computer-aided manufacturing (CAM) systems [6,15,21].

Voxel modeling or space partitioning is the third geometric modeling method employed by virtual machining [22,23,24,25,26,27]. This volumetric modeling method has an evident advantage over vector modeling since the fundamental operation involved in workpiece updates is the simple deactivation of any cubic element falling into the cutter-swept volume. Such an operation has a very low computational cost compared to the operation involved in updating other classes of models; hence, it is more efficient. A basic uniform grid, which means the entire volume is uniformly subdivided, is the most fundamental data structure for workpiece representation [22,23,24]. However, both memory consumption and computation time will increase dramatically when a high grid resolution is required to provide accurate enough modeling results. In order to mitigate such overhead, more advanced voxel grids are proposed for a milling simulation, which can localize and reduce the finest update calculations effectively. Yau et al. [25] utilized the octree to simulate the cutting process with a generalized automatically programmed tool (APT), while Joy et al. [26,27] computed the machined part geometry using the two-level grid. In general, voxel modeling, due to its low computational complexity and hierarchical space decomposition, is regarded as the most promising approach, as it can provide the best combination of accuracy, efficiency, and flexibility among all geometric modeling methods.

Even though the two-level grid has been widely used by many researchers, the limitation of such a fixed hierarchical structure is still obvious. Since the available grid resolution is measured by the product of the branching factors at both levels, a higher resolution suggests that either of the two branching factors should be increased. A large top-level branching factor increases the number of tiles sharply and hence the memory usage, while a large bottom-level branching factor indicates a large ratio of inactive-to-active voxels and redundant storage. More importantly, a large bottom-level subdivision factor will introduce a large quantity of repetitive fine-level traversal operations since numerous coarse voxels need to be accessed for an update multiple times due to the small sampling interval, leading to a longer computation time. Such a limitation reduces its adaptivity in geometry for extreme-scale volume modeling. Moreover, since a voxel model is a fully discrete representation of a solid object, the precision of the voxel-based CWE geometry is not satisfactory [22,23,24]. Thus, a voxel-based CWE determination method that can achieve the subvoxel resolution and precision also needs to be developed.

3. Representation of Workpiece

3.1. Basic Concepts

Voxel modeling, as a typical spatial occupancy enumeration type of a solid representation scheme, is mainly utilized in fields in which spatial data can be easily obtained from 3D scans of a surface or a volume, such as volumetric imaging in medicine and the representation of terrain in simulations. A voxel, essentially a small axis-aligned cube, is a 3D counterpart to a pixel, and a collection of such elementary cubes constitutes a voxel model of a solid object. Each voxel in the spatial grid is represented by a set of 3D indexes, indicating its positions along three orthogonal directions. Voxel modeling can also be understood as activating specific voxels that belong to the interior of the object. During this voxelization process, active voxels can include internal voxels (entirely inside the object) and surface voxels (partially inside the object), while external voxels (entirely outside the object) should stay inactive. In voxel modeling, the edge length of each voxel element or the finest edge length in a hierarchical data structure decides the available grid resolution of the volume. A higher resolution for a given space usually indicates the voxel model can describe the object with higher precision, but more time is required to manipulate the model.

3.2. Tri-Level Grid

Voxel modeling, with its enumerative nature through spatial indexing, is very suitable for simulating the material removal process in virtual machining as it simply deactivates all voxels inside the cutter-swept volume. It is also sufficiently robust to handle model updates from complex tool motions since it does not have to maintain topological information explicitly. However, for the simplest uniform grid, where the voxel space is uniformly subdivided, the total number of involved voxels is in the order of $\mathbf{O}\left(N^3\right)$, where $N$ is the available grid resolution. As a result, its memory usage and simulation time will rise drastically when a high modeling accuracy is required, and a uniform grid is thus rarely adopted by practical simulations.

In order to mitigate the high memory overhead of a dense uniform grid, an effective approach is to exploit adaptive refinement in spatial decomposition. The resulting hierarchical data structure is a two-level grid, also referred to as a tiled grid, where the volume is first decomposed into uniformly sized cells, also called blocks, bricks, or tiles, and each cell is then partitioned into a finer subgrid if it intersects the surface of the solid object. Note that to quickly reach a target subgrid for a workpiece update, every tile that has been partitioned should be recorded by a key-value container, such as a hash table and map, where its 3D index serves as the key. To simulate the cutting process, an activated voxel at the coarse level that has been swept by the cutting tool entirely needs to be deleted and removed from the container, while an active coarse voxel should be added into the container and partitioned with all finer voxels that fall into the cutter’s interior removed if it is partially cut by the tool for the first time.

The two-level grid manages to reduce both memory usage and computation time significantly; however, the aforementioned shortcoming restricts its further application in high-precision machining. To increase the adaptivity of the two-level grid, a tri-level grid, as shown in Figure 1, is employed, where each bottom-level voxel continues to be subdivided if it intersects with the object surface. Then, the available grid resolution is measured by the product of the branching factors at all three levels. All surface voxels at three levels constitute a surface voxel model of the object, and a solid voxel model, which is used to represent the workpiece model in this study, can be obtained if all internal voxels are also activated. By adding another level of adaptive refinement, the finest voxels that need to be updated can be approached very efficiently.

3.3. Surface Reconstruction

In order to visualize the in-process workpiece or machined part, a boundary representation, such as a triangle mesh, can be built by implementing the marching cubes algorithm [28] based on the voxel model. Triangle mesh has been widely used by many CAM systems due to its strong capability to represent free-form surfaces with simple triangles. Furthermore, by changing the triangle size according to the local curvature of the surface, the storage of the mesh model can be optimized for multi-axis milling simulations. The marching cubes algorithm is an effective approach for visualizing a conceptual surface called an isosurface by extracting its polygonal mesh from a 3D scalar field. An isosurface is formed from a set of 3D points satisfying the equation $f(x,y,z)=c$, where $c$ is a user-specified iso-value and remains constant for any point located on the isosurface. Therefore, an isosurface can be viewed as a surface within a cube where each point has the same parametric value. This algorithm works by iterating (“marching”) over a uniform grid of cubes superimposed over a region of the function, taking eight cube vertices at a time, and then determining the polygons needed to represent the part of the isosurface passing through this cube. The obtained polygons are then fused into the final mesh model. For the first step, an index to a precalculated array of 256 possible polygon configurations ($2^8=256$) is generated for each cube being considered, and each of the eight scalar values corresponds to one bit in an 8-bit integer. A scalar value higher than the user-specified iso-value sets the bit to one, while a lower scalar value sets the bit to zero. Then, the final integer would be the actual index to the polygon indices array, which guarantees any neighboring polygons can be connected appropriately. Finally, each vertex of the extracted polygons, which should be placed in the appropriate position along the voxel edge, is determined by linearly interpolating the two scalar values that are connected by that edge. In this research, the cutting tool is represented by an implicit function, which, however, is not for the computation of the scalar value of each vertex. The use of this function aims to identify every intersecting voxel as well as all voxel edges crossed by the cutter surface, as shown in Figure 2. The edge intersection point is then determined by calculating the real intersection between the voxel edge and the cutter surface.

During the reconstruction process, the polygons inside each cube can be built appropriately according to the lookup table of the marching cubes algorithm as long as all vertex states and edge intersection points are correctly computed. However, a simple data regularization procedure is still required to eliminate noises, inconsistent data, or missing data caused by numerical errors. Specifically, if the distance between a newly computed intersection point and an active voxel vertex is smaller than a specified value, this intersection point is deemed on the vertex whether it is inside or outside the cube. The coordinates of the intersection point are then replaced by those of the vertex, whose state is set to false in the meantime.

4. Representation of Generalized Cutter

A generalized cutting tool can be represented by an implicit surface function, which is used to determine the geometric relationship between an arbitrary 3D point and the cutter surface by simply putting its coordinates into the function. A negative value means the point is inside the cutter, a positive value means the point is outside the cutter, and zero means the point is on the cutter surface. The use of an implicit function is geometrically precise, and its programming is also very easy and straightforward.

According to the definition of an automatically programmed tool, the generalized cutter geometry, as shown in Figure 3, can be described fully using the following geometric parameters [25]: $d$ is the cutter diameter, twice the radial distance between the tool axis and the intersection between lower and upper line segments; $r$ is the radius of the corner torus; $e$ is the radial distance between the center of the corner circle and the tool axis; $f$ is the height of the center of the corner circle measured from the tool tip; $\alpha$ is the angle between the lower line segment and the cutter bottom line, $0°\le \alpha <90°$; $\beta$ is the angle between the cutter side and the tool axis, $-90°<\beta <90°$; $h$ is the overall cutting edge length.

Based on these user-specified independent parameters, several dependent parameters also need to be determined for the construction of the implicit function:

$$r_1=d/2+\left(h-\left(d\mathrm{t}\mathrm{a}\mathrm{n}\alpha \right)/2\right)\mathrm{t}\mathrm{a}\mathrm{n}\beta$$

(1)

$$\begin{array}{c}{r}_2=\left(u+\sqrt{u^2-4{\mathrm{c}\mathrm{o}\mathrm{s}}^2\alpha \left(e^2+f^2-r^2\right)}\right)/2\\ \mathrm{w}\mathrm{i}\mathrm{t}\mathrm{h} u=2{\mathrm{c}\mathrm{o}\mathrm{s}}^2\alpha \left(e+f\mathrm{t}\mathrm{a}\mathrm{n}\alpha \right)\end{array}$$

(2)

$$\begin{array}{c}{r}_3=\left\{\begin{array}{cc}d/2\hfill & \beta =0\hfill \\ e+{\displaystyle \frac{\left(v\mathrm{s}\mathrm{i}\mathrm{n}(2\beta )+\sqrt{v^2{\mathrm{s}\mathrm{i}\mathrm{n}}^2(2\beta )-4(v^2-r^2){\mathrm{s}\mathrm{i}\mathrm{n}}^2\beta }\right)}{2}}\hfill & \beta \ne 0\hfill \end{array}\right.\\ \mathrm{w}\mathrm{i}\mathrm{t}\mathrm{h} v=\left(r_1-e\right)\mathrm{c}\mathrm{o}\mathrm{t}\beta -h+f\end{array}$$

(3)

$$h_{uc}=\left\{\begin{array}{cc}h-f+\sqrt{r^2-{\left(r_3-e\right)}^2}\hfill & \beta =0\hfill \\ \left(r_1-r_3\right)\cot \beta \hfill & \beta \ne 0\hfill \end{array}\right.$$

(4)

$$h_{lc}=r_2\mathrm{t}\mathrm{a}\mathrm{n}\alpha$$

(5)

$$h_t=h-h_{uc}-h_{lc}$$

(6)

The implicit function of a universal cutter consists of four standard components with arbitrary axes. A lower cone, a corner torus, and an upper cone represent the cutting edge, and a cylinder represents the cutter shank. Then, the implicit functions of these four components can be expressed as follows:

$$\begin{array}{c}F(x, y, z)=\left\{\mathit{u}·\mathit{u}-\mathit{m}·\mathit{u}-\begin{array}{c}{D}_{cyl} \mathit{m}·\mathit{u}\ge h\\ D_{uc} { h}_{lc}+h_t<\mathit{m}·\mathit{u}<h \\ D_t h_{lc}<\mathit{m}·\mathit{u}\le h_{lc}+h_t\\ D_{lc} 0<\mathit{m}·\mathit{u}\le h_{lc}\\ D_{bottom} \mathit{m}·\mathit{u}=0\end{array}\right.\\ \mathrm{w}\mathrm{i}\mathrm{t}\mathrm{h}\left\{\begin{array}{c}{D}_{cyl}={r_1}^2\hfill \\ D_{uc}={\left(r_1-\left(h-\mathit{m}·\mathit{u}\right)\mathrm{t}\mathrm{a}\mathrm{n}\beta \right)}^2\hfill \\ D_t={\left(e+\sqrt{r^2-{\left(f-\mathit{m}·\mathit{u}\right)}^2}\right)}^2\hfill \\ D_{lc}=\left\{\begin{array}{cc}{r_2}^2\hfill & \alpha =0\hfill \\ {\left(\mathrm{c}\mathrm{o}\mathrm{t}\alpha \mathit{m}·\mathit{u}\right)}^2\hfill & \hfill \alpha \ne 0\end{array}\right.\hfill \\ D_{bottom}=\left\{\begin{array}{cc}0\hfill & \hfill \alpha \ne 0\\ \left(\left(\mathit{u}·\mathit{u}-\mathit{m}·\mathit{u}\right)+{r_2}^2-\left|\left(\mathit{u}·\mathit{u}-\mathit{m}·\mathit{u}\right)-{r_2}^2\right|\right)/2\hfill & \alpha =0\hfill \end{array}\right.\hfill \end{array}\right.\hfill \end{array}$$

(7)

where $\mathit{m}$ is the unit vector along the cutter axis, and $\mathit{u}$ is the vector from the tool tip to the 3D point being considered.

The use of the implicit surface function aims to detect the contact conditions of a target voxel with the cutter surface for workpiece update by inputting the coordinates of all its vertices into the implicit function. A voxel with all vertices inside the tool is an internal voxel that should be removed, a voxel with all vertices outside the tool is an external voxel that should be preserved, and a voxel with vertices partially inside the tool is an intersecting voxel, which should be decomposed for further updates.

5. Multi-Axis Milling Simulation

The present work employs the method of tool path sampling that takes tool instances along a tool path at a regular interval in order to construct the machined part surface and compute the instantaneous contact region. This method can be applied to a generalized milling tool and any kind of tool path, and it has been widely used for the generation of the cutter-swept volume in many works.

The simulation process proposed by this research is outlined in Figure 4. At each sampled cutter location, the in-process workpiece voxel model is updated by continuously subdividing coarser voxels that intersect the tool instance and removing any cell falling into the interior of the cutter. For each bottom-level intersecting voxel, the edges straddling the cutter surface need to be identified in order to compute and update their intersection points. Meanwhile, the states of all affected vertices need to be updated in time.

5.1. Solid Voxelization of Blank

The proposed algorithm takes the triangle mesh representation of the blank as input because it is relatively simple to create a water-tight solid voxel model, and the intersection between a voxel edge and a planar triangle is easy to compute. In addition, a triangle mesh can be readily obtained from a CSG or B-rep model and easily rendered for visualization. Note that a collection of discrete triangles is utilized to approximate a solid surface by the triangle mesh model; therefore, the approximation error is inevitable during transformation. However, such an error can be mitigated by reducing the average triangle size.

As illustrated in Figure 5, the solid voxelization of the input shape is accomplished by two major procedures: surface voxelization and 3D filling. For the first step, the level-1 voxels intersecting a target triangular face are identified by traversing all voxel candidates inside the axis-aligned bounding box of the triangle. Then, each intersecting voxel is subdivided into a level-2 subgrid, where each intersecting voxel also needs to be activated and subdivided into a level-3 subgrid, with the finest intersecting voxels being activated. For each bottom-level intersecting voxel, all intersection points between its edges and one or more triangles need to be computed for the first time. Once the surface voxel model of the blank is achieved, all internal voxels of the triangle mesh at three levels need to be activated by counting the intersections between the triangle mesh and the ray emitted from the interior of the voxel. An odd number of intersections indicates an internal voxel that needs to be activated, and all its finer voxels, if they exist, need further detection; meanwhile, an even number of intersection points indicates an external voxel.

5.2. Workpiece Geometry Computation

In order to reconstruct the surface of the final machined part, the in-process workpiece model needs to be updated accurately at each sampled cutter location. This multi-level update is accomplished by three levels of voxel updates and one level of edge intersection point update. The level-1 voxel update aims to carry out the fast bulk material removal from the workpiece, which needs to identify all coarse voxels that may possibly intersect the tool instance. A simple method is to traverse every voxel inside the axis-aligned bound box of the tool instance; however, this approach may introduce a large number of unnecessary detection computations, as many voxels are clearly far away from the cutter surface. Therefore, an oriented bounding cylinder is employed by this work, whose internal space can provide the minimal voxel candidates that need to be checked for updates. Note that the interior of the bounding cylinder can be extracted very efficiently by implementing the direct voxel tracing algorithm [24].

Theoretically, the geometric relationship between an affected coarse voxel and the cutter surface needs to be determined by checking all its vertices, which is practically inefficient since some voxels definitely inside or outside the cutter can be quickly ruled out by a relaxed proximity check [27]. Assuming $v_c$ is the coarse-level voxel size and the distance between the voxel center and the cutter surface is $d_{vc}$, then a voxel is deemed to be a near-field voxel, or a possible intersecting voxel, for the tool instance if $d_{vc}<\sqrt{3}{v}_c/2$, where $\sqrt{3}{v}_c$ is the thickness of this near-field region of the cutter surface which lies in the middle of this region, as shown in the dashed boundaries in Figure 6. Once a near-field voxel is recognized, its contact conditions are detected accurately by putting all its vertices into the implicit function of a generalized cutting tool.

Each possible intersecting voxel at the top level is subdivided into a finer level-2 subgrid, where a voxel definitely inside the tool instance is deleted directly, and an intersecting voxel is decomposed into the finest level-3 subgrid. Within each level-3 subgrid, an internal voxel of the cutter is removed, while both voxel edge intersection points and vertex states need to be updated immediately inside each intersecting voxel at the bottom level. Specifically, for an edge with both ends inside the cutter, its two vertex states are set to false, while for an edge with one end inside the cutter and the other outside the cutter, its intersection with the tool instance needs to be updated appropriately. The pseudocode of this multi-level update process is outlined below Algorithm 1:

Algorithm 1. Workpiece Geometry Computation

for each level-1 voxel v1 in oriented bounding cylinder of tool instance

if v1 inside cutter

deactivate v1

else if v1 intersects cutter

subdivide v1 into level-2 subgrid

for each level-2 voxel v2 in subgrid

if v2 inside cutter

deactivate v2

else if v2 intersect cutter

subdivide v2 into level-3 subgrid

for each level-3 voxel v3 in subgrid

if v3 inside cutter

deactivate v3

else if v3 intersects cutter

for each voxel edge of v3

if both vertices inside cutter

deactivate two vertices

else if one vertex inside cutter, the other one outside cutter

if both vertices outside cutter before

update edge intersection point

compute cutting depth and immersion angle

else if both vertex states stay the same as before

if edge intersection point inside cutter

update edge intersection point

compute cutting depth and immersion angle

5.3. CWE Extraction

During the update process of the workpiece at the subvoxel level, the intersection points between the edges of the finest intersecting voxels and the tool instance need to be computed and updated appropriately, which have managed to depict the instantaneous contact surface. Therefore, an accurate and efficient way to compute this region is to fuse these discrete intersection points directly, as shown in Figure 7. In order to obtain the CWE diagram suitable for cutting force prediction, the cutter is uniformly discretized into a group of cylindrical disks first, and their thickness is set equal to the finest voxel edge length for balanced computational performance. Then, for each newly computed intersection point, its cutting depth $d_c$ measured from the tool tip and immersion angle measured from the y-axis of the tool coordinate system (TCS) are calculated, and the newly determined angle is put into all elementary disks that fall into the range of $\left[d_c-\sqrt{3}{v}_f/2,d_c+\sqrt{3}{v}_f/2\right]$, where $v_f$ is the finest voxel size. After all intersection points are determined, the discrete immersion angles stored in each disk are fused by connecting any two neighboring angles with an interval smaller than $2{\tan }^{-1}\left(\sqrt{3}{v}_f/\left(2r_d-\sqrt{3}{v}_f\right)\right)$, where $r_d$ is the radius of the cutting disk. Finally, each group of connected immersion angles in a disk forms a CWE arc, whose first and last angles correspond to the entry and exit angles, respectively. This fusion process achieves the desired resolution and precision of the computed result at a subvoxel level while introducing a very limited number of extra computations.

5.4. CWE Trimming by Envelope Theory

In this work, the computed CWE arc may be slightly larger than its theoretical counterpart because only internal voxels are deleted for workpiece update. Therefore, the so-called feasible contact arc (FCA) is employed in this work to trim the computed result [3], as shown in Figure 8. The grazing curve is a collection of points on the cutter surface that remain on the enveloped surface during the milling process. The entire cutter surface is then partitioned into two parts by this grazing curve. The part facing the cutter moving direction is called the front-facing part, and it is responsible for removing the workpiece material, either in part or in whole; the other part will never be involved. Then, each sliced circle, which is obtained by slicing the cutter with a plane perpendicular to its axis, is divided into two arcs by its intersection points, namely grazing points, with the grazing curve. The arc within the front-facing part is called the feasible contact arc, which defines the largest arc that may engage with the workpiece.

Assume the trajectory of the tool tip can be described by a general curve $\mathit{P}\left(t\right)$ in 3D space, as illustrated by Figure 9; then, the tool coordinate system located at the tool tip can be determined as follows:

$$\left\{\begin{array}{c}{\mathit{o}}_t=\mathit{P}\left(t\right)\hfill \\ {\mathit{x}}_t={\mathit{y}}_t\times {\mathit{z}}_t\hfill \\ {\mathit{y}}_t={\displaystyle \frac{{\mathit{z}}_t\times {\mathit{P}}^{\prime }\left(t\right)}{‖{\mathit{z}}_t\times {\mathit{P}}^{\prime }\left(t\right)‖}}\hfill \\ {\mathit{z}}_t=\mathit{A}\left(t\right)\hfill \end{array}\right.$$

(8)

where $\mathit{A}\left(t\right)$ is the unit axis vector of the cutter and ${\mathit{P}}^{\prime }\left(t\right)$ is the tool tip velocity, which can be computed based on the CL data file:

$${\mathit{P}}^{\prime }\left(t_i\right)={\displaystyle \frac{\mathit{C}\left(t_{i+1}\right)-\mathit{C}\left(t_i\right)}{∆t}}$$

(9)

where $\mathit{C}\left(t_i\right)$ is the tool tip position, and $∆t$ is equal to one.

According to research [3], a 3D point located on the cutting tool surface can be expressed as follows:

$$\mathit{p}\left(z,\theta ,t\right)=\mathit{P}\left(t\right)+z\mathit{A}\left(t\right)+r\left(z\right)\sin \theta {\mathit{x}}_t+r\left(z\right)\cos \theta {\mathit{y}}_t$$

(10)

where $z$ is the axial height with respect to the tool tip, $r\left(z\right)$ represents the cutter radius along the cutter axis, and $\theta$ is the oriented angle measured from the y-axis of TCS.

Furthermore, the intersection point $\mathit{q}$ between the normal vector crossing point $\mathit{p}$ and the cutter axis is expressed as follows:

$$\begin{array}{c}\mathit{q}\left(z,t\right)=\mathit{P}\left(t\right)+l\left(z\right)\mathit{A}\left(t\right)\\ \mathrm{with} l\left(z\right)=r\left(z\right)r^{\prime }\left(z\right)+z\\ r^{\prime }\left(z\right)=\left\{\begin{array}{c}\left\{\begin{array}{c}0 \alpha =0\\ \cot \alpha \alpha \ne 0\end{array}\right. 0<z\le h_{lc}\\ {\displaystyle \frac{f-z}{\sqrt{r^2-{\left(f-z\right)}^2}}} h_{lc}<z\le h_{lc}+h_t\\ \tan \beta { h}_{lc}+h_t<z<h\end{array}\right.\end{array}$$

(11)

Thus, the unit normal vector crossing point $\mathit{p}$ can be written as

$$\mathit{n}\left(z,\theta ,t\right)={\displaystyle \frac{\mathit{p}\left(z,\theta ,t\right)-\mathit{q}\left(z,t\right)}{‖\mathit{p}\left(z,\theta ,t\right)-\mathit{q}\left(z,t\right)‖}}$$

(12)

According to the envelope theory in [13], a grazing point should meet the requirement below:

$${\mathit{v}}_q\left(z,t\right)·\mathit{n}\left(z,\theta ,t\right)=0$$

(13)

where ${\mathit{v}}_q$ is the velocity of point $\mathit{q}$ on the cutter axis, which can be computed by linearly interpolating the velocities of two endpoints on the cutter axis:

$${\mathit{v}}_q(z,t)={\mathit{v}}_1(t){\displaystyle \frac{h-l\left(z\right)}{h}}+{\mathit{v}}_2(t){\displaystyle \frac{l\left(z\right)}{h}}$$

(14)

where $h$ is the overall cutting-edge length, and the velocities ${\mathit{v}}_1$ and ${\mathit{v}}_2$ can be easily determined by connecting the endpoints at two adjacent cutter locations.

Finally, the equation below can be obtained by substituting Equations (12) and (14) into Equation (13):

$$\begin{array}{c}{A}_0\cos \theta +A_1\sin \theta =A_2\\ \mathrm{w}\mathrm{i}\mathrm{t}\mathrm{h} \left\{\begin{array}{c}{A}_0=r\left(z\right){\mathit{v}}_q\left(z,t\right)·{\mathit{y}}_t \\ A_1=r\left(z\right){\mathit{v}}_q\left(z,t\right)·{\mathit{x}}_t \\ A_2=r\left(z\right)r^{\prime }\left(z\right){\mathit{v}}_q\left(z,t\right)·\mathit{A}\left(t\right)\end{array}\right.\end{array}$$

(15)

Based on this equation, the grazing points, and hence the feasible contact arc, on each sliced circle can be represented by $\theta$:

$$\left\{\begin{array}{c}{\theta }_1=2{\tan }^{-1}{\displaystyle \frac{A_1+\sqrt{{A_0}^2+{A_1}^2-{A_2}^2}}{A_0+A_2}}\\ {\theta }_2=2{\tan }^{-1}{\displaystyle \frac{A_1-\sqrt{{A_0}^2+{A_1}^2-{A_2}^2}}{A_0+A_2}}\end{array}\right.$$

(16)

It should be noted that the above equation can be built only if ${A_2}^2\le {A_0}^2+{A_1}^2$. Then, any point on the sliced circle can possibly be a grazing point under the condition of $A_2<-\sqrt{\left({A_0}^2+{A_1}^2\right)}$, but cannot be a grazing point under the condition of $A_2>\sqrt{\left({A_0}^2+{A_1}^2\right)}$.

6. Case Studies

In this section, a series of practical machining cases with sufficient complexity have been conducted in order to compare the computational performance of the proposed method with that of the tri-dexel representation and the two-level grid in different aspects. The tri-dexel model is used as a comparison benchmark because it has been recognized as providing the best combination of computational accuracy, efficiency, and flexibility among all modeling methods.

Table 1. Case summary.

Cutter Type	Ball-End Mill	Bull-End Mill	Taper Ball-End Mill
Cutter parameter (mm)	$d=10$, $r=5$, $e=0$ $f=5$, $\alpha =\beta =0°,$ $h=50$	$d=10$, $r=2.5$, $e=2.5$ $f=2.5$, $\alpha =\beta =0°,$ $h=50$	$d=9.6569$, $r=5$, $e=0$ $f=5$, $\alpha =0°,$ $\beta =2°,$ $h=50$
No. of CL points	33,456	62,600	239,290
Finest voxel size (mm)	0.01	0.05	0.05

lists the geometric information needed to implement these two five-axis milling cases.

6.1. Accuracy Comparison

The computational performance of a voxel-based method highly depends on the voxel size at the bottom level or the ratio of the cutter diameter to the finest voxel size once the cutting tool is selected in the first place. To reveal the computational performance of the proposed method under the condition of extreme-scale resolution, the diameter of the ball-end mill in the first case is thus set to 1000 times the finest voxel size in the initial process. More importantly, a large ratio can help demonstrate how the simulation efficiency of a voxel-based method is influenced by its fine-level branching factor, which can be adjusted in a wider range.

To demonstrate the accuracy of the proposed method, four groups of CWE maps with IPW mesh models of the first case are shown in Figure 10, where the engagement regions computed by the proposed method are drawn in blue while the results determined by the tri-dexel method are drawn in red for comparison. The good similarity, as well as the high complexity of the contact areas, solidly demonstrates the high precision and robustness of the developed method.

The second case is accomplished by a bull-end mill and a taper ball-end mill for its roughing and finishing processes, respectively. Two groups of CWE maps from the roughing process and another two groups from the finishing process are shown in Figure 11, where the 3D contact regions can be seen clearly in close-up images. The good agreement further demonstrates the sufficient accuracy of the present method.

During the cutting process, the sampled tool instances along a tool path are used by this research to approximate the exact cutter-swept volume for workpiece update. As a result, some sampling scallops are left between adjacent tool instances, which can be observed in Figure 10. These scallops, however, are not evident since the distance between the original cutter location points has already been small enough.

Generally, the triangle meshes of the in-process workpiece and the machined part created by the proposed method are all of good quality and thus useful for real-time rendering and visual verification in virtual machining. More importantly, the obtained meshes will become useful when performing a quantitative comparison against their desired models for identifying potential machining errors such as gouging and undercuts. Such good quality, in this work, is due to the small bottom-level voxel size, which is directly associated with the available grid resolution and hence the memory usage. The required memory is mainly used for the storage of the bottom-level voxels with edge intersection points, the containers, and the obtained mesh model. Therefore, the largest grid resolution or workpiece dimensions allowed will be limited by the available memory space of the hardware platform in use.

6.2. Efficiency Comparison

The average processing time, obtained by dividing the overall simulation time by the total number of sampled cutter locations, is listed above each column in Figure 12 for both cases. To compare the modeling methods quantitatively, the execution time of the voxel-based method is divided by that of the tri-dexel method, and the computed result is listed in each column as the efficiency ratio. A smaller ratio indicates a larger improvement in efficiency. In the literature, even though vector modeling is able to accomplish workpiece updates efficiently by implementing clipping operations, it still lacks the ability to minimize these calculations since the workpiece vectors are always evenly spaced. This issue becomes more severe for a tri-dexel model because dexels along all three orthogonal directions need to be considered for possible calculations, which will greatly increase the processing time. Comparatively, because of the two levels of adaptive refinement in spatial decomposition, the proposed method is able to conduct the update process via the collective volume removal by batch processing at the coarse level before performing the finest computations, which effectively accelerates the bulk material removal simulation. With the fast removal of the workpiece material and the identification of near-field voxels, the coarse voxels in the vicinity of the machined part surface can be reached very quickly for finer updates. Then, the time-consuming voxel update at the finest level and the edge intersection point update at the subvoxel level are accomplished only inside the affected coarse-level surface voxels, which reduces the computation time significantly.

Figure 12 also illustrates how the computational efficiency of the two-level grid is affected by its level-2 branching factor. Note that the bottom-level voxel size remains unchanged after being set initially. Thus, different level-2 branching factors correspond to different level-1 voxel sizes. It can be observed that the computation speed of the two-level grid reaches its highest level when the level-2 branching factor is around 8, but it becomes unacceptable once the factor exceeds 32. In contrast, the fluctuations in the computation time of the tri-level grid are rather small, which solidly demonstrates its strong ability to provide high and stable computation speed without being affected by the level-2 branching factor. Essentially, a large level-2 branching factor indicates a large number of top-level voxels will be accessed multiple times for update since a tool path is sampled very closely, which will lead to a sharp increase in traversal operations at the bottom level. As illustrated in Figure 13, where edge intersection computation is neglected, to remove all cutter internal voxels inside level-1 cell $A$ at the current location based on the two-level grid (step 1), all its finer voxels need to be accessed for deactivation, and each access indicates a proximity check and possible accurate detection. By contrast, all level-3 voxels inside cell $B$ are excluded from deactivation since this mid-level voxel has already been deleted during the level-2 update process of the tri-level grid. Once the cutter moves to the next location (step 2), all finer voxels in cell $A$, whose vertices are still inside the tool instance partially, need to be accessed again for the two-level grid, while the finer voxels in cell $B$ and cell $C$, which just falls into the cutter’s interior are skipped when it comes to the tri-level grid. As can be seen, the computational efficiency of the traditional two-level grid is adversely affected by a large level-2 branching factor due to its repeated traversal operations, which, however, can be significantly reduced by adding another level of adaptive refinement.

Figure 14 shows the influence of the level-3 branching factor on the computational efficiency of the tri-level grid quantitatively with three different types of cutters. Since the level-3 voxel size remains unchanged in each case, the level-3 branching factor is reflected by the level-2 voxel size. It can be seen that the computation time of the tri-level grid reaches its lowest point when the level-3 branching factor becomes four or eight, and it rises quickly as the factor keeps increasing, also because of the rapidly increasing bottom-level traversal operations. Therefore, in order to obtain the maximum benefit from the tri-level grid, the level-3 branching factor should be set to four or eight in practice.

7. Conclusions and Future Work

This work presents an effective NC simulation method that can accurately compute machined part geometry and instantaneous CWE regions with higher and more stable computational efficiency under extreme-scale resolutions compared with the traditional two-level grid. The proposed method exploits its potential in high-precision machining by (1) eliminating the adverse impact of a large fine-level branching factor on the computational efficiency by adding another level of adaptive refinement in space decomposition and (2) calculating the immersion angles by fusing the edge intersection points, and trimming the computed results with feasible contact arcs determined by the envelope theory. Furthermore, in order to take full advantage of the tri-level grid, the mid-level subdivision factor should be selected appropriately according to its efficiency characteristic figures.

As discussed in case studies, minor sampling scallops are present on the workpiece surface due to the use of sampled tool instances. In order to obtain the mesh model with sufficient smoothness, the sampling interval needs to remain small, but this will hurt the overall simulation efficiency. Further study is thus needed to optimize the tool path sampling interval according to the fine voxel grid spacing or even eliminate the need to sample tool paths by directly calculating the intersection points between voxel edges and approximate linear cutter-swept volumes. Such a study will help researchers seek a balance between surface quality and computational efficiency in NC simulation.