Five-axis numerical control machine tool energy consumption modeling method and system based on energy flow footprint

Abstract

The invention belongs to the technical field of intelligent manufacturing twin modeling, and discloses a five-axis numerical control machine tool energy consumption modeling method and system based on energy flow footprint, and the method comprises the steps: building a five-axis numerical control machine tool three-dimensional geometric model, and obtaining part parameters; an electromechanical coupling dynamic model of the five-axis numerical control machine tool is built based on an SVPWM algorithm, and a machine tool machining process energy consumption model is built according to the electromechanical coupling dynamic model; a theoretical cutting force model is established based on an Altintas mechanical model, cutting loads of all axes are obtained according to a kinematical equation of the five-axis numerical control machine tool and applied to the load end of the machine tool machining process energy consumption model, and machine tool machining process energy consumption modeling is completed. According to the method, the cutting load model of each shaft of the machine tool is established, the cutting load can be estimated before actual machining, proper machining parameters are selected to meet the performance requirements of the machine tool, the energy consumption of the machine tool is reduced, and the problem that a method for obtaining the machine tool energy consumption model according to experimental measurement data is not applicable due to equipment change is solved.

Description

Technical Field

[0001] This invention belongs to the field of intelligent manufacturing twin modeling technology, and in particular relates to a method and system for energy consumption modeling of five-axis CNC machine tools based on energy flow footprint. Background Technology

[0002] With the development of science and technology and the constant changes in the international situation, the manufacturing industry, as the "ballast" for coping with various risks and challenges, has made high-quality development an important task. CNC machine tools are complex systems coupling multiple energy domains such as mechanics, electricity, hydraulics, pneumatics, and magnetism. Their energy consumption characteristics are affected by multiple factors, including component energy consumption, process parameters, and operating status, making modeling highly challenging. Currently, the closest existing technology for CNC machine tool energy consumption modeling is a statistical method based on experimental design. This technology obtains energy consumption measurement data during machine tool operation by designing specific experiments, and then fits and constructs an energy consumption model based on this data to meet basic energy consumption analysis needs.

[0003] Existing statistical methods rely heavily on scenario parameters set for specific experiments. The model parameters can only fit the energy consumption data under that scenario and lack general theoretical support that is independent of specific experimental scenarios. The multi-energy domain coupling characteristics of CNC machine tools mean that energy consumption is affected by multiple factors. Statistical methods can only capture local data correlations under specific scenarios and cannot reveal the inherent energy consumption laws under multi-energy domain coupling. They cannot adapt to the fluctuations in energy consumption characteristics caused by changes in scenarios. This method does not establish the inherent logical relationship between machine tool component parameters, dynamic characteristics and energy consumption. It only builds a model through the surface correlation of data, resulting in a lack of model transferability and an inability to cope with the modeling needs after changes in machine tools or optimization objectives.

[0004] Based on the above analysis, the problems and defects of the existing technology are as follows: when the energy consumption optimization target is adjusted or the machine tool itself is modified, the established energy consumption model will become completely invalid, and it is necessary to redesign the experiment, collect data and build the model, resulting in problems such as low modeling efficiency and poor scene adaptability. Summary of the Invention

[0005] To overcome the problems existing in related technologies, the present invention discloses an embodiment of a method and system for modeling the energy consumption of a five-axis CNC machine tool based on energy flow footprint. The technical solution is as follows:

[0006] This invention is implemented as follows: a five-axis CNC machine tool energy consumption modeling method based on energy flow footprint, comprising the following steps:

[0007] S1. Build a three-dimensional geometric model of a five-axis CNC machine tool and complete the assembly. Obtain the component parameters based on SolidWorks quality assessment, Hertz contact theory and Huth semi-empirical formula. The component parameters include moment of inertia, polar moment of inertia, joint stiffness and fastener compliance coefficient.

[0008] S2. Based on power bond graph theory and SVPWM algorithm, a five-axis CNC machine tool electromechanical coupling dynamics model is built, and an energy consumption model of the machine tool machining process is constructed based on the electromechanical coupling dynamics model.

[0009] S3. Based on the Altintas mechanical model, establish a theoretical cutting force model. Obtain the cutting load of each axis according to the kinematic equation of the five-axis CNC machine tool. Apply the cutting load of each axis to the load end of the machine tool machining process energy consumption model to complete the machine tool machining process energy consumption model.

[0010] In step S1, SolidWorks quality assessment is used to obtain the rotational inertia and polar moment of inertia of the component, including:

[0011] If the machine tool transmission system is a cylindrical rotating body, then the moment of inertia about the rotation center axis is:

[0012]

[0013] In the formula, J is the moment of inertia, m is the mass of the cylinder, and R is the radius of the cylinder;

[0014] For a torsion shaft, the angle through which the two cross sections rotate relative to each other during torsion is:

[0015]

[0016] In the formula, Let T be the torsion angle, T be the torque on the torsion shaft, l be the length of the torsion shaft, G be the shear modulus, and I be the torque on the torsion shaft. p The polar moment of inertia of the cross section;

[0017] The torsional stiffness K is:

[0018]

[0019] In the formula, K is the torsional stiffness.

[0020] In step S1, the stiffness of the joint at key locations is obtained based on Hertz contact theory, using simplified assumptions:

[0021] (1) The contacting object only produces elastic deformation and obeys Hooke's law;

[0022] (2) The load is perpendicular to the contact surface, and the friction between the contacting objects is negligible;

[0023] (3) The size of the contact surface is small compared to the radius of curvature of the contact object surface;

[0024] When two contacting objects are under the action of load Q, the contact point will expand into a contact surface. The projection of this contact surface onto the plane perpendicular to the contact normal is an ellipse. Within the contact area, the contact stress is distributed in a semi-elliptical pattern.

[0025] The contact stress and strain formulas derived from Hertz contact theory are as follows:

[0026]

[0027] In the formula, δ is the contact elastic approximation, E' is the equivalent elastic modulus, Q is the normal pressure, ∑ρ is the combined curvature of the two contacting bodies at the contact point, and K(e),m a All are Hertz coefficients;

[0028] The formula for calculating cosτ is obtained by looking up the auxiliary variable cosτ in a table:

[0029]

[0030] In the formula, ρ 11 ,ρ 12 Let ρ be the curvature of the two principal planes of the contact surface 1. 21 ,ρ 22 The curvatures of the two principal planes of the contact surface 2;

[0031] Then we have:

[0032]

[0033] In the formula, R 11 ,R 12 Let R be the radii of curvature of the two principal planes of the contact surface 1. 21 ,R 22 Let E' be the radius of curvature of the two principal planes of the contact surface 2; E' is obtained by the following formula:

[0034]

[0035] In the formula, E1 and E2 are the elastic moduli of the two contacting objects, and μ1 and μ2 are the Poisson's ratios of the two contacting objects.

[0036] To simplify the formula, let:

[0037]

[0038] Substituting into the formula, we get:

[0039] According to simplifying assumption 3, k h Approximately constant during the contact process, the contact stiffness K is obtained from the formula:

[0040]

[0041] In step S1, the fastener compliance coefficient is obtained based on Huth's semi-empirical formula;

[0042] For the fixed joint, the bolt connection stiffness between the harmonic reducer and the turntable is analyzed. The fastener compliance coefficient is calculated using Huth's semi-empirical formula. Then:

[0043]

[0044] In the formula, C is the fastener flexibility coefficient, t1 and t2 are the thickness of the connecting plate, E1 and E2 are the elastic modulus of the connected plate, E3 is the elastic modulus of the fastener, n is the fit coefficient, d is the fastener diameter, and a and b are empirical parameters.

[0045] In step S2, the electromechanical coupling dynamics model of the five-axis CNC machine tool includes power bond graph models of the spindle system, linear feed system, and rotary feed system. A state-space dynamics model is constructed using augmented bond graphs. The bond graph model of the linear feed system includes the mechanical properties of the coupling and bearing pairs, ball screws, ball screw-nut pairs, and the driven load. These mechanical properties include the mass coefficient, stiffness coefficient, damping coefficient, and converter module. The bond graph model of the rotary feed system includes the mechanical properties of the wave generator, rigid wheel, harmonic gear reducer, and rotating load.

[0046] Vector control of a three-phase permanent magnet synchronous motor is implemented based on the SVPWM algorithm, using i d =0 control method; selecting a mathematical model in the synchronous rotating coordinate system dq, the three-phase voltage is transformed from the three-phase natural coordinate system to the synchronous rotating coordinate system using the Clark transformation matrix and the Park transformation matrix, then:

[0047]

[0048] In the formula, T 3s / 2s T is the Clark coordinate transformation matrix. 2s / 2r Let θ be the Park coordinate transformation matrix. e For the mechanical angle of the motor;

[0049] The bond graph model of a three-phase permanent magnet synchronous motor in a synchronous rotating coordinate system includes the motor inductance resistance, electromagnetic torque coefficient, and motor rotor; by augmenting the bond graph, the augmented bond graph model of the motor is obtained.

[0050] Furthermore, the construction of the system dynamics model and the energy consumption model based on energy flow footprint includes:

[0051] (1) Determine the input variables and energy state variables;

[0052] (2) Construct the initial set of state equations;

[0053] (3) Simplify the initial system of equations into state-space form;

[0054] The construction of the state-space dynamic model of the linear feed system includes:

[0055] Based on the augmented bond diagram of the linear feed system and the permanent magnet synchronous motor, the system has nine energy storage elements, namely: inertial elements I2, I6, I... 12 ,I 15 ,I 20 ,I 26 and capacitive element C 18 C 23 C 26 The component numbers are determined according to the bond graph augmentation rules, proceeding from left to right and from top to bottom; among them, inertial components I2, I6, and I... 12 ,I 20 ,I 26 With C 18 C 26 There exists an integral causal relationship, I 15 With C 23 Differential causality exists;

[0056] The system's input variables are determined as follows:

[0057] U = [U d (t) U q (t) F(t)] T

[0058] In the formula, U d (t) represents the d-phase input voltage of the motor, U q (t) represents the input voltage of phase q of the motor, and F(t) represents the load force at the end of the workbench of the feed system.

[0059] The system's state variables are:

[0060] X = [p2 p6 p] 12 q 18 q 20 q 26 p 28 ] T

[0061] In the formula, p2, p6, p 12 ,p 20 ,p 28 For the generalized momentum of the corresponding numbered element, q 18 ,q 26 This refers to the generalized displacement of the corresponding numbered element;

[0062] Based on the properties of inertial elements with integral causality, we have:

[0063]

[0064] In the formula, f2, f6, f 12 ,f 20 ,f 28 For the corresponding numbered elements, I2, I6, I 12 ,I 20 ,I 28 These are the inertial parameters of the corresponding numbered components;

[0065] Based on the properties of capacitive elements with integral causality, we have:

[0066]

[0067] In the formula, e 18 ,e 26 For the potential variable of the corresponding numbered element, C 18 C 26 For the capacitive parameters of the corresponding numbered components;

[0068] Based on the existence of differential causality in inertial element I 15 The properties yield:

[0069] p 15 =f 15 ·I 15

[0070] In the formula, p 15 For the generalized momentum of the corresponding numbered element, f 15 For the flow variable corresponding to the numbered element, I 15 These are the inertial parameters of the corresponding numbered components;

[0071] From the system bond graph:

[0072]

[0073] Differentiating the above equation, we get:

[0074]

[0075] In the formula, The potential variable of the corresponding numbered element;

[0076] Based on the capacitive element C, which exhibits differential causality. 23 The properties yield:

[0077]

[0078] In the formula, For the flow variable corresponding to the numbered element, C 23 Here are the capacitive parameters for the corresponding numbered components, and TF is the converter module;

[0079] Based on resistive components R3, R7, R13 ,R 16 ,R 21 ,R 29 The properties yield:

[0080]

[0081] In the formula, e3, e7, e 13 ,e 16 ,e 21 ,e 29 For the potential variables of the corresponding numbered elements, R3, R7, R 13 ,R 16 ,R 21 ,R 29 For the resistive parameters of the corresponding numbered components, f3, f7, f 13 ,f 16 ,f 21 ,f 29 For the corresponding numbered element's flow variable;

[0082] Based on the causal relationships between the components and the direction of energy flow, we can conclude that:

[0083]

[0084] In the formula, e4,e8,e9,e 10 ,e 11 ,e 14 ,e 19 ,e 21 ,e 22 ,e 27 ,e 30 For the potential variable of the corresponding numbered element, MGY d MGY q The module of the gyroscope. f 17 ,f 19 ,f 25 ,f 27 For the corresponding numbered element's flow variable;

[0085] The simplified and integrated dynamic model of the linear feed system is as follows:

[0086]

[0087] Substituting the actual physical quantities into the system dynamics model, we get:

[0088]

[0089] In the formula, L d For d-phase inductance, i d Let R be the phase current (d). dFor the d-phase resistance, P n ω is the number of pole pairs of the motor. n L is the mechanical angular velocity of the motor. q For q-phase inductance, i q Let R be the phase q current. q Let ψ be the resistance of phase q. f For permanent magnet flux linkage, J z J is the moment of inertia of the motor rotor. lzq R is the moment of inertia of the bearing and coupling. b R is the rotor damping coefficient of the motor. lzq τ is the damping coefficient of the bearing and coupling. lzq C represents the coupling torque. lzq For the stiffness of the coupling, ω s J is the screw speed. s R is the moment of inertia of the lead screw. s τ is the damping coefficient of the lead screw. l To drive the load force, C l C represents the stiffness of the ball screw joint. s For the stiffness of the lead screw, ω l To drive the load speed, J l To drive the load mass, R l The damping coefficient for the driven load;

[0090] Let E be the sum of the quadratic forms on the left-hand side of each term in the system of differential equations, then:

[0091]

[0092] Differentiating E with respect to t, we get:

[0093]

[0094] Substitute the system of differential equations into:

[0095]

[0096] Rearranging the terms, we obtain the energy consumption equation P. x for:

[0097]

[0098] In the formula, P x E represents the input power of the linear feed system; E represents the energy stored in the system's energy storage element, with each term in the polynomial corresponding to the stored energy of the energy storage element; the remaining terms represent the power loss and processing power of the system's energy-consuming elements. This formula follows the law of conservation of energy and is used to characterize the direction and footprint of energy flow.

[0099] Construct the energy consumption equations for the remaining axes, and sum them to obtain the machine tool energy consumption equation:

[0100] P = P x +P y +P z +P b +P c

[0101] In the formula, P is the machine tool processing power, P x P represents the machining power of the X-axis feed system. y P is the machining power of the y-axis feed system. z P is the machining power of the z-axis feed system. b P is the machining power of the b-axis feed system. c The machining power of the C-axis feed system.

[0102] In step S3, establishing the theoretical cutting force model based on the Altintas mechanical model includes:

[0103] The formula for calculating the hysteresis angle ψ at an axial depth of cut z on a spiral end mill with a helix angle of β is:

[0104]

[0105] In the formula, D is the diameter of the milling cutter;

[0106] Establish a right-handed coordinate system with the feed direction of the helical end mill as the positive X-axis and the direction pointing to the spindle as the positive Z-axis; measure the instantaneous entry angle clockwise from the y-axis. When the instantaneous entry angle of the bottom point of the end mill's reference cutting edge is φ, the instantaneous entry angle of the cutting edge point located at an axial distance of z is φ-ψ.

[0107] A helical end mill with a helix angle of β, a diameter of D, and N cutting edges, has a tooth pitch angle of β. When the reference instantaneous angle of entry at the bottom of the cutting edge is φ, then the instantaneous angles of entry at the bottom of the other cutting edges are:

[0108] φ j (0)=φ+jφ p j = 0, 1, 2...(N-1)

[0109] Let the lag angle coefficient Then at the axial depth of tangency z, the hysteresis angle ψ = k β Then, the instantaneous entry angle of the cutting edge j at the axial depth of cut z is:

[0110] φ j (z)=φ+jφ p -k β z

[0111] In the formula, φ j (z) is the instantaneous entry angle of the cutting edge j at the axial depth of cut z, φ is the reference instantaneous entry angle at the bottom of the cutting edge, j is the cutting edge count, and z is the axial depth of cut;

[0112] The tangential, radial, and axial forces acting on the infinitesimal cutting edge at height dz are:

[0113]

[0114] In the formula, dF t,j (φ,z),dF r,j (φ,z),dF a,j (φ, z) represent the tangential force, radial force, and axial force acting on the infinitesimal cutting edge j, respectively. K tc ,K rc ,K ac The shear force coefficients h in the tangential, radial, and axial directions, respectively. j (φ j (z) represents the instantaneous cutting thickness, K te ,K re ,K ae denoted as the edge force coefficients in the tangential, radial, and axial directions, respectively, and dz is the height of the cutting micro-element;

[0115] Instantaneous cutting thickness h j (φ j (z))=c sinφ j (z) decomposes the force into three directions: feed, normal, and axial.

[0116]

[0117] In the formula, dF x,j (φ j (z)),dF y,j (φ j (z)),dF z,j (φ j (z) represents the feed, normal, and axial forces acting on the infinitesimal cutting edge j, respectively, and c is the milling cutter feed rate in millimeters per tooth angle;

[0118] Substituting the expressions for instantaneous cutting thickness, tangential force, radial force, and axial force into the equation:

[0119]

[0120] Therefore, the total cutting force generated by the cutting edge is:

[0121]

[0122] In the formula, z j,1 (φ j (z)) and z j,2 (φ j (z) represents the lower and upper limits of the axial cutting edge, respectively;

[0123] From φ j (z)=φ+jφ p -k β z dφ j (z)=-k β Substituting dz into the total cutting force equation, we get:

[0124]

[0125] In the formula, F x,j (φ j (z)),F y,j (φ j (z)),F z,j (φ j (z) represent the feed, normal, and axial cutting forces acting on the cutting edge j, respectively;

[0126] Based on the similarity triangle theorem and geometric relationships, the depth of cut z and the tool angle φ at that depth of cut are obtained during the inclination cutting process. j The relation for (z) is:

[0127]

[0128] In the formula, a is the deepest axial cutting depth, and R is the tool radius. The angle of inclination for cutting at an angle;

[0129] Substituting the tool angle formula for the axial depth of cut, we get:

[0130]

[0131] During climb milling, the tool approach angle is:

[0132]

[0133] In the formula, φ st Let a be the tool entry angle. e Radial depth of cut;

[0134] During climb milling, the cutter exit angle is:

[0135] φ ex =π

[0136] In the formula, φ ex Cut an angle with the tool;

[0137] When the spindle speed is n, taking the starting point of the cycle as the bottom of a certain cutting edge on the tool entering the cutting area, and after a time interval t after the bottom of that cutting edge enters the tool's entry angle, the instantaneous entry angle of the bottom of that cutting edge at this moment is:

[0138]

[0139] Then the instantaneous angle of entry at the bottom of the remaining cutting edges is:

[0140]

[0141] Then the instantaneous angle of entry of the cutting edge j at the axial depth of cut z is:

[0142]

[0143] Each cutting segment interacts with the cutting edge in five different ways to define the axial integral boundary, specifically including:

[0144] (1) If φ st <φ j (z=0)<φ ex Then z j,1 =0: If φ st <φ j (z=a)<φ ex Then z j,2 =a; if φ j (z=a)<φ st Then z j,2 =(φ+jφ) p -φ st ) / k β =(nt+jφ p ) / k β ;

[0145] (2) If φ j (z=0)>φ ex And φ j (z=a)<φ ex Then z j,1 =(φ+jφ) p -φ ex ) / k β =(nt+φ st +jφ p -φ ex ) / k β If φ j (z=a)<φ st Then z j,2 =(φ+jφ) p -φ st ) / k β =(nt+jφ p ) / k β If φ j (z=a)>φ st Then z j,2 =a;

[0146] (3) If φj (z=0)>φ ex And φ j (z=a)>φ ex or φ j (z=0)<φ st If so, the cutting edge does not cut;

[0147] Where, φ j (z=k) represents the instantaneous entry angle of the cutting edge j when the axial depth of cut is k, where k is a constant; z j,1 ,z j,2 These represent the lower and upper limits of axial cutting of cutting edge j, respectively;

[0148] Summing the cutting forces over all cutting edges yields the instantaneous total cutting force on the tool at the instantaneous entry angle φ. The cutting force acting on the milling cutter is:

[0149]

[0150] In the formula, F xj ,F yj ,F zj These represent the feed, normal, and axial cutting forces acting on the cutting edge j, respectively.

[0151] Furthermore, the cutting force parameters are obtained using the average value method, specifically including:

[0152] Since the cutting edge only cuts within the intrusion zone, the average milling force per tooth cycle is:

[0153]

[0154] In the formula, F is the average milling force per tooth cycle. q (φ) is the dynamic function of cutting force as a function of instantaneous angle of entry;

[0155] The derivation yields:

[0156]

[0157] In the formula, T is the tool rotation period, and F q (t) is the dynamic function of cutting force over time;

[0158] Therefore, in order to obtain the cutting force coefficient, a set of milling experiments were conducted using a milling cutter at different feed rates c, but the axial cutting depth a and radial cutting depth remained unchanged. The average milling force per tooth cycle was obtained, and then the average value method was used to evaluate the mechanism of the cutting force coefficient.

[0159] right Integrating, we get:

[0160]

[0161] In the formula, φ represents the average milling force per tooth cycle in the x, y, and z directions, respectively, and φ is the instantaneous angle of entry.

[0162] The entry angle and the removal angle for slotting milling are φ st =0 and φ ex =π, so the average force per tooth cycle simplifies to:

[0163]

[0164] The average cutting force is expressed as a linear function of the feed rate and the offset caused by the edge force:

[0165]

[0166] In the formula, The average shear force component per tooth cycle. This represents the average plowing force component per tooth cycle.

[0167] Calculate the average force at each feed rate, estimate the cutting force components through linear regression of the data, and obtain the cutting force coefficient:

[0168]

[0169] In the formula, These represent the average shear force components per tooth cycle in the x, y, and z directions, respectively. These represent the average plowing force components per tooth cycle in the x, y, and z directions, respectively.

[0170] Furthermore, based on the mechanical kinematics, the cutting load models for each axis are derived, including:

[0171] Perform kinematic modeling on a five-axis CNC machine tool. During the machine tool's operation, the absolute coordinate system O of the machine tool... M X M Y M Z M Let the B-axis coordinate system O remain unchanged. B X B Y B Z B The origin is fixed on the B-axis, and the origin of the coordinate system is O. B In machine tool coordinate system O M X M Y M Z M The coordinates in (L) mbx ,L mby ,L mbz When the B-axis rotates by α degrees, O B X B YB Z B Compared to O M X M Y M Z M The pose transformation matrix is:

[0172]

[0173] Since the rotation centers of the B and C rotary tables intersect, a table coordinate system O is established with the B axis as the Y-axis and the positive direction of the machine tool coordinate system's Y-axis as the positive direction; and the C axis as the Z-axis, with the direction away from the table surface as the positive direction. G X G Y G Z G It is stipulated that regardless of the rotation of axes B and C, the Y and Z axes of this coordinate system will always be collinear with the axes of axes B and C, and the origin of the coordinate system is O. G In machine tool coordinate system O B X B Y B Z B The coordinates in (L) bgx ,L bgy ,L bgz ), then O G X G Y G Z G Compared to O B X B Y B Z B The pose transformation matrix is:

[0174]

[0175] The workpiece is fixed to the worktable and rotates with the worktable around the C-axis. The workpiece coordinate system O W X W Y W Z W Fixed to the workpiece, with the origin O of the coordinate system W In the worktable coordinate system O G X G Y G Z G The coordinates in (L) gwx ,L gwy ,L gwz When the C-axis rotates by β degrees, O W X W YWZ W Compared to O G X G Y G Z G The pose transformation matrix is:

[0176]

[0177] With the tool tip as the origin, the tool feed direction as the positive X-axis, and the direction along the tool axis pointing to the spindle as the positive Z-axis, a right-handed coordinate system is established as the tool coordinate system. Based on the machining path and machining parameters, the X, Y, and Z axial cutting forces experienced by each cutting edge of the tool at that moment and position in the tool coordinate system are obtained. According to the translation theorem of forces, these three axial forces are translated to the tool tip. Ignoring the torque of these three axial forces, the X, Y, and Z axial cutting forces experienced by the tool tip at that moment and position in the tool coordinate system are obtained.

[0178] Given a system of parametric equations for the machining path with time as the independent variable, the derivative of each parametric equation with respect to time is used to obtain the tangent vector at any time point of the machining path. After normalization, the unit vector of the tool coordinate system X-axis in the workpiece coordinate system is obtained. In the machine coordinate system, the unit vector of the tool coordinate system Z-axis is always (0, 0, 1). After coordinate transformation and normalization of this unit vector, the unit vector of the tool coordinate system Z-axis in the workpiece coordinate system is obtained. The outer product of the two unit vectors is then normalized to obtain the unit vector of the tool coordinate system Y-axis in the workpiece coordinate system. The magnitude of the three-dimensional forces acting on the tool is multiplied by the unit vector and translated to the tool tip coordinates in the workpiece coordinate system to obtain the three-dimensional cutting forces acting on the tool in the workpiece coordinate system.

[0179] The load force model for each axis of a five-axis CNC machine tool involves six axis load models: the spindle, X, Y, and Z linear feed axes, and B and C rotary feed axes. The spindle load torque model is constructed by multiplying the theoretical cutting force x and y components by the tool radius. The load models for the other five axes are divided into rotary feed axis load models and linear feed axis load models according to their mechanical structure.

[0180] For the rotary feed axis, the tool tip interpolation point in the machine coordinate system is transformed into the tool tip interpolation point in the table coordinate system using a homogeneous coordinate transformation matrix. At the same time, based on the inverse kinematics method of the Jacobian matrix, the three-dimensional resultant force vector of the tool in the workpiece coordinate system is transformed into the three-dimensional resultant force vector in the table coordinate system. The resultant force vector is then multiplied by the force arm vector formed by the tool tip interpolation point and the origin of the table in the table coordinate system. The resulting J-axis component is the load torque of the B-axis, and the K-axis component is the load torque of the C-axis.

[0181] For the linear feed axis, the inverse kinematics method based on the Jacobian matrix converts the three-dimensional resultant force vector of the tool in the workpiece coordinate system into the three-dimensional resultant force vector in the machine coordinate system. This vector is then translated to the tool tip interpolation point in the machine coordinate system to obtain the resultant force of the tool in the machine coordinate system. The x-axis component obtained by decomposing this resultant force is the load force of the X-axis, the y-axis component is the load force of the Y-axis, and the z-axis component is the load force of the Z-axis.

[0182] Another objective of this invention is to provide a five-axis CNC machine tool energy consumption modeling system based on energy flow footprint. This system is used to implement the aforementioned five-axis CNC machine tool energy consumption modeling method based on energy flow footprint, comprising:

[0183] The parameter acquisition unit is configured to build a three-dimensional geometric model of a five-axis CNC machine tool and complete the assembly. It acquires component parameters based on SolidWorks quality assessment, Hertz contact theory and Huth semi-empirical formula. The component parameters include moment of inertia, polar moment of inertia, joint stiffness and fastener compliance coefficient.

[0184] The energy consumption model building unit is configured to build a five-axis CNC machine tool electromechanical coupling dynamics model based on power bond graph theory and SVPWM algorithm, and to build a machine tool machining process energy consumption model based on the electromechanical coupling dynamics model.

[0185] The load integration unit is configured to establish a theoretical cutting force model based on the Altintas mechanical model, obtain the cutting load of each axis according to the kinematic equation of the five-axis CNC machine tool, and apply the cutting load of each axis to the load end of the machine tool machining process energy consumption model to complete the machine tool machining process energy consumption model.

[0186] Combining all the above technical solutions, the beneficial effects of this invention are as follows:

[0187] First, this invention takes a five-axis CNC machine tool operating normally in a machining workshop as the analysis object and proposes a five-axis CNC machine tool energy consumption modeling method that considers the comprehensive energy flow mechanism of the machining process. This method takes the energy consumption of the spindle and each feed axis during the machining process as the analysis object, uses the energy flow electromechanical coupling dynamics model of the permanent magnet synchronous motor and ball screw transmission system of each axis as the main model, and uses the cutting force or cutting force of each axis after kinematic relationship processing as the load model to construct a relatively accurate energy consumption mechanism model based on the machining mechanism.

[0188] This invention employs a method that combines parameter identification and theoretical derivation to obtain system parameters, thereby identifying the essential factors affecting these parameters. It establishes cutting load models for each axis of the machine tool, enabling the estimation of cutting loads and selection of appropriate machining parameters before actual machining to meet machine tool performance requirements. This guides the selection of optimal parameters to reduce energy consumption, breaking away from the passive optimization mode of post-machining data back-calculation in existing technologies. Based on the electromechanical coupling dynamics model of a five-axis CNC machine tool, it constructs an energy consumption mechanism model for the machine tool machining process. Through energy flow footprints, it achieves full-path tracing of energy input—conversion—loss, clarifying the energy transfer logic under multi-energy domain coupling. This improves the problem of existing experimental statistical models failing due to equipment changes or changes in optimization objectives. It eliminates the need for extensive re-experimentation; only the corresponding energy flow unit parameters need to be updated to quickly adapt to new scenarios, breaking through the scenario-binding limitations of existing technologies at the mechanistic level.

[0189] Secondly, this invention designs a method for constructing a machining energy consumption mechanism model based on the energy flow footprint of the machine tool machining process. By utilizing the relevant mechanisms involved in the machine tool machining process, such as motor drive, mechanical motion, and cutting load, and combining the implementation conditions, operating parameters, and production parameters, a relatively accurate energy consumption mechanism model is obtained, which in turn guides the adjustment of machining parameters to reduce machine tool machining energy consumption and improve operating energy efficiency.

[0190] From the perspective of the overall technical solution, its core advantages and application value are reflected in the following aspects: First, compared with the existing data-driven "black box model," this invention makes the energy transfer process of each link traceable and quantifiable through energy flow footprints, and the physical meaning of the energy consumption results is clear and explicit, maintaining stable prediction accuracy in different working conditions and different models of five-axis machine tools. Second, the optimization is significantly more targeted. Energy flow footprints can accurately locate high-energy-consuming links such as motor no-load loss, transmission system friction loss, and cutting overload loss, clarify the influencing factors of each link, and guide users to carry out targeted optimization. Third, it has the value of large-scale promotion. The modeling logic is based on the general energy flow law of five-axis CNC machine tools. By establishing a standardized energy flow unit parameter library, it can be quickly adapted to various scenarios such as new equipment factory optimization and energy efficiency upgrade of existing equipment, without the need for customized experiments. It meets the industrial needs of batch energy efficiency optimization of workshop equipment in the context of intelligent manufacturing and has broad application prospects.

[0191] Third, this invention can significantly reduce the time and economic cost of energy consumption modeling for five-axis CNC machine tools, avoiding the ineffective investment of repeated experiments and modeling after changes in scenarios or machine tools; it helps machining enterprises to accurately optimize energy consumption, meet the green production needs under the "dual carbon" target, and enhance the market competitiveness of enterprises; the technology can be widely adapted to various five-axis CNC machine tools and different processing scenarios, has the potential for large-scale promotion, and can form commercial forms such as energy consumption modeling-related technical services and software products, creating continuous economic benefits.

[0192] This invention fills the gap in existing technologies that lack general theoretical support that is detached from specific experimental scenarios. By integrating theories from multiple fields such as three-dimensional geometric modeling, electromechanical coupling dynamics, and cutting force models through energy flow footprint, it establishes the inherent logical relationship between component parameters, dynamic characteristics, and energy consumption. This breaks through the limitation of traditional statistical methods that can only capture local data relationships and realizes cross-scenario and transferable energy consumption modeling.

[0193] For a long time, the industry has been eager to solve the problem of poor adaptability and insufficient transferability of energy consumption models caused by multi-energy domain coupling and changing scenarios of CNC machine tools. Traditional statistical methods have always been unable to break through the bottleneck of scenario binding. However, this invention has successfully achieved the effective reuse of models after scenario changes, machine tool modifications or optimization target adjustments by constructing a general modeling framework through the fusion of multiple theories, thus solving this long-standing technical problem.

[0194] There is a technical bias in the industry that CNC machine tool energy consumption modeling relies on fitting specific experimental data. Common methods construct models by collecting data through scenario-based experiments, neglecting the path of establishing a general model based on the essential characteristics and dynamic laws of components. This invention correlates component parameters, dynamic characteristics, and energy consumption based on energy flow footprints, without relying on specific experimental scenarios, verifying the feasibility of generalized energy consumption modeling and overcoming the limitations of traditional technical understanding. Attached Figure Description

[0195] The accompanying drawings, which are incorporated in and form part of this specification, illustrate embodiments consistent with this disclosure and, together with the description, serve to explain the principles of this disclosure;

[0196] Figure 1 This is a flowchart of the energy consumption modeling method for five-axis CNC machine tools based on energy flow footprint provided in this embodiment of the invention;

[0197] Figure 2 This is a schematic diagram of the energy consumption modeling method for five-axis CNC machine tools based on energy flow footprint provided in this embodiment of the invention;

[0198] Figure 3 This is a geometric model diagram of a five-axis CNC machine tool provided in an embodiment of the present invention;

[0199] In the diagram: 1 is the X-axis linear feed system, 2 is the Z-axis linear feed system, 3 is the spindle system, 4 is the C-axis rotary feed system, 5 is the Y-axis linear feed system, and 6 is the B-axis rotary feed system.

[0200] Figure 4 This is an interface diagram for obtaining the rotational inertia of components and the polar moment of inertia of cross sections provided in an embodiment of the present invention;

[0201] Figure 5 This is a stress distribution diagram of the point contact area and contact zone provided in an embodiment of the present invention;

[0202] Figure 6 This is a diagram of the energy flow model for five-axis CNC machine tool processing provided in an embodiment of the present invention;

[0203] Figure 7 This is a bonding diagram model of the spindle system provided in an embodiment of the present invention;

[0204] Figure 8 This is an augmented bond graph model diagram of the spindle system provided in an embodiment of the present invention;

[0205] Figure 9 This is a bond diagram model of the X-axis feed system provided in an embodiment of the present invention;

[0206] Figure 10 This is an augmented bond graph model diagram of the X-axis feed system provided in this embodiment of the invention;

[0207] Figure 11 This is a bond diagram model of the C-axis rotary feed system provided in an embodiment of the present invention;

[0208] Figure 12 This is an augmented bond graph model diagram of the C-axis rotary feed system provided in an embodiment of the present invention;

[0209] Figure 13 This is a bonding diagram model of a three-phase permanent magnet synchronous motor provided in an embodiment of the present invention;

[0210] Figure 14 This is an augmented bond graph model diagram of the motor provided in an embodiment of the present invention;

[0211] Figure 15 This is a vector control block diagram of a three-phase permanent magnet synchronous motor provided in an embodiment of the present invention;

[0212] Figure 16 This is a model diagram of the electromechanical coupling augmented bond graph of the X-axis feed system provided in this embodiment of the invention;

[0213] Figure 17 This is a diagram showing the interaction between the cutting edge and the workpiece provided in an embodiment of the present invention;

[0214] Figure 18 This is a flowchart of the rotary feed axis load torque calculation provided in an embodiment of the present invention;

[0215] Figure 19 This is a flowchart of the linear feed axis load force calculation provided in an embodiment of the present invention;

[0216] Figure 20 This is a comparison diagram of the milling force in the x-direction of the circular groove machining experiment provided in the embodiment of the present invention;

[0217] Figure 21This is a comparison diagram of the milling force in the y-direction of the circular groove machining experiment provided in the embodiment of the present invention;

[0218] Figure 22 This is a comparison diagram of the milling force in the z-direction of the circular groove machining experiment provided in the embodiment of the present invention;

[0219] Figure 23 This is a simulation power diagram of the X-axis linear feed system for machining S-shaped specimens provided in this embodiment of the invention;

[0220] Figure 24 This is a simulation power diagram of the spindle system for machining S-shaped specimens provided in an embodiment of the present invention. Detailed Implementation

[0221] To make the above-mentioned objects, features, and advantages of the present invention more apparent and understandable, specific embodiments of the present invention will be described in detail below with reference to the accompanying drawings. Many specific details are set forth in the following description to provide a thorough understanding of the present invention. However, the present invention can be practiced in many other ways different from those described herein, and those skilled in the art can make similar modifications without departing from the spirit of the present invention. Therefore, the present invention is not limited to the specific embodiments disclosed below.

[0222] The innovation of this invention lies in:

[0223] (1) Abandoning the traditional statistical modeling approach that relies on surface correlations of specific experimental data, this paper establishes an intrinsic logical link between the essential parameters of machine tool components, electromechanical coupling dynamic characteristics, and energy consumption, with energy flow footprint as the core. Core component attributes are obtained through SolidWorks quality assessment and Hertz contact theory, and then a dynamic model is constructed by combining power bond graphs and the SVPWM algorithm. This achieves a shift from data fitting to mechanism-driven approaches, giving the model theoretical universality and interpretability.

[0224] (2) A closed-loop modeling framework is formed by integrating technologies from multiple fields. The front end obtains component parameters accurately through 3D geometric modeling and semi-empirical formulas. The middle end builds an electromechanical coupling dynamic model based on power bond graph theory and derives an energy consumption model. The back end combines the Altintas mechanical model and the machine tool kinematic equations to calculate the cutting load and apply it to the energy consumption model. This system covers the entire process of component properties, dynamic characteristics, cutting load, and energy consumption, breaking through the limitations of traditional single technical paths that cannot cope with multi-energy domain coupling, and fully capturing the energy consumption formation mechanism of machine tools.

[0225] (3) In response to the problem of poor adaptability of traditional models due to scene binding, this invention constructs a model based on a general theoretical framework. When machine tool parts are modified, processing scenes are changed, or optimization targets are adjusted, there is no need to redesign experiments and modeling. Only the part parameters or load data need to be updated to reuse the model, which greatly improves the cross-scene adaptability and portability of the model and significantly reduces modeling costs and cycles.

[0226] Example 1, such as Figure 1 As shown in the embodiment of the present invention, the energy consumption modeling method for five-axis CNC machine tools based on energy flow footprint first constructs a three-dimensional geometric model of the machine tool based on the geometric parameters and geometric relationships of each system of the five-axis CNC machine tool to be analyzed, and obtains some parameters of each component by consulting relevant component manuals; then, SolidWorks quality assessment, Hertz contact theory and Huth semi-empirical formula are used to calculate the inertia and stiffness coefficient of each component, and the machine tool system is analyzed and simplified based on power bond graph theory to establish electromechanical coupling power bond graph models of each system; vector control of a three-phase permanent magnet synchronous motor is implemented based on the SVPWM algorithm, and an electromechanical coupling dynamic model of the five-axis CNC machine tool based on energy flow footprint is constructed; a cutting force model of the machining process is established based on the Altintas mechanical model, and a cutting load force model of each axis system is constructed based on the kinematic relationship of each axis of the machine tool; finally, the energy consumption model of the machine tool machining process is obtained to characterize the single-machine energy consumption of the machine tool under specific machining conditions. Figure 2 As shown, the method specifically includes the following steps:

[0227] S1. Build a three-dimensional geometric model of a five-axis CNC machine tool and complete the assembly. Obtain the component parameters based on SolidWorks quality assessment, Hertz contact theory and Huth semi-empirical formula. The component parameters include moment of inertia, polar moment of inertia, joint stiffness and fastener compliance coefficient.

[0228] S101: Geometric model of a five-axis CNC machine tool.

[0229] For the geometric dimensions of five-axis CNC machine tools, based on their geometric feature parameters and other information, an expandable geometric model covering the heterogeneous elements of the equipment is constructed. Currently, commonly used geometric model construction methods include direct modeling using 3D software, modeling using instrument measurement methods, and modeling using videos or images. Since the modeling object has a complete entity and drawings, SolidWorks 3D software is used to directly build the geometric model.

[0230] Using SolidWorks, 3D models of each part of the machine tool were created, and the models were then assembled. The resulting 5-axis machine tool geometric model is shown below. Figure 3 As shown in the figure. This model can be used in subsequent parameter solving processes and for building digital twin models of intelligent manufacturing production.

[0231] S102: Use SolidWorks quality assessment to obtain the rotational inertia and polar moment of inertia of the components.

[0232] To achieve dynamic simulation of the system, it is necessary to obtain the parameters of the spindle system and each feed system. Most parameters can be obtained by consulting the product manual, but some parameters, due to uncontrollable factors, need to be obtained using relevant theories and linearized to simplify the system complexity. Since the structures of various machine tool systems are similar, this invention takes the spindle, X-axis linear feed system, and C-axis rotary feed system as examples, and uses the built-in mass assessment module of SolidWorks software to obtain the rotational inertia and polar moment of inertia of the components.

[0233] In machine tool transmission systems, most components are cylindrical rotating bodies with a moment of inertia about their central axis of rotation:

[0234]

[0235] In the formula, J is the moment of inertia, m is the mass of the cylinder, and R is the radius of the cylinder;

[0236] For a torsion shaft, the angle through which the two cross sections rotate relative to each other during torsion is:

[0237]

[0238] In the formula, Let T be the torsion angle, T be the torque on the torsion shaft, l be the length of the torsion shaft, G be the shear modulus, and I be the torque on the torsion shaft. p The polar moment of inertia of the cross section;

[0239] Therefore, the torsional stiffness K can be obtained as:

[0240]

[0241] In the formula, K is the torsional stiffness.

[0242] To obtain the moment of inertia and polar moment of inertia of a component, the materials of each component are set in SolidWorks software. The moment of inertia and polar moment of inertia of the component are then queried through the Evaluation-Mass Properties dialog box. (See...) Figure 4 .

[0243] S103: Obtain the stiffness of the joint of key parts based on Hertz contact theory.

[0244] Hertz contact theory, under certain simplifying assumptions, can calculate the stress distribution on the contact surface. Its expression is simple, and the calculation results are in good agreement with reality. This invention adopts the following simplifying assumptions:

[0245] (1) The contacting object undergoes only elastic deformation, obeying Hooke's law. Under normal circumstances, all components of the machine tool operate below the elastic limit of the material, and the total plastic deformation at the contact point does not exceed one ten-thousandth of the radial diameter of the contacting object. Therefore, this assumption can be considered appropriate;

[0246] (2) The load is perpendicular to the contact surface, that is, it is assumed that the contact surface is completely smooth and the friction between the contacting objects is negligible. Since CNC machine tools are high-precision equipment, their minimum accuracy is determined by the minimum accuracy of the motion system. Therefore, the friction of the kinematic pairs at the joints of their key parts is very small. Based on this assumption, a sufficiently high accuracy can be obtained when calculating the maximum stress of the contact surface, the size of the contact surface, and the elastic approximation of the contacting objects.

[0247] (3) The size of the contact surface is small compared to the radius of curvature of the contact object surface. Strictly speaking, this assumption does not fully conform to the case of ball screws. However, after adopting this assumption, the calculation results are still relatively consistent with the experimental results. Therefore, for simplicity, this assumption is still adopted.

[0248] When two contacting objects are under a load Q, the contact point will expand into a contact surface. The projection of this contact surface onto a plane perpendicular to the contact normal is an ellipse. Within the contact area, the contact stress is distributed in a semi-elliptical pattern, as shown below. Figure 5 As shown.

[0249] The contact stress and strain formulas derived from Hertz contact theory are as follows:

[0250]

[0251] In the formula, δ is the contact elastic approximation, E' is the equivalent elastic modulus, Q is the normal pressure, ∑ρ is the combined curvature of the two contacting bodies at the contact point, and K(e),m a All are Hertz coefficients;

[0252] The formula for calculating cosτ is obtained by looking up the auxiliary variable cosτ in a table:

[0253]

[0254] In the formula, ρ 11 ,ρ 12 Let ρ be the curvature of the two principal planes of the contact surface 1. 21 ,ρ 22 The curvatures of the two principal planes of the contact surface 2;

[0255] Then we have:

[0256]

[0257] In the formula, R 11 ,R12 Let R be the radii of curvature of the two principal planes of the contact surface 1. 21 ,R 22 The radii of curvature of the two principal planes of the contact surface 2;

[0258] E' is obtained from the following formula:

[0259]

[0260] In the formula, E1 and E2 are the elastic moduli of the two contacting objects, and μ1 and μ2 are the Poisson's ratios of the two contacting objects.

[0261] To simplify the formula, let:

[0262]

[0263] Substituting into the formula, we get:

[0264]

[0265] Because the deformation generated during the contact process is small, and according to Hertz contact theory assumption 3, k h Approximately constant during the contact process, the contact stiffness K is obtained from the formula:

[0266]

[0267] The motion system of a CNC machine tool mainly consists of three linear feed systems (X-axis, Y-axis, and Z-axis) and two rotary feed systems (B-axis and C-axis), containing numerous motion engagement surfaces. Since the three linear feed systems (X-axis, Y-axis, and Z-axis) have the same structure, this invention takes the X-axis linear feed system and the B-axis rotary system as examples, considering the ball screw pair, linear guide slider pair, and bearing pair of the X-axis, and the bearing pair of the B-axis, and establishes their equivalent spring-mass model.

[0268] A ball screw pair is a helical transmission mechanism that converts linear motion into rotary motion, allowing the nut to rotate freely around the helix of the screw. Therefore, it can be considered equivalent to an axial tension spring, neglecting the influence of torsional stiffness. A linear guide slider pair achieves relative linear motion between the slider and the guide rail through ball rolling. A single slider bears force loads from two directions and torsional loads from three directions. Therefore, the rolling guide pair can be considered equivalent to four sets of tension springs evenly arranged along the guide rail direction, each set containing two tension springs along the contact normal direction. The X-axis feed system's screw uses a "one-end fixed, one-end floating" support method, while the B-axis rotary feed system's cradle uses a "one-end driven, one-end floating" support method. The fixed end uses a pair of angular contact bearings to restrict axial and radial displacement; the floating end uses a single deep groove ball bearing to restrict radial displacement. The pair of angular contact bearings at the fixed end can be considered equivalent to four radial tension springs and one axial tension spring; the deep groove ball bearing at the floating end can be considered equivalent to four radial tension springs.

[0269] Consult the product manual to obtain the product parameters of the main joint parts used in this invention, and substitute them into the above formula to calculate the stiffness of the key joint parts of the analysis object of this invention.

[0270] S104: Obtain the fastener compliance coefficient based on Huth's semi-empirical formula.

[0271] For the fixed joint, this invention needs to consider the bolt connection stiffness between the harmonic reducer and the turntable. In the engineering field, Huth's semi-empirical formula is often used to calculate the fastener compliance coefficient:

[0272]

[0273] In the formula, C is the fastener flexibility coefficient, t1 and t2 are the thickness of the connecting plate, E1 and E2 are the elastic modulus of the connected plate, E3 is the elastic modulus of the fastener, n is the fit coefficient (single shear: n = 1, double shear: n = 2), d is the fastener diameter, and a and b are empirical parameters, confirmed by Table 1.

[0274] Table 1 Empirical Parameters

[0275]

[0276] Substituting the parameters of the harmonic reducer and turntable bolt connection used in this invention into the above formula, the fastener connection stiffness can be obtained.

[0277] S2. Based on power bond graph theory and SVPWM algorithm, a five-axis CNC machine tool electromechanical coupling dynamics model is built, and an energy consumption model of the machine tool machining process is constructed based on the electromechanical coupling dynamics model.

[0278] S201: A dynamic model of electromechanical coupling of a five-axis CNC machine tool was built based on power bond graph theory and SVPWM algorithm.

[0279] As a complex multi-energy-domain coupled system, CNC machine tools have numerous energy-consuming components and complex energy flow patterns and consumption characteristics, mainly including electrical, hydraulic, pneumatic, and mechanical components. Therefore, when establishing the system's mechanistic model, it is necessary to simplify components with relatively small impact on the machine tool's energy consumption and dynamic characteristics (sensors, cables connecting moving parts). The main energy consumption is found to consist of six machining motion energy consumption systems: the spindle, three linear feed axes (X, Y, and Z axes), and two rotary feed axes (B and C axes). Based on the equipment composition of each machining motion energy consumption system, the spindle system can be broken down into the spindle motor, bearings, and spindle components; the three linear feed axes (X, Y, and Z axes) can be broken down into servo motors, couplings, bearings, ball screw nut pairs, and a load table; the two rotary feed axes (B and C axes) consist of servo motors, harmonic gear reducers, and a load table. This yields the five-axis CNC machine tool machining energy flow model as follows: Figure 6 As shown.

[0280] Ignoring some components that have a relatively small impact on machine tool energy consumption, a machine tool dynamics model is established based on bond graph theory. According to the energy flow model of a five-axis CNC machine tool, the mechanical characteristics of key components of the spindle system are shown in Table 2.

[0281] Table 2 Mechanical properties of key components of the spindle system

[0282]

[0283] Establish a bond graph model of the spindle system as follows Figure 7 As shown.

[0284] The augmentation of the bond graph is based on a bond graph that has been constructed and labeled with reasonable power flow directions. Additional information is added to make the model computable. The required additional conditions include: (1) numbering all bonds sequentially; and (2) specifying causal relationships for each bond based on potential and flow variables. By augmenting the bond graph, the augmented bond graph model of the spindle system is obtained as follows: Figure 8 As shown.

[0285] Based on the energy flow model of a five-axis CNC machine tool, since the X, Y, and Z linear feed systems have the same structure, this invention takes the X-axis feed system as an example, and obtains the mechanical characteristics of the key components of the X-axis feed system as shown in Table 3.

[0286] Table 3 Mechanical properties of key components of the X-axis feed system

[0287]

[0288] Establish the bond graph model of the X-axis feed system as follows: Figure 9 As shown, by augmenting the bond graph, the augmented bond graph model of the X-axis feed system is obtained as follows: Figure 10 As shown.

[0289] Based on the energy flow model of a five-axis CNC machine tool, since the B and C rotary feed systems have the same structure, the C-axis rotary feed system is taken as an example. The mechanical characteristics of the key components of the C-axis rotary feed system are shown in Table 4.

[0290] Table 4 Mechanical properties of key components of the C-axis rotary feed system

[0291]

[0292] Establish a bond graph model for the C-axis rotary feed system as follows: Figure 11 As shown. By augmenting the bond graph, the augmented bond graph model of the C-axis rotary feed system is obtained as follows. Figure 12 As shown.

[0293] For a three-phase permanent magnet synchronous motor, a mathematical model in the synchronous rotating coordinate system dq is chosen to facilitate controller design. The Clark transformation matrix and Park transformation matrix required to transform the three-phase voltage from the three-phase natural coordinate system to the synchronous rotating coordinate system are shown in the following equations:

[0294]

[0295] In the formula, T 3s / 2s T is the Clark coordinate transformation matrix. 2s / 2r Let θ be the Park coordinate transformation matrix. e For the mechanical angle of the motor;

[0296] Considering that the motor windings have inductance and current, and the rotor has rotational inertia, the electrical and mechanical characteristics of the key components of the motor are shown in Table 5.

[0297] Table 5 Electrical and Mechanical Characteristics of Key Components of the Motor

[0298]

[0299] Establish a bond diagram model of a three-phase permanent magnet synchronous motor in a synchronous rotating coordinate system, as follows: Figure 13 As shown. By augmenting the bond graph, the augmented bond graph model of the motor is obtained as follows. Figure 14 As shown.

[0300] Currently, common methods for vector control of traditional permanent magnet synchronous motors include i d =0 control and maximum torque-current ratio control are used. The former is mainly applicable to surface-mounted three-phase permanent magnet synchronous motors, while the latter is mainly used for built-in three-phase permanent magnet synchronous motors. Since this invention uses a surface-mounted three-phase permanent magnet synchronous motor, i is adopted. d The block diagram of the three-phase PMSM vector control based on the SVPWM algorithm with =0 control method is as follows: Figure 15 As shown.

[0301] S202: System dynamics modeling and energy consumption model establishment based on energy flow footprint.

[0302] After completing the augmented construction of the power bond graphs of each system in this invention, the system dynamics model can be constructed through the following general process:

[0303] (1) Determine the input variables and energy state variables;

[0304] (2) Construct the initial set of state equations;

[0305] (3) Simplify the initial system of equations into state-space form.

[0306] This invention takes a linear feed system as an example and establishes a dynamic model of the system in state-space form.

[0307] Based on the electromechanical coupling augmented bond diagram of the X-axis feed system, such as Figure 16 As shown. The system has nine energy storage elements, namely inertial elements I2, I6, I... 12 ,I 15 ,I 20 ,I 26 and capacitive element C 18 C 23 C 26 (The component numbers are arranged according to the bond diagram augmentation rules, from left to right and top to bottom). Among them, inertial components I2, I6, and I... 12 ,I 20 ,I 26 With C 18 C 26 There exists an integral causal relationship, I 15 With C 23 Differential causality exists.

[0308] First, determine the system's input variables as follows:

[0309] U = [U q (t) U q (t) F(t)] T

[0310] In the formula, U d (t) represents the d-phase input voltage of the motor, U q (t) represents the input voltage of phase q of the motor, and F(t) represents the load force at the end of the workbench of the feed system.

[0311] The system's state variables are:

[0312] X = [p2 p6 p] 12 q 18 q 20 q26 p 28 ] T

[0313] In the formula, p2, p6, p 12 ,p 20 ,p 28 For the generalized momentum of the corresponding numbered element, q 18 ,q 26 This refers to the generalized displacement of the corresponding numbered element;

[0314] Based on the properties of inertial elements with integral causality, we have:

[0315]

[0316] In the formula, f2, f6, f 12 ,f 20 ,f 28 For the corresponding numbered elements, I2, I6, I 12 ,I 20 ,I 28 These are the inertial parameters of the corresponding numbered components;

[0317] Based on the properties of capacitive elements with integral causality, we have:

[0318]

[0319] In the formula, e 18 ,e 26 For the potential variable of the corresponding numbered element, C 18 C 26 For the capacitive parameters of the corresponding numbered components;

[0320] Based on the existence of differential causality in inertial element I 15 The properties yield:

[0321] p 15 =f 15 ·I 15

[0322] In the formula, p 15 For the generalized momentum of the corresponding numbered element, f 15 For the flow variable corresponding to the numbered element, I 15 These are the inertial parameters of the corresponding numbered components;

[0323] From the system bond graph:

[0324]

[0325] Differentiating the above equation, we get:

[0326]

[0327] In the formula, The potential variable of the corresponding numbered element;

[0328] Similarly, based on the capacitive element C with differential causality... 23 The properties yield:

[0329]

[0330] In the formula, For the flow variable corresponding to the numbered element, C 23 Here are the capacitive parameters for the corresponding numbered components, and TF is the converter module;

[0331] Based on resistive components R3, R7, R 13 ,R 16 ,R 21 ,R 29 The properties yield:

[0332]

[0333] In the formula, e3, e7, e 13 ,e 16 ,e 21 ,e 29 For the potential variables of the corresponding numbered elements, R3, R7, R 13 ,R 16 ,R 21 ,R 29 For the resistive parameters of the corresponding numbered components, f3, f7, f 13 ,f 16 ,f 21 ,f 29 For the corresponding numbered element's flow variable;

[0334] Based on the causal relationships between the components and the direction of energy flow, we can conclude that:

[0335]

[0336] In the formula, e4,e8,e9,e 10 ,e 11 ,e 14 ,e 19 ,e 21 ,e 22 ,e 27 ,e 30 For the potential variable of the corresponding numbered element, MGY d MGY q The module of the gyroscope. f 17 ,f 19 ,f25 ,f 27 For the corresponding numbered element's flow variable;

[0337] The simplified and integrated dynamic model of the linear feed system is as follows:

[0338]

[0339] Substituting the actual physical quantities into the system dynamics model, we get:

[0340]

[0341] In the formula, L d For d-phase inductance, i d Let R be the phase current (d). d For the d-phase resistance, P n ω is the number of pole pairs of the motor. n L is the mechanical angular velocity of the motor. q For q-phase inductance, i q Let R be the phase q current. q Let ψ be the resistance of phase q. f For permanent magnet flux linkage, J z J is the moment of inertia of the motor rotor. lzp R is the moment of inertia of the bearing and coupling. b R is the rotor damping coefficient of the motor. lzq τ is the damping coefficient of the bearing and coupling. lzq C represents the coupling torque. lzq For the stiffness of the coupling, ω s J is the screw speed. s R is the moment of inertia of the lead screw. s τ is the damping coefficient of the lead screw. l To drive the load force, C l C represents the stiffness of the ball screw joint. s For the stiffness of the lead screw, ω l To drive the load speed, J l To drive the load mass, R l The damping coefficient for the driven load;

[0342] Let E be the sum of the quadratic forms on the left-hand side of each term in the system of differential equations, then:

[0343]

[0344] Differentiating E with respect to t, we get:

[0345]

[0346] Substitute the system of differential equations into:

[0347]

[0348] Rearranging the terms, we obtain the energy consumption equation P. x for:

[0349]

[0350] In the formula, P x E represents the input power of the linear feed system; E represents the energy stored in the system's energy storage element, with each term in the polynomial corresponding to the stored energy of the energy storage element; the remaining terms represent the power loss and processing power of the system's energy-consuming elements. This formula follows the law of conservation of energy and is used to characterize the direction and footprint of energy flow.

[0351] Construct the energy consumption equations for the remaining axes, and sum them to obtain the machine tool energy consumption equation:

[0352] P = P x +P y +P z +P b +P c

[0353] In the formula, P is the machine tool processing power, P x P represents the machining power of the X-axis feed system. y P is the machining power of the y-axis feed system. z P is the machining power of the z-axis feed system. b P is the machining power of the b-axis feed system. c The machining power of the C-axis feed system.

[0354] S3. Based on the Altintas mechanical model, establish a theoretical cutting force model. Obtain the cutting load of each axis according to the kinematic equation of the five-axis CNC machine tool. Apply the cutting load of each axis to the load end of the machine tool machining process energy consumption model to complete the machine tool machining process energy consumption model.

[0355] S301: Establish a theoretical cutting force model based on the Altintas mechanical model.

[0356] In actual machining, the loads on each system mainly come from the forces and torques generated by cutting. Therefore, a theoretical cutting force model based on the Altintas mechanical model is established as a load applied to the end of the machine tool table to establish a complete machine tool dynamics model.

[0357] The helix on a helical end mill causes a point on the cutting edge axis to lag behind the tool tip. Specifically, the lag angle ψ at a point with a helix angle of β and an axial depth of cut of z is calculated using the following formula:

[0358]

[0359] In the formula, D is the diameter of the milling cutter;

[0360] Establish a right-handed coordinate system with the feed direction of the helical end mill as the positive X-axis and the direction pointing to the spindle as the positive Z-axis; measure the instantaneous entry angle clockwise from the y-axis. When the instantaneous entry angle of the bottom point of the end mill's reference cutting edge is φ, the instantaneous entry angle of the cutting edge point located at an axial distance of z is φ-ψ.

[0361] A helical end mill with a helix angle of β, a diameter of D, and N cutting edges, has a tooth pitch angle of β. When the reference instantaneous angle of entry at the bottom of the cutting edge is φ, then the instantaneous angles of entry at the bottom of the other cutting edges are:

[0362] φ j (0)=φ+jφ p j = 0, 1, 2...(N-1)

[0363] Let the lag angle coefficient Then at the axial depth of tangency z, the hysteresis angle ψ = k β Then, the instantaneous entry angle of the cutting edge j at the axial depth of cut z is:

[0364] φ j (z)=φ+jφ p -k β z

[0365] In the formula, φ j (z) is the instantaneous entry angle of the cutting edge j at the axial depth of cut z, φ is the reference instantaneous entry angle at the bottom of the cutting edge, j is the cutting edge count, and z is the axial depth of cut;

[0366] The tangential, radial, and axial forces acting on the infinitesimal cutting edge at height dz are:

[0367]

[0368] In the formula, dF t,j (φ,z),dF r,j (φ,z),dF a,j (φ, z) represent the tangential force, radial force, and axial force acting on the infinitesimal cutting edge j, respectively. K tc ,K rc ,K ac The shear force coefficients h in the tangential, radial, and axial directions, respectively. j (φ j (z) represents the instantaneous cutting thickness, K te ,K re ,K ae denoted as the edge force coefficients in the tangential, radial, and axial directions, respectively, and dz is the height of the cutting micro-element;

[0369] Among them, the instantaneous cutting thickness h j (φj (z))=c sinφ j (z) decomposes the force into three directions: feed, normal, and axial.

[0370]

[0371] In the formula, dF x,j (φ j (z)),dF y,j (φ j (z)),dF z,j (φ j (z) represents the feed, normal, and axial forces acting on the infinitesimal cutting edge j, respectively, and c is the milling cutter feed rate in millimeters per tooth angle;

[0372] Substituting the expressions for instantaneous cutting thickness, tangential force, radial force, and axial force into the equation:

[0373]

[0374] Therefore, the total cutting force generated by the cutting edge is:

[0375]

[0376] In the formula, z j,1 (φ j (z)) and z j,2 (φ j (z) represents the lower and upper limits of the axial cutting edge, respectively;

[0377] From φ j (z)=φ+jφ p -k β z dφ j (z)=-k β Substituting dz into the total cutting force equation, we get:

[0378]

[0379] In the formula, F x,j (φ j (z)),F y,j (φ j (z)),F z,j (φ j (z) represent the feed, normal, and axial cutting forces acting on the cutting edge j, respectively;

[0380] Based on the similarity triangle theorem and geometric relationships, the depth of cut z and the tool angle φ at that depth of cut are obtained during the inclination cutting process. j The relation for (z) is:

[0381]

[0382] In the formula, a is the deepest axial cutting depth, and R is the tool radius. The angle of inclination for cutting at an angle;

[0383] Substituting the tool angle formula for the axial depth of cut, we get:

[0384]

[0385] During climb milling, the tool approach angle is:

[0386]

[0387] In the formula, φ st Let a be the tool entry angle. e Radial depth of cut;

[0388] During climb milling, the cutter exit angle is:

[0389] φ ex =π

[0390] In the formula, φ ex Cut an angle with the tool;

[0391] When the spindle speed is n, taking the starting point of the cycle as the bottom of a certain cutting edge on the tool entering the cutting area, and after a time interval t after the bottom of that cutting edge enters the tool's entry angle, the instantaneous entry angle of the bottom of that cutting edge at this moment is:

[0392]

[0393] Then the instantaneous angle of entry at the bottom of the remaining cutting edges is:

[0394]

[0395] Then the instantaneous angle of entry of the cutting edge j at the axial depth of cut z is:

[0396]

[0397] Each cutting segment interacts with the cutting edge in five different ways to define the axial integral boundary, specifically including:

[0398] (1) If φ st <φ j (z=0)<φ ex Then z j,1 =0: If φ st <φ j (z=a)<φ ex Then z j,2 =a; if φ j (z=a)<φ st Then zj,2 =(φ+jφ) p -φ st ) / k β =(nt+jφ p ) / k β ;

[0399] (2) If φ j (z=0)>φ ex And φ j (z=a)<φ ex Then z j,1 =(φ+jφ) p -φ ex ) / k β =(nt+φ st +jφ p -φ ex ) / k β If φ j (z=a)<φ st Then z j,2 =(φ+jφ) p -φ st ) / k β =(nt+jφ p ) / k β If φ j (z=a)>φ st Then z j,2 =a;

[0400] (3) If φ j (z=0)>φ ex And φ j (z=a)>φ ex or φ j (z=0)<φ st If so, the cutting edge does not cut;

[0401] Where, φ j (z=k) represents the instantaneous entry angle of the cutting edge j when the axial depth of cut is k, where k is a constant; z j,1 ,z j,2 These represent the lower and upper limits of axial cutting of cutting edge j, respectively;

[0402] Summing the cutting forces over all cutting edges yields the instantaneous total cutting force on the tool at the instantaneous entry angle φ. The cutting force acting on the milling cutter is:

[0403]

[0404] In the formula, F xj ,F yj ,F zj These represent the feed, normal, and axial cutting forces acting on the cutting edge j, respectively.

[0405] S302: Use the average value method to obtain cutting force parameters.

[0406] Since the cutting edge only cuts within the intrusion zone, the average milling force per tooth cycle is:

[0407]

[0408] In the formula, F is the average milling force per tooth cycle. q (φ) is the dynamic function of cutting force as a function of instantaneous angle of entry;

[0409] The derivation yields:

[0410]

[0411] In the formula, T is the tool rotation period, and F q (t) is the dynamic function of cutting force over time;

[0412] Therefore, in order to obtain the cutting force coefficient, a set of milling experiments were conducted using a milling cutter at different feed rates c, but the axial cutting depth a and radial cutting depth remained unchanged. The average milling force per tooth cycle was obtained, and then the average value method was used to evaluate the mechanism of the cutting force coefficient.

[0413] right Integrating, we get:

[0414]

[0415] In the formula, φ represents the average milling force per tooth cycle in the x, y, and z directions, respectively, and φ is the instantaneous angle of entry.

[0416] The entry angle and the removal angle for slotting milling are φ st =0 and φ ex =π, so the average force per tooth cycle simplifies to:

[0417]

[0418] The average cutting force is expressed as a linear function of the feed rate and the offset caused by the edge force:

[0419]

[0420] In the formula, The average shear force component per tooth cycle. This represents the average plowing force component per tooth cycle.

[0421] Calculate the average force at each feed rate, estimate the cutting force components through linear regression of the data, and obtain the cutting force coefficient:

[0422]

[0423] In the formula, These represent the average shear force components per tooth cycle in the x, y, and z directions, respectively. These represent the average plowing force components per tooth cycle in the x, y, and z directions, respectively.

[0424] S303: Derivation of cutting load models for each axis based on mechanical kinematics.

[0425] First, kinematic modeling is performed on the five-axis CNC machine tool used in this invention. During the operation of the machine tool, the absolute coordinate system O of the machine tool... M X M Y M Z M Let the B-axis coordinate system O remain unchanged. B X B Y B Z B The origin is fixed on the B-axis, and the origin of the coordinate system is O. B In machine tool coordinate system O M X M Y M Z M The coordinates in (L) mbx ,L mby ,L mbz When the B-axis rotates by α degrees, O B X B Y B Z B Compared to O M X M Y M Z M The pose transformation matrix is:

[0426]

[0427] Since the rotation centers of the B and C rotary tables intersect, a table coordinate system O is established with the B axis as the Y-axis and the positive direction of the machine tool coordinate system's Y-axis as the positive direction; and the C axis as the Z-axis, with the direction away from the table surface as the positive direction. G X G Y G Z G It is stipulated that regardless of the rotation of axes B and C, the Y and Z axes of this coordinate system will always be collinear with the axes of axes B and C, and the origin of the coordinate system is O. G In machine tool coordinate system O B X B Y B Z BThe coordinates in (L) bgx ,L bgy ,L bgz ), then O G X G Y G Z G Compared to O B X B Y B Z B The pose transformation matrix is:

[0428]

[0429] The workpiece is fixed to the worktable and rotates with the worktable around the C-axis. The workpiece coordinate system O W X W Y W Z W Fixed to the workpiece, with the origin O of the coordinate system W In the worktable coordinate system O G X G Y G Z G The coordinates in (L) gwx ,L gwy ,L gwz When the C-axis rotates by β degrees, O W X W YWZ W Compared to O G X G Y G Z G The pose transformation matrix is:

[0430]

[0431] With the tool tip as the origin, the tool feed direction as the positive X-axis, and the direction along the tool axis pointing to the spindle as the positive Z-axis, a right-handed coordinate system is established as the tool coordinate system. Based on the machining path and its machining parameters, the X, Y, and Z axial cutting forces of each cutting edge of the tool at that moment and position in the tool coordinate system can be obtained using methods S301 and S302. According to the translation theorem of forces, these three axial forces are translated to the tool tip. Considering the high rigidity of the tool and the small torque generated by the force translation, these three axial torques are ignored to reduce the amount of calculation. Thus, the X, Y, and Z axial cutting forces of the tool tip at that moment and position in the tool coordinate system can be obtained.

[0432] Given a system of parametric equations for the machining path with time as the independent variable, differentiating each parametric equation with respect to time yields the tangent vector at any given moment along the machining path. After normalization, this provides the unit vector of the tool coordinate system's X-axis relative to the workpiece coordinate system. In the machine coordinate system, the unit vector of the tool coordinate system's Z-axis is always (0, 0, 1). Performing the coordinate transformation described above on this unit vector, followed by normalization, yields the unit vector of the tool coordinate system's Z-axis relative to the workpiece coordinate system. Taking the outer product of these two unit vectors and normalizing the result, we obtain the unit vector of the tool coordinate system's Y-axis relative to the workpiece coordinate system. Multiplying the magnitude of the three-dimensional forces acting on the tool by this unit vector and translating it to the tool tip coordinates in the workpiece coordinate system, we obtain the three-dimensional cutting forces acting on the tool in the workpiece coordinate system.

[0433] For the load force (torque) model of each axis of a five-axis CNC machine tool, it mainly involves the load models of six axes: the spindle, the X, Y, and Z linear feed axes, and the B and C rotary feed axes. Among them, the spindle load torque model is relatively simple to build; it only requires multiplying the x and y components of the theoretical cutting force derived in methods S301 and S302 by the tool radius to obtain the spindle load torque. The load models of the other five axes can be divided into rotary feed axis load models and linear feed axis load models according to their mechanical structure.

[0434] For the rotary feed axis, the tool tip interpolation point in the machine coordinate system is transformed into the tool tip interpolation point in the table coordinate system using the homogeneous coordinate transformation matrix described above. Simultaneously, based on the inverse kinematics method using the Jacobian matrix, the three-dimensional resultant force vector of the tool in the workpiece coordinate system is transformed into the three-dimensional resultant force vector in the table coordinate system. This resultant force vector is then multiplied by the force arm vector formed by the tool tip interpolation point and the table origin in the table coordinate system. The resulting J-axis component is the load torque of the B-axis, and the K-axis component is the load torque of the C-axis. The solution process is as follows: Figure 17 As shown.

[0435] For linear feed axes, the inverse kinematics method based on the Jacobian matrix converts the three-dimensional resultant force vector of the tool in the workpiece coordinate system into a three-dimensional resultant force vector in the machine coordinate system. This vector is then translated to the tool tip interpolation point in the machine coordinate system to obtain the resultant force on the tool in the machine coordinate system. The x-component obtained from decomposing this resultant force is the X-axis load force, the y-component is the Y-axis load force, and the z-component is the Z-axis load force. The solution process is as follows: Figure 18 As shown.

[0436] Example 2: The five-axis CNC machine tool energy consumption modeling system based on energy flow footprint provided in this embodiment of the invention includes:

[0437] The parameter acquisition unit is configured to build a three-dimensional geometric model of a five-axis CNC machine tool and complete the assembly. It acquires component parameters based on SolidWorks quality assessment, Hertz contact theory and Huth semi-empirical formula. The component parameters include moment of inertia, polar moment of inertia, joint stiffness and fastener compliance coefficient.

[0438] The energy consumption model building unit is configured to build a five-axis CNC machine tool electromechanical coupling dynamics model based on power bond graph theory and SVPWM algorithm, and to build a machine tool machining process energy consumption model based on the electromechanical coupling dynamics model.

[0439] The load integration unit is configured to establish a theoretical cutting force model based on the Altintas mechanical model, obtain the cutting load of each axis according to the kinematic equation of the five-axis CNC machine tool, and apply the cutting load of each axis to the load end of the machine tool machining process energy consumption model to complete the machine tool machining process energy consumption model.

[0440] To further demonstrate the positive effects of the above embodiments, the present invention conducts the following experiments based on the above technical solutions.

[0441] For the identification of cutting force parameters in the Altintas mechanical model, 7050-T7451 aluminum alloy was selected as the machining material. To ensure cutting quality, a common single-flute helical end mill for aluminum alloys was selected as the cutting tool, with a main dimension of 6mm diameter and 45° helix angle. The feed rate gradient experiment was designed, and the cutting parameters and average milling force calculation results after tool coefficient calibration are shown in Table 6.

[0442] Table 6. Calculation results of cutting parameters and average milling force under tool coefficient calibration.

[0443]

[0444] Based on the calculation results of the average milling force, the milling force coefficient is obtained as: K tc =549.99 N / mm 2 K rc = 367.14 N / mm 2 K te =59.28 N / mm 2 K re = 45.94 N / mm 2 K ac = -101.65 N / mm 2 K ae =5.42 N / mm 2 This coefficient represents the milling force coefficient of the current milling cutter on 7050-T7451 aluminum alloy.

[0445] To verify the cutting force of the Altintas mechanical model, the experimental condition was the machining of a circular groove. The diameter of the groove was set to 46 mm, the width to 6 mm, and the spindle speed to 3000 r·min. -1 With an axial depth of 2mm and a tangential feed rate of 240mm / min, the experimental working conditions were simulated using the Altintas mechanical model to obtain simulation data. The actual three-dimensional cutting forces during the machining process were collected as real data. The three-dimensional cutting forces are as follows: Figure 20 , Figure 21 and Figure 22 As shown, based on the comparison of experimental and simulation results, it can be seen that some areas with larger errors may be affected by the linkage steering accuracy of the machine tool feed axis, resulting in more obvious vibration phenomena. At the same time, the vibration generated during the cutting process will also affect the instability of the three-phase force gauge measurement process. Considering the sources and effects of errors, the accuracy of this cutting force model basically meets the usage requirements.

[0446] This invention, using the X-axis linear feed system and spindle system during the machining of an S-shaped workpiece as examples, simulates the machining energy consumption in the MatLab-Simulink environment, obtaining the energy consumption curves as shown below. Figure 23 and Figure 24 As shown.

[0447] The above description is only a preferred embodiment of the present invention, but the scope of protection of the present invention is not limited thereto. Any modifications, equivalent substitutions and improvements made by those skilled in the art within the scope of the technology disclosed in the present invention and within the spirit and principles of the present invention should be covered within the scope of protection of the present invention.

Claims

1. A method for modeling energy consumption of a five-axis CNC machine tool based on energy flow footprint, characterized in that, The method includes the following steps: S1. Build a three-dimensional geometric model of a five-axis CNC machine tool and complete the assembly. Obtain the component parameters based on SolidWorks quality assessment, Hertz contact theory and Huth semi-empirical formula. The component parameters include moment of inertia, polar moment of inertia, joint stiffness and fastener compliance coefficient. S2. Based on power bond graph theory and SVPWM algorithm, a five-axis CNC machine tool electromechanical coupling dynamics model is built, and an energy consumption model of the machine tool machining process is constructed based on the electromechanical coupling dynamics model. S3. Based on the Altintas mechanical model, establish a theoretical cutting force model. Obtain the cutting load of each axis according to the kinematic equation of the five-axis CNC machine tool. Apply the cutting load of each axis to the load end of the machine tool machining process energy consumption model to complete the machine tool machining process energy consumption model.

2. The energy consumption modeling method for five-axis CNC machine tools based on energy flow footprint according to claim 1, characterized in that, In step S1, SolidWorks quality assessment is used to obtain the rotational inertia and polar moment of inertia of the component, including: If the machine tool transmission system is a cylindrical rotating body, then the moment of inertia about the rotation center axis is: In the formula, J is the moment of inertia, m is the mass of the cylinder, and R is the radius of the cylinder; For a torsion shaft, the angle through which the two cross sections rotate relative to each other during torsion is: In the formula, Let T be the torsion angle, T be the torque on the torsion shaft, l be the length of the torsion shaft, G be the shear modulus, and I be the torque on the torsion shaft. p The polar moment of inertia of the cross section; The torsional stiffness K is: In the formula, K is the torsional stiffness.

3. The energy consumption modeling method for five-axis CNC machine tools based on energy flow footprint according to claim 1, characterized in that, In step S1, the stiffness of the joint at key locations is obtained based on Hertz contact theory, using simplified assumptions: (1) The contacting object only produces elastic deformation and obeys Hooke's law; (2) The load is perpendicular to the contact surface, and the friction between the contacting objects is negligible; (3) The size of the contact surface is small compared to the radius of curvature of the contact object surface; When two contacting objects are under the action of load Q, the contact point will expand into a contact surface. The projection of this contact surface onto the plane perpendicular to the contact normal is an ellipse. Within the contact area, the contact stress is distributed in a semi-elliptical pattern. The contact stress and strain formulas derived from Hertz contact theory are as follows: In the formula, δ is the contact elastic approximation, E' is the equivalent elastic modulus, Q is the normal pressure, ∑δ is the combined curvature of the two contacting bodies at the contact point, and K(e),m a All are Hertz coefficients; The formula for calculating cosτ is obtained by looking up the auxiliary variable cosτ in a table: In the formula, ρ 11 ,ρ 12 Let ρ be the curvature of the two principal planes of the contact surface 1. 21 ,ρ 22 The curvatures of the two principal planes of the contact surface 2; Then we have: In the formula, R 11 ,R 12 Let R be the radii of curvature of the two principal planes of the contact surface 1. 21 ,R 22 The radii of curvature of the two principal planes of the contact surface 2; E' is obtained from the following formula: In the formula, E1 and E2 are the elastic moduli of the two contacting objects, and μ1 and μ2 are the Poisson's ratios of the two contacting objects. To simplify the formula, let: Substituting into the formula, we get: According to simplifying assumption 3, k h Approximately constant during the contact process, the contact stiffness K is obtained from the formula:

4. The energy consumption modeling method for five-axis CNC machine tools based on energy flow footprint according to claim 3, characterized in that, In step S1, the fastener compliance coefficient is obtained based on Huth's semi-empirical formula; For the fixed joint, the bolt connection stiffness between the harmonic reducer and the turntable is analyzed. The fastener compliance coefficient is calculated using Huth's semi-empirical formula. Then: In the formula, C is the fastener flexibility coefficient, t1 and t2 are the thickness of the connecting plate, E1 and E2 are the elastic modulus of the connected plate, E3 is the elastic modulus of the fastener, n is the fit coefficient, d is the fastener diameter, and a and b are empirical parameters.

5. The energy consumption modeling method for five-axis CNC machine tools based on energy flow footprint according to claim 1, characterized in that, In step S2, the electromechanical coupling dynamics model of the five-axis CNC machine tool includes power bond graph models of the spindle system, linear feed system, and rotary feed system. A state-space dynamics model is constructed using augmented bond graphs. The bond graph model of the linear feed system includes the mechanical properties of the coupling and bearing pairs, ball screws, ball screw-nut pairs, and the driven load. These mechanical properties include the mass coefficient, stiffness coefficient, damping coefficient, and converter module. The bond graph model of the rotary feed system includes the mechanical properties of the wave generator, rigid wheel, harmonic gear reducer, and rotating load. Vector control of a three-phase permanent magnet synchronous motor is implemented based on the SVPWM algorithm, using i d =0 control method; selecting a mathematical model in the synchronous rotating coordinate system dq, the three-phase voltage is transformed from the three-phase natural coordinate system to the synchronous rotating coordinate system using the Clark transformation matrix and the Park transformation matrix, then: In the formula, T 3s / 2s T is the Clark coordinate transformation matrix. 2s / 2r Let θ be the Park coordinate transformation matrix. e For the mechanical angle of the motor; The bond graph model of a three-phase permanent magnet synchronous motor in a synchronous rotating coordinate system includes the motor inductance resistance, electromagnetic torque coefficient, and motor rotor; by augmenting the bond graph, the augmented bond graph model of the motor is obtained.

6. The energy consumption modeling method for five-axis CNC machine tools based on energy flow footprint according to claim 5, characterized in that, The construction of the system dynamics model and the energy consumption model based on energy flow footprint includes: (1) Determine the input variables and energy state variables; (2) Construct the initial set of state equations; (3) Simplify the initial system of equations into state-space form; The construction of the state-space dynamic model of the linear feed system includes: Based on the augmented bond diagram of the linear feed system and the permanent magnet synchronous motor, the system has nine energy storage elements, namely: inertial elements I2, I6, I... 12 ,I 15 ,I 20 i 26 and capacitive element C 18 C 23 C 26 The component numbers are determined according to the bond graph augmentation rules, proceeding from left to right and from top to bottom; among them, inertial components I2, I6, and I... 12 ,I 20 ,I 26 With C 18 C 26 There exists an integral causal relationship, I 15 With C 23 Differential causality exists; The system's input variables are determined as follows: U=[U d (t) U q (t) F(t)] T In the formula, U d (t) represents the d-phase input voltage of the motor, U q (t) represents the input voltage of phase q of the motor, and F(t) represents the load force at the end of the workbench of the feed system. The system's state variables are: X=[p2 p6 p 12 q 18 q 20 q 26 p 28 ] T In the formula, p2, p6, p 12 ,p 20 ,p 28 For the generalized momentum of the corresponding numbered element, q 18 ,q 26 This refers to the generalized displacement of the corresponding numbered element; Based on the properties of inertial elements with integral causality, we have: In the formula, f2, f6, f 12 ,f 20 ,f 28 For the corresponding numbered elements, I2, I6, I 12 ,I 20 ,I 28 These are the inertial parameters of the corresponding numbered components; Based on the properties of capacitive elements with integral causality, we have: In the formula, e 18 ,e 26 For the potential variable of the corresponding numbered element, C 18 C 26 For the capacitive parameters of the corresponding numbered components; Based on the existence of differential causality in inertial element I 15 The properties yield: p 15 =f 15 ·I 15 In the formula, p 15 For the generalized momentum of the corresponding numbered element, f 15 For the flow variable corresponding to the numbered element, I 15 These are the inertial parameters of the corresponding numbered components; From the system bond graph: Differentiating the above equation, we get: In the formula, The potential variable of the corresponding numbered element; Based on the capacitive element C, which exhibits differential causality. 23 The properties yield: In the formula, For the flow variable corresponding to the numbered element, C 23 Here are the capacitive parameters for the corresponding numbered components, and TF is the converter module; Based on resistive components R3, R7, R 13 ,R 16 ,R 21 ,R 29 The properties yield: In the formula, e3, e7, e 13 ,e 16 ,e 21 ,e 29 For the potential variables of the corresponding numbered elements, R3, R7, R 13 ,R 16 ,R 21 ,R 29 For the resistive parameters of the corresponding numbered components, f3, f7, f 13 ,f 16 ,f 21 ,f 29 For the corresponding numbered element's flow variable; Based on the causal relationships between the components and the direction of energy flow, we can conclude that: In the formula, e4,e8,e9,e 10 ,e 11 ,e 14 ,e 19 ,e 21 ,e 22 ,e 27 ,e 30 For the potential variable of the corresponding numbered element, MGY d MGY q The module of the gyroscope. f 17 ,f 19 ,f 25 ,f 27 For the corresponding numbered element's flow variable; The simplified and integrated dynamic model of the linear feed system is as follows: Substituting the actual physical quantities into the system dynamics model, we get: In the formula, L d For d-phase inductance, i d Let R be the phase current (d). d For the d-phase resistance, P n ω is the number of pole pairs of the motor. n L is the mechanical angular velocity of the motor. q For q-phase inductance, i q Let R be the phase q current. q Let ψ be the resistance of phase q. f For permanent magnet flux linkage, J z J is the moment of inertia of the motor rotor. lzq R is the moment of inertia of the bearing and coupling. b R is the rotor damping coefficient of the motor. lzq τ is the damping coefficient of the bearing and coupling. lzq C represents the coupling torque. lzq For the stiffness of the coupling, ω s J is the screw speed. s R is the moment of inertia of the lead screw. s τ is the damping coefficient of the lead screw. l To drive the load force, C l C represents the stiffness of the ball screw joint. s For the stiffness of the lead screw, ω l To drive the load speed, J l To drive the load mass, R l The damping coefficient for the driven load; Let E be the sum of the quadratic forms on the left-hand side of each term in the system of differential equations, then: Differentiating E with respect to t, we get: Substitute the system of differential equations into: Rearranging the terms, we obtain the energy consumption equation P. x for: In the formula, P x E represents the input power of the linear feed system; E represents the energy stored in the system's energy storage element, with each term in the polynomial corresponding to the stored energy of the energy storage element; the remaining terms represent the power loss and processing power of the system's energy-consuming elements. This formula follows the law of conservation of energy and is used to characterize the direction and footprint of energy flow. Construct the energy consumption equations for the remaining axes, and sum them to obtain the machine tool energy consumption equation: P=P x +P y +P z +P b +P c In the formula, P is the machine tool processing power, P x P represents the machining power of the X-axis feed system. y P is the machining power of the y-axis feed system. z P is the machining power of the z-axis feed system. b P is the machining power of the b-axis feed system. c The machining power of the C-axis feed system.

7. The energy consumption modeling method for five-axis CNC machine tools based on energy flow footprint according to claim 1, characterized in that, In step S3, establishing the theoretical cutting force model based on the Altintas mechanical model includes: The formula for calculating the hysteresis angle ψ at an axial depth of cut z on a spiral end mill with a helix angle of β is: In the formula, D is the diameter of the milling cutter; Establish a right-handed coordinate system with the feed direction of the helical end mill as the positive X-axis and the direction pointing to the spindle as the positive Z-axis; measure the instantaneous entry angle clockwise from the y-axis. When the instantaneous entry angle of the bottom point of the end mill's reference cutting edge is φ, the instantaneous entry angle of the cutting edge point located at an axial distance of z is φ-ψ. A helical end mill with a helix angle of β, a diameter of D, and N cutting edges, has a tooth pitch angle of β. When the reference instantaneous angle of entry at the bottom of the cutting edge is φ, then the instantaneous angles of entry at the bottom of the other cutting edges are: φ j (0)=φ+jφ p ,j=0,1,2...(N-1) Let the lag angle coefficient Then at the axial depth of tangency z, the hysteresis angle ψ = k β Then, the instantaneous entry angle of the cutting edge j at the axial depth of cut z is: φ j (z)=φ+jφ p -k β With In the formula, φ j (z) is the instantaneous entry angle of the cutting edge j at the axial depth of cut z, φ is the reference instantaneous entry angle at the bottom of the cutting edge, j is the cutting edge count, and z is the axial depth of cut; The tangential, radial, and axial forces acting on the infinitesimal cutting edge at height dz are: In the formula, dF t,j (φ,z),dF r,j (φ,z),dF a,j (φ, z) represent the tangential force, radial force, and axial force acting on the infinitesimal cutting edge j, respectively. K tc ,K rc ,K ac The shear force coefficients h in the tangential, radial, and axial directions, respectively. j (φ j (z) represents the instantaneous cutting thickness, K te ,K re ,K ae denoted as the edge force coefficients in the tangential, radial, and axial directions, respectively, and dz is the height of the cutting micro-element; Instantaneous cutting thickness h j (φ j (z))=csinφ j (z) decomposes the force into three directions: feed, normal, and axial. In the formula, dF x,j (φ j (z)),dF y,j (φ j (z)),dF z,j (φ j (z) represents the feed, normal, and axial forces acting on the infinitesimal cutting edge j, respectively, and c is the milling cutter feed rate in millimeters per tooth angle; Substituting the expressions for instantaneous cutting thickness, tangential force, radial force, and axial force into the equation: Therefore, the total cutting force generated by the cutting edge is: In the formula, z j,1 (φ j (z)) and z j,2 (φ j (z) represents the lower and upper limits of the axial cutting edge, respectively; From φ j (z)=φ+jφ p -k β z dφ j (z)=-k β Substituting dz into the total cutting force equation, we get: In the formula, F x,j (φ j (z)),F y,j (φ j (z)),F z,j (φ j (z) represent the feed, normal, and axial cutting forces acting on the cutting edge j, respectively; Based on the similarity triangle theorem and geometric relationships, the depth of cut z and the tool angle φ at that depth of cut are obtained during the inclination cutting process. j The relation for (z) is: In the formula, a is the deepest axial cutting depth, and R is the tool radius. The angle of inclination for cutting at an angle; Substituting the tool angle formula for the axial depth of cut, we get: During climb milling, the tool approach angle is: In the formula, φ st Let a be the tool entry angle. e Radial depth of cut; During climb milling, the cutter exit angle is: f ex =π In the formula, φ ex Cut an angle with the tool; When the spindle speed is n, taking the starting point of the cycle as the bottom of a certain cutting edge on the tool entering the cutting area, and after a time interval t after the bottom of that cutting edge enters the tool's entry angle, the instantaneous entry angle of the bottom of that cutting edge at this moment is: Then the instantaneous angle of entry at the bottom of the remaining cutting edges is: Then the instantaneous angle of entry of the cutting edge j at the axial depth of cut z is: Each cutting segment interacts with the cutting edge in five different ways to define the axial integral boundary, specifically including: (1) If φ st < φ j (z = 0)< φ ex , then z j,1 = 0: If φ st < φ j (z = a)< φ ex , then z j,2 = a; If φ j (z = a)< φ st , then z j,2 =(φ + jφ p - φ st ) / k β =(nt + jφ p ) / k β ; (2) If φ j (z = 0) > φ ex and φ j (z = a) < φ ex , then z j,1 = (φ + jφ p - φ ex ) / k β = (nt + φ st + jφ p - φ ex ) / k β ): If φ j (z = a) < φ st , then z j,2 = (φ + jφ p - φ st ) / k β = (nt + jφ p ) / k β ; If φ j (z = a) > φ st , then z j,2 = a; (3) If φ j (z=0)>φ ex And φ j (z=a)>φ ex or φ j (z=0)<φ st If so, the cutting edge does not cut; Where, φ j (z=k) represents the instantaneous entry angle of the cutting edge j when the axial depth of cut is k, where k is a constant; z j,1 ,z j,2 These represent the lower and upper limits of axial cutting of cutting edge j, respectively; Summing the cutting forces over all cutting edges yields the instantaneous total cutting force on the tool at the instantaneous entry angle φ. The cutting force acting on the milling cutter is: In the formula, F xj ,F yj ,F zj These represent the feed, normal, and axial cutting forces acting on the cutting edge j, respectively.

8. The energy consumption modeling method for five-axis CNC machine tools based on energy flow footprint according to claim 7, characterized in that, The cutting force parameters are obtained using the average value method, specifically including: Since the cutting edge only cuts within the intrusion zone, the average milling force per tooth cycle is: In the formula, F is the average milling force per tooth cycle. q (φ) is the dynamic function of cutting force as a function of instantaneous angle of entry; The derivation yields: In the formula, T is the tool rotation period, and F q (t) is the dynamic function of cutting force over time; Therefore, in order to obtain the cutting force coefficient, a set of milling experiments were conducted using a milling cutter at different feed rates c, but the axial cutting depth a and radial cutting depth remained unchanged. The average milling force per tooth cycle was obtained, and then the average value method was used to evaluate the mechanism of the cutting force coefficient. right Integrating, we get: In the formula, φ represents the average milling force per tooth cycle in the x, y, and z directions, respectively, and φ is the instantaneous angle of entry. The entry angle and the removal angle for slotting milling are φ st =0 and φ ex =π, so the average force per tooth cycle simplifies to: The average cutting force is expressed as a linear function of the feed rate and the offset caused by the edge force: In the formula, The average shear force component per tooth cycle. This represents the average plowing force component per tooth cycle. Calculate the average force at each feed rate, estimate the cutting force components through linear regression of the data, and obtain the cutting force coefficient: In the formula, These represent the average shear force components per tooth cycle in the x, y, and z directions, respectively. These represent the average plowing force components per tooth cycle in the x, y, and z directions, respectively.

9. The energy consumption modeling method for five-axis CNC machine tools based on energy flow footprint according to claim 8, characterized in that, Based on the kinematic relationships of mechanics, the cutting load models for each axis are derived, including: Perform kinematic modeling on a five-axis CNC machine tool. During the machine tool's operation, the absolute coordinate system O of the machine tool... M X M Y M Z M Let the B-axis coordinate system O remain unchanged. B X B Y B Z B The origin is fixed on the B-axis, and the origin of the coordinate system is O. B In machine tool coordinate system O M X M Y M Z M The coordinates in (L) mbx ,L mby ,L mbz When the B-axis rotates by α degrees, O B X B Y B Z B Compared to O M X M Y M Z M The pose transformation matrix is: Since the rotation centers of the B and C rotary tables intersect, a table coordinate system O is established with the B axis as the Y-axis and the positive direction of the machine tool coordinate system's Y-axis as the positive direction; and the C axis as the Z-axis, with the direction away from the table surface as the positive direction. G X G Y G Z G It is stipulated that regardless of the rotation of axes B and C, the Y and Z axes of this coordinate system will always be collinear with the axes of axes B and C, and the origin of the coordinate system is O. G In machine tool coordinate system O B X B Y B Z B The coordinates in (L) bgx ,L bgy ,L bgz ), then O G X G Y G Z G Compared to O B X B Y B Z B The pose transformation matrix is: The workpiece is fixed to the worktable and rotates with the worktable around the C-axis. The workpiece coordinate system O W X W Y W Z W Fixed to the workpiece, with the origin O of the coordinate system W In the worktable coordinate system O G X G Y G Z G The coordinates in (L) gwx ,L gwy ,L gwz When the C-axis rotates by β degrees, O W X W YWZ W Compared to O G X G Y G Z G The pose transformation matrix is: With the tool tip as the origin, the tool feed direction as the positive X-axis, and the direction along the tool axis pointing to the spindle as the positive Z-axis, a right-handed coordinate system is established as the tool coordinate system. Based on the machining path and machining parameters, the X, Y, and Z axial cutting forces experienced by each cutting edge of the tool at that moment and position in the tool coordinate system are obtained. According to the translation theorem of forces, these three axial forces are translated to the tool tip. Ignoring the torque of these three axial forces, the X, Y, and Z axial cutting forces experienced by the tool tip at that moment and position in the tool coordinate system are obtained. Given a system of parametric equations for the machining path with time as the independent variable, the derivative of each parametric equation with respect to time is used to obtain the tangent vector at any time point of the machining path. After normalization, the unit vector of the tool coordinate system X-axis in the workpiece coordinate system is obtained. In the machine coordinate system, the unit vector of the tool coordinate system Z-axis is always (0, 0, 1). After coordinate transformation and normalization of this unit vector, the unit vector of the tool coordinate system Z-axis in the workpiece coordinate system is obtained. The outer product of the two unit vectors is then normalized to obtain the unit vector of the tool coordinate system Y-axis in the workpiece coordinate system. The magnitude of the three-dimensional forces acting on the tool is multiplied by the unit vector and translated to the tool tip coordinates in the workpiece coordinate system to obtain the three-dimensional cutting forces acting on the tool in the workpiece coordinate system. The load force model for each axis of a five-axis CNC machine tool involves six axis load models: the spindle, X, Y, and Z linear feed axes, and B and C rotary feed axes. The spindle load torque model is constructed by multiplying the theoretical cutting force x and y components by the tool radius. The load models for the other five axes are divided into rotary feed axis load models and linear feed axis load models according to their mechanical structure. For the rotary feed axis, the tool tip interpolation point in the machine coordinate system is transformed into the tool tip interpolation point in the table coordinate system using a homogeneous coordinate transformation matrix. At the same time, based on the inverse kinematics method of the Jacobian matrix, the three-dimensional resultant force vector of the tool in the workpiece coordinate system is transformed into the three-dimensional resultant force vector in the table coordinate system. The resultant force vector is then multiplied by the force arm vector formed by the tool tip interpolation point and the origin of the table in the table coordinate system. The resulting J-axis component is the load torque of the B-axis, and the K-axis component is the load torque of the C-axis. For the linear feed axis, the inverse kinematics method based on the Jacobian matrix converts the three-dimensional resultant force vector of the tool in the workpiece coordinate system into the three-dimensional resultant force vector in the machine coordinate system. This vector is then translated to the tool tip interpolation point in the machine coordinate system to obtain the resultant force of the tool in the machine coordinate system. The x-axis component obtained by decomposing this resultant force is the load force of the X-axis, the y-axis component is the load force of the Y-axis, and the z-axis component is the load force of the Z-axis.

10. A five-axis CNC machine tool energy consumption modeling system based on energy flow footprint, characterized in that, This system is used to implement the energy consumption modeling method for five-axis CNC machine tools based on energy flow footprint as described in any one of claims 1-9, including: The parameter acquisition unit is configured to build a three-dimensional geometric model of a five-axis CNC machine tool and complete the assembly. It acquires component parameters based on SolidWorks quality assessment, Hertz contact theory and Huth semi-empirical formula. The component parameters include moment of inertia, polar moment of inertia, joint stiffness and fastener compliance coefficient. The energy consumption model building unit is configured to build a five-axis CNC machine tool electromechanical coupling dynamics model based on power bond graph theory and SVPWM algorithm, and to build a machine tool machining process energy consumption model based on the electromechanical coupling dynamics model. The load integration unit is configured to establish a theoretical cutting force model based on the Altintas mechanical model, obtain the cutting load of each axis according to the kinematic equation of the five-axis CNC machine tool, and apply the cutting load of each axis to the load end of the machine tool machining process energy consumption model to complete the machine tool machining process energy consumption model.

Images