A method for obtaining the stability lobe map of a milling robot with variable posture

Abstract

The present invention provides a method for obtaining a stability lobe diagram of a milling robot with variable postures, including: gridding the robot's workspace, performing modal tests in the divided grids, identifying the joint stiffness and damping of the robot under different postures, establishing a dynamic model of the milling robot with variable postures, obtaining the robot's inherent vibration properties and the tool tip frequency response function, and drawing a stability lobe diagram of the robot's variable posture milling in a large range of workspaces based on a zero-order frequency domain stability algorithm. The method proposed in the present invention fully considers the impact of the robot's changing postures on milling stability during large-scale operations, combines the multi-body dynamics method with the zero-order frequency domain method to predict the milling stability domain, and can accurately and quickly select the process parameters (rotation speed, cutting depth) during stable milling, thereby avoiding vibration problems during processing.

Description

Technical Field

[0001] The invention belongs to the field of industrial robot milling stability prediction, and in particular relates to a method for acquiring a milling robot's variable posture stability lobe diagram. Background Art

[0002] With the booming development of the aerospace industry, industrial robots are required for in-situ processing in many work scenarios. As a milling processing carrier, robots are more flexible and efficient than traditional machine tools, and can perform milling operations on workpieces in complex processing sites. Therefore, they have broad development prospects in the fields of intelligent and flexible processing.

[0003] However, the rigidity of the serial structure of the robot body is poor. During milling, the end tool and the workpiece are prone to strong interaction, causing chatter in the robot-tool system. The violent vibration of the end tool produces machining noise and even damages the tool and workpiece, seriously affecting the machining accuracy and quality. Stability lobe diagram is an important means to prevent chatter in robot milling. It can avoid the occurrence of chatter by selecting stable milling parameters. Traditional stability lobe diagrams are mainly used for machine tool chatter, while research on robot chatter models is still not mature enough.

[0004] Said Mousavi et al. established a multi-body dynamic model of a serial robot based on beam elements in an article published in the academic journal International Journal of Advanced Manufacturing Technology, 2016, 88(9-12):3053-3065. They experimentally identified the geometric, elastic and damping parameters of the beam elements, used the regenerative chatter in the milling process as a function of the motion redundancy variable, and realized the drawing of stability lobe diagrams.

[0005] In an article published in the academic journal Robotics and Computer Integrated Manufacturing, 2019, 55: 11-18, Marcel Cordes et al. established a structural dynamics model of an articulated manipulator with a spindle and a tool, and analyzed the stability of the system using semi-discrete time domain and frequency domain methods. The accuracy of the stability lobe diagram was verified in milling experiments on aluminum alloy and titanium alloy parts.

[0006] The above research established a dynamic model of the robot and achieved a relatively accurate prediction of the robot's milling stability through time domain, frequency domain and other methods. However, its shortcoming is that it did not take into account the impact of the robot's posture changes on the modal parameters and stability state, and it is difficult to establish an accurate stability prediction model in large spaces and complex processing postures. Summary of the invention

[0007] Purpose of the invention: The technical problem to be solved by the present invention is to provide a method for obtaining a stability lobe diagram of a milling robot with variable posture, in view of the deficiencies of the prior art. Considering the influence of the robot posture change on the joint stiffness and damping, a robot variable posture milling dynamics model is established based on the multi-body system transfer matrix method, the system natural frequency and the tool tip frequency response function are predicted, and the vibration stability domain is analyzed and calculated in combination with the milling force model to draw a stability lobe diagram of the milling robot with variable posture. To achieve the above purpose, the present invention adopts the following technical solutions:

[0008] A method for obtaining a milling robot's variable posture stability lobe map comprises the following steps:

[0009] Step 1: Establish a robot dynamics model based on the multi-body system transfer matrix method, and at the same time establish a cutting force model of the milling process to form a milling robot dynamics simulation system;

[0010] Step 2: Standardize the robot milling workspace into a large cube, set the volume of the large cube to P, and further divide the space into Z small cube grids with side length h, Z = P / h 3 , Z is an integer, the grid number is i, i = 1, 2, 3...Z;

[0011] Step 3, set the initial mesh side length to h and the size change step length, and set the number of mesh divisions to r, where r = 0, 1, 2..., r = 0 represents the initial state when no mesh is divided;

[0012] Step 4: Control the center point of the robot end tool to be in the i-th grid to be sampled, obtain the joint angle under the corresponding posture of the grid to be sampled by calculating the inverse solution of the robot kinematics, refer to the robot joint stiffness and damping identification method published in the academic journal Applied Mathematical Modelling 2022, 105: 114-136 by Cui Guangyu et al., perform modal testing on the robot, and identify the robot joint stiffness matrix K in the grid to be sampled. i and the damping matrix C i ;

[0013] Step 5: transform the joint angle and joint stiffness matrix K i and the damping matrix C i Input the milling robot dynamics simulation system to calculate the system's natural frequency and tool tip frequency response function;

[0014] Step 6: Establish a flutter stability domain prediction model based on the zero-order frequency domain method, draw the stability lobe diagram of the i-th grid to be sampled, and record the critical cutting depth of the absolute stability zone;

[0015] Step 7: Count and record the average value of the critical cutting depth of all grid stability lobe diagrams CAr and maximum CM r , when r = 0, there is CA r =0.01 and CM r =0.01;

[0016] Step 8, reduce the grid size by X1 (usually 50 mm), re-divide the grid, and repeat steps 4 to 7;

[0017] Step 9: When r>0, calculate the average change rate of the critical cutting depth for the rth and r-1th times x=|CA r -CA r-1 | / CA r And the maximum rate of change y = |CM r -CM r-1 | / CM r When x and y are both less than X2 (generally 10%), output the stability lobe diagram under the current grid division; otherwise, repeat step 8 until the condition that x and y are both less than X2 is met, and output the variable attitude stability lobe diagram.

[0018] Step 1 includes:

[0019] Step 1-1, establishing a cutting force model for the milling process;

[0020] Step 1-2, measure the dynamic parameters of each component, splice the transfer matrix of each component, derive the transfer equation of the robot dynamic model, substitute the boundary conditions to solve the characteristic equation, and obtain the inherent vibration properties;

[0021] Step 1-3, according to the decoupled robot dynamic modal equation, the cutting force model of the milling process is input, and the vibration displacement in the physical space is obtained based on the modal superposition principle to complete the robot dynamic model establishment step;

[0022] Step 1-4: Input the cutting force model of the milling process into the robot dynamics model to form a milling robot dynamics simulation system, perform Fourier transform on the milling force and vibration displacement, and calculate the tool tip frequency response function.

[0023] Step 1-1 includes: The cutting force model of the instantaneous milling process of the robot with two degrees of freedom is expressed as:

[0024]

[0025] Among them, dF t 、dF r They represent the tangential milling force element and the radial milling force element, respectively, K tc , K rc Denote the tangential force coefficient and radial force coefficient respectively, a c is the milling thickness, da is the axial milling microelement, Kte , K re They represent the tangential edge force coefficient and the radial edge force coefficient respectively, and ds is the microelement of the cutting edge length;

[0026] The radial force F r , tangential force F t Projected to the workpiece coordinate system, the instantaneous milling force is obtained:

[0027]

[0028] Among them, F x 、F y They represent the milling force along the x direction and the milling force along the y direction in the workpiece coordinate system, φ is the contact angle of the milling tool, and φ ex ,φ st They represent the milling cut-out angle and cut-in angle, respectively, N is the number of milling tool edges, a is the axial cutting depth, and f is the feed per tooth;

[0029] Integrate equation (2) to obtain the average milling force in the x direction: Average milling force in y direction

[0030]

[0031] The average milling force under different milling parameters is obtained by orthogonal milling test, and K is identified by combining formula (3) tc , K rc .

[0032] Step 1-2 includes: taking an industrial six-degree-of-freedom serial robot as an example, each link of the robot is connected by gears at the joint, and the end flange carries an electric spindle and a tool capable of milling; the robot (the general embodiment of various robots is KUKA-KR210) is divided into body elements and hinge elements as a whole, numbered 0 to 14, among which No. 2, 4, 6, 8, 10, and 12 are robot body links, which are modeled as rigid body elements; No. 14 milling tool has poor rigidity and is modeled as a spatial vibration Euler-Bernoulli beam element; each body element has an independent body element that moves with the element. Coordinate system; 1, 3, 5, 7, 9, 11, 13 are 6-DOF elastic hinge elements added to the robot joints, which have linear elasticity along the x, y, z axes of the joint coordinate system and torsion elasticity around the x, y, z axes. The x, y, z axes of the hinge element coordinate system are consistent with the direction of the robot base coordinate system; Element 0 is the robot base, fixed to the ground, and is regarded as the fixed boundary of the robot dynamics model; Element 14 is the end tool, which is regarded as the free boundary of the system; The entire robot dynamics model is a tree structure, and its transmission direction is from the root 0 to the tip 14. The total transmission equation of the model is:

[0033] Z 14,0 =U1-14 Z 0,1 (4)

[0034] in, is the body element state vector, Z 0,1 , Z 14,0 They represent the state vector of the element connected to the earth and the state vector of the element connected to the free boundary, respectively. 1-14 is the total transfer matrix, which is obtained by concatenating the transfer matrices of each element in sequence. X, Y, and Z represent the linear displacements in the modal space in the x, y, and z directions in the body element coordinate system, respectively. x , Θ y , Θ z They represent the angular displacements in the x, y, and z directions of the modal space in the body element coordinate system, respectively. x 、M y 、M z They represent the internal moments in the modal space in the x, y, and z directions in the body element coordinate system, respectively, and Q x , Q y , Q z They represent the corresponding forces in the modal space in the x, y, and z directions in the body element coordinate system, respectively; the subscripts b and j are the numbers of the body element and the hinge element, respectively, and the number of the fixed boundary and the free boundary is 0;

[0035] The transfer matrices of each component of the robot are concatenated to obtain the total transfer matrix of the robot dynamics model:

[0036] U 1-14 =U 14 U 13 U 12 U 11 U 10 U9U8U7U6U5U4U3U2U1 (5)

[0037] Among them, U is the transfer matrix of each component in the model, U 14 is the transfer matrix of the 14th element;

[0038] According to the chain derivation rule of the transfer matrix method of multi-body system, element 1 is connected to the fixed boundary, and the linear displacement and angular displacement of element 1 are 0; element 14 is connected to the free boundary, and the internal force and internal moment of element 14 are 0. There are the following boundary conditions:

[0039]

[0040]

[0041] For rigid elements with b = 2, 4, 6, 8, 10, 12, C b is the mass center of the connecting rod numbered b, and the rigid element transfer matrix Ub At input point I b is the origin, O b The self-connection system for the output point is described as:

[0042]

[0043] In the formula,

[0044]

[0045]

[0046] Among them, I3 is the identity matrix of 3 rows and 3 columns, O 3×3 is a 3-row 3-column zero matrix, m b is the mass of the bth individual component, is the output point O in the coordinate system of the bth individual element b Relative to input point I b The coordinate antisymmetric matrix of is the center of mass point C in the coordinate system of the bth individual element b Relative to the output point O b The coordinate antisymmetric matrix of is the center of mass point C in the coordinate system of the bth individual element b Relative to input point I b The coordinate antisymmetric matrix, ω p is the p-th order natural frequency of the system, J b,I is the rigid body relative to the input point I b The inertia matrix, J b,x , J b,y , J b,z The input point I b The moment of inertia in the x, y, and z directions in the component coordinate system at the origin, J b,xy , J b,xz , J b,yz are the moments of inertia in the x, y, and z directions in space, respectively, (g b1 ,g b2 ,g b3 ) is the output point coordinate, (c b1 ,c b2 ,c b3 ) are the coordinates of the centroid;

[0047] For beam element b=14, the transfer matrix is:

[0048]

[0049] In the formula,

[0050] u 1,1 =u 10,10= cosβ x x b ,u 10,1 = β x EA b sinβ x x b ,u 4,4 = u 7,7 = cosγ θx x b , u 6,2 = u 11,9 = λ y V(λ y x b ), u 5,3 = u 12,8 = -λ z V(λ z x b ), u 8,3 = u 12,5 = -EI y λ z 2 U(λ z x b ), u 2,2 = u 6,6 = u 9,9 = u 11,11 = S(λ y x b ),u 8,5 = -EI y λ z V(λ z x b ),u 9,6 = EI z λ y V(λ y x b ),u 3,3 = u 5,5 = u 8,8 = u 12,12 = S(λ z x b ),

[0051] Among them, the calculation formulas for the intermediate variables λ y , γ θx , β x , λ z are as follows:

[0052]

[0053] Among them, the calculation formulas of the intermediate function variables S(u), T(u), U(u), and V(u) are:

[0054]

[0055] Where u is the function independent variable;

[0056] Among them, the transfer matrix U 14 The parameters in the self-joined coordinate system I 14 In the description, the subscript b indicates the component number, is the linear mass density, x b represents the length of the beam, ρ is the density of the beam, E is the elastic modulus, G is the shear modulus, A b is the cross-sectional area, I x is the moment of inertia of the cross section about the x-axis, I y is the moment of inertia of the cross section about the y-axis, I z is the moment of inertia of the cross section about the z-axis, J p is the polar moment of inertia of the cross section about the x-axis;

[0057] The transfer matrix parameters of all rigid elements and beam elements are obtained by measuring with 3D CAD software and do not change with the mesh changes;

[0058] For hinge elements j = 1, 3, 5, 7, 9, 11, and 13, the transfer matrix U j for:

[0059]

[0060] In the formula,

[0061]

[0062] Among them, K j Represents the linear stiffness coefficient matrix of the hinge element numbered j, K′ j Represents the torsional stiffness coefficient matrix of the hinge element numbered j, K″ j,x , K″ j,y , K″ j,z They represent the linear elastic stiffness of the hinge element numbered j in the x, y, and z directions, respectively, K′ j,x , K′ j,y , K′ j,z Respectively represent the torsional elastic stiffness in the x, y, and z directions, and the joint stiffness matrix K in the i-th grid i Contains the linear stiffness and torsional stiffness information of all hinge elements under the grid, which needs to be obtained through experimental identification;

[0063] Substituting the boundary conditions into the total transfer equation further derives the characteristic equation of the system:

[0064] det[U 1-14 (ω p )]=0 (11)

[0065] Among them, det[u] represents the determinant function of u, u is the function independent variable, ω p is the p-th order natural frequency of the robot; solving equation (11) yields the natural frequencies of each order of the robot.

[0066] Steps 1-3 include: performing modal tests to determine the damping matrix C under the i-th grid i After parameter identification and decoupling, the modal damping ratio is obtained and the vibration displacement response of the robot end is calculated according to formula (12):

[0067]

[0068] Among them, p is the p-th order modal damping ratio of the system, q p (t) are the pth order generalized displacement, generalized velocity and generalized acceleration, V p is the p-th order augmented feature vector, M p = <MV p ,V p > is the p-th order modal mass of the system, M is the system mass matrix, and f is the external force array.

[0069] Steps 1-4 include: using the cutting force model of the milling process as the external excitation input, superimposing the modal responses of each order based on the modal superposition principle, calculating the physical space vibration displacement output, performing Fourier transform on the input force and output displacement, respectively, and solving the tool tip frequency response function G(iω) and the chatter frequency ω d , the frequency response function matrix is:

[0070]

[0071] Among them, G xx (iω), G yy (iω) are the frequency response functions in the x direction and the y direction respectively, and i is an imaginary unit.

[0072] Step 5 includes: when the system vibrates, at the vibration frequency ω d The characteristic equation at is:

[0073]

[0074] Among them, I2 is a 2-row 2-column identity matrix, t is the time variable, and A0 is the first term of the Fourier series expansion of the total cutting force direction coefficient matrix that varies with time;

[0075] The eigenvalue Λ of the characteristic equation is:

[0076]

[0077] Where T is the cutting cycle, Λ R With Λ I are the real and imaginary parts of the system eigenvalues, respectively;

[0078] The eigenvalues ​​satisfy the following relationship:

[0079] b1Λ 2 +b2Λ+1=0 (16)

[0080] In the formula, b1 and b2 are intermediate parameters.

[0081] b1=G xx (iω d )G yy (iω d )(α xx α yy -α xy α yx ), b2=α xx G xx (iω d )+α yy G yy (iω d ),

[0082] Among them, α xx , α yy , α xy , α yx is the milling direction coefficient;

[0083] Rewrite equation (15) and solve equation (16) to obtain the critical milling depth a of the two-degree-of-freedom milling system: lim , for actual milling, the axial cutting depth of the tool is a real number, so we get:

[0084]

[0085] In the formula, κ is the intermediate parameter,

[0086] According to the relationship between the spindle speed n and the period T, the speed n corresponding to the axial cutting depth calculated by formula (17) in the lobe diagram is obtained as:

[0087]

[0088] Where, k = 1, 2, 3...10, is the number of lobes in the milling stability lobe diagram;

[0089] The number of lobes is selected according to the required speed range in the actual process, and the stability lobe diagram is drawn with the spindle speed as the horizontal axis and the critical axial milling depth as the vertical axis.

[0090] In step 3, the initial size of the grid and the change step size can be determined according to the actual workspace of the robot.

[0091] In step 4, the joint stiffness and damping of the robot change with the posture change. When the end moves in the same grid, the posture change is small. At this time, it is considered that the stiffness and damping are approximately unchanged in the grid.

[0092] In step 4, when performing a single-grid intra-modal test, the robot posture should be selected based on the robot posture during actual milling processing to identify the joint stiffness and damping.

[0093] The beneficial effects of the present invention are:

[0094] 1. The milling robot variable posture dynamics model proposed by the present invention, which is established by meshing the workspace, takes into account the changes in joint stiffness and damping caused by the robot posture changes, and realizes the modal characteristics and dynamic response (displacement, velocity, acceleration) of the robot under different postures. The model can save the test time cost while ensuring the required accuracy;

[0095] 2. The method for acquiring the variable posture stability lobe diagram of the milling robot proposed in the present invention takes into account the problem of variable milling postures of the robot. Compared with the existing stability prediction method, it can establish a variable posture milling stability prediction model more accurately and efficiently;

[0096] 3. The variable posture stability lobe diagram drawn by the present invention can provide appropriate milling process parameters in real time under different postures for robot milling operations in a large-scale complex workspace, effectively improving the robot milling processing stability and reducing chatter. BRIEF DESCRIPTION OF THE DRAWINGS

[0097] The present invention will be further described in detail below in conjunction with the accompanying drawings and specific embodiments, and the above and / or other advantages of the present invention will become more clear.

[0098] Figure 1 It is a flow chart of the method of the present invention.

[0099] Figure 2 It is a topological diagram of the dynamic model of the milling robot of the present invention.

[0100] Figure 3It is a schematic diagram of the gridded stiffness and damping identification process of the robot working space of the present invention.

[0101] Figure 4 It is a stability lobe diagram drawn by the robot of the present invention when the robot is in posture Φ1=(0,-90,90,0,0,180)°.

[0102] Figure 5 It is a stability lobe diagram drawn by the robot of the present invention when the robot is in posture Φ2=(0,-71,90,0,0,180)°.

[0103] Figure 6 It is a stability lobe diagram drawn by the robot of the present invention when the robot is in posture Φ3=(-11, -58, 95, -96, -95, 128)°. DETAILED DESCRIPTION

[0104] like Figure 1 The process shown in the figure takes the KUKA-210 milling robot as an example. The method for obtaining the stability lobe diagram of the milling robot according to the present invention comprises the following steps:

[0105] 1. Establish a milling robot dynamics simulation system;

[0106] Firstly, the cutting force model of the milling process is established to complete the identification of the milling force coefficient.

[0107] The cutting force model of the robot's two-degree-of-freedom instantaneous milling process can be expressed as:

[0108]

[0109] Where, dF t 、dF r They represent the tangential milling force element and the radial milling force element, respectively, K tc , K rc Denote the tangential force coefficient and radial force coefficient respectively, a c is the milling thickness, da is the axial milling microelement, K te , K re They represent the tangential edge force coefficient and the radial edge force coefficient respectively, and ds is the microelement of the cutting edge length;

[0110] The radial force F r , tangential force F t Projected to the workpiece coordinate system, the instantaneous milling force is obtained:

[0111]

[0112] In the formula, F x 、F yrepresents the milling force along the x and y directions in the workpiece coordinate system, φ is the contact angle of the milling tool, and φ ex ,φ st They represent the milling cut-out angle and cut-in angle, respectively, N is the number of milling tool edges, a is the axial cutting depth, and f is the feed per tooth;

[0113] Integrate equation (2) to obtain the average milling force in the x and y directions:

[0114]

[0115] The average milling force under different milling parameters is obtained by orthogonal milling test, and K is identified by combining formula (3) tc , K rc .

[0116] Secondly, a rigid-flexible coupling dynamic model of the robot-tool milling system is established to realize the calculation of the tool tip frequency response function.

[0117] Taking an industrial six-DOF serial robot as an example, the schematic diagram of the robot dynamics model is as follows: Figure 2 shown. Figure 2 Where o is the origin of the robot base coordinate system, x, y, z are the base coordinate axes, x, y, z with subscripts are the coordinate axes of the connection system of each component, and I, O, C, K′, K″, U with subscripts are the variables defined in the transfer matrix of each component.

[0118] The robot's connecting rods are connected by gears at the joints, and the end flange carries the electric spindle and tool capable of milling. The robot is divided into body elements and hinge elements, numbered 0 to 14, among which 2, 4, 6, 8, 10, and 12 are robot body connecting rods, which are modeled as rigid body elements. The 14th milling tool has poor rigidity and is modeled as a spatial vibration Euler-Bernoulli beam element. Each body element has an independent body coordinate system that moves with the element. 1, 3, 5, 7, 9, 11, and 13 are 6-DOF elastic hinge elements added to the robot joints, which have linear elasticity along the x, y, and z axes and torsional elasticity around the x, y, and z axes. The x, y, and z axes of the hinge element coordinate system are consistent with the robot base coordinate system. Element 0 is the robot base, fixed to the ground, and is regarded as the fixed boundary of the robot dynamics model. Element 14 is the end tool, which is regarded as the free boundary of the system. The entire robot dynamics model is a tree structure, and its transmission direction is from the root 0 to the tip 14. The total transmission equation of the model is:

[0119] Z 14,0 =U 1-14 Z 0,1 (4)

[0120] in, is the body element state vector, Z0,1 , Z 14,0 They represent the state vector of the element connected to the earth and the state vector of the element connected to the free boundary, respectively. 1-14 is the total transfer matrix, which is obtained by concatenating the transfer matrices of each element in sequence. X, Y, and Z represent the linear displacements in the modal space in the x, y, and z directions in the body element coordinate system, respectively. x , Θ y , Θ z They represent the angular displacements in the x, y, and z directions of the modal space in the body element coordinate system, respectively. x 、M y 、M z They represent the internal moments in the modal space in the x, y, and z directions in the body element coordinate system, respectively, and Q x , Q y , Q z They represent the corresponding forces in the modal space in the x, y, and z directions in the body element coordinate system, respectively; the subscripts b and j are the numbers of the body element and the hinge element, respectively, and the number of the fixed boundary and the free boundary is 0;

[0121] The transfer matrices of each component of the robot are concatenated to obtain the total transfer matrix of the robot dynamics model:

[0122] U 1-14 =U 14 U 13 U 12 U 11 U 10 U9U8U7U6U5U4U3U2U1 (5)

[0123] Among them, U is the transfer matrix of each component in the model, U 14 is the transfer matrix of the 14th element;

[0124] According to the chain derivation rule of the transfer matrix method of multi-body system, element 1 is connected to the fixed boundary, and the linear displacement and angular displacement of element 1 are 0; element 14 is connected to the free boundary, and the internal force and internal moment of element 14 are 0. There are the following boundary conditions:

[0125]

[0126]

[0127] For rigid elements with b = 2, 4, 6, 8, 10, 12, C b is the mass center of the connecting rod numbered b, and the rigid element transfer matrix U b At input point I b is the origin, O b The self-connection system for the output point is described as:

[0128]

[0129] In the formula,

[0130]

[0131]

[0132] Among them, I3 is the identity matrix of 3 rows and 3 columns, O 3×3 is a 3-row 3-column zero matrix, m b is the mass of the bth individual component, is the output point O in the coordinate system of the bth individual element b Relative to input point I b The coordinate antisymmetric matrix of is the center of mass point C in the coordinate system of the bth individual element b Relative to the output point O b The coordinate antisymmetric matrix of is the center of mass point C in the coordinate system of the bth individual element b Relative to input point I b The coordinate antisymmetric matrix, ω p is the p-th order natural frequency of the system, J b,I is the rigid body relative to the input point I b The inertia matrix, J b,x , J b,y , J b,z The input point I b The moment of inertia in the x, y, and z directions in the component coordinate system at the origin, J b,xy , J b,xz , J b,yz are the moments of inertia in the x, y, and z directions in space, respectively, (g b1 ,g b2 ,g b3 ) is the output point coordinate, (c b1 ,c b2 ,c b3 ) are the coordinates of the centroid;

[0133] For beam element b=14, the transfer matrix is:

[0134]

[0135] In the formula,

[0136] u 1,1 =u 10,10 =cosβ x x b ,u 10,1 =β x EA b sinβx x b , about 9,2 =u 11,6 =-NO z λ y 2 U(λ y x b ),u 7,4 =-γ θx (GJ p ) i you θx x b , about 6,2 =u 11,9 =λ y V(λ y x b ), about 5,3 =u 12,8 =-λ z V(λ z x b ), about 8,3 =u 12,5 =-NO y λ z 2 U(λ z x b ), about 2,2 =u 6,6 =u 9,9 =u 11,11 =S(λ y x b ),u 8,5 =-NO y λ z V(λ z x b ),u 9,6 =NO z λ y V(λ y x b )u 3,3 =u 5,5 =u 8,8 =u 12,12 =S(λ z x b )

[0137] among them,

[0138]

[0139] Intermediate variables defined for convenience of calculation

[0140] in,

[0141]

[0142] is the intermediate function variable defined, and u is the function independent variable;

[0143] Among them, the transfer matrix U 14 The parameters in the self-joined coordinate system I 14 In the description, the subscript b indicates the component number, is the linear mass density, x b represents the length of the beam, ρ is the density of the beam, E is the elastic modulus, G is the shear modulus, A b is the cross-sectional area, I x is the moment of inertia of the cross section about the x-axis, I y is the moment of inertia of the cross section about the y-axis, I z is the moment of inertia of the cross section about the z-axis, J p is the polar moment of inertia of the cross section about the x-axis;

[0144] The transfer matrix parameters of all rigid elements and beam elements are obtained by measuring with 3D CAD software and do not change with the mesh changes;

[0145] For hinge elements j = 1, 3, 5, 7, 9, 11, and 13, the transfer matrix U j for:

[0146]

[0147] In the formula,

[0148]

[0149] Among them, K j Represents the linear stiffness coefficient matrix of the hinge element numbered j, K′ j Represents the torsional stiffness coefficient matrix of the hinge element numbered j, K″ j,x , K″ j,y , K″ j,z They represent the linear elastic stiffness of the hinge element numbered j in the x, y, and z directions, respectively, K′ j,x , K′ j,y , K′ j,z Respectively represent the torsional elastic stiffness in the x, y, and z directions, and the joint stiffness matrix K in the i-th grid i Contains the linear stiffness and torsional stiffness information of all hinge elements under the grid, which needs to be obtained through experimental identification;

[0150] Substituting the boundary conditions into the total transfer equation further derives the characteristic equation of the system:

[0151] det[U 1-14 (ω p )]=0 (11)

[0152] Among them, det[u] represents the determinant function of u, u is the function independent variable, ω p is the p-th order natural frequency of the robot;

[0153] Solving equation (11) we can obtain the natural frequencies of each order of the robot;

[0154] The damping matrix C of the i-th grid is obtained by modal test. i After parameter identification and decoupling, the modal damping ratio is obtained and the vibration displacement response of the robot end is calculated according to formula (12):

[0155]

[0156] Among them, p is the p-th order modal damping ratio of the system, q p (t) are the pth order generalized displacement, generalized velocity and generalized acceleration, V p is the p-th order augmented feature vector, M p = <MV p ,V p > is the p-th order modal mass of the system, M is the system mass matrix, and f is the external force array.

[0157] The cutting force model of the milling process is used as the external excitation input. Based on the modal superposition principle, the modal responses of each order are superimposed to calculate the physical space vibration displacement output. The input force and output displacement are transformed by Fourier transform respectively to solve the tool tip frequency response function G(iω) and the chatter frequency ω d , the frequency response function matrix is:

[0158]

[0159] Among them, G xx (iω), G yy (iω) are the frequency response functions in the x direction and the y direction respectively, and i is an imaginary unit.

[0160] 2. Divide the working grid and draw the variable posture stability leaflet diagram;

[0161] like Figure 3 The figure shows the robot workspace grid division diagram. Figure 3 where i is the grid number, K i is the joint stiffness matrix of the i-th grid, C iis the damping matrix of the i-th grid, the robot milling workspace is standardized as a large cube with a side length of 2000 mm, and the entire space is further divided into 125 small cube grids with a side length of 400 mm, and the grid numbers are i = 1, 2, 3 ... 125. The initial grid side length is set to h = 400, the size change step is 50, and the number of grid divisions is set to r, where r = 0, 1, 2 ..., r = 0 represents the initial state when the grid is not divided.

[0162] The center point of the tool at the end of the robot is controlled to be in the 1st to 125th grids to be sampled. The joint rotation angle of the robot posture corresponding to the grid is obtained by calculating the inverse kinematics solution. At the same time, the robot is modally tested to identify the joint stiffness and damping of the robot in the grid. The rotation angle, stiffness, and damping information are input into the dynamics simulation system to calculate the system natural frequency and the tool tip frequency response function.

[0163] When the system vibrates, at the vibration frequency ω d The characteristic equation at is:

[0164]

[0165] Among them, I2 is a 2-row 2-column identity matrix, t is the time variable, and A0 is the first term of the Fourier series expansion of the total cutting force direction coefficient matrix that varies with time;

[0166] The eigenvalue Λ of the characteristic equation is:

[0167]

[0168] Where T is the cutting cycle, Λ R With Λ I are the real and imaginary parts of the system eigenvalues, respectively;

[0169] The eigenvalues ​​satisfy the following relationship:

[0170] b1Λ 2 +b2Λ+1=0 (16)

[0171] In the formula, b1 and b2 are intermediate parameters.

[0172] b1=G xx (iω d )G yy (iω d )(α xx α yy -α xy α yx ), b2=α xx G xx (iω d )+α yy Gyy (iω d ),

[0173] Among them, α xx , α yy , α xy , α yx is the milling direction coefficient;

[0174] Rewrite equation (15) and solve equation (16) to obtain the critical milling depth a of the two-degree-of-freedom milling system: lim , for actual milling, the axial cutting depth of the tool is a real number, so we get:

[0175]

[0176] In the formula, κ is the intermediate parameter,

[0177] According to the relationship between the spindle speed n and the period T, the speed n corresponding to the axial cutting depth calculated by formula (17) in the lobe diagram is obtained as:

[0178]

[0179] Where, k = 1, 2, 3...10, is the number of lobes in the milling stability lobe diagram;

[0180] The number of lobes is selected according to the required speed range in the actual process, and the stability lobe diagram is drawn with the spindle speed as the horizontal axis and the critical axial milling depth as the vertical axis.

[0181] Draw the stability lobe diagrams of the 1st to 125th grids, and count and record the average value CA of the critical cutting depth of all grid stability lobe diagrams r mm and maximum CM r mm, set CA r and CM r The initial value of is 0.01.

[0182] Reduce the grid size to 350 mm, re-divide the grid and repeat the stable lobe diagram drawing steps. Calculate the average change rate of the first and second critical cutting depths x = |CA r -CA r-1 | / CA r And the maximum rate of change y = |CM r -CM r-1 | / CM r , when x and y are both less than 10%, output the stability lobe diagram under the second grid division; otherwise, continue to reduce the grid size, repeat the stability lobe diagram drawing steps, and finally obtain the variable posture stability lobe diagram.

[0183] In this example, the grid side length is selected as 300 mm by the above method, and the stability lobe diagram of the KUKA-210 milling robot drawn by the present invention under the postures Φ1=(0,-90,90,0,0,180)°, Φ2=(0,-71,90,0,0,180)°, and Φ3=(-11,-58,95,-96,-95,128)° is as follows Figure 4 , Figure 5 , Figure 6 shown.

[0184] In a specific implementation, the present application provides a computer storage medium and a corresponding data processing unit, wherein the computer storage medium can store a computer program, and when the computer program is executed by the data processing unit, the invention content of the method for obtaining a lobe diagram of a milling robot with variable posture stability provided by the present invention and some or all of the steps in each embodiment can be executed. The storage medium can be a disk, an optical disk, a read-only memory (ROM) or a random access memory (RAM), etc.

[0185] Those skilled in the art can clearly understand that the technical solutions in the embodiments of the present invention can be implemented by means of computer programs and their corresponding general hardware platforms. Based on this understanding, the technical solutions in the embodiments of the present invention can be essentially or partly contributed to the prior art in the form of a computer program, i.e., a software product, which can be stored in a storage medium and includes several instructions for enabling a device including a data processing unit (which can be a personal computer, a server, a single-chip microcomputer, a MUU or a network device, etc.) to execute the methods described in various embodiments of the present invention or certain parts of the embodiments.

[0186] The present invention provides a method for obtaining a lobe diagram of a milling robot with variable posture stability. There are many methods and ways to implement the technical solution. The above is only a preferred embodiment of the present invention. It should be pointed out that for ordinary technicians in this technical field, several improvements and modifications can be made without departing from the principle of the present invention. These improvements and modifications should also be regarded as the protection scope of the present invention. All components not specified in this embodiment can be implemented using existing technologies.

Claims

1. A method for obtaining a milling robot's attitude stability lobe map, characterized in that: The following steps are involved: Step 1: Establish a robot dynamics model based on the multi-body system transfer matrix method, and at the same time establish a cutting force model of the milling process to form a milling robot dynamics simulation system; Step 2: Standardize the robot milling workspace into a large cube, set the volume of the large cube to P, and further divide the space into Z small cube grids with side length h, Z = P / h 3 , Z is an integer, the grid number is i, i = 1, 2, 3...Z; Step 3, set the initial mesh side length to h and the size change step length, and set the number of mesh divisions to r, where r = 0, 1, 2..., r = 0 represents the initial state when no mesh is divided; Step 4: Control the center point of the robot end tool to be in the i-th grid to be sampled, obtain the joint angle of the grid to be sampled by calculating the inverse solution of the robot kinematics, perform modal testing on the robot, and identify the robot joint stiffness matrix K in the grid to be sampled. i and the damping matrix C i ; Step 5: transform the joint angle and joint stiffness matrix K i and the damping matrix C i Input the milling robot dynamics simulation system to calculate the system's natural frequency and tool tip frequency response function; Step 6: Establish a flutter stability domain prediction model based on the zero-order frequency domain method, draw the stability lobe diagram of the i-th grid to be sampled, and record the critical cutting depth of the absolute stability zone; Step 7: Count and record the average value of the critical cutting depth of all grid stability lobe diagrams CA r and maximum CM r , when r = 0, there is CA r =0.01 and CM r =0.01; Step 8, reduce the grid size by X1, re-divide the grid, and repeat steps 4 to 7; Step 9: When r>0, calculate the average change rate of the critical cutting depth for the rth and r-1th times x=|CA r -CA r-1 | / CA r And the maximum rate of change y = |CM r -CM r-1 | / CM r , when x and y are both less than X2, output the stability lobe diagram under the current grid division; otherwise, repeat step 8 until the condition that x and y are both less than X2 is met, and output the variable posture stability lobe diagram.

2. The method according to claim 1, characterized in that Step 1 includes: Step 1-1, establishing a cutting force model for the milling process; Step 1-2, measure the dynamic parameters of each component, splice the transfer matrix of each component, derive the transfer equation of the robot dynamic model, substitute the boundary conditions to solve the characteristic equation, and obtain the inherent vibration properties; Step 1-3, according to the decoupled robot dynamic modal equation, the cutting force model of the milling process is input, and the vibration displacement in the physical space is obtained based on the modal superposition principle to complete the robot dynamic model establishment step; Step 1-4: Input the cutting force model of the milling process into the robot dynamics model to form a milling robot dynamics simulation system, perform Fourier transform on the milling force and vibration displacement, and calculate the tool tip frequency response function.

3. The method according to claim 2, characterized in that Step 1-1 includes: The cutting force model of the instantaneous milling process of the robot with two degrees of freedom is expressed as: Among them, dF t 、dF r They represent the tangential milling force element and the radial milling force element, respectively, K tc , K rc Denote the tangential force coefficient and radial force coefficient respectively, a c is the milling thickness, da is the axial milling microelement, K te , K re They represent the tangential edge force coefficient and the radial edge force coefficient respectively, and ds is the microelement of the cutting edge length; The radial force F r , tangential force F t Projected to the workpiece coordinate system, the instantaneous milling force is obtained: Among them, F x 、F y They represent the milling force along the x direction and the milling force along the y direction in the workpiece coordinate system, φ is the contact angle of the milling tool, and φ ex ,φ st They represent the milling cut-out angle and cut-in angle, respectively, N is the number of milling tool edges, a is the axial cutting depth, and f is the feed per tooth; Integrate equation (2) to obtain the average milling force in the x direction: Average milling force in y direction The average milling force under different milling parameters is obtained by orthogonal milling test, and K is identified by combining formula (3) tc , K rc .

4. The method according to claim 3, characterized in that: Step 1-2 includes: the robot's connecting rods are connected by gears at the joints, and the end flange carries an electric spindle and tool capable of milling; the robot is divided into body elements and hinge elements, numbered 0 to 14, among which No. 2, 4, 6, 8, 10, and 12 are robot body connecting rods, which are modeled as rigid body elements; No. 14 milling tool has poor rigidity and is modeled as a spatial vibration Euler-Bernoulli beam element; each body element has an independent connected coordinate system that moves with the element; 1, 3, 5, 7, 9, 11, 13 The number 0 is a 6-DOF elastic hinge element added to the robot joint, which has linear elasticity along the x, y, and z axes of the joint coordinate system and torsional elasticity around the x, y, and z axes. The x, y, and z axes of the hinge element coordinate system are consistent with the direction of the robot base coordinate system; the number 0 element is the robot base, fixed to the ground, and is regarded as the fixed boundary of the robot dynamics model; the number 14 element is the end tool, which is regarded as the free boundary of the system; the entire robot dynamics model is a tree structure, and its transmission direction is from the root 0 to the tip 14. The total transmission equation of the model is: WITH 14,0 =U 1-14 WITH 0,1 (4) in, is the body element state vector, Z 0,1 , Z 14,0 They represent the state vector of the element connected to the earth and the state vector of the element connected to the free boundary, respectively. 1-14 is the total transfer matrix, which is obtained by concatenating the transfer matrices of each element in sequence. X, Y, and Z represent the linear displacements in the modal space in the x, y, and z directions in the body element coordinate system, respectively. x , Θ y , Θ z They represent the angular displacements in the x, y, and z directions of the modal space in the body element coordinate system, respectively. x 、M y 、M z They represent the internal moments in the modal space in the x, y, and z directions in the body element coordinate system, respectively, and Q x , Q y , Q z They represent the corresponding forces in the modal space in the x, y, and z directions in the body element coordinate system, respectively; the subscripts b and j are the numbers of the body element and the hinge element, respectively, and the number of the fixed boundary and the free boundary is 0; The transfer matrices of each component of the robot are concatenated to obtain the total transfer matrix of the robot dynamics model: IN 1-14 =U 14 IN 13 IN 12 IN 11 IN 10 U9U8U7U6U5U4U3U2U1 (5) Among them, U is the transfer matrix of each component in the model, U 14 is the transfer matrix of the 14th element; According to the chain derivation rule of the transfer matrix method of multi-body system, element 1 is connected to the fixed boundary, and the linear displacement and angular displacement of element 1 are 0; element 14 is connected to the free boundary, and the internal force and internal moment of element 14 are 0. There are the following boundary conditions: For rigid elements with b = 2, 4, 6, 8, 10, 12, C b is the mass center of the connecting rod numbered b, and the rigid element transfer matrix U b At input point I b is the origin, O b The self-connection system for the output point is described as: In the formula, Among them, I3 is the identity matrix of 3 rows and 3 columns, O 3×3 is a 3-row 3-column zero matrix, m b is the mass of the bth individual component, is the output point O in the coordinate system of the bth individual element b Relative to input point I b The coordinate antisymmetric matrix of is the center of mass point C in the coordinate system of the bth individual element b Relative to the output point O b The coordinate antisymmetric matrix of is the center of mass point C in the coordinate system of the bth individual element b Relative to input point I b The coordinate antisymmetric matrix, ω p is the p-th order natural frequency of the system, J b,I is the rigid body relative to the input point I b The inertia matrix, J b,x , J b,y , J b,z The input point I b The moment of inertia in the x, y, and z directions in the component coordinate system at the origin, J b,xy , J b,xz , J b,yz are the moments of inertia in the x, y, and z directions in space, respectively, (g b1 ,g b2 ,g b3 ) is the output point coordinate, (c b1 ,c b2 ,c b3 ) are the coordinates of the centroid; For beam element b=14, the transfer matrix is: In the formula, you 1,1 =u 10,10 =cosβ x x b ,u 10,1 =b x EA b sinβ x x b , you 6,2 =u 11,9 =λ y V(λ y x b ), about 5,3 =u 12,8 =-λ z V(λ z x b ), about 8,3 =u 12,5 =-NO y λ z 2 U(λ z x b ), about 2,2 =u 6,6 =u 9,9 =u 11,11 =S(λ y x b ),u 8,5 =-NO y λ z V(λ z x b ),u 9,6 =NO z λ y V(λ y x b ), in 3,3 =in 5,5 =in 8,8 =in 12,12 =S(λ z x b ), Among them, the intermediate variable λ y , γ θx , β x , z The calculation formula is: Among them, the calculation formulas of the intermediate function variables S(u), T(u), U(u), and V(u) are: Where u is the function independent variable; Among them, the transfer matrix U 14 The parameters in the self-joined coordinate system I 14 In the description, the subscript b indicates the component number, is the linear mass density, x b represents the length of the beam, ρ is the density of the beam, E is the elastic modulus, G is the shear modulus, A b is the cross-sectional area, I x is the moment of inertia of the cross section about the x-axis, I y is the moment of inertia of the cross section about the y-axis, I z is the moment of inertia of the cross section about the z-axis, J p is the polar moment of inertia of the cross section about the x-axis; For hinge elements j = 1, 3, 5, 7, 9, 11, and 13, the transfer matrix U j for: In the formula, Among them, K j ″ represents the linear stiffness coefficient matrix of the hinge element numbered j, K j ′ represents the torsional stiffness coefficient matrix of the hinge element numbered j, K j ″ ,x , K j ″ ,y , K j ″ ,z They represent the linear elastic stiffness of the hinge element numbered j in the x, y, and z directions, respectively, K j ' ,x , K j ' ,y , K j ' ,z Respectively represent the torsional elastic stiffness in the x, y, and z directions, and the joint stiffness matrix K in the i-th grid i Contains the linear stiffness and torsional stiffness information of all hinge elements under the grid, which needs to be obtained through experimental identification; Substituting the boundary conditions into the total transfer equation further derives the characteristic equation of the system: it[U 1-14 (ω p )]=0 (11) Among them, det[u] represents the determinant function of u, u is the function independent variable, ω p is the p-th order natural frequency of the robot; solving equation (11) yields the natural frequencies of each order of the robot.

5. The method according to claim 4, characterized in that Steps 1-3 include: performing modal tests to determine the damping matrix C under the i-th grid i After parameter identification and decoupling, the modal damping ratio is obtained and the vibration displacement response of the robot end is calculated according to formula (12): Among them, p is the p-th order modal damping ratio of the system, q p (t) are the pth order generalized displacement, generalized velocity and generalized acceleration, V p is the p-th order augmented feature vector, M p = <MV p ,V p > is the p-th order modal mass of the system, M is the system mass matrix, and f is the external force array.

6. The method according to claim 5, characterized in that Steps 1-4 include: using the cutting force model of the milling process as the external excitation input, superimposing the modal responses of each order based on the modal superposition principle, calculating the physical space vibration displacement output, performing Fourier transform on the input force and output displacement, respectively, and solving the tool tip frequency response function G(iω) and the chatter frequency ω d , the frequency response function matrix is: Among them, G xx (iω), G yy (iω) are the frequency response functions in the x direction and the y direction respectively, and i is an imaginary unit.

7. The method according to claim 6, characterized in that Step 5 includes: when the system vibrates, at the vibration frequency ω d The characteristic equation at is: Among them, I2 is a 2-row 2-column identity matrix, t is the time variable, and A0 is the first term of the Fourier series expansion of the total cutting force direction coefficient matrix that varies with time; The eigenvalue Λ of the characteristic equation is: Where T is the cutting cycle, Λ R With Λ I are the real and imaginary parts of the system eigenvalues, respectively; The eigenvalues ​​satisfy the following relationship: b1Λ 2 +b2Λ+1=0 (16) In the formula, b1 and b2 are intermediate parameters. b1=G xx (iω d )G yy (iω d )(a xx a yy -a xy a yx ),b2=a xx G xx (iω d )+a yy G yy (iω d ), Among them, α xx , α yy , α xy , α yx is the milling direction coefficient; Rewrite equation (15) and solve equation (16) to obtain the critical milling depth a of the two-degree-of-freedom milling system: lim , for actual milling, the axial cutting depth of the tool is a real number, so we get: In the formula, κ is the intermediate parameter, According to the relationship between the spindle speed n and the period T, the speed n corresponding to the axial cutting depth calculated by formula (17) in the lobe diagram is obtained as: Where, k = 1, 2, 3...10, is the number of lobes in the milling stability lobe diagram; The number of lobes is selected according to the required speed range in the actual process, and the stability lobe diagram is drawn with the spindle speed as the horizontal axis and the critical axial milling depth as the vertical axis.

Images