Method for identifying radial stiffness of rotating shaft of machine tool

Abstract

The application provides a machine tool rotating shaft radial stiffness identification method, and relates to the technical field of rotating shaft system dynamics. First, the machine tool rotating shaft motion error is detected, and the Donaldson inversion method is used to separate the synchronous error in the combined measurement data; the 4-degree-of-freedom dynamics model of the rotating shaft system is established by applying the Lagrange equation; the relationship between the radial stiffness and damping of the supporting bearing at different angular velocities and the synchronous motion error of the rotating shaft is determined; and the change of the bearing radial stiffness and damping with the angular position of the rotating shaft is identified. By accurately identifying the radial stiffness of the rotating shaft, the machining process of the machine tool can be better controlled, thereby improving the surface quality of the machined parts. The identification method can provide important parameter basis for machine tool design, help to optimize the machine tool structure design, and improve the overall performance of the machine tool. By regularly detecting the stiffness change of the rotating shaft, potential faults can be predicted, maintenance or replacement of parts can be performed in advance, the waste rate can be reduced, and the production efficiency can be increased.

Description

Technical Field

[0001] This invention relates to the field of dynamics technology for rotating shaft systems, and in particular to a method for identifying the radial stiffness of a machine tool rotating shaft. Background Technology

[0002] Machine tools, as the "mother machines" of industry, play a crucial role in modern manufacturing. High stiffness and appropriate damping are essential for stable cutting. Identifying dynamic stiffness and dynamic damping helps to recognize and quantify the deformation and vibration characteristics of machine tools under dynamic loads, allowing for appropriate measures to reduce these deformations and vibrations, thereby improving machining accuracy, increasing machining efficiency, and extending tool life. Therefore, identifying the dynamic stiffness and dynamic damping of CNC lathe rotary axes is one of the key technologies for ensuring machining accuracy and improving production efficiency. Summary of the Invention

[0003] The technical problem to be solved by the present invention is to provide a method for identifying the radial stiffness of a machine tool rotating shaft, which addresses the shortcomings of the prior art.

[0004] To solve the above-mentioned technical problems, the technical solution adopted by the present invention is as follows:

[0005] A method for identifying the radial stiffness of a machine tool rotating shaft includes the following steps:

[0006] Step 1: Detection and separation of machine tool rotary axis motion error; The sensor readings are represented as a combination of the projection of the actual profile of the cross section onto the coordinate axis and the rotary axis motion error; If error sampling is performed continuously over M cycles, asynchronous errors can be eliminated when M is large enough; The workpiece profile error is periodic, and the average value of the sensor readings over M cycles is a combination of the workpiece profile error and the rotary axis motion synchronization error. The Donaldson inversion method is used to separate the synchronization error from the combined measurement data;

[0007] Step 2: Apply the Lagrange equations to establish a 4-DOF dynamic model of the rotating axis system;

[0008] Step 3: Determine the relationship between the radial stiffness and damping of the support bearing and the synchronous motion error of the rotating shaft at different angular velocities;

[0009] Step 4: Identify the changes in bearing radial stiffness and damping with the angular position of the rotating shaft.

[0010] Furthermore, in step 1, the method for detecting the motion error of the machine tool rotary axis is as follows:

[0011] In the machine tool rotary axis motion error detection system, when the rotary axis is stationary, O is the geometric center of the rotary axis, chuck, and workpiece assembly, and the z-axis is its centerline; the x-axis and y-axis lie on the workpiece cross-section. Along the x-axis and y-axis, two sets of eddy current displacement sensors S1 and S2, S3 and S4 detect the displacement signals at cross-sections a and b of the workpiece, respectively; cross-sections a and b, as well as the z-axis, intersect at points O' and O'', and the distances between O and points O' and O''' are represented by L1 and (L1+L2), respectively; the rotation angle of the rotary axis is... The actual profile of section a is The actual profile of section b is The radial motion error of the rotating axis at point O is and Swing error is and The sensor readings are expressed as a combination of the projection of the actual profile of the cross section onto the coordinate axes and the motion error of the rotation axis, as shown in Equation (1);

[0012]

[0013] in, Let be the angle between the workpiece radius at section a and the positive / negative x-axis; Let be the angle between the workpiece radius at section a and the positive / negative y-axis; Let be the angle between the workpiece radius at section b and the positive / negative x-axis; Let be the angle between the workpiece radius at section b and the positive / negative y-axis;

[0014] The sensor reading is a combination of workpiece contour error and rotation axis motion error; the workpiece's surface to be measured is uniformly divided into N parts in the circumferential direction, and the position of each sampling point in the circumferential direction is... i = 1, 2, ..., N; Since the motion error of the rotating shaft is much smaller than the workpiece radius, if the center of the detection surface is L away from the center of the rotating shaft, then the sensor reading S in equation (1) n (i) Simplified as follows:

[0015]

[0016] Where n = 1, 2, 3, 4; i = 1, 2, ..., N; α1 = α3 = 0, α2 = α4 = π / 2; P1 = P3 = T2 = T4 = x, P2 = P4 = T1 = T3 = y; H1 = H2 = L1, H3 = H4 = L1 + L2; when n = 1 or 2, k = a; when n = 3 or 4, k = b.

[0017] Furthermore, in step 1, the method for separating the motion error of the rotating axis is as follows:

[0018] The average value of the sensor readings over M cycles is a combination of the workpiece contour error and the rotation axis motion synchronization error, that is:

[0019]

[0020] Among them, S n,j (i) represents the reading of sensor n at position i during the j-th rotation cycle of the rotation axis, j = 1, 2, ..., M;

[0021] The synchronization error is separated from the measurement data of equation (4) using the Donaldson inversion method; the first measurement data of equation (4) is denoted as A, and the data measured after the workpiece is rotated by π is denoted as B, then:

[0022]

[0023] Solving equation (5), we obtain the synchronous motion error of the rotating shaft. and workpiece surface morphology as follows:

[0024]

[0025]

[0026]

[0027] Where, β n =2π-α n ;

[0028] Substituting equation (6) into equation (3), we can solve for the oscillation error of the rotating shaft. Right now:

[0029]

[0030] Where i = 1, 2, ..., N; n = 1, 2;

[0031] Substituting equation (9) into equation (3), we can solve for the radial motion error of the rotating shaft. Right now:

[0032]

[0033] Where i = 1, 2, ..., N; n = 1, 2.

[0034] Furthermore, in step 2, the establishment of the dynamic model is as follows:

[0035] In the dynamic model of the rotating axis system, the reference frame OXYZ is an inertial coordinate system, and point C is fixed on the rotating axis, coinciding with point O when the rotating axis is stationary; the motion of the rotating axis is regarded as the sum of translation in the x and y directions and rotation about point C.

[0036] During operation, assume that point C in the reference coordinate system OXYZ is r. C ={x,y,z} T Along vector r C Translating the coordinate system OXYZ yields the coordinate system Cx1y1z1; rotating the coordinate system Cx1y1z1 around the x1 and y1 axes by angles α and β yields the coordinate system Cx2y2z2; the coordinates of point P in the coordinate system Cx2y2z2 are expressed as:

[0037] r2 P ={ecosΩt esinΩt 0} T (11)

[0038] Where e is the radial distance between point P and point C; Ω is the angular velocity of the machine tool's rotating axis;

[0039] The coordinates of point P in the coordinate system OXYZ are expressed as follows:

[0040] r P =r C +T r r2 P (12)

[0041] Among them, T r Let be the rotation matrix; assuming angles α, β, and γ are small angles, the rotation matrix from Cx1y1z1 to Cx2y2z2 is:

[0042]

[0043] in,

[0044]

[0045] The kinetic energy T of the rotating shaft system is as follows:

[0046]

[0047] Where (·) represents d() / dt, M is the mass matrix of the rotation axis, J is the moment of inertia matrix of the rotation axis, and ω is the angular velocity vector of the rotation axis, as shown below:

[0048]

[0049] Among them, m and J d J pLet L be the mass, polar moment of inertia, and diametrical moment of inertia of the rotating shaft, respectively; let l2 and l1 be the distances between the front and rear support springs and point C, respectively; and let r be the connection point between the springs and the rotating shaft. bJ When the system is stationary, the coordinates of the connection point between the spring and the rotation axis in the OXYZ frame are:

[0050] r k10 ={l x1 ,0,-l1} T r k20 ={0,l y1 ,l1} T r k30 ={l x2 ,0,-l2} T r k40 ={0,l y2 ,l2} T (15)

[0051] Among them, l x1 Indicates the connection point r between the rear support spring and the rotating shaft. k1 x-coordinate; l y1 The point connection r between the rear support spring and the rotating shaft is indicated. k2 The ordinate; l x2 Indicates the connection point r between the front support spring and the rotating shaft. k3 x-coordinate; l y2 Indicates the connection point r between the front support spring and the rotating shaft. k4 The ordinate;

[0052] The inner and outer rings are connected by ball bearings. Therefore, assuming that the connection points between the spring and the rotation axis do not rotate around the velocity axis, the coordinates of these connection points in the coordinate system Cx2y2z2 are the same as those in equation (15). Therefore, the coordinates of the connection points in OXYZ are expressed as follows:

[0053] r J =r C +T r r J0 (16)

[0054] Where J = k1, k2, k3, k4;

[0055] Assume the connection point between the spring and the foundation is r. bJ Then the vector of the original spring is l. J0 =r J0 -r bJ The vector of the deformable spring is l. J =r J -r bJ The deformation of the spring is expressed as follows:

[0056] Δl J =l J -l J0 =r J -r J0 (17)

[0057] The potential energy V of the spring is as follows:

[0058]

[0059] The viscous dissipation function D of the rotating shaft system is described as follows:

[0060]

[0061] The stiffness matrix K of the spring J Represented as:

[0062]

[0063] Where J = k1, k2, k3, k4;

[0064] The damping matrix C of the spring J Represented as:

[0065]

[0066] Where J = k1, k2, k3, k4;

[0067] Choose q = {x, y, α, β} T Using generalized coordinates, the kinetic energy equation (13), potential energy equation (18), and dissipation energy equation (19) are introduced into the Lagrange equation:

[0068]

[0069] Where j represents the degree of freedom, j = 1, 2, ..., n; by neglecting higher-order terms, the differential equations of motion for the rotating axis system are established as follows:

[0070]

[0071] Where e is the eccentricity.

[0072] Furthermore, the specific method of step 3 is as follows:

[0073] The translation of the rotation axis in the x and y directions is the radial motion error, and the rotation of the rotation axis around the x and y axes is the oscillation error; the correspondence between the variables and the motion errors in the differential equation of motion (23) of the rotation axis system is as follows:

[0074]

[0075] In equation (23), the contact stiffness and damping depend on the surface profile between the rolling element and the raceway; assuming that the surface profile changes with the position of the rotating shaft system, that is, the contact stiffness and damping in equation (23) are functions of the position of the rotating shaft system, expressed as:

[0076]

[0077] Where i = 1, 2, ..., N;

[0078] The Newmark-β method is applied to the discrete five-DOF dynamic equations of a rotating axis; within the time interval [t, t+Δt], the following assumptions are made:

[0079]

[0080]

[0081] Where Δt is the time step; θ and These are parameters that can be adjusted according to the accuracy and stability requirements of the integration process. q t+1 These represent the acceleration, velocity, and displacement at time t+1, respectively. q t Let represent the acceleration, velocity, and displacement at time t, respectively.

[0082] Equation (26) can be further written as:

[0083]

[0084] Substituting equation (27) into equation (25), we get:

[0085]

[0086] Substituting equations (27) and (28) into the five-degree-of-freedom dynamic equations of the rotating axis and rearranging, we get:

[0087]

[0088] The expressions for each coefficient in equation (29) are shown below:

[0089]

[0090]

[0091]

[0092]

[0093]

[0094]

[0095]

[0096]

[0097]

[0098]

[0099]

[0100]

[0101]

[0102]

[0103]

[0104]

[0105]

[0106]

[0107]

[0108]

[0109]

[0110]

[0111]

[0112] in, This represents the displacement of the rotation axis at position i and time t+1; (·) represents the displacement of the rotation axis at position i and time t; (·) represents d(i) / dt; (··) represents d 2 () / dt 2 γ is the phase angle;

[0113] Equation (29) can be expressed in matrix form as follows:

[0114]

[0115] Where i = 1, 2, ..., N; l = 1, 2, ..., J l ;

[0116]

[0117]

[0118] Furthermore, the specific method of step 4 is as follows:

[0119] To determine the stiffness and damping parameters of the rotating shaft, it is necessary to measure the rotating shaft at different rotational speeds. Synchronization error under J l >2; In equations (29) and (30), the exponent l is the rotational speed Ω of the rotating shaft. l The experimental results, Ω l ∈[50,3000];

[0120] Assumption:

[0121] χ={η 1 η 2 …η N e} T ;

[0122]

[0123] All experimental results of equation (29) can be expressed in matrix form, i.e.:

[0124] Cχ=h (31)

[0125] in,

[0126] Solving for vector χ using the least squares method, i.e.

[0127] χ=(C T C) -1 C T h (32)

[0128] To determine the value of the phase angle, the Monte Carlo optimization method was used to obtain the most suitable phase γ. l l = 1, 2, ..., J l Assume the residual vector of equation (31) is:

[0129] v=Cχ-h (33)

[0130] The optimization process is to find the parameter γ. l Minimize the value of the following function expression:

[0131]

[0132] Furthermore, in step 4, the parameter γ that minimizes the value of the function expression is found. l The optimization process is as follows:

[0133] Step 4.1: Use the random number function rand in the interval The internal generation produces p phase angles, that is:

[0134]

[0135] Where l = 1, 2, ..., J l If i = 1, 2, ..., p, then we get p. J Combination of group parameters

[0136] Step 4.2: Substitute these combinations into the expressions for the coefficients in equation (29) to obtain the results in formula (32). A vector h; according to equations (32)-(34), the function expression is obtained. The values ​​are as follows:

[0137]

[0138] Response parameter f * Use respectively express;

[0139] Step 4.3: Update the upper and lower limits of the parameter using the following formula:

[0140]

[0141] The first formula represents the range δ of the phase angle γ. l Reduced to half its original size;

[0142] Step 4.4: Determine whether the optimal program has ended based on the stopping condition (38);

[0143]

[0144] In the formula, δ0 is a local minimum value;

[0145] If equation (38) is satisfied, the optimal process ends. At this point, the optimal phase angle between the excitation force and the rotation shaft response is calculated as follows:

[0146]

[0147] Use parameters Calculate the stiffness, damping, and eccentricity based on the expressions for each coefficient in equation (29) and equation (32); if equation (38) is not satisfied, proceed to step 4.1.

[0148] The beneficial effects of adopting the above technical solution are as follows: The radial stiffness identification method for machine tool rotary shafts provided by this invention, by accurately identifying the radial stiffness of the rotary shaft, can better control the machining process of the machine tool, thereby improving the surface quality of the machined parts. The identification method can provide important parameter basis for machine tool design, helping to optimize the machine tool structure design and improve the overall performance of the machine tool. By regularly detecting changes in the stiffness of the rotary shaft, potential faults can be predicted, allowing for early maintenance or replacement of parts, reducing scrap rates and increasing production efficiency. Attached Figure Description

[0149] Figure 1 Figure 1a is a schematic diagram of the rotary shaft motion error detection system provided in the embodiment of the present invention; Figure 1b is a schematic diagram of the geometric relationship between the displacement signal and the motion error at section a of the workpiece; Figure 1c is a schematic diagram of the geometric relationship between the displacement signal and the motion error at section b of the workpiece.

[0150] Figure 2 Figure 2a is a dynamic model diagram of the rotating shaft system provided in the embodiment of the present invention, and Figure 2b is a schematic diagram of the coordinate system of the rotating shaft system.

[0151] Figure 3 A flowchart of a machine tool rotary shaft radial stiffness identification method provided in an embodiment of the present invention;

[0152] Figure 4 A test system diagram provided for an embodiment of the present invention;

[0153] Figure 5 The diagram shows the identification results of dynamic stiffness provided in the embodiments of the present invention.

[0154] Figure 6 The diagram shows the identification results of dynamic damping provided in an embodiment of the present invention.

[0155] In the diagram: 1. Rotary shaft; 2. Chuck; 3. Workpiece; 4. Eddy current displacement sensor; 5. Bearing 1; 6. Bearing 2; 7. Data acquisition card; 8. Computer. Detailed Implementation

[0156] The specific embodiments of the present invention will be described in further detail below with reference to the accompanying drawings and examples. The following examples are for illustrative purposes only and are not intended to limit the scope of the invention.

[0157] The method for identifying the radial stiffness of the machine tool rotating shaft in this embodiment is described in detail below.

[0158] Step 1: Detection and Separation of Machine Tool Rotary Axis Motion Error. Sensor readings are represented as a combination of the projection of the actual profile of the cross-section onto the coordinate axes and the rotary axis motion error. If error sampling is performed continuously over M cycles, asynchronous errors can be eliminated when M is sufficiently large. The workpiece profile error is periodic; the average value of sensor readings over M cycles is a combination of the workpiece profile error and the rotary axis motion synchronization error. The synchronization error is separated from the measurement data of this combination using the Donaldson inversion method.

[0159] Step 1-1: Method for detecting spindle motion error.

[0160] The principle of the machine tool rotary axis motion error detection system is as follows: Figure 1 As shown in Figure (1a), when the rotating shaft 1 is stationary, O is the geometric center of the combination of rotating shaft 1, chuck 2, and workpiece 3, and the z-axis is its centerline. The x-axis and y-axis are located on the cross section of workpiece 3. Along the x-axis and y-axis, the eddy current displacement sensor 4 includes two sets (S1 and S2, S3 and S4) to detect the displacement signals at cross section a and cross section b of the workpiece, respectively. As shown in Figures (1b) and (1c), cross sections a and b and the z-axis intersect at points O' and O”, and the distances between O and points O' and O” are represented by L1 and (L1+L2), respectively. The rotation angle of the rotating shaft is... The actual profile of section a is The actual profile of section b is The radial motion error of the rotating axis at point O is and Swing error is and The reading of sensor 4 is expressed as a combination of the projection of the actual profile of the cross section onto the coordinate axis and the motion error of the rotation axis, as shown in equation (1);

[0161]

[0162] in, Let be the angle between the workpiece radius at section a and the positive / negative x-axis; Let be the angle between the workpiece radius at section a and the positive / negative y-axis; Let be the angle between the workpiece radius at section b and the positive / negative x-axis; Let be the angle between the workpiece radius at section b and the positive / negative half-axis of y.

[0163] The sensor reading is a combination of workpiece contour error and rotation axis motion error. The workpiece's surface to be measured is uniformly divided into N parts along the circumference, and the position of each sampling point along the circumference is... i = 1, 2, ..., N. Since the motion error of the rotating shaft is much smaller than the workpiece radius, if the center of the detection surface is L away from the center of the rotating shaft, then the sensor reading S in equation (1) is... n (i) Simplified as follows:

[0164]

[0165] Where n = 1, 2, 3, 4; i = 1, 2, ..., N; α1 = α3 = 0, α2 = α4 = π / 2; P1 = P3 = T2 = T4 = x, P2 = P4 = T1 = T3 = y; H1 = H2 = L1, H3 = H4 = L1 + L2; when n = 1 or 2, k = a; when n = 3 or 4, k = b.

[0166] Step 1-2: Method for separating spindle motion errors.

[0167] If the error sampling in equation (2) is performed continuously over M cycles, then when M is large enough, asynchronous errors can be eliminated because they are random. On the other hand, the workpiece contour error is periodic. Therefore, the average value of the sensor readings over M cycles is a combination of the workpiece contour error and the rotation axis motion synchronization error, i.e.:

[0168]

[0169] Among them, S n,j (i) represents the reading of sensor n at position i during the j-th rotation cycle of the rotating axis, where j = 1, 2, ..., M.

[0170] The synchronization error is separated from the measurement data in equation (4) using the Donaldson inversion method. Let A be the first measurement data in equation (4), and B be the data measured after the workpiece is rotated by π. Then:

[0171]

[0172] Solving equation (5), we obtain the synchronous motion error of the rotating shaft. and workpiece surface morphology as follows:

[0173]

[0174]

[0175]

[0176] Where, β n =2π-α n .

[0177] Substituting equation (6) into equation (3), we can solve for the oscillation error of the rotating shaft. Right now:

[0178]

[0179] Where i = 1, 2, ..., N; n = 1, 2.

[0180] Substituting equation (9) into equation (3), we can solve for the radial motion error of the rotating shaft. Right now:

[0181]

[0182] Where i = 1, 2, ..., N; n = 1, 2.

[0183] Step 2: Apply the Lagrange equations to establish a 4-DOF dynamic model of the rotating axis system;

[0184] The dynamic model of the rotating shaft system is as follows Figure 2 As shown in Figure (2a), the reference frame OXYZ is an inertial coordinate system. Point C is fixed on the rotation axis 1, and coincides with point O when the rotation axis 1 is stationary. The motion of the rotation axis 1 can be regarded as the sum of translation in the x and y directions and rotation about point C.

[0185] During operation, assume that point C in the reference coordinate system OXYZ is r. C ={x,y,z} T Along vector r C Translating the coordinate system OXYZ yields the coordinate system Cx1y1z1. Rotating Cx1y1z1 around the x1 and y1 axes by angles α and β yields the coordinate system Cx2y2z2. The coordinates of point P in coordinate system Cx2y2z2 are expressed as follows:

[0186] r2 P ={ecosΩt esinΩt 0} T (11)

[0187] Where e is the radial distance between point P and point C; Ω is the angular velocity of the machine tool's rotating axis.

[0188] The coordinates of point P in the coordinate system OXYZ are expressed as follows:

[0189] r P =r C +T r r2 P (12)

[0190] Among them, T rLet be the rotation matrix; assuming angles α, β, and γ are small angles, the rotation matrix from Cx1y1z1 to Cx2y2z2 is:

[0191]

[0192] in,

[0193]

[0194] The kinetic energy T of the rotating shaft system is as follows:

[0195]

[0196] Where (·) represents d() / dt, M is the mass matrix of the rotation axis, J is the moment of inertia matrix of the rotation axis, and ω is the angular velocity vector of the rotation axis, as shown below:

[0197]

[0198] Among them, m and J d J p Let L be the mass, polar moment of inertia, and radial moment of inertia of the rotating shaft, respectively. The distances between the front and rear support springs and point C are l2 and l1, respectively. Assume the connection point between the springs and the rotating shaft is r. bJ When the system is stationary, the coordinates of the connection point between the spring and the rotation axis in the OXYZ frame are:

[0199] r k10 ={l x1 ,0,-l1} T r k20 ={0,l y1 ,l1} T r k30 ={l x2 ,0,-l2} T r k40 ={0,l y2 ,l2} T (15)

[0200] Among them, l x1 Indicates the connection point r between the rear support spring and the rotating shaft. k1 x-coordinate; l y1 The point connection r between the rear support spring and the rotating shaft is indicated. k2 The ordinate; l x2 Indicates the connection point r between the front support spring and the rotating shaft. k3 x-coordinate; l y2 Indicates the connection point r between the front support spring and the rotating shaft. k4 The ordinate.

[0201] The inner and outer rings are connected by ball bearings. Therefore, assuming that the connection points between the spring and the rotation axis do not rotate around the velocity axis, the coordinates of these connection points in the coordinate system Cx2y2z2 are the same as those in equation (15). Therefore, the coordinates of the connection points in OXYZ are expressed as follows:

[0202] r J =r C +T r r J0 (16)

[0203] Where J = k1, k2, k3, k4.

[0204] Assume the connection point between the spring and the foundation is r. bJ Then the vector of the original spring is l. J0 =r J0 -r bJ The vector of the deformable spring is l. J =r J -r bJ The deformation of the spring is expressed as follows:

[0205] Δl J =l J -l J0 =r J -r J0 (17)

[0206] The potential energy V of the spring is as follows:

[0207]

[0208] The viscous dissipation function D of the rotating shaft system is described as follows:

[0209]

[0210] The stiffness matrix K of the spring J Represented as:

[0211]

[0212] Where J = k1, k2, k3, k4;

[0213] The damping matrix C of the spring J Represented as:

[0214]

[0215] Where J = k1, k2, k3, k4;

[0216] Choose q = {x, y, α, β} T Using generalized coordinates, the kinetic energy equation (13), potential energy equation (18), and dissipation energy equation (19) are introduced into the Lagrange equation:

[0217]

[0218] Where j represents the degree of freedom, j = 1, 2, ..., n. By neglecting higher-order terms, the differential equations of motion for the rotating axis system are established as follows:

[0219]

[0220] Where e is the eccentricity.

[0221] Step 3: Determine the relationship between the radial stiffness and damping of the support bearing and the synchronous motion error of the rotating shaft at different angular velocities.

[0222] The translation of the rotation axis in the x and y directions is the radial motion error, and the rotation of the rotation axis around the x and y axes is the oscillation error. The correspondence between the variables and the motion errors in the differential equation of motion (23) of the rotation axis system is as follows:

[0223]

[0224] In equation (23), the contact stiffness and damping depend on the surface profile between the rolling element and the raceway; assuming that the surface profile changes with the position of the rotating shaft system, that is, the contact stiffness and damping in equation (23) are functions of the position of the rotating shaft system, expressed as:

[0225]

[0226] Where i = 1, 2, ..., N. At this point, 4N stiffness values ​​and 4N damping values ​​need to be determined; in addition to the stiffness and damping parameters, the eccentricity e is also unknown. Therefore, 8N+1 parameters need to be identified.

[0227] The Newmark-β method is applied to the discrete five-DOF dynamic equations of a rotating axis. Within the time interval [t, t+Δt], the following assumptions are made:

[0228]

[0229]

[0230] Where Δt is the time step; θ and These are parameters that can be adjusted according to the accuracy and stability requirements of the integration process. q t+1 These represent the acceleration, velocity, and displacement at time t+1, respectively. q t Let represent the acceleration, velocity, and displacement at time t, respectively.

[0231] Equation (26) can be further written as:

[0232]

[0233] Substituting equation (27) into equation (25), we get:

[0234]

[0235] Substituting equations (27) and (28) into the five-degree-of-freedom dynamic equations of the rotating axis and rearranging, we get:

[0236]

[0237] The expressions for each coefficient in equation (29) are shown below:

[0238]

[0239]

[0240]

[0241]

[0242]

[0243]

[0244]

[0245]

[0246]

[0247]

[0248]

[0249]

[0250]

[0251]

[0252]

[0253]

[0254]

[0255]

[0256]

[0257]

[0258]

[0259]

[0260]

[0261] in, This represents the displacement of the rotation axis at position i and time t+1; (·) represents the displacement of the rotation axis at position i and time t; (·) represents d(i) / dt; (··) represents d 2 () / dt 2 γ is the phase angle.

[0262] Equation (29) can be expressed in matrix form as follows:

[0263]

[0264] Where i = 1, 2, ..., N; l = 1, 2, ..., J l .

[0265]

[0266]

[0267] Step 4: Identify the changes in bearing radial stiffness and damping with the angular position of the rotating shaft.

[0268] As shown in equation (30), 4N+1 equations can be obtained. Furthermore, the phase between the excitation force and the main shaft response is also unknown, therefore the number of identified parameters is greater than the number of equations. To determine the stiffness and damping parameters of the rotating shaft, it is necessary to measure the rotating shaft at different rotational speeds. Synchronization error under J l >2; In equations (29) and (30), the exponent l is the rotational speed Ω of the rotating shaft. l The experimental results, Ω l ∈[50,3000].

[0269] Assumption:

[0270] χ={η 1 η 2 …η N e} T ;

[0271]

[0272] All experimental results of equation (29) can be expressed in matrix form, i.e.:

[0273] Cχ=h (31)

[0274] in,

[0275] Solving for vector χ using the least squares method, i.e.

[0276] χ=(C T C) -1 C T h (32)

[0277] Because the phase angle γ between the excitation force and the spindle response is a nonlinear parameter, the parameter in (32) cannot be directly determined. To determine the value of the phase angle, the Monte Carlo optimization method is used to obtain the most suitable phase γ. l l = 1, 2, ..., J l Assume the residual vector of equation (31) is:

[0278] v=Cχ-h (33)

[0279] The optimization process is to find the parameter γ. l Minimize the value of the following function expression:

[0280]

[0281] like Figure 3 As shown, the optimization process is as follows:

[0282] Step 4.1: Use the random number function rand in the interval The internal generation produces p phase angles, that is:

[0283]

[0284] Where l = 1, 2, ..., J l Let i = 1, 2, ..., p. Then we get p. J Combination of group parameters

[0285] Step 4.2: Substitute these combinations into the expressions for the coefficients in equation (29) to obtain the results in formula (32). There are vectors h. Based on equations (32)-(34), the function expression is obtained. The values ​​are as follows:

[0286]

[0287] Response parameter f * Use respectively express.

[0288] Step 4.3: Update the upper and lower limits of the parameter using the following formula:

[0289]

[0290] The first formula represents the range δ of the phase angle γ. l Reduced to half its original size.

[0291] Step 4.4: Determine whether the optimal program has ended based on the stopping condition (38).

[0292]

[0293] In the formula, δ0 is a local minimum.

[0294] If equation (38) is satisfied, the optimal process ends. At this point, the optimal phase angle between the excitation force and the rotating shaft response is calculated as follows:

[0295]

[0296] Use parameters The stiffness, damping, and eccentricity are calculated based on the expressions for the coefficients in equation (29) and equation (32).

[0297] If equation (38) is not satisfied, proceed to step 4.1.

[0298] like Figure 4 As shown, the experimental machine tool is a commercial CNC 6130 lathe. Four eddy current sensors 4, model ML33-01-00-03, are installed on the CNC machine tool to detect the dynamic response of the rotating shaft 1 on the workpiece 3. Sensors S1 and S2 detect the displacement signal at section a of the workpiece 3, while S3 and S4 detect the displacement signal at section b of the workpiece 3. S1 and S3 detect the displacement signal in the y-direction, while S2 and S4 detect the displacement signal in the x-direction. The distance between the end face of the chuck 2 and sensor S1 is 75 mm, and the distance between sensors S1 and S3 is 50 mm. After the rotating shaft 1 is started, the data acquisition card 7 (NI-DAQ) inputs the sensor signals into the computer 8 at a sampling frequency of 25.6 kHz.

[0299] The dynamic stiffness and damping of a rotating shaft system, such as Figure 5 and Figure 6 As shown in the figure, the identified value of the eccentricity e of the rotating shaft system is 0.032 mm. The optimized values ​​of the phase difference obtained by the Monte Carlo method are shown in Table 1.

[0300] Table 1. Optimized values ​​for phase difference (rad)

[0301] rotational speed Workpiece 1 Workpiece 2 Workpiece 3 Workpiece Four 580r / min 0.008018 0.008114 0.007976 0.007937 600r / min 0.008323 0.008439 0.008173 0.008122 620r / min 0.008696 0.008762 0.008405 0.008465

[0302] The stiffness and damping of the bearing change periodically with the rotation angle, meaning the contact state between the rotating shaft and the bearing housing is constantly changing. The bearing provides strong nonlinear equivalent stiffness and damping to the rotating shaft. The constant variation in stiffness and damping in both directions, with different amplitudes, leads to different radial motion errors in the two directions. The horizontal direction is a machining-sensitive direction; drastic changes in stiffness and damping are detrimental to improving the quality of the machined surface. The masses of workpieces one and two are greater than those of workpieces three and four. When workpieces three and four are mounted on the chuck, the stiffness and damping of the bearings are consistent. When replaced with workpiece one, the stiffness and damping of the front support bearing increases by 3%, while the stiffness and damping of the rear support bearing decrease by 2.25%. When replaced with workpiece two, the stiffness and damping of the front support bearing increases by 6.98%, while the stiffness and damping of the rear support bearing decrease by 5.39%. The distance between the front and rear supports and the center of mass changes with the load, causing variations in the stiffness and damping of the bearings. Therefore, different machining parameters should be selected for workpieces of different masses.

[0303] Finally, it should be noted that the above embodiments are only used to illustrate the technical solutions of the present invention, and not to limit them; although the present invention has been described in detail with reference to the foregoing embodiments, those skilled in the art should understand that modifications can still be made to the technical solutions described in the foregoing embodiments, or equivalent substitutions can be made to some or all of the technical features therein; and these modifications or substitutions do not cause the essence of the corresponding technical solutions to deviate from the scope defined by the claims of the present invention.

Claims

1. A method for identifying the radial stiffness of a machine tool rotating shaft, characterized in that: Includes the following steps: Step 1: Detection and separation of machine tool rotary axis motion error; The sensor readings are represented as a combination of the projection of the actual profile of the cross section onto the coordinate axis and the rotary axis motion error; If error sampling is performed continuously over M cycles, asynchronous errors can be eliminated when M is large enough; The workpiece profile error is periodic, and the average value of the sensor readings over M cycles is a combination of the workpiece profile error and the rotary axis motion synchronization error. The Donaldson inversion method is used to separate the synchronization error from the combined measurement data; Step 2: Apply the Lagrange equations to establish a 4-DOF dynamic model of the rotating axis system; Step 3: Determine the relationship between the radial stiffness and damping of the support bearing and the synchronous motion error of the rotating shaft at different angular velocities; Step 4: Identify the changes in bearing radial stiffness and damping with the angular position of the rotating shaft.

2. The method for identifying the radial stiffness of a machine tool rotating shaft according to claim 1, characterized in that: In step 1, the method for detecting the motion error of the machine tool rotary axis is as follows: In the machine tool rotary axis motion error detection system, when the rotary axis is stationary, O is the geometric center of the rotary axis, chuck and workpiece assembly, and the z-axis is its centerline; the x-axis and y-axis are located on the workpiece cross-section, and along the x-axis and y-axis, two sets of eddy current displacement sensors S1 and S2, S3 and S4 respectively detect the displacement signals at cross-sections a and b of the workpiece; Cross sections a and b, and the z-axis intersect at points O' and O'', respectively. The distances between O and points O' and O'' are denoted by L1 and (L1+L2), respectively. The rotation angle of the rotation axis is φ, and the actual profile of cross section a is... The actual contour of section b is The radial motion error of the rotating axis at point O is: and The swing error is and The sensor readings are expressed as a combination of the projection of the actual profile of the cross section onto the coordinate axes and the motion error of the rotation axis, as shown in equation (1). (1)； Wherein, ϑ1(φ) is the angle between the workpiece radius at section a and the positive / negative x-axis; ϑ2(φ) is the angle between the workpiece radius at section a and the positive / negative y-axis; ϑ3(φ) is the angle between the workpiece radius at section b and the positive / negative x-axis; and ϑ4(φ) is the angle between the workpiece radius at section b and the positive / negative y-axis. The sensor reading is a combination of workpiece contour error and rotation axis motion error; the workpiece's surface to be measured is uniformly divided into N parts in the circumferential direction, and the position of each sampling point in the circumferential direction is... , Since the motion error of the rotating shaft is much smaller than the workpiece radius, if the center of the detection surface is L away from the center of the rotating shaft, then the sensor reading in equation (1) Simplified as follows: (2)； in, ; ; α1=α3=0, α2=α4=π / 2; P1=P3=T2=T4=x, P2=P4=T1=T3=y; H1=H2=L1, H3=H4=L1+L2; When n=1 or 2, k=a; When n=3 or 4, k=b.

3. The method for identifying the radial stiffness of a machine tool rotating shaft according to claim 2, characterized in that: In step 1, the method for separating the motion error of the rotating axis is as follows: The average value of the sensor readings over M cycles is a combination of the workpiece contour error and the rotation axis motion synchronization error, that is: (4)； in, This represents the reading of sensor n at position i during the j-th rotation cycle of the rotating axis. ; The synchronization error is separated from the measurement data of equation (4) using the Donaldson inversion method; the first measurement data of equation (4) is denoted as A, and the data measured after the workpiece is rotated by π is denoted as B, then: (5)； Solving equation (5), we obtain the synchronous motion error of the rotating shaft. and workpiece surface morphology as follows: (6)； (7)； (8)； in, ; Substituting equation (6) into equation (3), we can solve for the oscillation error of the rotating shaft. ,Right now: (9)； in, ; ; Substituting equation (9) into equation (3), we can solve for the radial motion error of the rotating shaft. ,Right now: (10)； in, ; .

4. The method for identifying the radial stiffness of a machine tool rotating shaft according to claim 3, characterized in that: In step 2, the establishment of the dynamic model is as follows: In the dynamic model of the rotating axis system, the reference frame OXYZ is an inertial coordinate system, and point C is fixed on the rotating axis, coinciding with point O when the rotating axis is stationary; the motion of the rotating axis is regarded as the sum of translation in the x and y directions and rotation about point C; During operation, assume that point C in the reference coordinate system OXYZ is... Along vector r C Translating the coordinate system OXYZ yields the coordinate system Cx1y1z1; rotating the coordinate system Cx1y1z1 around the x1 and y1 axes by angles α and β yields the coordinate system Cx2y2z2; the coordinates of point P in the coordinate system Cx2y2z2 are expressed as: (11)； Where e is the radial distance between point P and point C; Ω is the angular velocity of the machine tool's rotating axis; The coordinates of point P in the coordinate system OXYZ are expressed as follows: (12)； Among them, T r Let be the rotation matrix; assuming angles α, β, and γ are small angles, the rotation matrix from Cx1y1z1 to Cx2y2z2 is: ； in, ； The kinetic energy T of the rotating shaft system is as follows: (13)； in, express M is the mass matrix of the rotation axis, and J is the inertia matrix of the rotation axis. These are the angular velocity vectors of the rotation axis, represented as follows: ， ， (14)； Among them, m and J d J p Let L1 be the mass, polar moment of inertia, and diametrical moment of inertia of the rotating shaft; l2 and l1 be the distances between the front and rear support springs and point C, respectively; and assume that the connection point between the spring and the rotating shaft is L2. When the system is stationary, the coordinates of the connection point between the spring and the rotation axis in the OXYZ frame are: (15)； Among them, l x1 Indicates the connection point r between the rear support spring and the rotating shaft. k1 x-coordinate; l y1 The point connection r between the rear support spring and the rotating shaft is indicated. k2 The ordinate; l x2 Indicates the connection point r between the front support spring and the rotating shaft. k3 x-coordinate; l y2 Indicates the connection point r between the front support spring and the rotating shaft. k4 The ordinate; The inner and outer rings are connected by ball bearings. Therefore, assuming that the connection points between the spring and the rotation axis do not rotate around the velocity axis, the coordinates of these connection points in the coordinate system Cx2y2z2 are the same as the coordinates in equation (15). Therefore, the coordinates of the connection points in OXYZ are expressed as follows: (16)； in, ; Assume the connection point between the spring and the foundation is Then the vector of the original spring is The vector of the deformable spring is The deformation of the spring is expressed as follows: (17)； The potential energy V of the spring is as follows: (18)； The viscous dissipation function D of the rotating shaft system is described as follows: (19)； Spring stiffness matrix Represented as: (20)； in, ; , , , , ; Damping matrix of a spring Represented as: (21)； in, ; , , , , ; choose For generalized coordinates, the kinetic energy equation (13), potential energy equation (18), and dissipation energy equation (19) are introduced into the Lagrange equation: (22)； Where j represents the degree of freedom, By neglecting higher-order terms, the differential equations of motion for the rotating axis system are established as follows: (23)； Where e is the eccentricity.

5. The method for identifying the radial stiffness of a machine tool rotating shaft according to claim 4, characterized in that: The specific method for step 3 is as follows: The translation of the rotation axis in the x and y directions is the radial motion error, and the rotation of the rotation axis around the x and y axes is the oscillation error; the correspondence between the variables and the motion errors in the differential equation of motion (23) of the rotation axis system is as follows: ， ， ， ； In equation (23), the contact stiffness and damping depend on the surface profile between the rolling element and the raceway; assuming that the surface profile changes with the position of the rotating shaft system, that is, the contact stiffness and damping in equation (23) are functions of the position of the rotating shaft system, expressed as: (24)； in, ; The Newmark-β method is applied to the discrete five-DOF dynamic equations of a rotating axis; within a time interval Inside, assuming: (25)； (26)； Where Δt is the time step; θ and ς are parameters that can be adjusted according to the accuracy and stability requirements of the integration. , , These represent the acceleration, velocity, and displacement at time t+1, respectively. , , Let represent the acceleration, velocity, and displacement at time t, respectively. Equation (26) can be further written as: (27)； Substituting equation (27) into equation (25), we get: (28)； Substituting equations (27) and (28) into the five-degree-of-freedom dynamic equations of the rotating axis and rearranging, we get: (29)； The expressions for each coefficient in equation (29) are shown below: (B1)； (B2)； (B3)； (B4)； (B5)； (B6)； (B7)； (B8)； (B9)； (B10)； (B11)； (B12)； (B13)； (B14)； (B15)； (B16)； (B17)； (B18)； (B19)； (B20)； (B21)； (B22)； (B23)； Where, x t+1 ( ), y t+1 ( ), α t+1 ( ), β t+1 ( ) represents the displacement of the rotation axis at position i and time t+1; x t ( ), y t ( ), α t ( ), β t ( () represents the displacement of the rotation axis at position i and time t; express ; express ; It is the phase angle; Equation (29) can be expressed in matrix form as follows: (30)； in, ; ; ； ； ； 。 6. The method for identifying the radial stiffness of a machine tool rotating shaft according to claim 5, characterized in that: The specific method for step 4 is as follows: To determine the stiffness and damping parameters of the rotating shaft, it is necessary to measure the rotating shaft at different rotational speeds. Synchronization error under, In equations (29) and (30), the exponent l is the rotational speed of the rotating shaft. The experimental results ; Assumption: ； ； All experimental results of equation (29) can be expressed in matrix form, i.e.: (31)； in, ; Solving vectors using the least squares method ,Right now (32)； To determine the value of the phase angle, the Monte Carlo optimization method was used to obtain the most suitable phase. , Assume the residual vector of equation (31) is: (33)； The optimization process is to find the parameters. Minimize the value of the following function expression: (34)。 7. The method for identifying the radial stiffness of a machine tool rotating shaft according to claim 6, characterized in that: In step 4, the parameter that minimizes the value of the function expression is found. The optimization process is as follows: Step 4.1: Use the random number function rand in the interval The internal generation produces p phase angles, that is: (35)； in, ; Then we get Combination of group parameters ; Step 4.2: Substitute these combinations into the expressions for the coefficients in equation (29) to obtain equation (32). A vector h; according to equations (32)-(34), the function expression is obtained. The values ​​are as follows: (36)； Response parameters Use respectively express; Step 4.3: Update the upper and lower limits of the parameter using the following formula: (37)； The first formula represents the range δ of the phase angle γ. l Reduced to half its original size; Step 4.4: Determine whether the optimal program has ended based on the stopping condition (38); (38)； In the formula, δ0 is a local minimum value; If the stopping condition (38) is satisfied, the optimal process ends. At this time, the optimal phase angle between the excitation force and the response of the rotating shaft is calculated as follows: (39)； Use parameters Calculate the stiffness, damping and eccentricity according to the expressions of the coefficients in equation (29) and equation (32); if the stopping condition (38) is not satisfied, proceed to step 4.1.

Images