A turntable bearing-torque motor thermal vibration coupling modeling method

Abstract

The application provides a turntable bearing-torque motor thermal vibration coupling modeling method and relates to the technical field of machine tool reliability analysis. The method comprises the following steps: calculating bearing restoring force, torque and elastic stiffness; giving a main shaft system dynamics balance equation; calculating unbalanced magnetic pull and unbalanced mass excitation; solving a displacement field, transferring load changes caused by the displacement field to a thermal model; updating contact force borne by the bearing by considering displacement-related load effects, calculating bearing friction heat generation and motor part heat generation; calculating bearing thermal deformation and thermal displacement; calculating thermal resistance between nodes, giving a temperature field balance equation through a thermal resistance network to realize rapid thermal-mechanical coupling analysis; solving the thermal balance equation to calculate temperature response; introducing the latest temperature into a dynamics model; updating grease viscosity and bearing thermal deformation; repeating the above steps to finally obtain main shaft displacement field and temperature field response. The application has high accuracy, strong universality, low test cost and high test efficiency.

Description

Technical Field

[0001] This invention relates to the field of machine tool reliability analysis technology, and in particular to a method for thermal vibration coupling modeling of rotary table bearing-torque motor. Background Technology

[0002] The motor spindle of a precision CNC machine tool is the core power unit, and its performance directly affects machining accuracy. The spindle system must withstand the combined effects of static loads, cutting loads, and thermal loads, among which the interaction between mechanical vibration and temperature rise is particularly critical: vibration alters the heat generation mechanism, while thermal expansion affects the structural contact state, forming a thermo-mechanical coupling effect. This interaction leads to asymmetric thermal deformation, which, when superimposed with dynamic loads, may induce self-excited oscillations, ultimately reducing machining accuracy.

[0003] Torque motors (a type of permanent magnet synchronous motor), as energy conversion devices, generate heat during operation that can induce mechanical vibrations and affect system stability. Therefore, energy consumption analysis requires a multidisciplinary approach combining electromagnetics, thermodynamics, and other methods. Turntable bearings, as high-precision three-row roller bearings, directly affect the accuracy of the rotating shaft due to frictional heat, and are commonly found in precision equipment such as military radar. When components within the system expand due to heat, the contact forces dynamically adjust to counteract the deformation, creating a continuous thermo-mechanical interaction. By establishing a thermo-vibration coupling model, the complex correlation between temperature changes and mechanical motion in the bearing-spindle system can be revealed.

[0004] Existing methods use separate coupled models to analyze the displacement and temperature fields. However, because they neglect the complex feedback mechanism between heat generation and vibration, they cannot accurately predict the nonlinear response of the spindle system, resulting in inaccurate calculation results. Furthermore, calculating the heat generation of the torque motor using finite element simulation requires significant computational resources and time, and the simplified complex geometry ignores the load distribution and contact state in the bearings, leading to large discrepancies between simulation and actual results. Summary of the Invention

[0005] The technical problem to be solved by the present invention is to provide a thermal vibration coupling modeling method for turntable bearing-torque motor in response to the shortcomings of the prior art. The model calculation results of this method are highly accurate and fast, providing a new theoretical tool for thermal error compensation of precision spindle systems.

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

[0007] A method for thermal vibration coupling modeling of a rotary table bearing and torque motor is proposed. This method studies the thermal vibration characteristics of a machine tool during the cutting process, deriving the displacement response of each degree of freedom of the spindle and the temperature field of the bearing and stator / rotor structure. For the thermal vibration coupling analysis, a partitioned iterative algorithm is used to achieve bidirectional coupling modeling. A closed-loop system of thermal vibration coupling is constructed through the coordinated solution of the dynamic equilibrium equation and the thermal equilibrium equation. The method includes the following steps:

[0008] Step 1: Considering the initial conditions and the influence of bearing thermal deformation, derive the bearing restoring force, torque, and elastic stiffness;

[0009] Step 2: Give the dynamic equilibrium equations of the main shaft system;

[0010] Step 3: Calculate the unbalanced magnetic pull and the unbalanced mass excitation;

[0011] Step 4: Solve the dynamic equations using the Newmark-β method to calculate the dynamic response during the time interval Δt, and solve for the displacement field; then transfer the load changes caused by the displacement field to the thermal model;

[0012] Step 5: Consider the displacement-related load effects, update the contact force on the bearing, calculate the frictional heat generated by the bearing, and the heat generated by various parts of the motor;

[0013] Step 6: Calculate the bearing's thermal deformation and thermal displacement;

[0014] Step 7: Calculate the thermal resistance R between each node. ij Furthermore, rapid thermo-mechanical coupling analysis is achieved by deriving the temperature field equilibrium equation through a thermal resistance network.

[0015] Step 8: Solve the heat balance equation using the Runge-Kutta method to calculate the temperature response over time interval Δt; import the latest temperature into the kinetic model;

[0016] Step 9: Consider the influence of temperature-related parameters, and update the grease viscosity and bearing thermal deformation;

[0017] Step 10: Determine if the temperature change ΔT at the hot node is greater than 10. -3 If the value is greater than 10, repeat steps 1 to 9 until the node temperature change ΔT is less than 10. -3 The loop ends, and the final response of the principal shaft displacement field and temperature field is obtained.

[0018] The beneficial effects of adopting the above technical solution are as follows: The rotary table bearing-torque motor thermal vibration coupling modeling method provided by this invention innovatively constructs an energy transfer closed loop through a two-way coupling mechanism: Step 3 reflects the force-thermal coupling of vibration to the temperature field (involving nonlinear effects such as frictional heat generation), and Step 5 realizes the thermoelastic coupling of temperature to vibration (including thermal expansion effects). Compared with traditional co-simulation, the breakthrough of this invention lies in: establishing an improved thermal balance equation that includes nonlinear vibration characteristics, developing a physical bridge for the transfer of displacement-temperature dual-field parameters, and proposing a convergence criterion based on energy conservation. This method inherits and develops the multiphysics coupling modeling theory. Its core lies not in a general iterative framework, but in a customized equation system constructed for the special working conditions of the torque motor spindle system, providing a new theoretical tool for thermal error compensation of precision spindle systems. Attached Figure Description

[0019] Figure 1 A flowchart of the turntable bearing-torque motor thermal vibration coupling modeling method provided in an embodiment of the present invention;

[0020] Figure 2 This is a schematic diagram showing the position angle of the rollers of the turntable bearing in the turntable bearing-torque motor provided in an embodiment of the present invention.

[0021] Figure 3 This is a schematic diagram of the position angle of the angular contact bearing balls of the turntable bearing-torque motor provided in an embodiment of the present invention;

[0022] Figure 4 A schematic diagram of the inner and outer raceway curvature centers of the first angular contact bearing on the left side of the turntable bearing-torque motor provided in an embodiment of the present invention;

[0023] Figure 5 A schematic diagram of the inner and outer raceway curvature centers of the second left angular contact bearing of the turntable bearing-torque motor provided in an embodiment of the present invention;

[0024] Figure 6 This is a schematic diagram of the overall structure of the turntable bearing-torque motor spindle provided in an embodiment of the present invention;

[0025] Figure 7 This is a schematic diagram of stator and rotor eccentricity of a turntable bearing-torque motor provided in an embodiment of the present invention;

[0026] Figure 8 A schematic diagram of the electromagnetic phasors of a turntable bearing-torque motor provided in an embodiment of the present invention;

[0027] Figure 9 This is a schematic diagram of the nodal thermal grid of a turntable bearing-torque motor provided in an embodiment of the present invention; wherein the magnified view is the magnified thermal resistance perspective of the angular contact bearing and the turntable bearing;

[0028] Figure 10 This is a schematic diagram of the thermal resistance grid of the turntable bearing-torque motor provided in an embodiment of the present invention;

[0029] Figure 11 This is a schematic diagram of the thermal vibration coupling iteration of the turntable bearing-torque motor provided in an embodiment of the present invention. Detailed Implementation

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

[0031] A rotary table bearing-torque motor thermal vibration coupling modeling method is proposed to study the thermal vibration characteristics of machine tools during the cutting process, obtain the displacement response of each degree of freedom of the spindle, and the temperature field of the bearing and stator / rotor structure. For thermal vibration coupling analysis, a partitioned iterative algorithm is used to achieve bidirectional coupling modeling. A closed-loop system of thermal vibration coupling is constructed by synergistically solving the dynamic equilibrium equation and the thermal equilibrium equation. Figure 1 As shown, the method of this embodiment is described below.

[0032] Step 1: Considering the initial conditions and the influence of bearing thermal deformation, derive the bearing restoring force, torque, and elastic stiffness. The specific method is as follows.

[0033] Step 1.1: Calculate the position angle.

[0034] Figure 2 This is a schematic diagram of the position angle of the turntable bearing. Calculate using the following formula:

[0035] (1);

[0036] in, These are the bearing position angles of each row of turntables; , ω is the angular velocity of each row of cages in the rotary table bearing, ω is the angular velocity of the spindle rotation, and t is the time variable. It is the pitch circle diameter of the radial rollers. It is the radial roller diameter. This refers to the number of rollers in each row. The superscripts u / d / m represent the parameters corresponding to the upper row, lower row, or radial row of the turntable bearing.

[0037] Figure 3 This is a schematic diagram of the position angle of an angular contact bearing. Calculate using the following formula:

[0038] (2);

[0039] Where, ωc ω is the angular velocity of the angular contact bearing cage, ω is the angular velocity of the main spindle, and Z is the angular velocity of the main spindle. b D represents the number of bearing balls. i D o These refer to the inner and outer diameters of the bearing, respectively.

[0040] Step 1.2: Calculate the bearing deformation.

[0041] The axial variation δ of the center of curvature of the inner raceway groove of the i-th roller in each row of the rotary table bearing a Ru / Rd The radial variation value δ r m and rotational displacement θ u / d As shown in the following formula:

[0042] (3);

[0043] Among them, e R This is the axial distance between the center of the roller and the geometric center of the spindle. It is the pitch circle diameter of the upper and lower rows of rollers; where [x,y,z,θ] x ,θ y These represent the dynamic displacements in the x, y, and z directions of the principal axis and the angular displacements in the x and y directions, respectively.

[0044] For angular contact bearings, such as Figure 4 and Figure 5 As shown, the outer ring curvature center is fixed, and the radial displacement of the inner ring curvature center under dynamic load is... and axial displacement As shown in the following formula:

[0045] (4);

[0046] Among them, e L It is the axial distance between the center of the ball bearings and the geometric center of the spindle. It is the orbital radius of the center of curvature of the inner raceway groove; D m and D b These are the pitch circle diameter and ball diameter of the angular contact bearing, α0 is the initial contact angle, and f is the initial contact angle. i It is the inner raceway curvature radius coefficient.

[0047] Step 1.3: Calculate the normal approach and contact deformation between the two raceways.

[0048] The rollers are discretized using a slicing method, and the effective contact length l of each row of rollers is determined. u / d / m Divide into n equal parts k Each slice unit has an axial thickness expressed as: .

[0049] Normal approaching amount between the two raceway surfaces at the center of the k-th piece of the i-th roller in the upper row for:

[0050] (5);

[0051] Wherein, the superscript u represents the upper row of rollers, u a0 For the initial axial clearance, u a This is due to axial thermal expansion.

[0052] Normal approaching amount between the two raceway surfaces at the center of the k-th piece of the i-th roller in the lower row for:

[0053] (6);

[0054] The superscript d represents the lower row of rollers.

[0055] The normal approach of the two raceway surfaces at the center of the k-th piece of the i-th roller in the radial arrangement for:

[0056] (7);

[0057] Where the superscript m represents radial rollers, u r0 For the initial radial clearance, u r This is due to radial thermal expansion.

[0058] For angular contact bearings, such as Figure 4 and Figure 5 As shown, considering only the dynamic displacement of the inner ring, the axial distance between the centers of curvature of the inner and outer rings is... and radial distance for:

[0059] (8);

[0060] Among them, O i O o It is the distance between the centers of curvature of the inner and outer raceways, α F It is the bearing contact angle after preload.

[0061] Considering the dynamic displacement and centrifugal effect of the inner ring, as well as the thermal expansion of the bearing, the axial distance between the centers of curvature of the inner and outer rings. and radial distance Determined by the following formula:

[0062] (9);

[0063] Among them, the '+' and '-' of '±' correspond to the first and second angular contact bearings on the left end, respectively; δ jzTDand δ jrTD These are the axial and radial displacements of the inner ring curvature center, taking thermal expansion into account.

[0064] Based on the geometric compatibility relationship and considering the influence of external dynamic loads, the distance between the centers of curvature of the inner and outer rings on the j-th rolling element is obtained. and dynamic normal contact deformation for:

[0065] (10);

[0066] Considering the bearing expansion caused by temperature rise and centrifugal force, the actual distance A between the centers of curvature is... jT and the actual contact deformation δ at the j-th bearing ball jT With contact angle α jT It is given by the following formula:

[0067] (11);

[0068] Where, δ F It is pre-tightening deformation.

[0069] Step 1.4: Calculate the bearing's elastic restoring force, torque, and elastic stiffness.

[0070] The normal load on each roller of the turntable bearing is:

[0071] (12);

[0072] Considering dynamic loads and thermal expansion effects, the nonlinear stiffness characteristics of a bearing are characterized as a multivariable function of dynamic displacement, rotational speed, and temperature field. The nonlinear dynamic restoring force F of the bearing is obtained through the force balance equations of the inner ring. x R F y R F z R and torque M x R M y R for:

[0073] (13);

[0074] According to Hertzian contact theory, the dynamic normal contact force at the j-th ball of an angular contact bearing... for:

[0075] (14);

[0076] in, This indicates the Hertzian contact stiffness.

[0077] Actual normal contact force considering thermal effects for:

[0078] (15);

[0079] Finally, the nonlinear dynamic restoring force F of the angular contact bearing is obtained through the force balance equation of the inner ring. x L F y L F z L and torque M x L M y L :

[0080] (16);

[0081] in, It is the lever arm applied to the inner ring of the bearing.

[0082] Based on the load-deformation relationship, the nonlinear dynamic stiffness vector of the bearing is derived from the following formula:

[0083] (17);

[0084] Among them, only rigid body displacement is considered, and the elastic body displacement along the principal axis is ignored, [δ x , δ y , δ z , φ x , φ y ] T These are the dynamic displacements of the bearing in the x, y, and z directions and the angular displacements in the x and y rotational directions; then, the bearing stiffness matrix is ​​obtained from equation (17). for:

[0085] (18);

[0086] Among them, [k xL1 , k yL1 , k zL1 , k φxL1 , k φyL1 ]、[k xL2 , k yL2 , k zL2 , k φxL2 , k φyL2 ] and [k xR ,k yR , k zR , k φxR , k φyRThese are the displacement stiffness of the two angular contact bearings on the left and the rotary table bearing on the right in the x, y, and z directions, and the angular stiffness in the x and y rotation directions, respectively.

[0087] Step 2: Give the dynamic equilibrium equations of the main shaft system.

[0088] Figure 6 This is a model diagram of the motor. Considering the dynamic characteristics of the spindle system, a rigid body displacement assumption model is adopted, neglecting its elastic deformation components. The study focuses particularly on the thermal-vibration coupling effect of the bearing pair, positioning it as the core link in the system's energy transfer. To simplify engineering calculations, the displacement difference caused by the small distance between adjacent left bearings is ignored, and their displacement responses are assumed to be synchronous. Under external time-varying excitation loads, the system exhibits significant nonlinear vibration response characteristics. In summary, its dynamic equation can be expressed as:

[0089] (19);

[0090] Where m is the mass of the spindle system. Let the rotational inertia be the diameter of the system. R is the polar moment of inertia of the system. d The extreme diameter; F ij and M ij These represent the restoring force and restoring torque of the left and right bearings, respectively. i = x, y or z, j = 1, 2 or R, where 1, 2, and R represent the two angular contact bearings on the left and the turntable bearing, respectively. , , , , This indicates an unbalanced quality incentive; , , , , This indicates an unbalanced magnetic pull excitation; L1 represents the axial distance between the two bearings, and L2 represents the axial distance from the angular contact bearing to the center.

[0091] The above formula can be written in matrix form as follows:

[0092] (20);

[0093] Where q=[x,y,z,θ] x ,θ y ] T This represents the dynamic displacement of the principal axis, M is the mass matrix, and C = [c1 c2 c3 c4 c5]. T G is the damping matrix, F is the gyroscope matrix, and F is the damping matrix. b Q is the bearing restoring force matrix, and Q is the external force matrix.

[0094] Step 3: Calculate the unbalanced magnetic pull and the unbalanced mass excitation. The specific method is as follows.

[0095] Step 3.1: As Figure 7 and Figure 8 The diagram shows the unbalanced magnetic pull caused by stator-rotor eccentricity, and the phasor diagram of the motor's magnetomotive force. Rotor alignment deviation is mainly caused by spatial offset and can be divided into three types: static, dynamic, and mixed eccentricity. Static eccentricity is caused by systematic factors such as assembly errors; it is persistent and difficult to eliminate, inducing unbalanced magnetic pull (UMP) and weakening system stability. Dynamic eccentricity is coupled with the spindle vibration characteristics, and the resulting harmonics nonlinearly amplify the vibration response, forming a vicious cycle. Mixed eccentricity is a superposition of the first two, and its electromechanical coupling effect is the most complex.

[0096] The magnetomotive force (MOF) analysis of the stator windings of a permanent magnet synchronous motor (PMSM) is similar to that of an induction motor. PMSMs are typically powered by pulse-width modulation (PWM) voltage sources. However, the armature current is an approximately three-phase symmetrical sinusoidal waveform. According to motor winding theory, the fundamental magnetomotive force of the three-phase stator windings... Represented as:

[0097] (twenty one);

[0098] Among them, F sm The fundamental magnetomotive force amplitude of the excitation current for the stator armature reaction current is given by k, where N is the number of turns in series per pole per phase of the stator, and k is the number of turns in series per phase per pole of the stator. w It is the stator winding distribution factor, P is the number of pole pairs of the motor, and I is the stator winding distribution factor. smax It is the maximum value of the stator excitation current, ω cur α is the angular velocity of the rotating magnetic field, and α is the azimuth angle of the air gap space.

[0099] According to the magnetic circuit principle of permanent magnet motors, the permanent magnet is considered a constant source of magnetomotive force; the fundamental magnetomotive force of a surface permanent magnet motor... Represented as:

[0100] (twenty two);

[0101] Among them, F rm B represents the magnetomotive force amplitude of the permanent magnet. r For remanence of permanent materials, h m denoted as the magnetization thickness of the permanent magnet, and μ0 as the vacuum permeability.

[0102] Electromagnetic excitation is generated by the air gap magnetic field of the motor, which is determined by the magnetic force between the stator and rotor and the air gap permeability; the fundamental composite magnetomotive force is:

[0103] (twenty three);

[0104] When stationary, the rotor's axial center is the origin of the coordinate system; elastic deformation of the rotor is ignored and it is considered rigid. The initial static eccentricity of the rotor is [x0, y0, z0, θ]. x0 , θ y0 ] T The dynamic displacement of the rotor's axial center relative to the origin is [x, y, z, θ]. x ,θ y ] T Based on geometric relationships and considering the influence of angular displacement eccentricity, the dynamic displacement q of any point on the rotor is obtained. r for:

[0105] (twenty four);

[0106] Among them, L r It is the axial length of the rotor.

[0107] Based on the geometric compatibility relationship, the actual air gap length δ(α,t,Z) related to the rotor's dynamic displacement is derived from the following formula:

[0108] (25);

[0109] in, The air gap length when the stator and rotor coincide; Z is the rotor axial eccentricity, r0 is the rotor static eccentricity, β is the static eccentricity direction angle, and R... ro It is the outer diameter of the rotor.

[0110] The tilt angle of the rotor surface normal magnetic tensile stress relative to the stator radial direction is expressed by the following formula:

[0111] (26);

[0112] According to Ampere's circuital law, the permeability of the air gap per unit area is obtained by the following formula:

[0113] (27);

[0114] Where, k c It is the groove coefficient; for a closed groove, k c =1; Saturation coefficient k of the magnetic circuit for magnetic saturation u Consider the magnetic circuit length perpendicular to the rotor surface. for:

[0115] (28);

[0116] The air gap magnetic flux density is:

[0117] (29);

[0118] According to the magnetic stress formula, the normal magnetic tensile stress on the rotor surface is obtained as follows:

[0119] (30);

[0120] Projecting the magnetic tensile stress along the stator radial and axial directions and integrating it on the rotor surface, the unbalanced electromagnetic force and torque acting on the rotor are obtained as follows:

[0121] (31).

[0122] Step 3.2: Calculate the unbalanced mass excitation.

[0123] The centrifugal force caused by the eccentric mass projected along the radial and axial directions of the stator yields the unbalanced mass excitation as follows:

[0124] (32);

[0125] Where m is the mass of the spindle system, d e It is the mass eccentricity, and v is the spatial inclination of the centrifugal force relative to the stator radial direction.

[0126] Step 4: Solve the dynamic equations based on the Newmark-β method to calculate the dynamic response during the time interval Δt, and solve for the displacement field; then transfer the load changes caused by the displacement field to the thermal model.

[0127] The analytical functions of the time-varying unbalanced mass excitation, time-varying unbalanced magnetic pull, nonlinear dynamic bearing stiffness, and bearing restoring force in step 3 are related to the dynamic displacement q of the spindle. The dynamic equations are solved using the Newmark-β method within each time step Δt to obtain the spindle displacement field q. t+Δt The obtained displacement field is then used to update the heating of each component of the spindle system in subsequent steps.

[0128] Step 5: Consider the displacement-related load effect to update the contact force on the bearing, calculate the frictional heat generated by the bearing, and the heat generated by various parts of the motor, as shown below.

[0129] Step 5.1: Calculate motor heat generation; the internal losses of the motor are divided into stator losses (copper losses and iron losses, of which iron losses include hysteresis losses, eddy current losses and additional losses) and rotor losses (mainly permanent magnet eddy current losses and mechanical losses, iron losses can be ignored).

[0130] According to Joule's law, the copper loss P caused by the heating of the conductors in the stator windings... cus Represented as:

[0131] (33);

[0132] Where, mcus =3 indicates a three-phase motor, I represents the effective value of the current through each phase winding, and R represents the effective value of the resistance in each phase winding. The effect of temperature on the winding resistance is ignored.

[0133] In actual operating conditions, the non-ideal characteristics of three-phase power supply can induce a rotating magnetic field, resulting in iron losses having a dual source: the direct effect of the alternating magnetic field and the superimposed dynamic influence of the rotating magnetic field. Assuming the magnetic density of the motor core is uniformly distributed during operation, based on Bertotti's iron loss separation theory, the stator iron loss P... Fes The losses are categorized into three types: hysteresis loss, eddy current loss, and additional loss. The calculation model expression is as follows:

[0134] (34);

[0135] Among them, P h It is the unit mass hysteresis loss, P e It is the eddy current loss per unit mass, P exc It is the additional loss per unit mass; k h The value represents the hysteresis loss coefficient, f is the frequency, and B is the value of the frequency. m It is the magnetic flux density amplitude; the eddy current loss coefficient is The additional loss coefficient is σ is the electrical conductivity of the core material, d is the thickness of a single layer of silicon steel sheet, ρ is the core density, and G and V0 are performance parameters of the silicon steel sheet, where G is a dimensionless coefficient and V0 is a coefficient related to B. m The relevant statistical characteristic parameters are as follows: S is the cross-sectional area of ​​the silicon steel sheet.

[0136] Hysteresis loss refers to the energy consumed by the core material during repeated magnetization due to hysteresis, and is positively correlated with the area of ​​the hysteresis loop. Eddy current loss is caused by the induced circulation of alternating magnetic flux within the silicon steel sheet; the skin effect is neglected in the calculation. Additional loss refers to the remaining portion of the core material in an alternating magnetic field after removing hysteresis and eddy current losses. This portion of loss encompasses complex mechanisms such as domain wall motion and microscopic eddy currents, and is closely related to the material's crystal structure and dynamic magnetization process, making it difficult to accurately describe with simple physical models.

[0137] The rotor of the permanent magnet synchronous motor adopts a composite structure of permanent magnets and insulated silicon steel sheets. An insulation lamination process is used to block the transverse eddy current path, reducing iron losses (ignoring rotor core eddy current losses). Stator slotting causes periodic fluctuations in air gap permeability. According to the principle of electromagnetic coupling, the interaction between the stator harmonic magnetomotive force and the periodic permeability generates higher-order harmonic magnetic fields. When these harmonics penetrate the isotropic permanent magnets, they form closed eddy current loops, thus generating additional Joule losses. Based on the above electromagnetic mechanism, the permanent magnet eddy current loss P... wr The parsing expression is:

[0138] (35);

[0139] Where, k r ρ is the proportionality constant of the electromotive force. w It is the resistivity of permanent magnets. It is the volume of the permanent magnet, where h m L is the thickness of the permanent magnet. a L is the axial length. m This refers to the horizontal width.

[0140] The power loss P between the rotating air gap flow and the stator and rotor friction. w It is derived from the following formula:

[0141] (36);

[0142] in, It is the rotor rotation frequency, where n is the spindle speed, η air It is aerodynamic viscosity.

[0143] In summary, the stator heating P of the torque motor is calculated. s and rotor heating P r :

[0144] (37);

[0145] Where, m Fe It is the total mass of the stator silicon steel sheets, N w It refers to the number of permanent magnets in the rotor.

[0146] Step 5.2: Calculate bearing heat generation.

[0147] Under the influence of material elastic hysteresis, when the rollers of a turntable bearing roll along the inner and outer raceways, their contact area exhibits an asymmetric stress distribution: the rolling resistance torque generated by the normal contact stress in the leading edge region is significantly higher than the propulsion torque in the trailing edge region. This asymmetry in stress distribution directly leads to the generation of rolling friction resistance, and the energy loss between each row of rollers and the raceway is essentially caused by an energy dissipation mechanism dominated by the material elastic hysteresis effect. The following formula is used to derive:

[0148] (38);

[0149] in, It is the distance from the center of the k-th piece of the i-th roller to the center of the bearing. , ξ=1 is the material's elastic hysteresis coefficient, E is the combined elastic coefficient of the two contacting bodies, and ωI / E It refers to the angular velocity of the inner or outer ring.

[0150] During bearing operation, the difference in radii of gyration at the contact points between the rollers and the inner and outer races causes differential sliding at the contact interface. This relative slippage, induced by geometric kinematics, creates micro-friction on the contact surface. The corresponding energy dissipation mechanism is known in tribology as differential sliding friction. This frictional effect, together with rolling friction, constitutes an important component of the bearing's frictional torque. Its quantitative analysis is based on Hertz contact theory and an elastohydrodynamic lubrication model, and is derived from the following formula:

[0151] (39);

[0152] Where u=0.96 is the oil film drag coefficient. It is the relative sliding speed between the upper and lower rows of rollers and the inner raceway. It is the relative sliding speed between the upper and lower rows of rollers and the outer raceway. It is the rotational angular velocity of the upper and lower rows of rollers; for radial rollers, since the linear velocity along the contact line between the roller and the raceway remains constant, the relative sliding between the roller and the raceway is not considered.

[0153] The grease viscosity effect leads to frictional power loss in the contact area between the rolling elements and the raceway. for:

[0154] (40);

[0155] Where E0 is the equivalent elastic modulus and R0 is the equivalent radius of curvature; These are the speed parameters of the rolling elements and the inner and outer rings, where η0 is the dynamic viscosity of the grease. It is the average speed of the roller between the pitch circle and the inner raceway. It is the angular velocity of the radial rollers' rotation. It is the average speed between the roller and the outer raceway; These are the material parameters of the rolling elements and the inner and outer rings, among which... It is the viscosity-pressure coefficient of the lubricating grease; These are the load parameters for the rollers and the inner and outer rings.

[0156] In summary, the total heat generated by the inner rings of each row of rollers in the turntable bearing is... Total heat generation in the outer ring They are respectively:

[0157] (41);

[0158] To improve the heat generation of angular contact bearings, thermal expansion, nonlinear dynamic loads, and lubricant viscosity variations are considered. The total frictional torque of the bearing is calculated using the Palmgren formula. for:

[0159] (42);

[0160] Among them, M l and M v These are the frictional torques, M, caused by the applied load and the viscous friction of the lubricant, respectively. ω It is the rotational friction torque of the rolling element.

[0161] Frictional torque M caused by the applied load on the bearing l The analysis is as follows: the curvature compatibility of the rolling contact area (i.e., the geometric size effect), the elastohydrolubrication state of the Hertzian contact ellipse (i.e., the deformation energy storage), the Coulomb friction work in the micro-slip zone (i.e., the rolling to sliding friction ratio), and the energy dissipation field constituted by the load vector space distribution are calculated using the following formula:

[0162] (43);

[0163] Among them, f l F is a coefficient that takes into account the effects of load and bearing type. l This refers to the load on the bearing. For angular contact bearings, the calculation is obtained from the following formula:

[0164] (44);

[0165] Among them, F s It is the equivalent static load of the bearing. C s This is the basic static load rating; F r and F a These are radial and axial loads, calculated using the following formulas:

[0166] (45);

[0167] Frictional torque M caused by viscous friction of lubricant v Essentially, it is the energy loss generated by the shearing action of the oil film during bearing movement. Its magnitude is affected by various parameters such as lubricant dynamic viscosity, bearing speed, oil film thickness distribution, and effective friction area. For angular contact bearings, M v The calculation is as follows:

[0168] (46);

[0169] Among them, f v This is a coefficient that takes into account the bearing structure and lubrication method; for grease lubrication, it is taken as 2. , represents the relationship between the viscosity of the grease and temperature, where ν0 represents the kinematic viscosity of the grease at temperature T0, β is the temperature coefficient of the grease, and T is the temperature of the grease.

[0170] Rotational friction torque M of the rolling element ω This refers to the torque caused by friction between contact interfaces when a rolling element rolls along a raceway or other contact surface. This torque characteristic reflects the mechanical energy loss characteristics of the system due to surface friction during rotation, and the formula is as follows:

[0171] (47);

[0172] Where, μ ω It is the coefficient of rotational friction between the rolling elements and the inner ring raceway, which is taken as 0.0012~0.002 for ball bearings. a and b are the major and minor semi-axises of the rolling element deformation projected onto the inner ring surface, respectively. Then the total frictional heat generated by the bearing is P. b for:

[0173] (48);

[0174] To facilitate thermal network analysis, the heat distributed between the inner and outer raceways is converted into the heat distributed on the surfaces of the inner and outer raceways, expressed by the following formula:

[0175] (49);

[0176] To reduce friction, cages without guide structures are widely used. Ignoring slippage, the frictional loss P between the ball and the cage is considered. c We obtain it from the following formula:

[0177] (50);

[0178] Where, m c It is the cage mass, μ c It is the coefficient of sliding friction between the ball and the cage.

[0179] Step 6: Calculate the bearing's thermal deformation and thermal displacement, as shown below.

[0180] Step 6.1: Calculate the radial thermal expansion of the turntable bearing.

[0181] The inner ring of the rotary table bearing has an interference fit with the shaft, while the outer ring is fixed to a rigid bearing housing with no radial displacement. The rolling elements perform pure rolling between the inner and outer rings. The bearing is typically installed at ambient temperature, but its operating temperature may be higher than ambient temperature. Increased temperature will cause changes in the clearance. Assuming the material is isotropic and its thermal expansion is linearly distributed, the radial expansion u of the outer ring can be derived from the thermal expansion formula for a hollow cylinder. oand inner radial expansion u i and radial roller expansion for:

[0182] (51);

[0183] in, , and These are the coefficients of thermal expansion of the bearing's inner ring, rollers, and outer ring, respectively, ΔT. i and These are the temperature rises of the bearing inner ring and bearing rollers, respectively, ΔT o and ΔT h These are the temperature rises of the bearing outer ring and the bearing housing, respectively, d i and d o These are the inner and outer ring diameters of the turntable bearing, d h is the inner diameter of the bearing housing, and μ is Poisson's ratio.

[0184] According to equation (51), the radial thermal expansion u of the bearing is obtained. r for:

[0185] (52);

[0186] Step 6.2: Calculate the axial thermal expansion of the turntable bearing.

[0187] Axial expansion of bearing inner ring a i axial expansion of the outer ring a o and the axial expansion of the upper and lower rows of rollers a b u(d) Use the following formula:

[0188] (53);

[0189] Among them, l i and l o These represent the widths of the inner and outer rings of the bearing, respectively. It is the diameter of the upper and lower rows of rollers. It refers to the temperature rise of the upper and lower rollers.

[0190] According to equation (53), the axial thermal expansion u of the bearing is obtained. a for:

[0191] (54);

[0192] Step 6.3: The thermal deformation of the balls in an angular contact bearing is derived from the following formula:

[0193] (55);

[0194] Among them, D bIt is the diameter of the bearing balls, ΔT b It ensures uniform temperature rise of the ball bearings.

[0195] The difference in thermal expansion between the shaft and the inner ring leads to a redistribution of contact pressure. The assembly stress generated by the interference fit is superimposed on the thermal stress, jointly affecting the radial deformation of the inner ring raceway. Considering the thermal expansion effect between the shaft and the bearing inner ring, as well as the interference fit between the inner ring and the shaft journal, the radial thermal deformation U of the inner ring raceway of the bearing, as a hollow thin-walled circular ring structure, is... ir The calculation is performed using the following formula:

[0196] (56);

[0197] Among them, D i D is the inner diameter of the bearing. io U is the diameter of the bearing's inner raceway; io U i U ro The thermal deformations of the inner ring raceway, bearing inner diameter, and journal diameter are respectively calculated using the following formulas:

[0198] (57);

[0199] Among them, D ro The journal diameter is... μ is the coefficient of thermal expansion of the shaft. r For the axial Poisson's ratio, ΔT ro This ensures a uniform temperature rise for the journal.

[0200] Considering the thermal effects of the bearing housing and outer ring, as well as the clearance fit, the radial thermal deformation U of the outer ring raceway... oh Calculate using the following formula:

[0201] (58);

[0202] Among them, D o D is the outer diameter of the bearing. oi U is the outer raceway diameter, δ1 is the clearance fit between the outer ring and the bearing housing; oi U o U hi The thermal deformations of the outer raceway, bearing outer diameter, and bearing housing inner diameter are calculated using the following formulas:

[0203] (59);

[0204] Among them, D hi The inner diameter of the bearing housing. μ is the coefficient of thermal expansion of the bearing housing. h For the bearing housing Poisson's ratio, ΔT hi This ensures a uniform temperature rise in the bearing housing.

[0205] Step 6.4: The thermal displacement problem caused by the thermal expansion of the angular contact bearing is equivalent to a static load-displacement problem. Based on the geometric coordinate relationship, the radial displacement of the inner ring raceway curvature center caused by the thermal deformation of the bearing balls is obtained. and axial displacement They are respectively:

[0206] (60);

[0207] in, It refers to the bearing contact angle.

[0208] Radial displacement of the inner raceway curvature center caused by radial thermal deformation of the outer raceway and axial displacement They are obtained from the following formulas:

[0209] (61);

[0210] Considering the expansion of the inner ring due to centrifugal force, the radial displacement of the inner ring raceway curvature center caused by radial thermal deformation of the inner ring raceway is derived by the following formula:

[0211] (62);

[0212] Among them, U cent For the centrifugal expansion of the inner ring, These are the bearing pitch circle diameters, ρ, E, and μ. i These are the material density, elastic modulus, and Poisson's ratio of the inner ring, respectively.

[0213] like Figure 6 As shown, the angular contact bearings are arranged back-to-back. In the thermal displacement analysis, the center distance between the inner and outer ring curvatures of the left and right bearings shows opposite trends. Based on the assumption that the outer ring curvature center is fixed, the axial and radial displacement analysis of the inner ring curvature center is established by coupling the bearing thermal expansion effect with the centrifugal force of the inner ring, as shown in the following equation:

[0214] (63);

[0215] In this context, '-' corresponds to the first angular contact bearing on the left, and '+' corresponds to the second angular contact bearing on the left.

[0216] Step 7: Calculate the thermal resistance R between each node. ij Furthermore, a rapid thermo-mechanical coupling analysis is achieved by deriving the temperature field equilibrium equation through a thermal resistance network, as detailed below.

[0217] Step 7.1: Calculate thermal resistance.

[0218] Each node transfers heat through conduction, convection, and radiation. Neglecting the effect of radiation, calculate the axial thermal resistance R between each node. a Radial thermal resistance R r and convection thermal resistance R v As shown in the following formula:

[0219] (64);

[0220] Among them, R out and R in These are the outer radius and inner radius of the hollow cylinder, respectively, K. D A is the thermal conductivity of the solid, A is the convective heat transfer area, and h is the convective heat transfer coefficient, which is estimated by the Nusselt number N. u Estimate; L is the characteristic length of heat exchange, K l is the thermal conductivity of the liquid.

[0221] Step 7.2: Construct a thermal resistance network, give the temperature field equilibrium equation, and realize rapid thermal vibration coupling analysis.

[0222] like Figure 9 As shown, a 71-node thermal network model is constructed to analyze the temperature field of the torque motor spindle system, including 15 key heat-generating nodes and 59 nodes whose temperatures need to be calculated. Compared to the widely used finite element method, thermal network modeling significantly improves computational efficiency while maintaining temperature prediction accuracy: although the finite element method has the advantage of high accuracy, its complex modeling process and lengthy computation time limit its applicability to theoretical research. The model in this embodiment achieves rapid thermo-mechanical coupling analysis through a thermal resistance network, such as... Figure 10 As shown, the calculated thermal deformation of the components is imported into the dynamic model to study the thermally induced vibration response and its coupling effect with the system heating.

[0223] Thermal resistance networks consist of three core types of thermal resistance: thermal convection resistance, which characterizes the heat transfer efficiency between the solid surface and the cooling medium; thermal conduction resistance, which reflects the thermal conductivity characteristics inside the solid material; and thermal contact resistance, which describes the heat transfer mechanism of the contact interface and gap.

[0224] The model makes the following assumptions: the heat conduction unit is characterized by the central node to represent the average temperature field, the surface temperature of the heat convection unit is approximately equal to the central temperature, the contact surface temperature is dominated by the central node, the influence of thermal radiation and the temperature rise of the cooling medium is ignored, and the ambient temperature is kept constant as the boundary condition.

[0225] The thermal balance equations for each node in the angular contact bearing are as follows:

[0226] (64);

[0227] The heat balance equations for each node in the motor are as follows:

[0228] (65);

[0229] The heat balance equations for each node of the turntable bearing are as follows:

[0230] (66);

[0231] in, T i T is the temperature of the i-th node. j R is the temperature of the j-th node adjacent to the i-th node; i_j P is the total thermal resistance between the i-th node and the j-th node; i It is the heat generated per unit time at the i-th node, including P. i1 P c1 P o1 P i2 P c2 P o2 These are the heat generated by the inner ring, the frictional heat of the cage, and the heat generated by the outer ring of the two angular contact bearings, respectively. r P s P w These are the friction losses of the torque motor rotor, stator, and airflow, respectively, P I u P E u P I d P E d P I m P E m These refer to the heating of the inner and outer rings of the upper and lower rows and the radial row of the turntable bearings; c i It is the specific heat capacity of the i-th node, m i It is the quality of the i-th node. It is the rate of temperature change at the i-th node.

[0232] Step 8: Solve the heat balance equation using the Runge-Kutta method to calculate the temperature response over time Δt; import the latest temperature into the kinetic model.

[0233] Step 9: Consider the effects of temperature-related parameters and update the grease viscosity and bearing thermal deformation.

[0234] Step 10: Determine if the temperature change ΔT at the hot node is greater than 10. -3 If the value is greater than 10, repeat steps 1 to 9 until the node temperature change ΔT is less than 10. -3 The loop ends, and the final results are the principal shaft displacement field and temperature field responses. For example... Figure 11The process of thermal vibration coupling iteration is introduced.

[0235] This embodiment innovatively constructs a closed-loop energy transfer mechanism using a two-way coupling mechanism: step 5 reflects the force-thermal coupling of vibration to the temperature field (involving nonlinear effects such as frictional heat generation), and step 6 realizes the thermoelastic coupling of temperature to vibration (including thermal expansion effects). Compared with traditional co-simulation, the breakthrough of the algorithm in this embodiment lies in: establishing an improved thermal balance equation that includes nonlinear vibration characteristics, developing a physical bridge for the transfer of displacement-temperature dual-field parameters, and proposing a convergence criterion based on energy conservation. This method inherits and develops the multiphysics coupling modeling theory. Its core lies not in a general iterative framework, but in a customized equation system constructed for the special working conditions of torque motor spindle systems, providing a new theoretical tool for thermal error compensation of precision spindle systems.

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

Claims

1. A method for coupled thermal vibration modeling of a turntable bearing and a torque motor, characterized in that: The method described above studies the thermal vibration characteristics of machine tools during the cutting process, deriving the displacement response of each degree of freedom of the spindle and the temperature field of the bearing and stator / rotor structure. For the thermal vibration coupling analysis, a partitioned iterative algorithm is used to achieve bidirectional coupling modeling. A closed-loop system of thermal vibration coupling is constructed through the coordinated solution of the dynamic equilibrium equation and the thermal equilibrium equation. Specifically, the method includes the following steps: Step 1: Considering the initial conditions and the influence of bearing thermal deformation, derive the bearing restoring force, torque, and elastic stiffness; Step 2: Give the dynamic equilibrium equations of the main shaft system; Step 3: Calculate the unbalanced magnetic pull and the unbalanced mass excitation; Step 4: Solve the dynamic equations using the Newmark-β method to calculate the dynamic response during the time interval Δt, and solve for the displacement field; then transfer the load changes caused by the displacement field to the thermal model; Step 5: Consider the displacement-related load effects, update the contact force on the bearing, calculate the frictional heat generated by the bearing, and the heat generated by various parts of the motor; Step 6: Calculate the bearing's thermal deformation and thermal displacement; Step 7: Calculate the thermal resistance R between each node. ij Furthermore, rapid thermo-mechanical coupling analysis is achieved by deriving the temperature field equilibrium equation through a thermal resistance network. Step 8: Solve the heat balance equation using the Runge-Kutta method to calculate the temperature response over time interval Δt; import the latest temperature into the kinetic model; Step 9: Consider the influence of temperature-related parameters, and update the grease viscosity and bearing thermal deformation; Step 10: Determine if the temperature change ΔT at the hot node is greater than 10. -3 If the value is greater than 10, repeat steps 1 to 9 until the node temperature change ΔT is less than 10. -3 The loop ends, and the final response of the principal shaft displacement field and temperature field is obtained.

2. The method for coupled thermal vibration modeling of a turntable bearing and torque motor according to claim 1, characterized in that: The specific method for step 1 is as follows: Step 1.1: Calculate the position angle; Turntable bearing position angle Calculate using the following formula: (1)； in, These are the bearing position angles of each row of turntables; , ω is the angular velocity of each row of cages in the rotary table bearing, ω is the angular velocity of the spindle rotation, and t is the time variable. It is the pitch circle diameter of the radial rollers. It is the radial roller diameter. This refers to the number of rollers in each row. The superscripts u / d / m represent the parameters corresponding to the upper row, lower row, or radial row of the turntable bearing. Angular contact bearing position angle Calculate using the following formula: (2)； Where, ω c ω is the angular velocity of the angular contact bearing cage, ω is the angular velocity of the main spindle, and Z is the angular velocity of the main spindle. b D represents the number of bearing balls. i D o These are the inner and outer diameters of the bearing, respectively. Step 1.2: Calculate the bearing deformation; The axial variation δ of the center of curvature of the inner raceway groove of the i-th roller in each row of the rotary table bearing a Ru / Rd The radial variation value δ r m and rotational displacement θ u / d As shown in the following formula: (3)； Among them, e R This is the axial distance between the center of the roller and the geometric center of the spindle. It is the pitch circle diameter of the upper and lower rows of rollers; where [x,y,z,θ] x ,θ y These are the dynamic displacements in the x, y, and z directions of the principal axis and the angular displacements in the x and y directions, respectively. For angular contact bearings, with the outer ring curvature center fixed, the radial displacement of the inner ring curvature center under dynamic load... and axial displacement As shown in the following formula: (4)； Among them, e L It is the axial distance between the center of the ball bearings and the geometric center of the spindle. It is the orbital radius of the center of curvature of the inner raceway groove; D m and D b These are the pitch circle diameter and ball diameter of the angular contact bearing, α0 is the initial contact angle, and f is the initial contact angle. i It is the inner raceway curvature radius coefficient; Step 1.3: Calculate the normal approach and contact deformation between the two raceways; The rollers are discretized using a slicing method, and the effective contact length l of each row of rollers is determined. u / d / m Divide into n equal parts k Each slice unit has an axial thickness expressed as: ; Normal approaching amount between the two raceway surfaces at the center of the k-th piece of the i-th roller in the upper row for: (5)； Wherein, the superscript u represents the upper row of rollers, u a0 For the initial axial clearance, u a Axial thermal expansion; Normal approaching amount between the two raceway surfaces at the center of the k-th piece of the i-th roller in the lower row for: (6)； In this context, the superscript d represents the lower row of rollers; The normal approach of the two raceway surfaces at the center of the k-th piece of the i-th roller in the radial arrangement for: (7)； Where the superscript m represents radial rollers, u r0 For the initial radial clearance, u r Radial thermal expansion; For angular contact bearings, considering only the dynamic displacement of the inner ring, the axial distance between the centers of curvature of the inner and outer rings... and radial distance for: (8)； Among them, O i O o It is the distance between the centers of curvature of the inner and outer raceways, α F It is the bearing contact angle after preload; Considering the dynamic displacement and centrifugal effect of the inner ring, as well as the thermal expansion of the bearing, the axial distance between the centers of curvature of the inner and outer rings. and radial distance Determined by the following formula: (9)； Among them, the '+' and '-' of '±' correspond to the first and second angular contact bearings on the left end, respectively; δ jzTD and δ jrTD These are the axial and radial displacements of the inner ring curvature center, taking thermal expansion into account. Based on the geometric compatibility relationship and considering the influence of external dynamic loads, the distance between the centers of curvature of the inner and outer rings on the j-th rolling element is obtained. and dynamic normal contact deformation for: (10)； Considering the bearing expansion caused by temperature rise and centrifugal force, the actual distance A between the centers of curvature is... jT and the actual contact deformation δ at the j-th bearing ball jT With contact angle α jT It is given by the following formula: (11)； Where, δ F It is pre-tightening deformation; Step 1.4: Calculate the bearing's elastic restoring force, torque, and elastic stiffness; The normal load on each roller of the turntable bearing is: (12)； Considering dynamic loads and thermal expansion effects, the nonlinear stiffness characteristics of the bearing are characterized as a multivariable function of dynamic displacement, rotational speed, and temperature field; the nonlinear dynamic restoring force F of the bearing is obtained through the force balance equation of the inner ring. x R F y R F z R and torque M x R M y R for: (13)； According to Hertzian contact theory, the dynamic normal contact force at the j-th ball of an angular contact bearing... for: (14)； in, Indicates Hertzian contact stiffness; Actual normal contact force considering thermal effects for: (15)； Finally, the nonlinear dynamic restoring force F of the angular contact bearing is obtained through the force balance equation of the inner ring. x L F y L F z L and torque M x L M y L : (16)； in, It is the lever arm applied to the inner ring of the bearing; Based on the load-deformation relationship, the nonlinear dynamic stiffness vector of the bearing is derived from the following formula: (17)； Among them, only rigid body displacement is considered, and the elastic body displacement along the principal axis is ignored, [δ x , δ y , δ z , φ x , φ y ] T These are the dynamic displacements of the bearing in the x, y, and z directions and the angular displacements in the x and y rotational directions; then, the bearing stiffness matrix is ​​obtained from equation (17). for: (18)； Among them, [k xL1 , k yL1 , k zL1 , k φxL1 , k φyL1 ]、[k xL2 , k yL2 , k zL2 , k φxL2 , k φyL2 ] and [k xR , k yR ,k zR , k φxR , k φyR These are the displacement stiffness of the two angular contact bearings on the left and the rotary table bearing on the right in the x, y, and z directions, and the angular stiffness in the x and y rotation directions, respectively.

3. The method for coupled thermal vibration modeling of a turntable bearing and a torque motor according to claim 2, characterized in that: The specific dynamic equilibrium equations of the main shaft system in step 2 are as follows: (19)； Where m is the mass of the spindle system. Let the rotational inertia be the diameter of the system. R is the polar moment of inertia of the system. d The extreme diameter; F ij and M ij These represent the restoring force and restoring torque of the left and right bearings, respectively. i = x, y or z, j = 1, 2 or R, where 1, 2, and R represent the two angular contact bearings on the left and the turntable bearing, respectively. , , , , This indicates an unbalanced quality incentive; , , , , This indicates an unbalanced magnetic pull excitation; L1 represents the axial distance between the two bearings, and L2 represents the axial distance from the angular contact bearing to the center. The above formula can be written in matrix form as follows: (20)； Where q=[x,y,z,θ] x ,θ y ] T This represents the dynamic displacement of the principal axis, M is the mass matrix, and C = [c1 c2 c3 c4 c5]. T G is the damping matrix, F is the gyroscope matrix, and F is the damping matrix. b Q is the bearing restoring force matrix, and Q is the external force matrix.

4. The method for coupled thermal vibration modeling of a turntable bearing and torque motor according to claim 3, characterized in that: The specific method for step 3 is as follows: Step 3.1: According to motor winding theory, the fundamental magnetomotive force of the three-phase stator winding... Represented as: (21)； Among them, F sm The fundamental magnetomotive force amplitude of the excitation current for the stator armature reaction current is given by k, where N is the number of turns in series per pole per phase of the stator, and k is the number of turns in series per phase per pole of the stator. w It is the stator winding distribution factor, P is the number of pole pairs of the motor, and I is the stator winding distribution factor. smax It is the maximum value of the stator excitation current, ω cur It is the angular velocity of the rotating magnetic field, and α is the azimuth angle of the air gap space; According to the magnetic circuit principle of permanent magnet motors, the permanent magnet is considered a constant source of magnetomotive force; the fundamental magnetomotive force of a surface permanent magnet motor... Represented as: (22)； Among them, F rm B represents the magnetomotive force amplitude of the permanent magnet. r For remanence of permanent materials, h m denoted as the magnetization thickness of the permanent magnet, and μ0 as the vacuum permeability. Electromagnetic excitation is generated by the air gap magnetic field of the motor, which is determined by the magnetic force between the stator and rotor and the air gap permeability; the fundamental composite magnetomotive force is: (23)； When stationary, the rotor's axial center is the origin of the coordinate system; elastic deformation of the rotor is ignored and it is considered rigid. The initial static eccentricity of the rotor is [x0, y0, z0, θ]. x0 , θ y0 ] T The dynamic displacement of the rotor's axial center relative to the origin is [x, y, z, θ]. x , θ y ] T Based on geometric relationships and considering the influence of angular displacement eccentricity, the dynamic displacement q of any point on the rotor is obtained. r for: (24)； Among them, L r It is the axial length of the rotor; Based on the geometric compatibility relationship, the actual air gap length δ(α,t,Z) related to the rotor's dynamic displacement is derived from the following formula: (25)； in, The air gap length when the stator and rotor coincide; Z is the rotor axial eccentricity, r0 is the rotor static eccentricity, β is the static eccentricity direction angle, and R... ro It is the outer diameter of the rotor; The tilt angle of the rotor surface normal magnetic tensile stress relative to the stator radial direction is expressed by the following formula: (26)； According to Ampere's circuital law, the permeability of the air gap per unit area is obtained by the following formula: (27)； Where, k c It is the groove coefficient; for a closed groove, k c =1; Saturation coefficient k of the magnetic circuit used for magnetic saturation u Consider the magnetic circuit length perpendicular to the rotor surface. for: (28)； The air gap magnetic flux density is then: (29)； According to the magnetic stress formula, the normal magnetic tensile stress on the rotor surface is obtained as follows: (30)； Projecting the magnetic tensile stress along the stator radial and axial directions and integrating it on the rotor surface, the unbalanced electromagnetic force and torque acting on the rotor are obtained as follows: (31)； Step 3.2: Calculate the unbalanced mass excitation; The centrifugal force caused by the eccentric mass projected along the radial and axial directions of the stator yields the unbalanced mass excitation as follows: (32)； Where m is the mass of the spindle system, d e It is the mass eccentricity, and v is the spatial inclination of the centrifugal force relative to the stator radial direction.

5. The method for coupled thermal vibration modeling of a turntable bearing and torque motor according to claim 4, characterized in that: Step 4 specifically includes: The analytical functions of the time-varying unbalanced mass excitation, time-varying unbalanced magnetic pull, nonlinear dynamic bearing stiffness, and bearing restoring force in step 3 are related to the dynamic displacement q of the spindle. The dynamic equations are solved using the Newmark-β method within each time step Δt to obtain the spindle displacement field q. t+Δt The obtained displacement field is then used to update the heating of each component of the spindle system in subsequent steps.

6. The method for coupled thermal vibration modeling of a turntable bearing and a torque motor according to claim 5, characterized in that: Step 5 specifically includes: Step 5.1: Calculate motor heat generation; the internal losses of the motor are divided into stator losses and rotor losses; According to Joule's law, the copper loss P caused by the heating of the conductors in the stator windings... cus Represented as: (33)； Where, m cus =3 indicates a three-phase motor, I represents the effective value of the current through each phase winding, and R represents the effective value of the resistance in each phase winding. The effect of temperature on the winding resistance is ignored. Assuming the magnetic density of the motor core is uniformly distributed during operation, and based on Bertotti's theory of iron loss separation, the stator iron loss P... Fes The losses are categorized into three types: hysteresis loss, eddy current loss, and additional loss. The calculation model expression is as follows: (34)； Among them, P h It is the unit mass hysteresis loss, P e It is the eddy current loss per unit mass, P exc It is the additional loss per unit mass; k h The value represents the hysteresis loss coefficient, f is the frequency, and B is the value of the frequency. m It is the magnetic flux density amplitude; the eddy current loss coefficient is The additional loss coefficient is σ is the electrical conductivity of the core material, d is the thickness of a single layer of silicon steel sheet, ρ is the core density, and G and V0 are performance parameters of the silicon steel sheet, where G is a dimensionless coefficient and V0 is a coefficient related to B. m The relevant statistical characteristic parameters are as follows: S is the cross-sectional area of ​​the silicon steel sheet; The rotor of the permanent magnet synchronous motor adopts a composite structure of permanent magnets and insulating silicon steel sheets. Based on electromagnetic principles, the eddy current loss P of the permanent magnets is... wr The parsing expression is: (35)； Where, k r ρ is the proportionality constant of the electromotive force. w It is the resistivity of permanent magnets. It is the volume of the permanent magnet, where h m L is the thickness of the permanent magnet. a L is the axial length. m Horizontal width; The power loss P between the rotating air gap flow and the stator and rotor friction. w It is derived from the following formula: (36)； in, It is the rotor rotation frequency, where n is the spindle speed, η air It is aerodynamic viscosity; In summary, the stator heating P of the torque motor is calculated. s and rotor heating P r : (37)； Where, m Fe It is the total mass of the stator silicon steel sheets, N w It refers to the number of permanent magnets in the rotor; Step 5.2: Calculate bearing heat generation; Energy loss between rollers and raceways caused by energy dissipation mechanism dominated by material elastic hysteresis effect The following formula is used to derive: (38)； in, It is the distance from the center of the k-th piece of the i-th roller to the center of the bearing. , ξ=1 is the material's elastic hysteresis coefficient, E is the combined elastic coefficient of the two contacting bodies, and ω I / E It refers to the angular velocity of the inner or outer ring. During bearing operation, the difference in radii of gyration at the contact points between the rollers and the inner and outer races causes differential sliding at the contact interface. This relative slippage, induced by geometric kinematics, creates micro-friction on the contact surface. The corresponding energy dissipation mechanism is known in tribology as differential sliding friction. This frictional effect, together with rolling friction, constitutes an important component of the bearing's frictional torque. Its quantitative analysis is based on Hertz contact theory and elastohydrodynamic lubrication model for precise calculation, and is derived from the following formula: (39)； Where u=0.96 is the oil film drag coefficient. It is the relative sliding speed between the upper and lower rows of rollers and the inner raceway. It is the relative sliding speed between the upper and lower rows of rollers and the outer raceway. It is the rotational angular velocity of the upper and lower rows of rollers; for radial rollers, since the linear velocity along the contact line between the roller and the raceway remains constant, the relative sliding between the roller and the raceway is not considered. The grease viscosity effect leads to frictional power loss in the contact area between the rolling elements and the raceway. for: (40)； Where E0 is the equivalent elastic modulus and R0 is the equivalent radius of curvature; These are the speed parameters of the rolling elements and the inner and outer rings, where η0 is the dynamic viscosity of the grease. It is the average speed of the roller between the pitch circle and the inner raceway. It is the angular velocity of the radial rollers' rotation. It is the average speed between the roller and the outer raceway; These are the material parameters of the rolling elements and the inner and outer rings, among which... It is the viscosity-pressure coefficient of the lubricating grease; These are the load parameters for the rollers and the inner and outer rings; In summary, the total heat generated by the inner rings of each row of rollers in the turntable bearing is... Total heat generation in the outer ring They are respectively: (41)； To improve the heat generation of angular contact bearings, thermal expansion, nonlinear dynamic loads, and lubricant viscosity variations are considered. The total frictional torque of the bearing is calculated using the Palmgren formula. for: (42)； Among them, M l and M v These are the frictional torques, M, caused by the applied load and the viscous friction of the lubricant, respectively. ω It is the rotational friction torque of the rolling element; Frictional torque M caused by the applied load on the bearing l The analysis is as follows: the curvature compatibility of the rolling contact area (i.e., the geometric size effect), the elastohydrolubrication state of the Hertzian contact ellipse (i.e., the deformation energy storage), the Coulomb friction work in the micro-slip zone (i.e., the rolling to sliding friction ratio), and the energy dissipation field constituted by the load vector space distribution are calculated using the following formula: (43)； Among them, f l F is a coefficient that takes into account the effects of load and bearing type. l This refers to the load on the bearing. For angular contact bearings, the calculation is obtained from the following formula: (44)； Among them, F s It is the equivalent static load of the bearing. C s This is the basic static load rating; F r and F a These are radial and axial loads, calculated using the following formulas: (45)； Frictional torque M caused by viscous friction of lubricant v Essentially, it is the energy loss generated by the shearing action of the oil film during bearing movement. For angular contact bearings, M v The calculation is as follows: (46)； Among them, f v It is a coefficient that takes into account the bearing structure and lubrication method; , represents the relationship between the viscosity of the grease and temperature, where ν0 represents the kinematic viscosity of the grease at temperature T0, β is the temperature coefficient of the grease, and T is the temperature of the grease; Rotational friction torque M of the rolling element ω This refers to the torque caused by friction between contact interfaces when a rolling element rolls along a raceway or other contact surface. This torque characteristic reflects the mechanical energy loss characteristics of the system due to surface friction during rotation, and the formula is as follows: (47)； Where, μ ω Let a and b be the coefficient of rotational friction between the rolling elements and the inner raceway, respectively, and a and b be the major and minor semi-axises projected onto the inner race surface by the deformation of the rolling elements; then the total frictional heat generated by the bearing is P. b for: (48)； The heat distributed between the inner and outer raceways is converted into heat distributed on the surfaces of the inner and outer raceways, which can be expressed by the following formula: (49)； Frictional loss P between the ball and the cage c We obtain it from the following formula: (50)； Where, m c It is the cage mass, μ c It is the coefficient of sliding friction between the ball and the cage.

7. The method for coupled thermal vibration modeling of a turntable bearing and torque motor according to claim 6, characterized in that: Step 6 specifically includes: Step 6.1: Calculate the radial thermal expansion of the turntable bearing; Assuming the material is isotropic and its thermal expansion is linearly distributed, the radial expansion u of the outer ring can be derived from the thermal expansion formula of a hollow cylinder. o and inner radial expansion u i and radial roller expansion for: (51)； in, , and These are the coefficients of thermal expansion of the bearing's inner ring, rollers, and outer ring, respectively, ΔT. i and These are the temperature rises of the bearing inner ring and bearing rollers, respectively, ΔT o and ΔT h These are the temperature rises of the bearing outer ring and the bearing housing, respectively, d i and d o These are the inner and outer ring diameters of the turntable bearing, d h is the inner diameter of the bearing housing, and μ is Poisson's ratio; According to equation (51), the radial thermal expansion u of the bearing is obtained. r for: (52)； Step 6.2: Calculate the axial thermal expansion of the turntable bearing; Axial expansion of bearing inner ring a i axial expansion of the outer ring a o and the axial expansion of the upper and lower rows of rollers a b u(d) Use the following formula: (53)； Among them, l i and l o These represent the widths of the inner and outer rings of the bearing, respectively. It is the diameter of the upper and lower rows of rollers. It is the temperature rise of the upper and lower rollers; According to equation (53), the axial thermal expansion u of the bearing is obtained. a for: (54)； Step 6.3: The thermal deformation of the balls in an angular contact bearing is derived from the following formula: (55)； Among them, D b It is the diameter of the bearing balls, ΔT b The ball bearings experience uniform temperature rise. Considering the thermal expansion effect of the shaft and the inner ring of the bearing, as well as the interference fit between the inner ring and the journal, the inner and outer rings of the bearing are hollow thin-walled circular ring structures. The radial thermal deformation U of the inner ring raceway is... ir The calculation is performed using the following formula: (56)； Among them, D i D is the inner diameter of the bearing. io U is the diameter of the bearing's inner raceway; io U i U ro The thermal deformations of the inner ring raceway, bearing inner diameter, and journal diameter are respectively calculated using the following formulas: (57)； Among them, D ro The journal diameter is... μ is the coefficient of thermal expansion of the shaft. r For the axial Poisson's ratio, ΔT ro For uniform temperature rise of the journal; Considering the thermal effects of the bearing housing and outer ring, as well as the clearance fit, the radial thermal deformation U of the outer ring raceway... oh Calculate using the following formula: (58)； Among them, D o D is the outer diameter of the bearing. oi U is the outer raceway diameter, δ1 is the clearance fit between the outer ring and the bearing housing; oi U o U hi The thermal deformations of the outer raceway, bearing outer diameter, and bearing housing inner diameter are calculated using the following formulas: (59)； Among them, D hi The inner diameter of the bearing housing. μ is the coefficient of thermal expansion of the bearing housing. h For the bearing housing Poisson's ratio, ΔT hi To ensure uniform temperature rise in the bearing housing; Step 6.4: The thermal displacement problem caused by the thermal expansion of the angular contact bearing is equivalent to a static load-displacement problem. Based on the geometric coordinate relationship, the radial displacement of the inner ring raceway curvature center caused by the thermal deformation of the bearing balls is obtained. and axial displacement They are respectively: (60)； in, It is the bearing contact angle; Radial displacement of the inner raceway curvature center caused by radial thermal deformation of the outer raceway and axial displacement They are obtained from the following formulas: (61)； Considering the expansion of the inner ring due to centrifugal force, the radial displacement of the inner ring raceway curvature center caused by radial thermal deformation of the inner ring raceway is derived by the following formula: (62)； Among them, U cent For the centrifugal expansion of the inner ring, These are the bearing pitch circle diameters, ρ, E, and μ. i These are the material density, elastic modulus, and Poisson's ratio of the inner ring, respectively. Angular contact bearings are arranged back-to-back. In thermal displacement analysis, the center distance between the inner and outer ring curvatures of the left and right bearings shows opposite trends. Based on the assumption that the outer ring curvature center is fixed, the axial and radial displacement analysis of the inner ring curvature center is established by coupling the bearing thermal expansion effect with the centrifugal force of the inner ring, as shown in the following equation: (63)； in,' The '-' character corresponds to the first angular contact bearing on the left, and the '+' character corresponds to the second angular contact bearing on the left.

8. The method for coupled thermal vibration modeling of a turntable bearing and torque motor according to claim 7, characterized in that: Step 7 specifically includes: Step 7.1: Calculate thermal resistance; Each node transfers heat through conduction, convection, and radiation. Neglecting the effect of radiation, calculate the axial thermal resistance R between each node. a Radial thermal resistance R r and convection thermal resistance R v As shown in the following formula: (64)； Among them, R out and R in These are the outer radius and inner radius of the hollow cylinder, respectively, K. D A is the thermal conductivity of the solid, A is the convective heat transfer area, and h is the convective heat transfer coefficient, which is estimated by the Nusselt number N. u Estimate; L is the characteristic length of heat exchange, K l The thermal conductivity of the liquid; Step 7.2: Construct a thermal resistance network, derive the temperature field equilibrium equation, and achieve rapid thermal vibration coupling analysis; A 71-node thermal network model was constructed to analyze the temperature field of the torque motor spindle system, including 15 key heating nodes and 59 nodes whose temperatures need to be solved. A rapid thermo-mechanical coupling analysis was achieved through a thermal resistance network. The calculated thermal deformation of the components was imported into the dynamic model to study the thermally induced vibration response and its coupling with the system's heating. The thermal resistance network includes three core thermal resistances: thermal convection resistance, which characterizes the heat transfer efficiency between the solid surface and the cooling medium; thermal conduction resistance, which reflects the thermal conductivity of the solid material; and thermal contact resistance, which describes the heat transfer mechanism of the contact interface and gap. The model makes the following assumptions: the heat conduction unit is characterized by the central node to represent the average temperature field, the surface temperature of the heat convection unit is approximately equal to the central temperature, the contact surface temperature is dominated by the central node, the influence of heat radiation and the temperature rise of the cooling medium is ignored, and the ambient temperature is kept constant as the boundary condition. The thermal balance equations for each node in the angular contact bearing are as follows: （64）； The heat balance equations for each node in the motor are as follows: （65）； The heat balance equations for each node of the turntable bearing are as follows: （66）； in, T i T is the temperature of the i-th node. j R is the temperature of the j-th node adjacent to the i-th node; i_j P is the total thermal resistance between the i-th node and the j-th node; i It is the heat generated per unit time at the i-th node, including P. i1 P c1 P o1 P i2 P c2 P o2 These are the heat generated by the inner ring, the frictional heat of the cage, and the heat generated by the outer ring of the two angular contact bearings, respectively. r P s P w These are the friction losses of the torque motor rotor, stator, and airflow, respectively, P I u P E u P I d P E d P I m P E m These refer to the heating of the inner and outer rings of the upper and lower rows and the radial row of the turntable bearings; c i It is the specific heat capacity of the i-th node, m i It is the quality of the i-th node. It is the rate of temperature change at the i-th node.

Images