Stability of finish milling and surface position error comprehensive prediction method and device

Abstract

The application discloses a stability and surface position error comprehensive prediction method and equipment for milling finishing, and belongs to the technical field of milling finishing in mechanical processing. The method comprises the following steps: acquiring milling cutter information; calculating the cutting edge angle position, cutting radius, feed amount and total cutting force in X and Y directions at any time; constructing milling dynamics time-delay differential equation groups and transforming into space state equations; performing integral calculation on the space state equations to determine the spindle speed range, corresponding discrete numbers, axial cutting depth range, corresponding discrete numbers and single tooth cutting period discrete numbers; calculating state items, time-delay items and static force items in the space state equations, constructing the state transmission relationship of the milling finishing system in a single time period, and determining the stability of the milling finishing system; extracting corresponding steady-state coefficient vectors Y direction coefficients to predict dynamic machining errors; and obtaining stability petal diagrams and surface position error prediction diagrams. The application has good accuracy and timeliness.

Description

Technical Field

[0001] This invention relates to the field of milling finishing technology in machining, and more specifically to a method and equipment for comprehensively predicting the stability and surface position error of milling finishing. Background Technology

[0002] Milling technology is widely used in the manufacturing of complex parts in aerospace, shipbuilding, and automotive industries, and is one of the essential technologies for achieving high-precision machining quality. In actual machining processes, chatter, caused by improper selection of machining parameters, is a major reason for reduced surface quality and limited production efficiency. Therefore, stability analysis of the milling process is a crucial prerequisite for achieving chatter-free milling. However, due to the presence of forced vibration, chatter-free machining parameters are still insufficient to guarantee the desired machining accuracy. Especially in finishing processes, priority should be given to the surface quality of the machined parts, and maximum machining efficiency should be ensured only after meeting surface quality requirements.

[0003] Currently, some research has been conducted by relevant personnel in this field. For example, the literature "Li ZY, Jiang SL, Sun YW. Chatter stability and surface location error predictions in milling with mode coupling and process damping[J]. Proceedings of the Institution of Mechanical Engineers, Part B: Journal of Engineering Manufacture, 2019, 233(3): 686-698." discloses a milling dynamics modeling method that considers modal coupling and process damping, and uses a second-order semi-discrete method to simultaneously predict milling stability and surface location error. However, when considering factors such as tool runout, it cannot accurately predict the stability boundary and surface location error, and the prediction accuracy is limited. Another example is the literature "Wan M, Zhang WH, Dang JW, et al. A unified stability prediction method for milling process with multiple delays[J]. International Journal of Machine Tools and Manufacture, 2010, 50(1):29-41. A multi-delay milling stability prediction model considering tool runout is disclosed. It is only used to predict milling stability. Moreover, the method does not consider the modal coupling effect and process damping of the milling process, and the prediction result has a large error.

[0004] The patent with publication number "CN106808320A" discloses a milling force prediction method that considers the tool tooth angle, helix angle and eccentricity, which is only used to predict cutting force; the patent with publication number "CN108647413A" discloses a comprehensive prediction method for micro-surface position error and stability, which does not consider the modal coupling effect, process damping and tool runout during the milling process, and the prediction error is larger than the actual result.

[0005] In summary, most existing milling stability and surface position error prediction models only consider the regenerative effect during the milling process. However, in reality, the milling process is often accompanied by the influence of a variety of complex factors. Numerous studies have shown that there are significant discrepancies between the stability and surface position error predictions based on traditional dynamic models and experimental results. These discrepancies are largely due to the neglect of modal coupling effects, process damping, and tool runout.

[0006] Therefore, in view of the above-mentioned defects or improvement needs of the existing technology, how to provide a comprehensive prediction method and equipment for stability and surface position error of milling finishing is a problem that urgently needs to be solved by those skilled in the art. Summary of the Invention

[0007] In view of this, the present invention provides a method and device for comprehensively predicting the stability and surface position error of milling finishing, so as to solve the technical problems existing in the prior art.

[0008] To achieve the above objectives, the present invention provides the following technical solution:

[0009] On one hand, Embodiment 1 of the present invention discloses a method for comprehensively predicting the stability and surface position error of milling finishing, comprising the following steps:

[0010] Step 1: Obtain milling cutter information, which includes: tool geometry parameters, modal parameters, cutting force coefficients, and cutting conditions;

[0011] Step 2: Calculate the cutting edge angle position, cutting radius, feed rate, and total cutting force in the X and Y directions at any given time based on the milling cutter information;

[0012] Step 3: Construct a system of differential equations for milling dynamics and transform it into a spatial state equation;

[0013] Step 4: Divide the cutting cycle T of the cutting teeth into equal parts.

[0014] The cutting cycle of the cutting tooth is divided into several small intervals. The spatial state equation is integrated within any small interval to determine the spindle speed range and its corresponding discrete number, the axial depth of cut range and its corresponding discrete number, and the discrete number of a single cutting tooth cycle.

[0015] Step 5: Use third-order Newton's equation, Hermite's equation, and linear interpolation polynomial to approximate the state terms, time delay terms, and static force terms in the spatial state equation, construct the state transfer relationship of the milling and finishing system in a single time period, and determine the stability of the milling and finishing system;

[0016] Step 6: Extract the corresponding steady-state coefficient vector Y-direction coefficient to predict dynamic machining error;

[0017] Step 7: Based on Steps 5 and 6, obtain the stability lobe diagram and the surface position error prediction diagram.

[0018] Preferably, step 2, which calculates the cutting edge angle position, cutting radius, feed rate, and total cutting force in the X and Y directions at any given time based on the milling cutter information, comprises the following steps:

[0019] Step 2.1: Calculate the tool angle position of the cutting edge at any given time;

[0020] Step 2.2: Calculate the actual cutting radius and feed per tooth of the cutting edge under the condition of tool eccentricity;

[0021] Step 2.3: Calculate the cutting forces in the tangential and radial directions of the tool, and convert the cutting forces in the tangential and radial directions of the tool to the X and Y axes of the machine tool coordinate system through Cartesian coordinate transformation;

[0022] Step 2.4: Convert the cutting forces in the tangential and radial directions of the tool to the X and Y axes of the machine tool coordinate system through Cartesian coordinate system transformation. Calculate the total cutting forces acting on the milling cutter in the X and Y directions based on the cumulative cutting forces acting on the cutter teeth.

[0023] Preferably, the tool angle position of the cutting edge at any given time is calculated using the following formula:

[0024]

[0025] In the formula, Ω is the spindle speed; N is the number of cutter teeth; and j is the cutter tooth number.

[0026] Preferably, the specific steps for calculating the actual cutting radius and feed per tooth of the cutting edge under the condition of tool eccentricity are as follows:

[0027] The actual cutting radius of each tooth is:

[0028]

[0029]

[0030]

[0031] The actual feed per tooth is:

[0032]

[0033]

[0034]

[0035] In the formula, f t a is the nominal feed per tooth. r Radial cutting depth;

[0036] Based on geometric relationships, the angle between two adjacent cutting teeth can be obtained using the law of cosines:

[0037]

[0038]

[0039]

[0040] The entry and exit angles of the j-th tooth are expressed as:

[0041]

[0042]

[0043] The above provides the actual cutting radius and feed per tooth of the cutting edge under eccentric conditions.

[0044] Calculate the cutting forces in the tangential and radial directions of the tool, and then transform these cutting forces to the X and Y axes of the machine tool coordinate system using Cartesian coordinate transformation. The specific steps are as follows:

[0045] Calculation of cutting forces in the tangential and radial directions:

[0046] The milling forces acting on the cutting tooth j mainly include shearing force and process damping force, which are expressed as follows:

[0047] F t =F ts +F tp ,F r =F rs +F rp

[0048] In the formula, F ts and F rs F represents the tangential and normal forces caused by the shear effect. tp and F rp These are the tangential and normal forces caused by the process damping effect. Based on the shear force model, F ts and F rs They are respectively represented as

[0049] F ts =g(φ j (t))K t a p h(t),F rs =g(φ j (t))K n a p h(t)

[0050] In the formula, K t and K n These are the tangential and radial cutting force coefficients, respectively; ap is the axial cutting depth; and the window function g(φ) is... j (t) is defined as

[0051]

[0052] In the formula, φst and φ ex These are the entry and exit angles of the cutting teeth, respectively. Based on regenerative chatter theory, the instantaneous cutting thickness can be expressed as:

[0053] h(t) = [f t +x(tT)-x(t)]sin(φ j (t))+[y(tT)-y(t)]cos(φ j (t))

[0054] In the formula, f t T is the feed per tooth, and T is the time delay, which is equal to the cutting cycle of the cutting tooth, i.e., T = 60 / (NΩ).

[0055] Based on the process damping mechanism, the radial process damping force can be expressed as:

[0056] F rp =g(φ j (t))K sp a p S

[0057] In the formula, K sp Let S and S represent the indentation coefficient and area, respectively. The tangential force generated by the process damping effect can be obtained from the following equation:

[0058] F tp =μF rp

[0059] In the formula, μ is the coefficient of friction. The equivalent viscous damping can be expressed as:

[0060]

[0061] In the formula, v is the cutting speed, and W is the width of the tool wear band.

[0062] By transforming the Cartesian coordinate system, the cutting forces in the tangential and radial directions of the cutting tool are converted to the X and Y axes of the machine tool coordinate system. The total cutting force acting on the milling cutter in the X and Y directions is calculated based on the cumulative cutting forces acting on the cutting teeth.

[0063] Preferably, the specific method for constructing the milling dynamics time-delay differential equation system in step 3 and transforming it into a spatial state equation is as follows:

[0064] Step 3.1: Construct the set of time-delayed differential equations for milling dynamics:

[0065]

[0066] In the formula,

[0067]

[0068] Step 3.2: Let Using Cauchy transformation, the equations are converted into spatial state form, resulting in the following spatial state equation:

[0069]

[0070] In the formula,

[0071]

[0072] In the formula, M is the modal mass matrix, m xx m yy For the mass in the X and Y directions, m xy m yx The coupling mass in the X and Y directions in this coordinate system is given. C is the modal damping matrix, c xx c yy For damping in the X and Y directions, c xy c yx The X and Y axes represent the coupling damping in this coordinate system. K is the modal stiffness matrix, k xx k yy Let k be the stiffness in the X and Y directions. xy k yx C represents the coupling stiffness in the X and Y directions in this coordinate system. eq For equivalent viscous damping, K sp Let a be the compression factor. p denoted as C, where W is the axial depth of cut, W is the tool wear bandwidth, and v is the cutting speed. p H is the conversion coefficient matrix of cutting force, and we have:

[0073]

[0074]

[0075]

[0076]

[0077] f0(t) is the static force vector, f t This is the feed per tooth.

[0078] Preferably, the specific method for dividing the cutting cycle of the cutting teeth into several small intervals in step 4, and performing integral calculation of the spatial state equation within any small interval, is as follows:

[0079] Define κ(t) = A(t)x(t) and α(tT) = B(t)x(tT);

[0080] With x(t0) as the initial condition, the general solution of the space state equation is expressed as:

[0081]

[0082] Dividing the cutting cycle T of the cutting teeth into m equal time intervals, the time step τ = T / m; the above formula can be expressed as follows for the time interval kτ≤t≤(k+1)τ, (k=1,2,3,…m):

[0083]

[0084] In the formula, δ=ξ-kτ,δ∈[0,τ]; in the formula, x((k+1)τ) is denoted as x k+1 .

[0085] Dividing the cutting cycle T of the cutting teeth into m equal time intervals, the time step τ = T / m; the above formula can be expressed as follows for the time interval kτ≤t≤(k+1)τ, (k=1,2,3,…m):

[0086]

[0087] In the formula, δ=ξ-kτ,δ∈[0,τ]; in the formula, x((k+1)τ) is denoted as x k+1 .

[0088] The maximum permissible cutting speed is determined based on several major categories of constraints, including machine tool constraints (such as maximum permissible cutting power and maximum permissible cutting force), tool constraints (such as tool material and tool life), and part constraints (such as the accuracy and material of the workpiece). Then, the minimum permissible spindle speed n that the milling machining system can withstand is determined according to the maximum permissible cutting speed. min Maximum spindle speed n max Thus, the spindle speed range [n] can be obtained. min ,n max Discrete number k of the spindle speed range. n It is usually set to 200.

[0089] The maximum permissible axial depth of cut is determined based on machine tool rigidity, tool diameter, and machining accuracy requirements. Then, considering that a larger axial depth of cut is typically used during roughing to achieve a higher material removal rate and a smaller axial depth of cut is typically used during finishing to achieve higher machining accuracy, the acceptable minimum permissible spindle speed n is determined. min Maximum spindle speed n max Thus, the axial cutting depth range [w] can be obtained. min ,w max The discrete number k of the axial depth range. w It is usually set to 100.

[0090] The discrete number m of the cutting cycle T of the cutting tooth is usually set to 10 levels, such as 20, 30, 40, 50, 60, 70, 80, 90, 100, and 200.

[0091] In the process of solving chatter stability, the discrete number of spindle speed range, the discrete number of axial depth of cut range, and the discrete number of single tooth cutting cycle can be reduced or increased depending on the actual needs of calculation accuracy and efficiency.

[0092] Preferably, the specific steps for step 5—calculating the state terms, time delay terms, and static force terms in the spatial state equation, constructing the state transfer relationship of the milling and finishing system over a single time period, and determining the stability of the milling and finishing system—are as follows:

[0093] Step 5.1: Use third-order Newtonian, Hermite, and linear interpolation polynomials to approximate the state terms, time delay terms, and static force terms in the spatial state equations, respectively;

[0094] Step 5.2: Substitute the state term, time delay term, and static force term into the spatial state equation to obtain the following discrete dynamic mapping:

[0095] (F k+1 -I)x k+1 +(T1+F k )x k +F k-1 x k-1 +F k-2 x k-2 =F k-m x k-m +F k-m+1 x k-m+1 +F k-m+2 x k-m+2 -G k f k -G k+ 1f k+1

[0096] Step 5.3: Use the extended state vector z k = [x1,x2,…,x m+1 ] T y k =[f1,f2,…,f m+1 ] T Then the discrete dynamic mapping is expressed as:

[0097] D1z k =D2z k-m +D3y k

[0098] In the formula,

[0099]

[0100]

[0101]

[0102] Step 5:4: Establish the state transfer relationship of the milling finishing system in a single time cycle;

[0103] z m =Φz0+D3y k

[0104] In the formula,

[0105] Φ=(D1) -1 D2

[0106] Based on Floquet theory, the stability of a milling finishing system is determined by the modulus of the eigenvalues ​​of the transfer matrix Φ, as shown in the following formula:

[0107]

[0108] Preferably, the specific formula for approximating the state terms in the spatial state equation using a third-order Newton interpolation polynomial in step 5.1 is as follows:

[0109] κ(δ) is determined by the state vector at time t k-2 t k-1 t k and t k+1 Response κ at time k-2 κ k-1 κ k and κ k+1 Approximately:

[0110] κ(δ)=a1κ k+1 +b1κ k +c1κ k-1 +d1κ k-2

[0111] In the formula, κ k+1 =A k+1 x k+1 ,κ k =A k x k ,κ k-1 =A k-1 x k-1 ,κ k-2 =A k-2 x k-2 In the formula t k-2 =(k-2)τ, κ((k+1)τ) is denoted as κ k+1 =A((k+1)τ)x((k+1)τ), A((k+1)τ) is recorded as Ak+1 The interpolation coefficients a1, b1, c1, and d1 are respectively:

[0112]

[0113]

[0114]

[0115]

[0116] Preferably, the specific formula for approximating the time delay term in the spatial state equation using the Hermite interpolation polynomial in step 5.1 is as follows:

[0117] α(δ-T)=a2α k-m +b2α k-m+1 +c2α k-m+2

[0118] In the formula, α 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 In the formula, α(kmT) is denoted as α k-m =B(km)x(kmT), B(kτ) is denoted as B k The interpolation coefficients a2, b2, and c2 are respectively:

[0119]

[0120]

[0121]

[0122] Preferably, the specific method for approximating the static force term in the spatial state equation using a linear interpolation polynomial in step 5.1 is as follows:

[0123]

[0124] In the formula, f(kτ) is denoted as f k .

[0125] Preferably, step 6: extracting the corresponding steady-state coefficient vector Y-direction coefficient and predicting the dynamic machining error is specifically done as follows:

[0126] Step 6.1: For the multi-cycle case, the state transfer relationship between the i-th cycle and the (i-1)-th cycle is expressed as follows:

[0127] z i =Φz i-1 +D3y k

[0128] Step 6.2: During steady-state cutting, the steady-state coefficient is:

[0129] z k =z k-m =z*

[0130] Step 6.3: Substitute the steady-state coefficients into the state transitivity expression to obtain:

[0131] z * =(D1-D2) -1 D3y k

[0132] Step 6.4: Extract the corresponding steady-state coefficient vector Y-direction coefficient from the formula in Step 6.3 to predict dynamic processing error.

[0133] On the other hand, the present invention also provides an electronic device, including a memory, a processor, and a computer program stored in the memory and executable on the processor, characterized in that the processor, when executing the program, implements a method for comprehensively predicting the stability and surface position error of milling finishing.

[0134] As can be seen from the above technical solution, compared with the prior art, the present invention discloses a method and device for comprehensive prediction of stability and surface position error in milling finishing, which is used for comprehensive prediction of cutting stability and surface position error in milling finishing. It considers the influence of modal coupling, process damping and tool runout effect in the milling process, which improves the accuracy of the analytical model. At the same time, by establishing stability lobe diagram and surface position error diagram considering modal coupling, process damping and tool runout, it is possible to conveniently and quickly predict the machining state and machining error based on the graphics, which provides a theoretical basis for optimal selection of machining parameters with constraints of no chatter and dynamic machining error less than the expected value. Attached Figure Description

[0135] To more clearly illustrate the technical solutions in the embodiments of the present invention or the prior art, the drawings used in the description of the embodiments or the prior art will be briefly introduced below. Obviously, the drawings described below are only embodiments of the present invention. For those skilled in the art, other drawings can be obtained based on the provided drawings without creative effort.

[0136] Figure 1 This is a schematic diagram of the method flow of the present invention;

[0137] Figure 2Stability lobe diagram considering modal coupling, process damping, and tool runout for this invention;

[0138] Figure 3 This is a comprehensive prediction diagram of milling stability and surface position error in this invention;

[0139] Figure 4 This is a stability diagram of the leaflets with different eccentricities according to the present invention;

[0140] Figure 5 This is a diagram showing the stability of the leaflets at different eccentric angles according to the present invention. Detailed Implementation

[0141] The technical solutions of the embodiments of the present invention will be clearly and completely described below with reference to the accompanying drawings. Obviously, the described embodiments are only some embodiments of the present invention, and not all embodiments. Based on the embodiments of the present invention, all other embodiments obtained by those skilled in the art without creative effort are within the scope of protection of the present invention.

[0142] See appendix Figure 1 As shown in the figure, this invention discloses a method for comprehensively predicting the stability and surface position error of milling finishing, based on a milling finishing system, including the following steps:

[0143] Step 1: Obtain milling cutter information, which includes: tool geometry parameters, modal parameters, cutting force coefficients, and cutting conditions;

[0144] Step 2: Calculate the cutting edge angle position, cutting radius, feed rate, and total cutting force in the X and Y directions at any given time based on the milling cutter information;

[0145] Step 3: Construct a system of differential equations for milling dynamics and transform it into a spatial state equation;

[0146] Step 4: Divide the cutting cycle of the cutting tooth into several small intervals. Perform integral calculation on the spatial state equation within the time length of any small interval to determine the spindle speed range and its corresponding discrete number, the axial depth of cut range and its corresponding discrete number, and the discrete number of a single cutting tooth cycle.

[0147] Step 5: Use third-order Newton's equation, Hermite's equation, and linear interpolation polynomial to approximate the state terms, time delay terms, and static force terms in the spatial state equation, construct the state transfer relationship of the milling and finishing system in a single time period, and determine the stability of the milling and finishing system.

[0148] Step 6: Extract the corresponding steady-state coefficient vector Y-direction coefficient to predict dynamic machining error;

[0149] Step 7: Based on Steps 5 and 6, obtain the stability lobe diagram and the surface position error prediction diagram.

[0150] In one specific embodiment, step 2, which calculates the cutting edge angle position, cutting radius, feed rate, and total cutting force in the X and Y directions at any given time based on the milling cutter information, comprises the following steps:

[0151] Step 2.1: Calculate the tool angle position of the cutting edge at any given time;

[0152] Step 2.2: Calculate the actual cutting radius and feed per tooth of the cutting edge under the condition of tool eccentricity;

[0153] Step 2.3: Calculate the cutting forces in the tangential and radial directions of the tool, and convert the cutting forces in the tangential and radial directions of the tool to the X and Y axes of the machine tool coordinate system through Cartesian coordinate transformation;

[0154] Step 2.4: Convert the cutting forces in the tangential and radial directions of the tool to the X and Y axes of the machine tool coordinate system through Cartesian coordinate system transformation. Calculate the total cutting forces acting on the milling cutter in the X and Y directions based on the cumulative cutting forces acting on the cutter teeth.

[0155] Specifically, for a milling cutter with N cutting edges, the instantaneous radial contact angle φ on the j-th tooth is... j (j=1...N) can be expressed as a functional relationship between time t and spindle speed Ω:

[0156]

[0157] Manufacturing errors in cutting tools and installation errors on the machine tool can cause the tool axis to misalign with the machine tool spindle, resulting in tool eccentricity. Multi-tooth end mills exhibit tool eccentricity during rotary machining. This eccentricity will alter the actual cutting radius and feed per tooth.

[0158] Specifically, the steps for calculating the actual cutting radius and feed per tooth of the cutting edge under eccentric tool conditions are as follows:

[0159] The actual cutting radius of each tooth is:

[0160]

[0161] The actual feed per tooth is:

[0162]

[0163] In the formula, f t a is the nominal feed per tooth. r This represents the radial cutting depth.

[0164] Based on geometric relationships, the angle between two adjacent cutting teeth can be obtained using the law of cosines:

[0165]

[0166] The entry and exit angles of the j-th tooth are expressed as:

[0167]

[0168] Specifically, the cutting forces in the tangential and radial directions of the tool are calculated, and then transformed into the X and Y axes of the machine tool coordinate system using Cartesian coordinate transformation. The milling forces acting on the cutter tooth j mainly include shearing force and process damping force, which are expressed as follows:

[0169] F t =F ts +F tp ,F r =F rs +F rp (6)

[0170] In the formula, F ts and F rs F represents the tangential and normal forces caused by the shear effect. tp and F rp These are the tangential and normal forces caused by the process damping effect. Based on the shear force model, F ts and F rs They are represented as follows:

[0171] F ts =g(φ j (t))K t a p h(t),F rs =g(φ j (t))K n a p h(t) (7)

[0172] In the formula, K t and K n These are the tangential and radial cutting force coefficients, respectively, a p For the axial cutting depth, the window function g(φ) j (t) is defined as:

[0173]

[0174] In the formula, φ st and φ ex These are the entry and exit angles of the cutting teeth, respectively. Based on regenerative chatter theory, the instantaneous cutting thickness can be expressed as...

[0175] h(t) = [ft +x(tT)-x(t)]sin(φ j (t))+[y(tT)-y(t)]cos(φ j (t)) (9)

[0176] In the formula f t The feed per tooth is T, and the time delay is equal to the cutting cycle of the cutting tooth, i.e., T = 60. / (NΩ).

[0177] Based on the process damping mechanism, the radial process damping force can be expressed as:

[0178] F rp =g(φ j (t))K sp a p S (10)

[0179] In the formula, K sp S and S represent the indentation coefficient and area, respectively. The tangential force generated by the process damping effect can be obtained from the following equation.

[0180] F tp =μF rp (11)

[0181] In the formula, μ is the coefficient of friction. The equivalent viscous damping can be expressed as:

[0182]

[0183] In the formula, v is the cutting speed, and W is the width of the tool wear band.

[0184] Specifically, by accumulating the cutting forces acting on the cutter teeth and calculating the total cutting forces acting on the milling cutter in the X and Y directions:

[0185]

[0186] Substituting equations (6), (7), (11), and (12) into equation (13), the resultant cutting force can be rewritten as:

[0187]

[0188] In the formula,

[0189]

[0190]

[0191]

[0192]

[0193] In one specific embodiment, the tool dynamics equations that comprehensively consider regeneration effects, process damping, modal coupling effects, and tool runout are expressed as the following set of time-delay differential equations:

[0194]

[0195] In the formula,

[0196]

[0197] make With the aid of Cauchy transformation, formula (19) can be transformed into a spatial state form.

[0198]

[0199] In the formula,

[0200]

[0201] In a specific embodiment, to simplify the algorithm derivation process, κ(t) = A(t)x(t) and α(tT) = B(t)x(tT) are defined. With x(t0) as the initial condition, the general solution of formula (21) is expressed as:

[0202]

[0203] Divide the cutting cycle T into m equal time intervals, then the time step τ = T / m. Formula (23) can be expressed as follows for the time interval kτ≤t≤(k+1)τ, (k=1,2,3,…m)

[0204]

[0205] Dividing the cutting cycle T of the cutting teeth into m equal time intervals, the time step τ = T / m; the above formula can be expressed as follows for the time interval kτ≤t≤(k+1)τ, (k=1,2,3,…m):

[0206]

[0207] In the formula, δ=ξ-kτ,δ∈[0,τ]; in the formula, x((k+1)τ) is denoted as x k+1 .

[0208] The maximum permissible cutting speed is determined based on several major categories of constraints, including machine tool constraints (such as maximum permissible cutting power and maximum permissible cutting force), tool constraints (such as tool material and tool life), and part constraints (such as the accuracy and material of the workpiece). Then, the minimum permissible spindle speed n that the milling machining system can withstand is determined according to the maximum permissible cutting speed. min Maximum spindle speed nmax Thus, the spindle speed range [n] can be obtained. min ,n max Discrete number k of the spindle speed range. n It is usually set to 200.

[0209] The maximum permissible axial depth of cut is determined based on machine tool rigidity, tool diameter, and machining accuracy requirements. Then, considering that a larger axial depth of cut is typically used during roughing to achieve a higher material removal rate and a smaller axial depth of cut is typically used during finishing to achieve higher machining accuracy, the acceptable minimum permissible spindle speed n is determined. min Maximum spindle speed n max Thus, the axial cutting depth range [w] can be obtained. min ,w max The discrete number k of the axial depth range. w It is usually set to 100.

[0210] The discrete number m of the cutting cycle T of the cutting tooth is usually set to 10 levels, such as 20, 30, 40, 50, 60, 70, 80, 90, 100, and 200.

[0211] In the process of solving chatter stability, the discrete number of spindle speed range, the discrete number of axial depth of cut range, and the discrete number of single tooth cutting cycle can be reduced or increased depending on the actual needs of calculation accuracy and efficiency.

[0212] In one specific embodiment, the state term and the time delay term in formula (24) are approximated by third-order Newton and Hermite interpolation, respectively. For the static force term f(δ), linear interpolation is used for approximation.

[0213] For the state term κ(δ) in formula (24), the third-order Newton interpolation polynomial is used, that is, κ(δ) is obtained by the state vector in t k-2 t k-1 t k and t k+1 Response κ at time k-2 κ k-1 κ k and κ k+1 Approximately:

[0214] κ(δ)=a1κ k+1 +b1κ k +c1κ k-1 +d1κ k-2 (25)

[0215] In the formula, the interpolation coefficients a1, b1, c1 and d1 are respectively

[0216]

[0217]

[0218]

[0219]

[0220] κ k+1 =A k+1 x k+1 ,κ k =A k x k ,κ k-1 =A k-1 x k-1 ,κ k-2 =A k-2 x k-2 (30)

[0221] In one specific embodiment, the time delay term α(δ-T) in formula (24) is approximated using a Hermitian interpolation polynomial, as follows:

[0222] α(δ-T)=a2α k-m +b2α k-m+1 +c2α k-m+2 (31)

[0223] In the formula, the interpolation coefficients a2, b2, and c2 are respectively:

[0224]

[0225]

[0226]

[0227] α 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 (35)

[0228] In a specific embodiment, the static force term f(δ) in formula (24) is linearly approximated with the corresponding boundary values, i.e.:

[0229]

[0230] Substituting equations (25)-(36) into equation (24) yields the following discrete dynamic mapping:

[0231]

[0232] In the formula,

[0233]

[0234]

[0235]

[0236]

[0237]

[0238]

[0239]

[0240]

[0241]

[0242] Equation (37) can be rearranged to obtain:

[0243]

[0244] In one specific embodiment, an extended state vector z is used. k =[x1,x2,…,x m+1 ] T y k =[f1,f2,…,f m+1 ] T Then formula (47) can be expressed as

[0245] D1z k =D2z k-m +D3y k (48)

[0246] In the formula,

[0247]

[0248]

[0249]

[0250] The state transition relationship of the system over a single time period can be represented as:

[0251] z m =Φz0+D3y k (52)

[0252] In the formula,

[0253] Φ=(D1)- 1 D2 (53)

[0254] Based on Floquet theory, the stability of a system is determined by the modulus of the eigenvalues ​​of the transfer matrix Φ, that is:

[0255]

[0256] In one specific embodiment, the method for extracting the corresponding steady-state coefficient vector Y-direction coefficient and predicting dynamic processing error is as follows:

[0257] For the multi-cycle case, the state transfer relationship between the i-th cycle and the (i-1)-th cycle can be expressed as:

[0258] z i =Φz i-1 +D3y k (55)

[0259] During steady-state cutting, the steady-state coefficient can be obtained from the fixed point of formula (52):

[0260] z k =z k-m =z * (56)

[0261] Substituting formula (56) into formula (55) yields:

[0262] z * =(D1-D2)- 1 D3y k (57)

[0263] Extract the corresponding steady-state coefficient vector Y direction coefficient from formula (57) to predict dynamic processing error.

[0264] On the other hand, Embodiment 1 of the present invention also provides an electronic device, including a memory, a processor, and a computer program stored in the memory and executable on the processor, characterized in that the processor executes the program to realize a comprehensive prediction method for the stability of milling finishing and surface position error.

[0265] In summary, in one specific embodiment, the present invention converts the cutting radius and feed rate, and then the cutting forces acting on the tool in the tangential and radial directions, onto the X and Y axes of the machine tool coordinate system. The total cutting forces acting on the milling cutter in the X and Y directions are calculated by summing these forces. Then, within the framework of structural dynamics, a set of milling dynamic time-delay differential equations considering modal coupling effects, process damping, and tool runout is constructed. These equations are integrated within any small interval, and the state terms, time-delay terms, and static terms in the integral are approximated using third-order Newtonian, Hermite, and linear interpolation polynomials, respectively. This allows for the construction of the system's state transfer relationship over a single time period. Finally, based on Floquet theory and fixed-point theory, milling stability leaflet diagrams and surface position error prediction diagrams are simultaneously output. During finishing, surface quality should be prioritized, and maximum machining efficiency should be ensured while maintaining surface quality. The method of this invention can achieve comprehensive prediction of cutting stability and surface position error in milling finishing, thereby accurately predicting the surface quality during the finishing process, thus improving machining accuracy and quality, and exhibiting excellent accuracy and timeliness.

[0266] Example 2

[0267] Given the following parameters: number of milling cutter teeth N = 3, cutter diameter D = 25.4 mm; cutter wear band width W = 0.08 mm, and indentation coefficient K... sp =1.5×10 14 N / m 3 The Coulomb force coefficient μ = 0.3; the tangential and normal force coefficients are K and K, respectively. t =7×10 8 N / m 2 and K n =0.49×10 8 N / m 2 The tool eccentricity parameters are eccentricity ρ = 7.2 μm and eccentricity angle θ0 = 65.09°, as detailed in Table 1. The spindle cycle is divided into 100 small intervals, and the plane formed by the spindle speed and radial depth of cut is divided into a 200×200 grid.

[0268] Table 1 Milling-related parameters

[0269]

[0270] The above steps and parameters are programmed using Matlab software to draw stability lobe diagrams and surface position error diagrams. The stability diagrams are used to predict the stability during the milling process and the surface position error during machining. Figure 2 This is a stability lobe diagram constructed according to a preferred embodiment of the present invention. Figure 3The stability lobe diagram and the machining surface position error diagram are constructed simultaneously according to the preferred embodiment of the present invention. The stability lobe diagrams are obtained by taking different eccentricities and eccentric angles, as shown below. Figure 4 and Figure 5 As shown.

[0271] The various embodiments in this specification are described in a progressive manner, with each embodiment focusing on its differences from other embodiments. Similar or identical parts between embodiments can be referred to interchangeably. For the apparatus disclosed in the embodiments, since they correspond to the methods disclosed in the embodiments, the description is relatively simple; relevant parts can be referred to the method section.

[0272] The above description of the disclosed embodiments enables those skilled in the art to make or use the invention. Various modifications to these embodiments will be readily apparent to those skilled in the art, and the general principles defined herein may be implemented in other embodiments without departing from the spirit or scope of the invention. Therefore, the invention is not to be limited to the embodiments shown herein, but is to be accorded the widest scope consistent with the principles and novel features disclosed herein.

Claims

1. A method for comprehensively predicting stability and surface position error in milling finishing, based on a milling finishing system, characterized in that, Includes the following steps: Step 1: Obtain milling cutter information; Step 2: Calculate the cutting edge angle position, cutting radius, feed rate, and total cutting force in the X and Y directions at any given time based on the milling cutter information; Step 3: Construct a system of differential equations for milling dynamics and transform it into a spatial state equation; Step 4: Divide the cutting cycle T of the cutting tooth into several time intervals. Perform integral calculation on the spatial state equation within the time length of any time interval to determine the spindle speed range and its corresponding discrete number, the axial depth of cut range and its corresponding discrete number, and the discrete number of a single cutting tooth cycle. Step 5: Calculate the state terms, time delay terms, and static force terms in the spatial state equation, construct the state transfer relationship of the milling and finishing system in a single time period, and determine the stability of the milling and finishing system; Step 6: Extract the corresponding steady-state coefficient vector Y-direction coefficient to predict dynamic machining error; Step 7: Based on Steps 5 and 6, obtain the stability lobe diagram and the surface position error prediction diagram.

2. The method for comprehensively predicting stability and surface position error in milling finishing according to claim 1, characterized in that, The specific steps of step 2, which calculates the cutting edge angle position, cutting radius, feed rate, and total cutting force in the X and Y directions at any given time based on the milling cutter information, are as follows: Step 2.1: Calculate the tool angle position of the cutting edge at any given time; Step 2.2: Calculate the actual cutting radius and feed per tooth of the cutting edge under the condition of tool eccentricity; Step 2.3: Calculate the cutting forces in the tangential and radial directions of the tool, and convert the cutting forces in the tangential and radial directions of the tool to the X and Y axes of the machine tool coordinate system through Cartesian coordinate transformation; Step 2.4: Convert the cutting forces in the tangential and radial directions of the tool to the X and Y axes of the machine tool coordinate system through Cartesian coordinate system transformation. Calculate the total cutting forces acting on the milling cutter in the X and Y directions based on the cumulative cutting forces acting on the cutter teeth.

3. The method for comprehensively predicting stability and surface position error in milling finishing according to claim 1, characterized in that, The specific method for constructing the hindrance differential equations of milling dynamics in step 3 and transforming them into spatial state equations is as follows: Step 3.1: Construct the set of time-delayed differential equations for milling dynamics: In the formula, M is the modal mass matrix, m xx m yy For the mass in the X and Y directions, m xy m yx The coupling mass in the X and Y directions in this coordinate system is C; C is the modal damping matrix, c xx c yy For damping in the X and Y directions, c xy c yx The X and Y axes represent the coupling damping in this coordinate system; K is the modal stiffness matrix, k xx k yy Let k be the stiffness in the X and Y directions. xy k yx C represents the coupling stiffness in the X and Y directions in this coordinate system. eq For equivalent viscous damping, C eq =(K sp a p W 2 ) / (4v), K sp Let a be the compression factor. p Where C is the axial depth of cut, W is the tool wear bandwidth, and v is the cutting speed; p H is the conversion coefficient matrix of cutting force, and we have: f0(t) is the static force vector, f t This refers to the feed per tooth. Step 3.2: Let Using Cauchy transformation, the equations are converted into spatial state form, yielding the spatial state equation: In the formula, In the formula, A0, A(t), B(t), and f(t) are obtained by combining like terms. A0 is a constant coefficient matrix derived from x(t) and does not change with time, while A(t), B(t), and f(t) all change with time.

4. The method for comprehensively predicting stability and surface position error in milling finishing according to claim 1, characterized in that, Step 4 includes: The specific method for integrating the spatial state equation within any given time interval after dividing the cutting cycle T into several equal time intervals is as follows: Define κ(t) = A(t)x(t) and α(tT) = B(t)x(tT), with x(t0) as the initial condition, the general solution of the space state equation is expressed as: Dividing the cutting cycle T of the cutting teeth into m equal time intervals, the time step τ = T / m; the above formula is expressed as follows for the time interval kτ≤t≤(k+1)τ, (k=1,2,3,…m): In the formula, δ=ξ-kτ, δ∈[0,τ], and x((k+1)τ) is denoted as x k+1 。 5. The method for comprehensively predicting stability and surface position error in milling finishing according to claim 1, characterized in that, Step 5 involves calculating the state terms, time delay terms, and static force terms in the spatial state equation to construct the state transfer relationship of the milling and finishing system over a single time period. The specific steps for determining the stability of the milling and finishing system are as follows: Step 5.1: Calculate the state terms, time delay terms, and static force terms in the spatial state equation using third-order Newtonian, Hermite, and linear interpolation polynomials, respectively; Step 5.2: Substitute the state term, time delay term, and static force term into the spatial state equation to obtain the following discrete dynamic mapping: (F k+1 -I)x k+1 +(T1+F k )x k +F k-1 x k-1 +F k-2 x k-2 =F k-m x k-m +F k-m+1 x k-m+1 +F k-m+2 x k-m+2 -G k f k -G k+1 f k+1 Step 5.3: Use the extended state vector zk = [x1, x2, ..., x m+1 ] T yk = [f1, f2, ..., f m+1 ] T Then the discrete dynamic mapping is expressed as: D1z k =D2z k-m +D3y k In the formula, Step 5.4: Establish the state transfer relationship of the milling finishing system in a single time cycle; With m =Φz0+D3y k In the formula, Φ=(D1) -1 D2 Based on Floquet theory, the stability of a milling finishing system is determined by the modulus of the eigenvalues ​​of the transfer matrix Φ, as shown in the following formula:

6. The method for comprehensively predicting stability and surface position error in milling finishing according to claim 5, characterized in that, The specific formula for calculating the state terms in the spatial state equation using the third-order Newton interpolation polynomial in step 5.1 is as follows: κ(δ) is determined by the state vector at time t k-2 t k-1 t k and t k+1 Response κ at time k-2 κ k-1 κ k and κ k+1 Approximately: κ(δ)=a1κ k+1 +b1k k +c1k k-1 +d1k k-2 In the formula, κ k+1 =A k+1 x k+1 ,κ k =A k x k ,κ k-1 =A k-1 x k-1 ,κ k-2 =A k-2 x k-2 In the formula t k-2 =(k-2)τ, κ((k+1)τ) is denoted as κ k+1 =A((k+1)τ)x((k+1)τ), A((k+1)τ) is recorded as A k+1 The interpolation coefficients a1, b1, c1, and d1 are respectively:

7. The method for comprehensively predicting stability and surface position error in milling finishing according to claim 5, characterized in that, The specific formula for calculating the time delay term in the spatial state equation using Hermitian interpolation polynomials in step 5.1 is as follows: α(δ-T)=a2α k-m +b2a k-m+1 +c2a k-m+2 In the formula, α 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 In the formula, α(kmT) is denoted as α k-m =B(km)x(kmT), B(kτ) is denoted as B k The interpolation coefficients a2, b2, and c2 are respectively:

8. The method for comprehensively predicting stability and surface position error in milling finishing according to claim 5, characterized in that, The specific method for calculating the static force term in the spatial state equation using a linear interpolation polynomial in step 5.1 is as follows: In the formula, f(kτ) is denoted as f k .

9. The method for comprehensively predicting stability and surface position error in milling finishing according to claim 1, characterized in that, Step 6: Extracting the corresponding steady-state coefficient vector Y-direction coefficient, and the specific method for predicting dynamic machining errors is as follows: Step 6.1: For the multi-cycle case, the state transfer relationship between the i-th cycle and the (i-1)-th cycle is expressed as follows: With i =Φz i-1 +D3y k Step 6.2: During steady-state cutting, the steady-state coefficient is: With k =z k-m =z″ Step 6.3: Substitute the steady-state coefficients into the state transitivity expression to obtain: z * (D1-D2) -1 D3y k ; Step 6.4: Extract the corresponding steady-state coefficient vector Y-direction coefficient from the formula in Step 6.3 to predict dynamic processing error.

10. An electronic device comprising a memory, a processor, and a computer program stored in the memory and executable on the processor, characterized in that, When the processor executes the program, it implements the comprehensive prediction method for stability and surface position error of milling finishing as described in any one of claims 1 to 9.

Images