Milling stability region, surface roughness, position error and texture synchronization prediction method

Abstract

This invention discloses a method for synchronously predicting milling stability domain, surface roughness, position error, and texture. It establishes the dynamic equations of a five-axis milling system, constructs a discrete mapping relationship between two adjacent tool rotation cycles, and calculates the state transition matrix. Based on Floquet theory, it obtains the stability domain for five-axis milling. According to the fixed-point theorem, it transforms the modal displacements to physical space to obtain the relative vibration displacement between the tool and the workpiece, extracting the vibration displacements along the workpiece surface normal and along the tool feed direction. It selects some cutting edge trajectory points near the workpiece's forming surface to form a set of interpolated and densified cutting edge trajectory points. The densified point set is divided along the tool axial height to form the final workpiece surface topography point cloud, from which the surface roughness can be calculated. This invention enables high-efficiency and high-precision five-axis milling of weakly rigid complex components.

Description

Technical Field

[0001] This invention belongs to the field of advanced manufacturing, specifically relating to a method for synchronously predicting milling stability domain, surface roughness, position error, and texture. Background Technology

[0002] With the rapid development of manufacturing industries such as aerospace, machinery, and defense machinery, the performance requirements for some equipment are also increasing. Especially in the aerospace field, weakly rigid complex structural components are widely used, such as aircraft spars, wing ribs, bladed disks, missile bulkheads, and engine casings. These parts are characterized by thin walls, low stiffness, and complex structures. Furthermore, during five-axis milling, the movement of the cutting tool changes frequently, which seriously affects their machinability, resulting in low machining accuracy and poor surface quality. Therefore, achieving high-efficiency, high-precision, and high-reliability machining of such weakly rigid complex structural components and improving the stability of five-axis CNC machining plays a vital role in the modernization of my country's high-end precision manufacturing technology and national defense construction. To ensure that these weakly rigid components meet performance, accuracy, and reliability requirements, strict control of chatter and various errors during milling is essential. Therefore, in-depth research into the dynamics of five-axis CNC milling and early prediction of the milling chatter stability domain and the surface quality of machined parts are crucial.

[0003] Currently, Chinese patent CN107239603B proposes a method for modeling the leaf-lobe diagram of ball end mill chatter stability domain based on fine integration in five-axis CNC machine tool machining. This method performs high-precision time-domain numerical solutions on the dynamic system of the milling system to construct the leaf-lobe diagram of ball end mill chatter stability domain during five-axis CNC machine tool machining, but it does not consider the case of flexible workpieces. Chinese patent application CN113609611A proposes a method for predicting the surface position error of weakly rigid parts based on element nonlinearity technology. It utilizes finite element nonlinearity technology to establish a workpiece dynamic model considering the cutting process, achieving prediction of workpiece surface position error during stable cutting. This patented method only applies to three-axis milling and can determine in advance whether the selected machining conditions can meet the surface quality requirements. In the doctoral dissertation "Yuan Lei. Surface Quality Analysis and Machining Process Parameter Optimization of Thin-Walled Parts [D]. Wuhan: Wuhan University of Technology, 2016," the moving frame method of differential geometry theory was used to establish a mathematical model of the cutting edge point under the moving frame of machining feature points. For the milling of thin-walled parts, a zero-order frequency domain method for simultaneous prediction of surface morphology and surface position error was proposed. Furthermore, the milling force convolution model was combined with the Floquet-Nyquist method to obtain an analytical expression for stability analysis, which was applied to the stability prediction of the milling system. However, it did not achieve simultaneous prediction of multiple objectives such as stability, surface morphology, and surface position error. Summary of the Invention

[0004] To address the aforementioned problems, this invention proposes a method for simultaneous prediction of milling stability domain, surface roughness, position error, and texture. This method combines the prediction of chatter stability domain, which reflects machining stability constraints, with the overall surface position error, surface roughness, and surface texture, which reflect machining surface quality constraints. Simultaneously considering chatter-free operation, surface position error, surface roughness, and surface texture constraints is key to optimizing the selection of machining process parameters to achieve high-efficiency and high-precision five-axis milling of weakly rigid complex components.

[0005] To achieve the above objectives, the present invention adopts the following technical solution:

[0006] A method for synchronously predicting milling stability region, surface roughness, position error, and texture, the method specifically includes the following steps:

[0007] Step 1: Construct a dynamic model considering the five-axis milling system. Its dynamic equations can be expressed as follows:

[0008]

[0009] Where M, C, and K are the modal mass matrix, modal damping matrix, and modal stiffness matrix of the system, respectively, and q(t) represents the modal coordinates of the milling system; a p The axial cutting depth is represented by T; T represents the time delay, which is also the milling cycle, i.e., T = 60 / (ΩN). t N is the number of teeth on the milling cutter, and Ω is the spindle speed in rpm. This is a five-axis cutting force coefficient matrix, and H(t) = [h xx h xy h yx h yy R′ is the transformation matrix as follows:

[0010]

[0011] Where α and β are the tilt angle and lateral tilt angle of the machining tool, the tool-workpiece contact area is obtained by using the solid model method, and the tilt angle and lateral tilt angle of the milling cutter cutting edge are obtained.

[0012] By employing the Lagrange transformation, we make x′(t)=[q(t),Mdq(t) / dt+Cq(t) / 2] T Therefore, the system dynamics equations can be written as time-delay differential equations in modal space:

[0013] x′(t)=A0x(t)+B(t)x(t)-B(t)x(tT) (3)

[0014] Where A0 is a constant matrix; B(t) is a periodic coefficient matrix that varies with milling time, i.e., B(t+T)=B(t).

[0015]

[0016] Step 2: Based on the general solution of the spatial time-delay differential equation, a stability domain analysis is performed using the time-history fine integration method based on third-order Newton-Hermitian interpolation. According to the spatial time-delay differential equation obtained in Step 1, let v(t) = B(t)x(t) and θ(tT) = B(t)x(tT), and set the initial condition as x0 = x(0). Then, the general solution of the five-axis system can be expressed in integral form as follows:

[0017]

[0018] The milling cycle T is discretized into m infinitesimal time intervals, so the step size of each time interval is τ = T / m. In each time interval t ∈ [kτ, (k+1)τ], (k = 0, 1, ..., m), x... k+1 Let x(kτ+τ) be represented by x. k Let x(kτ) be the value, and simplifying formula (5) yields:

[0019]

[0020] Among them, δ=ξ-kτ, δ∈[0, τ].

[0021] Then, the state term and delay term in formula (6) are solved using the third-order Newton interpolation method and the Hermitian interpolation method. First, the state term is solved. The state term v(δ) on the time interval [kτ, (k+1)τ] is interpolated by the third-order Newton method through the interpolation nodes v(kτ+τ), v(kτ), v(kτ-τ), and v(kτ-2τ). The interpolation nodes can be written as v k+1 v k v k-1 and v k-2 The state term v(δ) is approximately expressed as:

[0022] v(δ)=a1v k+1 +b1v k +c1v k-1 +d1v k-2 (7)

[0023] Wherein, the coefficients a1, b1, c1, and d1 are represented as follows:

[0024]

[0025]

[0026]

[0027]

[0028] v k+1 =B k+1 x k+1 v k =B k x k v k-1 =B k-1 x k-1 v k-2 =B k-2 x k-2 (12)

[0029] Next, the delay term is approximated. The delay term θ(δ-T) on the time interval [kτ, (k+1)τ] is Hermitian interpolated through the interpolation nodes θ(kτ+2τ-T), θ(kτ+τ-T), and θ(kτ-T). The interpolation nodes can be represented as θ k-m+2 θ k-m+1 and θ k-m The delay term θ(δ-T) can be approximated as:

[0030] θ(δ-T)=a2θ k-m +b2θ k-m+1 +c2θ k-m+2 (13)

[0031] Wherein, the coefficients a2, b2, and c2 are respectively represented as:

[0032]

[0033]

[0034]

[0035] θ k-m =B k x k-m θ k-m+1 =B k+1 x k-m+1 θ k-m+2 =B k+2 x k-m+2 (17)

[0036] Finally, the time-delay differential equation can be written in the following form.

[0037]

[0038] in,

[0039]

[0040]

[0041]

[0042]

[0043]

[0044]

[0045]

[0046]

[0047] The exponent matrix T1 is calculated using the precise integration method. First, T1 is further rewritten as:

[0048]

[0049] Where Δt = τ / 2 n Generally, n = 20. When the time interval Δt is sufficiently small, T1 can be approximated using a Taylor expansion, i.e.

[0050]

[0051] Where, Ta = A0Δt + (A0Δt) 2 / 2! +(A0Δt) 3 / 3! +(A0Δt) 4 / 4!

[0052] Combining formulas (27) and (28), the expression for the exponential matrix T1 can be obtained as follows:

[0053]

[0054] Therefore, the increment Ta is obtained through n iterations, and the specific algorithm is written as follows:

[0055]

[0056] If k = 1, x on the left side of equation (18) k-1 and x k-2 Write as x0 and x -1 This corresponds to x in the previous period (tT). m-m and x m-m-1 To construct the state transition matrix of the milling cutter cutting edge in two adjacent milling cycles, when k = m, the time-delay differential equation can be further written in the following form:

[0057]

[0058] The following discrete mapping relation is finally constructed:

[0059]

[0060] in

[0061]

[0062]

[0063] Therefore, the state transition matrix Φ of the milling system over one time period can be calculated as Φ = (D1). -1 D2. Finally, based on Floquet theory, the stability of the milling system is determined by the modulus of the eigenvalues ​​of the state transition matrix Φ. If the modulus of all eigenvalues ​​of the state transition matrix Φ is less than 1, then the milling system is stable; otherwise, the milling system is unstable.

[0064] Step 3: During the milling process, the trajectory of the j-th cutting edge can be represented by the superposition of the relative motion of the milling cutter as follows:

[0065]

[0066] Where x(t) represents the displacement of the geometric center of the milling cutter, and f represents the feed rate.

[0067] Since the displacement of the geometric center of the milling cutter can be divided into static and dynamic displacement, we can obtain

[0068] x(t)=x p (t)+q(t) (36)

[0069] If the milling process is stable, the dynamic displacement q(t) of the milling system will converge to 0, and the milling cutter will generate periodic vibration x. p (t); If the milling process is unstable, the dynamic displacement q(t) of the milling system will diverge exponentially, and the static displacement x p q(t) will also diverge. During milling chatter, as q(t) continuously increases, the milling cutter will disengage from the workpiece, at which point the cutting force is 0, and milling will continue when the cutter returns to its equilibrium position. Therefore, in actual machining, q(t) will not increase indefinitely and has a finite amplitude.

[0070] In the time interval [t] i , t i+1 The discrete dynamical system on [the surface] is:

[0071]

[0072] Where F2 is the static force term, and the state vector z is defined. k = [x1, x2, ..., x m+1 ] T Therefore, the formula can be further expressed as

[0073] D1·z k =D2·z k-m +F2 (38)

[0074] The stability of a milling system is determined by the magnitudes of the eigenvalues ​​of the transition matrix Φ. When the magnitudes of all eigenvalues ​​of the transition matrix Φ are less than 1, indicating a stable milling system, the system is excited by a periodic cutting force, resulting in periodic motion. Using a method similar to the time-domain finite element analysis, let z be the fixed-point vector of the state vector in steady state. * ,make

[0075] z k =z k-m =z * (39)

[0076] Substituting the fixed-point vector into the common formula, we can obtain...

[0077] z * =(D1-D2) -1 F2 (40)

[0078] For different combinations of axial cutting depth and spindle speed, matrices D1, D2, and the static force term F2 can all be calculated, and the surface position error can be predicted using the following formula:

[0079]

[0080] Here, represents the coordinate value of the cutting edge trajectory in the y-direction. A positive SLE value indicates overcutting; a negative SLE value indicates undercutting.

[0081] Step 4: Select some cutting edge trajectory points close to the workpiece's machined surface to form the trajectory points to be interpolated and densified (u j，x_trim (t), u j，y_trim (t)), different milling methods can be selected using the following formulas respectively.

[0082]

[0083] Where γ is the adjustment coefficient, usually set to 0.9.

[0084] Next, the cutting edge trajectory points selected above are defined in the x-direction range as x. min =min{u j，x_trim (t)},x max =min{uj，x_trim (t)}, then for the interval [x min x max Discrete partitioning, where the number of discrete points is N. s The interval step size is δx, so the set of coordinate values ​​of the trajectory points in the x-direction is {x}. min x min +δx,x min +2δx,…,x max}

[0085] Then, when the axial height of the milling cutter is z, for each tooth j, respectively with (u j，x_trim (t), u j，y_trim Given nodes (t), the x-coordinates of the trajectory points are calculated using spline interpolation. min x min +δx,x min +2δx,…,x max The corresponding ordinate values ​​constitute the set of interpolation densification points for the cutting edge trajectory over a single time period (x). s (l), y s (l)), l = 1, 2, ..., N s .

[0086] Based on the periodic property of the dynamic response of chatter-free milling, i.e., x s (l) and x s (l)+n rev The y-coordinate value corresponding to T is y s (l), and then n rev (n rev The set of densified points (x) over ≥4 time periods s_n (l), y s_n (l)), l=1,2,…,N s_n , where N s_n For n rev The number of interpolation points over a time period.

[0087] Finally, x was chosen. min +T≤x s_n (l)≤n rev x max - The denser points of the milling cutter cutting edge trajectory during the intermediate cycle of -T constitute a new set of denser points (x s_n_trim (l), y s_n_trim (l)), l=1,2,…,N s_n_trim , where N s_n_trim This represents the number of interpolation points after clipping. This is achieved by using the new set of densified points (x...). s_n_trim (l), y s_n_trim(l) Divided according to the direction of the milling cutter axis, forming N a A subset of dense points (x) s_n_trim (i, l), y s_n_trim (i, l)), i=1, 2,…,N a .

[0088] When the axial height of the milling cutter is z, the y-coordinate values ​​corresponding to all cutter teeth j with the same x-coordinate value are compared, and the densification point closest to the workpiece forming surface is selected to form the final workpiece surface topography point cloud (x surf (i, l), y surf (i, l)), its expression is

[0089]

[0090] Step 5: Based on the point cloud (x) that contributes to the surface topography of the workpiece surf (i, l), y surf (i, l)) is used to calculate the surface roughness, and its expression is:

[0091]

[0092] Compared with the prior art, the beneficial effects of the present invention are as follows:

[0093] (1) The time history fine integration method based on third-order Newton-Hermit interpolation proposed in this invention can significantly improve the calculation efficiency and reduce the program calculation time while ensuring the calculation accuracy when predicting milling stability.

[0094] (2) The algorithm proposed in this invention combines the prediction of milling chatter stability domain with the prediction of surface position error, surface roughness and surface texture that reflect the constraints of machining surface quality, and realizes the synchronous prediction of five-axis machining chatter stability domain with surface position error, surface roughness and surface texture.

[0095] (3) This invention can provide more comprehensive theoretical support for optimizing the selection of machining process parameters to achieve high-efficiency and high-precision five-axis milling of weak rigidity complex components. Attached Figure Description

[0096] The accompanying drawings are provided to further illustrate the invention and constitute a part of this invention. The following illustrative embodiments and descriptions of the invention are used to explain the invention and do not constitute a limitation thereof.

[0097] Figure 1 This is a diagram of the dynamic model of a flexible tool and a flexible workpiece.

[0098] Figure 2(a) and Figure 2(b) are surface position error diagrams for different milling methods; where Figure 2(a) is climb milling and Figure 2(b) is conventional milling.

[0099] Figure 3 This is a flowchart of the algorithm for the synchronous prediction method of milling stability region, surface roughness, position error and texture of the present invention. Detailed Implementation

[0100] The embodiments of the present invention will now be described in detail with reference to the accompanying drawings, clearly indicating the purpose, technical solution, and advantages of the present invention. It should be noted that the following detailed descriptions of the embodiments are illustrative.

[0101] like Figure 3 As shown, the method for synchronous prediction of milling stability region, surface roughness, position error, and texture in this embodiment includes the following steps:

[0102] Step (1), as follows Figure 1 As shown, considering that the stiffness of the machine tool in the z-direction is much greater than that in the x and y directions, the dynamic displacement in the z-direction can be ignored, and thus the dynamic cutting force of the tool in the z-direction is 0. Therefore, the dynamic equations for two-degree-of-freedom five-axis milling are established:

[0103]

[0104] Where, q x (t) and q y (t) represents the vibration displacement of the tool in the x and y directions, q′ x (t) and q′ y (t) represents the vibration velocity of the tool in the x and y directions, q″ x (t) and q″ y (t) represents the vibration acceleration of the tool in the x and y directions, m x c x and k x Let m be the modal mass, damping, and stiffness of the system in the x-direction; y c y and k y Let a be the modal mass, damping, and stiffness of the system in the y-direction; p λ is the axial cutting depth; T is the time delay, which is usually the cutting cycle, i.e., T = 60 / (ΩN). t ), where Ω and N t These represent the spindle speed and the number of tool teeth, respectively. The cutting force coefficient matrix for five-axis milling is as follows:

[0105]

[0106] Where R′ is the coordinate transformation matrix, α is the tool rake angle, β is the tool side rake angle, and H(t) is the dynamic cutting force coefficient matrix.

[0107]

[0108] Where h xx (t), h xy (t), h yx (t) and h yy (t) is the dynamic cutting force coefficient, expressed as follows:

[0109]

[0110]

[0111]

[0112]

[0113] in, Let K be the angle of the j-th tooth of the cutting tool. tc and K rc The cutting force coefficients κ(z) are for the tangential and normal directions, respectively, and the axial contact angle is denoted by κ(z).

[0114] Then, through the Lagrange transformation, let x′(t) = [q(t), Mdq(t) / dt + Cq(t) / 2] T , where q(t)=[q x (t) q y (t)] T The superscript T denotes the transpose of the matrix. Therefore, formula (1) can be further rewritten in state-space form as follows:

[0115] x′(t)=A0x(t)+B(t)x(t)-B(t)x(tT) (8)

[0116] Where A0 is a constant matrix; B(t) is a periodic coefficient matrix that varies with milling time, i.e., B(t+T)=B(t), where

[0117]

[0118]

[0119]

[0120] Where w x and w y Let ζ be the natural frequencies of the tool in the x and y directions, respectively. x and ζ y These are the damping ratios of the tool in the x and y directions, respectively.

[0121] Step (2) solves the time-delay differential equation of the dynamic model of the two-degree-of-freedom five-axis milling described above, and calculates it using the time history fine integration method based on third-order Newton-Hermit interpolation.

[0122] Let v(t) = B(t)x(t) and θ(tT) = B(t)x(tT), and set the initial condition as x0 = x(0). Then the general solution of the system dynamic equation can be expressed in integral form as follows:

[0123]

[0124] First, the milling cycle T of the milling cutter cutting edge is discretized into m infinitesimal time intervals, so the step size of each infinitesimal time interval is τ = T / m. Therefore, in each time interval t∈[kτ, (k+1)τ], (k = 0, 1, ..., m), formula (12) can be further rewritten as:

[0125]

[0126] Then, for the convenience of subsequent calculations, x k+1 Represented as x(kτ+τ), x k Represented as x(kτ), the formula can be further simplified to obtain:

[0127]

[0128] Among them, δ=ξ-kτ, δ∈[0, τ].

[0129] Next, the state term v(δ) over the time interval [kτ, (k+1)τ] is approximated using the third-order Newton interpolation method. This can be achieved through the nodal values ​​v(kτ+τ), v(kτ), v(kτ-τ), and v(kτ-2τ) (abbreviated as v...). k+1 v k v k-1 and v k-2 ) is represented as:

[0130] v(δ)=a1v k+1 +b1v k +c1v k-1 +d1v k-2 (15)

[0131] The coefficients a1, b1, c1, and d1 are expressed as follows:

[0132]

[0133]

[0134]

[0135]

[0136] v k+1 =B k+1 x k+1 v k =B k x k v k-1 =B k-1 x k-1 v k-2 =B k-2 x k-2 (20)

[0137] Similarly, the delay term θ(δ-T) over the time interval [kτ, (k+1)τ] can be expressed as θ(kτ+2t-T), θ(kτ+τ-T), and θ(kt-T) (respectively θ). k-m+2 θ k-m+1 and θ k-m By performing Hermitian interpolation, we can approximately solve for the following:

[0138] θ(δ-T)=a2θ k-m +b2θ k-m+1 +c2θ k-m+2 (twenty one)

[0139] The coefficients a2, b2, and c2 are expressed as follows:

[0140]

[0141]

[0142]

[0143] θ k-m =B k x k-m θ k-m+1 =B k+1 x k-m+1 θ k-m+2 =B k+2 x k-m+2 (25)

[0144] Combining formulas (13), (14), and (21), the time-delay differential equation can finally be written in the following form:

[0145]

[0146] in,

[0147]

[0148]

[0149]

[0150]

[0151]

[0152]

[0153]

[0154]

[0155] Similarly, in order to construct the state transition matrix of the milling cutter cutting edge in two adjacent milling cycles, the responses on both sides of equation (26) under the conditions of k=1, k=2, and k=m can be expressed as follows:

[0156]

[0157]

[0158]

[0159] Combining the formulas, the following discrete mapping relation is finally constructed:

[0160]

[0161] in

[0162]

[0163]

[0164] Therefore, the state transition matrix Φ of the milling system over one time period can be calculated as follows:

[0165] Φ=(D1) -1 D2 (41)

[0166] Finally, according to Floquet's theory, the stability of the system is determined by the moduli of the eigenvalues ​​of the transition matrix Φ. If the moduli of all eigenvalues ​​of the transition matrix Φ are less than 1, then the milling system is stable; otherwise, the milling system is unstable.

[0167]

[0168] Step (3) During the milling process, the trajectory of the j-th cutting edge can be represented by the superposition of the relative motion of the milling cutter, as follows:

[0169]

[0170] Where x(t) represents the displacement of the geometric center of the milling cutter, f represents the feed rate, and D is the diameter of the milling cutter.

[0171] Since the displacement of the geometric center of the milling cutter can be divided into static and dynamic displacement, we can obtain

[0172] x(t)=x p (t)+q(t) (44)

[0173] If the milling process is stable, the dynamic displacement q(t) of the milling system will converge to 0, and the milling cutter will generate periodic vibration x. p (t); If the milling process is unstable, the dynamic displacement q(t) of the milling system will diverge exponentially, and the static displacement x p q(t) will also diverge. During milling chatter, as q(t) continuously increases, the milling cutter will disengage from the workpiece, at which point the cutting force is 0, and milling will continue when the cutter returns to its equilibrium position. Therefore, in actual machining, q(t) will not increase indefinitely and has a finite amplitude.

[0174] In the time interval [t] i , t i+1 The discrete dynamical system on [the surface] is:

[0175]

[0176] Where F2 is the static force term, and the state vector z is defined. k = [x1, x2, ..., x m+1 ] T Therefore, the formula can be further expressed as:

[0177] D1·z k =D2·z k-m +F2 (46)

[0178] The stability of a milling system is determined by the magnitudes of the eigenvalues ​​of the transition matrix Φ. When the magnitudes of all eigenvalues ​​of the transition matrix Φ are less than 1, indicating a stable milling system, the system is excited by a periodic cutting force, resulting in periodic motion. Using a method similar to the time-domain finite element analysis, let z be the fixed-point vector of the state vector in steady state. * ,make

[0179] z k =z k-m =z * (47)

[0180] Substituting the common fixed-point vector into the common vector, we get:

[0181] z * =(D1-D2) -1 F2 (48)

[0182] For different combinations of axial cutting depth and spindle speed, matrices D1, D2, and the static force term F2 can all be calculated, and the surface position error SLE can be predicted using the following formula:

[0183]

[0184] Among them, u j，y (t) represents the coordinate value of the cutting edge trajectory in the y direction. When SLE is positive, it indicates overcutting; when SLE is negative, it indicates undercutting, as shown in Figure 2(a) and Figure 2(b).

[0185] Step (4) Select some cutting edge trajectory points close to the workpiece machining surface to form the trajectory points to be interpolated and densified (u j，x_trim (t), u j，y_trim (t)), different milling methods can be selected using the following formulas:

[0186]

[0187] Where γ is the adjustment coefficient, usually set to 0.9.

[0188] Next, the cutting edge trajectory points selected above are defined in the x-direction range as x. min =min{u j，x_trim (t)},x max =min{u j，x_trim (t)}, then for the interval [x min x max Discrete partitioning, where the number of discrete points is N. s The interval step size is δx, so the set of coordinate values ​​of the trajectory points in the x-direction is {x}. min x min +δx,x min +2δx,…,x max}

[0189] Then, when the axial height of the milling cutter is z, for each tooth j, respectively with (u j，x_trim (t), u j，y_trim Given nodes (t), the x-coordinates of the trajectory points are calculated using spline interpolation. min x min +δx,x min +2δx,...,x max The corresponding ordinate values ​​constitute the set of interpolation densification points for the cutting edge trajectory over a single time period (x).s (l), y s (l)), l=1,2,…,N s .

[0190] Based on the periodic property of the dynamic response of chatter-free milling, i.e., x s (l) and x s (l)+n rev The y-coordinate value corresponding to T is y s (l), and then n rev (n rev The set of densified points (x) over ≥4 time periods s_n (l), y s_n (l)), l=1,2,…,N s_n , where N s_n For n rew The number of interpolation points over a time period.

[0191] Finally, x was chosen. min +T≤x s_n (l)≤n rev x max - The denser points of the milling cutter cutting edge trajectory during the intermediate cycle of -T constitute a new set of denser points (x s_n_trim (l), y s_n_trim (l)), l=1,2,…,N s_n_trim , where N s_n_trim This represents the number of interpolation points after clipping. This is achieved by using the new set of densified points (x...). s_n_trim (l), y s_n_trim (l) Divided according to the direction of the milling cutter axis, forming N a A subset of dense points (x) s_n_trim (i, l), y s_n_trim (i, l)), i=1, 2,…,N a .

[0192] When the axial height of the milling cutter is z, the y-coordinate values ​​corresponding to all cutter teeth j with the same x-coordinate value are compared, and the densification point closest to the workpiece forming surface is selected to form the final workpiece surface topography point cloud (x surf (i, l), y surf (i, l)), its expression is:

[0193]

[0194] Based on the point cloud (x) that contributes to the surface morphology of the workpiece surf (i, l), y surf (i, l)) are used to calculate the surface roughness Ra(z), and its expression is:

[0195]

[0196] While some embodiments of the overall technical concept of this disclosure have been shown and described, those skilled in the art will understand that changes may be made to these embodiments without departing from the principles and spirit of the overall technical concept, the scope of which is defined by the claims and their equivalents.

Claims

1. A method for synchronously predicting milling stability domain, surface roughness, positional error, and texture, characterized in that, Includes the following steps: Step 1: Construct a dynamic model of the five-axis milling system; Step 2: A time-history fine integration method based on third-order Newton-Hermitian interpolation is proposed for stability domain analysis; Step 3: Calculate the relative vibration displacement between the tool and the workpiece and the running trajectory of the milling cutter cutting edge at the discrete points of the tool rotation cycle using fixed point theory, and simultaneously predict the surface position error. Step 4: Use spline interpolation to densify the cutting edge trajectory generated on the workpiece surface, and finally construct the surface topography point cloud; Step 5: Calculate the surface roughness based on the point cloud that constitutes the surface morphology of the workpiece. Step 1 includes: establishing the dynamic equations of a five-axis milling system considering flexible cutting tools and flexible workpieces, expressed as: Where M, C, and K are the modal mass matrix, modal damping matrix, and modal stiffness matrix of the system, respectively, and q(t) represents the modal coordinates of the milling system; a p The axial cutting depth is represented by T; T represents the time delay, which is also the milling cycle, i.e., T = 60 / (ΩN). t N is the number of teeth on the milling cutter, and Ω is the spindle speed in rpm. This is a five-axis cutting force coefficient matrix, and H(t)=[h xx h xy h yx h yy R′ is the transformation matrix as follows: Where α and β are the tilt angle and side tilt angle of the machining tool, the tool-workpiece contact area is obtained by using the solid model method, and the tilt angle and side tilt angle of the milling cutter cutting edge are obtained; By employing the Lagrange transformation, we make x′(t)=[q(t),Mdq(t) / dt+Cq(t) / 2] T Therefore, the system dynamics equations can be written as time-delay differential equations in modal space: x′(t)=A0x(t)+B(t)x(t)-B(t)x(tT) (8) Where A0 is a constant matrix; B(t) is a periodic coefficient matrix that varies with milling time, i.e., B(t+T)=B(t); 2. The method according to claim 1, characterized in that, In step 2, the stability region analysis based on the time-history fine integration method using third-order Newton-Hermitian interpolation is performed as follows: Based on the spatial time-delay differential equation, let v(t) = B(t)x(t) and θ(tT) = B(t)x(tT), and set the initial condition as x0 = x(0), the general solution of the five-axis system can be expressed in integral form as follows: The milling cycle T is discretized into m infinitesimal time intervals, so the step size of each time interval is τ = T / m; in each time interval t∈[kτ, (k+1)τ], where k = 0, 1, ..., m, x k+1 Let x(kτ+τ) be represented by x. k Let x(kτ) be the value, and simplify equation (12) to obtain Among them, δ = ξ-kτ, δ∈[0, τ]; Then, the state term and delay term in formula (14) are solved using the third-order Newton interpolation method and the Hermitian interpolation method. First, the state term is solved. The state term v(δ) at the time interval [kτ, (k+1)τ] is interpolated by the third-order Newton method through the interpolation nodes v(kτ+τ), v(kτ), v(kτ-τ), and v(kτ-2τ). The interpolation nodes v(kτ+τ), v(kτ), v(kτ-τ), and v(kτ-2τ) are simplified to v k+1 v k v k-1 and v k-2 The state term v(δ) is approximately expressed as: v(δ)=a1v k+1 +b1v k +c1v k-1 +d1v k-2 (15) Wherein, the coefficients a1, b1, c1, and d1 are represented as follows: v k+1 =B k+1 x k+1 v k =B k x k v k-1 =B k-1 x k-1 v k-2 =B k-2 x k-2 (20) Next, the delay term is approximated. The delay term θ(δ-T) on the time interval [kτ, (k+1)τ] is Hermitian interpolated through the interpolation nodes θ(kτ+2τ-T), θ(kτ+τ-T), and θ(kτ-T). The interpolation nodes θ(kτ+2τ-T), θ(kτ+τ-T), and θ(kτ-T) are respectively denoted as θ k-m+2 θ k-m+1 and θ k-m The delay term θ(δ-T) is approximately expressed as: θ(δ-T)=a2θ k-m +b2θ k-m+1 +c2θ k-m+2 (21) Wherein, the coefficients a2, b2, and c2 are respectively represented as: i k-m =B k x k-m i k-m+1 =B k+1 x k-m+1 i k-m+2 =B k+2 x k-m+2 (25) Finally, the time-delay differential equation can be written in the following form: in, The exponent matrix T1 is calculated using the precise integration method. First, T1 is further rewritten as: Where Δt = τ / 2 n Let n = 20; when the time interval Δt is sufficiently small, T1 is approximated using a Taylor expansion, i.e. Where, Ta = A0Δt + (A0Δt) 2 / 2! +(A0Δt) 3 / 3! +(A0Δt) 4 / 4! ; Combining the two formulas above, the expression for the exponential matrix T1 is as follows: Therefore, the increment Ta is obtained through n iterations, and the specific algorithm is written as follows: for(i=1;i≤n;i++) Ta = 2Ta + Ta × Ta; T1 = I + Ta; If k = 1, x on the left side of equation (26) k-1 and x k-2 Write as x0 and x -1 This corresponds to x in the previous period (tT). m-m and x m-m-1 To construct the state transition matrix of the milling cutter cutting edge in two adjacent milling cycles, when k = m, the time-delay differential equation can be further written in the following form: The following discrete mapping relation is finally constructed. in Therefore, the state transition matrix Φ of the milling system over one time period is calculated as Φ=(D1) -1 D2. Finally, based on Floquet theory, the stability of the milling system is determined by the modulus of the eigenvalues ​​of the state transition matrix Φ. If the modulus of all eigenvalues ​​of the state transition matrix Φ is less than 1, then the milling system is stable; otherwise, the milling system is unstable.

3. The method according to claim 1, characterized in that, Step 3 includes: In the milling process, the trajectory of the j-th cutting edge can be represented by the superposition of the relative motion of the milling cutter as follows: Where x(t) represents the displacement of the geometric center of the milling cutter, f represents the feed rate, and D is the diameter of the milling cutter; since the displacement of the geometric center of the milling cutter can be divided into static and dynamic displacements, we obtain: x(t)=x p (t)+u(t) (44) If the milling process is stable, the dynamic displacement u(t) of the milling system will converge to 0, and the milling cutter will generate periodic vibration x. p (t); If the milling process is unstable, the dynamic displacement u(t) of the milling system will diverge exponentially, and the static displacement x p u(t) will also diverge; when milling chatter occurs, as u(t) increases continuously, the milling cutter will disengage from the workpiece, at which point the cutting force is 0, and when the cutter returns to the equilibrium position, milling will continue; therefore, in actual machining, u(t) will not increase indefinitely, but has a finite amplitude. In the time interval [t] i , t i+1 The discrete dynamical system on [the surface] is: Where F2 is the static force term, and the state vector z is defined. k = [x1, x2, ..., x m+1 ] T Therefore, the formula can be further expressed as: D1·z k =D2·z k-m +F2 (46) The stability of a milling system is determined by the magnitudes of the eigenvalues ​​of the transition matrix Φ. When the magnitudes of all eigenvalues ​​of the transition matrix Φ are less than 1, indicating a stable milling system, the system is excited by a periodic cutting force, resulting in periodic motion. Using a method similar to the time-domain finite element analysis, let z be the fixed-point vector of the state vector in steady state. * ,make With k =z k-m =z * (47) Substituting the common fixed-point vector into the common vector, we get: z * (D1-D2) -1 F2 (48) For different combinations of axial cutting depth and spindle speed, matrices D1, D2, and the static force term F2 are all calculated, and the surface position error is predicted using the following formula: Here, SLE represents the coordinate value of the cutting edge trajectory in the y direction; and when SLE is positive, it indicates overcutting; when SLE is negative, it indicates undercutting.

4. The method according to claim 1, characterized in that, Step 4 includes: Selecting some cutting edge trajectory points close to the workpiece's machined surface to form the trajectory points to be interpolated and densified (u j，x_trim (t), u j，y_trim (t)), different milling methods can be selected using the following formulas respectively. Where γ is the adjustment coefficient, let γ = 0.9; Next, the cutting edge trajectory points selected above are defined in the x-direction range as x. min =min{u j，x_trim (t)},x max =min{u j，x_trim (t)}, then for the interval [x min x max Discrete partitioning, where the number of discrete points is N. s The interval step size is δx, so the set of coordinate values ​​of the trajectory points in the x-direction is {x}. min x min +δx,x min +2δx,...,x max }; Then, when the axial height of the milling cutter is z, for each tooth j, respectively with (u j，x_trim (t), u j，y_trim Given nodes (t), the x-coordinates of the trajectory points are calculated using spline interpolation. min x min +δx,x min +2δx,...,x max The corresponding ordinate values ​​constitute the set of interpolation densification points for the cutting edge trajectory over a single time period (x). s (l), y s (l)), l = 1, 2, ..., N s ; Based on the periodic property of the dynamic response of chatter-free milling, i.e., x s (l) and x s (l)+n rev The y-coordinate value corresponding to T is y s (l), and then n rev The set of densification points over a time period (x) s_n (l), y s_n (l)), l = 1, 2, ..., N s_n , where n rev ≥4, N s_n For n rev The number of interpolation points over a time period; Finally, x was chosen. min +T≤x s_n (l)≤n rev x max - The denser points of the milling cutter cutting edge trajectory during the intermediate cycle of -T constitute a new set of denser points (x s_n_trim (l), y s_n_trim (l)), l = 1, 2, ..., N s_n_trim , where N s_n_trim The number of interpolation points after clipping; by adding a new set of densified points (x s_n_trim (l), y s_n_trim (l) Divided according to the direction of the milling cutter axis, forming N a A subset of dense points (x) s_n_trim (i, l), y s_n_trim (i, l)), i=1, 2,...,N a ; When the axial height of the milling cutter is z, the y-coordinate values ​​corresponding to all cutter teeth j with the same x-coordinate value are compared, and the densification point closest to the workpiece forming surface is selected to form the final workpiece surface topography point cloud (x surf (i, l), y surf (i, l)), its expression is:

5. The method according to claim 4, characterized in that, In step 5, based on the point cloud (x) that participates in constituting the surface morphology of the workpiece surf (i, l), y surf (i, l)) is used to calculate the surface roughness. The formula for calculating surface roughness is:

Images