Curvilinear space-filling curves for five-axis machining

Computer-Aided Design 40 (2008) 350–367 www.elsevier.com/locate/cad

Curvilinear space-filling curves for five-axis machining Weerachai Anotaipaiboon a , Stanislav S. Makhanov b,∗ a Electrical Engineering Department, Faculty of Engineering, Thammasat University, Thailand b School of Information and Computer Technology, Sirindhorn International Institute of Technology, Thammasat University, Thailand

Received 20 March 2007; accepted 18 November 2007

Abstract The paper presents a new combination of two methods for tool path generation for five-axis machining proposed earlier by the authors. The first method is based on the grid generation technologies whereas the second method exploits the space-filling curve approach. Combination of the two techniques is superior with regard to the conventional methods and with regard to the case when the two methods are applied independently. In particular the algorithm allows us to generate tool paths for workpieces with complex boundaries as well as when the scallop and gouging constraints are changing sharply and irregularly. In this case the conventional methods are inefficient, whereas the proposed algorithms construct the required tool path and reduce the length of the path and the time of the machining. The numerical experiments are complemented by the real machining as well as by the test simulations on the Unigraphics 18. c 2007 Elsevier Ltd. All rights reserved.

Keywords: CNC machining; Optimization; Grid generation; Space-filling curves

1. Introduction Milling machines are programmable mechanisms for cutting complex industrial parts. The machines are designed in such a way that the cutting device (the tool) is capable of approaching the desired surface at a given point with a required orientation. The machines consist of several moving parts designed to establish the required coordinates and orientations of the tool during the cutting process. The movements of the machine parts are guided by a controller which is fed with a so-called NC program comprising commands carrying three spatial coordinates of the tool-tip and a pair of rotation angles needed to rotate the machine parts to establish the orientation of the tool. The tool path Π = {Π0 , Π1 , . . . , Πm } is a sequence of positions Π p possibly arranged into a structured spatial pattern such as the zigzag or the spiral pattern. However, modern tool path generation procedures may employ combinations of unconventional patterns and include tool retractions. Let M be a set of parameters defining the configuration of the machine and T parameters which characterize the tool such as the tool length, diameter, shape, etc. ∗ Corresponding author. Tel.: +66 2501 3505; fax: +66 2501 3524.

E-mail address: [email protected] (S.S. Makhanov). c 2007 Elsevier Ltd. All rights reserved. 0010-4485/$ - see front matter doi:10.1016/j.cad.2007.11.007

Furthermore, let S(u, v) = (x(u, v), y(u, v), z(u, v)) be the required surface. The general optimization problem is then formulated by minimize(C), Π ,M ,T

where the criteria vector C includes one or several estimates of the difference between the required and the actual surface. The vector may or may not include the length of the path, the negative of the machining strip (strip maximization), the machining time, etc (see, for instance, [15,17]). Besides, the optimization could be subjected to constraints [16,18,26] the most important of which are: - The scallop height constraints. The scallops between the successive tool tracks must not exceed a prescribed tolerance. - The local accessibility constraints. The constraint ensures against the removal of an excess material when the tool comes in contact with the desired surface due to the so-called curvature interference. - The global accessibility constraints. The constraint ensures against the tool coming in contact with either machine parts (collision detection) or unwanted parts of the desired surface.

W. Anotaipaiboon, S.S. Makhanov / Computer-Aided Design 40 (2008) 350–367

It should be noted that minimization of a particular component of the criteria vector may be replaced by the corresponding constraint and vice versa. For instance, the scallop height could be either minimized or constrained depending on the application. Application of the space-filling curves (SFC) to the NC tool path generation has been first reported in [10,8]. Griffiths [10] proposed the use of the Hilbert’s curve as a tool path while Cox et al. [8] used various forms of space-filling curves such as the Moore’s curve. Both tool path generation methods have been developed for three-axis NC machining with the ball-end cutter. However, the space-filling curves have not been particularly popular in the five-axis machining community due to large inaccuracies produced by sharp angular turns which characterize the standard SFC patterns. Nonetheless, tool path generation based on the SFC has a number of attractive features such as the possibility to locally adapt the curve without changing the global structure. In addition, the entire surface is cut in one path eliminating the need of tool retractions. Recently the authors proposed a new tool path generation method using the adaptive SFC which overcomes the above mentioned drawbacks while keeping the above mentioned advantages [1,2]. The proposed SFC tool path generation method requires four steps: -

Construction of a basic rectangular grid; Generation of the adaptive space-filling tool path on the grid; Correction of the tool path; Inserting additional points along the path to reduce the kinematics error.

First of all, the corresponding parametric region in the u, v-plane is covered by two zigzag tool paths. One path is constructed in the u-direction and another one in the vdirection. The step between two adjusted parallel tracks of the zigzag paths is evaluated in such a way that that the machining strips overlap and the scallops are within the allowable margins. The evaluation step requires as an input the surface curvature, the cutter shape and size, and the maximum allowable scallop height. The output of the procedure is the minimal tool inclination angle and the size of the maximum possible machining strip which defines an offset to the adjusted parallel segment of the path (track). Furthermore, the grid created by overlaying the two zigzag tool paths is called the basic grid. The SFC generation step treats the cells of the basic grid as nodes and connects them by arcs. The SFC is constructed as a Hamiltonian path on the such a grid-like graph using a cover and merge algorithm [2,20]. Finally, the kinematics error is reduced by inserting additional points along the resulting space- filling curve until the error does not exceed a prescribed tolerance. It is plain that once the maximum achievable machining strip is evaluated, constructing the basic grid for the SFC generation is trivial. However, the rectangular grid is often inefficient since a small step between the tracks could be required only in certain areas. The grid is also inefficient in treating complex geometries appearing in the case of the so-called trimmed surfaces having the boundaries created by intersections with other surfaces. The

351

complex boundaries also occur in the case of pocket milling when the parametric region includes internal boundaries around one or several pockets (islands). Furthermore, the above geometrical complexities and sharp variations of the surface curvature have been proved to be successfully treated by numerically generated curvilinear zigzag tool paths obtained from adaptive grids topologically equivalent to the rectangular grids. In [21,23] a modification of a classic grid generation method based on Euler–Lagrange equations for the Winslow functional [25] has been adapted to the zigzag tool path generation. The zigzag tool path is constructed by solving numerically Euler–Lagrange equations for a functional representing desired properties of the grid such as smoothness, adaptivity to the boundaries and to a certain weight (control) function [25,5]. A similar idea to use a Laplacian curvilinear grid was suggested independently by Bieterman and Sandstrom in 2003 [6]. In the framework, proposed in [21] the grid is adapted to the kinematics error subject to constrains relevant to the heights of scallops between the successive tool tracks. However, the above techniques have a number of drawbacks. In particular, they may converge slowly for complicated constraints. Besides, the approach requires equal number of the cutter contact points on each track of the tool. Therefore, if the kinematics error changes sharply from track to track, the method may require an excessive number of points. This paper introduces a new modification of the grid refinement which fits better in the framework of tool path optimization and is designed specifically for the SFC generation. The method does not require an equal number of points on each track since the grid generation does not try to minimize or to reduce the kinematics error along the tracks as in [21]. At the grid generation stage the kinematics error is ignored and the grid is generated with regard to the scallop height constraints. Next, the grid can be converted into a curvilinear zigzag path or into a space- filling curve. Only after that the cutter location points along the resulting curve are distributed using a standard method such as the bisection. Additionally, the algorithm automatically evaluates the number of the required grid lines. As opposed to the preceding approach, where the weight function represents either the kinematics error or an estimate of the kinematics error (such as the surface curvature or the rotation angles), the proposed algorithm iteratively constructs an adaptive control function designed to represent the scallop height constraints. This important modification makes it possible to consider an arbitrary number of points along the tool tracks. The kinematics error is reduced by means of inserting additional points along the resulting curvilinear coordinates. In other words, this approach replaces the scallop constraints by a weight function and then treats the kinematics error independently. Additionally, instead of the Winslow functional the new optimization is based on the harmonic functional derived from the theory of harmonic maps [12,13]. The functional not only provides the smoothness and the adaptivity but under certain conditions guarantees the numerical convergence [13]. Finally, this approach merges with the SFC techniques. In this case,

352

W. Anotaipaiboon, S.S. Makhanov / Computer-Aided Design 40 (2008) 350–367

the grid is not converted to the tool path directly. Instead, it becomes the basic grid required for the SFC generation. With this modification, the SFC tool path can be constructed for surfaces with complex irregular boundaries, cuts off, pockets, islands, etc. Besides, the adaptive grid allows to efficiently treat complex spatial variability of the constraints in such a way that the SFC is created on a grid having the small cells only where necessary. Finally, the combination of the two techniques is superior with regard to the case when the two methods are applied independently. We present a variety of examples when the conventional methods are inefficient whereas the proposed algorithms allow for constructing the required tool path and reduce the length of the path and the time of the machining. We also show that the combination of the two techniques could be superior with regard to the case when the two methods are applied independently. The numerical experiments are complemented by the real machining as well as by the test simulations on the Unigraphics 18. 2. An introductory example of a curvilinear space-filling curve Let us demonstrate that the combination of the two techniques may decrease the tool path. Consider the following minimization with regard to the length of the tool path min L , Π

subject to h < h max , ε < εmax ,

(1)

where h is the scallop height, measured between each pair of the adjusted tracks, and h max is the maximum allowed scallop height. If the machining strips do not overlap, the remaining areas are considered as the scallops as well. Therefore, the first constraint controls the scallops and provides that the machining strips cover the entire surface. The second constraint requires that the difference ε between the actual and the desired trajectory extracted in a certain way from the desired surface, does not exceed a prescribed tolerance ε max . Furthermore, we shall call ε the kinematics error and define D D D D it as follows. Let W p, p+1 (t) ≡ (x p, p+1 (t), y p, p+1 (t), z p, p+1 (t)) ∈ S(u, v) be a curve between two tool positions W p and W p+1 extracted from the machined surface S(u, v), t is the parametric coordinate along the curve. The curve is extracted in such a way that it represents the desired tool trajectory. The kinematics error is defined as a norm of the difference D between W p, p+1 (t) and the actual trajectory W p, p+1 (t) ≡ (x p, p+1 (t), y p, p+1 (t), z p, p+1 (t)), namely, XZ 1 2 D |W p, ε= p+1 − W p, p+1 | dt p

=

0

1

XZ p

0

D 2 D 2 (x p, p+1 − x p, p+1 ) + (y p, p+1 − y p, p+1 )

2 + (z D p, p+1 − z p, p+1 ) dt.

Fig. 1. The unimodal surface which exponential peak along a line. D The actual trajectory W p, p+1 (t) is generated using the inverse kinematics equations derived for the specific machine configuration (see, for instance, [3]). Finally, the two conditions in (1) are nothing else than a numerical approximation of a general constraint that the machined surface does not deviate from the actual surface by more than the prescribed tolerance Tmax = max(h max , εmax ). However, Tmax characterizes only an approximation to the actual cut, suffice to say that h max , εmax are measured at different points. As a matter of fact, the actual error might exceed Tmax , however, h max and εmax are the conventional error estimates which should be minimized first. Consider an introductory example. Let us perform minimization (1) by choosing the positions and inclinations Π p ∈ Π in such a way that the machining strip width is maximized. We will show that the resulting curvilinear SFC (CSFC) makes it possible to cut a surface by travelling along a shorter path. Consider a unimodal surface which an exponential peak along a line shown in Fig. 1. Let us use a flat end tool with the radius r = 3 mm. In this case the surface requires a vertical zigzag (isoparametric) tool path of about 5952 mm (Fig. 2(a)). The SFC tool path requires a basic grid with a small spatial step. Consequently it leads to a SFC of about 4144 mm (Fig. 2(b)). The zigzag tool path based on the adaptive grid is shown in Fig. 2(c). The length of the tool path is 4348 mm. Finally, the CSFC on that grid requires a shorter path of about 3814 mm (Fig. 2(d)). The step between two adjusted tracks on the rectangular grid is set to the minimal step ensuring the required scallop. As opposed to that, the curvilinear grid makes it possible to vary the distance between the tracks making the step small only where necessary. Therefore, the CSFC constructed within a new basic grid leads to better results.

3. Grid generation method Recall that S(u, v) ≡ (x(u, v), y(u, v), z(u, v)) denotes the surface to be machined. As usual u and v are the parametric variables. Consider a set of cutter location points {u i, j , vi, j } arranged as a curvilinear grid. Mathematically, it means that (u i, j , vi, j ), 0 ≤ i ≤ Nξ , 0 ≤ j ≤ Nη is a discrete analogy of a mapping from the computational region {0 ≤ ξ ≤ Nξ , 0 ≤ η ≤ Nη } onto a parametric region defined in the parametric

W. Anotaipaiboon, S.S. Makhanov / Computer-Aided Design 40 (2008) 350–367

353

Fig. 2. (a) Iso-parametric tool path. (b) SFC tool path from two iso-parametric toolpaths, (c) zigzag tool path from adaptive grid. (d) SFC tool path from adaptive grid.

3.1. The harmonic functional

Fig. 3. Coordinate transformations and curvilinear grids.

coordinates u, v. In other words, there exists a pair of functions {u(ξ, η), v(ξ, η)} such that the rectangular grid {i, j} being fed to {u(ξ, η), v(ξ, η)} becomes {u i, j , vi, j }, see Fig. 3.

The required (unknown) grid is a discretized solution of the minimization problem given in Box I, where subscripts u, v, ξ, η denote partial derivatives and f the control function. The harmonic functional I is a generalization the Winslow functional to the case of the grids lying on the surface of control function f (u, v). The harmonic functional is derived from the theory of harmonic maps [13]. It has been proved that the functional minimizes an “energy of mapping” [24] and produces a grid adapted to the regions of large gradients of f . Note that if f u = f v ≡ 0, then the harmonic functional becomes the Winslow functional [25], however, it is important that I adapts the grid to the gradients of f rather than to f itself as in the previous version [21]. It is known that the minimization could be computationally expensive as compared with minimization of the Winslow functional [7]. However, it has many points in its favour. In particular, it is possible to construct a computational procedure which, under certain conditions, converges to a non degenerate grid [24], that is, the grid without twisted or non convex cells.

354

W. Anotaipaiboon, S.S. Makhanov / Computer-Aided Design 40 (2008) 350–367

Z (u 2 + u 2 )(1 + f 2 ) + (v 2 + v 2 )(1 + f 2 ) + 2 f f (u v + u v ) u v ξ η η ξ η u η v ξ ξ p dξ dη, min I ≡ min 2 2 u,v u,v (u ξ vη − u η vξ ) 1 + f u + f v Box I.

The constraint minimization of the functional in Box I can be performed by using efficient penalty type techniques similar to those presented in [21]. Finally, functional is more reliable and converges for sharp variations of the input data whereas the Winslow functional often produces degenerated grids. Note, that the functional in Box I is just one from many possibilities offered by grid generation methods chosen for its convenience, efficiency and robustness. Other grid generation schemes could be also utilized for problems of tool path planning. 3.2. The data structure We introduce the following data structure suitable for regular, block structured grids or even irregular grids. Instead of double indexing {i, j}, the grid nodes are characterized by the global node numbers n = 1, . . . , Nn . A grid element, be it triangle, quadrilateral or other shape, is characterized by the element number N = 1, . . . , Ne . Finally, the nodes inside each element are numbered using local node numbers. Function C 0 (N , k) defines the correspondence between the local and the global node numbers as follows C 0 (N , k) = n,

where k = 1, . . . ke .

(2)

We will consider quadrilateral grids, ke = 4. The element (cell) vertices are numbered from 1 to 4 in the clockwise direction. Additionally, each vertex is associated with a triangle: vertex 1 with triangle 421, vertex 2 with triangle 123, and so on. It can be shown that the Jacobian of the mapping Jk , k = 1, 2, 3, 4 at node k can be approximated by two times the area of the corresponding triangle. Finally, we need function N 0 (n, d) which for each node n returns its neighbours in the direction d, where d ∈ {left, right, up, down}. An example of such data structure is shown in Fig. 4.

Fig. 4. Correspondence of node numbers for a mapping of the unit squre in the (ξ, η) plane on to the quadrilateral cell 1 of the grid in the (u, v) plane.

Here ( f u )k and ( f v )k are the derivatives of the control function in the u and v directions, respectively, computed at node k of the cell N . The requirement of the grid cells convexity introduced by [12] corresponds to the positivness of the Jacobian [Jk ] N > 0.

(4)

A grid for which (4) is satisfied for every N and k is called the convex grid. 3.4. The solution method

3.3. Approximation The functional in Box I is approximated by a discrete functional as follows: I ≈ Ih =

Ne X 4 X 1 N =1 k=1

4

[Fk ] N ,

(3)

where Fk =

D1 [1 + ( f u )2k ] + D2 [1 + ( f v )2k ] + 2D3 ( f u )k ( f v )k Jk [1 + ( f u )2k + ( f v )2k ]1/2

,

D1 = (u k−1 − u k )2 + (u k+1 − u k )2 , D2 = (vk−1 − vk )2 + (vk+1 − vk )2 , D3 = (u k−1 − u k )(vk−1 − vk ) + (u k+1 − u k )(vk+1 − vk ).

Discrete functional (3) is minimized numerically using an initial grid for the first iteration. The minimization is performed in such a way that the minimum is attained on a grid composed of convex quadrilaterals. This property improves the stability of the grid generation procedure. It can be shown that if the set of convex grids corresponding to the mapping is not empty, the system of algebraic equations Ru =

∂Ih = 0, ∂u n

Rv =

∂Ih =0 ∂vn

has at least one solution which minimizes (locally) discrete functional (3). Suppose that the grid at the l-th step of the iterations is calculated. Then, the grid at the (l + 1)-th step is computed

355

W. Anotaipaiboon, S.S. Makhanov / Computer-Aided Design 40 (2008) 350–367

using a quasi-Newtonian procedure as follows:   ∂ Rv ∂ Ru l+1 l u n = u n − τ Ru − Rv ∂vn ∂vn   ∂ Rv ∂ Ru −1 ∂ Ru ∂ Rv − , × ∂u n ∂vn ∂u n ∂vn   ∂ Ru ∂ Rv vnl+1 = vnl − τ Rv − Ru ∂u n ∂u n   ∂ Rv ∂ Ru −1 ∂ Ru ∂ Rv , × − ∂u n ∂vn ∂u n ∂vn

(5)

where τ is the iteration parameter, which is chosen so that the grid remains convex. Namely, at each step l, condition (4) is verified. If (4) is not satisfied, then τ is reduced by half and the grid at (l + 1)-th iteration is re-computed. 3.5. Control function for tool path optimization Let us now describe the use of the harmonic functional for tool path generation. Our control function f represents the distance between the tool tracks so that the resulting grid guarantees the required height of the scallop. The maximum allowable distance (the machining strip) T between the cutter location points is evaluated by finding numerically the intersection of the tool effective cutting shape and the machined surface. The cutting shape of the tool depends on the tool inclination which must be chosen so that the curvature interference is eliminated (see the details in Appendices A and B. Since the tool path is a discrete set of points, the derivatives of the control functions f u (u, v) and f v (u, v) for a given surface can not be explicitly evaluated. Therefore, we generate these derivatives “artificially” as follows. First ( f u )0n ≡ ( f v )0n ≡ 0. Next,  ( f u )ln + λ+ , if su (n) > 0, = ( f u )l+1 n ( f u )ln − λ− , otherwise,  ( f v )ln + λ+ , if sv (n) > 0, ( f v )l+1 = n ( f v )ln − λ− , otherwise , where λ+ and λ− is the prescribed increment and decrement, respectively, and sd (n) is the difference between the actual distance between the tracks and the machining strip evaluated at node n as follows: su (n) = sv (n) =

max

W (n, N 0 (n, d)) − T (n, N 0 (n, d)),

max

W (n, N 0 (n, d)) − T (n, N 0 (n, d)),

d∈{left,right} d∈{up,down}

where W (n, m) is the distance between nodes n and m given by W (n, m) = |S((u, v)n ) − S((u, v)m )|, and T (n, m) is an
estimate of the machining strip at midpoint m vn +vm S u n +u , (see Appendices A and B). 2 2

Fig. 5. Partitioning of block-structured grid for refinement.

points must be inserted for convergence. For a structured grid having initially nr rows and n c columns, the number of rows and columns at the next step is evaluated as follows. nr,new = nr + n add,r , n add,r max

1≤i≤nr

=

nP c −1

W n i, j , n i, j+1 − T n i, j , n i, j+1


!


j=1

2r

(6)

n c,new = n c + n add,c , n add,c max

=

1≤ j≤n c

nP r −1

!

W n i, j , n i+1, j − T n i, j , n i+1, j

i=1

2r

,

where n i, j is the grid node and r is the tool radius. If the parametric region in the computational domain is not rectangular (Fig. 5(a)), then the grid is partitioned vertically as shown in Fig. 5(b) and horizontally as shown in Fig. 5(c). It should be noted that if the grid is constructed to produce a curvilinear zigzag tool path in one direction, then only one from the two formulas (6) must be applied. However, if the grid is needed for the CSFC generation, then (6) must be applied in the both directions. Finally, (6) may overestimate the number of the required tracks. Consequently, it can be replaced by nr,c,new = nr,c + αrel n add,r , where αrel < 1 is a “the rate of release” of the additional curves. Such a procedure may lead to a decrease in the number of the zigzag curves, thus, improving the efficiency of the machining. The “rate of release” is determined experimentally. 1 An inexperienced user will be safe with αrel = n add,r which, however, may lead to the largest computational cost. 3.7. The grid generation algorithm

3.6. Inserting additional tracks

The grid generation algorithm consists of the following steps:

The initial grid does not (and should not) satisfy the scallop height constraint. However, it is often the case where additional

1. Generate an initial convex grid. The grid is generated manually or by interpolating. Note that interpolation may

356

W. Anotaipaiboon, S.S. Makhanov / Computer-Aided Design 40 (2008) 350–367

Fig. 6. Adaptive harmonic tool path generation (a) the surface, (b) the control function, (c) the grid.

generate a grid with the nodes outside the boundary of the region. In that case, several iterations can be performed by the classic Brackbill’s and Saltzman’s method [5] which will move the nodes back inside the region. The initial grid should not satisfy the scallop constraints, because, if it does, the adaptation is not necessary. It also means than the number of grid nodes could be reduced. 2. Insert additional grid lines (tool tracks) using (6). 3. Adapt the grid by minimizing functional (3) until all the grid points satisfy the scallop constraint or until a prescribed number of iterations has been exceeded. 4. If the scallop constraint has not been satisfied for all the points, goto the refinement stage 2. Fig. 6 illustrates the adaptive harmonic grid generation applied to a simple surface characterized by a large gradients along a sinus-shaped zone (Fig. 6(a)). The required small machining strips generate the control function depicted in Fig. 6(b) which in turn produces a grid adapted to the control function depicted in Fig. 6(c). Finally, our curvilinear grid is nothing else than two overlapped curvilinear zigzag isoscallop paths. We may stop at that point and use one of them or go further and move inside them using the space-filling curve technology presented next. 4. Tool path generation using curvilinear space-filling curve method Minimization of the length of the tool path is formulated as follows. Given the grid of points, connect them by constructing a continuous curve with the minimal length on the design surface. In this chapter we will present a solution of the problem based on space-filling curve technologies. A space-filling curve is defined as a continuous mapping of a unit line segment onto the unit square. In the limit, the SFC passes through every point in the two-dimensional unit square, even though it is a one-dimensional entity. The most popular is the recursive Hilbert’s curve [14] considered for numerous applications including the tool path planning [10]. Hilbert’s curve is particularly appealing in tool path planning as its local refinement property can be used to adaptively increase the density of the path only where necessary. However, each

local refinement of the tool path based on the Hilbert’s curve increases the tool path density in the refined region by a factor of 2 resulting in lower machining efficiency due to the increased total path length. Besides, the Hilbert’s curve has an undesirable property that it leads to a path where the tool is constantly changing directions which slows down the machining process and produces large kinematics errors. To overcome these drawbacks [2] proposed to use adaptive SFC characterized by the following features. First of all, the adaptive SFC always follows the local optimal direction. Second, as opposed to the conventional SFC, the adaptive SFC turns only when necessary, in other words, only when the optimal direction changes. Third, the adaptive SFC eliminates the large kinematics errors and the overcuts appearing due to the sharp angular turns. Finally, local refinement of the adaptive SFC is accomplished in exactly the same fashion that the conventional SFC is refined. In this paper we will adapt this algorithm to the case of curvilinear block structured grid. 4.1. Space filling curve generation As an example consider a basic curvilinear grid generated using the algorithm developed in Section 3 (Fig. 7(a)). In order to construct the CSFC, each grid cell is replaced by a vertex in the middle of the cell (Fig. 7(b)). Every pair of adjacent vertices is then connected by an edge as shown in Fig. 7(c) and (d) to create an initial set of small circuits. Note that vertices of the graph correspond to an initial set of CC points on the required surface. Therefore, the distance between two connected vertices is defined as the distance between the corresponding CC points on the surface in R 3 . A cut along the path between any two connected vertices satisfies the scallop height constraint. This feature allows for the tool path optimization by means of the SFC. The SFC tool path generation algorithm is presented next. The tool path generation on the grid-like graph is a particular case of the travelling salesman problem called the Hamiltonian path problem [19]. Since the problem is NP-hard [11], the algorithms for finding the optimal solution are slow and inefficient. We present a simple and computationally efficient algorithm for producing the Hamiltonian path based on the cover and

W. Anotaipaiboon, S.S. Makhanov / Computer-Aided Design 40 (2008) 350–367

357

Fig. 7. (a)–(c) Construction of an undirected graph for the CSFC tool path generation. (d) Generation of an initial set of circuits for constructing the CSFC.

merge algorithm developed by [20] for 2-dimentional image scanning. In this paper we extend the algorithm for non rectangular domains and block structured grids and apply it to CSFC tool path generation. The CSFC tool path generation works as follows. First, all vertices are covered by small disjoint circuits. The circuits are then merged into a single Hamiltonian circuit. The initial circuits are created by constructing small rectangular cyclic paths over every four adjacent vertices, i.e. by connecting the vertices on even rows and columns with the vertices on odd rows and columns, respectively, as shown in Figs. 7(d) and 8(a). For structured grid, if the grid size is odd, we construct the virtual circuits to cover the vertices along the boundaries as shown by the dashed lines in Fig. 8(b). Any two adjacent circuits can be merged into one bigger circuit. The cost of merging is defined by Cost(A, B) = |s| + |t| − |e| − | f | ,

(7)

Fig. 8. (a) Undirected graph G covered by initial small circuits (b) dual graph G0.

where |e| denotes the distance between two vertices connected by edge e.

358

W. Anotaipaiboon, S.S. Makhanov / Computer-Aided Design 40 (2008) 350–367

Note that the distance | | along the edges is measured on the surface. Therefore, by minimizing the cost function, for every pair of the adjustent circuits the algorithm attempts to minimize the total length of the tool path on the surface. The cost of merging of two virtual circuits is set to −∞, i.e. all the virtual circuits are initially merged. This is to ensure that there is no discontinuity of the tool path after removing the virtual edges from the Hamiltonian path. Furthermore, a nonvirtual circuit A can be merged with a virtual circuit D only if A is merged with a nonvirtual circuit C located on the opposite side. To enforce this merging dependency, the cost of merging A and D is set equal to that of merging A and C. The merging dependency is used to eliminate a possibility of an inappropriate narrow zigzag tool path with a large number of turns along the boundaries. To merge all small circuits, we construct a dual graph as shown in Fig. 8(b) and a minimum spanning tree is constructed by iteratively merging circuits according to the cost defined by (7). After all the circuits are merged into a Hamiltonian circuit, the CSFC is generated by removing all the virtual edges. Recall that the tool orientation is changing dynamically in order to eliminate the curvature interference. To ensure against the curvature interference a modification of a standard gouging avoidance algorithm is used (see Appendix B). In turn, the cost functions f u (u, v) and f v (u, v) in Box I include the machining strip calculated at each point of the basic curvilinear grid (Appendix A). Therefore, the CSFC is constructed within the grid for which the scallop height constraint h < h max is satisfied for every cutter location point. In order to improve the CSFC at the corners a correction (smoothing) procedure (Appendix C) is used. The details can be found in [2]. The procedure introduces two corrections at the corners of the CSFC. The first correction eliminates undercut areas, the second smoothes the tool orientation when the tool is sharply changing the direction. Observe that the direction of the CSFC at a given point can be selected using a variety of criteria such as the speed of the tool, acceleration of the linear or/and rotation axis (for instance, for high speed machining), etc. In this case the definition of the distance between the two cutter location points must be modified accordingly. Of course, having angular accelerations constrained could turn the SFC into a curvilinear zigzag (which is a simple SFC as well). In this case only the grid generation part will lead to an improvement. In some cases it is possible to perform optimization directly with regard to the machining time, however, for the general machining tool it is still an open problem. Furthermore, the cost function may be also adapted to the case when the direction of the SFC is fixed at certain areas. For example, the blades of industrial impellers or propellers may be produced given preferable directions of the scallops at some prescribed areas to provide better hydrodynamics conditions. Finally, our method has been implemented and verified for conventional (“slow”) five-axis machines which are still important in manufacturing. An adaptation and efficient implementation of our algorithm to high speed machining is still an open problem.

Table 1 Comparison of the methods in terms of the tool path length Method used for constructing the basic grid

Tool path length (mm) Vertical zigzag SFC

Isoparametric tool path Curvilinear grid generation

5952.47 4348.10

4131.18 3813.57

5. Examples and practical machining This section demonstrates the efficiency and advantages of the use of curvilinear grid in tool path generation by examples and practical machining. Example 1. A unimodal surface The first example demonstrates the efficiency of the CSFC with the reference to the traditional isoparametric tool path method. Consider a unimodal surface which exponential peak along a line in the parametric domain (u, v) given by (see Fig. 1) x = 100u − 50, y = 100v − 50,

(8)

z = 10 exp(−40(2u − 0.5 − v) ) − 15. 2

For flat-end tool of radius 3 mm and machined surface tolerance Tmax = 0.01 mm the final curvilinear grid and SFC are shown in Fig. 9. The comparison of the zigzag and SFC tool paths generated from traditional isoparametric tool path and curvilinear grid is presented in Table 1. The length of the zigzag and SFC tool paths based on the adaptive grid are shorter by 27% and 8%, respectively, when compared with the zigzag and SFC based on isoparametric tool path method (Fig. 2). Besides the SFC curvilinear tool path is shorter than the curvilinear zigzag by approximately 12%. Of course, in five axis machining the shortest path does not necessarily mean the shortest time. Therefore, Table 2 compares the machining time of surface (8) for different angular and linear feed rates. Symbol * in the table indicates the time measured in a real machining experiment. For other combinations of the feed rates we use a simple estimate based on the time required for the slowest (for this step) axis. Given the spatial and angular increments ∆x M , ∆y M , ∆z M , ∆a, ∆b and the feed rate F, the machining time t M is calculated by t M = max{t0 , tx , t y , tz , t A , t B }, where t0 = ∆FL ,where p xM ∆L = (∆x M )2 + (∆y M )2 + (∆z M )2 and where tx = v∆x,max , yM yM A B , tz = v∆z,max , t A = v∆ , t B = v∆ , where tx = v∆y,max A,max B,max vmax denotes the maximum speed in the corresponding axis. For our experiments on MAHO600E the angular speed is given by v A,max = 235 ◦ /s, v B,max = 162 ◦ /s and vx,max , v y,max , vz,max are the varying feed rates given in Table 2. Since SFCs require many turns, a slow angular feed rate may lead to a time longer than that provided by the curvilinear grid. However, the SFC always overperforms the rectangular grid. Besides, for fast angular feed rates and short, slow linear motions the curvilinear SFC is a very good choice. Tables 1 and 2

359

W. Anotaipaiboon, S.S. Makhanov / Computer-Aided Design 40 (2008) 350–367

Fig. 9. Curvilinear grid adapted to the unimodal surface which exponential peak along a line in (a) u − v domain and (b) grid on the surface in the workpiece coordinate systems, (c) the corresponding curvilinear space-filling curve.

Table 2 Comparison of the methods in terms of the machining time Tool path

Feedrate 100 mm/min

Vertical zigzag based on isoparametric tool path SFC based on isoparametric tool path Vertical zigzag based on curvilinear grid generation SFC based on curvilinear grid generation

224.5 min∗ 152.0 min 134.5 min∗ 121.0 min∗

also show that the use of the curvilinear grid always produces the machining time better than the conventional method. Unfortunately, constructing a SFC which is designed to minimize the machining time or a distance in the angular space (a, b) for an arbitrary tool cannot be done by merely replacing cost function (7). The direction of the tool path at a given point is changing dynamically during the optimization and there is no a priori knowledge in which direction the SFC is going to turn. However, such optimization can be done for the ball nose tool which can be positioned irrespectively of the tool path direction.

500 mm/min

1000 mm/min

4000 mm/min

45.00 min 30.5 min 27.00 min 24.2 min

22.5 min 15.0 min 13.5 min 12.0 min

5.6 min 3.8 min 3.4 min 3.0 min

Example 2. Curvilinear boundaries and pocket milling. The second example demonstrates the use of the CSFC to construct tool paths to machine surfaces with complex irregular boundaries, cuts off and islands. Consider a surface shown in Fig. 10(a). Fig. 10(b) shows the basic curvilinear grid constructed using the proposed method. Fig. 10(c) shows the CSFC and Fig. 10(d) the CSFC on the surface. Finally, Fig. 10(e) shows the machining result obtained with the use of the solid modelling engine of the Unigraphics. The surface has been machined by a flat-end tool of radius 3 mm with surface tolerance

360

W. Anotaipaiboon, S.S. Makhanov / Computer-Aided Design 40 (2008) 350–367

Fig. 10. Constructing the curvilinear space- filling curve. (a) The surface, (b) the curvilinear grid, (c) the space-filling curve, (d) the space-filling curve on the surface, (e) virtual machining by Unigraphics 18.

W. Anotaipaiboon, S.S. Makhanov / Computer-Aided Design 40 (2008) 350–367 D1

Fk =

D2

361

|vk,k−1 ||vk,k+1 | [1 + ( f u )k ] + |vk,k−1 ||vk,k+1 | [1 + ( f v )k ] + 2( f u )k ( f v )k cos ϑ . Jk 2 + ( f )2 ]1/2 [1 + ( f ) u v k k |v ||v | 2

k,k−1

2

k,k+1

Box II.

Fig. 11. Convergence of the method with and without the “angular” term for Example 1.

Tmax = 0.05 mm. Consequently, the method is capable of creating tool path for surfaces with complex nonrectangular boundaries and islands.

Finally, although, the improvement in the convergence rate is not evidently large, the complexity of the algorithm is reduced when the “angular” term is eliminated.

Remark. Observe that the algorithm converges slightly faster and exhibits better stability if the third term in the nominator of functional given in Box I is neglected. In order to interpret this term geometrically for each Fk in (3) consider two vectors from local node k to two adjacent local nodes k − 1 and k + 1 on the same cell given by

Example 3. Point milling of an impeller blade.

vk,k−1 = [u k−1 − u k , vk−1 − vk ], vk,k+1 = [u k+1 − u k , vk+1 − vk ].

(9)

Clearly, vk,k−1 · vk,k+1 = (u k−1 − u k )(u k+1 − u k ) + (vk−1 − vk )(vk+1 − vk ) = D3 = vk,k−1 vk,k+1 cos ϑ (see page 8), where ϑ is the angle between the two vectors. Thus, we have the equation given in Box II. Therefore, omitting the last “angular” term leads to an adaptation with less impact on the angles between the coordinate curves. Figs. 11 and 12 display the convergence of the method given in terms of the number of the “bad points”, that is, the points where the scallop height requirements are not satisfied. The dynamics of the bad points is given in terms of a number of the points and the percentage. Note that at certain iterations the absolute number of the bad points may increase due to the deformations of the grid as well as due to inserting additional curves at step 2.

The impeller blade is one of the most challenging benchmarks of 5-axis machining. Many industrial parts such as the car body, various moulds and dies which seem to be more complicated than the impeller can be produced in the 3axis mode. However, the impeller blade can be produced in one setup only in the 5-axis mode. Therefore, the applicability to the impeller blades is a very good test for five-axis algorithms. Frequently, the blades are produced by the so-called fiveaxis swarf milling made by a side of the tool. In this case the contact between the workpiece and the cutter is characterized by a contact line rather than a contact point. However, the cutter contact line may lead to large errors. Our example demonstrates machining of the impeller by 5-axis point milling with the use of curvilinear adaptive space-filling curves. Fig. 13 shows a blade produced from the hard aluminium by the proposed method on MAHOO600E. In order to demonstrate the advantages of the proposed method, we will additionally increase the geometrical complexity of the example blade as follows. Consider a blade depicted in Fig. 14(a). After long-hours serves in a harsh environment the blades may suffer from a variety of defects, such as distortion, cracks, nicks and dents. Suppose that the blade is broken as shown in Fig. 14(a) (the dashed line) and requires a restoration. The missing part is shown in

362

W. Anotaipaiboon, S.S. Makhanov / Computer-Aided Design 40 (2008) 350–367

Fig. 12. Convergence of the method with and without the “angular” term. Example 2.

Fig. 14. (a) A broken blade. (b) The missing part of the blade.

Fig. 13. The impeller blade machined from a hard aluminum using the proposed method.

Fig. 14(b). Note that similar (but smaller in size) restorations through the reverse engineering techniques are described, in [9]. The tool path must be generated using the shape and the boundary of the repair volume to reduce the machining time. We will show that the CSFC method generates the tool path which follows exactly the boundary of the region being restored, produces a shorter tool path and reduces the time. Our basic curvilinear grid adapted to the shape of the blade is shown in Fig. 15(a). Unfortunately, the machine angular limits do not allow us to do that using the flat end cutter. Therefore, we used a ball-end tool of the radius 3 mm. The SFC tool path constructed for Tmax = 0.05 mm is shown in Fig. 15(b). The corresponding grid and the SFC tool path in the parametric region are shown in Fig. 15(c) and (d). The length of the curvilinear tool path is 6724.32 mm in the direction across the tip of the blade and 3991.86 mm in the direction along the tip. The reason for such a remarkable

difference is that the small machining strip on the tip of the blade effects distance between all the curves across the tip of the blade. However, few curves must be adapted when the tool path goes along the tool tip. The length of the isoparametric tool path across the tip is slightly larger then the curvilinear tool path whereas an isoparametric tool path is not applicable across the tip. Finally the length of the CSFC is 3286.41 mm which means a 17% decrease with regard to the best curvilinear grid in terms of the tool path length. The CSFC also gives a 7.5% decrease in terms of machining time for the linear feed rate 100 mm/min. The virtual cutting using the proposed CSFC tool path is shown in Fig. 16 whereas real machining (wood) is shown in Fig. 17. For demonstration purposes the size of scallops has been chosen so that the CSFC is visible on the blade surface. 6. Conclusions Numerically-generated adaptive-curvilinear grid is introduced to replace the rectangular grid used for construction of the space filling tool path for 5-axis machining. With this modification the SFC can be constructed for surfaces with complex irregular boundaries, cuts off, pockets, islands, etc. Besides, the adaptive grid allows us to efficiently treat complex spatial

W. Anotaipaiboon, S.S. Makhanov / Computer-Aided Design 40 (2008) 350–367

363

Fig. 15. Curvilinear grid adapted to part of a surface of the blade in (a) on the surface (c) in the u–v domain, (b) CSFC on the surface, (d) CSFC in the u–v domain.

variability of the constraints. The combination of the two techniques is superior with regard to the case when the two methods are applied independently. Acknowledgements We acknowledge a sponsorship of the Thailand Research Fund, grant BRG50800012. We also thank the National Electronics and Computer Technology Center of Thailand, and the National Science and Technology Development Agency of Thailand for their support. Finally, we wish to thank Ir. E. Bohez for fruitful discussions and excellent cooperation and the anonymous referees of the paper for their insightful remarks and attention to details. Appendix A. Machining strip evaluation Fig. 16. The virtual blade restoration. Simulation in the unigraphics.

Given the maximum allowable scallop height h max , the distance between the tool tracks is found by computing the

364

W. Anotaipaiboon, S.S. Makhanov / Computer-Aided Design 40 (2008) 350–367

Appendix B. The tool orientation The effective cutter radius re is given by 2 3/2 , where a = sin λ cos ω, b = re = ra 2 ( a1+b 2 +b2 ) tan λ sin ω, [17]. To optimize the machining strip width, λ and ω are usually set so that re is the best match to the radius of curvature at the CC point. For convex or planar surfaces, the tool inclination angle λ is set to a small default angle or zero and the tilt angle ω is set to zero as well. If the surface is nonconvex, a nonzero tool inclination angle λ is needed to avoid gouging. Consider a flat-end cutter shown in Fig. B.1. Gouging occurs whenever a point on the circle touches or goes inside the surface. Let G be a gouging point (Fig. B.1(a)). The line connecting the two points, Ol and G, forms a chord on the circle. Denote the angle between this line and O1 Oc by φ (see Fig. B.1(a) and (b)). Let λφ be the tool inclination angle that corresponds to a specific φ. The minimum tool inclination angle to avoid gouging is then λmin = Fig. 17. Blade restoration using CSFC, real machining.

max

−π/2≤φ≤π/2

λφ .

It is not hard to demonstrate that for a nonconvex surface λmin = sin−1 (r kmax ),

machining strip width. Let us introduce a local coordinate system (Ol , xl , yl , z 1 ) at the CC (cutter contact) point Ol shown in Fig. A.1, where xl denotes the normalized projection of the tool cutting direction onto the tangent plane, zl denotes the surface normal vector, and yl = zl × xl . The tool is rotated by an inclination angle λ about the yl axis, then by a tilt angle ω about the zl axis. The projected bottom edge of a flat-end cutter with radius r onto the (yl , zl )-plane becomes an ellipse called the effective cutting shape. In order to evaluate the machining strip, the surface cross-section perpendicular to the tool cutting direction xl is approximated by a circular arc, for which radius R y is equal to the radius of the normal curvature1 of the surface in the yl direction as shown in Fig. A.1. Suppose that h = h max . The maximum machined surface error is represented by a virtual circular arc with radius R y −h as shown in Fig. A.2. The machining strip width is then obtained by finding intersections of the effective cutting shape with the virtual arc. Let P be an arbitrary point on the cutter bottom edge (see Fig. B.1). We consider an angle θ required to turn the yc axis around the z c axis in such a way that the negative yc -axis passes through P. Furthermore, angles corresponding to the left and the right intersections Pl and Pr are denoted by θl and θr respectively. It is not hard to demonstrate that the left and the right machining strip wl and wr are then given by wµ = r | cos ω cos θµ − cos λ sin ω(1 − sin θµ )|, where µ = l or µ = r . The entire machining strip width is then w = wl + wr .

1 The radius of normal curvature is positive for concave surface and negative for convex surface.

where kmax is the maximum surface curvature at the CC point. Clearly, for a convex surface an inclination is not required, so λmin = 0, however, from technological viewpoints a small inclination angle is still recommended. Furthermore, for ∀λ ≥ λmin gouging will be avoided as well. It also can be shown [17] that for the flat end (cylindrical) cutter the orientation maximizes the machining strip. If gouging cannot be eliminated by inclining the tool alone or the inclination angle λ requires rotations which exceed the limit of the machine, the tilt angle ω can be optimized or a smaller tool size must be used. Appendix C. Correction of the CSFC tool paths The CSFC requires two modifications. First, the tool path should be modified to eliminate undercut areas. Second, the tool orientation needs to be carefully set when the tool is changing the direction. The tool path adjustment is needed since the SFC tool path contains turns that cause the tool to miss some areas of the surface when the tool is changing the direction. At each turn, the machining strips of the two adjacent tool paths travelling in different directions might not overlap leaving an undercut or the overlap is insufficient and produces a large scallop height (see Fig. C.1). To eliminate the error, the tool path at the turn is modified so that the machining strip overlaps with the machining strip of the adjacent tool path. Fig. C.2 displays the two types of tool path alterations at the corner (left and right), and at the U turns. At these turns, the tool path is extended so that it overlaps with the machining strip of the adjacent tool path. The minimal extension length is L ext = L − w,

W. Anotaipaiboon, S.S. Makhanov / Computer-Aided Design 40 (2008) 350–367

Fig. A.1. Geometrical analysis of the cutter.

Fig. A.2. Machining strip width estimation.

Fig. B.1. Tool gouging.

365

366

W. Anotaipaiboon, S.S. Makhanov / Computer-Aided Design 40 (2008) 350–367

Fig. C.1. Machining strips (dashed lines) on adjacent tool paths generated by using space-filling curve.

Fig. C.3. Trajectories of the cutter’s effective cutting edge (projected onto the x–y plane) (a) before and (b) after the tool path correction.

Fig. C.3 shows the trajectories of the effective cutting edge of the tool projected onto the x–y plane before and after applying the tool path correction. References

Fig. C.2. Tool path trajectory alteration at (a) corner turn (left or right) and (b) U turn.

where L and w denote the tool path interval and the side machining strip width (left or right), respectively, as shown in Fig. C.1 The second modification is applied to the tool orientation. The tool orientation is usually set by inclining the tool by λ in the tool cutting direction.2 At sharp turns, the tool orientation changes abruptly creating large kinematics error. This kinematics error could not be reduced by merely inserting more points as usually done for tool path segmentation [4,22]. Additionally, the tool orientation of the newly inserted CC point needs to be adjusted by interpolating the tool orientations at the two adjacent CC points. Consider the sharp turn o–p–q shown in Fig. C.2. To make a correct turn, the feed direction at the turning point p is first aligned with the feed direction at the previous point o. To reduce the kinematics error when going from point p to point q, a new point p0 is inserted and the feed direction is set to fp + fq , fp0 = fp + fq where fp is the feed direction at point p. 2 The tool cutting direction is the direction from the current CC point to the next CC point.

[1] Anotaipaiboon W, Makhanov SS. Tool path generation for five-axis NC machining using space-filling curves. In: Proceeding of the third asian conference on industrial automation and robotics, vol. 1; 2003. p 23–8. [2] Anotaipaiboon W, Makhanov SS. Tool path generation for five-axis NC machining using adaptive space-filling curves. Int J Prod Res 2005; 1643–65. [3] Anotaipaiboon W, Makhanov SS, Bohez EJ. Optimal setup for 5-axis machining. Int J Mach Tools Manu 2006;964–77. [4] Bohez E, Makhanov SS, Sonthipermpoon K. Adaptive nonlinear tool path optimization for five-axis machining. Int J Prod Res 2000;4329–43. [5] Brackbill JU, Saltzman JS. Adaptive zoning for singular problems in two dimensions. J Comput Phys 1982;342–68. [6] Bieterman MB, Sandstrom DR. A curvilinear tool-path method for pocket machining. J Mater Process Technol 2003;125(4):709–15. [7] Charakhch’yan AA, Ivanenko SA. A variational form of the Winslow grid generator. J Comput Phys 1997;136(2):385–98. [8] Cox JJ, Takezaki T, Ferguson HRP, Kohkonen KE, Mulkay EL. Spacefilling curves in tool-path applications. Comput Aided Design 1994;26: 215–24. [9] Gao J, Chen X, Zheng D, Yilmaz O, Gindy N. Adaptive restoration of complex geometry parts through reverse engineering application. Adv Eng Softw 2006;37(9):592–600. [10] Griffiths JG. Tool path based on Hilbert’s curve. Comput-Aided Design 1994;26:839–44. [11] Hopcroft JE, Ullman JD. Introduction to automata theory, languages, and computation. New York: Addison-Wesley; 1979. [12] Ivanenko SA. Generation of non-degenerate meshes. USSR Comput Math Math Phys 1988;28:141–6. [13] Ivanenko SA. Harmonic mappings. In: Thompson JF, Soni BK, Weatherill NP, editors. Handbook of grid generation, CRC Press LLC; 1999. pp. 8 (1–43). [14] Lam WM, Shapiro JH. A class of fast algorithms for the Peano–Hilbert space-filling curve. Proceedings ICIP-94, vol. 1. IEEE Computer Society; 1994. p. 638–41. [15] Lauwers B, Dejonghe P, Kruth JP. Optimal and collision free tool posture in five-axis machining through the tight integration of tool path generation and machining simulation. Comput Aided Design 2003;35:421–32. [16] Lee Y-S. Admissible tool orientation control of gouging avoidance for 5axis complex surface machining. Comput Aided Design 1997;29:507–21. [17] Lo C-C. Efficient cutter-path planning for five-axis surface for machining with a flat-end cutter. Comput Aided Design 1999;31:557–66. [18] Rao A, Sarma R. On local gouging in five-axis sculptured surface machining using flat end tools. Comput Aided Design 2000;32:409–20. [19] Robert FS. Applied combinatorics. Englewood CLiffs (NJ): Prentice-Hall; 1984. [20] Dafner R, Cohen-Or D, Matias Y. Context-based space filling curves. Comput Graph Forum 2000;19:C209–C217. [21] Makhanov SS. An application of the grid generation techniques to optimize a tool-path of industrial milling robots. J Comput Math Math Phys 1999;39:1589–600.

W. Anotaipaiboon, S.S. Makhanov / Computer-Aided Design 40 (2008) 350–367 [22] Makhanov SS, Batanov D, Bohez E, Sonthipaumpoon K, Anotaipaiboon W, Tabucanon M. On the tool-path optimization of a milling robot. Comput Ind Eng 2002;43:455–72. [23] Makhanov SS, Ivanenko SA. Grid generation as applied to optimize cutting operations of a five-axis milling machine. Appl Numer Math 2003; 46:353–77.

367

[24] Thomson JF, Soni B, Weatherill N. Handbook of grid generation. CRC Press; 1999. [25] Winslow AM. Numerical solution of the quasilinear Poisson equation in a nonuniform triangle mesh. J Comput Phys 1966;1(2):149–72. [26] Yoon J-H. Tool tip gouging avoidance and optimal tool positioning for 5-axis sculptured surface machining. Int J Prod Res 2003;41:2125–42.