Improving the dynamics of five-axis machining through optimization of workpiece setup and tool orientations

Computer-Aided Design 43 (2011) 1693–1706

Contents lists available at SciVerse ScienceDirect

Computer-Aided Design journal homepage: www.elsevier.com/locate/cad

Improving the dynamics of five-axis machining through optimization of workpiece setup and tool orientations Pengcheng Hu, Kai Tang ∗ Department of Mechanical Engineering, Hong Kong University of Science and Technology, Hong Kong

article

info

Article history: Received 21 April 2011 Accepted 8 September 2011 Keywords: Five-axis NC machining Inverse kinematics Angular acceleration Workpiece setup Tool orientation

abstract Existing works in optimization of five-axis machining mainly focus on the machining efficiency and precision, while the dynamic performance of the machine tools has not been fully addressed, especially in high-speed machining, in which the rotary actuators have limited dynamic ability. In this paper, a study is reported on how to generate a tool path so that the maximal angular accelerations of the rotary axes of the five-axis machine can be reduced. Two independent methods are proposed for this task: (1) by optimizing the setup of the workpiece on the machine’s table, and (2) by finding better tilt and yaw angles for the tool orientations. In this paper, the setup parameters of the workpiece are incorporated into the inverse kinematic equations, and angular acceleration functions are established according to the numerical solutions of those equations. While varying the tool orientations unquestionably would affect the surface quality of the machining, we introduce the so called Domain of Geometric Constraints that will restrict the allowable tilt and yaw angle of the tool at the cutter contact points on the part surface, so to ensure the satisfaction of the requirement of both local-gouging-free and cusp-height. For the first method – finding the optimal workpiece setup – a heuristic-based approach, i.e., the Genetic Algorithm (GA), is adopted, whereas for the second method – the constrained optimization of tool orientations – we present an elaborate algorithm based on the results from the analysis conducted by the authors. At the end of the paper, computer simulation experiments are reported that demonstrate the effectiveness of our proposed methods and algorithms. © 2011 Elsevier Ltd. All rights reserved.

1. Introduction Five-axis Numerical Controlled (NC) machining has been widely used in industry for processing sculptured surfaces, such as ship impellers, turbine blades and propellers, high-precision lenses, etc. Five-axis NC machining is much more complex as compared to traditional three-axis machining because of the two additional rotary axes. Although the research in five-axis NC machining has been extensive, with its focus spreading over a variety of topics, such as the tool trajectory generation [1–5], the tool orientation optimization [6–10], the kinematics error reduction [11,12], the tool path interpolation [13], the singular configuration avoidance [14,15] and so on, there are still many pending problems, and one of them is the dynamic performance of the machine tool’s rotary axes. Specifically, a tool path that is deemed good or even optimal according to some existing criteria might demand too high angular acceleration of the rotary axes and thus exert too severe – often infeasible – dynamic loading on the

∗

Corresponding author. Tel.: +852 2358 8656. E-mail address: [email protected] (K. Tang).

0010-4485/$ – see front matter © 2011 Elsevier Ltd. All rights reserved. doi:10.1016/j.cad.2011.09.005

machine tool. This is especially of great concern in the case of highspeed machining, where actuation of the rotary axes is always the bottleneck for improving the machining efficiency. Our objective is thus to conduct an in-depth study of this dynamic loading problem and provide practical and efficient solutions to it. In the following, we first review some background and past works on the general subject of tool path generation in five-axis NC machining, starting with the inverse kinematics of five-axis machining. 1.1. Inverse kinematics of the five-axis machining The existing research on inverse kinematics related to fiveaxis machining has mainly focused on the kinematic model, the kinematic errors and singularity configurations. Bohez et al. [16] gave a systematic introduction of different types of five-axis milling machines, and analyzed their respective disadvantages and advantages. In a further work [17], they summarized the related errors and built an error model in five-axis machining. Lee and She [18] established an inverse kinematic model of the three most popular types of five-axis milling machines, but no analysis of kinematic errors was given. Sorby [14] extended the

1694

P. Hu, K. Tang / Computer-Aided Design 43 (2011) 1693–1706

work of [18] by establishing a mathematical model of a fiveaxis machine tool with non-orthogonal rotary axes, and proposed a solution dealing with the singularities in inverse kinematics. However, the robustness of his solution to singularities is achieved at the cost of requiring very small tool orientation changes, which inevitably is at risk of causing larger accumulated numerical errors. Affouard et al. [15] introduced an interesting concept, called singular cone, to represent a singular point at which the tool path is deformed according to the singular cone and interpolated with B-spline curves to avoid the cone area. In the research of reducing inverse kinematics errors, Munlin et al. [12] minimized the kinematics errors near a stationary point of the machined surface by formulating the problem as a shortest path searching problem with the constraint region of a combination of the angles of A and B axis, and the optimization process is carried out by a greedy discrete method—the Dijkstra’s algorithm. Later, Anotaipaiboon et al. [11] minimized the kinematics errors by choosing an appropriate workpiece setup in five-axis machining, which includes no more than 6 parameters, i.e., the position and orientation of the workpiece with respect to the machine coordinate system and the machine’s configuration, and the optimization process is carried out by the least squares method. 1.2. Optimization of five-axis machining Over the last 15 years or so, there has been a surge of research interest in various kinds of optimization in five-axis machining. Those optimizations, however, mainly aim at improving machining efficiency and/or machining precision from the point of view of part coordinate system. In the planning of the tool trajectory (also referred to as the cutter contact curves), Suresh and Yang [19] proposed the concept of iso-scallop that aims at maintaining equal scallop height in a tool path. Lo [20] extended the idea further by adaptively selecting the cutter-inclination for a flat-end cutter, while at the same time trying to minimize the total length of the tool path. Agrawal et al. [21] minimized the length of the tool path by choosing proper primary path direction using the Genetic Algorithm (GA). Wang and Tang [22] presented a tool path generation algorithm based on the so called iso-conic concept in which a special direction on the cone is found in the part surface to generate the cutter contact curves, reducing the angular velocity and angular acceleration of the A-axis into zero. Unfortunately, their algorithm only applies to the specific five-axis machine with two rotary axes on the spindle. Whereas in the planning of tool orientations in a tool path, most existing methods opt to the fixed-orientation practice, i.e., the tool orientation is fixed relative to the surface normal vector at the Cutter Contact (CC) point [23]. This inflexibility suffers from the serious drawback that, for a part surface with high curvature, the tool orientation often has to undergo drastic change that in turn causes large angular velocity and acceleration on the rotary axes of the machine. Morishige et al. [24] relied on the powerful algebraic tool of C -space to determine the tool orientation, wherein the thus obtained tool orientations are further smoothed with the collisionfree condition upheld; however, some important local geometric properties, such as local gouging, is not taken into consideration in their algorithm. Jun et al. [8] improved this by adding the local-gouging-free constraint into the C -space; in addition, to address the issue of dramatic change in tool orientation, they proposed to smooth the tool orientation in the C -space. In their recent work, Wang and Tang [9] went one step further by also taking into account the limit on the maximum angular velocity of the tool and presented an elaborate tool path generation algorithm that entails a systematic framework knitting together all the necessary considerations such as collision avoidance, local gouging, and the angular velocity of the tool. Nevertheless, all the

works in [24,8,9] are done with respect to the Part Coordinate System (PCS), independent of the specific machines, and therefore are susceptible to causing heavy dynamic loading on the machine’s rotary axes—it is well known that a tool path smooth in tool orientation (but only measured in PCS) may deceivingly impose severe angular acceleration on the A, B, or C axis of a specific fiveaxis machine. Aside from the abovementioned, there are also many attempts at optimizing tool orientations, with different objectives and different techniques. Lavernhe et al. [6] optimized the tool orientation of every CC point, aiming at finding the best kinematic parameters, thus reducing the difference between the programmed and the actual feed rate in high-speed machining. Ye and Xiong [25] improved the machining accuracy by optimizing the geometric machining parameters such as the tool orientation, the cutter shape, and the setup of workpiece. Ho et al. [7] used the quaternion interpolation scheme to smooth the tool orientation, aiming at reducing the cutting error, though only in the PCS. Fleisig and Spence [13] used continuous splines to approximate tool position and orientation, which yields near constant feed speed and reduces angular acceleration; unfortunately again, the interpolation is carried out in the PCS. Chiou and Lee [26] searched for better tool orientations based on the swept envelope approach, considering both the local gouging and global collision. Kim and Sarma [27] modeled a tool path as the streamlines of a vector field, which is constructed by selecting the best directions of CC points subject to the constraints of the machine’s rotary motor speed limits and the cusp-height limit, and the best directions for performance envelope are found by a greedy approximation method; however, this very complicated system is not implemented and no specific machine-dependent solution is offered. Chiou [28] introduced the concepts of floor, wall and ceiling, for the determination of optimal local-gouge-free tool orientations, albeit in PCS. Recently, Castagnetti et al. [29] introduced an interesting concept of Domain of Admissible Orientation (DAO) which delineates the allowable tool orientations at each CC point. The DAO is originally defined in PCS as a rectangle, which then is transformed into an irregular area in the specific Machine Coordinate System (MCS) and is subsequently approximated by a quadrilateral. Every point in the quadrilateral represents a valid tool orientation for the associated CC point. Optimizations can then be performed directly in these quadrilaterals to find the best tool orientations for the desired objective. Debout et al. [30] utilized the concept of DAO in the smoothing of tool orientations for a particular manufacturing process – the automated fiber placement – which though cannot be extended to general five-axis machining. In summary, existing five-axis tool path generation or optimization algorithms are not able to address the following fundamental question: given a free-form part surface, a specific type of five-axis machine, and a required limit of cusp-height, how to generate a tool path that will cut the surface with the cusp-height requirement upheld while at the same time minimize the dynamic loading on the rotary axes of the machine? While the dynamics of a fiveaxis machine could be measured in various different ways, in this paper we focus on the angular accelerations of the rotary axes (the A, B, or C axis) of the machine, i.e., how to minimize the maximal magnitudes of angular accelerations, individual or combined. 1.3. Contribution of the paper In this paper, we present two optimization methods for the prescribed problem of minimizing the maximal (magnitude of the) angular acceleration of either an individual rotary axis, e.g., the A-axis, or the quadratic sum of two axes, e.g., the A and B axis. In the first method, we take as input an already computed tool

P. Hu, K. Tang / Computer-Aided Design 43 (2011) 1693–1706

1695

Fig. 1. Table-tilt type five-axis machine. Fig. 2. Procedure of inverse kinematics.

path that meets the cusp-height requirement – which though is computed in PCS as all the current tool path generation algorithms do – and find, among all the possible setups of the workpiece, a best workpiece setup with respect to the given five-axis machine (i.e., the rigid body transformation from PCS to MCS) that will minimize the maximum of the angular acceleration on the specific rotary axis. In the second method, with respect to a fixed workpiece setup, we compute a tool path that will guarantee the satisfaction of the cusp-height requirement while at the same time strive to minimize the maximal angular acceleration. These two methods can also be combined together in an iterative manner. The paper is organized as follows. First of all, the inverse kinematics of five-axis machines is analyzed, and parameters for the setup of the workpiece are introduced and incorporated into the inverse kinematics. Next, we present a simple iso-parametric tool path generation algorithm with which an initial tool path (that is, the CC curves and their associated tilt and yaw angles) will be computed. The Domain of Geometric Constraint (DGC) is then defined. Details of the algorithms of the two minimization methods are presented next. The experiment results of two representative surfaces are then presented with discussions. Finally, we conclude the paper with some comments on the future research in this subject. 2. Modeling of the setup of the workpiece in inverse kinematics In this section, we first introduce the parameters that define the setup of the workpiece in terms of the inverse kinematics of the five-axis machine, and then establish the relationship between the angular accelerations of the rotary axes and these parameters. In this paper we focus exclusively on the table-tilt type five-axis machine, but the general idea presented can be easily generalized to other types of five-axis machine. Fig. 1 shows a table-tilt type machine. 2.1. Inverse kinematics of the five-axis machine Although inverse kinematics has been well studied, the analysis given here is somewhat different from the traditional one for the following reasons. (1) The inverse kinematics we consider here contains the procedure directly from the placement of the tool to the rotary axes of the machine—that means, if the position and orientation of the tool is given, we can explicitly obtain its corresponding solutions in A and B. (2) The setup of the workpiece is incorporated into the inverse kinematics, and the parameters related to the setup will be included in the kinematic equations.

(3) Only the solutions for the two rotary angles of A and B are needed. We will exclusively assume that our part surface is a parametric surface S (u, v) defined in the PCS, whereas the situation if the part surface is given in the form of a mesh surface will be discussed in the end of the paper. The PCS is installed relative to the Table Coordinate System (TCS) of the machine which then is further defined with respect to the MCS. In Fig. 2, the tilt angle α and the yaw angle β are defined in Local Coordinate System (LCS) at each CC point. For a table-tilt type five-axis machine, the procedure for the inverse kinematics is shown in Fig. 2. In the process of solving angle ϕ A and ϕ B , which are respectively the angles of A and B axis, the equations are constructed according to the fact that the tool axis of the cutter in MCS is parallel with zm . Thus, the relationship is obtained by transforming the vector T in LCS into zm in MCS. 1. The tool orientation (T ) in LCS: T = T L = Rot (n, −β) ∗ Rot (k , −α) ∗ n

(1)

where: n is the surface normal vector measured in LCS, and in homogeneous format, n = [0, 0, 1, 0], matrices Rot (n, −β) and Rot (k , −α) have the same expressions as the standard matrices Rot (z , −β) and Rot (x, −α), respectively. 2. The tool orientation (TP ) in PCS: TP = TL→P ∗ TL TL→P = Trans (xP , yP , zP ) ∗

(2)

[

T

k 0

T

f 0

T

n 0

0 1

]

where: xP , yP , zP are the coordinates of the CC point in PCS. f

T

 =

∂S ∂v

 T  ∂S    ,  ∂v 

T

n =



∂S ∂S × ∂ u ∂v

T   ∂S ∂ S    ∂ u × ∂v  ,

k T = (f × n)T . 3. The tool orientation (TT ) in TCS TT = TP →T ∗ TP

(3)

TP →T represents the setup of the workpiece on the machine table, i.e., the transformation from the PCS to TCS, and its mathematical expression is: TP →T = Trans (x0 , y0 , z0 ) ∗ Rot (z , r2 ) ∗ Rot (x, r1 ) ∗ Rot (z , r3 ) where, x0 , y0 , z0 are the coordinates of the origin of PCS in TCS, r1 , r2 , r3 are the three Euler angles from PCS to TCS, and Rot (z , r2 ) , Rot (x, r1 ), Rot (z , r3 ) are standard rotation matrices.

1696

P. Hu, K. Tang / Computer-Aided Design 43 (2011) 1693–1706

Fig. 4. CC points on one CC curve.

2.2. Numerical solution of aA and aB We first need to define the time parameter t. For the convenience of calculating the angular acceleration, the tool is assumed to maintain a constant feed rate f on the CC curve, which though is different from the feed rate F defined in most of the CNC [31]; in our setting, the CC curve arc length l is proportional to time t. In other words, we use l to represent t when calculating aA and aB . The exact analytical expressions could be obtained from Eqs.

Fig. 3. Coordinate transformation from PCS to TCS.

The transformation is shown is Fig. 3.

d2 ϕ A

4. The tool orientation (TM ) in MCS TM = TT →M ∗ TT

(4)

where, TT →M = Trans (x1 , y1 , z1 ) ∗ Rot y, ϕ B ∗ Rot x, ϕ A . ϕ A and ϕ B are respectively the angles of A and B axis; x1 , y1 , and z1 are the coordinates of the crossing point of A and B measured in MCS. As on a table-tilt type machine, the orientation of the tool is constant in MCS, for the machine chosen in this paper, shown in Fig. 1, TM = [0, 0, 1, 0], in homogeneous coordinates. From Eqs. (1)–(4), the final equation for solving ϕ A and ϕ B is:









TM = TT →M ∗ TP →T ∗ TL→P ∗ Rot (n, −β) ∗ Rot (k, −α) ∗ T .

(5)

Assuming: TP →T ∗ TL→P ∗ Rot (n, −β) ∗ Rot (k, −α) ∗ T

= [m1 , m2 , m3 , 0]

(6)

m1 , m2 and m3 are functions of parameters u, v, α, β, r1 , r2 and r3 . Substitute Eqs. (1)–(4) and Eq. (6) into Eq. (5), we obtain two inverse solutions:

ϕ A = arctan(m2 /m3 )

(7)

ϕ = −arctan(m1 /(m2 sin(ϕ ) + m3 cos(ϕ ))). B

A

A

(8)

One important fact of the above two expressions is that ϕ A and ϕ B do not contain the parameters (x0 , y0 , z0 ) and (x1 , y1 , z1 ), which implies that in the inverse kinematics, the position of the workpiece on the machine table, which is represented by (x0 , y0 , z0 ), has no influence on ϕ A and ϕ B , and the same for (x1 , y1 , z1 ) and (xP , yP , zP ). This is due to the fact that, in the derivation leading to ϕ A and ϕ B , all the parameters related to translation are vanished because they are multiplied by the fourth element (equal to 0) of vector T . It is necessary to state that the solutions in Eqs. (7) and (8) in general have two branches, since trigonometric functions are involved. In the current study, only one is considered. Moreover, singularities in the solutions are not considered either. This nonuniqueness of the solutions, in particular when coupled with singularities, is by itself an interesting research topic. Next, we will study the influence of the three Euler angles r1 , r2 and r3 on aA and aB .

d2 ϕ B

(7) and (8) by taking aA = f 2 dl2 and aB = f 2 dl2 . The resultant expressions are however quite complicated and hence time consuming to calculate. In our implementation, numerical methods are adopted to calculate aA and aB instead. As illustrated in Fig. 4, at each CC point, there is one associated pair of ϕ A and ϕ B . The angular accelerations at CCi can be obtained by numerical differentiation:

 aAi

aBi

= 2f

2

= 2f

2

ϕiA+1 − ϕiA



li

  ϕiB+1 − ϕiB li

−

 A   ϕi − ϕiA−1 li−1

 −

ϕiB − ϕiB−1 li−1

 

(li + li−1 )

(9)

(li + li−1 )

(10)

where, i = 2, 3, 4, . . . . In the numerical calculation of aA and aB at CCi , the angles of the point itself and its two adjacent points are used. ϕiA , ϕiB and li are functions of parameters u, v, α, β, r1 , r2 and r3 , and from Eqs. (9) and (10), we know that aA and aB are also functions of u, v, α, β, r1 , r2 and r3 (assuming f to be a constant). It needs to be pointed out that, theoretically, the curve length l in Eqs. (9) and (10) should be along the trajectory of tool tip point which is different than the CC point. However, due to the proximity between CCi and CCi+1 , the difference between ‖CCi − CCi+1 ‖ and its counter-part on the trajectory of the tip point should be negligible. 3. Initial tool path generation Our tool path generation algorithm can be outlined into two phases. In the first phase, an initial tool path is generated in PCS. This initial tool path is computed based on the popular idea of maintaining equal maximum scallop height along the tool path, which is widely used in machining of complicated sculptured surfaces because of its simplicity and computational efficiency. In the second phase, after the initial tool path is available, with the CC curves kept unchanged, we then modify the tilt and yaw angles at the sample CC points of the CC curves to minimize the maximal angular acceleration of the A and/or B axis, while in all the time satisfying the cusp-height constraint. The details of Phase 1 are given next while that of Phase 2 will be given in the next section. On the initial tool path, the tool orientation is constant at every CC point in LCS. The generation scheme itself is very simple: starting from u1 = 0 (assuming the parameter domain of S (u, v) is [0, 1] × [0, 1]), the current CC curve is the iso-parametric curve S (ui , v) which is further sampled into discrete CC points Ci =

P. Hu, K. Tang / Computer-Aided Design 43 (2011) 1693–1706

1697

{S (ui , v1 = 0) , S (ui , v2 ) , S (ui , v3 ) , . . . , S (ui , vn = 1)} for some n; a constant tilt angle α and yaw angle β are assigned to the tool orientation at all the CC points in Ci ; a 1ui is then sought for the next CC curve Ci+1 at ui+1 = ui + 1ui , such that the cusp-height due to the two adjacent CC curves Ci and Ci+1 is under the given limit; this process then continues until the end condition u = 1 is reached. 3.1. Initial tool orientation At the beginning, the tilt angle α and yaw angle β are fixed in LCS and their initial value α0 and β0 are assigned based on the geometrical proprieties of the surface and the radius of the flat-end cutter. In our current setting, β0 is set to 0, which is a common choice for flat-end cutter. The initial value of the tilt angle α0 is decided by considering both the machining efficiency and the local-gougingfree requirement. To improve machining efficiency, the effective cutting radius r e should be as large as possible, since a larger effective cutting radius will generate a wider machining strip width [32]. However, on the other hand, due to the requirement of no-local-gouging, r e must be smaller than the local curvature of the surface in the direction perpendicular to the feed direction (which is the partial derivative ∂ S /∂v ) in concave regions of the surface, which is expressed as r ∗ as defined in [19]. Considering altogether the two conflicting conditions above, the minimum, min(r ∗ ), of all the r ∗ in the concave regions of the surface will be chosen as the (initial) effective cutting radius, as described below:

 r0∗ = min r ∗ (u, v) : (u, v) ∈ [0, 1] 

and

II (u, v) > 0

where, II is the second fundamental form of the part surface. The r0∗ is taken as the (initial) effective cutting radius of the whole surface. Then, the initial tilt angle α0 is obtained by solving Eq. (A.1) (see Appendix A) with β0 = 0. 3.2. Initial tool path generation In the equal maximum cusp-height tool path generation process on a parametric surface, we assume that u = constant and v varies in one CC curve on the sculptured surface S (u, v), and all the calculations are based on this assumption henceforth in this paper. The mathematics of the calculations is based on the works in [19,20]. Our tool path generation process is divided into two steps. Step 1. Forward-step calculation: decide the sample CC points on the ith (i = 1, 2, 3, . . .) CC curve Ci . In this step, point CCi,j (the jth CC point on the ith CC curve) is obtained with ui = constant and vi,j obtained according to the method given in [19] by iteration, starting from vi,1 = 0, until the entire ith CC curve is sampled. Step 2. Side-step calculation: calculate 1ui,j at each CCi,j point (j = 1, 2, 3, . . .) using the method in principle from [19]. In this step, after all the 1ui,j are calculated, their minimum – 1ui = min(1ui,j ), (j = 1, 2, 3, . . .) – is obtained and the next CC curve will be the iso-parametric curve S (ui+1 , v) with ui+1 = ui + 1ui . Step 1 and 2 are repeated until the final CC curve with u = 1 is reached. Fig. 5 illustrates our scheme of generating the initial isoparametric tool path. The CCi,j (j = 1, 2, 3, . . .) point on the ith CC curve is calculated by forward-step calculation along v direction; at each CCi,j , 1ui,j is obtained by the side-step calculation to meet the specified cusp height constraint.

Fig. 5. Initial iso-parametric tool path generation.

3.3. Construction of domain of geometrical constraint After the initial tool path is computed, we next construct the Domain of Geometrical Constraint (DGC) of the tilt angle α and the yaw angle β for every sample CC point on the path. The DGC is a region in the α –β plane at each CC point in which any (α, β ) is admissible under the geometrical constraint of the surface. More specifically, the DGC is constructed according to the requirement of local-gouging-free and the specified cusp-height limit. To define DGC, let rie,j denote the effective cutting radius of the cutter at the CC point CCi,j . Local gouging could happen only when CCi,j is locally concave, and the gouging-free condition can be characterized by the following equation: rie,j ≤ Ri,j

(11)

where

 Ri,j =

ri∗,j ,

+∞,

IIi,j > 0 IIi,j ≤ 0.

The ri∗,j is the radius of local curvature of the geodesic path at CCi,j . Thus, as long as Eq. (11) is satisfied, CCi,j will be free of local gouging. Next, we define the lower limit of rie,j , in accordance with the specified cusp height limit. When the CC curve Ci of the initial tool path is generated, the increment 1ui = min(1ui,j ) (j = 1, 2, 3, . . .) is conservative except at those points (in general only one) where 1ui,j = 1ui . We divide the CC points into two groups, the Cusp-Height Critical Point (CCP) and None-Cusp-Height Critical Point (NCCP): CCP = {CCi,j |1ui,j = 1ui } NCCP = CCi,j |1ui,j ̸= 1ui .





On the ith CC curve Ci , let hi,j denote the cusp-height at CCi,j ; in most cases it should be smaller than the given cusp-height limit h. That is, there are two cases: Case 1. CCi,j ∈ CCP. Case 2. CCi,j ∈ NCCP. In the first case, the CC point is a CCP and hence its cusp-height is h; however, barring any rarity, in general exactly only a single point will fall into Case 1, while the rest belong to the second case. Refer to Fig. 6. At CCi,j , rie0 ,j is the initial effective cutting radius calculated by α0 , β0 and Eq. (A.1) (Appendix A). The 1ui,j is originally decided by h according to the method given in [19], and u′i+1 = ui + 1ui,j . However, after the final adjustment of ui+1 in the side-step (Fig. 5), ui+1 = ui + 1ui < u′i+1 . As a result, if the effective cutting radius rie0 ,j is kept unchanged, then the real cusp height is h0i,j , which is clearly smaller than h. Therefore, for those NCCP CC points, their cusp-heights can become larger from their initial values, as long as not exceeding h. As shown in Fig. 6, the smaller the effective cutting radius, the larger the cusp height. In other words, we can decrease (by

1698

P. Hu, K. Tang / Computer-Aided Design 43 (2011) 1693–1706

the yaw angle is set as β0 = 0, which belongs to sub-region RIi,j , we hence choose RIi,j to be the domain for the DGC of CCi,j , as illustrated in Fig. 7(b). 4. Optimization algorithms In our current study, the optimization mainly focuses on three objectives: the dynamics of A-axis, B-axis, and the combination of the two. The three can be represented by the angular acceleration of A, B, and their quadratic sum, respectively. Our objective is to improve the dynamic performance, which means to reduce the maximum of the (magnitude of the) desired angular acceleration. Specifically, the following three types of minimizations are to be solved:



  A   min max ai,j i ,j





 

min max aBi,j 

Fig. 6. Cusp height determination.

i ,j

min

−  (aAi,j )2 + (aBi,j )2

(14)

(15) (16)

i ,j

Fig. 7. (a) DGC with two disconnected sub-regions; (b) the chosen DGC Domain.

changing the tilt and yaw angle) the effective cutting radius at CCi,j until the corresponding cusp-height reaches h. Let rie1 ,j denote the smallest effective cutting radius allowed at CCi,j , i.e., it enforces the equation h0i,j = h.

From equations in [19], we can get the value of rie1 ,j with specified cusp-height h. rie1 ,j is the lower limit of the effective cutting radius at CCi,j point, in accordance with the requirement of the specified cusp-height limit, i.e., rie,j ≥ r i,j . e1

(12)

Combining both Eqs. (11) and (12), we get: e

rie1 ,j ≤ r i,j ≤ Ri,j .

(13)

From Eq. (A.1) (Appendix A), we know that rie,j is a function of α and β . Therefore, the region of DGCi,j should be defined according to Eq. (A.1) (Appendix A) and Eq. (13). As an illustration, in Fig. 7, the DGC at one CC point from a test example in our paper is shown. As seen in Fig. 7(a), due to the cyclic nature of sinusoidal functions involved in the definition of rie,j , each DGC is actually divided into two disconnected sub-regions separated at the special yaw angle β = π /2, denoted by RIi,j and RIIi,j respectively, and they are defined as follows:

RIi,j = {(α, β)|DGCi,j ∩ αϵ(0, π /2) ∩ βϵ(−π /2, π /2)}; RIIi,j = {(α, β)|DGCi,j ∩ αϵ(0, π /2) ∩ βϵ(π /2, 3π /2)}. Apparently, from the concern of both numerical stability and smooth continuity of the cutting motions, the point (α, β) must not cross between the two sub-regions. Since the initial value of

where aAi,j and aBi,j are defined and calculated by Eqs. (9) and (10) respectively. In the above three, the first two are of min–max type optimization, while the third is an overall minimization integrated over the entire surface. As already alluded in the beginning, the optimization comprises two independent modules: the Optimal Workpiece Setup module and the Optimal Tool Orientation Assignment module. In the first module, with respect to a given and fixed tool path (such as a one from Section 3.2) that though is computed in PCS independent of the MCS, we look for the best transformation from the PCS to MCS that will achieve the specified objective from the three equations (14)–(16). In the second module, with respect to a given and fixed transformation from the PCS to MCS, we iteratively modify the tilt–yaw angle pair (α, β) within the DGC of certain CC points, and continuously reduce the maximum (magnitude of the) angular acceleration of the A or B axis, or their quadratic sum, until a termination condition is reached. 4.1. The optimal workpiece setup module Owing to the nonlinearity of the inverse kinematics, per Eqs. (7) and (8), different transformations from PCS to MCS will yield different ϕ A and ϕ B , and hence their second order derivatives aA and aB as well. As this transformation is uniquely determined by the three Euler angles r1 , r2 , and r3 , with their domains r1 ∈ [0, π /2] and r2 , r3 ∈ [0, 2π ], we need to search within these domains to find a best solution, pertaining to the given optimization objective from Eqs. (14)–(16). In our current approach, the heuristic-based method GA – which is implemented by us using MATLAB – is utilized to realize this task. In the following, we present one example of our preliminary experiments to demonstrate the influence of the workpiece setup on those three objective functions Eqs. (14)–(16). For better elucidation, only the data of one particular CC curve on the part surface are shown here, while that of the experiments on an entire part surface will be presented in the next section. Refer to Fig. 8 and Table 1. In Fig. 8, four different cases are shown: (1) r = (r1 , r2 , r3 ) = r0 ; r0 is an arbitrary reference setup (i.e., the original one), and in our case r0 = (0, 0, 0). (2) r = (r1 , r2 , r3 ) = rA∗ ; rA∗ is the optimized result minimizing the maximal angular acceleration of A axis, i.e., corresponding to the objective function Eq. (14).

P. Hu, K. Tang / Computer-Aided Design 43 (2011) 1693–1706

1699

Fig. 8. (a) ϕ A ; (b) |aA |; (c) ϕ B ; (d) |aB |. Table 1 Results of three different optimizations. Unit: rad/s2 or (rad/s2 )2 .

max(|aA |) B max ∑ (|Aa 2|) ((a ) + (aB )2 )

r = r0

r = rA∗

r = rB∗

r = rA∗+B

15.8382 17.7772 2.3167 × 104

10.5391 20.5939 1.8775 × 104

20.8390 7.0233 1.6823 × 104

18.5879 14.7009 1.5349 × 104

(3) r = (r1 , r2 , r3 ) = rB∗ ; rB∗ is the optimized result minimizing the maximal angular acceleration of B axis, i.e., corresponding to the objective function Eq. (15). (4) r = (r1 , r2 , r3 ) = rA∗+B ; rA∗+B is the optimized result corresponding to the objective function Eq. (16). Fig. 8(a) and (c) show the angles ϕ A and ϕ B , respectively, each corresponding to the four cases above; whereas the corresponding (magnitudes of the) angular accelerations of the A and B axis, as well as their norm, are displayed in Fig. 8(b) and (d), respectively. Table 1 lists the relevant data from Fig. 8. From both the figure and table, it shows that: when individually optimized, the reduction in max(|aA |) or max(|aB |) is obvious, sometimes tremendous; however, the reduction on one axis may cause increase on the other. If the goal is the improvement in overall dynamics rather than on one particular axis, then perhaps the combined optimization (Eq. (16)) might suit better. Regardless, this simple example clearly demonstrates that changing the workpiece setup indeed will have direct effect on the angular accelerations of the rotary axes; the experiment also validates the effectiveness of GA in the optimization. 4.2. The optimal tool orientation assignment module With the transformation from the PCS to MCS (hence the workpiece setup) fixed, we now proceed to work on how to modify the tilt–yaw angle pair (α, β) inside the DGC of a CC point. Since the variation of (α, β) is restricted within the DGC, the cusp-height constraint will never be violated. On the other hand, plausibly, its effect on aA and aB is apparent. The crux here is to devise a computationally efficient numerical algorithm to fulfill the task.

Because of the extremely high algebraic complexity of aA and aB , let alone the highly nonlinear constraints imposed by the DCGs, it is virtually impossible to simultaneously move the (α, β) for more than one CC point. Instead, a plausible approach is to work on the CC points one at a time, in an iterative manner. Let CCi,j be the CC point whose angular acceleration (in magnitude) is currently the highest—either aA or aB , depending on the specific optimization objective. Take aA for example. Per Eq. (9), it is a function of three A-axis angles, namely the ϕiA,j itself plus its two neighbors ϕiA,j−1 and ϕiA,j+1 on the CC curve Ci . (Note that if higher order numerical differentiations are used in Eqs. (9) and (10), more neighboring CC points will be involved.) A simple analysis of Eq. (9) tells that, when ϕiA,j varies, the moving direction of aAi,j is opposite of that of aAi,j−1 and aAi,j+1 . A tempting local modification idea would naturally

arise: decrease or increase ϕiA,j (depending on the sign of aAi,j ), but   within its DGC, to reduce aAi,j  until it equals the new aAi,j−1 or aAi,j+1 , whichever occurs first. However, this strategy immediately runs into a quandary since, after this initial adjustment, none of the three angles can be modified without increasing their respective |aA |. We instead propose a reduction strategy based on the idea called flooring and ceiling (FC) operation. Suppose we change ϕiA,j (by modifying its α and β angle within its DGC) so that the new aAi,j is reduced by a small percentage, say 5%. Because of the change of ϕiA,j , however, aAi,j−1 is raised. If the new aAi,j−1 is still smaller than the new aAi,j , the current iteration ends. Otherwise, which we refer to be a ceiling situation, we modify angle ϕiA,j−2 (but not ϕiA,j−1 ) to pull down aAi,j−1 to be equal to the new aAi,j . As ϕiA,j−1 is kept unchanged,

so will aAi,j . On the other hand, aAi,j−2 has now been raised. If the

1700

P. Hu, K. Tang / Computer-Aided Design 43 (2011) 1693–1706

new aAi,j−2 is smaller than the new aAi,j−1 , no action will be taken and a new iteration starts. Otherwise, we modify angle ϕiA,j−3 (but

not ϕiA,j−2 ) to pull down aAi,j−2 to be equal to the new aAi,j−1 . This process continues until either no ceiling situation is met or the first CC point on the CC curve is reached. Note that this treatment is also applied to the other neighboring CC point ϕiA,j+1 in the exactly same manner, except this time the advancement moves forward along the CC curve. We next give the details of For better    the procedure. lucidity of description, assume the min maxi,j aAi,j  objective, and all the subscripts ‘‘A’’ are dropped in the variables in the procedure. Procedure Flooring_and_Ceiling (a CC curve, a CC point) /* The input are a CC curve and a CC point on the curve whose absolute angular acceleration is the largest. Without loss of generality, assume that the input CC point is CC i . */ BEGIN Step 1. M = |ai |, s = 1 if ai > 0 and –1 otherwise. Step 2. Reduce |ai | of the CC i point. In this step, according to Eq. (9), we tentatively try to reduce ai to λM (for some λ < 1 and very close to 1) by modifying the angle ϕi while keeping ϕi−1 and ϕi+1 unchanged. This step is divided into 3 sub-steps. 2.1. Solve the minimization problem. min(−sϕi (αi , βi )), subject to (αi , βi ) ∈ DGC i ; Let (αi∗ , βi∗ ) be the optimal point found. 2.2. Checkif it is feasible to reduce |ai | to λM.  IF sϕi αi∗ , βi∗ < sϕiλ Exit. /* current |ai | cannot be lowered further (to λM), exit with failure. */ ELSE Starting from (αi , βi ) = (αi0 , βi0 ), search for new (αi , βi ) that will satisfy ϕ i (αi , βi ) = ϕiλ . Where, ϕiλ is defined in Eq. (17). 2.3. Reassign (αi0 , βi0 ).

(αi0 , βi0 ) ← (αi , βi )

Step 3. Initialize index p. p ← 1, this index is used for the flooring and ceiling operation on the CC j (j < i) before CC i on the CC curve. Step  4. Check the ceiling condition at CC i−p (1 ≤ p ≤ i − 2). IF ai−p  ≤ λM, Go to Step 6. /* No pulling-down is needed; go to process CC j (j > i). */ /* Pulling-down is necessary at CC i−p . */  Step 5. Pull-down ai−p  to λM . Here we try to lower ai−p  to λM by modifying the angle ϕi−p−1 (instead of ϕ i−p ). This step is divided into 5 sub-steps. 5.1. Solve the minimization problem. min  (sϕi−p−1 (αi− p−1 , βi−p−1 )), subject to αi−p−1 , βi−p−1 ∈ DGC i−p−1 ; Let (αi∗−p−1 , βi∗−p−1 ) be the optimal point found.





5.2. Check if it’s feasible to reduce the ai−p  to λM. IF sϕi−p−1 (αi∗−p−1 , βi∗−p−1 ) > sϕiλ−p−1





Exit. /* Impossible to lower ai−p  to λM, exit with failure. */ ELSE   From αi−p−1 , βi−p−1 = (αi0−p−1 , βi0−p−1 ), search for new (αi−p−1 , βi−p−1 ) that will satisfy ϕi−p−1 (αi−p−1 , βi−p−1 ) = ϕi∗−p−1 .





Go to Step 6. 5.5. Move to the next CC point. p ← p + 1, and go to Step 4. Step 6. Initialize p. p ← 1, this index is used for the flooring and ceiling operation on the CC j (i < j) after CC i on the CC curve. Step  7. Check the ceiling condition at CC i+p (1 ≤ p ≤ n − i − 1). IF ai+p  ≤ λM, Exit. /* No pulling-down is needed; exit with success. */ /* Pulling-down is necessary at CC i+p . */  Step 8. Pull-down ai+p  to λM . The action here is identical to that in Step 5 except this time the advancement moves forward along the CC curve. The details are omitted. END In the procedure above we use a sign variable s to indicate the sign of ai : if s is ‘‘+1’’, then ai should be reduced, whereas a ‘‘−1’’ s means that ai should actually be increased (so to reduce its magnitude). The adoption of s helps put the treatment of both positive and negative angular acceleration in a single and more succinct form. In Step 2, firstly, at sub-step 2.1, we find the best value of ϕi for reducing |ai |. This is actually a nonlinear constrained optimization problem; to solve it, we first use the interior-point method to transform the original form ‘‘minx f (x), subject to g (x) ≤ 0’’∑ into its approximation form ‘‘minx,s fµ (x) = minx,s f (x) − µ i ln(si ) , (µ > 0), subject to g (x) + si = 0’’, by introducing more variables si . The approximated problem is then solved by the direct step search method and the conjugate gradient step search method. Next, at sub-step 2.2, we ensure that the reduction of|ai |  to λM is indeed realizable by checking the condition sϕi αi∗ , βi∗ < sϕiλ , where ϕiλ is given by (refer to Eq. (9)):

ϕiλ =

li li−1 li + li−1

(α

,β

  ) ← αi−p−1 , βi−p−1

5.4. Check if the start of the CC curve is reached IF p = i − 2

li

+

ϕi−1 li−1

−

sλ M 2f 2

] (li + li−1 ) .

(17)





pulling-down operation for ai−p  at Step 5 is similar to that of |ai | at Step 2 except for the calculation of the target ϕiλ−p−1 . Refer to Eq. (9), substituting the index ‘‘i’’ therein with ‘‘i − p’’, and  keeping the two angles ϕi−p and ϕi−p+1 constant, to force ai−p  to equal λM, the new ϕi−p−1 should satisfy the following equation:

ϕiλ−p−1 =

5.3. Reassign (αi0−p−1 , βi0−p−1 ). 0 i−p−1

ϕi+1

The angle ϕi (αi∗ , βi∗ ) represents the extreme case for ϕi to reduce |ai |. If even ϕi (αi∗ , βi∗ ) is not able to reduce |ai | to λM, that would mean that, albeit conservatively, no (αi , βi ) in DGCi can fulfill the reduction either.In this case, we simply exit the procedure (with failure). If sϕi αi∗ , βi∗ < sϕiλ , then there is a (αi , βi ) in DGCi that satisfies ϕi (αi , βi ) = ϕiλ . Our numerical algorithm of searching for (αi , βi ) is given in Appendix B. From Step 3 to Step 5, the CC points on the CC curve prior to CCi are processed,   one at a time. Basically, for the current CC point CCi−p , if ai−p  does not exceed λM, there will be no need to check the remaining CC points CCj (j < i − p) since their |a| is not affected or less than λM, and hence we  go  to process those CC points after CCi (Step 6–8). Otherwise, ai−p  is pulled down to λM by modifying    the angle  ϕi−p−1 . Because ϕi−p is kept unchanged, all the ai−p+1  , ai−p+2  , . . . , |ai | remain the same value λM. The

sλMli−p−1 (li−p + li−p−1 ) 2f 2

Where, ϕiλ−p−1 is defined in Eq. (18). 0 i−p−1

[

− 



  li−p−1 ϕi−p+1 − ϕi−p li−p

+ ϕi−p .

(18)

After ai−p  is pulled down, we then move backward to check the next CC point CCi−p−1 , CCi−p−2 , . . . , until either the first CC point

P. Hu, K. Tang / Computer-Aided Design 43 (2011) 1693–1706

1701

Fig. 9. (a) Results of FC (min–max|aA |); (b) results of FC (min–max|aB |).

Fig. 10. Test part surfaces: (a) surface S1 ; (b) surface S2 .

∗ ∗ Fig. 11. Setup of surface S1 in TCS when: (a) r = rA1 , and (b) r = rB1 .

is reached or a ‘‘No pulling-down’’ condition is detected (Step 4). The procedure now proceeds to Step 6–8 to process the CC points after CCi on the CC curve. The treatment is identical to that of Step 3–5 except this time the advancement moves forward along the CC curve. The details are omitted. In Fig. 9, an illustrative example of the floor and ceiling operations on a CC curve is given. In this particular example, all the aA (and aB ) on the entire CC curve eventually converges to a single value. This is however not the general case. As exemplified by our experiments given in the next section, in most cases the |aA | or |aB | on a CC curve after the optimization with FC operations exhibits a jagged form, exactly due to the constraints imposed by the DGCs. The choice of the value for λ deserves examination; λ should be some value between 0 and 1. Ideally, it should be as close to 1 as possible, so that |aA | (or |aB |) can be gradually ‘‘floored’’ without running into the danger of violating some DGC and having to stop. On the other hand, a too small 1 − λ will inevitably increase the computing time and the accumulated numerical errors. After a great deal of experimenting, in our current implementation, a λ = 0.91 is used. With the procedure Flooring_and_Ceiling available, our final algorithm for optimizing the tool orientations becomes straightforward. In a nutshell, given an initial tool path (Sections 3.1 and 3.2), we recursively do the following: Step 1. In the current tool path, find the CC point and the CC curve it belongs that has the largest maximal angular acceleration (depending on which of Eqs. (14)–(15) is used). Step 2. Call procedure Flooring_and_Ceiling with the CC point and curve found in Step 1. Step 3. If no improvement is made at Step 2, the entire optimization terminates. Otherwise, go to Step 1. Note that, since our main objective from Eqs. (14)–(15) is min–max over the entire part surface, the optimization stops if the maximum cannot be reduced. However, if our concern is some other type of dynamics, then the Flooring_and_Ceiling operation can be continually applied to other CC points, not necessarily the one with the maximum.

In our current system, the optimizations of Eqs. (14)–(15) are dealt with using the proposed Flooring_and_Ceiling operation. However, for Eq. (16), which not only involves two axes but also requires integration over the entire surface, a heuristic algorithm based on GA is employed. 5. Experiential results and discussions The optimization algorithms presented in Sections 3 and 4 are implemented using MATLab and a number of preliminary experiments are conducted. In this section, we report the results of some conducted experiments on two parametric surfaces S1 and S2 shown in Fig. 10. As S2 is symmetric in both in u and v at u = 0.5 and v = 0.5, its optimization is restricted to only one quadrant, i.e., the partial parametric domain u × v = [0.5, 1] × [0.5, 1]. A flat-end tool with the diameter of 15 units is used, whereas the size of the two surfaces each is about 100 units by 150 units—i.e., the curve lengths of the iso-curves at u = 0 and v = 0 are roughly about 100 units and 150 units respectively. In terms of the computing time, for both surfaces, the Optimal Workpiece Setup module and the optimizations of Eqs. (14) and (15) (including both the initial tool path generation and the subsequent tool orientation optimization) all run in one minute or so, while the GA algorithm for solving Eq. (16) takes several minutes to finish before achieving the data given in Tables 2–5. 5.1. Results of optimizations of workpiece setup To test the Optimal Workpiece Setup module (Section 4.1), a reference setup r = (r1 , r2 , r3 ) = r0 = (0, 0, 0) is used for comparison purpose. For surface S1 , the two optimized workpiece setups ∗ ∗ r = rA1 and r = rB1 are shown in Fig. 11, each corresponding to the min–max of |aA | and |aB | respectively. The distributions of |aA | and |aB | on surface S1 , before and after the optimization with ∗ ∗ r = r0 , r = rA1 and r = rB1 , are displayed in Fig. 12, while the maximum data are given in Table 2. Similarly, the Optimal

1702

P. Hu, K. Tang / Computer-Aided Design 43 (2011) 1693–1706

∗ ∗ Fig. 12. Angular accelerations on surface S1 : (a) original |aA | and |aB | with r = r0 ; (b) |aA | and |aB | with r = rA1 ; (c) |aA | and |aB | with r = rB1 .

∗ ∗ Fig. 13. Setup of surface S2 in TCS when: (a) r = rA2 , and (b) r = rB2 .

Workpiece Setup minimization is performed on surface S2 and the ∗ ∗ two found orientations in TCS, rA2 and rB2 , are shown in Fig. 13, A B with the distributions of |a | and |a | shown in Fig. 14, and the maximum data given in Table 3. The experimental results convincingly validate the idea of changing the workpiece setups to improve the dynamics of the rotary axes in five-axis machining. The improvement is often substantial; for instance, as seen in Table 2, the reductions in maximum |aA | and |aB | on surface S1 are respectively 29% and 43%. In the case of S2 , we purposely make it concave and have a near constant curvature in one parametric direction. Ideally, to such a surface, due to the constant feed rate assumption, it could be oriented in TCS in such a way that |aA | is near zero at all the CC points; our optimization algorithm clearly has captured this overall optimal result, as revealed by Fig. 14(b). 5.2. Results of optimizations of tool orientations Next, with the workpiece setup fixed at r = (r1 , r2 , r3 ) = r0 = (0, 0, 0), we apply the procedure Flooring_and_Ceiling to the two surfaces to find better tool orientations at the CC points on the entire surface, in terms of the objectives Eqs. (14)–(15) respectively. GA is adopted for the third optimization Eq. (16). Let α(S ) and β(S ) stand respectively for the tilt and yaw angle on surface S. Then, let (α∗A1 , β∗A1 ) and (α∗B1 , β∗B1 ) denote the (α, β) assignments on S1 after the application of Flooring_and_Ceiling with respect to Eqs. (14) and (15) respectively. The |aA | and |aB | distributions on surface S1 before and after the optimizations are shown in Fig. 15, while Table 4 lists the maximum data before and after the optimizations. Similarly, let (α∗A2 , β∗A2 ) and (α∗B2 , β∗B2 ) denote the (α, β) assignments on S2 after the application of Flooring_and_Ceiling with respect to Eqs. (14) and (15) respectively, and the corresponding optimization results are shown in Fig. 16, and Table 5. Results of GA for the objective Eq. (16) are also shown in Tables 4–5.

The test results also evidently demonstrate the effectiveness of the proposed Flooring_and_Ceiling algorithm in reducing the maximum of |aA |, |aB |, or their norm, mostly more than 40% on the two surfaces. It is though necessary to point out the simple fact that the improvement on one axis may unavoidably degrade the other. For instance, as revealed in Table 4, after our min–max |aA | on S1 , while the maximal |aA | is reduced from 15.84 to 9.33 (41.07%), the maximal |aB | increases from 17.78 to 20.03 (12.67%). However, this caveat does not simply rule out the usefulness and need of single axis optimization. This is because, on many five-axis machine tools, the physical requirements on the dynamics of different rotary axes are different, with different conditions, and often with priority. Also, for part surfaces with some special or biased shapes, the dynamic loading on one particular axis is usually much larger and critical than the others. Taking S2 for instance (see Table 5), after the minimization on maximal |aA |, the maximal |aA | is reduced from 31.88 to 14.17 (a 56% reduction), while in meantime the maximal |aB | has risen from 2.03 to 3.09, which though is still very small compared to |aA |. In this case, improving the dynamics of Aaxis is more critical and hence necessary. One point that needs to be re-addressed is that in our algorithms, the feed rate f is assumed to be constant and defined on the CC curves. However, the assumption of a constant feed rate is not necessary—it only serves to be convenient for comparing the results before and after the proposed optimizations. In real fiveaxis machining, the CNC usually applies the inverse time feed rate F , signifying that the move between the current position and the destination should be completed in a certain period of time [31]. Inverse time feed rate F will result in different feed rate f in each segment of the CC curve, which will affect the expressions of Eqs. (9) and (10). In our setting, the CC curves on the surface do not change before and after the optimizations. Conceivably, on any equal bases, regardless the varying feed rate f , our optimization algorithms always reduce the maximum angular acceleration on A, B axis or combined. 6. Summary Traditional five-axis tool path generation and optimization algorithms for complex surface milling machining are all computed in Part Coordinate System and in general disregard the final dynamic behaviors of the specific five-axis machine to be used. In this paper, we study the particular dynamic problem of how to reduce the maximum angular acceleration imposed on a rotary

P. Hu, K. Tang / Computer-Aided Design 43 (2011) 1693–1706

1703

Table 2 Results of the optimizations of workpiece setups for surface S1 . Unit: rad/s2 or (rad/s2 )2 .

maxi,j (|a |) B max ∑ i,j (|Aa 2|) (( a ) + (aB )2 ) i,j A

r = r0

r = rA∗

r = rB∗

r = rA∗+B

15.8382 17.7772 1.8751 × 105

11.2795 20.5849 2.1459 × 105

19.4060 10.1421 1.6318 × 105

19.4164 10.2371 1.6211 × 105

∗ ∗ Fig. 14. Angular accelerations on surface S2 : (a) original |aA | and |aB | with r = r0 ; (b) |aA | and |aB | with r = rA2 ; (c) |aA | and |aB | with r = rB2 .

Table 3 Results of the optimizations of workpiece setups for surface S2 . Unit: rad/s2 or (rad/s2 )2 .

maxi,j (|aA |) maxi,j (|aB |) ∑ A 2 B 2 i,j ((a ) + (a ) )

r = r0

r = rA∗

r = rB∗

r = rA∗+B

31.8786 2.0265 1.8751 × 105

1.6852 30.5686 1.5355 × 105

31.1462 1.7239 1.5545 × 105

2.3797 30.5162 1.4357 × 105

Fig. 15. Angular accelerations on surface S1 : (a) original |aA | and |aB |; (b) the |aA | and |aB | when (α(S1 ), β(S1 )) = α∗A1 , β∗A1 ; (c) the |aA | and |aB | when (α(S1 ), β(S1 )) =





 ∗  αB1 , β∗B1 .

Table 4 Results of the optimizations after the Flooring_and_Ceiling operations on surface S1 . Unit: rad/s2 or (rad/s2 )2 .

maxi,j (|a |) maxi,j (|aB |) ∑ A 2 B 2 i,j ((a ) + (a ) ) A

((aAi,j )2 + (aBi,j )2 )

Original

min–max|aA |

min–max|aB |

min

15.8382 17.7772 1.8751 × 105

9.3333 20.0287 1.9541 × 105

16.7981 11.0847 2.1325 × 105

11.0035 17.3265 1.5751 × 105

∑

i,j

1704

P. Hu, K. Tang / Computer-Aided Design 43 (2011) 1693–1706

Fig. 16. Angular accelerations on surface S2 : (a) original |aA | and |aB |; (b) the |aA | and |aB | when (α(S2 ), β(S2 )) = α∗A2 , β∗A2 ; (c) the |aA | and |aB | when (α(S2 ), β(S2 )) =







αB2 , βB2 . ∗

∗



Table 5 Results of the optimizations after the Flooring_and_Ceiling operations on surface S2 . Unit: rad/s2 or (rad/s2 )2 .

maxi,j (|a |) B max ∑ i,j (|Aa 2|) (( a ) + (aB )2 ) i ,j A

((aAi,j )2 + (aBi,j )2 )

Original

min–max|aA |

min–max|aB |

min

31.8786 2.0265 1.1924 × 105

14.1689 3.0924 1.0007 × 105

30.8587 0.9471 1.6161 × 105

15.8671 2.652 0.8521 × 105

axis of the machine by a tool path. Two approaches are proposed. In the first approach, one takes as input an already computed tool path (though in PCS) and, by using the Genetic Algorithm technique, finds the best workpiece setup on the machine’s table to minimize the maximum angular acceleration. In the second, with the workpiece setup fixed, we first compute an initial tool path and subsequently, with the CC curves kept the same, modify the tilt and yaw angle of the tool at the CC points to minimize the maximum angular acceleration. To ensure the machined surface quality, i.e., the surface must be free of local gouging and under the maximum cusp-height limit, a concept called Domain Geometric Constraint is introduced, which defines the allowable domain of the tilt and yaw angle. An iterative algorithm is then proposed to realize this constrained optimization. Preliminary experiments in computer simulation have shown the effectiveness of both approaches in the prescribed optimization problem, sometimes achieving more than 50% reduction in the maximum angular acceleration. Regarding the issues and potential future research subjects of the presented solutions, the following are of particular interest. Part surface representation. Currently, our initial path generation algorithm (which relies on the parameters of the surface to advance the CC curves) is limited to only parametric surfaces. In most practical situations, e.g., automotive parts, a part surface is often a composite of many parametric surfaces, or even a triangular mesh. For such a surface, a strategy must be determined for deciding the first CC curve and the subsequent CC curves. One plausible idea is to, starting from the boundary curve, use the geodesic offset curves on the surface as subsequent CC curves. A careful analysis however must be conducted on how to advance the CC curves and also how a DGC should be constructed this time, in order to maintain the maximum cusp-height. Modification of CC curves. Angular accelerations depend on both CC curves and the tool orientations assigned along them. While

∑

i ,j

conceivably the influence of the tool orientations on the angular acceleration may be more direct and larger than that of CC curves, it is still worth investigating the effectiveness of modifying the CC curves to reduce the angular acceleration, especially if other types of optimization, e.g., minimizing the total length of the CC curves, are desired along with the optimization on dynamics. Nevertheless, due to the complex geometric relationship between the CC point, the tool axis, the surface normal, and the tool tip (reference) point, most existing tool path generation algorithms decouple the planning of CC curves from their associated tool orientations—when optimizing the tilt and yaw angles, the CC curves are fixed and, conversely, when optimizing the CC curves, the tilt and yaw angles are constants. How to consider the two at the same time, especially when dynamics of the machine is involved, remains a challenging task. Domain of constraint of tilt and yaw angles. In our current solutions, the allowable domain for the modification of the tilt and yaw angle at a CC point is defined per the constraints of local-gouging and cusp-height, which are purely geometric (thus the name ‘‘domain of geometric constrain’’). It is natural to extend it to include other non-geometric considerations, particularly those pertinent to machining, e.g., the cutting force. In addition, the global collision is currently not considered and should be included as well. Finally, methods for avoiding singular configurations of the five-axis machine tools need to be addressed and integrated with the proposed methods, before the application of our methods in real machining. Another point is that, a better initial tool path generation algorithm should be explored, especially in the way how CC curves are generated. Advancing along the iso-parametric curves makes it ultra-convenient for planning CC curves. However, in order to handle non-parametric part surfaces, and also for further improvement in the stipulated optimization objectives, it ought be asked if a more intelligent non-iso-parametric strategy for planning CC curves can be devised, akin to the idea in [22] for spindle–tilt type five-axis machines.

P. Hu, K. Tang / Computer-Aided Design 43 (2011) 1693–1706

1705

Fig. A. Effective cutting ellipse. Fig. B. (a) (αi , βi ) ∈ DGCi . (b) (αi , βi ) ̸∈ DGCi and the modified search routes.

Acknowledgement This work is supported by Hong Kong RGC General Research Fund 620409. Appendix A. Mathematical formulation for calculation effective cutting radius for flat-end cutter [20] In the LCS of the CC point, the effective cutting radius of the flatend cutter can be expressed as the function of the tilt angle α and yaw angle β . In the plane of n − k, the mapping of the end of the cutter is an ellipse, as displayed in Fig. A, and it is called the effective cutting ellipse. The effective cutting radius of the flat-end cutter is the curvature of the ellipse of the CC point, which can be described as re = t 2 r 2



1 + tan2 θ

3/2 (A.1)

t 2 + r 2 tan θ 2

where t = r sin α · cos β,

θ = tan−1 (tan α · sin β), r is the radius of the cutter.

and

Appendix B. Algorithm for searching (αi , βi ) satisfying

ϕi (αi , βi ) = ϕiλ

WHILE (|ϕi (αi , βi ) − ϕiλ | > E ) and (running time limit has not been reached) m = 0.5(αi + αi∗ ) n = 0.5(βi + βi∗ ) λ IF sϕ  i (∗m, n∗ ) > sϕi αi , βi = (m, n) ELSE (αi , βi ) = (m, n) END END In the above algorithm, simple bisection method is used, and the termination criteria are the running time limit and a desired tolerance E (10−6 in our current implementation). In our tests, it is found that in most cases indeed a (αi , βi ) is returned within DGCi , e.g., Fig. B(a). In case this search fails, either because the found (α  i′, βi) ̸∈ DGCi or it is not found at all, we insert a new point αi , 0 in the middle of DGCi , and then search along the two lines as shown in Fig. B(b). These two lines are virtually always inside DGCi (unless the part surface is extremely complicated) and in general the (αi , βi ) satisfying ϕi (αi , βi ) = ϕiλ will be found on them.

References [1] Han Z, Yang DCH. Iso-phote based tool-path generation for machining freeform surfaces. Journal of Manufacturing Science and Engineering 1999;121: 656. [2] Tournier C, Duc E. A surface based approach for constant scallop height tool-path generation. The International Journal of Advanced Manufacturing Technology 2002;19(5):318–24. [3] Chen T, Ye P. A tool path generation strategy for sculptured surfaces machining. Journal of Materials Processing Technology 2002;127(3):369–73. [4] Anotaipaiboon W, Makhanov SS. Tool path generation for five-axis NC machining using adaptive space-filling curves. International Journal of Production Research 2005;43(8):1643–65. [5] He W, Lei M, Bin H. Iso-parametric CNC tool path optimization based on adaptive grid generation. The International Journal of Advanced Manufacturing Technology 2009;41(5):538–48. [6] Lavernhe S, Tournier C, Lartigue C. Optimization of 5-axis high-speed machining using a surface based approach. Computer-Aided Design 2008; 40(10–11):1015–23. [7] Ho MC, Hwang YR, Hu CH. Five-axis tool orientation smoothing using quaternion interpolation algorithm. International Journal of Machine Tools and Manufacture 2003;43(12):1259–67. [8] Jun CS, Cha K, Lee YS. Optimizing tool orientations for 5-axis machining by configuration-space search method. Computer-Aided Design 2003;35(6): 549–66. [9] Wang N, Tang K. Automatic generation of gouge-free and angular-velocity compliant five-axis tool path. Computer-Aided Design 2007;39(10):841–52. [10] Chen KH. Investigation of tool orientation for milling blade of impeller in five-axis machining. The International Journal of Advanced Manufacturing Technology 2011;1–10. [11] Anotaipaiboon W, Makhanov SS, Bohez ELJ. Optimal setup for five-axis machining. International Journal of Machine Tools and Manufacture 2006; 46(9):964–77. [12] Munlin M, Makhanov SS, Bohez ELJ. Optimization of rotations of a five-axis milling machine near stationary points. Computer-Aided Design 2004;36(12): 1117–28. [13] Fleisig R, Spence A. A constant feed and reduced angular acceleration interpolation algorithm for multi-axis machining. Computer Aided Design 2001;33(1):1–16. [14] Sorby K. Inverse kinematics of five-axis machines near singular configurations. International Journal of Machine Tools and Manufacture 2007;47(2):299–306. [15] Affouard A, Duc E, Lartigue C, Langeron JM, Bourdet P. Avoiding 5-axis singularities using tool path deformation. International Journal of Machine Tools and Manufacture 2004;44(4):415–25. [16] Bohez E. Five-axis milling machine tool kinematic chain design and analysis. International Journal of Machine Tools and Manufacture 2002;42(4):505–20. [17] Bohez ELJ. Compensating for systematic errors in 5-axis NC machining. Computer-Aided Design 2002;34(5):391–403. [18] Lee RS, She CH. Developing a postprocessor for three types of five-axis machine tools. The International Journal of Advanced Manufacturing Technology 1997; 13(9):658–65. [19] Suresh K, Yang C. Constant scallop-height machining of free-form surfaces. Transactions of the ASME. Journal of Engineering for Industry USA 1994; 116(2):253–9. [20] Lo C. Efficient cutter-path planning for five-axis surface machining with a flatend cutter. Computer-Aided Design 1999;31(9):557–66. [21] Agrawal RK, Pratihar D, Roy Choudhury A. Optimization of CNC isoscallop free form surface machining using a genetic algorithm. International Journal of Machine Tools and Manufacture 2006;46(7–8):811–9. [22] Wang N, Tang K. Five-axis tool path generation for a flat-end tool based on iso-conic partitioning. Computer-Aided Design 2008;40(12):1067–79. [23] Lee YS. Admissible tool orientation control of gouging avoidance for 5-axis complex surface machining. Computer-Aided Design 1997;29(7):507–21.

1706

P. Hu, K. Tang / Computer-Aided Design 43 (2011) 1693–1706

[24] Morishige K, Takeuchi Y, Kase K. Tool path generation using C -space for 5-axis control machining. Journal of Manufacturing Science and Engineering 1999; 121:144. [25] Ye T, Xiong CH. Geometric parameter optimization in multi-axis machining. Computer-Aided Design 2008;40(8):879–90. [26] Chiou JCJ, Lee YS. Optimal tool orientation for five-axis tool-end machining by swept envelope approach. Journal of Manufacturing Science and Engineering 2005;127:810. [27] Kim T, Sarma SE. Toolpath generation along directions of maximum kinematic performance: a first cut at machine-optimal paths. Computer-Aided Design 2002;34(6):453–68.

[28] Chiou JCJ. Floor, wall and ceiling approach for ball-end tool pocket machining. Computer-Aided Design 2005;37(4):373–85. [29] Castagnetti C, Duc E, Ray P. The domain of admissible orientation concept: a new method for five-axis tool path optimisation. Computer-Aided Design 2008;40(9):938–50. [30] Debout P, Chanal H, Duc E. Tool path smoothing of a redundant machine: application to automated fiber placement. Computer-Aided Design 2011; 43(2):122–32. [31] Karlo Apro. Secrets of five-axis machining. New York: Industrial Press; 2008. [32] Lee YS. Non-isoparametric tool path planning by machining strip evaluation for 5-axis sculptured surface machining. Computer-Aided Design 1998;30(7): 559–70.