Milling chatter stability prediction method

Abstract

The invention relates to the field of milling, in particular to a method for predicting milling chatter stability, which comprises the following steps of: 1, establishing a milling process dynamic differential equation containing a regeneration effect; 2, converting a time-delay kinetic equation of the milling system into a state-space equation; 3, solving through a third-order Hermite interpolation polynomial approximation state item and a second-order Newton interpolation polynomial approximation time delay item; and 4, constructing a Floquet transition matrix according to the solution of the state-space equation, and judging the stability of the milling system. The invention provides a stable lobe graph drawing method based on an improved full discrete method. According to the method, the calculation precision is remarkably improved, the calculation efficiency is effectively improved, and the requirements in practical application can be better met.

Description

Technical Field

[0001] The invention relates to the technical field of milling processing, and in particular to a method for predicting milling chatter stability. Background Art

[0002] Milling technology is widely used in the processing of complex large parts such as wind power manufacturing, aerospace, and marine propellers due to its advantages of high efficiency, high yield, high precision, and low cost. It has become one of the key technologies to ensure high-precision and high-quality workpieces. However, in the milling process, regenerative chatter is a common self-excited instability phenomenon, which seriously affects the surface quality of the workpiece, processing efficiency, and tool life. Therefore, it is crucial to take effective measures to avoid or suppress chatter. An effective method to suppress regenerative chatter is to avoid the occurrence of chatter by selecting appropriate process parameters based on the stability lobe diagram of the milling process. Based on this, it is of great practical significance to analyze and study the stability of the milling process.

[0003] At present, the stability prediction methods of the milling process mainly include numerical simulation method, frequency domain method and time domain method. The numerical simulation method sets the initial conditions for the time-delay differential equation of the milling process, uses the step-by-step numerical integration method to simulate the time domain of several cycles, and then uses the simulation signal to determine whether chatter occurs. For example, the existing literature "An improved time domain simulation for dynamic milling at small radial immersions" established an improved numerical simulation model and proposed to determine the conditions for the occurrence of chatter by dimensionless coefficients, that is, the ratio of the maximum dynamic to static cutting thickness. The numerical simulation method has high accuracy in the stability prediction of the milling process, but this method can only simulate under one set of cutting parameters at a time, so it takes a long time to draw a complete stability domain diagram, and the calculation efficiency is low. The zero-order frequency domain method has high calculation efficiency, but due to the use of zero-order approximate cutting force, its prediction accuracy has large errors, especially under the condition of small radial cutting depth.

[0004] In the frequency domain model of stability research, the existing literature "Analytical prediction of stability lobes in milling" proposed a method to decompose the periodic cutting force coefficient matrix using Fourier series expansion technology, extract the zero-order term of the Fourier series as an approximate estimate of the periodic cutting force, and derive the analytical expression of the critical axial cutting depth when no chatter occurs by solving the system characteristic equation.

[0005] The time domain method directly discretizes the time delay differential equation in the time domain and transforms it into a set of ordinary differential equations. The infinite dimensional characteristic operator in the time delay differential equation is approximated by the finite dimensional transfer matrix of adjacent tooth periods, and then the stability of the milling system is judged according to the Floquet theorem. The existing literature "Semi-discretization method for delayed systems" proposes a semi-discretization method, which transforms the time delay differential equation into a set of ordinary differential equations by averaging the time delay terms and periodic coefficient terms in the time delay differential equation to judge the stability of the system; while the literature "A full-discretization method for prediction of milling stability" expands the semi-discretization method into a full-discretization method to improve the accuracy and computational efficiency of stability prediction. The semi-discretization method has good applicability in the milling process and can provide high prediction accuracy and numerical stability. However, this method needs to calculate the exponential matrix when scanning the cutting depth, resulting in low computational efficiency, and its accuracy still has room for improvement. In contrast, the fully discrete method greatly improves the computational efficiency by decomposing the coefficient matrix into a constant matrix and a periodic variation matrix, without calculating the exponential matrix when scanning the cutting depth. However, there is still room for further improvement in the accuracy of the fully discrete method. Summary of the invention

[0006] To achieve the above purpose, the present invention proposes a stability lobe diagram drawing method based on an improved full discrete method. This method not only significantly improves the calculation accuracy, but also effectively improves the calculation efficiency, and can better meet the needs in practical applications. The present invention provides the following technical solutions:

[0007] A method for predicting milling chatter stability, comprising:

[0008] In the first step, the differential equation of the milling process dynamics including the regeneration effect is established;

[0009] The second step is to transform the time-delay dynamics equation of the milling system into the state-space equation;

[0010] The third step is to solve the problem by approximating the state term with a third-order Hermite interpolation polynomial and the delay term with a second-order Newton interpolation polynomial.

[0011] In the fourth step, the Floquet transition matrix is ​​constructed based on the solution of the state-space equation to determine the stability of the milling system.

[0012] As a further solution of the present invention: the first step is to establish a milling dynamics model based on the milling regenerative chatter principle under the framework of structural dynamics. The milling dynamics model is shown in Figure 2The milling process dynamics differential equations are the following differential equations that contain both parameter excitation and time delay terms: Where M, C, and K are the modal mass matrix, damping matrix, and stiffness matrix of the tool end, respectively; q(t) is the tool modal coordinate, and the vibration coefficient is normalized at the tool tip, that is, q(t) = [x(t), y(t)] T ; T is the time delay and is equal to the cutting cycle of the cutter teeth, the matrix K c (t) is expressed as: (2), where K t and K r are the tangential and radial cutting force coefficients, The function used to determine whether the cutting edge is within the effective cutting thickness is defined as: and are the starting cutting angle and the ending cutting angle of milling, respectively, which can be determined according to the radial intrusion ratio of the cutter teeth of the milling system and the milling direction as shown in the following formula (4):

[0013] As a further solution of the present invention: In the second step, the time-delay dynamics equation of the milling system is converted into a state space equation: using the conversion formula in the Hamiltonian system, let x(t)=[q(t)p(t)] T , formula (1) can be expressed in the following state space form: A0 is a constant matrix indicating that the system does not change with time, and B(t) is a periodic coefficient matrix determined by the dynamic cutting force considering the regeneration effect, and satisfies B(t) = B(t + T). T is the time period and is equal to the time delay τ. The tooth cutting cycle is divided into m time periods, and the interval of each time period is Δt = τ / m.

[0014] In each time period t∈[t i ,t i+1 ], (i=1,2,3…m), formula (5) can be expressed as: Formula (6) can be equivalently expressed as: Considering the special form of equation (5), x(t) at t = t i At t = t i+1 The derivative of can be expressed as:

[0015] As a further solution of the present invention: in the third step, the state item x(t) is in the time interval [t i ,t i+1 ], the state term x(t i ), x(t i+1 ) and derivatives The two-point cubic Hermite interpolation is expressed as: x(t) = a1x i +b1x i+1 +c1x i-m +d1x i-m+1 (10);

[0016] a, b, c, d are coefficient terms, ξ is the integral of the time term; use linear interpolation to find the periodic coefficient term B(ξ) in the time interval [t i ,t i+1 ] is linearly approximated with the corresponding boundary value, that is, (15); Use the second-order Newton interpolation polynomial to approximate the time delay term, node t i-m ,t i+1-m ,t i+2-m The corresponding value is x i-m ,x i+1-m ,x i+2-m Based on the second-order Newton interpolation polynomial, the time delay term x(t-τ) can be expressed as: x(t-τ) = a2x i-m +b2x i+1-m +c2x i+2-m (16); In formula (16), Substituting formula (10)-(19) into formula (5), we get the total iteration formula: Q i x i+1 =H i x i +H i-m x i-m +H i-m+1 x i-m+1 +H i-m+2 x i-m+2 (20); where Q i =(I-(K1B i +K2B i+1 ))(twenty one); H i-m =K5B i +K6B i+1 (23); H i-m+1 =K7B i +K8B i+1 (24); H i-m+2 =K9B i +K 10 B i+1 (25), Q i , H i , H i-m , H i-m+1 , H i-m+2is the time period matrix B i and B i+1 Related coefficient matrix; I represents the identity matrix; K1-K 10 is a constant matrix.

[0017] As a further solution of the present invention: In the fourth step, according to formula (20), the following discrete mapping can be obtained: i+1 =D i y i (26); where y i =[y i ,y i-1 ,y i-2 …y i-m ] T (27); If the matrix Q i is non-singular, then D i It can be expressed as The state transition matrix Φ of the system in a single time period can be obtained by the matrix sequence D i (i=1,2,…m), constructed, that is, y m =Φy0(29); where Φ = D m-1 D m-2 …D1D0(30); if Q i If is a singular matrix, the generalized inverse of the matrix can be used instead of the inverse. According to Floquet theory, the stability of the system can be determined by the eigenvalues ​​of the state transfer matrix Φ: when the modulus of any eigenvalue is greater than 1, the system is unstable; when the modulus of all eigenvalues ​​of the state transfer matrix Φ is less than 1, the system is stable. The flutter stability lobe diagram of the milling flutter stability prediction method based on the third-order Hermite second-order Newton interpolation is drawn using MATLAB software as shown in the figure. Figure 3 As shown in the figure, the flutter stability lobe diagram of the first-order full discrete method is as follows Figure 4 As shown. For reference, the stability boundary calculated by the first-order semi-discrete method when the number of discrete intervals m = 200 is used as the benchmark and is given as a black thick line in the figure. Figure 3 and Figure 4 It can be seen that the milling chatter stability prediction method based on third-order Hermite and second-order Newton interpolation can more accurately describe the relationship between the limiting stability cutting depth and the tool speed. The convergence rate is an important indicator for evaluating the quality of the chatter stability lobe diagram calculation method, which can reflect the local error |μ| between the absolute value of the maximum critical eigenvalue of the state transfer matrix and the exact eigenvalue |μ0|. The exact eigenvalue μ0 is calculated using the first-order semi-discrete method when the discrete interval number m=200 as the modulus of the approximate exact critical eigenvalue. In order to clearly demonstrate the convergence of different methods, the present invention sets the same processing parameters and compares the convergence of the proposed method with the first-order fully discrete method under different discrete numbers. The relevant comparison results are shown in the figure. Figure 5As shown, from Figure 5 It can be seen that compared with the existing method, the method of the present invention has a faster convergence rate. The analysis results show that under the premise of ensuring the calculation accuracy, the proposed method can further improve the calculation efficiency.

[0018] Compared with the prior art, the present invention has the following beneficial effects:

[0019] The present invention has important practical significance in actual application in milling processing. The verified test results of this method show that the accuracy and calculation efficiency of vibration stability prediction are significantly improved, which means that when setting processing parameters, it is possible to more accurately predict whether vibration will occur in a shorter time, which not only helps to ensure the quality and accuracy of the processing, but also effectively improves the overall production efficiency and further optimizes the manufacturing process. BRIEF DESCRIPTION OF THE DRAWINGS

[0020] Figure 1 Flow chart of a method for predicting milling chatter stability in an embodiment of the present invention.

[0021] Figure 2 This is a milling regenerative chatter dynamics model diagram obtained by the method for predicting milling chatter stability in an embodiment of the present invention.

[0022] Figure 3 It is a stability lobe diagram of the third-order Hermite and second-order Newton interpolation method in the prediction method of milling chatter stability in an embodiment of the present invention.

[0023] Figure 4 It is a stability lobe diagram of the first-order full discrete method in the prediction method of milling chatter stability in an embodiment of the present invention.

[0024] Figure 5 This is a graph of convergence results obtained by the prediction method for milling chatter stability in an embodiment of the present invention.

[0025] Figure 6 It is a stability lobe diagram obtained by the prediction method of milling chatter stability in an embodiment of the present invention.

[0026] Figure 7 This is a part surface diagram obtained by performing milling processing by selecting the corresponding spindle speed and axial cutting depth in the stable area below the black curve in an embodiment of the present invention.

[0027] Figure 8 This is a part surface diagram obtained by selecting the spindle speed and axial cutting depth for milling in the unstable area above the black curve in an embodiment of the present invention. DETAILED DESCRIPTION

[0028] The following will be combined with the drawings in the embodiments of the present invention to clearly and completely describe the technical solutions in the embodiments of the present invention. Obviously, the described embodiments are only part of the embodiments of the present invention, not all of the embodiments. Based on the embodiments of the present invention, all other embodiments obtained by ordinary technicians in this field without creative work are within the scope of protection of the present invention.

[0029] Milling stability analysis usually solves the dynamic equations of the milling process and combines the dynamic characteristics of the cutting force, tool, workpiece and machine tool to optimize the processing parameters (such as spindle speed, cutting depth, etc.) to avoid the occurrence of chatter. However, the existing milling stability prediction algorithms still have deficiencies in computational efficiency and prediction accuracy, especially when faced with complex processing conditions and high-precision requirements. To solve this problem, the present invention proposes a stability lobe diagram drawing method based on an improved full discrete method, which can more effectively perform milling chatter stability analysis and prediction. Through the improved numerical discretization technology, this method significantly improves the computational efficiency of the algorithm while improving the computational accuracy, thereby providing a more accurate and efficient solution for stability analysis in actual processing.

[0030] The specific implementation of the present invention is described in detail below in conjunction with specific embodiments.

[0031] Example 1, see Figure 1-Figure 8 , set the corresponding cutting parameters, the number of milling cutter teeth is 4, the radial immersion ratio is set to 1, the tangential cutting force coefficient and the radial cutting force coefficient are K t =6×10 8 N / m 2 and K n =6=2×10 8 N / m 2 , the natural frequency f of the milling cutter n =ω n / 2π=922Hz, damping ratio ζ=0.011, modal mass is m=0.03993kg. The discrete numbers of spindle speed and axial cutting depth are 400 and 200. The above process is written into a program using MATLAB software, and the given cutting parameters are input into the program for calculation, and the following is obtained: Figure 6 The stability diagram is shown in Figure 6 In the example, when using an end mill to mill thin-walled aluminum alloy, the corresponding spindle speed and axial cutting depth are selected in the stable area below the black curve to perform milling, and vibration-free cutting can be obtained. After cutting, the surface of the part that meets the requirements is obtained, as shown in the attached figure. Figure 7 As shown in the figure; if the spindle speed and axial cutting depth are selected for milling in the unstable area above the black curve, cutting vibration will occur, and the unsatisfactory part surface will be obtained after cutting. Figure 8shown.

[0032] In addition, it should be understood that although the present specification is described according to implementation modes, not every implementation mode contains only one independent technical solution. This description of the specification is only for the sake of clarity. Those skilled in the art should regard the specification as a whole. The technical solutions in each embodiment may also be appropriately combined to form other implementation modes that can be understood by those skilled in the art.

Claims

1. A method for predicting milling chatter stability, characterized in that: The following steps are involved: In the first step, the differential equation of the milling process dynamics including the regeneration effect is established; The second step is to transform the time-delay dynamics equation of the milling system into the state-space equation; The third step is to solve the problem by approximating the state term with a third-order Hermite interpolation polynomial and the delay term with a second-order Newton interpolation polynomial. In the fourth step, the Floquet transition matrix is ​​constructed based on the solution of the state-space equation to determine the stability of the milling system.

2. The method for predicting milling chatter stability according to claim 1, characterized in that: The expression of the milling process dynamic differential equation is as follows: M, C, and K are the modal mass matrix, damping matrix, and stiffness matrix of the tool end, respectively; q(t) is the tool modal coordinate, and the vibration coefficient is normalized at the tool tip, that is, q(t) = [x(t), y(t)] T ; T is the time delay and is equal to the cutting cycle of the cutter teeth, K c (t) is the coefficient matrix.

3. The method for predicting milling chatter stability according to claim 2, characterized in that: The K c (t) is expressed as: K t and K r are the tangential and radial cutting force coefficients, The function used to determine whether the cutting edge is within the effective cutting thickness is defined as: and are the starting cutting angle and the ending cutting angle of milling, respectively, which can be determined according to the radial intrusion ratio of the cutter teeth of the milling system and the milling direction as shown in formula (4):

4. The method for predicting milling chatter stability according to claim 2, characterized in that: The method of converting the time-delay dynamics equation of the milling system into the state space equation comprises the following steps: using the conversion formula in the Hamiltonian system, let x(t)=[q(t)p(t)] T , formula (1) is expressed in the following state space form: A0 is a constant matrix indicating that the system does not change with time. B(t) is a periodic coefficient matrix determined by the dynamic cutting force considering the regeneration effect, and satisfies B(t) = B(t + T). T is the time period and is equal to the time delay τ. The tooth cutting cycle is divided into m time periods, and the interval between each time period is Δt = τ / m. In each time period t∈[t i ,t i+1 ], (i=1,2,3…m), formula (5) is expressed as: Formula (6) can be expressed equivalently as 5. The method for predicting milling chatter stability according to claim 4, characterized in that: The state term is approximated by the third-order Hermite interpolation polynomial as follows: The state term x(t) is in the time interval [t i ,t i+1 ] using the state term x(t i ), x(t i+1 ) and derivatives The two-point cubic Hermite interpolation is expressed as: x(t) = a1x i +b1x i+1 +c1x i-m +d1x i-m+1 (10); a, b, c, d are coefficient terms, and ξ is the integral of the time term.

6. The method for predicting milling chatter stability according to claim 5, characterized in that: The second-order Newton interpolation polynomial approximates the time lag term to solve the problem as follows: use the linear interpolation method to calculate the periodic coefficient term B(ξ) in the time interval [t i ,t i+1 ] make a linear approximation with the corresponding boundary value, The time delay term is approximated by a second-order Newton interpolation polynomial, node t i-m ,t i+1-m ,t i+2-m The corresponding value is x i-m ,x i+1-m ,x i+2-m Based on the second-order Newton interpolation polynomial, the time delay term x(t-τ) is expressed as: x(t-τ) = a2x i-m +b2x i+1-m +c2x i+2-m (16); In formula (16), At this time, the total iteration of formula (5) is Q i x i+1 =H i x i +H i-m x i-m +H i-m+1 x i-m+1 +H i-m+2 x i-m+2 (20); where Q i =(I-(K1B i +K2B i+1 ))(twenty one); H i-m =K5B i +K6B i+1 (23); H i-m+1 =K7B i +K8B i+1 (24); H i-m+2 =K9B i +K 10 B i+1 (25), Q i , H i , H i-m , H i-m+1 , H i-m+2 is the time period matrix B i and B i+1 Related coefficient matrix; I represents the identity matrix; K1-K 10 is a constant matrix.

7. The method for predicting milling chatter stability according to claim 6, characterized in that: The Floquet transition matrix is ​​constructed based on the solution of the state space equation to determine the stability of the milling system. Specifically, the following discrete mapping is obtained according to formula (20): i+1 =D i y i (26); where y i =[y i ,y i-1 ,y i-2 …y i-m ] T (27); If the matrix Q i is non-singular, then D i Expressed as The state transition matrix Φ of the system in a single time period is obtained by the matrix sequence D i (i=1,2,…m), constructed, that is, y m =Φy0(29); where Φ = D m-1 D m-2 …D1D0(30); if Q i If it is a singular matrix, the generalized inverse of the matrix is ​​used instead of inverting it.

Images