Motorized spindle thermal characteristic analysis modeling method based on multi-patch isogeometric analysis

Abstract

The invention discloses a motorized spindle thermal characteristic analysis modeling method based on multi-patch isogeometric analysis, which comprises the following steps: S10, carrying out geometric modeling of a high-speed motorized spindle: constructing a multi-patch solid model comprising a motorized spindle stator, a rotor, a bearing and a cooling pipeline of the high-speed motorized spindle by adopting a three-dimensional NURBS patch splicing technology, geometric continuity and parameterization consistency are ensured; s20, performing heat conduction analysis: constructing a transient heat conduction control equation based on a Fourier heat conduction differential equation in combination with high-order smoothness characteristics of the IGA, and solving a temperature evolution process in a time domain through a backward difference method; and S30, performing thermal error modeling: introducing an elastic mechanical constitutive relation and a fixed support boundary, and establishing a thermal-force coupling solving framework.

Description

Technical Field

[0001] This invention belongs to the field of precision CNC machine tool technology, specifically relating to a method for analyzing and modeling the transient thermal characteristics of a high-speed electric spindle system based on geometric analysis such as multi-faceted plates. Background Technology

[0002] As a core functional component of precision CNC machine tools, the thermal characteristics of high-speed electric spindles directly affect machining accuracy and equipment reliability. Under high-speed rotation conditions, the uneven distribution of the transient temperature field generated by heat sources such as motor losses and bearing friction in the electric spindle system leads to non-uniform thermal expansion and thermal stress accumulation in key components, causing thermally induced errors such as axial elongation and radial runout of the spindle. Therefore, how to efficiently and accurately analyze and predict the temperature field and thermal deformation behavior of high-speed electric spindle systems has become one of the current research hotspots in the field of precision manufacturing.

[0003] In recent years, various models for predicting spindle thermal error have been proposed in existing technologies. These models can be broadly classified into two categories: data-driven models and mechanistic analysis models. Data-driven modeling is a black-box model, such as multivariate regression analysis, neural networks, support vector machines, and Bayesian networks. For example, an encoder and multi-head attention mechanism were used to explore key information in thermal images, and a spindle thermal error model based on a visual transformer was proposed. However, the generalization ability of data-driven models is limited by the quality of training data, and it is difficult to reveal the physical mechanisms of heat conduction.

[0004] Research on the thermal characteristic modeling of electric spindles based on thermal mechanism analysis can elucidate the logic of heat generation, transfer, and conversion mechanisms within the spindle. Among the various methods employed, the Finite Element Method (FEA) is the most popular and widespread method for temperature field mechanism modeling due to its high simulation accuracy and ease of modeling. For example, using the optimized heat transfer coefficient of a spiral cooling system as the boundary condition of the finite element model, a predictive model of the electric spindle temperature field is established. However, EFA modeling and solving typically rely on mesh reconstruction, which can easily lead to geometric distortion and data inconsistencies during the conversion from CAD to CAE models. The thermal network method based on heat transfer theory is also a classic method for simulating heat transfer in electric spindles. While its simplified model is computationally efficient, it relies excessively on empirical parameters and struggles to characterize complex three-dimensional heat conduction paths. Analytical methods based on physics fields, when comprehensively considering multiple influencing factors of the electric spindle temperature field, result in highly complex physics field governing equations. Summary of the Invention

[0005] In view of the above-mentioned problems, the present invention provides a transient thermal characteristic analysis and modeling method for high-speed electric spindle systems based on geometric analysis such as multi-faceted plates, so as to accurately characterize complex heat conduction paths and predict thermally induced errors.

[0006] To solve the above-mentioned technical problems, the present invention adopts the following technical solution:

[0007] A method for transient thermal characteristic analysis and modeling of a high-speed electric spindle system based on geometric analysis such as multi-faceted patches includes the following steps:

[0008] S10, Perform geometric modeling of the high-speed electric spindle: Use 3D NURBS patch splicing technology to construct a multi-patch solid model of the high-speed electric spindle, including the spindle stator, rotor, bearings, and cooling pipes, to ensure geometric continuity and parametric consistency;

[0009] S20, heat conduction analysis: Based on the Fourier heat conduction differential equation and combined with the high-order smoothness characteristics of IGA, a transient heat conduction control equation is constructed, and the temperature evolution process in the time domain is solved by the backward difference method.

[0010] S30, Perform thermal error modeling: Introduce the constitutive relation of elasticity and the fixed support boundary to establish a thermo-mechanical coupling solution framework.

[0011] In one possible implementation, S10 includes: using multi-faceted modeling technology to transform the entire electric spindle system into a structure assembled from multiple 3D non-uniform rational B-spline NURBS facets, each facet having a different thermal effect.

[0012] The following formula applies when stitching together individual 3D NURBS patches:

[0013]

[0014] Where P i a and Let A be the control points, belonging to facet a and facet b respectively, and N be the set of contact surfaces of adjacent facets. Then (a, b) represents that facet a and facet b are adjacent.

[0015] During the establishment of the integration domain, the parameter space of the NURBS patch is defined, and the existing overlapping or gap regions are corrected. The control points of each patch are globally numbered, and a mapping rule from local numbering to global numbering is adopted.

[0016] In one possible implementation, S10 further includes: isogeometric spatial mapping relations under a three-dimensional NURBS entity; constructing a three-node index space from a given node vector; and constructing a corresponding parameter space based on the units formed by the non-zero node intervals in the index space, thereby calculating the NURBS basis functions; the units are defined in the parameter space, and the Gaussian integral domain is defined in the parent space. Therefore, the entity in the spatial discretization process of IGA includes two direct mappings. In the first mapping, the mapping originates from the parent space supporting Gaussian integrals. Transform to parameter space unit Through the Jacobian matrix Completed; in the second mapping, from the parameter space unit To physical space unit Ω e The transformation is achieved through the Jacobian matrix. Completed, among which Composed of the gradients of the basis functions, the Jacobian matrix J corresponding to the indirect mapping from the parent space to the physical space in 3D NURBS entities is obtained from these two direct mapping relationships. v In this way, PDEs in the physical model are transformed from strong form to weak form.

[0017] In one possible implementation, S20 includes:

[0018] Fourier's law is the foundation of heat conduction. Based on the three-dimensional heat transfer domain Ω formed by the boundary Γ, its transient heat conduction governing differential equation is as follows:

[0019]

[0020] In the formula, T is the unknown temperature field function, Q is the intensity of the internal heat source, ρ is the material density, c is the specific heat capacity of the material, and λ is the specific heat capacity of the material. x , λ y , λ z These are the thermal conductivity of the material in the x, y, and z directions, respectively;

[0021] The initial temperature conditions, isothermal boundary conditions, isothermal boundary conditions, and convective heat transfer boundary conditions of the transient heat transfer system are as follows:

[0022]

[0023] In the formula, T0 is the initial temperature of the given region, and T λ Here, q is the constant temperature of a given region, and h is the heat flux density of that region. c Γ is the convective heat transfer coefficient for a given region, and n is the unit outward normal vector on the boundary Γ.

[0024] The temperature field function T at any point within a NURBS solid element is represented by the temperature value T at the control point. i,j,k and three-variable NURBS basis functions The interpolation form is as follows:

[0025]

[0026] In the formula, n1 = p b +1,n2=q b +1,n3=r b +1, n n=n1·n2·n3 represents the total number of control points or basis functions associated with this NURBS entity cell;

[0027] Transforming equation (12) into the standard Galerkin weak form and combining it with equation (14) yields the discrete system equations based on isogeometric analysis, in the following form:

[0028]

[0029] In the formula, N T K T , {P} and {P} are respectively composed of NURBS solid units Ω e It is assembled from correlation matrices or vectors, and its unit correlation matrix or vector is described as follows: For the element heat capacity matrix, For the unit heat conduction matrix, P is the element heat conduction matrix related to the convective heat transfer boundary conditions. e Let be the element temperature load vectors. They can be written in Gaussian integral form over the standard integration domain in the parent space, as follows:

[0030]

[0031] In the formula, To be compatible with unit Ω e The relevant basis function matrix, NURBS unit Ω e The temperature gradient interpolation matrix within the element Ω, and related to the element Ω. e The temperature gradient interpolation matrix B corresponding to the k-th control point Tk (k = 1, 2, ..., n) n The format is as follows:

[0032] Then, based on the chain rule and the Jacobian matrix... Further results were obtained:

[0033]

[0034] In addition to employing isogeometric discretization in the spatial domain, the transient heat conduction system based on isogeometric analysis also uses the backward difference method in the time domain. The system discretization equation (15) for transient heat conduction is further written as:

[0035]

[0036] In the formula, Δt is the time step.

[0037] In one possible implementation, S20 includes: in numerical calculation, employing h-refinement and p-refinement strategies, where h-refinement increases spatial resolution by increasing the density of control points, and p-refinement reduces numerical error by increasing the order of basis functions, and solving using a global heat conduction matrix assembly formula.

[0038] In one possible implementation, S30 includes: defining convection boundary conditions and heat sources at relevant locations on the electric spindle, including forced convection between the shaft and air, the natural convection coefficient of air (marked with dark green lines), the convection region between the stator cooling jacket and coolant, the convection region for cooling the bearings, forced convection between the stator and rotor, the heating regions of the stator and rotor, and the heating regions of the front and rear bearings; estimating the natural convection heat transfer coefficient of ambient air to obtain the heat dissipation caused by local convection under different thermal conditions.

[0039] The present invention has the following beneficial effects:

[0040] (1) In terms of geometric modeling, by introducing multi-faceted NURBS surface structure and global node numbering strategy, the geometric accuracy modeling and parameter consistency control of key components of the spindle are realized, laying the foundation for the simulation of heat conduction path continuity.

[0041] (2) In terms of thermal field modeling, a three-dimensional transient heat conduction geometric discrete model was established. Combined with the backward time difference method and the local-global coupling assembly method, it can stably solve the temperature evolution process under the action of complex non-uniform heat sources in the spindle system, avoiding the temperature gradient distortion problem caused by mesh distortion in the traditional FEA method.

[0042] (3) In terms of thermal error modeling, the thermo-mechanical coupling mechanism is integrated, and the advantages of high-order continuity in IGA are combined to capture the influence of small thermal deformations on the overall structural accuracy in the parameter space. Simulation results show that the maximum deviation between the simulated thermal expansion error in the Z-axis direction and the experimental data is only 1.1 μm, which is significantly better than the FEA method (deviation 3.9 μm).

[0043] (4) By comparing the experimental temperature measurement data with the simulation results of IGA and FEA models, the accuracy and robustness of the proposed method in predicting temperature field and thermal deformation are verified. It is especially suitable for the transient thermal performance simulation of electric spindle systems with complex structure, multiple heat sources and multiple material distributions. Attached Figure Description

[0044] Figure 1 This is a flowchart illustrating the steps of a transient thermal characteristic analysis and modeling method for a high-speed electric spindle system based on geometric analysis such as multi-faceted patches, according to an embodiment of the present invention.

[0045] Figure 2 This is a schematic diagram of the basic structure of the electric spindle;

[0046] Figure 3 This is a schematic diagram illustrating the relationship between B-spline curves and NURBS curves.

[0047] Figure 4 This relates to the isogeometric model transformation and spatial mapping relationship based on NURBS;

[0048] Figure 5 This is a schematic diagram illustrating a multi-faceted model.

[0049] Figure 6 This is a schematic diagram illustrating the segmentation and refinement of solid surfaces.

[0050] Figure 7 This is a schematic diagram showing the local and global numbering of a double-sided patch model;

[0051] Figure 8 A simplified schematic diagram of the geometric model of the multifaceted electric spindle.

[0052] Figure 9 A schematic diagram of the multifaceted assembly of the geometric model of the electric spindle;

[0053] Figure 10 Schematic diagram of heat source and heat dissipation path analysis inside the electric spindle;

[0054] Figure 11 Schematic diagram of finite element analysis and isogeometric analysis of the transient temperature field of an electric spindle;

[0055] Figure 12 The temperature convergence characteristic curves are shown for different mesh densities.

[0056] Figure 13 This is a schematic diagram of the finite element analysis and isogeometric analysis of the axial thermal deformation field of the electric spindle.

[0057] Figure 14 This is a schematic diagram of the temperature and displacement measurement device for a high-speed electric spindle system. Detailed Implementation

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

[0059] Reference Figure 1 The image shows an embodiment of the present invention of a method for transient thermal characteristic analysis and modeling of a high-speed electric spindle system based on geometric analysis such as multi-faceted patches, comprising the following steps:

[0060] S10, Perform geometric modeling of the high-speed electric spindle: Use 3D NURBS patch splicing technology to construct a multi-patch solid model of the high-speed electric spindle, including the spindle stator, rotor, bearings, and cooling pipes, to ensure geometric continuity and parametric consistency;

[0061] like Figure 2 The diagram shows the internal structure of the electric spindle of a three-axis precision CNC machine tool in a specific application example. The physical model of this electric spindle, designed for a maximum speed of 20,000 revolutions per minute (RPM) on the experimental three-axis machine tool, includes a central spindle, motor stator, motor rotor, front and rear bearings, and front and rear bearing housings. During high-speed operation, the electric spindle generates a significant amount of heat due to power losses in the stator and rotor, as well as frictional heat generation from the front and rear bearings. Therefore, it is essential to circulate coolant through the cooling pipes in the stator cooling jacket and the cooling pipes in the front and rear bearing housings. Simultaneously, the airflow around the spindle continuously carries away the heat.

[0062] S20, heat conduction analysis: Based on the Fourier heat conduction differential equation and combined with the high-order smoothness characteristics of IGA, a transient heat conduction control equation is constructed, and the temperature evolution process in the time domain is solved by the backward difference method.

[0063] S30, Perform thermal error modeling: Introduce the constitutive relation of elasticity and the fixed support boundary to establish a thermo-mechanical coupling solution framework.

[0064] In a specific application example, an embodiment of the present invention provides a method for transient thermal characteristic analysis and modeling of a high-speed electric spindle system based on geometric analysis such as multi-faceted patches. S10 includes:

[0065] The entire electric spindle system is transformed into a structure composed of multiple 3D non-uniform rational B-spline NURBS facets using multi-facet modeling technology. Each facet has a different thermal effect.

[0066] The following formula applies when stitching together individual 3D NURBS patches:

[0067]

[0068] Where P i a and Let A be the control points, belonging to facet a and facet b respectively, and N be the set of contact surfaces of adjacent facets. Then (a, b) represents that facet a and facet b are adjacent.

[0069] During the establishment of the integration domain, the parameter space of the NURBS patch is defined, and the existing overlapping or gap regions are corrected. The control points of each patch are globally numbered, and a mapping rule from local numbering to global numbering is adopted.

[0070] S10 further includes: isogeometric spatial mapping relations under 3D NURBS entities, constructing a three-node index space from given node vectors, and constructing the corresponding parameter space based on the units formed by non-zero node intervals in the index space, thereby calculating NURBS basis functions; the units are defined in the parameter space, and the Gaussian integration domain is defined in the parent space. Then, the entity contains two direct mappings in the spatial discretization process of IGA. In the first mapping, from the parent space supporting Gaussian integration... Transform to parameter space unit Through the Jacobian matrix Completed; in the second mapping, from the parameter space unit To physical space unit Ω e The transformation is achieved through the Jacobian matrix. Completed, among which Composed of the gradients of the basis functions, the Jacobian matrix J corresponding to the indirect mapping from the parent space to the physical space in 3D NURBS entities is obtained from these two direct mapping relationships. v In this way, PDEs in the physical model are transformed from strong form to weak form.

[0071] The following further explains the implementation process of S10. Non-uniform rational B-splines (NURBS) have advantages such as a unified mathematical model and flexible shape control in the geometric definition of industrial products. For example... Figure 3 As shown, taking NURBS as an example, we first need to base our work on the node vector ξ={ξ1,ξ2,...,ξ...} m+p+1 Construct a B-spline curve where the elements of ξ are monotonically increasing in the parameter space. The basis functions of the B-spline can be constructed using the following formula:

[0072]

[0073] Parameter description: i = 1, 2, ..., m + p b +1, ξ i It is the i-th element in the node vector, m is the number of basis functions, and p b It is the highest order of the basis functions, connecting the basis functions with the control point P. i b Combining these, the functional form of the B-spline curve is as follows:

[0074]

[0075] It's important to note that a node vector containing equally spaced elements is called a uniform node vector; otherwise, it's called a non-uniform node vector. Furthermore, if the first and last elements of a node vector have p... b If a node vector is repeated +1 times, it is called an open node vector.

[0076] Non-uniform rational B-spline (NURBS) curves are obtained by projecting a higher-dimensional B-spline curve. The expression for the NURBS basis functions is:

[0077]

[0078] Where i = 1, 2, ..., m, ω i are basis functions The weighting coefficients. The basis functions and control points P. i b In summary, the functional form of the NURBS curve is as follows:

[0079]

[0080] To meet the requirements of solving two-dimensional problems, the one-dimensional NURBS curve is extended into a two-dimensional NURBS surface, which can be done using tensor products. The basis functions of the NURBS surface can be expressed as:

[0081]

[0082] Where i = 1, 2, ..., n, j = 1, 2, ..., m, and It has p in the ξ direction and η direction respectively. b and q b Basis functions of order. Connecting basis functions to control points. In combination, the functional form of the NURBS surface is:

[0083] Similarly, three-variable NURBS basis functions It can be represented as:

[0084] Therefore, the NURBS entity function R(ξ,η,β) can be expressed as:

[0085]

[0086] For most physical problems with complex geometries, the partial differential equations (PDEs) used to describe the physical model are difficult to solve analytically directly. Typically, the PDEs in the physical model are transformed from strong to weak forms, and then the physical model is discretized using a large number of minimum elements to obtain a mathematical analysis model similar to the PDEs of the physical model. The numerical solution of the mathematical model is approximately equal to the analytical solution under the same conditions. This is a common idea in isogeometric analysis (IGA) and traditional finite element analysis (FEA).

[0087] Figure 4 This paper demonstrates the NURBS-based isogeometric model transformation, aiming to illustrate that the core idea of ​​the IGA method is to unify CAD and CAE. Specifically, it adopts NURBS basis functions with precise geometric representation as shape functions in numerical computation, and these NURBS basis functions are defined based on the node vectors given in the parameter space. Furthermore, due to the non-standard nature of the elements in the parameter space, these elements need to be further mapped to a standard Gaussian integration domain for integration operations. Therefore, numerical integration in isogeometric analysis involves mappings between the parameter space, physical space, and parent space. Figure 4 The isogeometric spatial mapping relationship under a certain 3D NURBS entity is given. A three-node index space is constructed from the given node vectors Ξ1={0,0,0,0.25,0.25,0.5,0.5,0.75,0.75,1,1,1}, Ξ2={0,0,0.25,0.5,0.75,1,1}, and Ξ3={0,0,1,1}. The corresponding parameter space is then constructed based on the "cells" formed by the non-zero node intervals in the index space, thereby calculating the NURBS basis functions. Defined within the parameter space, the Gaussian integration domain within the parent space is defined as follows: This entity involves two direct mappings during the spatial discretization process of the IGA. In the first mapping, the mapping originates from the parent space that supports Gaussian integrals. Transform to parameter space unit Through the Jacobian matrix Completed can be represented as follows:

[0088]

[0089] In the second mapping, from the parameter space unit To physical space unit Ω e The transformation is achieved through the Jacobian matrix. Completed, among which Composed of the gradients of the basis functions, it can be expressed in the following form:

[0090]

[0091] From these two direct mapping relationships, we can obtain the Jacobian matrix J corresponding to the indirect mapping from the parent space to the physical space in 3D NURBS entities. v According to the chain rule, it can be expressed as:

[0092]

[0093] In this way, the PDEs in the physical model are transformed from strong form to weak form, making them easier to solve.

[0094] Furthermore, S20 includes:

[0095] Fourier's law is the foundation of heat conduction. Based on the three-dimensional heat transfer domain Ω formed by the boundary Γ, its transient heat conduction governing differential equation is as follows:

[0096]

[0097] In the formula, T is the unknown temperature field function, Q is the intensity of the internal heat source, ρ is the material density, c is the specific heat capacity of the material, and λ is the specific heat capacity of the material. x , λ y , λ z These are the thermal conductivity of the material in the x, y, and z directions, respectively;

[0098] The initial temperature conditions, isothermal boundary conditions, isothermal boundary conditions, and convective heat transfer boundary conditions of the transient heat transfer system are as follows:

[0099]

[0100] In the formula, T0 is the initial temperature of the given region, and T λ Here, q is the constant temperature of a given region, and h is the heat flux density of that region. c Γ is the convective heat transfer coefficient for a given region, and n is the unit outward normal vector on the boundary Γ.

[0101] The temperature field function T at any point within a NURBS solid element is represented by the temperature value T at the control point. i,j,k and three-variable NURBS basis functions The interpolation form is as follows:

[0102]

[0103] In the formula, n1 = p b +1,n2=q b +1,n3=r b +1, n n =n1·n2·n3 represents the total number of control points or basis functions associated with this NURBS entity cell;

[0104] Transforming equation (12) into the standard Galerkin weak form and combining it with equation (14) yields the discrete system equations based on isogeometric analysis, in the following form:

[0105]

[0106] In the formula, N T K T , {P} and {P} are respectively composed of NURBS solid units Ω e It is assembled from correlation matrices or vectors, and its unit correlation matrix or vector is described as follows: For the element heat capacity matrix, For the unit heat conduction matrix, P is the element heat conduction matrix related to the convective heat transfer boundary conditions. e Let be the element temperature load vectors. They can be written in Gaussian integral form over the standard integration domain in the parent space, as follows:

[0107]

[0108] In the formula, To be compatible with unit Ω e The relevant basis function matrix, NURBS unit Ω e The temperature gradient interpolation matrix within the element Ω, and related to the element Ω. e The temperature gradient interpolation matrix B corresponding to the k-th control point Tk (k = 1, 2, ..., n) n The format is as follows:

[0109] Then, based on the chain rule and the Jacobian matrix... Further results were obtained:

[0110]

[0111] In addition to employing isogeometric discretization in the spatial domain, the transient heat conduction system based on isogeometric analysis also uses the backward difference method in the time domain. The system discretization equation (15) for transient heat conduction is further written as:

[0112]

[0113] In the formula, Δt is the time step.

[0114] The following provides a further explanation of the multi-faceted solid structure. Assume the existence of a solid... It consists of κ non-overlapping faces, then:

[0115]

[0116] Where Γ i,j It's a sheet of dough Ω i Kneading dough Ω j The interface between these surfaces is called the interface. For each facet in ν, it can be defined by the parameter domain [-1, 1]. 3 The tensor product NURBS entity R(ξ,η,β) defined above is obtained from formula (8).

[0117] like Figure 5 As shown, the interface between Ω1 and Ω2

[0118] Typically, each facet in a solid ν is defined independently, so it is necessary to consider whether the geometry and parameterization of the interfaces between faces match. Here, two properties of multi-facet structures are defined: (1) Geometric consistency: The interfaces between two faces are perfectly aligned, otherwise they are geometrically inconsistent. (2) Parameterization matching: The premise is that there must be consistency between faces. If the parameters of each facet are the same, then they are matched, otherwise they are parameterized mismatched.

[0119] For example, from Figure 5 As can be seen, pieces Ω1 and Ω3 are in The two points intersect, and the corresponding intersection surface is Γ. 1,3 However, this interface only captures a portion of the boundary surface of Ω3; therefore, Ω1 and Ω3 are defined as geometrically inconsistent patches. Although Ω1 and Ω2 share the boundary surface Γ... 1,2 They are consistent, but because the nodes on the interface are different, they are defined as parameterized mismatched patches.

[0120] Before analyzing the multi-faceted solid structure, a foundation must be laid by handling segmentation and refinement operations. These two operations ensure consistent physical properties are established at the interfaces of adjacent segments. Here, we will still use... Figure 5 Taking the multi-faceted model as an example, for geometrically inconsistent scenarios, we simply divide the solid facet into sub-facets along the boundary surface to obtain a geometrically consistent boundary surface for any pair of adjacent facets.

[0121] Figure 6 First, the process of dividing the Ω3 facet was demonstrated, transforming the entire model from a 3-facet model into a 4-facet model. We divided the Ω3 facet along the isoparametric lines in the XY plane (shown as the yellow gradient surface in the figure), and then aligned the dividing surfaces with the interface Γ. 1,2 On the same plane. This divides the Ω3 facet into two sub-facets. and They are at the interface The area clearly exhibits geometric consistency. Furthermore, the newly segmented patch... At the boundary surface, it exhibits geometric consistency with patch Ω1. We use sub-patterns and It replaces the original patch Ω3 and is used for the patch representation of the geometric domain ν.

[0122] Each facet in ν is independently defined, and their parameterizations are not interdependent or constrained. For a pair of non-parameterized matched faces, we refine them by performing order increases and node insertions in a parameterization direction aligned with the interface. The aim of this process is to establish uniformity in order and nodes, ensuring that the two faces share the same control points on their corresponding boundary faces. This refinement process... Figure 6 As explained in the text, the original isoparametric lines of each facet in the 4-faceted cuboid are represented by solid black lines, while newly inserted isoparametric lines are represented by red lines. After refinement, we replace the original faces in ν with the refined faces.

[0123] Through the above process, the C between the surface pieces 0 - Continuity can be easily handled by simply unifying the degrees of freedom (DOF) of a patch onto another patch, without requiring complex algorithms or additional coupling operations. For example... Figure 7 As shown, a cube model with two faces is illustrated. The red and blue dots represent the control points of the two faces. While the control points of each face are independent in their local numbering, they are integrated into a global numbering system through unified numbering, ensuring connectivity and consistency between the faces. These preprocessing steps only need to be performed once, and the relationships between the faces are established in a simple way.

[0124] When calculating the discrete equations of the global system of transient temperature field with matched patches (Equation (19)), the temperature matrix T contains the temperature control coefficients on all patches, as well as the interface Γ between adjacent patches S1 and S2. 1,2 Shared degrees of freedom on (denoted as Γ). These degrees of freedom of the face are divided into two parts: shared degrees of freedom (T). Γ The heat conduction matrix K is defined as follows: K represents the heat transfer matrix K. The heat transfer matrix ... T convective heat transfer matrix and the term of rate of change over time Taking into account:

[0125] So This represents the combined thermal conduction coupling matrix between the facets.

[0126] The heat conduction matrix K in formula (19) T and The temperature load vector P is calculated based on the local degrees of freedom of each patch, and then they are assembled into a global system, in which K in equation (21) is used. t The term is used to represent the comprehensive heat conduction matrix. Similarly, the temperature load vector for each patch is divided into two parts: P1 (for S1) and P2 (for S2), and the temperature load vector of the interface P... Γ The heat conduction matrix is ​​divided into multiple submatrices: related to shared degrees of freedom. And related to the unique degrees of freedom of each facet. and It also includes the thermal conduction coupling term between the two surfaces. and Finally, combining the external heat source terms P1, P2, and P... Γ Combined with the temperature value from the previous time step, the temperature field at the current time step is calculated. The assembly of the heat conduction matrices for the two surfaces is expressed as follows:

[0127] By solving this equation, we can obtain the temperature distribution of each patch at the current moment and calculate the heat conduction process and temperature change of the entire system.

[0128] Based on 3D NURBS model entities Figure 8 To illustrate, the control point mesh is hidden from the overall main axis model. To better observe the dimensions of the geometric model, two perspectives are added to showcase the details of the control point mesh. The red dots represent the control points. This geometric model is composed of multiple 3D NURBS model entities. Figure 8 The bottom right corner shows the control point mesh inside one of the patches. It can be seen that the mesh is arranged very regularly, which will help to capture the temperature gradient changes more accurately in later simulations.

[0129] When establishing a multi-faceted electric spindle model, the first step is to simplify the electric spindle's geometry appropriately. Based on boundary conditions, key heat source regions are identified, and features with minimal impact on heat transfer, such as cooling channel details, are ignored. Multi-faceted modeling technology is then used to transform the entire electric spindle system into a structure assembled from 49 NURBS facets. Figure 9 The assembly of the multi-faceted panels is illustrated using a half-axial section view. Each panel has a different thermal effect; the bearing area generates significant frictional heat, the motor stator and cooling sleeves are responsible for heat dissipation, the high-speed rotation of the rotor generates the main heat source, and the outer casing provides thermal insulation. The material properties of each panel are given in Table 1.

[0130] Table 1 Material properties of the faceplate

[0131]

[0132] During the assembly of a multi-faceted geometric model, it is necessary to ensure the C0 continuity between the faces to guarantee the stability of heat conduction calculations. Therefore, the following formula must be satisfied when assembling the various 3D NURBS faces:

[0133]

[0134] Where P i a and Let A be the control points, belonging to facet a and facet b respectively, and let N be the set of contact surfaces of adjacent facets. Then (a, b) represents that facet a and facet b are adjacent.

[0135] During the establishment of the integration domain, it is necessary to define the parameter space of the NURBS patches and correct for any overlapping or gap regions. To ensure the correct assembly of the heat conduction matrix, the control points of each patch need to be globally numbered using a mapping rule from local to global numbers.

[0136] In numerical computation, to improve computational accuracy and ensure mesh independence, S20 employs h-refinement and p-refinement strategies. h-refinement increases spatial resolution by increasing the density of control points, while p-refinement reduces numerical errors by increasing the order of basis functions. Simultaneously, a global heat conduction matrix assembly formula is used for solving the problem to ensure the accuracy of numerical computation for the multi-faceted system.

[0137] Another embodiment of the present invention provides a transient thermal characteristic analysis and modeling method for a high-speed electric spindle system based on geometric analysis such as multi-faceted plates. S30 includes: defining convection boundary conditions and heat sources at relevant locations of the electric spindle, including forced convection between the shaft and air, natural air convection coefficient (marked with dark green lines), convection region between the stator cooling jacket and coolant, convection region for cooling bearings, forced convection between the stator and rotor, heating regions of the stator and rotor, and heating regions of the front and rear bearings; estimating the natural convection heat transfer coefficient of ambient air to obtain heat dissipation caused by local convection under different thermal conditions.

[0138] In a specific application example, such as Figure 10 As shown, convection boundary conditions and heat sources are defined at relevant locations on the electric spindle. For ease of differentiation, different convection surfaces are represented by lines of different colors in the axial section view based on the electric spindle model. h1 represents forced convection between the shaft and air, h2 represents the natural air convection coefficient (marked with a dark green line), h3 represents the convection region between the stator cooling jacket and the coolant, h4 and h6 show the convection regions used for cooling the bearings, and h5 represents forced convection between the stator and rotor. Simultaneously, each heat source region is represented by a rectangular block of a different color. s and qr These represent the stator and rotor heating regions, respectively. The heating of the front and rear bearings is analyzed as a single unit. Figure 10 The middle is represented as q a and q b Calculate the thermal power and convective cooling power, and estimate the natural convective heat transfer coefficient of the ambient air in order to obtain the heat dissipation caused by local convection under different thermal conditions.

[0139] When the electric spindle is in an unloaded state and the speed is constant at 15000 rpm, the magnitude of the boundary conditions used for numerical analysis of the electric spindle system is shown in Table 2.

[0140] Table 2 Boundary Condition Setting Table

[0141]

[0142] To verify the accuracy of the transient thermal analysis and thermal error model based on isogeometric analysis, experiments were conducted on a three-axis precision machine tool. The experimental subject was a high-speed milling electric spindle, and the actual performance of temperature field and thermal deformation was tested. The entire testing process followed the ISO230-3 standard.

[0143] The experimental system consists of two main parts: temperature measurement and spindle deformation measurement. Temperature measurement utilizes a magnetically attached PT100 temperature sensor, with eight temperature measurement points arranged on the outer wall of the electric spindle for real-time monitoring of the spindle and its components. Two additional temperature measurement points are also arranged to monitor the ambient temperature. The actual arrangement is as follows: Figure 14 As shown in Table 3, the specific locations of the 10 temperature measurement points are displayed. In order to ensure the accuracy of the measurement data at different locations of the spindle, two temperature sensors are arranged in each temperature measurement area.

[0144] In addition to temperature measurement, three KEYENCE LK-H008W laser displacement sensors were used to measure the thermal expansion of the electric spindle in the X, Y, and Z directions. The laser displacement sensors accurately capture the deformation of the electric spindle caused by thermal expansion during high-speed operation.

[0145] Temperature and displacement sensors work together to monitor the thermal performance of the electric spindle in real time at high speeds. This experiment simulated the electric spindle operating at 150,000 rpm, while simultaneously activating the spindle's cooling system to regulate the temperature.

[0146] The experimental results will be compared with the designed isogeometric analysis model and the simulation results of the commercial finite element software ANSYS 2021R1 to verify the accuracy and effectiveness of the model in predicting the temperature field and thermal deformation of the electric spindle system.

[0147] Table 3. Precise Location of Temperature Measurement Points

[0148]

[0149] In addition to temperature measurement, three KEYENCE LK-H008W laser displacement sensors are used to measure the thermal expansion of the electric spindle in the X, Y, and Z directions. The laser displacement sensors accurately capture the deformation of the electric spindle caused by thermal expansion during high-speed operation.

[0150] Temperature and displacement sensors work together to monitor the thermal performance of the electric spindle in real time at high speeds. The experiment simulated the electric spindle operating at 150,000 rpm, while the spindle's cooling system was activated to regulate the temperature.

[0151] The experimental results will be compared with the designed isogeometric analysis model and the simulation results of the commercial finite element software ANSYS 2021R1 to verify the accuracy and effectiveness of the model in predicting the temperature field and thermal deformation of the electric spindle system.

[0152] The analysis and discussion of the spindle thermal performance simulation results continue. The isogeometric analysis method is applied to the specific thermal conditions described above, and the results are compared with those obtained using commercial FEA software. The dimensions, materials, and boundary conditions of the physical model used in the established FEA analysis model are consistent with the isogeometric analysis model presented in this paper, to verify the accuracy and effectiveness of the proposed analytical thermal performance characterization method. Assuming a constant ambient temperature of 21.5℃ and a coolant temperature controlled at 20℃, these are used as the simulated thermal conditions for the spindle under no-load operation. The spindle speed is constant at 15000 rpm, and the transient temperature field and transient thermal deformation field of the spindle from 0 to 4800 s are simulated with a simulation step size of 60 s.

[0153] Furthermore, the simulation results of the transient temperature field were analyzed. The simulation results of the transient temperature field using the EFA and IGA methods are as follows: Figure 11 As shown in the figure. The figure displays the simulated temperature field results at 900s and 4800s, where... Figure 11 Image (a) shows the spindle temperature field simulated using the commercial FEA software ANSYS 2021R1, with 48,839 nodes. Figure 11 (b) shows the spindle temperature field calculated using the proposed IGA method, implemented independently in the MATLAB R2021 b platform, with a total of 45275 control points. The results show that the temperature of the central shaft structure, especially the part near the spindle motor rotor, is extremely high, and the temperature rise in the front bearing region is higher than that in the rear bearing region. Further quantitative analysis shows that the direct mapping between the NURBS parameter space and physical space in the boundary condition handling, and the global node numbering strategy of the multi-faceted model, ensure the C-value between adjacent NURBS faces. 0Continuity makes the temperature gradient transition at the geometric connection interface of the IGA model smoother, thus avoiding the problem of mesh mismatch between the physical model and the analysis model in FEA, which would otherwise lead to inconsistent simulated temperature field results.

[0154] To quantify the robustness of the FEA model and the multi-faceted IGA model, Figure 12 Temperature convergence characteristics under different grid densities were compared. Taking the highest temperature at a simulation time of 4800s as the comparison object, for FEA, the temperature change was less than 0.2℃ when the number of nodes in the grid was greater than 37452, while for IGA, the temperature change was less than 0.2℃ when the number of control points in the grid was greater than 35026. Both are less than the measurement accuracy.

[0155] Further analysis of the transient thermal deformation field simulation results was conducted. The thermal deformation of the principal shaft is caused by changes in its temperature field. Therefore, after obtaining the principal shaft temperature field, the IGA method can be used to calculate the thermal deformation. The strain ε caused by linear thermal expansion is... th Calculated using the following formula:

[0156] ε th =α(T-T0)[1 1 1 0 0 0] T (twenty four)

[0157] This formula indicates that the strain caused by thermal expansion is uniform and varies only in the translational direction (i.e., the effect of shear strain is not considered). Regarding the specific analysis process, an isogeometric analysis framework for thermo-mechanical coupling problems was independently developed in the MATLAB R2021 b platform, consisting of the following steps:

[0158] (1) The elastic matrix D is calculated based on the elastic modulus (E) and Poisson's ratio (ν) of the material, and then the strain-displacement matrix B is constructed using the derivative of the NURBS shape function.

[0159] (2) Assemble the stiffness matrix K and thermal load F by Gaussian integration. th .

[0160] (3) Since the front of the spindle is fixed by the flange, fixed boundary conditions are applied.

[0161] (4) Finally, solve the linear equation system Ku = F to obtain the displacement field.

[0162] To compare with the IGA analysis results, the transient thermal deformation of the spindle was calculated using the commercial FEM software ANSYS 2021R1. Figure 13 Numerical simulation results show that the axial thermal deformation trend of the spindle calculated by the IGA modeling method is consistent with the results of simulations by commercial EFA software. Figure 13In (a), corresponding to the FEA modeling method, the maximum deviation of thermal deformation at the Z-axis tool tip is 12.7 μm. Figure 13 In Figure (b), the IGA modeling method is represented. Its maximum deviation of thermal deformation at the Z-axis tool tip is 16.2 μm, which means that the FEA simulation result is smaller than the IGA simulation result. This also corresponds to the fact that in temperature field analysis, the maximum temperature of the FEA simulation result is smaller than the maximum temperature of the IGA simulation result.

[0163] It should be understood that the exemplary embodiments described herein are illustrative and not restrictive. Although one or more embodiments of the invention have been described in conjunction with the accompanying drawings, those skilled in the art will understand that various changes in form and detail may be made without departing from the spirit and scope of the invention as defined by the appended claims.

Claims

1. A method for analyzing and modeling the thermal properties of an electric spindle based on geometric analysis such as multi-faceted patches, characterized in that, Includes the following steps: S10, Perform geometric modeling of the high-speed electric spindle: Use 3D NURBS patch splicing technology to construct a multi-patch solid model of the high-speed electric spindle, including the spindle stator, rotor, bearings, and cooling pipes, to ensure geometric continuity and parametric consistency; S20, heat conduction analysis: Based on the Fourier heat conduction differential equation and combined with the high-order smoothness characteristics of IGA, a transient heat conduction control equation is constructed, and the temperature evolution process in the time domain is solved by the backward difference method. S30, Perform thermal error modeling: Introduce the constitutive relation of elasticity and the fixed support boundary to establish a thermo-mechanical coupling solution framework.

2. The method for transient thermal characteristic analysis and modeling of a high-speed electric spindle system based on multi-faceted geometric analysis as described in claim 1, characterized in that, S10 includes: The entire electric spindle system is transformed into a structure composed of multiple 3D non-uniform rational B-spline NURBS facets using multi-facet modeling technology. Each facet has a different thermal effect. The following formula applies when stitching together individual 3D NURBS patches: Where P i a and Let A be the control points, belonging to facet a and facet b respectively, and N be the set of contact surfaces of adjacent facets. Then (a, b) represents that facet a and facet b are adjacent. During the establishment of the integration domain, the parameter space of the NURBS patch is defined, and the existing overlapping or gap regions are corrected. The control points of each patch are globally numbered, and a mapping rule from local numbering to global numbering is adopted.

3. The method for transient thermal characteristic analysis and modeling of a high-speed electric spindle system based on multi-faceted geometric analysis as described in claim 2, characterized in that, S10 further includes: isogeometric spatial mapping relations under the three-dimensional NURBS entity, constructing a three-node index space from the given node vectors, and constructing the corresponding parameter space based on the units formed by the non-zero node intervals in the index space, thereby calculating the NURBS basis functions; the units are defined in the parameter space, and the Gaussian integration domain is defined in the parent space, then the entity contains two direct mappings in the spatial discretization process of IGA. In the first mapping, from the parent space supporting Gaussian integration... Transform to parameter space unit Through the Jacobian matrix Completed; in the second mapping, from the parameter space unit To physical space unit Ω e The transformation is achieved through the Jacobian matrix. Completed, among which Composed of the gradients of the basis functions, the Jacobian matrix J corresponding to the indirect mapping from the parent space to the physical space in 3D NURBS entities is obtained from these two direct mapping relationships. v In this way, PDEs in the physical model are transformed from strong form to weak form.

4. The transient thermal characteristic analysis and modeling method for high-speed electric spindle systems based on multi-faceted geometric analysis as described in claim 3, characterized in that, S20 includes: Fourier's law is the foundation of heat conduction. Based on the three-dimensional heat transfer domain Ω formed by the boundary Γ, its transient heat conduction governing differential equation is as follows: In the formula, T is the unknown temperature field function, Q is the intensity of the internal heat source, ρ is the material density, c is the specific heat capacity of the material, and λ is the specific heat capacity of the material. x , λ y , λ z These are the thermal conductivity of the material in the x, y, and z directions, respectively; The initial temperature conditions, isothermal boundary conditions, isothermal boundary conditions, and convective heat transfer boundary conditions of the transient heat transfer system are as follows: In the formula, T0 is the initial temperature of the given region, and T λ Here, q is the constant temperature of a given region, and h is the heat flux density of that region. c Γ is the convective heat transfer coefficient for a given region, and n is the unit outward normal vector on the boundary Γ. The temperature field function T at any point within a NURBS solid element is represented by the temperature value T at the control point. i,j,k and three-variable NURBS basis functions The interpolation form is as follows: In the formula, n1 = p b +1,n2=q b +1,n3=r b +1, n n =n1·n2·n3 represents the total number of control points or basis functions associated with this NURBS entity cell; Transforming equation (12) into the standard Galerkin weak form and combining it with equation (14) yields the discrete system equations based on isogeometric analysis, in the following form: In the formula, N T K T , {P} and {P} are respectively composed of NURBS solid units Ω e It is assembled from correlation matrices or vectors, and its unit correlation matrix or vector is described as follows: For the element heat capacity matrix, For the unit heat conduction matrix, P is the element heat conduction matrix related to the convective heat transfer boundary conditions. e Let be the element temperature load vectors. They can be written in Gaussian integral form over the standard integration domain in the parent space, as follows: In the formula, To be compatible with unit Ω e The relevant basis function matrix, NURBS unit Ω e The temperature gradient interpolation matrix within the element Ω, and related to the element Ω. e The temperature gradient interpolation matrix B corresponding to the k-th control point Tk (k = 1, 2, ..., n) n The format is as follows: Then, based on the chain rule and the Jacobian matrix... Further results were obtained: In addition to employing isogeometric discretization in the spatial domain, the transient heat conduction system based on isogeometric analysis also uses the backward difference method in the time domain. The system discretization equation (15) for transient heat conduction is further written as: In the formula, Δt is the time step.

5. The method for transient thermal characteristic analysis and modeling of a high-speed electric spindle system based on multi-faceted geometric analysis as described in claim 1, characterized in that, S20 includes: in numerical calculation, h-refinement and p-refinement strategies are adopted. h-refinement improves spatial resolution by increasing the density of control points, and p-refinement reduces numerical error by increasing the order of basis functions. The global heat conduction matrix assembly formula is used for solving.

6. The method for transient thermal characteristic analysis and modeling of a high-speed electric spindle system based on geometric analysis such as multi-faceted patches as described in claim 1, characterized in that, S30 includes: defining convection boundary conditions and heat sources at relevant locations of the electric spindle, including forced convection between the shaft and air, the natural convection coefficient of air (marked with dark green lines), the convection region between the stator cooling jacket and the coolant, the convection region used to cool the bearings, forced convection between the stator and rotor, the heating regions of the stator and rotor, and the heating regions of the front and rear bearings; estimating the natural convection heat transfer coefficient of ambient air in order to obtain the heat dissipation caused by local convection under different thermal conditions.

Images