Five-axis machining center rapid finite element modeling method based on g-code instruction driving

Abstract

The application discloses a quick finite element modeling method of a five-axis machining center driven by G code instructions, and belongs to the technical field of machine tool simulation and structural dynamics. Step 1: independent meshing of parts and pre-computation of element matrices; step 2: establishment of a machine tool kinematics model and construction of a joint surface node index table; step 3: G code instruction analysis and node coordinate updating based on inverse kinematics solution; step 4: automatic pairing of joint surface nodes based on the K nearest neighbor algorithm; step 5: master-slave node freedom degree deduplication and preprocessing based on a spherical hinge connection pair; step 6: overall matrix assembly based on matrix reuse and rotation transformation; and step 7: multiple constraint application and system matrix condensation processing. The application can realize dynamic updating of a machine tool finite element model, reflect time-varying structural characteristics driven by G code instructions in a machining process, provide a more accurate mechanical model for dynamic simulation, and support forward design iteration and machining process optimization of a machine tool structure.

Description

Technical Field

[0001] This invention belongs to the fields of machine tool simulation technology and structural dynamics technology, and relates to a rapid finite element modeling method for five-axis machining centers based on G-code instruction drive. It is applicable to the analysis of time-varying structural dynamic characteristics caused by component movement and rotation during machining, as well as the dynamic performance evaluation of the entire domain during the forward design stage of machine tool structure. Background Technology

[0002] Five-axis machining centers are core equipment in modern precision manufacturing, and their dynamic characteristics directly affect the machining accuracy and surface quality of workpieces. During machining, as linear motion components and angle adjustment components move and rotate along preset trajectories, the overall geometry of the machine tool continuously changes, leading to changes in the overall stiffness and mass distribution, which in turn causes significant time-varying structural characteristics in the natural frequencies and mode shapes. Especially under five-axis linkage conditions, changes in tool posture and machining position cause significant differences in the dynamic characteristics of the machine tool between different workstations, inducing cutting chatter under specific positioning conditions, threatening the surface integrity of the workpiece and accelerating tool wear.

[0003] In the forward design phase of machine tools, traditional finite element analysis methods typically only perform static modeling and analysis when the machine tool is in a fixed posture, failing to efficiently evaluate the dynamic performance changes at different machining stations across the entire travel range. Designers often need to manually adjust the model configuration and repeatedly submit calculation tasks, resulting in long design iteration cycles and high computational costs, making it difficult to meet the demands of modern machine tool products for rapid development and structural optimization requiring comprehensive dynamic characteristic evaluation. In the actual machining phase, when dynamic instability phenomena such as chatter occur, on-site operators typically adjust process parameters by stopping the machine for inspection, offline simulation, or repeated trial cuts. This approach not only significantly reduces production efficiency but also makes it difficult to accurately locate problematic workstations and lacks predictive and control methods tailored to specific machining tasks.

[0004] The existing technology has the following shortcomings: (1) Slow finite element model update: After the movement position of the machining center parts changes, it is necessary to manually remodel or adjust the connection relationship. The model is difficult to reuse and cannot keep in line with the actual running state of the machine tool driven by G code instructions. It is difficult to effectively obtain the dynamic characteristic changes in the machining process. (2) Low computational efficiency: Each time the working configuration of the machining center changes, the model structure needs to be readjusted, the mesh divided, and the element stiffness matrix and mass matrix recalculated. This results in a long modeling cycle and high computational cost, making it difficult to meet the needs of rapid simulation and iterative optimization. (3) Complex joint surface treatment: The joint surface of the slider guide rail and the mechanical connection joint surface of the machining center are usually simulated by spring-damping elements. However, the node pairing needs to be manually selected and cannot be automatically adjusted with the movement of the parts, resulting in a cumbersome and inefficient modeling process.

[0005] To address the above issues, scholars both domestically and internationally have conducted extensive research and proposed several solutions. For example, Chinese invention patent application number 201910908032.4 provides a multi-pose finite element modeling method for a five-axis moving beam gantry milling machine. This method pre-plans hard points at joints, uses offline manually input pose files to drive the local coordinate system transformation of various machine tool components, and then uses a spring-damped model to connect the hard points to generate a finite element model corresponding to the pose. Chinese invention patent 202610162135.0 provides an automatic optimization method and system for machining parameters of five-axis CNC machine tools. It uses discrete analytical techniques to convert G-code instructions into tool position sequence, and performs spatial interpolation by calling pre-calibrated and stored global static stiffness field data of the machine tool to calculate the instantaneous stiffness of the current pose and predict regenerative chatter. While existing technical solutions can model and simulate the dynamic characteristics of machine tools in specific positions, they rely on offline manual setting of position parameters or calling preset static databases for calculation. They cannot complete the automatic pairing of mating surface nodes and the rapid reconstruction of the underlying matrix of the whole machine, which are directly driven by G-code instructions.

[0006] Therefore, there is an urgent need for a finite element analysis method driven by G-code instructions and capable of automatic finite element modeling. This method should be able to automatically identify and couple interface nodes, quickly reuse pre-calculated matrices to update the overall machine model, and provide mechanical model support for efficient evaluation of the dynamic performance and structural optimization of machine tools across their entire travel range. Summary of the Invention

[0007] To address the problems existing in the prior art, this invention provides a rapid finite element modeling method for five-axis machining centers based on G-code command-driven approaches. This invention integrates inverse kinematics driven by G-code commands, reuse of component stiffness and mass matrices, and automatic node pairing processing at interface surfaces based on the K-nearest neighbor algorithm. This enables dynamic updating of the machine tool finite element model, obtaining the overall stiffness and mass matrices of the machine tool configuration under the current G-code command-driven approach. This reflects the time-varying structural characteristics driven by G-code commands during machining, providing a more accurate mechanical model for dynamic simulation and supporting the forward design iteration of the machine tool structure and optimization of the machining process.

[0008] To achieve the above objectives, the technical solution adopted by the present invention is as follows: A rapid finite element modeling method for five-axis machining centers based on G-code instructions includes the following steps: Step 1: Component-independent mesh generation and element matrix pre-calculation.

[0009] Finite element mesh discretization was used to discretize the components of the five-axis machining center, establishing finite element mesh models of each component under the initial configuration. Numerical integration was used to calculate the element stiffness matrix and element mass matrix of each element under the initial configuration, and an element-matrix index table was established. Specifically: Step 1.1: Establish geometric models for each component of the five-axis machining center. Constrain each component of the geometric model to the zero point position, i.e., the initial configuration of the machine tool. The initial configuration of the machine tool is when each motion axis is at zero stroke. The components include the bed, slide, column, slider, guide rail, and swivel head or rotary table. Using the mesh generation method, the geometric model is meshed with three-dimensional solid elements including hexahedral elements and tetrahedral elements to establish the finite element mesh model of each component. The model file containing node coordinates and element node numbers is obtained, and the material properties of each element are defined.

[0010] Step 1.2: Based on the finite element mesh model established in Step 1.1, the element stiffness matrix and element mass matrix of each element under the initial configuration are calculated using the Gaussian numerical integration method. The calculation formulas are as follows: (1) (2) In the formula, Indicates the unit number; Representation unit The element stiffness matrix under the initial configuration; Representation unit The unit mass matrix under the initial configuration; Representation unit The strain-displacement matrix; Representation unit The elasticity matrix; Representation unit The shape function matrix; Representation unit Material density; Representation unit The volume integral region; superscript This represents the matrix transpose operation.

[0011] Step 1.3 converts the element stiffness matrix and element mass matrix of each element obtained in Step 1.2 under the initial configuration into a sparse format. This is performed once during the initialization phase to establish and generate an element-matrix index table, which can be directly called when updating the model later.

[0012] Step 2: Machine Tool Kinematic Model Establishment and Mating Surface Node Index Table Construction. This step includes three parts, performed concurrently with Step 1: The first part is establishing a machine tool kinematic model describing the direction of the entire machine's series shaft system; the second part is automatically identifying boundary contact areas based on geometric features and constructing a mating surface node index table for each component; the third part is establishing a list of ball joint connections describing the internal fixed connection relationships of the machine tool, as well as an external fixed constraint table describing the external fixed constraints of the machine tool. Specifically: Step 2.1: Based on the axis sequence and transmission structure of the five-axis machining center, each axis of the machine tool is divided into two independent motion chains starting from the machine tool bed: the first independent motion chain from the machine tool bed to the workpiece coordinate system, and the second independent motion chain from the machine tool bed to the tool system.

[0013] A global machine tool coordinate system is established with the bed as the global reference datum. For rotating parts, a local coordinate system is established at their rotation center. For slide parts, a local coordinate system is established at the guide rail reference point under their initial configuration. For columns, a local coordinate system is established at the reference reference point where their bottom surface connects to the bed. A workpiece coordinate system describing the workpiece position and orientation is established at the workpiece design reference point.

[0014] Construct a library of third-order basic rotation matrices for the A, B, and C axes of a machine tool. The general formulas for the third-order fundamental rotation matrices of each axis are as follows: The third-order fundamental rotation matrix of the A-axis about the X-axis The formula is: (3) The third-order fundamental rotation matrix of the B-axis about the Y-axis The formula is: (4) The third-order fundamental rotation matrix of the C-axis about the Z-axis The formula is: (5) In the formula, Indicates the rotation angle of axis A; Indicates the rotation angle of the B-axis; This indicates the rotation angle of the C-axis.

[0015] Based on the actual axis series topology of the five-axis machining center, the first independent motion chain in the first independent motion chain... The third-order fundamental rotation matrix corresponding to each rotation axis and the second independent kinetic chain The third-order fundamental rotation matrix corresponding to each rotation axis The following mapping rules are used to obtain the third-order basic rotation matrix from the library. Select from: (6) (7) In the formula, This represents the selected fundamental rotation matrix. ; Then, the third-order basic rotation matrices of each rotating component in the first independent motion chain are multiplied and cascaded in order from the bed to the workpiece coordinate system to establish the transformation matrix of the first rotating component. The third-order basic rotation matrices of each rotating component in the second independent kinematic chain are multiplied and cascaded in order from the bed to the tool system to establish the transformation matrix of the second rotating component. The formula for its calculation is: (8) (9) In the formula, This represents the transformation matrix of the first rotating component; This indicates the total number of rotational axes contained in the first independent kinematic chain; Indicates the cascade index number of the rotating axis in the first independent kinematic chain; In the first independent kinematic chain, the first The third-order basic rotation matrix corresponding to each rotation axis; This represents the transformation matrix of the second rotating component; This indicates the total number of rotational axes contained in the second independent kinematic chain; Indicates the cascade index number of the rotation axis in the second independent kinematic chain; Indicating the second independent kinetic chain, the first The third-order basic rotation matrix corresponding to each rotation axis.

[0016] Based on the first independent kinematic chain, the second independent kinematic chain, the global machine tool coordinate system, the component local coordinate system, and the workpiece coordinate system, a machine tool kinematic model is established for inverse kinematics calculation under the tool tip following control function. This machine tool kinematic model describes the geometric mapping relationship between the tool tip target position vector in the workpiece coordinate system, the origin translation vector of the workpiece coordinate system in the global machine tool coordinate system, the position vector of the rotation center of the rotating component relative to the global machine tool coordinate system in the current configuration, the first rotating component transformation matrix, the second rotating component transformation matrix, and preset structural constants. The current configuration of the machine tool is defined as the geometric pose state of each component after spatial translation and rotation transformation relative to the initial configuration following the execution of the current G-code instruction.

[0017] Specifically, the target position vector of the tool tip in the workpiece coordinate system is transformed by the first rotating component transformation matrix, and then superimposed with the origin translation vector of the workpiece coordinate system in the global machine tool coordinate system to obtain the expression of the target position vector of the tool tip in the global machine tool coordinate system. The preset structural constant is then transformed by the second rotating component transformation matrix to obtain the geometric vector pointing from the rotation center of the rotating component to the tool tip in the current configuration. In machining mode with the tool tip following control function enabled, the expression of the target position vector of the tool tip in the global machine tool coordinate system is consistent with the actual tool tip position of the tool system. The actual tool tip position of the tool system is determined by the position vector of the rotation center of the rotating component relative to the global machine tool coordinate system in the current configuration and the preset structural constant after transformation by the second rotating component transformation matrix. Their geometric relationship is as follows: (10) By organizing the above geometric relationships, we obtain the position vector of the rotation center of the rotating component relative to the global machine tool coordinate system under the current configuration: (11) In the formula, This represents the position vector of the rotation center of the rotating component relative to the global machine tool coordinate system under the current configuration; This represents the transformation matrix of the first rotating component; This represents the target position vector of the tool tip in the workpiece coordinate system; This represents the translation vector of the origin of the workpiece coordinate system in the global machine tool coordinate system; This represents the transformation matrix of the second rotating component; This represents the preset structural constant, which is the initial geometric vector from the rotation center of the rotating component to the tool tip under the initial configuration, and is related to the length of the tool and the tool holder.

[0018] Step 2.2: Perform global sorting and independent numbering on the mesh nodes in the finite element mesh model established in Step 1.1 to obtain a list containing the global numbers of all nodes. Establish and generate a component-node index table, and mark the feature node areas of the slider bottom surface and the guide rail contact surface that participate in the contact of the mating surfaces in the component-node index table.

[0019] Step 2.3: For the slider bottom surface node area and guide rail contact surface feature node area marked in Step 2.2, based on the planar or cylindrical surface geometric features of the component and user-specified feature nodes, automatic identification processing is performed using the spatial distance tolerance method to obtain a subset of planar contact nodes and a subset of cylindrical surface contact nodes. The planar contact node subset and the cylindrical surface contact node subset are then merged and deduplicated to establish and generate a mating surface node index table. The user-specified feature nodes include vertex nodes used to define the spatial boundary of the planar contact area of ​​the component, and circumferential nodes used to fit the geometric features of the cylindrical contact surface of the component. The mating surface node index is obtained using feature nodes, specifically as follows: Planar node extraction: Based on three or more non-collinear specified feature nodes, fit a three-dimensional spatial plane equation, traverse all nodes on the corresponding component one by one, calculate the spatial perpendicular distance from the node to the plane, and classify the nodes whose spatial perpendicular distance is not greater than a preset surface tolerance into the planar contact node subset. The value range of the preset surface tolerance is taken as a fraction of the average edge length of the finite element mesh in the contact area of ​​the component. to .

[0020] Cylindrical surface node extraction: Based on the cylinder axis direction, a reference point on the axis, and the cylinder radius, the radial distance from the node to the cylinder axis is calculated. Nodes whose absolute value of the difference between the radial distance and the cylinder radius is not greater than a preset circular tolerance are included in the cylindrical surface contact node subset. The value range of the preset circular tolerance is also taken as a fraction of the average edge length of the finite element mesh in the component contact area. to The direction of the cylinder axis, the reference point on the axis, and the cylinder radius are determined by the cylindrical surface geometric parameters in the geometric model.

[0021] Merging and deduplication: Merge the subset of planar contact nodes and the subset of cylindrical contact nodes in different mating areas on the same component surface, remove duplicate node numbers, and determine the mating surface node index table.

[0022] Step 2.4: Obtain ball joint connection pair information describing the fixed connection relationship between various components inside the machine tool, and establish a ball joint connection pair list.

[0023] Step 2.5: Obtain the fixed boundary conditions connecting the bottom surface of the machine tool bed to the external foundation, the set of degrees of freedom subject to external fixed constraints, and establish the external fixed constraint table.

[0024] Step 3: G-code instruction parsing and node coordinate update based on inverse kinematics. The G-code instructions are read using the G-code instruction parsing method, and the machining trajectory contained in the G-code instructions is converted into the relative translational displacement and rotational attitude of each machine tool component. The current configuration of the machine tool is defined as the geometric pose of each component after spatial translation and rotation transformation relative to the initial configuration following the execution of the current G-code instruction. The spatial coordinates of the mesh nodes of each component are updated using the rigid body spatial coordinate transformation formula. Specifically: Step 3.1: Parse the G-code instructions line by line. The parsing method involves scanning the function words representing the machine tool's motion mode, including linear interpolation or rapid positioning. Based on the identified function words, extract the coordinate address characters representing the spatial geometric position and rotation axis attitude in the current line, including X, Y, Z, A, B, C, and their suffixes. Convert the extracted characters into real-valued numerical parameters to obtain the target position vector of the tool tip in the workpiece coordinate system. and the target rotation angle values ​​of each rotation axis , , .

[0025] Step 3.2, based on the tool tip tracking control principle, the target position vector of the tool tip in the workpiece coordinate system calculated in Step 3.1 is... The target rotation angle values ​​of each rotation axis are input into the machine tool kinematic model established in step 2.1 for inverse kinematic calculation, and the indexed values ​​are determined. Spatial translation vector of machine tool components and rotation matrix .

[0026] If the machine tool G-code instruction is in the machining state with the tool tip point following control function enabled, the machine tool kinematic model established in step 2.1 based on the global machine tool coordinate system, workpiece coordinate system, and component local coordinate system is invoked. Its inverse kinematics calculation formula is as follows: (12) In the formula, This represents the position vector of the rotation center of the rotating component relative to the global machine tool coordinate system under the current configuration; This represents the transformation matrix of the first rotating component; This represents the target position vector of the tool tip in the workpiece coordinate system; This represents the translation vector of the origin of the workpiece coordinate system in the global machine tool coordinate system; This represents the transformation matrix of the second rotating component; This represents the preset structural constants.

[0027] Using the calculated position vector The spatial translation vector of the rotating component is obtained by calculating the difference between the initial position vector of the rotating component's rotation center and the initial position vector of the rotating component under the initial configuration. : (13) In the formula, Indicates having an index number The spatial translation vector of the machine tool component in the global machine tool coordinate system; This indicates the index number of the machine tool component. For purely translational components, the stroke displacement of each axis relative to its initial configuration, obtained from the inverse kinematics operation of the current G-code instruction, is used as the spatial translation vector of the corresponding purely translational component.

[0028] Step 3.3, analyze the target rotation angle values ​​of each rotation axis obtained in Step 3.1. , , Substitute the values ​​into the third-order basic rotation matrix library constructed in step 2.1 respectively. In the process, the basic rotation matrix of each rotation axis under the current configuration is calculated; then, according to the series sequence of the machine tool's shaft system, the third-order basic rotation matrices of each rotation axis under the current configuration are multiplied by matrix multiplication to calculate the indexed matrix. Rotation matrix corresponding to machine tool components .

[0029] Step 3.4: Update the node coordinates in the finite element mesh models of each component from Step 1.1 using the following rigid body space coordinate transformation formula: For a purely translational component, the nodal coordinate update formula is as follows: (14) For rotating parts: (15) In the formula, Indicates having an index number Finite element node numbering inside machine tool components; Indicates having an index number The internal parts of the machine tool The global machine tool coordinate vector of each finite element node in the current configuration; Indicates having an index number The internal parts of the machine tool The global machine tool coordinate vector of each finite element node under the initial configuration; Indicates having an index number The spatial translation vector of the machine tool component in the global machine tool coordinate system; Indicates having an index number The rotation matrix corresponding to the machine tool component; This represents the initial position vector of the rotation center of the rotating component in the global machine tool coordinate system under the initial configuration.

[0030] Using the rigid body space coordinate transformation formula, the global machine tool coordinate vector of all mesh nodes of each component in the whole machine finite element model under the current configuration is calculated and obtained. Each component includes sliders and guide rails. The current configuration of the machine tool is determined by the updated global machine tool coordinate vector of the mesh nodes under the current configuration.

[0031] Step 4: Automatic node pairing at the mating surface based on the K-nearest neighbor algorithm. Using the spatial search method and the current machine tool configuration calculated in Step 3.4, matching nodes are automatically found among the nodes of the components that need to be connected. The contact state after the movement of each axis is determined, and a node pairing table, Pairs, is established to describe the connection between the slider and guide rail mating surfaces. Specifically: Step 4.1: For the slider-guide rail pair, extract the set of slider nodes involved in the contact based on the mating surface node index table established in Step 2.3. With guide rail node set .

[0032] Step 4.2, denote the global machine tool coordinate vector of the corresponding slider node obtained in step 3.4 as... The global machine tool coordinate vector of the corresponding guide rail node is denoted as ; Step 4.3: Assign the global machine coordinate vectors of the slider node and guide rail node under the current configuration. and As the retrieval benchmark, the K-nearest neighbor method with a search count of 1 is used for spatial search and distance calculation to obtain the slider nodes. With guide rail node Spatial Euclidean distance between : (16) In the formula, Represents slider node With guide rail node The spatial Euclidean distance between them; Represents slider node The global machine tool coordinate vector under the current configuration; Indicates guide rail node Global machine tool coordinate vector in the current configuration; superscript This represents the matrix transpose operation.

[0033] The calculated spatial Euclidean distance The comparison is performed with a preset interface threshold; if the spatial Euclidean distance is... If the distance between the nodes is not greater than a preset joint surface threshold, then the two nodes are determined to form a joint surface connection, and the corresponding pair is added to the node pairing table Pairs. The preset joint surface threshold is set based on the nominal assembly gap between the slider and the guide rail and the average size of the finite element mesh of the contact surface, and its value range is [range missing]. to If the same guide rail node is matched repeatedly by multiple slider nodes, only a unique correspondence will be retained according to the principle of minimizing spatial Euclidean distance.

[0034] Step 4.4: Collect all corresponding relationships that satisfy the connection of the mating surfaces and establish the node pairing table Pairs.

[0035] Step 5: Deduplication and preprocessing of master-slave node degrees of freedom based on ball joint connections. This step is an auxiliary process performed in parallel with steps 3 and 4 in the overall calculation. Based on deformation compatibility conditions, the constraints of the ball joint connections between internal machine tool components are processed to fix the two internal components and maintain their common motion relationship, establishing a transformation relationship to eliminate the slave node degrees of freedom. Specifically: Step 5.1: Read the list of ball joint connections established in Step 2.4; perform uniqueness verification and deduplication on the master node degrees of freedom and slave node degrees of freedom corresponding to the ball joint connection pairs to ensure that each slave node degree of freedom in the whole machine finite element model corresponds to only one master node degree of freedom; divide the total degree of freedom of the whole machine finite element model into the set of eliminated slave node degrees of freedom and the set of retained degrees of freedom composed of master node degrees of freedom and unconstrained degrees of freedom; where master node degrees of freedom and slave node degrees of freedom respectively represent the finite element node degrees of freedom that play the main control role and the following role in the constraint equation formed by the ball joint connection.

[0036] Step 5.2: For the deduplicated linear constraint equations of the ball joint connection, each slave node degree of freedom is written as an algebraic linear function of the corresponding master node degree of freedom, thus establishing a multi-point constraint equation system whose block algebraic form satisfies: (17) Based on the algebraic coefficients of the multi-point constraint equation system, a master-slave constraint transformation matrix for the reduction of the degree of freedom of the ball joint connection constraint is established by expanding and combining the equations. Its explicit parsing expression is: (18) In the formula, This represents the constraint coefficient matrix corresponding to the degrees of freedom of each node; This represents the vector of slave node degrees of freedom, which consists of the slave node degrees of freedom. This represents the constraint coefficient matrix corresponding to the retained degrees of freedom; This represents the vector of retained degrees of freedom, consisting of the master node degrees of freedom and the unconstrained degrees of freedom. Represents the master-slave constraint transformation matrix; This represents the identity matrix of the same order as the number of degrees of freedom retained. The number of rows in this matrix equals the original total number of degrees of freedom of the system, and the number of columns equals the number of degrees of freedom retained after eliminating the degrees of freedom of slave nodes. It is used in subsequent steps to shrink the overall matrix to an algebraic space containing only the degrees of freedom of master nodes and unconstrained degrees of freedom.

[0037] Step 6: Overall matrix assembly based on matrix reuse and rotation transformation. Utilizing the element-matrix index table generated in Step 1.3 and the rotation matrix calculated in Step 3.3... Based on the energy equivalence principle of coordinate transformation, the element matrix is ​​reused, assembled, and spatially rotated to form the unconstrained overall stiffness matrix and the unconstrained overall mass matrix of the entire machine. Then, based on the matrix order of the unconstrained overall stiffness matrix and the unconstrained overall mass matrix, the unconstrained overall damping matrix of the entire machine is established. Specifically: Step 6.1, targeting the moving parts inside the machine tool Obtain its rotation matrix from step 3.3. If it has an index number If the machine tool component is a rotating component, then its third-order rotation matrix is ​​used. As basic unit items, they are laid out flat along the main diagonal without overlap. This establishes a size scale of block diagonal rotation transformation matrix ,in Indicates having an index number The total number of nodes in the machine tool components; if it has an index number. If a machine tool component is a purely translational component, then its rotation matrix... It is a third-order identity matrix, and its corresponding block diagonal rotation transformation matrix. for An identity matrix of order 1.

[0038] Step 6.2: Extract components from the unit-matrix index table established in Step 1.3. The element stiffness matrix and element mass matrix under the initial configuration are assembled into the initial component stiffness matrix. and the initial component mass matrix For rotating components, a spatial rotational transformation correction is applied using the derivation process based on the energy equivalence principle of continuum mechanics, as described in steps 6.2.1 to 6.2.2, and the global stiffness matrix of the current configuration after the rotational transformation is calculated. With the current configuration global mass matrix Specifically: Step 6.2.1: Based on the equivalence principle of system strain energy under coordinate transformation, when a machine tool component undergoes rigid body rotation, its inherent elastic characteristics remain unchanged relative to the initial configuration. Introduce the global displacement vector under the current configuration. and the nodal displacement vectors in the component's local coordinate system It satisfies the following strain energy invariance formula: (19) Based on the block diagonal rotation transformation matrix The orthogonal geometric characteristics, since the global displacement vector under the current configuration is the product of the block diagonal rotation transformation matrix and the initial displacement vector, satisfy the relation... Then the initial displacement vector can be expressed as Substituting the aforementioned correspondence into the strain energy invariance formula, the transformation formula is obtained by expansion: (20) Based on this, the rotational transformation formula for the stiffness matrix is ​​derived as follows: (twenty one) Step 6.2.2: Based on the equivalence principle of system kinetic energy under coordinate transformation, when a component undergoes rigid body rotation, its inherent inertial characteristics remain unchanged relative to the initial configuration. To ensure consistency in the expression of system kinetic energy, a global velocity vector under the current configuration is introduced. and the nodal velocity vectors in the component's local coordinate system It satisfies the following kinetic energy invariance formula: (twenty two) Since the initial velocity vector is equivalent to the product of the transpose of the block diagonal rotation transformation matrix and the global velocity vector, it satisfies the following relationship: Substituting the aforementioned correspondence into the kinetic energy invariance formula, the rotation transformation formula for the mass matrix is ​​derived as follows: (twenty three) In the formula, Indicates having an index number The initial component stiffness matrix of the machine tool component under the initial configuration; Indicates having an index number The initial component mass matrix of the machine tool components under the initial configuration; Indicates having an index number The dimensions of the machine tool components are The block diagonal rotation transformation matrix; Indicates having an index number The global stiffness matrix of the current configuration of the machine tool components after rotational transformation; Indicates having an index number The global mass matrix of the current configuration of the machine tool components after rotational transformation. For purely translational components, due to their block diagonal rotation transformation matrix... Since it is an identity matrix, no spatial rotation transformation correction is needed. The initial matrix, i.e., the global stiffness matrix of the current configuration, can be directly reused according to the above formula. Equal to the initial component stiffness matrix The current configuration global quality matrix equal to the initial component mass matrix .

[0039] Step 6.3: Obtain the global number of each mesh node using the component-node index table established in Step 2.2; based on the matrix assembly principle of the finite element direct stiffness method, perform global addressing and accumulation processing of the stiffness matrix and mass matrix: For the overall unconstrained stiffness matrix of the whole machine Based on the global number of the mesh nodes, the global stiffness matrix of each component's current configuration is determined. The matrix elements in the matrix are accumulated into the overall unconstrained stiffness matrix of the entire machine. In the process, the unconstrained overall stiffness matrix of the entire machine is formed through assembly. Its assembly mathematical formula is: (twenty four) For the unconstrained overall mass matrix of the whole machine Based on the global number of the grid nodes, the global mass matrix of each component's current configuration is calculated. The matrix elements in the matrix are accumulated into the unconstrained overall mass matrix of the whole machine. In the process, the unconstrained overall mass matrix of the entire machine is formed through assembly. Its assembly mathematical formula is: (25) In the formula, Represents the overall unconstrained stiffness matrix of the entire machine. The Middle line, number Matrix elements at column positions; ; represents the unconstrained overall mass matrix of the entire machine. The Middle line, number Matrix elements at column positions; Indicates the total number of independent moving parts contained in the machine tool; indicates the index number of the machine tool parts; Indicates having an index number The machine tool component in the global stiffness matrix under the current configuration is the first... line, number Matrix elements at column positions; Indicates having an index number The machine tool component in the global mass matrix under the current configuration is the [number]th [unit]. line, number Matrix elements at column positions; and These represent the global row number and global column number of the overall matrix, respectively; and These represent the internal row number and internal column number of the current configuration global stiffness matrix or the current configuration global mass matrix of the corresponding component, respectively; Since the finite element mesh model uses three-dimensional solid elements with 3 translational degrees of freedom per node, if the global number of a certain mesh node is... Then the global row number corresponding to the grid node in the overall matrix or global column number The set of values ​​is This enables the mapping from the global number of grid nodes to the row and column positions of the overall matrix.

[0040] Step 6.4: Based on the matrix orders of the unconstrained overall stiffness matrix and the unconstrained overall mass matrix of the whole machine assembled in Step 6.3, establish an unconstrained overall damping matrix of the same order, and initialize the unconstrained overall damping matrix of the whole machine to a zero matrix; the unconstrained overall damping matrix of the whole machine is used to receive the accumulated terms of the damping matrix of the damping unit corresponding to the joint surface connection in subsequent steps, and its initialization expression is: (26) In the formula, This represents the unconstrained overall damping matrix of the entire machine after initialization. This represents the matrix order of the unconstrained overall stiffness matrix and the unconstrained overall mass matrix of the entire machine; express A zero matrix of order 1.

[0041] Step 7: Applying multiple constraints and system matrix shrunk processing. This involves combining the node pairing table (Pairs) established in Step 4.4 with the master-slave constraint transformation matrix established in Step 5.2. The unconstrained overall stiffness matrix of the whole machine assembled in step 6.3 Unconstrained overall mass matrix of the whole machine And perform multi-constraint processing on the unconstrained overall damping matrix of the whole machine established in step 6.4; through the unconstrained overall stiffness matrix of the whole machine With respect to the unconstrained overall damping matrix of the whole machine A flexible connection is applied to the joint surface, and the degrees of freedom from the nodes and those subject to external fixed constraints are eliminated sequentially to obtain the final overall stiffness matrix after eliminating the influence of rigid body displacements in the system. Final overall quality matrix and the final overall damping matrix Specifically: Step 7.1 involves connecting the mating surfaces of the nodes in the Pairs table established in Step 4.4, i.e., the corresponding slider nodes. With guide rail node By introducing spatial spring elements that simulate contact stiffness characteristics and spatial damping elements that simulate energy dissipation characteristics, the corresponding spring element stiffness matrix is ​​established. With the damping matrix of the damping unit The damping matrix of the damping unit is used to accumulate in the overall unconstrained damping matrix of the whole machine established in step 6.4. (27) (28) In the formula, This represents the stiffness matrix of the spring element established by the translational stiffness coefficient of the mating surface; , , , These are all stiffness sub-matrices used to describe the stiffness characteristics of the physical connection between paired nodes; specifically, equivalent translational stiffness coefficients are introduced from the mating surface in the three spatial translational directions X, Y, and Z of the global machine tool coordinate system. , , This corresponds to the Pairs slider node in the node pairing table. The third-order self-stiffness submatrix With the corresponding guide rail node The third-order self-stiffness submatrix and slider nodes With guide rail node The coupling stiffness submatrix between and Its matrix representation is as follows: (29) (30) In the formula, This represents the damping matrix of the damping element established by the translational damping coefficient of the interface; , , , These are all damping sub-matrices used to describe the energy dissipation damping characteristics between paired nodes; specifically, equivalent translational damping coefficients of the joint surface in the three spatial translational directions X, Y, and Z of the global machine tool coordinate system are introduced. , , This corresponds to the slider node in the node pairing table Pairs. The third-order self-damped matrix With the corresponding guide rail node The third-order self-damped matrix and the slider node With guide rail node The coupling damping sub-matrix between and Its matrix representation is as follows: (31) (32) Among them, the equivalent translational stiffness coefficient preset at the interface , , With equivalent translational damping coefficient , , It can be obtained by conducting experimental modal tests and frequency response function parameter identification on the machine tool mating surface, or by analytical calculation based on contact mechanics theory combined with surface micro-morphology parameters, or by directly setting the rated contact stiffness and damping parameters provided by the supplier of moving parts such as linear guides.

[0042] Based on the slider nodes in the Pairs node pairing table With guide rail node The degree of freedom numbering is used to determine the stiffness matrix of the spring element. With the damping matrix of the damping unit The sub-block matrices in the matrix are summed to the overall unconstrained stiffness matrix of the whole machine. With respect to the unconstrained overall damping matrix of the whole machine middle.

[0043] Specifically, let the slider nodes in the node pairing table Pairs be... The global number is Guide rail node The global number is Since the finite element mesh model uses three-dimensional solid elements with three translational degrees of freedom per node, the slider node... The set of row and column indices corresponding to the overall matrix is: (33) Guide rail node The set of row and column indices corresponding to the overall matrix is: (34) Based on the above index set, perform sub-block addressing and accumulation processing of matrix elements: This involves processing the third-order self-stiffness submatrix... Accumulated into the overall unconstrained stiffness matrix of the machine The third-order self-damped submatrix is ​​positioned using row and column indices to determine its submatrix position. Accumulated into the unconstrained overall damping matrix of the whole machine The submatrix positions are determined by row and column indices; the coupling stiffness submatrix will be used. Accumulated into the overall unconstrained stiffness matrix of the machine As a row index and by The position of the submatrix determined by the column index will be used to couple the damping submatrix. Accumulated into the unconstrained overall damping matrix of the whole machine As a row index and by The position of the submatrix is ​​determined by the column index; another coupling stiffness submatrix... Accumulated into the overall unconstrained stiffness matrix of the machine As a row index and by The position of the submatrix determined by the column index will be used to locate another coupling damping submatrix. Accumulated into the unconstrained overall damping matrix of the whole machine As a row index and by The submatrix position is determined by the column index; the third-order self-stiffness submatrix corresponding to the guide rail node is... Accumulated into the overall unconstrained stiffness matrix of the machine The submatrix positions determined by the row and column indices will correspond to the third-order self-damped submatrix of the guide rail nodes. Accumulated into the unconstrained overall damping matrix of the whole machine The position of the submatrix is ​​determined by the row and column indices.

[0044] Therefore, the overall stiffness matrix of the unconstrained machine is... With respect to the unconstrained overall damping matrix of the whole machine The application of flexible connection at the bonding surface is completed in the middle.

[0045] Step 7.2, using the master-slave constraint transformation matrix calculated in step 5.2. The overall stiffness matrix of the unconstrained machine Unconstrained overall mass matrix of the whole machine and the overall unconstrained damping matrix of the whole machine By performing a condensation transformation to eliminate the slave node degrees of freedom constrained by the ball joint connection, a fixed connection of the internal components of the machine tool is achieved. Its condensation analytical expression is: (35) (36) (37) In the formula, This represents the overall stiffness matrix after the constraint condensation of the ball joint connection; This represents the overall mass matrix after the constraint condensation of the ball joint connection; This represents the overall damping matrix after the ball joint connection is constrained and condensed; This represents the master-slave constraint transformation matrix used for the condensation of the degree of freedom of the ball joint connection constraint; This represents the overall unconstrained stiffness matrix of the entire machine. This represents the unconstrained overall mass matrix of the entire machine. This represents the unconstrained overall damping matrix of the entire machine.

[0046] Step 7.3: Determine the new row and column indices based on the external fixed constraint table established in Step 2.5, and use the elimination method to refine the condensed global stiffness matrix. Overall quality matrix and the overall damping matrix Boundary condition processing is performed to fix the bottom surface of the machine tool bed to the external foundation, eliminating the degrees of freedom subject to external constraints, i.e., eliminating their corresponding rows and columns. This reduces the order of the system matrix and eliminates the matrix singularity caused by the rigid body displacement of the system, obtaining the final overall stiffness matrix after eliminating the influence of the rigid body displacement. Final overall quality matrix and the final overall damping matrix .

[0047] Step 8: Application of the dynamic model. Apply the final global stiffness matrix calculated in Step 7.3. Final overall quality matrix and the final overall damping matrix The data is imported into the solver for dynamic calculations. Following the sequence of G-code instructions, steps 3 through 8 are executed iteratively to obtain the evolution of the machine tool's dynamic characteristics along the machining trajectory. This method provides mechanical model support for chatter prediction during machine tool machining, optimization of cutting process parameters, and forward design of machine tool structures oriented towards dynamic performance.

[0048] The beneficial effects of this invention are as follows: (1) The inverse kinematics solution and rigid body space coordinate transformation method of machine tool driven by G code instructions proposed in this invention can realize the dynamic update of the coordinates of the finite element mesh nodes of each component, and effectively solve the problem that the finite element model is difficult to keep in line with the actual operating state of the machine tool, and can reflect the time-varying structural characteristics in the processing process.

[0049] (2) The element matrix reuse and spatial rotation transformation method proposed in this invention, compared with the traditional modeling method that requires re-grid and recalculate the element stiffness matrix and mass matrix for each configuration change, can avoid repeated numerical integration while ensuring the accuracy of the whole machine mechanical model, and significantly improve the computational efficiency of assembling the whole machine unconstrained overall stiffness matrix and overall mass matrix.

[0050] (3) The present invention uses a spatial search method based on the K-nearest neighbor algorithm to establish a node pairing table, and introduces spatial spring elements and spatial damping elements to apply flexible connections of the joint surface in the overall matrix. Compared with the traditional method of manually selecting joint surface node pairing, the method of the present invention has a higher degree of automation, which can effectively avoid the problems of cumbersome and inefficient modeling process, and realize automatic identification and rapid coupling of the contact state after each axis moves. Attached Figure Description

[0051] Figure 1 This is the overall flowchart of the method of the present invention; Figure 2 This is a schematic diagram of the initial configuration of a five-axis machining center; Figure 3 This is a schematic diagram of the component mass matrix reuse principle; Figure 4 This is a schematic diagram of the component stiffness matrix reuse principle; Figure 5 This is a schematic diagram of the inverse kinematics solution of a five-axis machining center. Figure 6 This is a schematic diagram illustrating the KNN algorithm for identifying slider-rail node pairings. Figure 7 This is a schematic diagram of the finite element model of the five-axis machining center output by this invention; Figure 8 This is a comparison chart of the efficiency and speed of finite element modeling in this invention; Figure 9 This is a schematic diagram of the finite element modeling efficiency test model of the present invention; In the diagram, 1 is a slider; 2 is a guide rail. Detailed Implementation

[0052] The following specific embodiments will further illustrate the present invention in detail.

[0053] A rapid finite element modeling (FEA) of a certain model of a five-axis swivel-type mill-turning machining center is used as an analysis example. In this five-axis swivel-type mill-turning machining center, the machine bed is fixed, the A-axis rotary table is mounted on the machine bed, the X-axis slide moves on the machine bed, the Y-axis slide moves on the X-axis slide, the Z-axis slide moves on the column of the Y-axis slide, and the B-axis swivel head is mounted on the Z-axis slide. This embodiment uses this mill-turning machining center as an example, but the method is applicable to other five-axis machining centers; it is only necessary to establish the corresponding kinematic model of the machine tool according to the specific machine tool structure and maintain consistent component division.

[0054] Based on the specific machine tool structure described above, the overall process of the rapid finite element modeling method for a five-axis machining center driven by G-code instructions is as follows: Figure 1 As shown, it includes the following steps: Step 1: The finite element mesh discretization method is used to generate meshes for each component of the five-axis machining center, establishing a finite element mesh model for each component under the initial configuration. Numerical integration is then used to calculate the element stiffness matrix and element mass matrix of each element under the initial configuration, and an element-matrix index table is established. Specifically: Step 1.1: A geometric model is established for the bed, X-axis slide, Y-axis slide, Z-axis slide, B-axis swivel head, slider, and guide rails of the five-axis swivel-head milling and turning machining center in this embodiment. The A-axis rotary table mainly serves as the moving component carrying the workpiece, and its structural flexibility has a relatively small impact on the overall dynamic characteristics of the machine tool. The A-axis structure is simplified in the finite element solid model of this embodiment; the actual rotation of the A-axis will be calculated through coordinate transformation using the subsequent inverse kinematics model. Each component of the geometric model is constrained to the zero point position of the machine tool's spatial coordinate system, with each motion axis of the machine tool in a zero-stroke state as the initial configuration of the machine tool, such as... Figure 2 As shown. Using finite element preprocessing software, the geometric model is meshed using hexahedral elements to establish finite element mesh models for each component, resulting in model files containing node spatial coordinates and element node numbers. The material properties of each element are then defined. In this embodiment, the main body of the component is made of isotropic cast iron, and its elastic modulus is set to [value missing]. Poisson's ratio is The material density is The completed finite element mesh model of the whole machine contains 144,661 nodes and a total of 114,483 elements.

[0055] Step 1.2: Based on the finite element mesh model established in Step 1.1, the element stiffness matrix and element mass matrix of each element under the initial configuration are calculated using the Gaussian numerical integration method. The calculation formula is as follows: (1) (2) In the formula, Indicates the unit number; This represents the element stiffness matrix under the initial configuration; This represents the element mass matrix under the initial configuration; This represents the strain-displacement matrix constructed based on shape functions; Represents the elasticity matrix of an isotropic material; The shape function matrix representing the element; Indicates the material density of the unit cell; This represents the volume integral region of the element. In this embodiment, the element mass matrix is ​​calculated using a lumped mass matrix method.

[0056] Step 1.3 converts the element stiffness matrix and element mass matrix of each element obtained in Step 1.2 under the initial configuration into a sparse matrix format. This step is performed only once during system initialization to establish and generate an element-matrix index table. In this embodiment, the element number, the corresponding eight finite element node numbers, the sparse storage data of the element stiffness matrix, and the sparse storage data of the element mass matrix are written into this element-matrix index table for direct reading and calling during subsequent G-code-driven model updates, eliminating the need to repeatedly perform numerical integration calculations of the element matrices.

[0057] Step 2: Machine Tool Kinematic Model Establishment and Mating Surface Node Index Table Construction. This step includes three parts, performed concurrently with Step 1: The first part is establishing a machine tool kinematic model describing the direction of the entire machine's series shaft system; the second part is automatically identifying boundary contact areas based on geometric features and constructing a mating surface node index table for each component; the third part is establishing a list of ball joint connections describing the internal fixed connection relationships of the machine tool, as well as an external fixed constraint table describing the external fixed constraints of the machine tool. Specifically: Step 2.1: Based on the axis sequence and transmission structure of the five-axis swivel-head milling and turning machining center, each axis of the machine tool is divided into two independent motion chains starting from the machine tool bed: a first independent motion chain from the machine tool bed to the workpiece coordinate system, and a second independent motion chain from the machine tool bed to the tool system. In this embodiment, the first independent motion chain includes the path from the machine tool bed to the A-axis rotary table; the second independent motion chain includes the path from the machine tool bed to the X-axis slide, Y-axis slide, Z-axis slide, and B-axis swivel head in sequence.

[0058] A global machine tool coordinate system is established with the machine tool bed as the global reference datum; a component local coordinate system is established at the rotation center of the B-axis swivel head; a component local coordinate system is established at the guide rail reference point under the initial configuration of the X-axis slide, Y-axis slide, and Z-axis slide; and a workpiece coordinate system describing the workpiece position and orientation is established at the workpiece design reference point.

[0059] Based on the actual shaft series topology of the machine tool in this embodiment, the first independent kinematic chain only includes the A-axis rotary table rotating around the X-axis, and the second independent kinematic chain only includes the B-axis oscillating head rotating around the Y-axis. Therefore, the corresponding transformation relationship is extracted from the third-order basic rotation matrix library to establish the specific model of this embodiment: the transformation matrix of the first rotating component. That is, the third-order fundamental rotation matrix of the A-axis and the transformation matrix of the second rotating component. This is the third-order fundamental rotation matrix along the B-axis, and its specific instance matrix form is as follows: (8) (9) In the formula, This indicates the target rotation angle of the A-axis driven by G-code instructions; This indicates the target rotation angle of the B-axis driven by G-code instructions.

[0060] Based on the first independent kinematic chain, the second independent kinematic chain, the global machine tool coordinate system, the component local coordinate system, and the workpiece coordinate system, a specific machine tool kinematic model is established for inverse kinematics calculation under the tool tip tracking control function. Its inverse kinematics geometric relationships are as follows: Figure 5 As shown. In machining mode with the tool tip point following control function enabled, the matrix from this embodiment is substituted into the geometric mapping relationship to establish a specific instance model: (10) By organizing the specific geometric relationships described in this embodiment, we obtain the position vector calculation model of the B-axis oscillating head rotation center relative to the global machine tool coordinate system under the current configuration: (11) In the formula, This represents the position vector of the rotation center of the rotating component relative to the global machine tool coordinate system under the current configuration; This represents the target position vector of the tool tip in the workpiece coordinate system; This represents the translation vector of the origin of the workpiece coordinate system in the global machine tool coordinate system; This represents the preset structural constants, which are the initial geometric vectors from the rotation center of the rotating component to the tool tip under the initial configuration, and are related to the length of the tool and tool holder. The current configuration of the machine tool is defined as the geometric pose state of each component after spatial translation and rotation transformation relative to the initial configuration following the execution of the current G-code instructions.

[0061] Step 2.2 involves sorting and independently numbering the mesh nodes in the finite element mesh model established in Step 1.1 to obtain a list containing the global numbers of all nodes. A component-node index table is then created, and the feature node regions of the slider bottom surface and the guide rail contact surface involved in the mating surface contact are marked in the component-node index table. In this embodiment, the entire finite element model has 144,661 nodes. The nodes of the B-axis oscillating head, Z / Y / X-axis sliding plates, column, slider, guide rail, and bed are sorted and recorded sequentially according to the globally unified numbering to construct the component-node index table shown in Table 1.

[0062] Table 1: Component-Node Index Table

[0063] Step 2.3: For the slider bottom surface node area and guide rail contact surface feature node area marked in Step 2.2, based on the planar or cylindrical surface geometric features of the component and user-specified feature nodes, automatic identification processing is performed using the spatial distance tolerance method to obtain a subset of planar contact nodes and a subset of cylindrical surface contact nodes. The planar contact node subset and the cylindrical surface contact node subset are then merged and deduplicated to establish and generate a mating surface node index table. The user-specified feature nodes include vertex nodes used to define the spatial boundary of the planar contact area of ​​the component, and circumferential nodes used to fit the geometric features of the cylindrical contact surface of the component. The mating surface node index is obtained using feature nodes, specifically as follows: Planar Node Extraction: A three-dimensional spatial plane equation is fitted based on three or more non-collinear specified feature nodes. In this embodiment, spatial reference node features such as globally numbered 1 and 17181, and 55687 and 111305 are used as the plane fitting benchmark. All nodes on the corresponding component are traversed one by one, and the spatial vertical distance from the node to the plane is calculated. Nodes whose spatial vertical distance is not greater than a preset surface tolerance are included in the planar contact node subset. The preset surface tolerance ranges from 10% to 30% of the average edge length of the finite element mesh in the component contact area, and the tolerance is set to 1.375 mm.

[0064] Cylindrical Surface Node Extraction: Based on the cylinder axis direction, reference point on the axis, and cylinder radius, the radial distance from the node to the cylinder axis is calculated. Nodes whose absolute value of the difference between the radial distance and the cylinder radius is not greater than a preset circular tolerance are included in the cylindrical surface contact node subset. The preset circular tolerance ranges from 10% to 30% of the average edge length of the finite element mesh in the component contact area, with a tolerance set to 1.375 mm. The cylinder axis direction, reference point on the axis, and cylinder radius are determined by the cylindrical surface geometric parameters in the geometric model. Merging and Deduplication: The planar contact node subsets of different mating areas on the same component surface are merged with the cylindrical surface contact node subset, duplicate node numbers are removed, and the mating surface node index table is determined.

[0065] Step 2.4: Obtain ball joint connection pair information describing the fixed connection relationships between various components inside the machine tool, and establish a ball joint connection pair list. In this embodiment, for internal component connections that do not require relative sliding, the corresponding coincident nodes are extracted as ball joint connection pairs. This embodiment extracts and constructs a one-to-one correspondence such as global node numbers 16874 and 15279, 16875 and 15190, which together form a ball joint connection pair list reflecting the fixed assembly constraints of the components.

[0066] Step 2.5: Obtain the fixed boundary conditions connecting the machine tool bed bottom surface to the external foundation, the set of degrees of freedom subject to external fixed constraints, and establish an external fixed constraint table. In this embodiment, for the anchor bolt installation positions where the machine tool bed is fixed to the foundation, extract the global numbers of the finite element mesh nodes in the corresponding area. In this embodiment, the constrained nodes, including the bed bottom surface nodes with global node numbers 135115, 135222, and 135113, are extracted to completely restrict the external fixed constraint table for the three translational degrees of freedom of this part of the node space.

[0067] Step 3: G-code instruction parsing and node coordinate update based on inverse kinematics. This step follows directly from Step 2. The G-code instructions are read using the G-code instruction parsing method, and the machining trajectory contained in the G-code instructions is converted into the relative translational displacement and rotational attitude of each machine tool component. The current configuration of the machine tool is defined as the geometric pose of each component after spatial translation and rotation transformation relative to the initial configuration following the execution of the current G-code instruction. The spatial coordinates of the mesh nodes of each component are updated using the rigid body spatial coordinate transformation formula. Specifically: Step 3.1: Parse the G-code instructions line by line. The parsing method involves scanning the function words representing the machine tool's motion mode, including linear interpolation or rapid positioning. Based on the identified function words, extract the coordinate address characters representing the spatial geometric position and rotation axis attitude in the current line. For the five-axis swivel-head milling and turning machining center in this embodiment, extract the data words including X, Y, Z, A, B, and their suffixes. Convert the extracted characters into real-number numerical parameters to obtain the target position vector of the tool tip in the workpiece coordinate system. And the target rotation angle values ​​of the A-axis turntable and the B-axis oscillating head. , .

[0068] Step 3.2, based on the tool tip tracking control principle, the target position vector of the tool tip in the workpiece coordinate system calculated in Step 3.1 is... The target rotation angle values ​​of each rotation axis are input into the specific machine tool kinematic model established in step 2.1, and inverse kinematics calculations are performed to determine the indexed values. Spatial translation vector of machine tool components and rotation matrix In this embodiment, when the machine tool is in machining mode with the tool tip tracking control function enabled, the kinematic model established in step 2.1 is invoked, and the specific physical dimensional features of this embodiment are extracted and substituted. The preset structural constant in this embodiment is the initial geometric vector from the B-axis oscillating head rotation center to the tool tip under the initial configuration. A three-dimensional space vector The unit is millimeters. This constant is formed by superimposing the tool length vector (with a negative direction) and the spindle offset vector; the origin translation vector of the workpiece coordinate system in the global machine tool coordinate system. A three-dimensional space vector The specific calculation model for its inverse kinematics solution is as follows: (12) In the formula, This represents the position vector of the rotation center of the rotating component, i.e., the B-axis oscillating head in this embodiment, relative to the global machine tool coordinate system under the current configuration; This represents the target position vector of the tool tip in the workpiece coordinate system. Here, we are referring to the A-axis rotary table, which carries and drives the workpiece to rotate; its equivalent rotation angle in the global machine tool coordinate system is taken as... To maintain alignment with the positive direction of rotation of the machine tool rotary table.

[0069] Using position vectors The spatial translation vector of the B-axis oscillating head, which is a rotating component, is obtained by calculating the difference between the initial position vector of the B-axis oscillating head and the initial position vector of the rotation center under the initial configuration. : (13) In the formula, Indicates having an index number The spatial translation vector of the machine tool component in the global machine tool coordinate system; This represents the initial position vector of the rotating component, i.e., the B-axis oscillating head in this embodiment, relative to the global machine tool coordinate system under the initial configuration. The specific coordinate data of the initial position vector in this embodiment are as follows: ; This indicates the index number of the machine tool component. For the pure translation components in this embodiment, namely the X-axis slide, Y-axis slide, and Z-axis slide, the stroke displacement of each axis relative to its initial configuration, obtained from the inverse kinematics operation of the current G-code instruction, is directly used as the spatial translation vector of the pure translation component.

[0070] Step 3.3, analyze the target rotation angle values ​​of each rotation axis obtained in Step 3.1. , Substitute these values ​​into the third-order basic rotation matrix formula. Based on the shaft series sequence of the machine tool in this embodiment, calculate the rotation matrix corresponding to each component. For the A-axis rotary table, its corresponding rotation matrix is... This is the third-order fundamental rotation matrix for rotation around the X-axis; for the B-axis rotation, its corresponding rotation matrix is... This refers to the third-order basic rotation matrix for rotation around the Y-axis; for the machine tool bed, column, X-axis slide, Y-axis slide, and Z-axis slide, since they do not undergo spatial orientation rotation, their corresponding rotation matrices are... All are third-order identity matrices.

[0071] Step 3.4: Update the node coordinates in the finite element mesh models of each component in Step 1.1 using the following rigid body space coordinate transformation formula: For purely translational components, namely the machine tool bed, column, X-axis slide, Y-axis slide, and Z-axis slide, the node coordinate update formula is as follows: (14) For the rotating components, namely the A-axis turntable and the B-axis oscillating head, the nodal coordinate update formula is as follows: (15) In the formula, Indicates machine tool components Internal finite element node numbering; Indicates machine tool components Internal The global machine tool coordinate vector of each finite element node in the current configuration; Indicates machine tool components Internal The global machine tool coordinate vector of each finite element node under the initial configuration; Indicates having an index number The spatial translation vector of the machine tool component in the global machine tool coordinate system; Indicates having an index number The rotation matrix corresponding to the machine tool component; This represents the initial position vector of the rotation center of the corresponding rotating component in the global machine tool coordinate system under the initial configuration; This represents the position vector of the rotation center of the corresponding rotating component relative to the global machine tool coordinate system. Here, for the A-axis rotary table, the position vector of the rotation center of the aforementioned rotating component relative to the global machine tool coordinate system and the initial position vector are respectively the geometric vectors of the A-axis rotary table under its own configuration.

[0072] Using the rigid body space coordinate transformation formula, the global machine tool coordinate vectors of all mesh nodes of each component in the whole machine finite element model under the current configuration are calculated and obtained. Considering the mesh node interval distribution in this embodiment, when specifically performing node coordinate updates: firstly, for the B-axis swaying head mesh nodes within the range of global node numbers 1 to 13845, the coordinates are updated around the three-dimensional coordinate system according to the aforementioned rotating component coordinate update formula. The rotation center is rotated in millimeters to update the spatial rotation and translation coordinates; then, for the Z-axis slide mesh nodes in the range of global node numbers 13846 to 21242, the Y-axis slide mesh nodes in the range of numbers 21243 to 62678, and the X-axis slide mesh nodes in the range of numbers 62679 to 79859, the corresponding component translation vectors are applied according to the above pure translation component coordinate update formula to complete the spatial coordinate update; the current configuration of the machine tool is finally determined by the global machine tool coordinate vector of the updated mesh nodes under the current configuration.

[0073] Step 4: Automatic node pairing at the interface based on the K-nearest neighbor algorithm. Using the spatial search method and the current machine tool configuration calculated in Step 3.4, matching nodes are automatically found among the nodes of the components that need to be connected. The contact state after the movement of each axis is determined, and a node pairing table, Pairs, is established to describe the connection between the slider and guide rail interface. A schematic diagram of the KNN algorithm identifying slider-guide rail node pairing is shown below. Figure 6 As shown, Figure 6 In the diagram, 1 represents the slider and 2 represents the guide rail. Specifically: Step 4.1: For the slider-guide rail pair, extract the set of slider nodes involved in the contact based on the mating surface node index table established in Step 2.3. With guide rail node set In this embodiment, the slider node set is extracted. The data sources include the underlying index records of each axis slider and the node records of the B-axis oscillating head connection area; the guide rail node set is extracted. The data sources include the underlying index records of the guide rail and the node records of the backplane. The above data are combined to form a complete set of nodes participating in the pairing calculation.

[0074] Step 4.2, denote the global machine tool coordinate vector of the corresponding slider node obtained in step 3.4 as... The global machine tool coordinate vector of the corresponding guide rail node is denoted as .

[0075] Step 4.3: Assign the global machine coordinate vectors of the slider node and guide rail node under the current configuration. and As the retrieval benchmark, the K-nearest neighbor method with a search count of 1 is used for spatial search and distance calculation to obtain the slider nodes. With guide rail node Spatial Euclidean distance between : (16) In the formula, Represents slider node With guide rail node The spatial Euclidean distance between them; Represents slider node The global machine tool coordinate vector under the current configuration; Indicates guide rail node Global machine tool coordinate vector in the current configuration; superscript This represents the matrix transpose operation. It applies to the calculated spatial Euclidean distance. The comparison is performed with a preset interface threshold; if the spatial Euclidean distance is... If the distance between two nodes is not greater than the preset joint surface threshold, the two nodes are determined to form a joint surface connection, and the corresponding relationship is added to the node pairing table Pairs. The preset joint surface threshold is set to 0.15mm based on the actual mesh accuracy. If the same guide rail node is repeatedly matched by multiple slider nodes, only a unique correspondence is retained according to the principle of minimum spatial Euclidean distance.

[0076] Step 4.4: Collect all corresponding relationships that satisfy the mating surface connection and establish a node pairing table Pairs. In this embodiment, the system merges the slider node array index obtained based on the nearest neighbor matching algorithm with the corresponding guide rail node array index to form a node pairing table Pairs that reflects the mating surface connection state of the machine tool in the current machining pose.

[0077] Step 5: Deduplication and preprocessing of master-slave node degrees of freedom based on ball joint connections. The constraints of the ball joint connections between internal machine tool components are processed according to deformation compatibility conditions to fix the two internal components and maintain their common motion relationship, establishing a transformation relationship to eliminate the slave node degrees of freedom. Specifically: Step 5.1: Read the list of ball joint connections established in Step 2.4; perform uniqueness verification and deduplication on the master node degrees of freedom and slave node degrees of freedom corresponding to the ball joint connection pairs to ensure that each slave node degree of freedom in the whole machine finite element model corresponds to only one master node degree of freedom; divide the total degree of freedom of the whole machine finite element model into the set of eliminated slave node degrees of freedom and the set of retained degrees of freedom composed of master node degrees of freedom and unconstrained degrees of freedom.

[0078] In this embodiment, for the fixed assembly constraints inside the machine tool, a list of ball joint connection pairs containing corresponding relationships such as global node numbers 16874 and 15279 is extracted in sequence. Redundant duplicate constraint pairs are identified and removed. The total of 433983 translational degrees of freedom of the whole machine is divided into slave node degrees of freedom that will be condensed and master node and unconstrained degrees of freedom that are retained for calculation.

[0079] Step 5.2: For the deduplicated linear constraint equations of the ball joint connection, each slave node degree of freedom is written as an algebraic linear function of the corresponding master node degree of freedom, thus establishing a multi-point constraint equation system whose block algebraic form satisfies: (17) Based on the algebraic coefficients of the multi-point constraint equation system, a master-slave constraint transformation matrix for the reduction of the degree of freedom of the ball joint connection constraint is established by expanding and combining the equations. Its explicit parsing expression is: (18) In the formula, This represents the constraint coefficient matrix corresponding to the degrees of freedom of each node; This represents the vector of slave node degrees of freedom, which consists of the slave node degrees of freedom. This represents the constraint coefficient matrix corresponding to the retained degrees of freedom; This represents the vector of retained degrees of freedom, consisting of the master node degrees of freedom and the unconstrained degrees of freedom. Represents the master-slave constraint transformation matrix; This represents the identity matrix of the same order as the number of degrees of freedom retained. The number of rows in this matrix is ​​equal to the original total number of degrees of freedom of the system. In this specific embodiment, it has 433,983 rows, and the number of columns is equal to the number of degrees of freedom retained after eliminating the degrees of freedom of slave nodes. It is used in subsequent steps to shrink the overall matrix to an algebraic space containing only the degrees of freedom of master nodes and unconstrained degrees of freedom.

[0080] Step 6: Overall matrix assembly based on matrix reuse and rotation transformation. Utilizing the element-matrix index table generated in Step 1.3 and the rotation matrix calculated in Step 3.3... Based on the energy equivalence principle of coordinate transformation, the element matrix is ​​reused, assembled, and spatially rotated to form the unconstrained overall stiffness matrix and the unconstrained overall mass matrix of the entire machine. Then, based on the matrix order of the unconstrained overall stiffness matrix and the unconstrained overall mass matrix, the unconstrained overall damping matrix of the entire machine is established. Specifically: Step 6.1, targeting the moving parts inside the machine tool Obtain its rotation matrix from step 3.3. If components If it is a rotating component, then its third-order rotation matrix is ​​used. As basic unit items, they are laid out flat along the main diagonal without overlap. This establishes a size scale of block diagonal rotation transformation matrix ,in Indicates components The total number of nodes. In this embodiment, for the B-axis oscillating head, which is a rotating component, the total number of grid nodes it contains is... The value is 13845. The system uses its corresponding third-order fundamental rotation matrix around the Y-axis to tile and combine along the main diagonal, constructing a block diagonal rotation transformation matrix of order 41535; if the component If it is a purely translational component, then its rotation matrix... It is a third-order identity matrix, and its corresponding block diagonal rotation transformation matrix. for An identity matrix of order 1.

[0081] Step 6.2: Extract components from the unit-matrix index table established in Step 1.3. The element stiffness matrix and element mass matrix under the initial configuration are assembled into the initial component stiffness matrix. and the initial component mass matrix For rotating components, a spatial rotational transformation correction is applied using the derivation process based on the energy equivalence principle of continuum mechanics as described in steps 6.2.1 to 6.2.2, and the global stiffness matrix of the current configuration after the rotational transformation is calculated. With the current configuration global mass matrix The reuse principles of the mass matrix and stiffness matrix of the components are as follows: Figure 3 , Figure 4 As shown. Specifically: Step 6.2.1: Based on the equivalence principle of system strain energy under coordinate transformation, when a component undergoes rigid body rotation, its inherent elastic characteristics remain unchanged relative to the initial configuration. Introduce the global displacement vector under the current configuration. and the nodal displacement vectors in the component's local coordinate system It satisfies the following strain energy invariance formula: (19) Based on the block diagonal rotation transformation matrix The orthogonal geometric characteristics, since the global displacement vector under the current configuration is the product of the block diagonal rotation transformation matrix and the initial displacement vector, satisfy the relation... Then the initial displacement vector can be expressed as Substituting the aforementioned correspondence into the strain energy invariance formula, the transformation formula is obtained by expansion. (20) Based on this, the rotation transformation formula for the stiffness matrix is ​​derived as follows: (twenty one) Step 6.2.2 Introduce the global velocity vector under the current configuration and the nodal velocity vectors in the component's local coordinate system It satisfies the following kinetic energy invariance formula: (twenty two) Since the initial velocity vector is equivalent to the product of the transpose of the block diagonal rotation transformation matrix and the global velocity vector, it satisfies the following relationship: Substituting the corresponding relationships into the kinetic energy invariance formula, the rotation transformation formula for the mass matrix is ​​derived as follows: (twenty three) In the formula, Indicates components Initial component stiffness matrix under initial configuration; Indicates components Initial component mass matrix under the initial configuration; Indicates components Size is The block diagonal rotation transformation matrix; Indicates components The global stiffness matrix of the current configuration after rotational transformation; Indicates components Global mass matrix of the current configuration after rotational transformation; superscript This represents the matrix transpose operation; This indicates the index number of the machine tool component. For purely translational components, due to their block diagonal rotation transformation matrix... Since it is an identity matrix, no spatial rotation transformation correction is needed. The initial matrix, i.e., the global stiffness matrix of the current configuration, can be directly reused according to the above formula. Equal to the initial component stiffness matrix The current configuration global quality matrix equal to the initial component mass matrix .

[0082] Step 6.3: Obtain the global number of each mesh node using the component-node index table established in Step 2.2; based on the matrix assembly principle of the finite element direct stiffness method, perform global addressing and accumulation processing of the stiffness matrix and mass matrix: for the unconstrained overall stiffness matrix of the whole machine Based on the global number of the mesh nodes, the global stiffness matrix of each component in its current configuration is determined. The matrix elements in the matrix are accumulated into the overall unconstrained stiffness matrix of the entire machine. In the process, the unconstrained overall stiffness matrix of the entire machine is formed through assembly. Its assembly mathematical formula is: (twenty four) For the unconstrained overall mass matrix of the whole machine Based on the global number of the grid nodes, the global mass matrix of each component's current configuration is calculated. The matrix elements in the matrix are accumulated into the unconstrained overall mass matrix of the whole machine. In the process, the unconstrained overall mass matrix of the entire machine is formed through assembly. Its assembly mathematical formula is: (25) In the formula, Represents the overall unconstrained stiffness matrix of the entire machine. The Middle line, number Matrix elements at column positions; Represents the unconstrained overall mass matrix of the entire machine. The Middle line, number Matrix elements at column positions; This indicates the total number of independent moving parts contained in the machine tool; Indicates the index number of the machine tool component; Indicates having an index number The machine tool component in the global stiffness matrix under the current configuration is the first... line, number Matrix elements at column positions; Indicates having an index number The machine tool component in the global mass matrix under the current configuration is the [number]th [unit]. line, number Matrix elements at column positions; and These represent the global row number and global column number of the overall matrix, respectively; and These represent the internal row and column numbers of the global stiffness matrix or global mass matrix of the corresponding component for the current configuration, respectively. Since the finite element mesh model uses three-dimensional solid elements with 3 translational degrees of freedom per node, if the global number of a certain mesh node is... Then the global row number corresponding to the grid node in the overall matrix or global column number The set of values ​​is , , This achieves the mapping from the global numbering of grid nodes to the row and column positions of the overall matrix. In this embodiment, the B-axis oscillating head global stiffness matrix containing 41,535 internal degrees of freedom, along with the matrix elements of other independent moving parts such as the skateboard, are cumulatively added to the unconstrained overall stiffness matrix of the whole machine with an order of 433,983 × 433,983 according to this mapping rule. Unconstrained overall mass matrix of the whole machine middle.

[0083] Step 6.4: Based on the matrix orders of the unconstrained overall stiffness matrix and the unconstrained overall mass matrix of the whole machine assembled in Step 6.3, establish an unconstrained overall damping matrix of the same order, and initialize the unconstrained overall damping matrix of the whole machine to a zero matrix; the unconstrained overall damping matrix of the whole machine is used to receive the accumulated terms of the damping matrix of the damping unit corresponding to the joint surface connection in subsequent steps, and its initialization expression is: (26) In the formula, This represents the unconstrained overall damping matrix of the entire machine after initialization. This represents the order of the unconstrained overall stiffness matrix and the unconstrained overall mass matrix of the entire machine, combined with the order of this ultra-large matrix in this embodiment. The strict setting is 433983; express A zero matrix of order 1.

[0084] Step 7: Applying multiple constraints and system matrix shrunk processing. This involves combining the node pairing table (Pairs) established in Step 4.4 with the master-slave constraint transformation matrix established in Step 5.2. The unconstrained overall stiffness matrix of the whole machine assembled in step 6.3 Unconstrained overall mass matrix of the whole machine and the unconstrained overall damping matrix of the whole machine established in step 6.4. Multiple constraint processing is performed; through the overall unconstrained stiffness matrix of the whole machine. With respect to the unconstrained overall damping matrix of the whole machine A flexible connection is applied to the joint surface, and the degrees of freedom from the nodes and those subject to external fixed constraints are eliminated sequentially to obtain the final overall stiffness matrix after eliminating the influence of rigid body displacements in the system. Final overall quality matrix and the final overall damping matrix Specifically: Step 7.1 involves connecting the mating surfaces of the nodes in the Pairs table established in Step 4.4, i.e., the corresponding slider nodes. With guide rail node By introducing spatial spring elements that simulate contact stiffness characteristics and spatial damping elements that simulate energy dissipation characteristics, the corresponding spring element stiffness matrix is ​​established. With the damping matrix of the damping unit The damping matrix of the damping unit is used to accumulate in the overall unconstrained damping matrix of the whole machine established in step 6.4. (27) (28) In the formula, This represents the stiffness matrix of the spring element established by the translational stiffness coefficient of the mating surface; , , , These are all stiffness sub-matrices used to describe the stiffness characteristics of the physical connection between paired nodes; specifically, equivalent translational stiffness coefficients are introduced from the mating surface in the three spatial translational directions X, Y, and Z of the global machine tool coordinate system. , , In this embodiment, the stiffness coefficients in the three translational directions are all set to 1. Corresponding to the node pairing table Pairs slider node The third-order self-stiffness submatrix With the corresponding guide rail node The third-order self-stiffness submatrix and the slider node With guide rail node The coupling stiffness submatrix between and Its matrix representation is as follows: (29) (30) in, This represents the damping matrix of the damping element established by the translational damping coefficient of the interface; , , , These are all damping sub-matrices used to describe the energy dissipation damping characteristics between paired nodes; specifically, equivalent translational damping coefficients of the joint surface in the three spatial translational directions X, Y, and Z of the global machine tool coordinate system are introduced. , , In this embodiment, its corresponding value is set to This corresponds to the slider node in the node pairing table Pairs. The third-order self-damped matrix With the corresponding guide rail node The third-order self-damped matrix and the slider node With guide rail node The coupling damping sub-matrix between and Its matrix representation is as follows: (31) (32) Based on the slider nodes in the Pairs node pairing table With guide rail node The degree of freedom numbering is used to determine the stiffness matrix of the spring element. With the damping matrix of the damping unit The sub-block matrices in the matrix are summed to the overall unconstrained stiffness matrix of the whole machine. With respect to the unconstrained overall damping matrix of the whole machine In detail, let the slider node be in the node pairing table Pairs. The global number is Guide rail node The global number is Since the finite element mesh model uses three-dimensional solid elements with three translational degrees of freedom per node, the slider node... The set of row and column indices corresponding to the overall matrix is: (33) Guide rail node The set of row and column indices corresponding to the overall matrix is: (34) Sub-block addressing and accumulation of matrix elements based on the index set: This involves processing the third-order self-stiffness submatrix... Accumulated into the overall unconstrained stiffness matrix of the machine The third-order self-damped submatrix is ​​positioned using row and column indices to determine its submatrix position. Accumulated into the unconstrained overall damping matrix of the whole machine The submatrix positions are determined by row and column indices; the coupling stiffness submatrix will be used. Accumulated into the overall unconstrained stiffness matrix of the machine As a row index and by The position of the submatrix determined by the column index will be used to couple the damping submatrix. Accumulated into the unconstrained overall damping matrix of the whole machine As a row index and by The position of the submatrix is ​​determined by the column index; another coupling stiffness submatrix... Accumulated into the overall unconstrained stiffness matrix of the machine As a row index and by The position of the submatrix determined by the column index will be used to locate another coupling damping submatrix. Accumulated into the unconstrained overall damping matrix of the whole machine As a row index and by The submatrix position is determined by the column index; the third-order self-stiffness submatrix corresponding to the guide rail node is... Accumulated into the overall unconstrained stiffness matrix of the machine The submatrix positions determined by the row and column indices will correspond to the third-order self-damped submatrix of the guide rail nodes. Accumulated into the unconstrained overall damping matrix of the whole machine The submatrix positions are determined by row and column indices. This leads to the overall unconstrained stiffness matrix of the entire machine. With respect to the unconstrained overall damping matrix of the whole machine The application of flexible connection at the bonding surface is completed in the middle.

[0085] Step 7.2, using the master-slave constraint transformation matrix calculated in step 5.2. The overall stiffness matrix of the unconstrained machine Unconstrained overall mass matrix of the whole machine and the overall unconstrained damping matrix of the whole machine By performing a condensation transformation to eliminate the slave node degrees of freedom constrained by the ball joint connection, a fixed connection of the internal components of the machine tool is achieved. Its condensation analytical expression is: (35) (36) (37) In the formula, This represents the overall stiffness matrix after the constraint condensation of the ball joint connection; This represents the overall mass matrix after the constraint condensation of the ball joint connection; This represents the overall damping matrix after the ball joint connection is constrained and condensed; This represents the master-slave constraint transformation matrix used for the condensation of the degree of freedom of the ball joint connection constraint; This represents the overall unconstrained stiffness matrix of the entire machine. This represents the unconstrained overall mass matrix of the entire machine. This represents the unconstrained overall damping matrix of the entire machine.

[0086] Step 7.3: Determine the new row and column indices based on the external fixed constraint table established in Step 2.5, and use the elimination method to refine the condensed global stiffness matrix. Overall quality matrix and the overall damping matrix Boundary condition processing is performed to fix the bottom surface of the machine tool bed to the external foundation, eliminating the degrees of freedom subject to external constraints, i.e., eliminating their corresponding rows and columns. This reduces the order of the system matrix and eliminates the matrix singularity caused by the rigid body displacement of the system, obtaining the final overall stiffness matrix after eliminating the influence of the rigid body displacement. Final overall quality matrix and the final overall damping matrix .

[0087] Step 8: Application of the dynamic model. Apply the final global stiffness matrix calculated in Step 7.3. Final overall quality matrix and the final overall damping matrix The data is imported into the solver for dynamic calculations. Following the sequence of G-code instructions, steps 3 through 8 are executed iteratively to obtain the evolution of the machine tool's dynamic characteristics along the machining trajectory. This method provides mechanical model support for chatter prediction during machine tool machining, optimization of cutting process parameters, and forward design of machine tool structures oriented towards dynamic performance.

[0088] As a verification of the modeling effect in this embodiment, when the G-code instruction is X531.3684Y-294.9790Z248.5775B45.000A307.000, the finite element model of the five-axis machining center output by the current G-code configuration is as follows: Figure 7 As shown. Figure 7 The connection annotations on the mating surfaces show the node matching relationship between the guide rail and the slider component under this specific configuration. This result demonstrates that the spatial distance matching method described in this embodiment can obtain the paired nodes of the mating surfaces and apply corresponding flexible connection constraints after the machine tool's pose transformation, ensuring the integrity of the overall finite element assembly topology.

[0089] In addition, to evaluate the efficiency of this method, six processing trajectory points were extracted for multi-configuration reconstruction tests, and the corresponding whole-machine finite element model is as follows: Figure 9 As shown. Figure 9 The connection states of the mating surfaces presented in each model verify the applicability of the aforementioned node pairing algorithm under the variable configuration state of the machine tool. The comparison results of the finite element modeling calculation time for the above six configurations are as follows: Figure 8 As shown in the test data, the method improves the efficiency of modeling the whole machine finite element model.

[0090] The above embodiments are merely illustrative of the implementation methods of the present invention, but should not be construed as limiting the scope of the present invention. It should be noted that those skilled in the art can make various modifications and improvements without departing from the concept of the present invention, and these modifications and improvements all fall within the protection scope of the present invention.

Claims

1. A rapid finite element modeling method for five-axis machining centers based on G-code instruction drive, characterized in that, Includes the following steps: Step 1: Component-specific mesh generation and element matrix pre-calculation; The finite element mesh discretization method is used to mesh each component of the five-axis machining center to establish the finite element mesh model of each component under the initial configuration; the element stiffness matrix and element mass matrix of each element under the initial configuration are calculated using the numerical integration method, and an element-matrix index table is established. Step 2: Establishment of machine tool kinematic model and construction of mating surface node index table; Establish a kinematic model of the machine tool describing the direction of the series shaft system of the whole machine; automatically identify the boundary contact area based on geometric features and construct an index table of mating surface nodes of each component; establish a list of ball joint connections describing the internal fixed connection relationship of the machine tool, and an external fixed constraint table describing the external fixed constraint of the machine tool. Step 3: G-code instruction parsing and node coordinate update based on inverse kinematics; The G-code instructions are read according to the G-code instruction parsing method, and the machining trajectory contained in the G-code instructions is converted into the relative translational displacement and rotational attitude of each component of the machine tool. The current configuration of the machine tool is defined as the geometric pose of each component after spatial translation and rotation transformation relative to the initial configuration following the execution of the current G-code instruction; the spatial coordinates of the mesh nodes of each component are updated using the rigid body spatial coordinate transformation formula. Step 4: Automatic pairing of interface nodes based on the K-nearest neighbor algorithm; Based on the spatial search method, using the current configuration of the machine tool, matching nodes are automatically found among the nodes of the components that need to be connected, the contact state after each axis moves is determined, and a node pairing table Pairs is established to describe the connection of the slider-guide rail mating surface. Step 5: Deduplication and preprocessing of master-slave node degrees of freedom based on ball joint connection pairs; Based on the deformation compatibility conditions, the ball joint connection between the internal components of the machine tool is constrained to fix the two internal components and maintain their common motion relationship. A transformation relationship to eliminate the degree of freedom of the slave node is established, and the master-slave constraint transformation matrix is ​​obtained. Step 6: Overall matrix assembly based on matrix reuse and rotation transformation; The element matrix is ​​reused, assembled, and spatially rotated to form the unconstrained overall stiffness matrix and the unconstrained overall mass matrix of the whole machine, and then the unconstrained overall damping matrix of the whole machine is established. Step 7: Applying multiple constraints and performing system matrix condensation processing; By combining the node pairing table Pairs and the master-slave constraint transformation matrix, multiple constraint processing is applied to the overall unconstrained stiffness matrix, overall unconstrained mass matrix, and overall unconstrained damping matrix of the whole machine. By applying flexible connections to the joint surfaces and successively eliminating the degrees of freedom of the slave nodes and the degrees of freedom subject to external fixed constraints, the final overall stiffness matrix, final overall mass matrix, and final overall damping matrix after eliminating the influence of rigid body displacement of the system are calculated. Step 8: Application of the dynamic model.

2. The rapid finite element modeling method for five-axis machining centers based on G-code instruction drive according to claim 1, characterized in that, Specifically, step 1 is as follows: Step 1.1: Establish geometric models for each component of the five-axis machining center, constrain each component of the geometric model to the zero point position, and take the zero stroke of each motion axis of the machine tool as the initial configuration of the machine tool. The components include the bed, slide, column, slider, guide rail, swivel head or rotary table; use the mesh generation method to mesh the geometric model with three-dimensional solid elements, establish the finite element mesh model of each component, obtain the model file containing node coordinates and element node numbers, and define the material properties of each element; Step 1.2: Based on the finite element mesh model established in Step 1.1, the element stiffness matrix and element mass matrix of each element under the initial configuration are calculated using the Gaussian numerical integration method. The calculation formulas are as follows: In the formula, Indicates the unit number; Representation unit The element stiffness matrix under the initial configuration; Representation unit The unit mass matrix under the initial configuration; Representation unit The strain-displacement matrix; Representation unit The elasticity matrix; Representation unit The shape function matrix; Representation unit Material density; Representation unit The volume integral region; Step 1.3: Convert the element stiffness matrix and element mass matrix of each element under the initial configuration into a sparse format. Perform this once during the initialization phase to establish and generate the element-matrix index table.

3. The rapid finite element modeling method for five-axis machining centers based on G-code instruction drive according to claim 2, characterized in that, Step 2 specifically includes: Step 2.1: Divide each axis of the machine tool into two independent motion chains starting from the machine tool bed: the first independent motion chain pointing from the machine tool bed to the workpiece coordinate system, and the second independent motion chain pointing from the machine tool bed to the tool system. A global machine tool coordinate system is established with the bed as the global reference datum; for rotating parts, a local coordinate system is established at their rotation center; for slide parts, a local coordinate system is established at the guide rail reference point under their initial configuration; for columns, a local coordinate system is established at the reference reference point where their bottom surface connects to the bed; and a workpiece coordinate system describing the workpiece position and orientation is established at the workpiece design reference point. Construct a library of third-order basic rotation matrices for the A, B, and C axes of a machine tool. Based on the actual series topology of the axis system in a five-axis machining center, the first independent motion chain... The third-order fundamental rotation matrix corresponding to each rotation axis and the second independent kinetic chain The third-order fundamental rotation matrix corresponding to each rotation axis The third-order basic rotation matrix library is obtained through mapping rules. Select from: In the formula, This represents the selected fundamental rotation matrix. ; Then, the third-order basic rotation matrices of each rotating component in the first independent motion chain are multiplied and cascaded in order from the bed to the workpiece coordinate system to establish the transformation matrix of the first rotating component. The third-order basic rotation matrices of each rotating component in the second independent kinematic chain are multiplied and cascaded in order from the bed to the tool system to establish the transformation matrix of the second rotating component. The formula for its calculation is: In the formula, This represents the transformation matrix of the first rotating component; This indicates the total number of rotational axes contained in the first independent kinematic chain; Indicates the cascade index number of the rotating axis in the first independent kinematic chain; In the first independent kinetic chain, the first The third-order basic rotation matrix corresponding to each rotation axis; This represents the transformation matrix of the second rotating component; This indicates the total number of rotational axes contained in the second independent kinematic chain; Indicates the cascade index number of the rotation axis in the second independent kinematic chain; Indicating the second independent kinetic chain, the first The third-order basic rotation matrix corresponding to each rotation axis; A machine tool kinematic model is established for inverse kinematics calculation under the tool tip tracking control function; the machine tool kinematic model is used to describe the geometric mapping relationship between the tool tip target position vector in the workpiece coordinate system, the origin translation vector of the workpiece coordinate system in the global machine tool coordinate system, the position vector of the rotation center of the rotating part relative to the global machine tool coordinate system under the current configuration, the first rotating part transformation matrix, the second rotating part transformation matrix, and the preset structural constants. In machining mode with the tool tip tracking control function enabled, the expression of the tool tip target position vector in the global machine tool coordinate system is consistent with the actual tool tip position of the tool system. The actual tool tip position of the tool system is determined by the position vector of the rotation center of the rotating component relative to the global machine tool coordinate system under the current configuration and the preset structural constant after transformation by the second rotating component transformation matrix. The geometric relationship is as follows: Based on formula (10), the position vector of the rotation center of the rotating component relative to the global machine tool coordinate system under the current configuration is obtained: In the formula, This represents the position vector of the rotation center of the rotating component relative to the global machine tool coordinate system under the current configuration; This represents the transformation matrix of the first rotating component; This represents the target position vector of the tool tip in the workpiece coordinate system; This represents the translation vector of the origin of the workpiece coordinate system in the global machine tool coordinate system; This represents the transformation matrix of the second rotating component; This represents the preset structural constant, which is the initial geometric vector from the rotation center of the rotating component to the tool tip under the initial configuration, and is related to the length of the tool and the tool holder. Step 2.2: Perform global sorting and independent numbering on the mesh nodes in the finite element mesh model established in Step 1.1 to obtain a list containing the global numbers of all nodes. Establish and generate a component-node index table, and mark the feature node areas of the bottom surface of the slider and the feature node areas of the guide rail contact surface that participate in the contact of the mating surfaces in the component-node index table. Step 2.3: For the slider bottom surface node area and guide rail contact surface feature node area marked in Step 2.2, based on the planar or cylindrical surface geometric features of the component, and based on the user-specified feature nodes, the spatial distance tolerance method is used for automatic identification processing to obtain the planar contact node subset and the cylindrical surface contact node subset. The planar contact node subset and the cylindrical contact node subset are merged and deduplicated to establish and generate a joint surface node index table; Step 2.4: Obtain ball joint connection pair information describing the fixed connection relationship between various components inside the machine tool, and establish a ball joint connection pair list; Step 2.5: Obtain the fixed boundary conditions connecting the bottom surface of the machine tool bed to the external foundation, the set of degrees of freedom subject to external fixed constraints, and establish the external fixed constraint table.

4. The rapid finite element modeling method for five-axis machining centers based on G-code instruction drive according to claim 3, characterized in that, In step 2: In step 2.1, the general formulas for the third-order fundamental rotation matrices of the A-axis, B-axis, and C-axis are as follows: The third-order fundamental rotation matrix of the A-axis about the X-axis The formula is: The third-order fundamental rotation matrix of the B-axis about the Y-axis The formula is: The third-order fundamental rotation matrix of the C-axis about the Z-axis The formula is: In the formula, Indicates the rotation angle of axis A; Indicates the rotation angle of the B-axis; Indicates the rotation angle of the C-axis; In step 2.3, the user-specified feature nodes include vertex nodes used to define the spatial boundary of the planar contact area of ​​the component, and circumferential nodes used to fit the geometric features of the cylindrical contact surface of the component; the mating surface node index is obtained using the feature nodes, specifically as follows: Planar node extraction: Based on three non-collinear specified feature nodes, fit the three-dimensional spatial plane equation, traverse all nodes on the corresponding component one by one, calculate the spatial vertical distance from the node to the plane, and classify the nodes whose spatial vertical distance is not greater than the preset surface tolerance into the planar contact node subset. The range of the preset surface tolerance is taken from the average edge length of the finite element mesh in the contact area of ​​the component. to ; Cylindrical surface node extraction: Based on the cylinder axis direction, reference points on the axis, and cylinder radius, the radial distance from the node to the cylinder axis is calculated. Nodes whose absolute value of the difference between the radial distance and the cylinder radius is not greater than a preset circular tolerance are included in the cylindrical surface contact node subset; the value range of the preset circular tolerance is also taken from the average edge length of the finite element mesh in the component contact area. to ; Merging and deduplication: Merge the subset of planar contact nodes and the subset of cylindrical contact nodes in different mating areas on the same component surface, remove duplicate node numbers, and determine the mating surface node index table.

5. The rapid finite element modeling method for five-axis machining centers based on G-code instruction drive according to claim 4, characterized in that, Step 3 specifically includes: Step 3.1: Parse the G-code instructions line by line; using the identified function words as a reference, extract the coordinate address characters representing the spatial geometric position and rotation axis attitude in the current line; convert the extracted characters into real-number numerical parameters to obtain the target position vector of the tool tip in the workpiece coordinate system. and the target rotation angle values ​​of each rotation axis , , ; Step 3.2, calculate the target position vector of the tool tip in the workpiece coordinate system obtained in Step 3.

1. The target rotation angle values ​​of each rotation axis are input into the machine tool kinematic model established in step 2.1 for inverse kinematic calculation, and the indexed values ​​are determined. Spatial translation vector of machine tool components and rotation matrix ; If the machine tool's G-code instructions are in a machining state with the tool tip tracking control function enabled, the machine tool kinematic model established in step 2.1 is invoked, and its inverse kinematics calculation formula is as follows: In the formula, This represents the position vector of the rotation center of the rotating component relative to the global machine tool coordinate system under the current configuration; This represents the transformation matrix of the first rotating component; This represents the target position vector of the tool tip in the workpiece coordinate system; This represents the translation vector of the origin of the workpiece coordinate system in the global machine tool coordinate system; This represents the transformation matrix of the second rotating component; This represents a preset structural constant; Using the calculated position vector The spatial translation vector of the rotating component is obtained by calculating the difference between the initial position vector of the rotating component's rotation center and the initial position vector of the rotating component under the initial configuration. : In the formula, Indicates having an index number The spatial translation vector of the machine tool component in the global machine tool coordinate system; Indicates the index number of the machine tool component; for pure translation components, the stroke displacement of each axis relative to its initial configuration obtained by the inverse kinematics operation of the current G-code instruction is used as the spatial translation vector of the corresponding pure translation component; Step 3.3, analyze the target rotation angle values ​​of each rotation axis obtained in Step 3.

1. , , Substitute the values ​​into the third-order basic rotation matrix library constructed in step 2.1 respectively. In the process, the basic rotation matrix of each rotation axis under the current configuration is calculated; Based on the series sequence of the machine tool's shafts, the third-order basic rotation matrices of each rotating axis under the current configuration are multiplied using matrix multiplication to calculate the result with index numbers. Rotation matrix corresponding to machine tool components ; Step 3.4: Update the node coordinates in the finite element mesh models of each component from Step 1.1 using the following rigid body space coordinate transformation formula: For a purely translational component, the nodal coordinate update formula is as follows: For rotating parts: In the formula, Indicates having an index number Finite element node numbering inside machine tool components; Indicates having an index number The internal parts of the machine tool The global machine tool coordinate vector of each finite element node in the current configuration; Indicates having an index number The internal parts of the machine tool The global machine tool coordinate vector of each finite element node under the initial configuration; Indicates having an index number The spatial translation vector of the machine tool component in the global machine tool coordinate system; Indicates having an index number The rotation matrix corresponding to the machine tool component; This represents the initial position vector of the rotation center of the rotating component in the global machine tool coordinate system under the initial configuration; Using the rigid body space coordinate transformation formula, the global machine tool coordinate vector of all mesh nodes of each component in the whole machine finite element model under the current configuration is calculated and obtained. Each component includes sliders and guide rails. The current configuration of the machine tool is determined by the updated global machine tool coordinate vector of the mesh nodes under the current configuration.

6. The rapid finite element modeling method for a five-axis machining center based on G-code instruction drive according to claim 5, characterized in that, Step 4 specifically includes: Step 4.1: For the slider-guide rail pair, extract the set of slider nodes involved in the contact based on the mating surface node index table established in Step 2.

3. With guide rail node set ; Step 4.2, denote the global machine tool coordinate vector of the corresponding slider node obtained in step 3.4 as... The global machine tool coordinate vector of the corresponding guide rail node is denoted as ; Step 4.3: Assign the global machine coordinate vectors of the slider node and guide rail node under the current configuration. and As the retrieval benchmark, the K-nearest neighbor method with a search count of 1 is used for spatial search and distance calculation to obtain the slider nodes. With guide rail node Spatial Euclidean distance between : In the formula, Represents slider node With guide rail node The spatial Euclidean distance between them; Represents slider node The global machine tool coordinate vector under the current configuration; Indicates guide rail node Global machine tool coordinate vector in the current configuration; superscript This represents the matrix transpose operation; The calculated spatial Euclidean distance The comparison is performed with a preset interface threshold. If the Euclidean distance in space If the distance between the nodes is not greater than a preset mating surface threshold, then the two nodes are determined to form a mating surface connection, and the node pairing table Pairs is added accordingly. The preset mating surface threshold is set based on the nominal assembly gap between the slider and the guide rail and the average size of the finite element mesh of the contact surface, and its value range is [range missing]. to If the same guide rail node is matched repeatedly by multiple slider nodes, only a unique correspondence will be retained according to the principle of minimum spatial Euclidean distance. Step 4.4: Collect all corresponding relationships that satisfy the connection of the mating surfaces and establish the node pairing table Pairs.

7. The rapid finite element modeling method for five-axis machining centers based on G-code instruction drive according to claim 6, characterized in that, Step 5 specifically includes: Step 5.1: Based on the list of ball joint connections in Step 2.4, perform uniqueness verification and deduplication on the master node degrees of freedom and slave node degrees of freedom of the ball joint connections; divide the total degrees of freedom of the whole machine finite element model into the set of eliminated slave node degrees of freedom and the set of retained degrees of freedom composed of master node degrees of freedom and unconstrained degrees of freedom. Step 5.2: For the deduplicated linear constraint equations of the ball joint connection, each slave node degree of freedom is written as an algebraic linear function of the corresponding master node degree of freedom, establishing a multi-point constraint equation system whose block algebraic form satisfies: Based on the algebraic coefficients of the multi-point constraint equation system, a master-slave constraint transformation matrix for the reduction of the degree of freedom of the ball joint connection constraint is established by expanding and combining the equations. Its explicit analytical expression is: In the formula, This represents the constraint coefficient matrix corresponding to the degrees of freedom of each node; This represents the vector of slave node degrees of freedom, which consists of the slave node degrees of freedom. This represents the constraint coefficient matrix corresponding to the retained degrees of freedom; This represents the vector of retained degrees of freedom, consisting of the master node degrees of freedom and the unconstrained degrees of freedom. Represents the master-slave constraint transformation matrix; Represents the identity matrix of the same order as the number of degrees of freedom retained.

8. The rapid finite element modeling method for a five-axis machining center based on G-code instruction drive according to claim 7, characterized in that, Specifically, step 6 includes: Step 6.1, targeting the moving parts inside the machine tool Obtain its rotation matrix from step 3.

3. If it has an index number If the machine tool component is a rotating component, then its third-order rotation matrix is ​​used. As basic unit items, they are laid out flat along the main diagonal without overlap. This establishes a size scale of block diagonal rotation transformation matrix ,in Indicates having an index number The total number of nodes in the machine tool components; if it has an index number. If a machine tool component is a purely translational component, then its rotation matrix... It is a third-order identity matrix, and its corresponding block diagonal rotation transformation matrix. for An identity matrix of order 1; Step 6.2: Extract components from the unit-matrix index table established in Step 1.

3. The element stiffness matrix and element mass matrix under the initial configuration are assembled into the initial component stiffness matrix. and the initial component mass matrix For rotating components, calculate the global stiffness matrix of the current configuration after rotational transformation. With the current configuration global mass matrix Specifically: Step 6.2.1: When a machine tool component undergoes rigid body rotation, its inherent elastic characteristics remain unchanged relative to the initial configuration; introduce the global displacement vector under the current configuration. and the nodal displacement vectors in the component's local coordinate system It satisfies the following strain energy invariance formula: The initial displacement vector is expressed as Substituting into the strain energy invariance formula, we obtain the transformation formula: Based on this, the rotational transformation formula for the stiffness matrix is ​​derived as follows: Step 6.2.2: Introduce the global velocity vector under the current configuration. and the nodal velocity vectors in the component's local coordinate system It satisfies the following kinetic energy invariance formula: Since the initial velocity vector is equivalent to the product of the transpose of the block diagonal rotation transformation matrix and the global velocity vector, it satisfies the following relationship: Substituting this into the kinetic energy invariance formula, the rotation transformation formula for the mass matrix is ​​derived as follows: In the formula, Indicates having an index number The initial component stiffness matrix of the machine tool component under the initial configuration; Indicates having an index number The initial component mass matrix of the machine tool components under the initial configuration; Indicates having an index number The dimensions of the machine tool components are The block diagonal rotation transformation matrix; Indicates having an index number The global stiffness matrix of the current configuration of the machine tool components after rotational transformation; Indicates having an index number The global mass matrix of the current configuration after the rotational transformation of the machine tool components; Step 6.3: Obtain the global number of each mesh node using the component-node index table established in Step 2.2; based on the matrix assembly principle of the finite element direct stiffness method, perform global addressing and accumulation processing of the stiffness matrix and mass matrix: For the overall unconstrained stiffness matrix of the whole machine Based on the global number of the mesh nodes, the global stiffness matrix of each component's current configuration is determined. The matrix elements in the matrix are accumulated into the overall unconstrained stiffness matrix of the entire machine. In the process, the unconstrained overall stiffness matrix of the entire machine is formed through assembly. Its formula is: For the unconstrained overall mass matrix of the whole machine Based on the global number of the grid nodes, the global mass matrix of each component's current configuration is calculated. The matrix elements in the matrix are accumulated into the unconstrained overall mass matrix of the whole machine. In the process, the unconstrained overall mass matrix of the entire machine is formed through assembly. Its formula is: In the formula, Represents the overall unconstrained stiffness matrix of the entire machine. The Middle line, number Matrix elements at column positions; Represents the unconstrained overall mass matrix of the entire machine. The Middle line, number Matrix elements at column positions; Indicates the total number of independent moving parts contained in the machine tool; indicates the index number of the machine tool parts; Indicates having an index number The machine tool component in the global stiffness matrix under the current configuration is the first... line, number Matrix elements at column positions; Indicates having an index number The machine tool component in the global mass matrix under the current configuration is the [number]th [unit]. line, number Matrix elements at column positions; and These represent the global row number and global column number of the overall matrix, respectively; and These represent the internal row number and internal column number of the current configuration global stiffness matrix or the current configuration global mass matrix of the corresponding component, respectively; Implement the mapping from the global number of grid nodes to the row and column positions of the overall matrix; Step 6.4: Based on the matrix order of the unconstrained overall stiffness matrix and the unconstrained overall mass matrix of the whole machine assembled in Step 6.3, establish an unconstrained overall damping matrix of the same order, and initialize the unconstrained overall damping matrix of the whole machine to a zero matrix. The initialization expression for the unconstrained overall damping matrix of the entire machine is: In the formula, This represents the unconstrained overall damping matrix of the entire machine after initialization. This represents the matrix order of the unconstrained overall stiffness matrix and the unconstrained overall mass matrix of the entire machine; express A zero matrix of order 1.

9. The rapid finite element modeling method for a five-axis machining center based on G-code instruction drive according to claim 8, characterized in that, Specifically, step 7 is as follows: Step 7.1: For the mating surface connections of the Pairs node pairing table established in Step 4.4, introduce spatial spring elements simulating contact stiffness characteristics and spatial damping elements simulating energy dissipation characteristics to establish the corresponding spring element stiffness matrix. With the damping matrix of the damping unit ;in, The damping matrix of the damping unit is used to accumulate in the overall unconstrained damping matrix of the whole machine established in step 6.4: In the formula, This represents the stiffness matrix of the spring element established by the translational stiffness coefficient of the mating surface; , , , These are all stiffness sub-matrices used to describe the stiffness characteristics of the physical connection between paired nodes; specifically, equivalent translational stiffness coefficients are introduced from the mating surface in the three spatial translational directions X, Y, and Z of the global machine tool coordinate system. , , This corresponds to the Pairs slider node in the node pairing table. The third-order self-stiffness submatrix With the corresponding guide rail node The third-order self-stiffness submatrix and slider nodes With guide rail node The coupling stiffness submatrix between and Its matrix representation is as follows: In the formula, This represents the damping matrix of the damping element established by the translational damping coefficient of the interface; , , , These are all damping sub-matrices used to describe the energy dissipation damping characteristics between paired nodes; specifically, equivalent translational damping coefficients of the joint surface in the three spatial translational directions X, Y, and Z of the global machine tool coordinate system are introduced. , , This corresponds to the slider node in the node pairing table Pairs. The third-order self-damped matrix With the corresponding guide rail node The third-order self-damped matrix and the slider node With guide rail node The coupling damping sub-matrix between and Its matrix representation is as follows: Among them, the equivalent translational stiffness coefficient preset at the interface , , With equivalent translational damping coefficient , , It can be obtained by conducting experimental modal tests and frequency response function parameter identification on the machine tool mating surface, or by analytical calculation based on contact mechanics theory combined with surface micro-morphology parameters, or by directly setting the rated contact stiffness and damping manual parameters provided by suppliers of moving parts such as linear guides. Based on the slider nodes in the Pairs node pairing table With guide rail node The degree of freedom numbering, and the stiffness matrix of the spring element. With the damping matrix of the damping unit The sub-block matrices in the matrix are summed to the overall unconstrained stiffness matrix of the whole machine. With respect to the unconstrained overall damping matrix of the whole machine middle; Suppose there are slider nodes in the node pairing table Pairs. The global number is Guide rail node The global number is Since the finite element mesh model uses three-dimensional solid elements with 3 translational degrees of freedom per node, the slider node... The set of row and column indices corresponding to the overall matrix is: Guide rail node The set of row and column indices corresponding to the overall matrix is: Sub-block addressing and accumulation of matrix elements are performed based on the index set, thereby obtaining the overall unconstrained stiffness matrix of the entire machine. With respect to the unconstrained overall damping matrix of the whole machine The application of flexible connection at the interface is completed in the middle; Step 7.2, using the master-slave constraint transformation matrix calculated in step 5.

2. The overall stiffness matrix of the unconstrained whole machine Unconstrained overall mass matrix of the whole machine and the overall unconstrained damping matrix of the whole machine By performing a condensation transformation to eliminate the slave node degrees of freedom constrained by the ball joint connection, a fixed connection of the internal components of the machine tool is achieved. Its condensation analytical expression is: In the formula, This represents the overall stiffness matrix after the constraint condensation of the ball joint connection; This represents the overall mass matrix after the constraint condensation of the ball joint connection; This represents the overall damping matrix after the ball joint connection is constrained and condensed; This represents the master-slave constraint transformation matrix used for the condensation of the degree of freedom of the ball joint connection constraint; This represents the overall unconstrained stiffness matrix of the entire machine. This represents the unconstrained overall mass matrix of the entire machine. This represents the overall unconstrained damping matrix of the entire machine; Step 7.3: Determine the new row and column indices based on the external fixed constraint table established in Step 2.5, and use the elimination method to refine the condensed global stiffness matrix. Overall quality matrix and the overall damping matrix Boundary condition processing is performed to fix the bottom surface of the machine tool bed to the external foundation, eliminating the degrees of freedom subject to external constraints. This reduces the order of the system matrix and eliminates the matrix singularity caused by rigid body displacements, resulting in the final global stiffness matrix after eliminating the influence of rigid body displacements. Final overall quality matrix and the final overall damping matrix .

10. The rapid finite element modeling method for a five-axis machining center based on G-code instruction drive according to claim 9, characterized in that, Specifically, step 8 includes: The result calculated in step 7.3 , , Import the data into the solver for dynamic calculations; according to the G-code instruction sequence, execute steps 3 to 8 in a loop to obtain the evolution law of the dynamic characteristics of the entire machine tool along the machining trajectory.

Images