Analysis and Prediction of Surface Topographyin Peripheral Milling Based on Workpiece Vibration and Milling-Tool Structure

LIU Yang;QIN Guohua;WU Zhuxi;LOU Weida;LAI Xiaochun

Acta Armamentarii ›› 2023, Vol. 44 ›› Issue (7) : 2132-2146. DOI: 10.12382/bgxb.2022.0289

Analysis and Prediction of Surface Topographyin Peripheral Milling Based on Workpiece Vibration and Milling-Tool Structure

Author information +
History +

Abstract

In the milling process, vibrations can lead to variations in the structure of machined parts, significantly influencing the machining quality of the finished surface. Therefore, investigating the milling formation method of the machined surface morphology is of great significance for predicting surface roughness and selecting reasonable machining parameters. First, by equivalently transforming the contact mode between workpiece and fixture into a linear spring damping system, the vibration differential equation of the workpiece is derived for the milling process by using the energy method. By using the coordinate transformation method, the modal analysis solution method is proposed for the vibration differential equations with constraint conditions to obtain the workpiece position deviation. Secondly, a method for determining the tool contact point and establishing the motion trajectory equation for the cutting edge in peripheral milling is presented. This enables the development of a discrete algorithm for surface topography by incorporating the workpiece position deviation. Experimental validations are conducted to verify the proposed algorithm for peripheral milling topography simulation. The results show that the maximum surface roughness error is 7.56% without considering the fixturing layout vibration, whereas the maximum contact force error is no more than 12% with fixturing layout vibration. Compared with calculations without clamping vibration considered, the calculation accuracy of contact force and surface roughness considering clamping vibration are improved by 3% and 4.72% (up milling)/3.00% (down milling), respectively.

Key words

milling / workpiece position deviation / milling vibration / surface topography

QR code of this article

Cite this article

Download Citations
LIU Yang , QIN Guohua , WU Zhuxi , LOU Weida , LAI Xiaochun. Analysis and Prediction of Surface Topographyin Peripheral Milling Based on Workpiece Vibration and Milling-Tool Structure. Acta Armamentarii. 2023, 44(7): 2132-2146 https://doi.org/10.12382/bgxb.2022.0289

0 Introduction

The dynamic cutting process is bound to cause vibration of the workpiece, which affects the surface morphology of the workpiece, and then affects the physical and mechanical properties and practical life of the parts.
Based on the free-form surface equation, Zhao Jun et al calculated the kinematic envelope surface of the tool edge line, and modeled and simulated the surface topography machined by a ball-end milling cutter. The results show that the selection of the swing angle of the tool axis has a great influence on the surface topography when the workpiece surface is a plane[1]. When the machined surface is a free-form surface, the selection of the swing angle of the cutter shaft has little effect on the surface topography. Liu et al. Developed a simulation system for milling the machined surface, which makes it possible to dynamically track the tool-workpiece interaction with the tool motion. By modifying the tool center motion and the shape of the cutting edge, the effects of tool runout and wear are included in the developed model, respectively, so that the microscopic characteristics of the machined surface can be described more accurately[2]. Lavernhe et al. Established a prediction model of surface topography by using N-buffer method, and the analysis of the prediction model showed that the tool inclination angle had a great influence on the three-dimensional topography[3]. On the basis of dimensional analysis and finite element analysis, Fernando et al. Carried out numerical simulation on the surface topography of AISI H13 steel ball end milling, and verified the effectiveness of the model by milling experiments[4]. Schmitz et al. Studied the influence of tool installation error caused by tool position offset and tool tilt on the surface topography of the machined workpiece in the face milling process, and verified the relationship between runout, surface roughness, stability and surface positioning error through simulation[5]. Chen et al. Established a surface generation model based on the geometric features of a micro-milling cutter, and studied the effects of axial runout and radial runout on the surface topography[6]. Zhang et al. Proposed a new surface topography model by using the dynamic characteristic function obtained by the receiver coupling method and combining with the milling force model to calculate the dynamic deflection of the tool.By introducing the tool deflection into the micro-flank milling edge trajectory equation, it is concluded that the most important factors affecting the surface topography in micro-flank milling are the spindle speed, feed rate, radial depth of cut, and milling method[7]. Arizmendi et al. Proposed a new method to predict the three-dimensional topography of face milling surface, and the results showed that the axial runout between face milling cutter inserts significantly affected the surface topography[8]. Wang et al. Considered the influence of tool parallel axis offset and position angle, changed the position of the intersection point between the ellipse center and the cutting edge, and analyzed the factors affecting the development trend of roughness by establishing a surface topography prediction model[9]. The above research work is mainly based on the cutting edge geometry of the tool, and further considers the influence of the tool runout and deflection on the machined surface topography.
Zeroudi et al. Established the cutting force prediction model of the tool in free milling, and established the reconstruction model of the surface topography of the milled workpiece by considering the influence of the cutting force on the surface topography of the machined workpiece[10]. Chen et al. Proposed a numerical model to predict the surface topography of thin-walled parts after elastic deformation, established the cutting force model of thin-walled parts according to the geometric characteristics in the cutting process, and studied the meshing relationship between the tool and the workpiece[11]. Peng et al. Established a plastic deformation model considering tool tilt, chatter, runout, cutting force and material, and proved that the parameters recommended by the model can well control the three-dimensional surface topography[12]. Song Wengang et al. Used the arithmetic average height and the maximum height of the surface to characterize the surface topography, and established a cutting layer model considering the initial surface geometry of the workpiece.The cutting layer model is verified by two-step milling experiments, and the results show that the initial surface geometry of the workpiece is positively correlated with the milling force and surface roughness parameters of the second step milling[13]. Taking the tool tip as the research object, Bai Lijuan et al. Used the simulated tool tip trajectory obtained by vibration-assisted milling to generate biological surface topography, and compared with the experimental results, it showed that vibration-assisted milling surface topography was feasible in the field of bionic surface processing[14]. Lu et al. Established a surface topography simulation model to predict the surface roughness based on the actual cutting trajectory and flexible deformation of the micro-milling cutter, and verified the accuracy of the model through experiments[15]. Liang Xinguang et al. Used the three-dimensional trochoid trajectory of the ball-end milling cutter blade to read the scallop height of the workpiece surface and the vibration response of the tool at equal intervals and for a long time, and realized the nominal topography modeling of the machined surface and the dynamic response trajectory simulation of the system[16]. Zhang Guohua et al. Used metal cutting theory to analyze the elliptical vibration cutting process and three-dimensional cutting model, and the results showed that the different phase difference eigenvalues between two adjacent revolutions in the cutting process had an important impact on the final surface morphology[17]. Costes et al. Used a non-contact displacement sensor to measure the high-frequency vibration signal of the tool in the milling process, and established the mechanical model of tool deformation by extracting the step length and the tool center point, so as to obtain the surface topography characteristics of the workpiece[18].
Although the above research work considers the changes of milling topography caused by tool deformation, vibration and other factors, it does not study the formation law of workpiece surface topography from the perspective of fixture layout.
In this paper, the kinematics model of the fixture layout is established according to the actual condition of one-way contact between the workpiece and the fixture, and the vibration differential equation of the workpiece in the milling process is derived by using the energy method. Then the decoupling method of vibration differential equation is proposed to calculate the deviation of the workpiece position. Finally, the motion trajectory equation of the cutting edge in the peripheral milling process is established, and the simulation method of the workpiece surface topography is proposed combined with the workpiece position deviation.

1 Vibration differential equation of fixture layout

In the machining process, the rigid workpiece will be subjected to external loads such as gravity Fg, machining force Fc, and torque Mc. In order to resist the destructive effect of these external loads, a clamping force must be applied to the workpiece.
Suppose that after the workpiece is positioned by m positioning elements, n clamping forces Fj(1≤j≤n) are applied. The work-locating element is equivalent to a semi-elastic contact model, then the ith locating element can be regarded as a damping spring with neglected mass in the normal direction ni=[nxi, nyi, nzi]T and two tangential directions ti=[txi, tyi, tzi]T, ai=[axi, ayi, azi]T,The stiffness and damping coefficients are kni, kti,kai and cni, cti, cai. Record rw=[xw, yw, zw]T, Θw=[αw, βw, γw]T as the position and direction of the workpiece coordinate system Owxwywzw relative to the global coordinate system Oxyz, respectively, as shown in Fig. 1. Assuming that rw=[xw, yw, zw]T is the coordinate of any point P on the workpiece in the Owxwywzw coordinate system, the attitude of point P at any time in the Oxyz coordinate system is r=rw+T(Θw)rw, where
T(Θw)= cos βwcos γw-cos αwsin γw+sin αwsin βwcos γwsin αwsin γw+cos αwsin βwcos γwcos βwsin γwcosαw cos γw+sin αwsin βw sin γw-sin αwcos γw+cos αwsiin βwsin γw-sin βwsin αwcos βwcos αwcos βw
Fig.1 Workpiece vibration model

图1 工件振动模型

Full size|PPT slide

Therefore, the vibration differential equation of the workpiece in the fixture layout system can be obtained, that is[19],
Mx··+ Cx·+Kx=F(t)
(1)
Where: x· = dqwdt and x·· = d2qwd2t represent the velocity and acceleration of the workpiece, respectively, and dqw represents the position offset of the workpiece; M = , M = m   m   m, m denotes the mass of the workpiece, J = Jx-Jxy-Jxz-JxyJy-Jyz-Jxz-JyzJz,Jx= (y2+z2)dm 、Jy= (x2+z2)dm 、Jz= U UN y2+x2)dm moment,Jxy= U Nxydm,Jyz= K m J6 yzdm,Jzx= U xdm;K= i=1m UKi,Ki= N krikrθikTrθikθi+kθθiWhere the tension-compression stiffness kri= U kxi   kyi   kzi U N ni,ti,aikni   kti   kaiTorsional stiffness kθi= U kαi   kβi   kγiFlexural rigidity krθi= U 0kxβi-kxγi-kyαi0kyγikzαi-kzβi0 U N kθθi= K 0-kαβi-kγαi-kαβi0-kβγi-kγαi-kβγi0And kαi=kyi U zi2 N +kzi K yi2kβi=kzi U xi2 N+kxi K zi2kγi=kxi U yi2 N+kyi K xi2kxβi=kxizi,kxγi=kxiyi,kyγi=kyixi,kyαi=kyizi,kzαi=kziyi,kzβi=kzixi,kβγi=kxiyizi,kγαi=kyizixi,kαβi=kziyixi;C= U i=1m NCi,Ci= K CriCrθiCTrθiCθi+CθθiWhere the tension-compression damping Cri= U Cxi   Cyi   CziTorsional damping Cθi= U Cαi   Cβi   CγiBending damping Crθi= U 0Cxβi-Cxγi-Cyαi0CyγiCzαi-Czβi0 UN Cθθi= K 0-Cαβi-Cγαi-Cαβi0-Cβγi-Cγαi-Cβγi0And Cαi=Ciy U zi2 N +Ciz K yi2Cβi=Ciz U xi2 N+Cix K zi2Cγi=Cix U yi2 N+Ciy K xi2 K6Cxβi=Cixzi,Cxγi=Cixyi,Cyγi=Ciyxi,Cyαi=Ciyzi,Czαi=Cizyi,Czβi=Cizxi,Cβγi=Cixyizi,Cγαi=Ciyzixi,Cαβi=Cizyixi
It is worth noting that the purpose of workpiece clamping is to apply a clamping force to the workpiece to ensure that the position of the workpiece relative to the tool is unchanged during positioning. To this end, the support reaction of the ith positioning element in equation (1) shall satisfy
nTikridri≥0, i=1,2,…,m
(2)
To describe equation (2) more simply, it is further expressed in matrix form as

kNψdqw≥0

(3)
Where :kN= kn,,knm, kn=[kni,0,0];ψ= I3×3I3×3I3×3I3×3Ω1Ω2ΩiΩmT,I3×3 is the identity matrix and Ωi= 0zi-yi-zi0xiyi-xi0.

2 Decoupling and Solution of Vibration Equation

In some cases, M, C, and K are all diagonal matrices, so it is easy to solve equation (1), but it is not universal. In the actual machining process, the differential equation of the vibration system will appear coupling terms. Therefore, the decoupling method of vibration equation is the key to solve the position deviation of the workpiece.

2.1 Orthogonality of mode shape vectors

From equation (1), the expression of the free vibration response without damping can be expressed as
Mx··+Kx=0
(4)
In the free vibration of the system, assuming that all the masses move harmonically, the solution of the equation can be expressed as

xs=A(s)sin(ωst+φs)

(5)
Where the :xs represents the array of displacements of the s-th mode shape; The A(s) represents the maximum displacement or amplitude vector of each point in the s-th mode shape; ωs represents the natural frequency of the s-th mode shape; φs is the phase angle of the s-th mode shape.
Substituting formula (5) into formula (4), we can get
(K- ωs2M)A(s)=0
(6)
The following two equations for the principal modes A(s) and A(f) corresponding to the natural frequencies ωs and ωs are obtained:
KA(s)=ωs2MA(s)
(7)
KA(f)=ωf2MA(f)
(8)
Left multiplying equation (7) by (A(f))T, then right multiplying equation (8) by A(s), we can get
(A(f))TKA(s)=ωs2(A(f))TMA(s)
(9)
(A(f))TKA(s)=ωf2(A(f))TMA(s)
(10)
In this way, equation (10) is subtracted from equation (9) to obtain
(ωs2- ωf2)(A(f))TMA(s)=0
(11)
Dividing both sides of equation (9) by ωs2 and subtracting both sides of equation (10) by ωf2 gives
(1ωs2- 1ωf2)(A(f))TKA(s)=0
(12)
When s s f and the eigenvalues are ωs≠ωf, to satisfy equations (11) (12), the following relations must hold:

(A(f))TMA(s)=(A(s))TMA(f)=0

(13)

(A(f))TKA(s)=(A(s))TKA(f)=0

(14)
Equations (13) and (14) show that there is orthogonality between two principal modes with unequal natural frequencies with respect to the mass matrix M and the stiffness matrix K.
When s = f, equations (11) and (12) are true for any value, so that

(A(s))TMA(s)=Mps

(15)

(A(s))TKA(s)=Kps

(16)
Where :Mps and Kps are the principal mass and principal stiffness of the s-th order, respectively, both constants, and they depend on how the eigenvector A(s) is normalized.
Normalize the A(s) and arrange them into columns in order to obtain the modal matrix Ap, which can be expressed in the following form:
Ap=[A(1),…,A(s),…,A(λ)]=A1(1)A1(s)A1(λ)A2(1)A2(s)A2(λ)A6(1)A6(s)A6(λ)
(17)
Where: λ is the maximum mode order.
Substituting formula (17) into formula (15) and formula (16), we can obtain
ATpMAp=Mp
(18)
ATpKAp=Kp
(19)
Where :Mp and Kp are the principal mass matrix and principal stiffness matrix, respectively, both of which are diagonal matrices.

2.2 The decoupling of Eq.

The properties of the coupled equations do not depend on the inherent characteristics of the system, but on the choice of generalized coordinates, so the principal coordinates Xp of the decoupled equations must be found. Left multiplying equation (4) by ATp and inserting I=Ap Ap-1 in front of x· and x··, we have
ATpMAp Ap-1x··+ ATpKAp Ap-1x=0
(20)
Substituting equations (18) and (19) into equation (20) by using the orthogonality relationship, we can obtain
Mp X··p+KpXp=0
(21)
Where:
Xp=Ap-1x
(22)

X··p= Ap-1x··

(23)
Since the principal mass matrix Mp and the principal stiffness matrix Kp are diagonal matrices, there is no coupling in the differential equations of motion of the system described by using principal coordinates.
For the convenience of operation, each principal mode is regularized, that is, a group of specific principal modes, that is, regular modes, are taken and expressed by the column matrix AN(s), s = 1,2, …, λ, so there is
(AN(s))TM AN(s)=1
(24)
The AN(s) can be found from the A(s) of the principal modes.
AN(s)=1μsA(s)
(25)
Where :μs is an undetermined constant called the regularization factor.
Substituting equation (25) into equation (24), we obtain
μs=Mps=(A(s))TMA(s)
(26)
Solving the equation (25) derived from the μs, the regular mode shape matrix AN is obtained as
AN=[A(1),…,A(s),…,A(λ)1μ1          1μi          1μλ=AN1(1)AN1(s)AN1(λ)AN2(1)AN2(s)AN2(λ)AN6(1)AN6(s)AN6(λ)
(27)
The canonical mass matrix MN, calculated according to equation (18) with AN, is an identity matrix I, that is
ATNMAN=MN=I
(28)
It can be seen from formulas (7) and (8) that
ATKA=ATMA ωs2
(29)
Substituting the regular mode shape array AN(s) into equation (29), we can obtain
ωs2=(AN(s))TKAN(s)(AN(s))TMAN(s)=KNs1=KNs
(30)
The canonical stiffness KNs is equal to the square of the natural frequency, i.e., ωs2. Therefore, using the canonical mode AN, the diagonal elements of the canonical stiffness KN calculated according to Equation (19) are the square values of the natural frequencies of each order, respectively, that is,
KN=KN1000KN20000KNs=ω12000ω220000ωs2
(31)
For an equation of motion with damping C, the canonical damping matrix CN can be expressed as
CN=ATNCAN=CN11CN12CN16CN21CN22CN26CN61CN62CN66
(32)
In practical engineering applications, the CN is not a diagonal matrix in most cases, but the damping of most vibration systems in engineering is small, and because of the complexity of various damping, the simplest way to diagonalize the regular mode matrix is to change the value of its off-diagonal elements to 0, that is,
CNCN11000CN22000CN66
(33)
Where :CNss is the damping coefficient of the s-th regular mode.
In practical vibration analysis, the damping ratio ξs of each mode is usually given by the damping coefficient CNss by experiment or measurement. Therefore, combined with equation (32), we can get
CNCN11000CN22000CNnn=2ξ1ω10002ξ2ω20002ξsωs
(34)

2.3 Equation solving

Let AN be the regular mode shape matrix of equation (1), and let Xp be the principal coordinates of the vibration system. Then we have ATN MAN=I6×6, ATN KAN=Γ=diag [ω12, ω22 ,…, ω62 ], ATN CAN=Z=diag[2ξ1ω1, 2ξ2ω2,…,2ξ6ω6] 。 The Xp= AN-1 δqw can be obtained by linear transformation of the principal coordinates. Denote G (t) = ATN F (t), then equation (1) can be further transformed into
X··p+ ZX·p+ΓXp=G(t)
(35)
The decoupled form of (35) is
X··ps+2ξsωs X·ps+ ωs2Xps=Gs(t)
(36)
Where: 1 ≤ s ≤ 6.
Each principal coordinate can be solved separately by using the Duhamel integral, that is
Xps(t)=1ωs1-ξs20t e-ξsωs(t-η)·sin ωs 1-ξs2(t-η)Gs(η)dη
(37)
Where: t is the system time; η is any time in the system time.
Therefore, the solution of the vibration differential equation is

dqw=ANXp
s.t.
kNψdqw≥0

(38)

3 Blade Motion Trajectory

In addition to the dynamic clamping layout of the workpiece, another major factor affecting the machined surface topography of the workpiece is the geometry of the cutting edge of the tool. The key point of modeling is to solve the motion trajectory of the blade in space.
The peripheral milling process is shown in Fig. 2. Oxyz and Owxwywzw are the global coordinate system and workpiece coordinate system respectively, while Ocxcyczc and Osxsyszs are the tool coordinate system and spindle coordinate system respectively. In general, for convenience, the coordinate system should be established so that Oxyz coincides with the Owxwywzw; R is the milling cutter radius, γ is the helix angle, e is the cutter eccentricity, α is the rotation angle, vf is the tool feed rate, φj is the initial angle of the jth cutter tooth, w is the tool rotation speed, w is positive for down milling and negative for up milling, P is the cutting point, and ABCD is the cutting area.
Fig.2 Diagram of peripheral milling cutting edge movement

图2 周铣切削刃运动示意图

Full size|PPT slide

If the cutter has Z cutter teeth, the initial angle of the jth cutter tooth can be expressed as
φj=2π(j-1)z
(39)
In this way, the coordinates of any point P on the jth cutting edge of the milling cutter in the Ocxcyczc coordinate system should be

rcP= xcPycPzcP= Rcos(φj+α)Rsin(φj+α)Rα/tanγ

(40)
Assuming that the xs axis of the Osxsyszs coordinate system coincides with the xc axis of the Ocxcyczc coordinate system,Then the position of the Ocxcyczc coordinate system with respect to the Osxsyszs at any time t,The directions are rsc =[e,0,0]T, Θsc =[wt,wt,0]T, respectively. Then the coordinates of point P in the Osxsyszs coordinate system should be
rsP=rsc+T(Θsc) rcP
(41)
Moreover, at any time instant t, the spindle has moved by vft in the feed direction y. If Osxsyszs coincides with Oxyz at the initial time, then at time t, the position and direction of Osxsyszs relative to Oxyz coordinate system are respectively rs=[0, vft, 0]Ts= [0,0,0]T, then the coordinate of point P on the jth blade in Oxyz coordinate system is
rP=rs+T(Θs) rsP
(42)

4 Workpiece surface topography

Taking the derivative of the position R of any point on the workpiece, we know that at time t, if the position of the workpiece changes dqw due to the vibration of the fixture layout, the position change of any point on the workpiece should be

dr=Edqw

(43)
Where:
E={[1000zy010z0x001yx0],三维[10y01x],二维.
(44)
At time t, if the point P on the cutter tooth is in the cutting area ABCD along the back engagement direction, then if and only if the coordinate of the point P is between the theoretical cutting plane and the plane to be cut, it indicates that the point P is involved in the cutting, as shown in Figure 3. ap in Fig. 3 is the axial depth of cut.
Fig.3 Diagram of surface formation in the coordinate plane Oxy

图3 坐标平面Oxy内的表面成形示意图

Full size|PPT slide

According to equations (42) and (43), the coordinate rP=[xP, yP, zP] of point P on the cutting plane should be
rP+EPdqw(t)=rs+T(Θs)(rsc+T(Θsc) rcP)
(45)
Where:
EP=1000zP-yP010-zP0xP001yP-xP0, 10yP01-xP, 
(46)
Obviously, equation (45) is difficult to solve, so the discrete method is used to solve it. The steps of the surface topography discrete algorithm are as follows:
Step 1 Discretize the plane to be cut.
The workpiece is uniformly divided into m × n grids with step sizes Δy and Δz, respectively, and the coordinate values of each grid node are stored in a matrix as X (u, e), y (u, e), and Z (U, E), respectively, where u = 1,2, …, m + 1 and e = 1,2, …, n + 1.
Step 2 Discrete cutting time and cutting edge.
The time discrete step Δt and the number of discrete points of the cutting edge l are set according to the mesh accuracy of the workpiece.Generally, Δt=min{Δy, Δz}/(wR),l=hT/min{Δy, Δz},hT is taken as the effective length of the cutting part of the milling cutter to ensure that the projection of the discrete infinitesimal element of the cutting edge in the machining plane does not exceed the mesh spacing of the workpiece and sweeps through at most one mesh point of the workpiece in a unit time step.
Step 3 Initialize the cutting time t = 0s.
Step 4 calculates the fixture layout vibration at the current time.
The coordinate r(u,e)=[x(u,e), y(u,e), z(u,e)]T of each mesh node on the cutting plane is obtained according to equation (43).
Step 5 Initialize the cutting edge number J = 1.
Step 6 Initialize the discrete point number K = 1.
Step 7: Determine whether the discrete point is in the processing area?
According to equation (42), calculate the coordinate value (xk, yk, zk) of the discrete point Pk on the cutting edge at the current time, and judge whether it is in the plane to be machined. If yes, go to step 8, otherwise go to step 10.
Step 8 searches for the plane mesh node to be machined nearest to the discrete point of the cutting edge. A discrete point corresponds to a grid point when its position (yk, zk) coincides with or is closest to the position of the grid point (y (u, e), Z (u, e)).
Step 9 Judge the cutting condition of the discrete points of the cutting edge. Compare the relationship between the nearest grid node X (u, e) and the xk of the discrete point Pk of the cutting edge, if xk<x(u,e), it indicates that the tool has cut into the workpiece, and update the stored value of X (u, e) with the value of xk; Otherwise, no treatment will be made.
Step 10 Determine if it is the last discrete point? If yes, go to step 11; Otherwise, calculate the next discrete point K = K + 1, and go to step 7.
Step 11 Determine if it is the last cutting edge? If yes, go to step 12; Otherwise, calculate the next cutting edge J = J + 1, go to step 6.
Step 12 determines whether it is the last moment? If so, the calculation process is ended, and the stored X (u, e), y (u, e) and Z (u, E) form the cutting surface topography; Otherwise, calculate the next moment and go to step 4.

5 Application example

Two typical examples are used to illustrate the calculation of surface topography in detail, and the correctness of the calculation results of surface topography is verified. In the first example, the workpiece is fastened on the worktable, the influence of clamping vibration is ignored, and the surface roughness value is directly measured for verification. The second example is from reference [20], which indirectly verifies the validity of the calculation process by verifying the magnitude of the contact force, and gives the surface topography considering the clamping vibration.

5.1 Surface topography without consideration of vibration

Fig. 4 is a field diagram and an experimental result diagram of an end mill peripheral milling workpiece. The material of the workpiece is aluminum alloy, and the Young's modulus is Ew=207GPa; The material of the tool is cemented carbide, the Young's modulus is ET=530GPa, and other parameters are shown in Table 1. According to the clamping scheme and tool parameters, the rigidity of the tool and the workpiece in the peripheral milling process is strong, so the vibration of the tool and the workpiece in the machining process is ignored.
Fig.4 Peripheral milling experiment

图4 周铣实验

Full size|PPT slide

Table 1 Peripheral milling parameters

表1 周铣加工参数

刀具与工艺 参数 数值
齿数z 2
刀具 半径R/mm 5
螺旋角β/ (°) 45
刃长L/mm 10
工艺 主轴转速ω/(r·min-1) 500
进给速度vf/(mm·min-1) 600
The effectiveness of the peripheral milling algorithm is verified separately without considering the influence of the workpiece fixture layout on the milled surface topography of the workpiece in the machining process, that is, dqw=0.
Set the milling parameters according to Table 1. Let Δy = Δz = 0.05 mm, and divide the milled surface into 120 × 120 meshes uniformly. The coordinate values of each grid node are stored with matrices X (u, e), y (u, e), Z (u, E). After that, the discrete time interval Δt=min{Δy, Δz}/(wR)=9.55×10-5s and the discrete number of cutting edge points l=hT/min{Δy, Δz}=120 are calculated.
Since the vibration of the workpiece is not considered, the discrete algorithm of the surface topography can be simplified without considering step 4 in the algorithm, and the result of the surface topography after calculation is shown in Figure 5.
Fig.5 Surface morphology without considering vibration

图5 不考虑振动的表面形貌

Full size|PPT slide

Fig. 6 gives the actual machined surface by up milling and climb milling in the case of steady state machining. By comparing the results of the actual machined surface and the simulation topography, it can be seen that the calculation results of the discrete algorithm of surface topography are very consistent with the actual machining situation.
Fig.6 Actual machined surface in stable machining(the left represents the actual machining table, the middle represents the black and white simulated surface, and the right represents the colored simulated surface)

图6 稳态加工中的实际加工表面(左为实际加工表面,中为黑白仿真表面,右为彩色仿真表面)

Full size|PPT slide

For further numerical comparison, the roughness of the milled surface was measured on a roughness profiler JB-6 CA manufactured by Shanghai Precision Instrument Co., Ltd., China, as shown in Fig. 7, and the measurement results are shown in Fig. 4 (B). The Ra of the arithmetic mean deviation of the profile of the measured and calculated results are listed in Table 2. By comparison, the results obtained by the discrete algorithm of up milling surface topography have better accuracy.
Fig.7 Roughness measurement

图7 粗糙度测量

Full size|PPT slide

Table 2 Experimental verification of surface morphology

表2 周铣形貌仿真实验验证结果

铣削方式 实测粗糙
/μm
仿真粗糙
/μm
相对
误差/%
逆铣 2.428 2.294 5.52
顺铣 3.057 2.826 7.56

5.2 Surface topography considering vibration

A milling cutter with radius R = 9 mm was used to mill a groove on the workpiece surface with axial depth of cut ap=3mm, feed per tooth f = 0.05 mm, and milling speed V = 100 mm/min, as shown in Fig. 8[19-20]. In Fig. 8, F1, F2 are the clamping force, and L0~L5 is the positioning element. The size of the workpiece is 220 mm × 122 mm × 112 mm, the material is 7075-T6 aluminum alloy, and the Young's modulus and Poisson's ratio of the workpiece are Ew=70GPa and bw=0.334, respectively. The gravity of the workpiece is Fgrav=[0N, 59.73N, 0N]T, the position of the center of mass is rgrav=[61mm]T, and the machining force and torque are Fmach=[-131N, -55N, 232N]T,Mmach=[0N·m, 2.77N·m, 0N·m]T, respectively[20].
Fig.8 3-2-1 fixturing layout

图8 3-2-1装夹方案

Full size|PPT slide

To ensure the smooth realization of the machining process, the clamping forces F1= 640N and F2=670N are applied at r1=[60mm,0mm]T and r2=[60mm]T, respectively. The position and unit normal vector of each positioning element are shown in Table 3. The flat-head positioning element is made of hard steel, the Young's modulus and Poisson's ratio of the positioning element are Ef=207GPa and bf=0.292 respectively, the radius of the positioning element is Rf=6mm, the length is Hf=12mm, and the damping coefficient is cf=7000N·s/m.
Table 3 Position and normal vector of locators

表3 定位元件的位置与法向量

According to the data of workpiece mass and elastic modulus of positioning element, the stiffness matrix, mass matrix and damping matrix are easily obtained as follows:
K=1.2042×106000-1.2042×106-4.8169×10603.6127×1090-5.6598×10700002.4084×1099.6338×106000-5.6598×1079.6338×1068.9534×10600-1.2042×1060002.4086×1074.8169×105-4.8169×1060004.8169×1052.4104×107
(47)
M=6.09490000006.09490000006.09490000000.01390000000.03210000000.0310
(48)
C=7000000-7-280210000-32.900001400056000-3295652.04500-7000140.0070.028-280000.028140.112
(49)
Combining the mass, stiffness and damping matrices and equation (1) and decoupling, six second-order differential equations can be obtained, where the dynamic external load is F(t)=[801N, 114.73N, -872N, -9.1024N·m, 22.08-0.3867tN·m, -0.0917t-0.3242N·m]T, and the cutting time is 0s ≤ t ≤ 132s. According to equation (37), the change of the position offset of the workpiece with time t can be obtained, and the results are shown in Figure 9 and Figure 10.
By comparison, it can be found that the position deviation trend of the workpiece is roughly the same in the two cases of considering damping and not considering damping. However, when the damping is not considered, the change of the workpiece position deviation is unstable. When damping is considered, the workpiece position deviation is stable.
Fig.9 Workpiece position error without damping

图9 不考虑阻尼的工件位置变化

Full size|PPT slide

Fig.10 Workpiece position error with damping

图10 考虑阻尼的工件位置变化

Full size|PPT slide

Because the vibration deviation of the workpiece position is caused by the position deviation of each positioning element, the calculated value of the workpiece vibration can be indirectly guaranteed by measuring the supporting reaction force at the positioning element. A diagram of the experimental setup is shown in Fig. 11[20]. The experimental measured value of the support reaction force at each positioning element is shown in Fig. 12 (a), the calculated value of the static model is shown in Fig. 12 (B), and the dynamic support reaction force calculated according to formula (3) is shown in Figure 12 (C), in which the numbers 0 to 5 represent the L0~L5 of the positioning element respectively[20][21]. It can be seen from the analysis and comparison in Fig. 12 that the error of the calculated value of the static model of the support reaction at the positioning element is 15%, and the calculated value of the dynamic model is closer to the experimental result, with an error of about 12%, and the accuracy is improved by 3%[20]. Obviously, both the change trend of the support reaction and the magnitude of the force are very consistent.
Fig.11 Experimental device

图11 实验装置

Full size|PPT slide

Fig.12 Supporting reaction force

图12 支撑反力

Full size|PPT slide

According to the discrete algorithm of surface topography proposed in this paper, any side of the groove is selected to simulate the topography, and the stiffness G=πER2/H is calculated by the formula, where H is the length of the positioning element, and the stiffness of the positioning element is about half of that of the tool, so the influence of clamping layout vibration on the surface topography is mainly considered, and the results are shown in Figure 13Figure 14.
Fig.13 Morphology of up milling

图13 逆铣形貌

Full size|PPT slide

Fig.14 Morphology of down milling

图14 顺铣形貌

Full size|PPT slide

It can be seen from Fig. 13 and Fig. 15 that the three-dimensional topography without considering the clamping layout is roughly the same as that with considering the clamping layout. However, the overall cutting depth is deeper without considering the clamping layout, because the workpiece is rigid, and the vibration of the workpiece will cause the overall position deviation of the workpiece when considering the clamping layout.
Fig.15 Profile comparison

图15 截面轮廓对比

Full size|PPT slide

From the surface quality point of view, the simulated surface is significantly more uneven when considering the effect of fixture layout on the machined surface topography. Through the calculation, the arithmetic average profile deviation Ra of up milling and down milling without considering the fixture layout is 5. 956 μm and 6. 358 μm, respectively, and the arithmetic average profile deviation Ra of up milling and down milling with considering the fixture layout is 6. 237 μm and 6. 549 μm, respectively, and the calculation accuracy is improved by 4. 72% and 3. 00%, respectively.

6 Conclusion

The main results are as follows: 1) The contact mode between the workpiece and the fixture is equivalent to a linear spring-damper system, and the workpiece fixture layout model is established. The vibration differential equation of the workpiece is established by using the energy method, and the single contact constraint condition of the workpiece-fixture layout is derived by using Hooke's law according to the working condition that the contact force must be pressure in the machining process.
2) According to the decoupling method in mechanical vibration, the modal analysis solution method of vibration differential equation under constraint condition is proposed by using coordinate transformation method, so as to realize the prediction of workpiece position deviation.
3) Based on the transformation matrix and vector operation rule of the motion homogeneous coordinate, the motion trajectory equation of the peripheral milling cutting edge is established, and the discrete simulation algorithm of the peripheral milling surface topography is given.
4) Without considering the workpiece vibration caused by fixture layout, the topography simulation algorithm is verified by experiments, and the relative errors of up milling and down milling roughness are 5.52% and 7.56%, respectively. In the case of considering the clamping layout, the validity of the positioning element position offset is verified by the support reaction, and the error of the model is within about 12%, compared with the static model, the accuracy is improved by 3%, and the topography of the two cases is simulated and compared.
5) In the actual machining process, if the cutting part of the tool is long, the influence of tool vibration on the surface topography is more serious than that of the clamping scheme, which is also the focus of future research.

References

[1]
赵军, 方文勇, 韩式国. 自由曲面数控加工表面形貌仿真[J]. 工具技术, 2007, 41(8): 17-19.
ZHAO J, FANG W Y, HAN S G. Topography simulation of NC milling of free-form surface[J]. Tool Engineering, 2007, 41(8): 17-19. (in Chinese)
[2]
LIU X B, SOSHI M, SAHASRABUDHE A, et al. A geometrical simulation system of ball end finish milling process and its application for the prediction of surface micro features[J]. Journal of Manufacturing Science & Engineering, 2006, 128(1):74-85.
[3]
LAVERNHE S, QUINSAT Y, LARTIGUE C. Model for the prediction of 3D surface topography in 5-axis milling[J]. International Journal of Advanced Manufacturing Technology, 2010, 51(9/10/11/12): 915-924.
[4]
FERNANDO R P, SOLDANI X, LOYA J, et al. A new approach for time-space wear modeling applied to machining tool wear[J]. Wear, 2017, 390-391:125-134.
[5]
SCHMITZ T L, COUEY J, MARSH E, et al. Runout effects in milling: Surface finish, surface location error, and stability[J]. International Journal of Machine Tools & Manufacture, 2007, 47(5): 841-851.
[6]
CHEN W Q, SUN Y Z, HUO D H, et al. Modelling of the influence of tool runout on surface generation in micro milling[J]. Chinese Journal of Mechanical Engineering, 2019, 32(1): 1-9.
[7]
ZHANG X, PAN X D, WANG G L. Influence factors of surface topography in micro-side milling[J]. International Journal of Advanced Manufacturing Technology, 2019, 105(12): 5239-5245.
[8]
ARIZMENDI M, JIM NEZ A. Modelling and analysis of surface topography generated in face milling operations[J]. International Journal of Mechanical Sciences, 2019, 163: 105061-1-9.
[9]
WANG L P, GE S Y, SI H, et al. Elliptical model for surface topography prediction in five-axis flank milling[J]. Chinese Journal of Aeronautics, 2020, 33(4): 1361-1374.
[10]
ZEROUDI N, FONTAINE M. Prediction of machined surface geometry based on analytical modelling of ball-end milling[J]. Procedia CIRP, 2012, 1: 108-113.
[11]
CHEN Z T, YUE C X, LIU X L, et al. Surface topography prediction model in milling of thin-walled parts considering machining deformation[J]. Materials, 2021, 14(24): 7679.
With the continuous improvement of the performance of modern aerospace aircraft, the overall strength and lightweight control of aircraft has become a significant feature of modern aerospace parts. With the wide application of thin-walled parts, the requirements for dimensional accuracy and surface quality of workpieces are increasing. In this paper, a numerical model for predicting surface topography of thin-walled parts after elastic deformation is proposed. In view of the geometric characteristics in the cutting process, the cutting force model of thin-walled parts is established, and the meshing relationship between the tool and the workpiece is studied. In addition, the influence of workpiece deformation is considered based on the beam deformation model. Cutting force is calculated based on deformed cutting thickness, and the next cutting–meshing relationship is predicted. The model combines the radial deflection of the workpiece in the feed direction and the changing meshing relationship of the tool–workpiece to determine the three-dimensional topography of the workpiece. The error range between the experimental and the simulation results of surface roughness is 7.45–13.09%, so the simulation three-dimensional morphology has good similarity. The surface topography prediction model provides a fast solution for surface quality control in the thin-walled parts’ milling process.
[12]
PENG Z X, JIAO L, YAN P, et al. Simulation and experimental study on 3D surface topography in micro-ball-end milling[J]. International Journal of Advanced Manufacturing Technology, 2018, 96(5): 1943-1958.
[13]
宋文刚, 刘战强, 吕文军. 工件初始表面几何状态对铣削力及表面形貌的影响[J]. 工具技术, 2021, 55(7): 28-32.
SONG W G, LIU Z Q, W J. Influence of initial surface geometric states of workpiece on milling force and surface topography[J]. Tool Engineering, 2021, 55(7): 28-32. (in Chinese)
[14]
白利娟, 张建华, 陶国灿, 等. 振动辅助铣削加工仿生表面研究[J]. 中国机械工程, 2016, 27(9): 1229-1233,1242.
BAI L J, ZHANG J H, TAO G C, et al. Vibration assisted milling for bionic surface manufacturing[J]. China Mechanical Engineering, 2016, 27(9): 1229-1233,1242. (in Chinese)
[15]
LU X H, HU X C, JIA Z Y, et al. Model for the prediction of 3D surface topography and surface roughness in micro-milling Inconel 718[J]. The International Journal of Advanced Manufacturing Technology, 2018, 94(5): 2043-2056.
[16]
梁鑫光, 姚振强. 基于动力学响应的球头刀五轴铣削表面形貌仿真[J]. 机械工程学报, 2013, 49(6): 171-178.
Abstract
针对自由曲面球头刀五轴铣削中刀具与工件复杂的位姿关系,利用球头铣削刃的三维次摆线轨迹,提出一种在工件表面等距离间隔缓存残留高度、在刀具端等时间间隔缓存振动响应的双缓存离散机制,分别实现五轴铣削中的名义加工表面形貌建模与系统动态响应轨迹仿真。在此基础上,以切削刃经过工件表面残留区域的时间信息为纽带,利用系统瞬时动态位移响应的插值信息对已加工表面法向残留高度进行修正,建立考虑铣削系统动力学响应的已加工表面形貌预测模型。以变切深、变转速、恒定的刀具倾角与每齿进给量所进行的验证试验表明,提出的球头刀五轴铣削表面形貌建模方法可有效预测颤振工况与稳定铣削工况的加工表面形貌与纹理特征。
LIANG X G, YAO Z Q. Dynamic-based simulation for machined surface topography in 5-axis ball-end milling[J]. Journal of Mechanical Engineering, 2013, 49(6): 171-178. (in Chinese)
Focusing on the complexity of the relative motion between cutter and workpiece in five-axis milling, a 3D trochoid mechanism is adopted to describe the actual tooth trajectory of the ball-end cutting edge. A kind of two-buffer mechanism is proposed to calculate residual height. One is established to determine each spatial buffer point which is discretized by a constant interval on workpiece surface, the other is built in cutter centre to simulate system dynamic response for each time interval. Thus, the nominal machined surface topography and system vibration trajectory can be calculated respectively. According to the exact time when cutting edge passing by the discretized point of the residual surface, the corresponding instantaneous vibration response can be interpolated and used to update the nominal residual height so that the prediction model of surface topography can be established with the consideration of milling system dynamic response. Validation tests are carried out with the varying cutting depth, increasing spindle speed, constant lead angle and fixed feed rate. It indicates that the proposed method can be used to predict surface topography accurately and distinguish stable milling status from chatter effectively.
[17]
张国华, 李咚咚, 李茂伟, 等. 超声椭圆振动车削三维形貌形成研究[J]. 兵工学报, 2017, 38(10): 2002-2009.
Abstract
为研究椭圆振动车削表面三维几何形貌形成规律及影响因素,采用金属切削理论对椭圆振动切削过程和三维切削模型进行了分析。研究结果发现:切削过程中相邻两转之间不同的相位差特征值对最终表面形貌有着重要的影响;通过已有的试验结果表明,相位差特征值分别为0和0.5时所形成的表面形貌截然不同,与理论分析结果相吻合。
ZHANG G H, LI D D, LI M W, et al. Research on formation of 3-D micro-surface in ultrasonic elliptical vibration cutting[J]. Acta Armamentarii, 2017, 38(10): 2002-2009. (in Chinese)
[18]
COSTES J P, MOREAU V. Surface roughness prediction in milling based on tool displacements[J]. Journal of Manufacturing Processes, 2011, 13(2): 133-140.
[19]
秦国华, 何志芬, 王华敏, 等. 基于动力学与遗传算法的工件位置偏离预测与控制方法[J]. 机械工程学报, 2017, 53(1): 110-120.
Abstract
装夹布局的合理规划是影响工件加工质量与生产安全的关键因素,为此系统地提出了工件位置偏离的动力学预测模型及其控制方法。首先通过将工件与夹具之间的接触模式等价地转化为线性弹簧-阻尼系统,利用运动学和弹性力学理论推导出工件的振动微分方程,以确定工件在装夹布局中的位置偏离。根据工件与夹具之间的接触为单向这一实际情况,运用胡克定律建立了工件位置偏离的约束条件。通过坐标转换方法,提出在约束条件下振动微分方程的模态分析求解方法,实现工件位置偏离的预测,试验结果表明所提出预测结果与试验结果完全吻合。其次,为使工件位置偏离达到最小,构建工件位置偏离的最小2-范数为目标的装夹布局优化模型,通过以合理装夹布局中工件位置偏离的2-范数为自变量构造个体适应度,提出装夹布局优化模型的遗传算法求解技术。提出的工件位置偏离预测与控制方法,能够避免工件处于非稳定状态下优化模型的求解过程,为复杂工件装夹布局方案的合理设计提供了基础理论。
QIN G H, HE Z F, WANG H M, et al. Prediction and control of workpiece position deviation based on dynamics and genetic algorithm[J]. Journal of Mechanical Engineering, 2017, 53(1): 110-120. (in Chinese)
The reasonable planning of the fixturing layout is crucial to guarantee the machining quality and safety manufacturing. Therefore, the prediction and control method is systematically proposed for the workpiece position error. Above all, by transforming the contacts between the workpiece and fixture into the linear spring-damp system, the vibration differential equation of the workpiece, which can determine the position movement of the workpiece in the fixturing layout, is concluded to obtain the workpiece position deviation according to the theory of kinematics and elastic mechanics. Based on the fact that the contacts between the workpiece and fixture are unidirectional, the hook law is used to formulate the constraint conditions of the workpiece position deviation. According to coordinate transformation method, the modal analysis solution technology of the vibration differential equation is proposed to predict the workpiece position error. A numerical test is demonstrated to validate the prediction results are in good agreement with the experimental data in the reference. And then, in order to minimize the workpiece position deviation, the optimal model with the objective of minimizing the 2-norm of the workpiece position deviation is proposed for the fixturing layout. According to the workpiece position deviation in the reasonable fixturing layout, the individual fitness is defined by taking its 2-norm as independent variable. Thus, the genetic algorithm can be skillfully developed to solve the optimal model of fixturing layout. The presented prediction and control method of workpiece position error can avoid the solution of the optimal model under the condition of the workpiece in unstable fixturing. It can provide a basic theory of fixturing layout design for the complex workpiece.<p>&nbsp;</p>
[20]
TAO Z J, KUMAR A S, NEE A Y C, et al. Modelling and experimental investigation of a sensor-integrated workpiece-fixture system[J]. International Journal of Computer Applications in Technology, 1997, 10(3/4): 236-250.
[21]
QIN G H, ZHANG W H, WU Z X, et al. Systematic modeling of workpiece-fixture geometric default and compliance for the prediction of workpiece machining error[J]. ASME Journal of Manufacturing Science and Engineering, 2007, 129: 789-801.
Control of workpiece machining error (WME) is a key concern in the design of a fixture system. In this paper, source errors, which are categorized into workpiece-fixture geometric default and workpiece-fixture compliance, are systematically investigated to reveal their effects upon the WME. The underlying mechanism is that source errors lead to the workpiece position error (WPE), the workpiece elastic deformations (WED), and the inconsistent datum error (IDE), and all of them will contribute together to the WME. Here, the IDE refers to the dimension deviation of the processing datum from the locating datum once two references do not coincide. An overall quantitative formulation is proposed for the computing of WME in terms of WPE, WED, and IDE for the first time. In detail, the WPE raised in the workpiece-locating and clamping process is evaluated based on the geometric defaults and local deformations of workpiece-fixture in the contact region. The WED relative to the workpiece-clamping process is determined by solving a nonlinear mathematical programming problem of minimizing the total complementary energy of the frictional workpiece-fixture system. Some numerical tests are finally demonstrated to validate the proposed approach on the basis of both theoretical and experimental data given in the references.

62

Accesses

0

Citation

Detail

Sections
Recommended

/