# Modeling of Drilling Forces Based on Twist Drill Point Angles Using Multigene Genetic Programming

Authors: Myong-Il Kim; Ping Zou
Journal: Mathematical Problems in Engineering
 (2016)

---

## Abstract

The mathematical model was developed for predicting the influence of the drill point angles on the cutting forces in drilling with the twist drills, which was used to optimize those angles for reducing drilling forces. The approach was based on multigene genetic programming, for the training data, the grinding tests of twist drill were firstly conducted for the different drill point angles in Biglide parallel machine, and then drilling tests were performed on carbon fiber reinforced plastics using the grinded drills. The effectiveness of the proposed approach was verified through comparing with published data. It was found that the proposed model agreed well with the experimental data and was useful for improving the performance of twist drill.


---

## Body

## 1. Introduction

The drilling process is one of the most important and basic metal cutting operations because of its large portion of overall machining operations. Moreover, the drilling problems can result in costly production waste because many drilling operations are usually in the final steps for fabricating a part. In drilling process, the drilling forces are one of the primary responses; therefore, predicting the influences of the drill point geometry on the cutting forces will provide many platforms needed for the evaluation and comparison of new drill designs. The geometry and cutting mechanics of twist drill have been well studied over many years; many models have been developed by many researchers for predicting torque and thrust force in drilling [1–4].Early drilling models were developed based on the analytical mathematical approach using the mechanics of cutting analysis and the shear plane approach, where the drilling mechanism was analytically described by complicated equations in 3D space; 2D projected geometry was used instead in many cases, in order to reduce the amount of calculation [5, 6]. In these studies, the modeling of drilling process was conducted using a series of oblique cutting slices; after that, it was expanded to the more detailed modeling of material removal in both the cutting lip and chisel edge regions for predicting the drill temperatures and simulating arbitrary drill geometries [7, 8]. More recent developments in drilling models have utilized either a mechanistic or a finite element method. The finite element method was used for determining drilling torque and thrust force and for predicting drill heat flux, temperatures, and the thermal distortion of drill holes based on the Lagrangian and Eulerian methods [3, 9–11]. Chandrasekharan et al. [7, 12] and Hamade et al. [13] developed mechanistic drilling models, where an extensive amount of experiments took place and the results were stored in databases so that different parameters could be used for experimentally derived equations. Mechanistic approach typically involves determining the forces acting along the tool rake face from which the resultant force projected along any desired orientation may be calculated [14–16].As stated above, although a number of approaches for predicting the drilling forces were developed using the different cutting conditions, there is little research to evaluate the influences of the drill point geometry on the cutting forces.The drilling forces were found to be mainly dependent on not only the cutting conditions but also the drill point geometry [12, 17]. Moreover, as the point geometry of the twist drill is completely featured with the drill point angles such as the semipoint angle, the chisel edge angle, and the relief angle, the drilling forces depend naturally on the drill point angles. Therefore, it is very significant to formulate the mathematical model of drilling forces in terms of drill point angles for designing a high performance twist drill.For formulating this model, various advanced mathematical modeling methods could be applied such as support vector regression (SVR), artificial neural network (ANN), and multigene genetic programming (MGGP), which are comprehensively used in engineering problems because of their powerful prediction ability [18, 19]. Amongst these methods, MGGP possesses the ability to evolve the model structure and its coefficients, which shows the advantages in predicting accuracy and reliability compared to the others in many applications [20, 21].In this paper, the experimental model is presented to predict the influence of twist drill point angles on the cutting forces to be presented based on MGGP method. In order to acquire the training data, the grinding tests and measurements of the drill point angles are conducted on Biglide parallel machine [22] and the drilling tests are performed using carbon fiber reinforced plastics (CFRP) work piece, which is widely used in aerospace industry due to its ability to resist corrosion and withstand high loads while reducing the weight of the structural parts [23–25]. The twist drills with different drill point angles are used in this experiment. The model is verified by comparing predicted drilling force and torque with published data and measured.


## 2. Preparing of Training Data

In order to prepare the training input data, firstly, the experiments were conducted to grind twenty twist drills with different drill point angles with Biglide parallel machine; the structure is as shown Figure1 [26].Figure 1
Structure chart of Biglide parallel machine: (1) sliding carriage, (2) ball screw, (3) slider, (4) rod, (5) drill, (6) moving platform, (7) grinding wheel, (8) work table, and (9) worm gearing.Biglide parallel machine has three translations, which is comprised of parallelogram links to connect the moving platform with the slide units by revolute joint. The work table moves up and down by the unit of worm and worm wheel. The linear movement of each slider is provided with the lead ball screw system and actuated by a stepping motor. The moving platform moves along the lengthwise direction by another ball screw and the sliding carriage under the condition of driving motor, which is a decoupled motion. The different drill point angles are obtained by controlling the velocity of the slide, the angular velocity of the drill, the angle between the axis of drill and the axis of grinding wheel, and the distance between two slides. After the twist drills are grinded out with the designated different angles, the main drill point angles are measured for each twist drill as proposed by Kim and Zou [22]; the measurement values of the drill point angles are shown in Table 1.Table 1
Measured values of drill point angles.
Number
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
rd
3
3
3
3
3
3
3
3
3
3
3
3
3
3
3
3
3
3
3
3
α
7
8
9
12
12
15
16
16
16
19
23
31
33
34
34
37
37
38
44
50
ψ
54
53
67
73
51
61
53
55
57
65
40
41
60
36
43
45
49
39
55
48
ρ
59
59
59
59
59
59
59
59
59
59
59
59
59
59
59
59
59
59
59
59In Table1, rd means the drill radius mm, α the relief angle °, ψ the chisel edge angle °, and ρ the semipoint angle °; the values of the semipoint angle are constant 59°, because the angle between the drill axis and the grinding wheel axis is fixed as a constant angle in Biglide parallel machine during grinding of drill point.The drilling experiments are carried out by Z3040 radical drilling machine in which the thrust forces and the torque during the drilling operations are measured by a Kistler model 9257B quartz 3-component dynamometer. In this series of tests, the spindle speed is fixed at 400 rpm with the feed rate of 0.1 mm/rev; the material of the work piece is CFRP. All drilling tests are conducted under dry conditions, that is, without any coolant. The thrust force and torque are given by DynoWare software. For example, those in drilling test using the twist drill corresponding to number 7 in Table1 are shown in Figures 2 and 3.Figure 2
The thrust force measured for number 7 twist drill.Figure 3
The torque measured for number 7 twist drill.The mean values of the thrust force and torque calculated by DynoWare software areFzm=97.6N and Mzm=0.86Nm, respectively.For convenience [27], a resultant quantity FM is introduced instead of the thrust force and torque; FM=Fzm+0.5Mzm=98.03. The relevant FM values in drilling tests using twist drills in Table 1 are as shown in Table 2; like this, the training data are ready for MGGP model.Table 2
F
M values.
Number
1
2
3
4
5
6
7
8
9
10
r
3
3
3
3
3
3
3
3
3
3
α
7
8
9
12
12
15
16
16
16
19
ψ
54
53
67
73
51
61
53
55
57
65
ρ
59
59
59
59
59
59
59
59
59
59
F
M
120.96
86.22
170.12
308.88
145.12
151.39
98.03
186.49
239.82
171.76
Number
11
12
13
14
15
16
17
18
19
20
r
3
3
3
3
3
3
3
3
3
3
α
23
31
33
34
34
37
37
38
44
50
ψ
40
41
60
36
43
45
49
39
55
48
ρ
59
59
59
59
59
59
59
59
59
59
F
M
97.64
87.73
107.28
186.16
75.98
85.03
103.39
113.91
303.79
128.05


## 3. MGGP-Based Modeling

As shown in Table2, the relationship of the drill point angles of the twist drill and drilling forces is very irregular and complex to make the analytical model; it is needed to use the powerful predicting methods for empirical modeling. The genetic programming (GP) method has been widely used to predict the behavior [28–30], of which structure is shown in Figure 4.Figure 4
Structure of GP model.The performance of a GP model is evaluated on the fitness function, that is, root mean square error (RMSE):(1)RMSE=∑i=1NPi-Ai2N,where Pi is the value of ith sample data predicted by the GP model, Ai is the actual value, and N is the number of sample data. The subtree of the crossover operation is shown in Figure 5, where two branches are selected randomly to be swapped.Figure 5
Subtree of crossover operation in GP model.Figure6 shows the subtree of the mutation operation, where the node (terminal or functional) is selected randomly to be replaced by branch of the new generated random tree.Figure 6
Subtree of mutation operation in GP model.MGGP combines the ability of selecting the model structure in GP with the parameter estimation power of the classical regression by using “multigene.” The MGGP model is as shown in Figure7, where d0 is the bias of the offset term and d1 and d2 are the weights for gene1 and gene2, respectively.Figure 7
MGGP model.The MGGP algorithm is outlined as follows:(a)
Set the initial parameters such as the function and terminal set, the number of generations, the population size, and the maximum depth of gene.(b)
Generate the initial population of genes randomly.(c)
Construct the models by combining set of genes using least squares method.(d)
Evaluate the performance of models based on the fitness function.(e)
Conduct the genetic operations and construct the new population.(f)
Cross-check the models performance with the termination criterion; if not satisfied then go to step (e); else select the model as the best one.In this paper, our goal is to get a relationship model from the experimental data in Table2 by using MGGP; the relief angle and chisel edge angle are set as input variables and FM is set as output. The fitness function is defined as RMSE in (1); MGGP-based model is selected based on minimum fitness function and its implementation and performance are discussed in the next section.


## 4. Application of the Proposed Modeling Method

### 4.1. Comparing Proposed Modeling Method with Published Data

For evaluating the performance of MGGP-based method, the MGGP model is firstly obtained using the experimental data proposed by Chandrasekharan et al. [12] in Table 3.Table 3
Predicted (bold) and experimental (italic) thrust forces.
Feed (mm/rev)
Speed (rpm)
Diameter (mm)
Point angle
Web thickness (mm)
Pilot hole diameter (mm)
⁢Entire drill thrustforce (N)
0.229
400
15.9
118
2.3
3.2
2247.01
2228.64
0.102
200
15.9
118
2.3
3.2
1279.97
1352.70
0.102
800
15.9
118
2.3
3.2
1157.96
1149.73
0.254
400
12.7
118
2.3
2.8
2001.38
1873.23
0.102
400
9.5
118
1.5
2.4
781.99
811.97
0.178
400
9.5
118
1.5
2.4
1184.69
1131.71
0.178
400
12.7
80
1.9
2.8
1221.16
1376.90
0.102
400
12.7
160
1.9
2.8
1765.58
1893.29
0.102
400
15.9
135
2.3
4.4
997.64
868.69
0.178
400
15.9
135
2.3
4.4
1528.14
1542.99The parameter setting is very important for the generalization ability of the MGGP model; the parameters selected by trial-and-error approach are shown in Table4.Table 4
Parameters for MGGP (1).
Parameters
Setting values
Population size
400
Number of generations
150
Lexicographic selection
True
Max tree depth
4
Max nodes per tree
Inf
Functional set
times, minus, plus, pdivide, sin, cos, exp, mypower
Terminal set
Number of inputs: 6, constants range: [−10 10]
Max genes
3
Using fitness function
regressmulti_fitfun.m
Crossover probability
0.85
Reproduction probability
0.10
Mutation probability
0.05As Table4 shows, the elements in the functional set evolve a variety of nonlinear mathematical models and those in the terminal set consist of the 6 input variables (i.e., feed, speed, diameter, point angle, web thickness, and pilot hole diameter) and random constants are in the range -1010 so as to take into account the experimental errors (±10%) in the data collection.The best fitness function is given at generation 145 and the predicted function based on MGGP for the thrust force is as follows:(2)fx1,x2,x3,x4,x5,x6=4.82391×102-3.58052×102×x5expx2-cos⁡x1-x1x3-1.83528×102×cos⁡2.687335+x4×x1x5+5.16001×102×cos⁡x3x4x5-x6.The thrust force calculated on each data point by (2) is shown in Figure 8; the approximation proportion of the predicted values to the experimental ones is 0.985906.Figure 8
Thrust forces predicted by MGGP model.Figure9 shows the results of the thrust force predicted by Chandrasekharan and MGGP model, respectively; it is clearly seen that the latter is closer to the ideal line than the former.Figure 9
Comparison of thrust forces predicted by Chandrasekharan and MGGP model.For comparing the performance of the proposed method with Chandrasekharan data, the three following error statistics are introduced:(3)RMSPE=1N∑i=1NPFi-AFi2(4)MAPE=1N∑i=1NPFi-AFi(5)R=∑PFi-PFi¯AFi-AFi¯∑PFi-PFi¯2∑AFi-AFi¯2,where AFi and PFi are actual and predicted values of the thrust force, respectively, and AFi¯ and PFi¯ are the average values of them, respectively. MAPE can express the estimation accuracy generally, and RMSPE is the fundamental scale for comparing the accuracies of different interpolation methods. Square of correlation coefficient R is used to evaluate the performance of the predicting method. The values of them are shown in Table 5; the results show significantly that MGGP method has better performance.Table 5
Assessment of two models by error statistics.
Model
RMSPE
MAPE
R
Chandrasekharan
91.0488
71.932
0.97857
MGGP
52.2156
39.6891
0.99293


### 4.2. Modeling of Drilling Forces Based on Drill Point Angles

For the improved new type of twist drill, the minimization of the thrust forces and torque leads to an improvement in performance by reducing drill deflection due to thrust loading, by lowering the power required for the drilling operation, and by improving drill life. In order to make the model of drill point angles for minimizing drilling forces, in this paper, MGGP method is adopted using the measurement data of Table1 as training data. As shown in Table 6, MGGP parameters are selected by trial-and-error approach in previous section; the difference is to select two input variables (i.e., relief angle and chisel edge angle) for terminal set.Table 6
Parameters for MGGP (2).
Parameters
Setting values
Population size
300
Number of generations
200
Tournament size
4
Lexicographic selection
True
Max tree depth
8
Max nodes per tree
Inf
Functional set (F)
TIMES MINUS PLUS PDIVIDE SIN COS EXP MYPOWER
Terminal set (T)
Number of inputs: 2, constants range: [−10 10]
Max genes
5
Using fitness function
regressmulti_fitfun.m
Crossover probability
0.85
Reproduction probability
0.10
Mutation probability
0.05Figure10 shows the convergence course of the fitness function; the best fitness function is given at 146th generation.Figure 10
The convergence of the fitness function.The predicting model is given as(6)fx1,x2=-7.10886×103+3.50915×104×3x2+2×8.326996x22+1.97076×102×x2+sin⁡x1-2.85913×10×x2+sin⁡x1-cos⁡8.256093-x2×x2+sin⁡x1-8.02450×10-1×x2+sin⁡sin⁡sin⁡x1×8.256093+sin⁡sin⁡x1+x2-7.92810×10×cos⁡cos⁡sin⁡x2-x2-sin⁡x1-2.79330×8.228307+sin⁡sin⁡x1×x2-sin⁡sin⁡sin⁡x1sin⁡x2.Equation (6) shows the formulated model to predict FM as function of drill point relief angle x1 and chisel edge angle x2 based on MGGP. The predicted and actual values of FM on data point are shown in Figure 11; it is clearly seen that the predicted values agree very well with the actual measurement ones.Figure 11
The actual and predicted values ofFM.For evaluating the generalization performance of the proposed model, 4 sets of data are ready for checking as shown in Table7, which have never been used for training data and are almost out of collection range for training data.Table 7
F
M values for checking data.
Number
1
2
3
4
r
3
3
3
3
α
28
61
67
41
ψ
25
35
83
52
ρ
59
59
59
59
F
M
910.32
220.46
380.32
176.38Compared to the actual values ofFM for checking data in Table 7, the corresponding values predicted using (6) are shown in Figure 12. In this case, the value of estimation accuracy MAPE is 32.3940 (the original 6.5653), which is calculated for using (4).Figure 12
The actual and predicted values ofFM for performance check.Consequently, the estimation accuracy in checking data is less than in the original, but it is well known that the proposed model could be very useful in large range.In the optimization of the drill point relief angle and chisel edge angle for reducing drilling forces, the objective is limited to minimizing the thrust forces and torque. The objective function is chosen asU=Ththrust  force+0.5Ttorque [27] for minimizing the thrust forces and torque; therefore, the optimization above can be just defined as a problem to find the optimal relief angle and chisel edge angle for minimizing of FM.Mathematically, this optimization problem is formulated as(7)minimizefXoverX=x1,x2Tsubject  tog1:0°<x1<90°g2:0°<x2<90°.This optimization problem is easily solved using (6); the best solution is found as(8)X=34°,46°T.


## 4.1. Comparing Proposed Modeling Method with Published Data

For evaluating the performance of MGGP-based method, the MGGP model is firstly obtained using the experimental data proposed by Chandrasekharan et al. [12] in Table 3.Table 3
Predicted (bold) and experimental (italic) thrust forces.
Feed (mm/rev)
Speed (rpm)
Diameter (mm)
Point angle
Web thickness (mm)
Pilot hole diameter (mm)
⁢Entire drill thrustforce (N)
0.229
400
15.9
118
2.3
3.2
2247.01
2228.64
0.102
200
15.9
118
2.3
3.2
1279.97
1352.70
0.102
800
15.9
118
2.3
3.2
1157.96
1149.73
0.254
400
12.7
118
2.3
2.8
2001.38
1873.23
0.102
400
9.5
118
1.5
2.4
781.99
811.97
0.178
400
9.5
118
1.5
2.4
1184.69
1131.71
0.178
400
12.7
80
1.9
2.8
1221.16
1376.90
0.102
400
12.7
160
1.9
2.8
1765.58
1893.29
0.102
400
15.9
135
2.3
4.4
997.64
868.69
0.178
400
15.9
135
2.3
4.4
1528.14
1542.99The parameter setting is very important for the generalization ability of the MGGP model; the parameters selected by trial-and-error approach are shown in Table4.Table 4
Parameters for MGGP (1).
Parameters
Setting values
Population size
400
Number of generations
150
Lexicographic selection
True
Max tree depth
4
Max nodes per tree
Inf
Functional set
times, minus, plus, pdivide, sin, cos, exp, mypower
Terminal set
Number of inputs: 6, constants range: [−10 10]
Max genes
3
Using fitness function
regressmulti_fitfun.m
Crossover probability
0.85
Reproduction probability
0.10
Mutation probability
0.05As Table4 shows, the elements in the functional set evolve a variety of nonlinear mathematical models and those in the terminal set consist of the 6 input variables (i.e., feed, speed, diameter, point angle, web thickness, and pilot hole diameter) and random constants are in the range -1010 so as to take into account the experimental errors (±10%) in the data collection.The best fitness function is given at generation 145 and the predicted function based on MGGP for the thrust force is as follows:(2)fx1,x2,x3,x4,x5,x6=4.82391×102-3.58052×102×x5expx2-cos⁡x1-x1x3-1.83528×102×cos⁡2.687335+x4×x1x5+5.16001×102×cos⁡x3x4x5-x6.The thrust force calculated on each data point by (2) is shown in Figure 8; the approximation proportion of the predicted values to the experimental ones is 0.985906.Figure 8
Thrust forces predicted by MGGP model.Figure9 shows the results of the thrust force predicted by Chandrasekharan and MGGP model, respectively; it is clearly seen that the latter is closer to the ideal line than the former.Figure 9
Comparison of thrust forces predicted by Chandrasekharan and MGGP model.For comparing the performance of the proposed method with Chandrasekharan data, the three following error statistics are introduced:(3)RMSPE=1N∑i=1NPFi-AFi2(4)MAPE=1N∑i=1NPFi-AFi(5)R=∑PFi-PFi¯AFi-AFi¯∑PFi-PFi¯2∑AFi-AFi¯2,where AFi and PFi are actual and predicted values of the thrust force, respectively, and AFi¯ and PFi¯ are the average values of them, respectively. MAPE can express the estimation accuracy generally, and RMSPE is the fundamental scale for comparing the accuracies of different interpolation methods. Square of correlation coefficient R is used to evaluate the performance of the predicting method. The values of them are shown in Table 5; the results show significantly that MGGP method has better performance.Table 5
Assessment of two models by error statistics.
Model
RMSPE
MAPE
R
Chandrasekharan
91.0488
71.932
0.97857
MGGP
52.2156
39.6891
0.99293


## 4.2. Modeling of Drilling Forces Based on Drill Point Angles

For the improved new type of twist drill, the minimization of the thrust forces and torque leads to an improvement in performance by reducing drill deflection due to thrust loading, by lowering the power required for the drilling operation, and by improving drill life. In order to make the model of drill point angles for minimizing drilling forces, in this paper, MGGP method is adopted using the measurement data of Table1 as training data. As shown in Table 6, MGGP parameters are selected by trial-and-error approach in previous section; the difference is to select two input variables (i.e., relief angle and chisel edge angle) for terminal set.Table 6
Parameters for MGGP (2).
Parameters
Setting values
Population size
300
Number of generations
200
Tournament size
4
Lexicographic selection
True
Max tree depth
8
Max nodes per tree
Inf
Functional set (F)
TIMES MINUS PLUS PDIVIDE SIN COS EXP MYPOWER
Terminal set (T)
Number of inputs: 2, constants range: [−10 10]
Max genes
5
Using fitness function
regressmulti_fitfun.m
Crossover probability
0.85
Reproduction probability
0.10
Mutation probability
0.05Figure10 shows the convergence course of the fitness function; the best fitness function is given at 146th generation.Figure 10
The convergence of the fitness function.The predicting model is given as(6)fx1,x2=-7.10886×103+3.50915×104×3x2+2×8.326996x22+1.97076×102×x2+sin⁡x1-2.85913×10×x2+sin⁡x1-cos⁡8.256093-x2×x2+sin⁡x1-8.02450×10-1×x2+sin⁡sin⁡sin⁡x1×8.256093+sin⁡sin⁡x1+x2-7.92810×10×cos⁡cos⁡sin⁡x2-x2-sin⁡x1-2.79330×8.228307+sin⁡sin⁡x1×x2-sin⁡sin⁡sin⁡x1sin⁡x2.Equation (6) shows the formulated model to predict FM as function of drill point relief angle x1 and chisel edge angle x2 based on MGGP. The predicted and actual values of FM on data point are shown in Figure 11; it is clearly seen that the predicted values agree very well with the actual measurement ones.Figure 11
The actual and predicted values ofFM.For evaluating the generalization performance of the proposed model, 4 sets of data are ready for checking as shown in Table7, which have never been used for training data and are almost out of collection range for training data.Table 7
F
M values for checking data.
Number
1
2
3
4
r
3
3
3
3
α
28
61
67
41
ψ
25
35
83
52
ρ
59
59
59
59
F
M
910.32
220.46
380.32
176.38Compared to the actual values ofFM for checking data in Table 7, the corresponding values predicted using (6) are shown in Figure 12. In this case, the value of estimation accuracy MAPE is 32.3940 (the original 6.5653), which is calculated for using (4).Figure 12
The actual and predicted values ofFM for performance check.Consequently, the estimation accuracy in checking data is less than in the original, but it is well known that the proposed model could be very useful in large range.In the optimization of the drill point relief angle and chisel edge angle for reducing drilling forces, the objective is limited to minimizing the thrust forces and torque. The objective function is chosen asU=Ththrust  force+0.5Ttorque [27] for minimizing the thrust forces and torque; therefore, the optimization above can be just defined as a problem to find the optimal relief angle and chisel edge angle for minimizing of FM.Mathematically, this optimization problem is formulated as(7)minimizefXoverX=x1,x2Tsubject  tog1:0°<x1<90°g2:0°<x2<90°.This optimization problem is easily solved using (6); the best solution is found as(8)X=34°,46°T.


## 5. Conclusions

MGGP-based modeling method was proposed for predicting the drilling forces. The effectiveness of the proposed method was verified by comparing its performance with the published data [12]. The grinding tests were conducted using Biglide parallel machine and the drill point angles were measured to prepare the input values. The thrust forces and torques were measured during drilling on work piece of CFRP material to prepare the output values. Using these input and output values as the training data, the model for predicting thrust forces and torque was proposed based on MGGP to obtain the optimal drill point angles for minimizing the thrust forces and torque. This work can be further extended to analysis and design of the reasonable geometry of twist drill under the different cutting conditions.



---

*Source: 6749182-2016-01-19.xml*