Quick prediction method for variable-speed milling stability domain of variable-pitch variable-helix cutter

Abstract

The invention provides a rapid prediction method for a variable-rotating-speed milling stability domain of a variable-tooth-pitch variable-spiral cutter, and the method comprises the steps: equally dividing the axial direction of the variable-tooth-pitch variable-spiral cutter into a plurality of discrete intervals, calculating the overall time lag of the small intervals according to the tooth pitches at the middle height positions of the small intervals and the rotating speed at the current time, equally dividing the cutter tooth period into a plurality of discrete small intervals, and calculating the overall time lag of the small intervals; performing two-point numerical differential processing on the differential factor in each cell, performing left and right two-point linear interpolation on a state item and a periodic coefficient item thereof, performing left and right three-node value Hermite interpolation on a time delay item and a periodic coefficient item thereof, and finally performing arrangement to obtain an inter-cell transfer function. And according to the Floquet theory, a variable-rotating-speed milling stable lobe graph of the variable-tooth-pitch variable-spiral cutter is obtained. According to the method, the calculation efficiency can be effectively improved, so that the milling stability domain is predicted more efficiently, and a theoretical basis is provided for cutting technicians to select reasonable cutting parameters to improve the surface quality of parts.

Description

Technical Field

[0001] The invention belongs to the field of advanced manufacturing technology, and in particular relates to a method for quickly predicting a variable speed milling stability domain for a variable pitch and variable helix tool. Background Art

[0002] Milling is a highly efficient machining process with the advantages of high machining accuracy, high material removal rate and low machining cost. Therefore, it is widely used in the fields of aviation, aerospace, shipbuilding and automotive key parts processing. Chatter often exists in the milling process and is the main factor causing the instability of the machining process.

[0003] In addition, chatter in the milling process not only affects the surface quality of the machined parts, but also reduces the life of the CNC machine tool. Through the prediction of chatter stability in the milling process, a chatter-free milling parameter combination can be selected to avoid the occurrence of chatter, improve machining efficiency, and reduce machining costs. Variable pitch variable helix tool and spindle variable speed modulation can disturb the regeneration effect and suppress chatter by changing the phase of the dynamic cutting thickness during cutting. At present, the existing stability prediction methods (full discrete method, semi-discrete method, etc.) often select fewer time points within the cycle and their responses, and cannot take into account both the prediction accuracy and the calculation speed, which increases the calculation cost and time. In addition, the variable pitch variable helix tool and spindle variable speed modulation conditions further increase the amount of calculation. Therefore, the efficient stability prediction method for variable pitch variable helix tool and spindle variable speed modulation is of great significance to improving the high-performance machining characteristics of high-end CNC equipment.

[0004] At present, the literature "Nie WY, Zheng ML, Zhang W, et al. Analytical prediction of chatter stability with the effect of multiple delays for variable pitch endmills and optimization of pitch parameters[J]. International Journal of Advanced Manufacturing Technology, 2023, 124(7-8): 2645-2658." introduces a stability prediction method for variable pitch end mills using frequency domain and semi-discrete. This method has high calculation accuracy. However, since only the time delay term is discretized, the calculation efficiency of this method is relatively low. The literature "Sims ND, Mann B, Huyanan S. Analytical prediction of chatter stability for variable pitch and variable helix milling tools [J]. Journal of Sound and Vibration, 2008, 317 (3-5): 664-686." and the literature "Xie Qizhi. Research on Hopf bifurcation and semi-discrete algorithm of machine tool chatter [D]. Tianjin University, 2013." also use semi-discrete stability prediction methods for variable pitch and variable helix tools. In addition, the analytical modeling of the tools further reduces the computational efficiency of stability prediction of variable pitch and variable helix tools.

[0005] At present, the patent with application number "201510068454.7" requires the transfer matrix to be iterated m times at discrete intervals when obtaining the stability lobe diagram, which reduces the efficiency of calculation; the patent with application number "201811022089.6" proposes a milling chatter suppression method based on multi-frequency variable speed, but the description of the stability prediction method of multi-frequency variable speed is less. The patent with application number "202010784217.1" establishes a variable pitch milling cutter cutting model considering the helix angle effect. Its stability prediction method has high accuracy, but it does not consider the variable helix situation, and its calculation process is relatively complicated. The patent with application number "202010881385.2" uses the implicit exponential time difference multi-step method to solve the discrete values ​​at discrete points in each time interval, and the transfer matrix needs to be inverted when calculating the matrix, which has a large amount of calculation. Patent application number: "202110442947.8" solves the milling dynamics differential equation through the variable step-size Runge-Kutta algorithm. The step-size is changed according to the relationship between the given accuracy and the deviation size to ensure the calculation accuracy, but it greatly increases the amount of calculation. Patent application number: "202410632692.5" involves a large number of terms when performing discrete mapping, and the amount of calculation is large. Summary of the invention

[0006] In order to solve the problem of the variable speed milling stability domain prediction method for a variable pitch variable helix tool. The present invention provides a variable speed milling stability domain fast prediction method for a variable pitch variable helix tool. Compared with other methods for predicting milling stability, the method according to the embodiment of the present invention has higher computational efficiency.

[0007] The technical solution of the present invention is implemented as follows: a method for quickly predicting the stable region of variable speed milling for a variable pitch and variable helix tool, comprising:

[0008] Step 1: Divide the axial cutting height of the variable pitch variable helix tool into several small height intervals, and calculate the overall time lag of the small height intervals using the pitch and variable speed model at the middle height position of the small height intervals;

[0009] Step 2: Substitute the overall time delay of the height interval in step 1 into the milling dynamics delay differential equation considering the regeneration effect;

[0010] Step 3: Perform variable transformation on the time-delay differential equation to obtain the spatial state equation of the time-delay differential equation, and decompose the spatial state equation into three parts, including differential factors, state terms and time-delay terms;

[0011] Step 4: Divide the tooth period into several small time intervals, perform two-point numerical differentiation on the differential factors of the spatial state equation in the small time intervals, perform left and right two-point linear interpolation approximation on the state terms and periodic coefficient terms of the spatial state equation, and perform left and right three-point Hermitian interpolation approximation on the time lag terms and periodic coefficient terms of the spatial state equation;

[0012] Step 5: Substitute the approximate expressions of the differential factor, state term and time lag term in step 4 into the spatial state equation, sort out the recursive formula of the variables, and then transform the variables of the recursive formula to obtain the transfer function between the adjacent tooth time cells;

[0013] Step 6: Multiply the transfer functions of all adjacent tooth time cells in the tooth cycle to obtain the transfer function Fi of the tooth cycle;

[0014] Step 7: Calculate the characteristic value λ(Fi) of the transfer function of the tooth period, and judge the stability of the milling system according to the Floquet theory. The specific judgment is shown in Formula 1:

[0015]

[0016] Further, in step 1, obtaining the overall time lag between high-altitude cells includes:

[0017] First, the actual axial cutting height a of the variable pitch variable spiral cutter p Divide the tool into L equal height intervals, and then divide the tool again by the tool edge line. At this time, there are N teeth and L heights, a total of N*L tool intervals (N is the number of tool teeth), and the height of each height interval is δz=a p / L, the length of the tool cell is idealized as the length at the height of the tool cell δz / 2; time lag τ j,l (t) is the time required for tooth j-1 in the lth layer to rotate to tooth j at time t (j∈(0,N) and j is an integer, l∈(0,L) and l is an integer); the distance between tooth j-1 and tooth j in the lth layer is ψ j,l =p j +p l ,p j is the pitch from the jth tooth to the j-1th tooth (obtained from the tool parameters), p l is the change in tooth pitch caused by axial height at layer l, and:

[0018] p l =dz(l-0.5)(tanβ j-1 -tanβ j ) / L*360° (2)

[0019] In the formula, β jis the helix angle of the jth tooth, and the mathematical modeling of the variable speed of the spindle is:

[0020]

[0021] Where, Ω0 is the average spindle speed; Ω A is the modulation amplitude; P is the spindle speed modulation period; RVA is the modulation amplitude ratio; RVF is the modulation frequency ratio; RVA=Ω A / Ω0; RVF=60 / (pΩ0). Time lag τ j,l (t) is the time it takes for tooth j-1 in layer l to rotate to tooth j, and we have:

[0022]

[0023] In the formula, t0 is the time for the lth layer to rotate to tooth j. Substituting formula (3) into formula (4) and integrating it, we get:

[0024]

[0025] If Ω A Relative to Ω0, τ j,l (t) will be approximately

[0026]

[0027] In the formula,

[0028] Furthermore, in step 2, the overall time delay between tool cells is introduced into the milling dynamics delay differential equation considering the regeneration effect:

[0029] φ j,l (t) is the angular position at time t, jth tooth, and lth layer of height cell, and we have:

[0030]

[0031] According to the regeneration effect, at time t and angular position φ j,l (t) cutting thickness h(φ j,l (t)) has:

[0032] h(φ j,l (t)) = [x(t)-x(t-τ j,l (t))]sin(φ j,l (t))+[y(t)-y(t-τ j,l (t))]cos(φ j,l (t)) (8)

[0033] In the formula, x(t) is the function of the displacement in the x-degree of freedom direction and the time relationship, and y(t) is the function of the displacement in the y-degree of freedom direction and the time relationship. After being substituted into formula (8) and rearranged, the time-delay differential equations of the two degrees of freedom of x and y are:

[0034]

[0035] In the formula, is the first-order derivative of x(t) with respect to time t, is the second-order derivative of x(t) with respect to time t, is the first-order derivative of y(t) with respect to time t, is the second-order derivative of y(t) with respect to time t; M, C and K are the mass matrix, damping matrix and stiffness matrix corresponding to the milling system in the x and y directions respectively. If modal coupling is not considered, then:

[0036]

[0037] In the formula, m xx and m yy are the modal masses in the x and y directions respectively; c xx and c yy are the modal damping in the x and y directions respectively; k xx and k yy are the modal stiffness in the x and y directions respectively; K t and K n are the tangential and radial cutting force coefficients respectively; g(φ j,l ) represents the angular position φ j,l The coefficient restriction function of whether the tool cell is involved in cutting is 1 if it is involved in cutting, otherwise it is 0;

[0038] Furthermore, in step 3, the spatial state equation of the delay differential equation is derived:

[0039] Let: p(t) = [x(t) y(t)] T , substituting p(t) into equation (9) yields:

[0040]

[0041] In the formula, is the first-order derivative of p(t) with respect to time t, is the second-order derivative of p(t) with respect to time t;

[0042] make: The first derivative of q(t) with respect to time t have:

[0043]

[0044]

[0045] Substituting equation (12) and equation (13) into equation (11), we get:

[0046]

[0047] In the formula,

[0048]

[0049]

[0050]

[0051] is the differential factor, q(t) is the state term, A j,l (t) is the periodic coefficient term of q(t), q(t-τ j,l (t)) is the time lag term, B j,l (t) is q(t-τ j,l (t)) periodic coefficient term.

[0052] Furthermore, in step 4, the tooth period T is equally divided into several time intervals, and the various terms of the spatial state equation are approximated:

[0053] The tooth period T is

[0054]

[0055] Divide the tooth period T into k equal parts, and the time step Δt is:

[0056]

[0057] For the i-th hour interval (t i ,t i+1 ), T = kΔt, t i =iΔt, let q i =q(t i ), i∈(0,k-1) and i is an integer, then the i-th hour interval (t i ,t i+1 ) within the average time lag for:

[0058]

[0059] make:

[0060]

[0061] In the formula, int() is a rounding function that tends to 0, m i,j,lThe average time lag The time step difference of the corresponding small time interval is the i-th hour time interval of the l-th height small interval of the j-th tooth (t i ,t i+1 ) corresponds to the time series with a time lag, then the im i,j,l Hour interval is the time interval of the last tooth action. Take m i,j,l The maximum value of m max .

[0062] Differentiation Factor Discretization is performed by two-point numerical differentiation (t∈(t i ,t i+1 ),δ=t i -t):

[0063]

[0064] For the state item q(t), pass through its left and right nodes q i and q i+1 The linear difference is (q i =q(t i )):

[0065]

[0066] For the periodic coefficient term A of the state term q(t) j,l (t), through its left and right node values ​​A j,l (t i ) and A j,l (t i+1 ) Using first-order linear interpolation:

[0067]

[0068] For the delay term q(t-τ j,l (t)) through its left and right three node values and Perform Hermitian interpolation, for the delay term q(t-τ j,l (t)) of the periodic coefficient term B j,l (t) through its left and right three node values ​​B j,l (t i ), B j,l (t i+1 ) and B j,l (t i+2 ) performs Hermitian interpolation:

[0069]

[0070] B j,l(t) = a1B j,l (t i )+b1B j,l (t i+1 )+c1B j,l (t i+2 ) (26)

[0071] Among them, a1, b1 and c1 are difference coefficients:

[0072]

[0073] Furthermore, in step 5, the transfer function D between adjacent blade tooth time cells is i :

[0074] Substituting equations (22), (23), (24), (25) and (26) into equation (14) and rearranging them, we obtain:

[0075]

[0076] In the formula,

[0077]

[0078] P i,j,l =a1B j,l (t i )+b1B j,l (t i+1 )+c1B j,l (t i+2 ) (33)

[0079] If the above H i Existence, order And put it into formula (30):

[0080] Z i+1 =D i Z i (34)

[0081] In the formula,

[0082]

[0083] Furthermore, in step 6, the transfer function Fi of the adjacent tooth period is:

[0084]

[0085] Further, in step 7, according to Floquet theory, the system stability is given by the transfer matrix Fi The maximum value μ of the modulus of the eigenvalue is determined by:

[0086]

[0087] Beneficial effects:

[0088] In an embodiment of the present invention, a time-delay differential equation for the milling dynamics of a variable pitch and variable spiral tool and a variable speed is established through steps 1 and 2, and the time-delay differential equation is converted into a spatial state equation through step 3; the tooth period is equally divided into a number of time intervals through step 4, and the differential factor of the spatial state equation is subjected to two-point numerical differentiation processing within the time interval, the integration step is removed, the state term and the periodic coefficient term of the spatial state equation are subjected to left and right two-point linear interpolation approximation, and the time-delay term and the periodic coefficient term of the spatial state equation are subjected to left and right three-point Hermite interpolation approximation, thereby reducing the fitting error of the calculation method and the amount of interpolation calculation, and significantly improving the calculation efficiency (>=30%) at the expense of a small amount of calculation accuracy (<=5%). BRIEF DESCRIPTION OF THE DRAWINGS

[0089] Figure 1 The stability lobe diagram of the 4-tooth variable pitch variable spiral cutter of the present invention

[0090] Figure 2 This is the convergence diagram of the present invention with an axial cutting depth of 4 mm and a spindle speed of 5760 rpm

[0091] Figure 3 It is the overall flow chart of the present invention. DETAILED DESCRIPTION

[0092] 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.

[0093] like Figure 3 As shown, an embodiment of the present invention provides a method for quickly predicting the stable region of variable speed milling for a variable pitch and variable helix tool, and the method may include the following steps:

[0094] Step 1: Divide the axial cutting height of the variable pitch variable helix tool into several small height intervals, and calculate the overall time lag of the small height intervals using the pitch and variable speed model at the middle height position of the small height intervals;

[0095] Step 2: Substitute the overall time delay of the height interval in step 1 into the milling dynamics delay differential equation considering the regeneration effect;

[0096] Step 3: Perform variable transformation on the time-delay differential equation to obtain the spatial state equation of the time-delay differential equation, and decompose the spatial state equation into three parts, including differential factors, state terms and time-delay terms;

[0097] Step 4: Divide the tooth period into several small time intervals, perform two-point numerical differentiation on the differential factors of the spatial state equation in the small time intervals, perform left and right two-point linear interpolation approximation on the state terms and periodic coefficient terms of the spatial state equation, and perform left and right three-point Hermitian interpolation approximation on the time lag terms and periodic coefficient terms of the spatial state equation;

[0098] Step 5: Substitute the approximate expressions of the differential factor, state term and time lag term in step 4 into the spatial state equation, sort out the recursive formula of the variables, and then transform the variables of the recursive formula to obtain the transfer function between the adjacent tooth time cells;

[0099] Step 6: Multiply the transfer functions of all adjacent tooth time cells in the tooth cycle to obtain the transfer function Fi of the tooth cycle;

[0100] Step 7: Calculate the characteristic value λ(Fi) of the transfer function of the tooth period, and judge the stability of the milling system according to the Floquet theory. The specific judgment is shown in Formula 1:

[0101]

[0102] Further, in step 1, obtaining the overall time lag between high-altitude cells includes:

[0103] First, the actual axial cutting height a of the variable pitch variable spiral cutter p Divide the tool into L equal height intervals, and then divide the tool again by the tool edge line. At this time, there are N teeth and L heights, a total of N*L tool intervals (N is the number of tool teeth), and the height of each height interval is δz=a p / L, the length of the tool cell is idealized as the length at the height of the tool cell δz / 2; time lag τ j,l (t) is the time required for tooth j-1 in the lth layer to rotate to tooth j at time t (j∈(0,N) and j is an integer, l∈(0,L) and l is an integer); the distance between tooth j-1 and tooth j in the lth layer is ψ j,l =p j +p l ,p j is the pitch from the jth tooth to the j-1th tooth (obtained from the tool parameters), p l is the change in tooth pitch caused by axial height at layer l, and:

[0104] p l =dz(l-0.5)(tanβj-1 -tanβ j ) / L*360° (39)

[0105] In the formula, β j is the helix angle of the jth tooth, and the mathematical modeling of the variable speed of the spindle is:

[0106]

[0107] Where, Ω0 is the average spindle speed; Ω A is the modulation amplitude; P is the spindle speed modulation period; RVA is the modulation amplitude ratio; RVF is the modulation frequency ratio; RVA=Ω A / Ω0; RVF=60 / (pΩ0). Time lag τ j,l (t) is the time it takes for tooth j-1 in layer l to rotate to tooth j, and we have:

[0108]

[0109] In the formula, t0 is the time for the lth layer to rotate to tooth j. Substituting formula (3) into formula (4) and integrating it, we get:

[0110]

[0111] If Ω A Relative to Ω0, τ j,l (t) will be approximately

[0112]

[0113] In the formula,

[0114] Furthermore, in step 2, the overall time delay between tool cells is introduced into the milling dynamics delay differential equation considering the regeneration effect:

[0115] The present invention mainly constructs the x and y two-degree-of-freedom time-delay differential equations.

[0116] φ j,l (t) is the angular position at time t, jth tooth, and lth layer of height interval, and we have:

[0117]

[0118] According to the regeneration effect, at time t and angular position φ j,l (t) cutting thickness h(φ j,l (t)) has:

[0119] h(φ j,l (t)) = [x(t)-x(t-τ j,l (t))]sin(φj,l (t))+[y(t)-y(t-τ j,l (t))]cos(φ j,l (t)) (45)

[0120] In the formula, x(t) is the function of the displacement in the x-degree of freedom direction and the time relationship, and y(t) is the function of the displacement in the y-degree of freedom direction and the time relationship. After being substituted into formula (8) and rearranged, the time-delay differential equations of the two degrees of freedom of x and y are:

[0121]

[0122] In the formula, is the first-order derivative of x(t) with respect to time t, is the second-order derivative of x(t) with respect to time t, is the first-order derivative of y(t) with respect to time t, is the second-order derivative of y(t) with respect to time t; M, C and K are the mass matrix, damping matrix and stiffness matrix corresponding to the milling system in the x and y directions respectively. If modal coupling is not considered, then:

[0123]

[0124] In the formula, m xx and m yy are the modal masses in the x and y directions respectively; c xx and c yy are the modal damping in the x and y directions respectively; k xx and k yy are the modal stiffness in the x and y directions respectively; K t and K n are the tangential and radial cutting force coefficients respectively; g(φ j,l ) represents the angular position φ j,l The coefficient restriction function of whether the tool cell is involved in cutting is 1 if it is involved in cutting, otherwise it is 0;

[0125] Furthermore, in step 3, the spatial state equation of the delay differential equation is derived:

[0126] Let: p(t) = [x(t) y(t)] T , substituting p(t) into equation (9) yields:

[0127]

[0128] In the formula, is the first-order derivative of p(t) with respect to time t, is the second-order derivative of p(t) with respect to time t;

[0129] make: The first derivative of q(t) with respect to time t have:

[0130]

[0131]

[0132] Substituting equation (12) and equation (13) into equation (11), we get:

[0133]

[0134] In the formula,

[0135]

[0136]

[0137] is the differential factor, q(t) is the state term, A j,l (t) is the periodic coefficient term of q(t), q(t-τ j,l (t)) is the time lag term, B j,l (t) is q(t-τ j,l (t)) periodic coefficient term.

[0138] Furthermore, in step 4, the tooth period T is equally divided into several time intervals, and the various terms of the spatial state equation are approximated:

[0139] The tooth period T is

[0140]

[0141] Divide the tooth period T into k equal parts, and the time step Δt is:

[0142]

[0143] For the i-th hour interval (t i ,t i+1 ), T = kΔt, t i =iΔt, let q i =q(t i ), i∈(0,k-1) and i is an integer, then the i-th hour interval (t i ,t i+1 ) within the average time lag for:

[0144]

[0145] make:

[0146]

[0147] In the formula, int() is a rounding function that tends to 0, m i,j,l The average time lag The time step difference of the corresponding small time interval is the i-th hour time interval of the l-th height small interval of the j-th tooth (t i ,t i+1 ) corresponds to the time series with a time lag, then the im i,j,l Hour interval is the time interval of the last tooth action. Take m i,j,l The maximum value of m max .

[0148] Differentiation Factor Discretization is performed by two-point numerical differentiation (t∈(t i ,t i+1 ),δ=t i -t):

[0149]

[0150] For the state item q(t), pass through its left and right nodes q i and q i+1 The linear difference is (q i =q(t i )):

[0151]

[0152] For the periodic coefficient term A of the state term q(t) j,l (t), through its left and right node values ​​A j,l (t i ) and A j,l (t i+1 ) Using first-order linear interpolation:

[0153]

[0154] For the delay term q(t-τ j,l (t)) through its left and right three node values and Perform Hermitian interpolation, for the delay term q(t-τ j,l (t)) of the periodic coefficient term B j,l (t) through its left and right three node values ​​B j,l (t i ), B j,l (t i+1 ) and B j,l (t i+2 ) performs Hermitian interpolation:

[0155]

[0156] B j,l (t) = a1B j,l (t i )+b1B j,l (t i+1 )+c1B j,l (t i+2 ) (63)

[0157] Among them, a1, b1 and c1 are difference coefficients:

[0158]

[0159] Furthermore, in step 5, the transfer function D between adjacent blade tooth time cells is i :

[0160] Substituting equations (22), (23), (24), (25) and (26) into equation (14) and rearranging them, we obtain:

[0161]

[0162] In the formula,

[0163]

[0164] P i,j,l =a1B j,l (t i )+b1B j,l (t i+1 )+c1B j,l (t i+2 ) (70)

[0165] If the above H i Existence, order And put it into formula (30):

[0166] Z i+1 =D i Z i (71)

[0167] In the formula,

[0168]

[0169] Furthermore, in step 6, the transfer function Fi of the adjacent tooth period is:

[0170]

[0171] Furthermore, in step 7, according to Floquet theory, the system stability is determined by the maximum value μ of the modulus of the eigenvalue of the transfer matrix Fi:

[0172]

[0173] This embodiment is programmed using Matlab software, using a vertical milling cutter with a 4-tooth helix of [30° 35° 30° 35°] and a pitch of [70° 110° 70° 110°], a tooth period discrete number k=144, an axial cutting height discrete number L=5, and an axial cutting depth a p The discrete number of the average spindle speed Ω0 is 50X50, and the variable speed parameters RVA=0.1, RVF=0.5, Figure 1 The stability lobe diagram under this condition. Axial cutting depth 4mm, spindle speed 5760rmp, μ precise value u0=0.418, μ calculated value u, Figure 2 This is the convergence diagram of |u-u0| changing with the number of discrete parts of the tooth period k.

[0174] The above description is only a preferred embodiment of the present invention and is not intended to limit the present invention. Any modifications, equivalent substitutions, improvements, etc. made within the spirit and principle of the present invention should be included in the protection scope of the present invention.

Claims

1. A method for quickly predicting the stable region of variable speed milling for a variable pitch and variable helix tool, characterized in that: The specific steps include: Step 1: Divide the axial cutting height of the variable pitch variable helix tool into several small height intervals, and calculate the overall time lag of the small height intervals using the pitch and variable speed model at the middle height position of the small height intervals; Step 2: Substitute the overall time delay of the height interval in step 1 into the milling dynamics delay differential equation considering the regeneration effect; Step 3: Perform variable transformation on the time-delay differential equation to obtain the spatial state equation of the time-delay differential equation, and decompose the spatial state equation into three parts, including differential factors, state terms and time-delay terms; Step 4: Divide the tooth period into several small time intervals, perform two-point numerical differentiation on the differential factors of the spatial state equation in the small time intervals, perform left and right two-point linear interpolation approximation on the state terms and periodic coefficient terms of the spatial state equation, and perform left and right three-point Hermitian interpolation approximation on the time lag terms and periodic coefficient terms of the spatial state equation; Step 5: Substitute the approximate expressions of the differential factor, state term and time lag term in step 4 into the spatial state equation, sort out the recursive formula of the variables, and then transform the variables of the recursive formula to obtain the transfer function between the adjacent tooth time cells; Step 6: Multiply the transfer functions of all adjacent tooth time cells in the tooth cycle to get the transfer function of the tooth cycle. Fi ; Step 7: Calculate the characteristic value λ(Fi) of the transfer function of the tooth period, and judge the stability of the milling system according to the Floquet theory. The specific judgment is shown in Formula 1:

2. The method according to claim 1, characterized in that The obtaining of the overall time lag between the height cells includes: First, the actual axial cutting height a of the variable pitch variable spiral cutter p Divide the tool into L equal height intervals, and then divide the tool again by the tool edge line. At this time, there are N teeth and L heights, a total of N*L tool intervals (N is the number of tool teeth), and the height of each height interval is δz=a p / L, the length of the tool cell is idealized as the length at the height of the tool cell δz / 2; time lag τ j,l (t) is the time required for tooth j-1 in the lth layer to rotate to tooth j at time t (j∈(0,N) and j is an integer, l∈(0,L) and l is an integer); the distance between tooth j-1 and tooth j in the lth layer is ψ j,l =p j +p l ,p j is the pitch from the jth tooth to the j-1th tooth (obtained from the tool parameters), p l is the change in tooth pitch caused by axial height at layer 1, and: <h2 style=";text-align:left;direction:ltr">p<h2 style=";text-align:left;direction:ltr"> l <h2 style=";text-align:left;direction:ltr"> =dz(l-0.5)(tanβ<h2 style=";text-align:left;direction:ltr"> j-1 <h2 style=";text-align:left;direction:ltr"> -tanβ<h2 style=";text-align:left;direction:ltr"> j <h2 style=";text-align:left;direction:ltr"> ) / L*360° (2) In the formula, β j is the helix angle of tooth j, and the mathematical modeling of the variable speed of the spindle is: Where, Ω0 is the average spindle speed; Ω A is the modulation amplitude; P is the spindle speed modulation period; RVA is the modulation amplitude ratio; RVF is the modulation frequency ratio; RVA=Ω A / Ω0; RVF = 60 / (pΩ0). Due to the time lag τ j,l (t) is the time it takes for tooth j-1 in layer l to rotate to tooth j at time t, then: Where t0 is the time lag τ j,l (t) corresponds to the time when the lth layer rotates to tooth j. Substituting equation (3) into equation (4) and integrating it, we get: If Ω A Relative to Ω0, τ j,l (t) will be approximately In the formula, 3. The method according to claim 2, characterized in that The overall time delay between tool cells is introduced into the milling dynamics delay differential equation considering the regeneration effect: The present invention mainly constructs the x and y two-degree-of-freedom time-delay differential equations. φ j,l (t) is the angular position at time t, jth tooth, and lth layer of height cell, and we have: According to the regeneration effect, at time t and angular position φ j,l (t) cutting thickness h(φ j,l (t)) has: h(φ j,l (t))=[x(t)-x(t-τ j,l (t))]sin(φ j,l (t))+[y(t)-y(t-τ j,l (t))]cos(φ j,l (t))(8) In the formula, x(t) is the function of the displacement in the x-degree of freedom direction and the time relationship, and y(t) is the function of the displacement in the y-degree of freedom direction and the time relationship. After being substituted into formula (8) and rearranged, the time-delay differential equations of the two degrees of freedom of x and y are: In the formula, is the first-order derivative of x(t) with respect to time t, is the second-order derivative of x(t) with respect to time t, is the first-order derivative of y(t) with respect to time t, is the second-order derivative of y(t) with respect to time t; M, C and K are the mass matrix, damping matrix and stiffness matrix of the milling system in the x and y directions respectively. If the modal coupling is not considered, then: In the formula, m xx and m yy are the modal masses in the x and y directions respectively; c xx and c yy are the modal damping in the x and y directions respectively; k xx and k yy are the modal stiffness in the x and y directions respectively; K t and K n are the tangential and radial cutting force coefficients respectively; g(φ j,l ) represents the angular position φ j,l The coefficient limit function of whether the tool cell participates in cutting is 1 if it participates in cutting, otherwise it is 0.

4. The method according to claim 3, characterized in that The spatial state equation of the delay differential equation: Let: p(t) = [x(t) y(t)] T , substituting p(t) into equation (9) yields: In the formula, is the first-order derivative of p(t) with respect to time t, is the second-order derivative of p(t) with respect to time t; make: The first derivative of q(t) with respect to time t have: Substituting equation (12) and equation (13) into equation (11), we get: In the formula, is the differential factor, q(t) is the state term, A j,l (t) is the periodic coefficient term of q(t), q(t-τ j,l (t)) is the time lag term, B j,l (t) is q(t-τ j,l (t)) periodic coefficient term.

5. The method according to claim 4, characterized in that Divide the tooth period T into several small time intervals and approximate each term of the spatial state equation: The tooth period T is Divide the tooth period T into k equal parts, and the time step Δt is: For the i-th hour interval (t i ,t i+1 ), T = kΔt, t i =iΔt, let q i =q(t i ), i∈(0,k-1) and i is an integer, then the i-th hour interval (t i ,t i+1 ) within the average time lag for: make: In the formula, int() is a rounding function that tends to 0, m i,j,l The average time lag The time step difference of the corresponding small time interval is the i-th hour time interval of the l-th height small interval of the j-th tooth (t i ,t i+1 ) corresponds to the time series with a time lag, then the im i,j,l Hour interval is the time interval of the last tooth action. Take m i,j,l The maximum value of m max . Differentiation Factor Discretization is performed by two-point numerical differentiation (t∈(t i ,t i+1 ),δ=t i -t): For the state item q(t), pass through its left and right nodes q i and q i+1 The linear difference is (q i =q(t i )): For the periodic coefficient term A of the state term q(t) j,l (t), through its left and right node values ​​A j,l (t i ) and A j,l (t i+1 ) Using first-order linear interpolation: For the delay term q(t-τ j,l (t)) through its left and right three node values and Perform Hermitian interpolation, for the delay term q(t-τ j,l (t)) of the periodic coefficient term B j,l (t) through its left and right three node values ​​B j,l (t i ), B j,l (t i+1 ) and B j,l (t i+2 ) performs Hermitian interpolation: B j,l (t)=a1B j,l (t i )+b1B j,l (t i+1 )+c1B j,l (t i+2 ) (26) Among them, a1, b1 and c1 are difference coefficients:

6. The method according to claim 5, characterized in that The transfer function Fi of adjacent tooth periods: Substituting equations (22), (23), (24), (25) and (26) into equation (14) and rearranging them, we obtain: In the formula, P i,j,l =a1B j,l (t i )+b1B j,l (t i+1 )+c1B j,l (t i+2 ) (33) If the above H i Existence, order And put it into formula (30) to get: Z i+1 =D i Z i (34) In the formula, The transfer function Fi of adjacent tooth periods: According to Floquet theory, the stability of the system is determined by the maximum value μ of the modulus of the eigenvalue of the transfer matrix Fi:

Images