Multi-physical field coupling simulation method for motorized spindle

Abstract

The application belongs to the technical field of simulation analysis, and particularly relates to a multi-physical field coupling simulation method of an electric spindle. The multi-physical field coupling simulation method of the electric spindle comprises the following steps: S100, a multi-physical field coupling model is established through finite element software, including: importing a geometric model, setting a material, setting a rotor dynamics module, setting a solid heat transfer module and mesh division; S200, solving is performed through the finite element software, including: performing multi-physical field coupling of a bearing nonlinear restoring force model, a motor unbalanced magnetic pull model and a spindle thermal model with the electric spindle rotor model, and obtaining dynamic performance by solving. The simulation model for the dynamic characteristic research of the electric spindle in the technical scheme can meet the demand of the dynamic characteristic research of the electric spindle under actual complex working conditions.

Description

TECHNICAL FIELD

[0001] The application belongs to the technical field of simulation analysis, and particularly relates to a multi-physical field coupling simulation method for an electric spindle. BACKGROUND

[0002] The simulation models for the dynamic characteristics of the electric spindle in the prior art mainly include a bearing-rotor-thermal coupling model, a bearing-rotor-motor (electromagnetic) coupling model, and a motor (electromagnetic)-rotor-thermal coupling model. In a complex working condition, the electric spindle has strong nonlinear characteristics during operation. The simulation models for the dynamic characteristics of the electric spindle in the prior art cannot meet the research needs of the dynamic characteristics of the rotor in the actual complex working condition.

[0003] Therefore, it is urgent to provide a multi-physical field coupling simulation method for an electric spindle that is closer to the actual complex working condition to solve the above problems. SUMMARY

[0004] The application aims to provide a multi-physical field coupling simulation method for an electric spindle that can realize the research of the dynamic characteristics in a complex working condition. The application achieves the goal by the following technical scheme:

[0005] The first aspect of the application provides a multi-physical field coupling simulation method for an electric spindle, including the following steps:

[0006] S100, a rotor model of an electric spindle is established by using finite element software, including importing a geometric model, setting a material, setting a rotor dynamics module, setting a solid heat transfer module, and meshing, wherein the geometric model includes a rotor structure model;

[0007] The rotor dynamics module is set, including applying a nonlinear restoring force and an unbalanced magnetic pull to the rotor structure model;

[0008] The solid heat transfer module is set, including applying a motor-acting rotor surface temperature T motor and a bearing-acting rotor surface temperature T bearing to the surface of the rotor structure model;

[0009] The nonlinear restoring force is solved by a bearing nonlinear restoring force model, the unbalanced magnetic pull is solved by a motor unbalanced magnetic pull model, and the motor-acting rotor surface temperature T motor and the bearing-acting rotor surface temperature T bearing are solved by a spindle thermal model;

[0010] S200, solving by the finite element software, comprising: coupling the bearing nonlinear restoring force model, the motor unbalanced magnetic pull model and the spindle thermal model with the motorized spindle rotor model in multiple physical fields, and solving to obtain dynamic performance.

[0011] The motorized spindle simulation model in the technical solution is a multi-physical field coupling simulation model, that is, a bearing-motor-rotor-spindle thermal coupling dynamic model is integrated at the same time, the model can more accurately reflect the spindle time domain rotation accuracy, critical speed and order distribution of the spindle in actual complex working conditions, and provide a theoretical basis for subsequent high-performance spindle optimization design.

[0012] In addition, the motorized spindle multi-physical field coupling simulation method of the application can have the following additional technical features:

[0013] In some embodiments of the application, the geometric model further comprises a tool head structure model, a tool shank structure model and a motor rotor structure model, the tool head structure model and the rotor structure model are connected to two ends of the tool shank structure model respectively, and the motor rotor structure model is sleeved on the outer periphery of the rotor structure model.

[0014] In some embodiments of the application, the rotor structure model has a cavity through the axial direction.

[0015] In some embodiments of the application, in S100, the unbalanced magnetic pull is applied to the rotor structure model, comprising the following steps:

[0016] S111, dividing the motor rotor structure model into a plurality of shaft segments along the axial direction thereof;

[0017] S112, establishing a coordinate system with the center point of the end face away from the motor on each shaft segment as the origin, and acquiring the position information of the selected node in real time by a probe built in the finite element software;

[0018] S113, taking the position information as an input variable of the motor unbalanced magnetic pull model;

[0019] S114, applying each unbalanced magnetic pull to the corresponding shaft segment.

[0020] In some embodiments of the application, in S100, the bearing nonlinear restoring force model solving comprises the following steps:

[0021] S121, setting bearing parameters and convergence accuracy, the bearing parameters comprising a contact angle of a ball and an inner and outer raceway , a ball diameter ; and obtaining the relative radial displacement of the bearing inner and outer rings in the X-axis direction The relative radial displacement of the inner and outer rings of the bearing in the Y-axis direction The relative axial displacement of the inner and outer rings of the bearing The relative angular displacement of the inner and outer rings of the bearing around the X-axis The relative angular displacement of the inner and outer rings of the bearing around the Y-axis direction The X-axis and the Y-axis are perpendicular

[0022] S122, calculating the curvature center of the outer raceway groove of the single ball O The axial distance between the final position of the curvature center of the inner raceway groove The radial distance between the final position of the curvature center of the inner raceway groove The curvature center of the outer raceway groove of the single ball O The final position of the curvature center of the inner raceway groove The radial distance between the final position of the curvature center of the inner raceway groove ;

[0023]

[0024] wherein, is the distance between the curvature centers of the inner and outer raceway grooves , , are the groove curvature radius coefficients of the inner and outer raceways, respectively and satisfy the following relationship:

[0025]

[0026] In the formula, is the azimuth angle of the ball is the radius of the curvature center trajectory of the inner raceway, which satisfies the following relationship:

[0027]

[0028] In the formula, is the pitch diameter

[0029] The constraint , is the initial value and is substituted into the following formula:

[0030]

[0031] In the formula, , are the axial and radial projections of the distance between the final position of the ball center and the curvature center of the outer raceway, respectively

[0032] According to the Hertz contact theory, the contact deformation of the ball and the inner ring raceway and the contact deformation of the ball and the outer ring raceway can be calculated by the following formula:

[0033]

[0034] where represents the normal pressure of the contact deformation, , , , , are calculated by the following equations, respectively:

[0035]

[0036] where , is the elastic modulus of the two contacts; , is the Poisson's ratio of the two contacts;

[0037] and are the radii of curvature of the two contacts in the two principal planes, respectively;

[0038] Ball and inner raceway contact:

[0039]

[0040] Ball and outer raceway contact:

[0041]

[0042] where ;

[0043] Since the rolling bearing is lubricated, the effect of the lubricating oil film thickness should also be considered:

[0044]

[0045] where ， are the distances between the centers of curvature of the inner and outer raceway grooves and the center of the ball, respectively; , are the central oil film thicknesses between the ball and the inner and outer raceway grooves, respectively, and can be calculated by the following equations:

[0046]

[0047] where ; e is a natural constant; other unknown quantities are determined by the following equations:

[0048]

[0049] According to the position relationship between the curvature center of the inner and outer raceway groove and the ball center, the following geometric relationships exist:

[0050]

[0051] wherein, , are the contact angles between the ball and the inner and outer raceway respectively, wherein is the ball number, satisfying , is the total number of balls;

[0052] S123, the force relationship between a single ball and the inner and outer raceway under the action of centrifugal force and gyroscopic moment is established and solved:

[0053]

[0054] wherein, , are the friction coefficients between the ball and the inner and outer raceway respectively, , are the contact angles between the ball and the inner and outer raceway respectively, is the ball diameter, , represents the normal pressure of the contact deformation between the ball and the inner and outer raceway; centrifugal force and gyroscopic moment are determined by the following formula:

[0055]

[0056] wherein, is the rotation angular velocity of the ball, is the revolution angular velocity of the ball, is the rotational inertia of the ball, is the bearing inner ring rotation speed, is the included angle between the rotation axis and the revolution axis, is the mass of the ball;

[0057] According to the raceway control hypothesis of Jones, i.e., the bearing is closer to the outer raceway control when operating at high speed. The following relationship can be established:

[0058]

[0059] wherein, , ;

[0060] The bearing inner ring is simultaneously subjected to the force and moment of the rotation shaft and the ball. The force analysis is performed to obtain the non-linear restoring force calculation formula of a single ball:

[0061]

[0062] wherein, ;

[0063] S124, repeat steps S2-S3 until all balls are calculated;

[0064] S125, summing up to obtain the non-linear restoring force.

[0065] In some embodiments of the present application, in S100, the motor unbalanced magnetic pull model solving comprises the following steps:

[0066] S131, setting motor parameters, the motor parameters including initial air gap length δ, stator slot number Z;

[0067] S132, calculating eccentric air gap length :

[0068]

[0069] S133, calculating air gap permeance :

[0070]

[0071] wherein, as follows:

[0072]

[0073] S134, calculating motor magnetic potential ;

[0074]

[0075] wherein, is the stator winding fundamental magnetic potential, is the electrical frequency; is the current maximum value; is the number of turns in series per phase winding; is the number of motor pole pairs; is the power factor angle; is the stator winding coefficient, expressed as follows:

[0076]

[0077]

[0078] wherein, is the number of slots per pole per phase; are the slot pitch angles related to the pole number respectively, expressed as follows:

[0079]

[0080] in, This refers to the number of stator slots;

[0081]

[0082] In the formula, The winding pitch ratio;

[0083] S135, Calculate the air gap magnetic flux density :

[0084]

[0085] S136. Calculate Maxwell stress σ :

[0086]

[0087] S137. Calculate the unbalanced magnetic pull:

[0088] ;

[0089] ;

[0090] in, , , , The following relationship must be satisfied:

[0091] .

[0092] In some embodiments of the present invention, solving the spindle thermal model includes solving the bearing thermal model, the motor thermal model, and the heat transfer model.

[0093] In some embodiments of the present invention, solving the bearing thermal model includes the following steps:

[0094] Calculate the frictional heat generated by the ball and the inner and outer raceways. :

[0095]

[0096] In the formula, where It is the coefficient of rotational friction between the ball and the inner raceway. For ball bearings, ; It is the angular velocity of the ball's spin; It is contact resilience; ;

[0097] Calculate the heat generated by friction between the ball and the cage. :

[0098]

[0099] In the formula, Sphere diameter; The bearing pitch circle is equal to the average of the bearing's inner and outer diameters. ; Initial junction angle; To maintain rack quality; The coefficient of friction between the ball and the cage; To maintain the frame angular velocity, ;

[0100] Calculate the total heat generated by the bearing :

[0101] .

[0102] In some embodiments of the present invention, solving the motor thermal model includes the following steps:

[0103] Calculate the copper loss of an asynchronous motor :

[0104]

[0105] Stator copper loss and rotor copper loss Calculated by the following formula:

[0106]

[0107]

[0108] In the formula, This refers to the stator winding current. For winding resistance; This refers to the rotor winding current; For rotor bar resistance;

[0109] Calculate the iron loss of an asynchronous motor :

[0110]

[0111] Stator or rotor yoke losses Calculated by the following formula:

[0112]

[0113] In the formula This is an empirical coefficient, when the capacity hour, Otherwise, it is 1.3; The weight of the yoke; is the yoke loss coefficient, and the expression is as follows:

[0114]

[0115] where is the loss per unit weight of silicon steel sheet, with the unit of W / kg; is the magnetic induction of the idle yoke;

[0116] Stator or rotor tooth iron loss calculated by the following formula:

[0117]

[0118] where is the weight of the yoke; is an empirical coefficient, and is taken as 1.8 for an induction motor; is the tooth loss coefficient;

[0119]

[0120] where is the loss per unit weight of silicon steel sheet, with the unit of W / kg; is the magnetic induction of the idle tooth;

[0121] Rotor and air friction loss

[0122]

[0123] where is the rotation frequency; is the rotor outer diameter; is the rotor axial length; is the air dynamic viscosity; is the average air gap length;

[0124] Asynchronous motor stray loss

[0125]

[0126] where the copper bar is taken as 0.005, is the output power;

[0127] Total heat generation power

[0128] .

[0129] ​​​​​In some embodiments of the present application, the solving of the heat transfer model comprises the following steps:

[0130] The convective heat transfer amount of the stator or rotor core to the air gap :

[0131]

[0132] wherein, is the Stefan-Boltzmann constant; is the relative emissivity, usually taken as 0.85; is the heat outflow surface temperature; is the heat absorption surface temperature; is the heat sink surface area;

[0133] The surface temperature rise of the stator or rotor core T is calculated by the following formula:

[0134]

[0135] wherein, is the stator or rotor iron loss; is the stator or rotor heat dissipation area, is the heat transfer coefficient per unit surface area to the external air;

[0136] The convective heat transfer amount of the stator or rotor core to the air gap :

[0137]

[0138] The heat transfer coefficient of the stator or rotor core to the air gap is calculated by the following formula:

[0139]

[0140] wherein, is the air thermal conductivity; the air thermal conductivity is 0.0267 W / m-K at a normal temperature of 20°C; the air thermal conductivity is 0.0251 W / m-K at 0°C; the air thermal conductivity is 0.0321 W / m-K at 100°C;

[0141] Nusselt number Nu is calculated by the following formula:

[0142]

[0143] wherein, ;

[0144] Taylor number Ta is calculated by the following formula:

[0145]

[0146] In the formula, The mass density of the fluid; This refers to the rotor angular velocity; The average radius of the stator and rotor; This is the air gap length; For fluid dynamic viscosity;

[0147] Calculate the convective heat transfer from the rotor core ends to the air. P 4:

[0148]

[0149] Natural convection heat transfer coefficient Calculated by the following formula:

[0150]

[0151] In the formula, Temperature of the surface from which heat flows out; Temperature of the heat absorption surface; The outer diameter of the rotor;

[0152] Calculate the convective heat transfer from the stator core to the coolant. P 5:

[0153]

[0154] Stator core heat transfer coefficient to coolant Calculated by the following formula:

[0155]

[0156] In the formula, For coolant velocity; This refers to the length of the motor body;

[0157] Calculate the convective heat transfer from the stator core to the air. :

[0158]

[0159] Natural convection heat transfer coefficient Calculated by the following formula:

[0160]

[0161] In the formula, The temperature of the radiant surface; Temperature of the absorption surface; Stator outer diameter;

[0162] Calculating the amount of heat radiated by the stator core to the air :

[0163]

[0164] wherein, is the Stefan-Boltzmann constant; is the relative emissivity, usually taken as 0.85; is the radiating surface temperature; is the absorbing surface temperature; is the radiating body surface area;

[0165] Calculating the heat dissipated by the bearing :

[0166]

[0167] Calculating the residual heat W :

[0168]

[0169] Calculating the surface temperature of the shaft acted upon by the motor :

[0170]

[0171] wherein, is the surface area of the shaft acted upon by the motor heat generation temperature;

[0172] Calculating the surface temperature of the shaft acted upon by the bearing :

[0173]

[0174] wherein, is the surface area of the shaft acted upon by the bearing heat generation temperature. BRIEF DESCRIPTION OF DRAWINGS

[0175] In the drawings:

[0176] Figure 1 a structural schematic diagram of a geometric model according to an embodiment of the present application is shown schematically;

[0177] Figure 2 a sectional view of a shaft structural model according to an embodiment of the present application is shown schematically;

[0178] Figure 3 a flowchart of the solution of a multi-physical field coupling model of an electric spindle according to an embodiment of the present application is shown schematically;

[0179] Figure 4The position relationship between the curvature center of the inner and outer raceway groove and the ball center of the bearing according to the embodiment of the application is schematically shown;

[0180] Figure 5 The schematic diagram of the force on the bearing ball according to the embodiment of the application is schematically shown;

[0181] Figure 6 The motor stator-rotor double-cylinder model considering eccentricity according to the embodiment of the application is schematically shown;

[0182] Figure 7 The flowchart of solving the bearing thermal model according to the embodiment of the application is schematically shown;

[0183] Figure 8 The flowchart of solving the motor thermal model according to the embodiment of the application is schematically shown;

[0184] Figure 9 The curve diagram of the nonlinear restoring force of the front bearing varying with time according to the embodiment of the application is schematically shown;

[0185] Figure 10 The curve diagram of the nonlinear restoring force of the rear bearing varying with time according to the embodiment of the application is schematically shown;

[0186] Figure 11 The curve diagram of the unbalanced magnetic pull varying with time according to the embodiment of the application is schematically shown;

[0187] Figure 12 The curve diagram of the temperature of the front bearing varying with time according to the embodiment of the application is schematically shown;

[0188] Figure 13 The curve diagram of the temperature of the rear bearing varying with time according to the embodiment of the application is schematically shown;

[0189] Figure 14 The curve diagram of the temperature of the motor varying with time according to the embodiment of the application is schematically shown;

[0190] Figure 15 The curve diagram of the radial displacement of the rotating shaft varying with time according to the embodiment of the application is schematically shown;

[0191] Figure 16 The Campbell diagram of the dynamic performance analysis result of Example One is schematically shown;

[0192] Figure 17 The trajectory diagram of the dynamic performance analysis result of Example One is schematically shown;

[0193] Figure 18A spectrum diagram of a dynamic performance analysis result of the embodiment one is schematically shown.

[0194] The reference signs in the drawings represent the following:

[0195] 100, tool head structure model; 200, tool holder structure model; 300, shaft structure model; 301, cavity; 400, motor rotor structure model. DETAILED DESCRIPTION

[0196] Figure 1 A structural schematic diagram of a geometric model according to an embodiment of the present application is schematically shown. Figure 2 A structural schematic diagram of a shaft structure model according to an embodiment of the present application is schematically shown. Figure 3 A flowchart of a solution of a motor spindle multi-physical field coupling model according to an embodiment of the present application is schematically shown. Figures 1 to 3 As shown in the figure, the present application proposes a motor spindle multi-physical field coupling simulation method, including the following steps:

[0197] S100, establishing a motor spindle rotor model through finite element software, including: importing a geometric model, setting a material, setting a rotor dynamics module, setting a solid heat transfer module and meshing, the geometric model including a shaft structure model 300;

[0198] Setting the rotor dynamics module includes applying a nonlinear restoring force (NRF) and an unbalanced magnetic pulling force (UMP) to the shaft structure model 300;

[0199] Setting the solid heat transfer module includes applying a motor acting shaft surface temperature T motor and a bearing acting shaft surface temperature T bearing to the surface of the shaft structure model 300;

[0200] Solving the nonlinear restoring force through a bearing nonlinear restoring force model, solving the unbalanced magnetic pulling force through a motor unbalanced magnetic pulling force model, and solving the motor acting shaft surface temperature T motor and the bearing acting shaft surface temperature T bearing through a spindle thermal model;

[0201] S200, solving through finite element software, including: coupling the bearing nonlinear restoring force model, the motor unbalanced magnetic pulling force model and the spindle thermal model with the motor spindle rotor model in a multi-physical field, and solving to obtain dynamic performance.

[0202] It can be understood that, in the finite element software, the rotating shaft structure model 300 is equivalent to the rotating shaft in the actual working condition. As shown in Figure 3 , the non-linear restoring force is calculated by the bearing non-linear restoring force model and is applied to the front and rear rotating shafts. At the same time, after the rotating shaft rotates to generate eccentricity, the bearing node is acted on, thereby affecting the bearing non-linear restoring force model. Similarly, the unbalanced magnetic pull model is calculated by the motor unbalanced magnetic pull model, and is applied to the rotating shaft. At the same time, after the rotating shaft rotates to generate eccentricity, the motor unbalanced magnetic pull model is acted on. As shown in Figure 4 and Figure 5 , the motor temperature T motor and the bearing temperature T bearing are applied to the surface of the rotating shaft. The motorized spindle simulation model in the technical solution is a multi-physical field coupling simulation model, that is, a bearing-motor-rotating shaft-main shaft thermal coupling dynamic model is integrated, which can more accurately reflect the main shaft time domain rotation accuracy, critical speed and order distribution of the main shaft in the actual complex working condition, and provide a theoretical basis for subsequent high-performance main shaft optimization design.

[0203] Further, referring to Figure 1 , the geometric model further includes a tool head structure model 100, a tool shank structure model 200, and a motor rotor structure model 400, the tool head structure model 100 and the rotating shaft structure model 300 are respectively connected to two ends of the tool shank structure model 200, and the motor rotor structure model 400 is sleeved on the outer periphery of the rotating shaft structure model 300.

[0204] In the prior art, the geometric model is only a rotating shaft model without other additional structures, which cannot meet the research needs of the dynamic characteristics of the rotating shaft in the actual complex working condition. The geometric model in the technical solution simultaneously includes a tool shank, a tool head and a motor rotor, and can accurately reflect the dynamic performance of the main shaft in the actual complex working condition. Optionally, the tool head can be a standard bar or other cutters.

[0205] Further, referring to Figure 2 , the rotating shaft structure model 300 has a cavity 301 penetrating in the axial direction.

[0206] In the prior art, the geometric model is usually a solid rotating shaft model, but in actual situations, the rotating shaft is a hollow structure. Therefore, the technical solution sets the cavity 301 penetrating in the axial direction in the rotating shaft structure model 300, which is closer to the real working condition, so that more accurate analysis results can be obtained.

[0207] Further, continuing to refer to Figure 1 , in S100, the unbalanced magnetic pull is applied to the rotating shaft structure model 300, including the following steps:

[0208] S111. Divide the motor rotor structure model 400 into multiple axial segments along its axis evenly;

[0209] S112. Establish a coordinate system with the center point of the end face far from the motor on each axial segment as the origin, and obtain the position information of the selected nodes in real time through the probe built into the finite element software;

[0210] S113. Use the position information as the input variable of the motor unbalanced magnetic pull force model;

[0211] S114. Apply each unbalanced magnetic pull force to the corresponding axial segment.

[0212] Among them, the probe is a built-in function of the finite element software, which can be used to obtain the position information of the selected nodes of the finite element model in real time.

[0213] During the actual operation of the rotating shaft, its radial eccentricity will vary due to different axial positions. However, in the existing technology, the unbalanced magnetic pull force modeling method only applies a total unbalanced magnetic pull force on the overall structure of the motor rotor, making it difficult to reflect the actual distribution of the unbalanced magnetic pull force at different axial positions along the rotating shaft. In this technical solution, by dividing the motor rotor part into multiple parts along the axis and applying the unbalanced magnetic pull forces to the corresponding motor rotor axial segments respectively, the difference in the unbalanced magnetic pull force along the axis of the rotating shaft under different eccentricities can be reflected, so as to accurately reflect the dynamic performance of the rotating shaft under actual complex working conditions. Exemplarily, in this technical solution, the motor rotor part is divided into four parts, and the length of each part along the axis of the geometric model is equal. In other embodiments, it can be divided into three parts, five parts, six parts, etc. according to needs, and the length of each part can be equal or unequal.

[0214] Further, in this technical solution, the finite element software used in the model establishment and solution process is COMSOL Multiphysics.

[0215] COMSOL Multiphysics is a leading multi-physics simulation platform in the industry, providing functions for simulating single physics fields and flexibly coupling multiple physics fields, and can be used to accurately analyze equipment, processes and flows in various engineering fields. The built-in model developer of the software contains a complete modeling workflow, which can implement all simulation steps from geometric modeling, material parameter and physics field setting, solution to result processing.

[0216] Further, as Figure 3 shown, in S100, the solution of the bearing nonlinear restoring force model includes the following steps:

[0217] S121. Set the bearing parameters and convergence accuracy. The bearing parameters include the contact angle between the ball and the inner and outer raceways Diameter of ball Obtain relative radial displacement of inner and outer rings of bearing in X-axis direction Obtain relative radial displacement of inner and outer rings of bearing in Y-axis direction Obtain relative axial displacement of inner and outer rings of bearing Obtain relative angular displacement of inner and outer rings of bearing around X-axis Obtain relative angular displacement of inner and outer rings of bearing around Y-axis direction X-axis and Y-axis are perpendicular

[0218] S122, as shown in Figure 4 Calculate axial distance between outer raceway groove curvature center of single ball and final position of inner raceway groove curvature center O Calculate radial distance between outer raceway groove curvature center of single ball and final position of inner raceway groove curvature center O

[0219]

[0220] Wherein, Distance between inner and outer raceway groove curvature centers , , Respectively, the radius of curvature coefficient of inner and outer raceway groove And Satisfy the following relationship:

[0221]

[0222] In the formula, Azimuth angle of ball Radius of curvature center track of inner raceway, satisfying the following relationship:

[0223]

[0224] In the formula, Pitch diameter

[0225] Constraint ,Initial value of and substitute into the following formula:

[0226]

[0227] In the formula, ,Respectively, the distance between the final position of the ball center and the outer raceway curvature center in the axial and radial projection

[0228] ​​​​​According to the Hertz contact theory, the contact deformation of the ball and the inner ring raceway and the contact deformation of the ball and the outer ring raceway can be calculated by the following equations:

[0229]

[0230] wherein denotes the normal pressure of the contact deformation, , , , , are calculated by the following equations, respectively:

[0231]

[0232] wherein , is the modulus of elasticity of the two contact bodies; , is the Poisson's ratio of the two contact bodies;

[0233] and are the radii of curvature of the two contact bodies in the two main planes, respectively;

[0234] Contact of the ball and the inner raceway:

[0235]

[0236] Contact of the ball and the outer raceway:

[0237]

[0238] wherein ;

[0239] Since the rolling bearing is lubricated, the influence of the lubricating oil film thickness should also be considered:

[0240]

[0241] wherein ， are the distances of the centers of curvature of the inner and outer raceway grooves and the center of the ball, respectively; , are the central oil film thicknesses between the ball and the inner and outer ring raceways, respectively, and can be calculated by the following equations:

[0242]

[0243] wherein ; e is a natural constant; other unknown quantities are determined by the following equations:

[0244]

[0245] As Figure 4 shown, according to the position relationship between the curvature center of the inner and outer raceway groove and the ball center, there are also the following geometric relationships:

[0246]

[0247] In the formula, , are the contact angles of the ball and the inner and outer raceway, respectively, wherein is the ball number, satisfying , is the total number of balls;

[0248] S123, as Figure 5 shown, the force relationship between a single ball and the inner and outer raceway under the action of centrifugal force and gyroscopic moment is established and solved:

[0249]

[0250] wherein, , are the friction coefficients between the ball and the inner and outer raceway, respectively, , are the contact angles of the ball and the inner and outer raceway, respectively, is the ball diameter, , represent the normal pressure of the contact deformation of the ball and the inner and outer raceway; the centrifugal force and the gyroscopic moment are determined by the following formula:

[0251]

[0252] wherein, is the rotation angular velocity of the ball, is the revolution angular velocity of the ball, is the moment of inertia of the ball, is the bearing inner ring rotation speed, is the included angle between the rotation axis and the revolution axis, is the mass of the ball;

[0253] According to the raceway control hypothesis of Jones, that is, the bearing is closer to the outer raceway control when operating at high speed, the following relationship can be established:

[0254]

[0255] wherein, , ;

[0256] The bearing inner ring is simultaneously subjected to the force and moment of the shaft and the ball. The non-linear restoring force calculation formula of a single ball can be obtained by stress analysis:

[0257]

[0258] In the formula, ;

[0259] S124, repeat steps S2-S3 until all balls are calculated;

[0260] S125, sum to obtain the non-linear restoring force.

[0261] Further, as shown in Figure 6 , the motor unbalanced magnetic pull model adopts a motor stator-rotor double-cylinder model. In S100, the motor unbalanced magnetic pull model solution includes the following steps:

[0262] S131, set motor parameters, including initial air gap length δ, stator slot number Z;

[0263] S132, calculate eccentric air gap length :

[0264]

[0265] S133, calculate air gap permeance :

[0266]

[0267] Wherein, as follows:

[0268]

[0269] S134, calculate motor magnetic potential ;

[0270]

[0271] Wherein, is the stator winding fundamental magnetic potential, is the electrical frequency; is the current maximum value; is the number of turns in series per phase winding; is the number of motor pole pairs; is the power factor angle; is the stator winding coefficient, the expression is as follows:

[0272]

[0273] wherein is the number of slots per pole per phase; are the slot pitch angles associated with the number of poles, respectively, expressed as follows:

[0274]

[0275] wherein is the number of stator slots;

[0276]

[0277] wherein is the winding pitch ratio;

[0278] S135, calculating the air gap flux density :

[0279]

[0280] S136, calculating the Maxwell stress σ :

[0281]

[0282] S137, calculating the unbalanced magnetic pull

[0283]

[0284] wherein , , , satisfy the following relationships:

[0285]

[0286] Further, the solving of the main shaft thermal model comprises solving a bearing thermal model, a motor thermal model and a heat transfer model.

[0287] Further, as shown in Figure 7 , the solving of the bearing thermal model comprises the following steps:

[0288] calculating the frictional heat generation of the balls and the inner and outer raceways :

[0289]

[0290] wherein, wherein is the rotational friction coefficient between the balls and the inner raceway, for ball bearings, ; is the spin angular velocity of the balls; is the contact restoring force; ;

[0291] Frictional heating of balls and cage :

[0292]

[0293] where is the ball diameter; is the bearing pitch diameter, which is equal to the average of the inner and outer bearing diameters ; is the initial contact angle; is the cage mass; is the friction factor between the balls and the cage; is the cage angular velocity, ;

[0294] Calculating the total power generated by the bearing :

[0295]

[0296] Further, as shown in Figure 8 , the solution of the motor thermal model comprises the following steps:

[0297] Calculating the asynchronous motor copper losses :

[0298]

[0299] The stator copper losses and the rotor copper losses are calculated by the following equations:

[0300]

[0301] where is the stator winding current; is the winding resistance; is the rotor winding current; is the rotor bar resistance;

[0302] Calculating the asynchronous motor iron losses :

[0303]

[0304] The stator or rotor yoke iron losses are calculated by the following equation:

[0305]

[0306] where is an empirical coefficient, which is 1.3 when the capacity is less than 100 kW, , otherwise 1.3; is the weight of the yoke; is the yoke loss coefficient, expressed as follows:

[0307]

[0308] where, is the loss per unit weight of silicon steel sheet, in W / kg, when is the magnetic induction of the yoke at no load;

[0309] Stator or rotor tooth iron loss is calculated from the following equation:

[0310]

[0311] where, is the weight of the yoke; is an empirical coefficient, taken as 1.8 for an induction motor; tooth loss coefficient;

[0312]

[0313] where, is the loss per unit weight of silicon steel sheet, in W / kg, when is the magnetic induction of the tooth at no load;

[0314] Rotor and air friction loss

[0315]

[0316] where, is the rotational frequency; is the rotor outer diameter; is the rotor axial length; is the air dynamic viscosity; is the average air gap length;

[0317] Asynchronous motor stray loss

[0318]

[0319] where, the copper bar is taken as 0.005, is the output power;

[0320] Total heat generation power

[0321]

[0322] ​​​​​Further, the solution of the heat transfer model includes the following steps:

[0323] The convective heat transfer amount of the stator or rotor core to the air gap

[0324]

[0325] wherein, is the Stefan-Boltzmann constant; is the relative emissivity, usually taken as 0.85; is the heat flow surface temperature; is the heat absorption surface temperature; is the heat sink surface area;

[0326] The surface temperature rise of the stator or rotor core T is calculated by the following formula:

[0327]

[0328] wherein, is the stator or rotor iron loss; is the stator or rotor heat dissipation area, is the heat transfer coefficient per unit surface area to the external air;

[0329] The convective heat transfer amount of the stator or rotor core to the air gap

[0330]

[0331] The heat transfer coefficient of the stator or rotor core to the air gap is calculated by the following formula:

[0332]

[0333] wherein, is the air thermal conductivity; the air thermal conductivity is 0.0267 W / m-K at a normal temperature of 20°C; at 0°C, the air thermal conductivity is 0.0251 W / m-K; at 100°C, the air thermal conductivity is 0.0321 W / m-K;

[0334] Nusselt number Nu is calculated by the following formula:

[0335]

[0336] wherein,

[0337] Taylor number Ta is calculated by the following formula:

[0338] ​​​

[0339] where, is the mass density of the fluid; is the angular velocity of the rotor; is the average radius of the stator and rotor; is the air gap length; is the dynamic viscosity of the fluid;

[0340] Calculating the convective heat transfer from the rotor core to air P 4:

[0341]

[0342] Natural convection heat transfer coefficient is calculated from the equation:

[0343]

[0344] where, is the temperature of the heat leaving surface; is the temperature of the heat absorbing surface; is the outer diameter of the rotor;

[0345] Calculating the convective heat transfer from the stator core to coolant P 5:

[0346]

[0347] Stator core heat transfer coefficient to coolant is calculated from the equation:

[0348]

[0349] where, is the velocity of the coolant; is the length of the motor body;

[0350] Calculating the convective heat transfer from the stator core to air :

[0351]

[0352] Natural convection heat transfer coefficient is calculated from the equation:

[0353]

[0354] where, is the temperature of the radiating surface; is the temperature of the absorbing surface; is the outer diameter of the stator;

[0355] Calculating the amount of heat radiated by the stator core to the air :

[0356]

[0357] where, is the Stefan-Boltzmann constant; is the relative emissivity, usually taken as 0.85; is the radiating surface temperature; is the absorbing surface temperature; is the radiating surface area;

[0358] Calculating the heat dissipated by the bearing :

[0359]

[0360] Calculating the residual heat W :

[0361]

[0362] Calculating the surface temperature of the rotating shaft acted on by the motor :

[0363]

[0364] where, is the surface area of the rotating shaft acted on by the motor heat generation temperature;

[0365] Calculating the surface temperature of the rotating shaft acted on by the bearing :

[0366]

[0367] where, is the surface area of the rotating shaft acted on by the bearing heat generation temperature.

[0368] Example One

[0369] In this example, the rotating shaft material is structural steel, and the rotating speed is taken as 15000 rpm. The simulation obtains the nonlinear restoring force of the front and rear bearings, the unbalanced magnetic pull, the temperature of the front and rear bearings, and the radial component of the rotating shaft, as shown in T bearing , the motor temperature T motor , and the dynamic performance analysis results, as shown in Figures 9 to 15 Figures 16 to 18

[0370] ​​The above merely describes the preferred embodiments of the present application, but the protection scope of the present application is not limited thereto, and any person skilled in the art can easily think of changes or replacements within the technical scope disclosed by the present application, which should be covered in the protection scope of the present application. Therefore, the protection scope of the present application should be subject to the protection scope of the claims.

Claims

1. A multiphysics coupling simulation method for an electric spindle, characterized in that, Includes the following steps: S100. Establish an electric spindle rotor model using finite element software, including: importing the geometric model, setting the material, setting the rotor dynamics module, setting the solid heat transfer module, and mesh generation. The geometric model includes the shaft structure model (300). The rotor dynamics module is configured to apply nonlinear restoring force and unbalanced magnetic pull to the shaft structure model (300); Setting up the solid heat transfer module includes controlling the surface temperature of the motor shaft. T motor The surface temperature of the shaft due to the action of the bearing T bearing The action is applied to the surface of the rotating shaft structure model (300); The nonlinear restoring force is solved using a bearing nonlinear restoring force model, the unbalanced magnetic pull is solved using a motor unbalanced magnetic pull model, and the surface temperature of the motor-acting shaft is solved using a spindle thermal model. T motor The surface temperature of the shaft acting on the bearing T bearing ; S200. Solving using the finite element software includes: coupling the bearing nonlinear restoring force model, the motor unbalanced magnetic pull model, and the spindle thermal model with the electric spindle rotor model using multiphysics fields to obtain the dynamic performance. In S100, solving the bearing nonlinear restoring force model includes the following steps: S121. Set bearing parameters and convergence accuracy, wherein the bearing parameters include the contact angle between the ball and the inner and outer raceways. sphere diameter ; Obtain the relative radial displacement of the inner and outer rings of the bearing along the X-axis. The relative radial displacement of the inner and outer rings of the bearing in the Y-axis direction The relative axial displacement of the inner and outer rings of the bearing The relative angular displacement of the inner and outer rings of the bearing about the X-axis The relative angular displacement of the inner and outer rings of the bearing about the Y-axis direction The X-axis and the Y-axis are perpendicular; S122. Calculate the center of curvature of the outer raceway groove of a single ball. O Final position of the inner raceway groove curvature center axial distance between Calculate the center of curvature of the outer raceway groove of the single ball. O Final position of the inner raceway groove curvature center radial distance between ; in, The distance between the centers of curvature of the inner and outer raceway grooves; , , These are the groove curvature radius coefficients of the inner and outer raceways, respectively; and The following relationship must be satisfied: In the formula, The azimuth angle of the ball; Let be the radius of the trajectory of the inner raceway curvature center, satisfying the following relationship: In the formula, The diameter of the pitch circle; constraint , Find the initial value and substitute it into the following formula: In the formula, , These are the axial and radial projections of the distance between the final position of the ball center and the center of curvature of the outer raceway, respectively. According to Hertz contact theory, the contact deformation between the ball and the inner raceway... and the amount of contact deformation between the ball and the outer raceway It can be calculated using the following formula: In the formula, The normal force representing contact deformation. , , , , Calculated using the following formulas respectively: In the formula, , The elastic modulus of the two contacting bodies; , The ratio of the two contacting bodies is Poisson's ratio. and Let be the radii of curvature of the two contacting bodies on the two principal planes, respectively; Ball in contact with inner raceway: Ball in contact with outer raceway: In the formula, ; Since rolling bearings are lubricated, the effect of lubricating oil film thickness should also be considered: In the formula, ， These are the distances from the center of curvature of the inner and outer raceway grooves to the center of the ball, respectively. , The thicknesses of the center oil film between the ball and the inner and outer raceways are respectively, and can be calculated using the following formula: in, ; e Here, is a natural constant; other unknowns are determined by the following formula: Based on the positional relationship between the curvature centers of the inner and outer raceways and the center of the sphere, the following geometric relationship also exists: In the formula, , These are the contact angles between the ball and the inner and outer raceways, respectively. Number the balls, satisfying , The total number of balls; S123, Establishing a single sphere subjected to centrifugal force With gyro torque The force relationship between the inner and outer raceways under the action of the action is determined and solved: in, , These are the coefficients of friction between the ball and the inner and outer raceways, respectively. , These are the contact angles between the ball and the inner and outer raceways, respectively. The diameter of the sphere, , The normal force representing the deformation of the ball upon contact with the inner and outer raceways; centrifugal force. With gyro torque Determined by the following formula: in, Let be the angular velocity of the ball's rotation. Let be the angular velocity of the ball's revolution. Let be the moment of inertia of the ball. This refers to the rotational speed of the bearing inner ring. The angle between the axis of rotation and the axis of revolution. The mass of the ball; Based on Jones's raceway control assumption, namely that bearings are more closely controlled by the outer raceway at high speeds, the following relationship can be established: in, , ; The inner ring of the bearing is simultaneously subjected to the forces and torques of the shaft and the balls. By performing a force analysis, the nonlinear restoring force of a single ball can be calculated: In the formula, ; S124. Repeat steps S2 to S3 until all balls have been calculated. S125. Summation yields the nonlinear restoring force.

2. The multiphysics coupling simulation method for electric spindles according to claim 1, characterized in that, The geometric model also includes a tool head structure model (100), a tool holder structure model (200), and a motor rotor structure model (400). The tool head structure model (100) and the shaft structure model (300) are respectively connected to the two ends of the tool holder structure model (200), and the motor rotor structure model (400) is fitted around the outer periphery of the shaft structure model (300).

3. The multiphysics coupling simulation method for electric spindles according to claim 1, characterized in that, The rotating shaft structure model (300) has a cavity (301) that runs through the axial direction.

4. The multiphysics coupling simulation method for electric spindles according to claim 1, characterized in that, In S100, applying the unbalanced magnetic pull force to the rotating shaft structure model (300) includes the following steps: S111. Divide the motor rotor structure model (400) into multiple shaft segments along its axial direction; S112. Establish a coordinate system with the center point of the end face away from the motor on each shaft segment as the origin, and obtain the position information of the selected node in real time through the probe built into the finite element software. S113. Use the position information as the input variable for the motor unbalanced magnetic pull model; S114. Apply each unbalanced magnetic force to the corresponding shaft segment.

5. The multiphysics coupling simulation method for electric spindles according to claim 1, characterized in that, In S100, the solution to the motor unbalanced magnetic pull model includes the following steps: S131. Set motor parameters, including initial air gap length δ and stator slot number Z; S132. Calculate the eccentric air gap length. : S133, Calculate the air gap permeability : in, as follows: S134, Calculate the motor magnetomotive force ; in, For the fundamental magnetomotive force of the stator winding, It is the electrical frequency; This represents the maximum current value. This refers to the number of turns connected in series in each phase winding; This represents the number of pole pairs of the motor. Power factor angle; The stator winding coefficient is expressed as follows: In the formula, The number of slots per pole per phase; These are the slot pitch angles related to the number of poles, expressed as follows: in, This refers to the number of stator slots; In the formula, This refers to the winding pitch ratio; S135, Calculate the air gap magnetic flux density : S136. Calculate Maxwell stress σ : S137. Calculate the unbalanced magnetic pull: ； ； in, , , , The following relationship must be satisfied: 。 6. The multiphysics coupling simulation method for electric spindles according to claim 1, characterized in that, Solving the spindle thermal model includes solving the bearing thermal model, the motor thermal model, and the heat transfer model.

7. The multiphysics coupling simulation method for electric spindles according to claim 6, characterized in that, Solving the bearing thermal model includes the following steps: Calculate the frictional heat generated by the ball and the inner and outer raceways. : In the formula, where It is the coefficient of rotational friction between the ball and the inner raceway. For ball bearings, ; It is the angular velocity of the ball's spin; It is contact resilience; ; Calculate the heat generated by friction between the ball and the cage. : In the formula, Sphere diameter; The bearing pitch circle is equal to the average of the bearing's inner and outer diameters. ; Initial junction angle; To maintain rack quality; The coefficient of friction between the ball and the cage; To maintain the frame angular velocity, ; Calculate the total heat generated by the bearing : 。 8. The multiphysics coupling simulation method for electric spindles according to claim 6, characterized in that, Solving the thermal model of the motor includes the following steps: Calculate the copper loss of an asynchronous motor : Stator copper loss and rotor copper loss Calculated by the following formula: In the formula, This refers to the stator winding current. For winding resistance; This refers to the rotor winding current; For rotor bar resistance; Calculate the iron loss of an asynchronous motor : Stator or rotor yoke losses Calculated by the following formula: In the formula This is an empirical coefficient, when the capacity hour, Otherwise, it is 1.3; The weight of the yoke; The yoke loss coefficient is expressed as follows: In the formula For when The loss per unit weight of silicon steel sheet, expressed in W / kg; The magnetic flux density of the unloaded yoke; Stator or rotor tooth iron loss Calculated by the following formula: In the formula, The weight of the yoke; This is an empirical coefficient, taken as 1.8 for induction motors; This is the tooth loss coefficient; In the formula, For when The loss per unit weight of silicon steel sheet, expressed in W / kg; The magnetic flux density of the unloaded tooth; Calculate the friction loss between the rotor and the air. : In the formula, For frequency conversion; The outer diameter of the rotor; This refers to the axial length of the rotor. Aerodynamic viscosity; This represents the average air gap length. Calculate stray losses of asynchronous motors : In the formula, copper conductor bar Take 0.005, This refers to the output power. Calculate the total heat generation power : 。 9. The multiphysics coupling simulation method for electric spindles according to claim 6, characterized in that, Solving the heat transfer model includes the following steps: Calculate the convective heat transfer radiated from the stator or rotor core into the air gap. : In the formula, It is the Stefan-Boltzmann constant; The relative emissivity is typically taken as 0.85; Temperature of the surface from which heat flows out; Temperature of the heat absorption surface; This refers to the surface area of ​​the heat sink. Temperature rise on the surface of stator or rotor core T Calculated by the following formula: In the formula, For stator or rotor iron loss; For stator or rotor heat dissipation area, The heat transfer coefficient per unit surface area to the outside air; Calculate the convective heat transfer from the stator or rotor core to the air gap. : Heat transfer coefficient from stator or rotor core to air gap Calculated by the following formula: In the formula, The thermal conductivity of air is 0.0267 W / mK at room temperature (20℃); 0.0251 W / mK at 0℃; and 0.0321 W / mK at 100℃. Nusel number Nu Calculated by the following formula: In the formula, ; Taylor number Ta Calculated by the following formula: In the formula, The mass density of the fluid; This refers to the rotor angular velocity; The average radius of the stator and rotor; This is the air gap length; For fluid dynamic viscosity; Calculate the convective heat transfer from the rotor core ends to the air. P 4: Natural convection heat transfer coefficient Calculated by the following formula: In the formula, Temperature of the surface from which heat flows out; Temperature of the heat absorption surface; The outer diameter of the rotor; Calculate the convective heat transfer from the stator core to the coolant. P 5: Stator core heat transfer coefficient to coolant Calculated by the following formula: In the formula, For coolant velocity; This refers to the length of the motor body; Calculate the convective heat transfer from the stator core to the air. : Natural convection heat transfer coefficient Calculated by the following formula: In the formula, The temperature of the radiant surface; Temperature of the absorption surface; Stator outer diameter; Calculate the amount of heat radiation from the stator core to the air. : In the formula, It is the Stefan-Boltzmann constant; The relative emissivity is typically taken as 0.85; The temperature of the radiant surface; Temperature of the absorption surface; The surface area of ​​the radiating body; Calculate bearing heat dissipation : Calculate the remaining heat W : Calculate the surface temperature of the motor shaft : In the formula, The surface area of ​​the shaft is affected by the heat generated by the motor. Calculate the surface temperature of the bearing shaft : In the formula, The surface area of ​​the shaft affected by the bearing's heat generation temperature.

Images