Nonlinear vibration coupling analysis method for aerostatic spindle under sliding mode control

Abstract

The invention provides a nonlinear vibration coupling analysis method for an aerostatic spindle under sliding mode control, and relates to the technical field of machine tool precision machining. A multi-degree-of-freedom system dynamics model is constructed by performing nonlinear coupling calculation on gas film force, unbalanced magnetic pulling force and mass eccentric excitation. Time-varying nonlinearity of bearing rigidity is accurately described through calculation of dynamic air film force, the two-way coupling effect between an electromagnetic field and mechanical motion is revealed through introduction of unbalanced magnetic pulling force, and complex combined frequency vibration can be generated through interaction of mass eccentricity serving as classical external synchronous excitation and the two kinds of internal excitation. According to the method, the controller design can consider and compensate the nonlinear time variation of the gas film force, the state dependence disturbance of the magnetic pulling force and the eccentric force excitation in advance, the reverse influence of the control action on the gas film stability and the electromagnetic field can be evaluated, and the real mechanical-electric-gas-control collaborative optimization is realized.

Description

Technical Field

[0001] This invention relates to the field of precision machining technology for machine tools, and in particular to a nonlinear vibration coupling analysis method for an air hydrostatic spindle under sliding diaphragm control. Background Technology

[0002] Air hydrostatic spindles are widely used in optical material processing due to their excellent stability. During operation, the spindle bearing is subjected to load excitations from multiple sources, resulting in complex and diverse vibration behaviors. Various physical fields are coupled with each other through displacement amplitude responses, exerting a combined effect. Nonlinear vibrations cause instantaneous changes in different load characteristic parameters, thus affecting motion stability. Under extreme conditions, excessive vibration may damage the bearing structure and reduce machining accuracy. Therefore, implementing active real-time control can not only prevent spindle wear but also improve machining accuracy and efficiency. Establishing a coupled model integrating active control to study the interaction between various load excitations can enable vibration behavior prediction. Simulation analysis based on this model can compare the vibration suppression effect under different operating conditions, providing a theoretical basis for the selection of control parameters.

[0003] Currently, most studies simplify the air film force into a linear spring-damped model for analysis purposes. However, the nature of air hydrostatic bearings is highly nonlinear. A strong nonlinear relationship exists between the air film pressure distribution and rotor displacement, especially under large disturbances or high-speed conditions, where linear models cannot accurately describe the stiffness and damping characteristics of the air film. The tiny air gap between the spindle rotor and stator creates a complex fluid-structure interaction field. Rotor vibration disturbs the air film pressure distribution, and the changing pressure distribution reacts back to the rotor, forming a tightly coupled nonlinear system. Existing models often ignore or simplify this strong coupling mechanism, leading to a significant decrease in model prediction accuracy at critical speeds or under external disturbances. To suppress spindle vibration, traditional linear control strategies such as PID controllers are commonly used. Although PID controllers are simple in structure and easy to implement, they exhibit significant limitations when dealing with a controlled object like an air hydrostatic spindle, which possesses nonlinearity, time-varying parameters, and uncertainties related to external disturbances. Sliding mode control (SMC), as a robust nonlinear control method, is theoretically well-suited to address the uncertainties and nonlinearities of air hydrostatic spindles.

[0004] Current research lacks a unified modeling method that can accurately describe the coupling relationship between sliding mode active control and the nonlinear vibration of the air hydrostatic spindle. This deficiency leads to blind spots in the design of existing control systems, failing to fundamentally understand and solve the complex vibration behavior of the spindle under strong robust control, thus limiting the potential of the air hydrostatic spindle to achieve ultimate accuracy under extreme conditions. Summary of the Invention

[0005] To address the shortcomings of existing technologies, a coupled analysis method for nonlinear vibration of an air hydrostatic master shaft under sliding diaphragm control is proposed. This method not only improves computational efficiency and reduces costs, but also provides a theoretical basis for predicting nonlinear vibration phenomena under active control.

[0006] On one hand, the present invention provides a nonlinear vibration coupling analysis method for an air hydrostatic master shaft under sliding diaphragm control, comprising the following steps:

[0007] Step 1: Construct an air static pressure spindle system and obtain real-time active control excitation under sliding diaphragm control based on the position error signal and the closed-loop feedback principle;

[0008] Based on the closed-loop control of the air static pressure spindle, a special tool post fixed to the tool holder was customized. The special tool post integrates a pair of deep groove ball bearings, four piezoelectric actuators, and a tool post housing fixed with bolts. The inner ring of the deep groove ball bearing rotates at high speed with the tool, while the outer ring contacts the piezoelectric actuator to transmit real-time control force. The dynamic amplitude signal of the spindle is detected by two sets of radially arranged displacement sensors. The acquisition card converts the analog signal into digital data and transmits it to the computer. The computer calculates the compensation voltage in real time through the active sliding diaphragm control program. The calculated compensation voltage digital signal is converted into an analog signal and sent to the piezoelectric controller. After signal amplification, it is applied to the four piezoelectric actuators to realize closed-loop instantaneous active compensation control.

[0009] Specifically: Based on the core idea of ​​sliding membrane control, the sliding membrane surface s is calculated using the following formula: ;

[0010] In the formula, λ is the speed error, λ is the slope of the sliding surface, and e is the position error;

[0011] Calculate the control function u of the active slicker control program: ;

[0012] In the formula, b is the control gain, and f is the nonlinear function of the air static pressure spindle system. It's a speed error. λ is the second derivative of the desired position, K is the set gain, s is the sliding surface, λ is the sliding surface parameter, u is the control gain, and sgn is the switching control law.

[0013] Considering the four mutually perpendicular radial directions, calculate the control gain u experienced by the piezoelectric actuator in each of the four directions. x1 u x2 u y1 and u y2 :

[0014] ;

[0015] In the formula, u x1It is the active control gain in the positive x-axis direction; u x2 It is the active control gain in the negative x-axis direction; u y1 It is the active control gain in the positive y-axis direction; u y2 It is the active control gain in the negative y-axis direction, k x1 It is the stiffness coefficient in the positive x-axis direction, k x2 It is the stiffness coefficient in the negative x-axis direction, k y1 It is the stiffness coefficient in the positive y-axis direction, k y2 It is the stiffness coefficient in the negative y-axis direction, c x1 It is the damping coefficient in the positive x-axis direction, c x2 It is the damping coefficient in the negative x-axis direction, c y1 It is the damping coefficient in the positive y-axis direction, c y2 It is the damping coefficient in the negative y-axis direction. It is the second derivative of the desired position in the positive x-axis direction. It is the second derivative of the desired position in the negative x-axis direction. It is the second derivative of the desired position in the positive y-axis direction. It is the second derivative of the desired position in the negative y-axis direction. It is the positive velocity error along the x-axis. It is the negative velocity error along the x-axis. It is the positive y-axis velocity error. It is the negative velocity error along the y-axis. It is the positive x-axis sliding surface. It is the negative x-axis sliding surface. It is the positive y-axis sliding surface. It is the negative y-axis sliding surface.

[0016] Based on the determination of the control function, the active control excitation for real-time compensation of the sliding membrane control is calculated:

[0017] ;

[0018] In the formula, F controlx1 It is the active control force in the positive x-axis direction; F controlx2 It is the active control force in the negative x-axis direction; F controly1 It is the active control force in the positive y-axis direction; F controly2 It is the active control force in the negative y-axis direction; M controlx It is the active sliding membrane control bending moment in the x-direction, M controly It is an active sliding membrane control of bending moment in the y-direction; L T It is the axial length of the piezoelectric actuator from the center of the spindle system; d 33 It is the piezoelectric constant; K s It refers to the stiffness of piezoelectric ceramics.

[0019] Step 2: Based on the initial air film conditions, considering the changes in air film under the influence of vibration, derive the dynamic air film force excitation of the bearing;

[0020] Step 2.1: Calculate the instantaneous radial air film force excitation experienced by the spindle during operation;

[0021] Calculate the radial instantaneous air film thickness h J :

[0022] ;

[0023] In the formula, Δx is the vibration displacement in the radial x-direction; Δy is the vibration displacement in the radial y-direction; Δφ x0 It is the deflection angle in the radial x-direction; Δφ y0 It is the deflection angle in the radial y-direction; It is the initial eccentricity deflection angle; c0 is the film gap of the radial bearing; e0 is the initial eccentricity of the radial bearing; φ x0 It is the initial deflection state in the radial x-direction; φ y0 It represents the initial deflection state in the radial y-direction; θ is the angular position coordinate; z is the axial position coordinate; z m is the axial position coordinate of the midpoint of the main shaft; h0 is the air film thickness considering vibration at the initial position;

[0024] Substituting the obtained gas film thickness into the two-dimensional finite difference equation, the gas film pressure at each point is calculated:

[0025] ;

[0026] In the formula, the dimensionless parameters are defined as follows:

[0027] ;

[0028] In the formula, p a It is atmospheric pressure; c a It is the radial initial film gap; L is the bearing length; p s ρ is the gas supply pressure; μ is the gas viscosity; u1 is the operating speed; p is the actual pressure at each grid point; z is the axial position of each grid point; h is the actual gas film thickness at each grid point.

[0029] Based on the principle of flow conservation, an equilibrium region is defined at the orifice location, and Newton's iterative method is applied within this region to solve for the air pressure at each point. The termination condition for the calculation is set as follows:

[0030] ;

[0031] In the formula, The pressure at node (i, j) after N iterations represents the air film pressure; i and j represent the node coordinates; N is the number of iterations.

[0032] Based on the determined pressure conditions, calculate the instantaneous radial film force excitation experienced by the spindle during operation:

[0033] ;

[0034] In the formula, It is the film force in the radial x-direction; It is the film force in the radial y-direction; It is the bending moment of the air film force in the radial x direction; It is the bending moment of the air film force in the radial y direction; It is the dimensionless air film pressure at each node;

[0035] Step 2.2: Calculate the instantaneous axial film excitation force experienced by the spindle during operation;

[0036] Calculate the thickness h of the air film on both sides of the axis. T1 and h T2 :

[0037] ;

[0038] In the formula, Δz is the vibration displacement in the axial z-direction; Δφ x It is the deflection angle of the thrust bearing in the x-direction; Δφ y It is the deflection angle of the thrust bearing in the y-direction; c T1 and c T2 It is a bidirectional closed thrust film gap, φ zx0 and φ zy0 It represents the initial deflection state of the thrust bearing, where r is the radius of the circle containing the differential grid point;

[0039] Substituting the obtained film thickness into the two-dimensional finite difference equation expressed in polar coordinates, the film pressure at each point is calculated:

[0040] ;

[0041] In the formula, the dimensionless parameter is defined as:

[0042] ;

[0043] In the formula, p a It is atmospheric pressure; c T R2 is the initial axial film gap; R2 is the thrust disc radius; p s It is the supply pressure; μ is the gas viscosity; w is the operating angular velocity; w s is the rated angular velocity; p is the actual pressure at each grid point; r is the radius of the circle containing the differential grid point; h is the actual air film thickness at each grid point.

[0044] Considering the boundary conditions of the inner and outer rings, the air film within the bidirectional closed thrust bearing is divided into finite difference grid points. Flow balance elements are set based on each throttling orifice, and the pressure distribution at each point is finally determined using Newton's iteration method, with the termination condition being:

[0045] ;

[0046] In the formula, The gas film pressure at node (r, θ) after N iterations represents the gas film pressure; r and θ represent the polar coordinates of the node; N is the number of iterations.

[0047] Based on the determined pressure conditions, calculate the instantaneous axial film force excitation experienced by the spindle during operation:

[0048] ;

[0049] In the formula, It is the dimensionless gas film pressure at each node. Dimensionless polar radius.

[0050] Step 3: Considering the relative positional relationship between the stator and rotor and the change in the center of mass deflection angle, derive the unbalanced magnetic pull excitation and mass eccentricity excitation under the influence of nonlinear vibration;

[0051] Step 3.1: Calculate the unbalanced magnetic pull excitation experienced by the spindle during operation;

[0052] Based on Kirchhoff's voltage and current laws, the stator and rotor voltages, electromotive force, and core winding excitation are given by the following formula:

[0053] ;

[0054] In the formula, , ,and These are the excitation current, stator current, and rotor phase current, respectively, α Fe It is the stator iron loss angle, φ1 and φ2 are the power factor angles of the stator and rotor, respectively, I em I s and The root mean square values ​​of the excitation current, stator current, and rotor current, respectively. It is the power supply voltage applied to each phase of the stator. It is a single-phase induced electromotive force in the stator. It is the equivalent rotor induced electromotive force, Z s It is the stator leakage reactance, Z em It is the magnetizing impedance. The equivalent impedance of the rotor is expressed as:

[0055] ;

[0056] In the formula, R s and X s These are the stator single-phase resistance and reactance, and These are the rotor equivalent resistance and reactance, R. em and X em These are the motor's excitation resistance and excitation reactance, respectively. R le is the rotor equivalent additional resistance, j is the phase coordinate, and s is the slip of the induction motor.

[0057] Calculate the slip rate of the induction motor through speed analysis: ;

[0058] In the formula, n is the spindle speed, n s =60f cur / p is the synchronous rotational speed of the rotating magnetic field, f cur p and p are the stator current frequency and the number of pole pairs, respectively, and the additional resistance R is... le The power dissipated above represents the mechanical power P generated by the electric spindle. e .

[0059] Considering the actual mechanical output power in the electric spindle, calculate the additional resistance R in the equivalent circuit. le : ;

[0060] An equivalent circuit for the induction motor was established based on the root mean square values ​​of voltage and current. Given the motor parameters and power supply voltage, the excitation current, stator current, and rotor phase current were calculated.

[0061] ;

[0062] Furthermore, the rotor equivalent current is obtained. With the output mechanical power P of the electric spindle e The relationship between them:

[0063] ;

[0064] Solving the above equation yields the current. Then the stator current is obtained. Excitation current and motor slip s;

[0065] Calculate the three-phase AC power required for the built-in motor:

[0066] ;

[0067] In the formula, ω cur =2πf cur It is the angular velocity of the rotating magnetic field, θ. cur=pα is the electrical angle, p is the pole pair number, and α is the spatial angle. It is the equivalent phase current, i s It is the stator phase current, i em It is the excitation phase current, t is the operating time, and α is the excitation phase current. Fe It is the stator iron damage angle. It is the stator power factor angle. It is the rotor power factor angle. It is the stator space phase angle. This is the stator spatial phase angle. The fundamental magnetomotive force F of the rotor under symmetrical load is calculated using the following formula. rf The fundamental magnetomotive force F of the stator sf The fundamental magnetomotive force F of the air gap m :

[0068] ;

[0069] In the formula, m1=3 is the number of stator phases, N1 is the number of turns of the series coil per pole per phase of the stator, and k w1 It is the stator winding distribution factor. It is the equivalent phase current.

[0070] In the air gap magnetic field between the rotor and stator, the expression for the normal magnetic flux density is:

[0071] ;

[0072] In the formula, δ n It is the normal magnetic gap length on the rotor surface; It is the normal magnetic flux density on the outer surface of the rotor; It is the air gap permeability per unit area of ​​the rotor's outer surface; μ0 is the vacuum permeability; φ3 is the spatial azimuth angle; f m It is magnetic potential; ω cur It is the angular velocity of the rotating magnetic field; θ cur It is the electric field angle; α is the rotor power factor angle; k is the air gap azimuth angle; c It is the stator winding distribution coefficient.

[0073] When the magnetic material on the rotor comes into contact with the magnetic field in the air gap, a magnetic force is generated. Calculation of conventional magnetic force... :

[0074] ;

[0075] In the formula, It is the magnetic flux density perpendicular to the rotor surface; It is the magnetic pull perpendicular to the rotor surface; α Fe It is the stator iron loss angle; It is the spatial phase of σ.

[0076] 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:

[0077] ;

[0078] In the formula, It is the average 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... si It is the rotor outer diameter, R ro θ is the stator inner diameter, r is the instantaneous linear distance between the rotor cross-section considering vibration and the rotor cross-section not considering vibration, Δx is the vibration displacement in the x-direction, Δy is the vibration displacement in the y-direction, and θ is the vibration displacement in the y-direction. x It is the vibration deflection angle in the x-direction, θ y θ is the vibration deflection angle in the y direction, θ is the dynamic eccentricity azimuth angle, and l is the effective axial length of the rotor.

[0079] Calculate the instantaneous magnetic circuit length considering the vibration response effects in five degrees of freedom. :

[0080] ;

[0081] In the formula, D si =2R si and D ro =2R ro These are the inner diameter of the stator and the outer diameter of the rotor of the built-in motor; (R−R) ro This indicates the change in air gap length caused by the change in the cross-sectional shape of the built-in motor rotor; It is the rotor radial deflection angle. It refers to the axial relative position.

[0082] Based on geometric coordinate relationships, calculate the horizontal, vertical, and axial projected components of the magnetic pull force acting at any point on the rotor surface, as well as the transverse magnetic moment:

[0083] ;

[0084] In the formula, It is the horizontal projection component of the unbalanced magnetic pull. It is the vertical projection component of the unbalanced magnetic pull. It is the axial projection component of the unbalanced magnetic pull. It is the transverse magnetic moment component of the unbalanced magnetic pull.

[0085] Furthermore, by integrating the horizontal, vertical, and axial projection components of the magnetic pull force in the rotor surface normal direction, the projections of the unbalanced magnetic pull force excitation in the X, Y, and Z axes and the components of its additional transverse electromagnetic torque in the X and Y axes are calculated:

[0086] ;

[0087] In the formula, It is an unbalanced magnetic pull effect horizontal excitation. It is an unbalanced magnetic pull effect with vertical excitation. It is an axial excitation caused by unbalanced magnetic pull effect. It is the transverse magnetic moment in the x-direction due to the unbalanced magnetic pull effect. It is the transverse magnetic moment in the y-direction due to the unbalanced magnetic pull effect.

[0088] Step 3.2: Calculate the unbalanced mass eccentric excitation experienced by the spindle during operation;

[0089] The centrifugal force caused by the eccentric mass is projected along the radial and axial directions of the stator, resulting in unbalanced mass excitation. Specifically, the projections of the unbalanced mass excitation onto the X, Y, and Z axes, and the components of the additional lateral moment in the X and Y axes, are calculated:

[0090] ;

[0091] In the formula, m is the mass of the spindle system, w is the angular velocity of the spindle, v is the spatial tilt of the centrifugal force relative to the stator radial direction, and e is the spindle eccentricity. It is the initial eccentricity angle in the x-direction. It is the initial eccentricity angle in the y-direction. It is a quality eccentricity level excitation, It is a mass eccentric vertical excitation. It is axial excitation due to mass eccentricity. It is the eccentric bending moment of mass in the x-direction. It is the eccentric bending moment of mass in the y-direction.

[0092] Step 4: Due to the coupling relationship between multiple physics fields, construct a multi-degree-of-freedom dynamic model under complex excitation influence;

[0093] Considering the coupling effect of multiple excitation loads, the specific dynamic equations of the air hydrostatic principal shaft under sliding film control are as follows:

[0094] ;

[0095] In the formula M s G s and C s These are the mass matrix, gyroscope matrix, and damping matrix of the unit node, respectively, Fair It is based on the radial and axial air film force excitation calculated by dynamic air film thickness, F ump It is the unbalanced magnetic excitation generated by the built-in motor under the condition of air gap change, F e It is an unbalanced mass eccentric excitation under the action of instantaneous deflection angle, F control It is an active sliding membrane control excitation determined based on vibration signals, where q is the vibration displacement of each element node. It is the first derivative of the vibration displacement of each element node. It is the second derivative of the vibration displacement of each unit node.

[0096] Based on the nonlinear finite element dynamics model, the displacement and excitation distribution of each node are shown in the following matrix:

[0097] ;

[0098] ;

[0099] ;

[0100] ;

[0101] ;

[0102] In the formula, x1 is the x-direction displacement of node 1 of the mass element in the dynamic model, y1 is the y-direction displacement of node 1 of the mass element in the dynamic model, z1 is the z-direction displacement of node 1 of the mass element in the dynamic model, and θ x1 θ is the x-direction deflection angle of node 1 of the mass element in the dynamic model. y1 It is the y-direction deflection angle of node 1 of the mass element in the dynamic model, x 19 It is the x-direction displacement of node 19 of the mass element in the dynamic model, y 19 It is the y-direction displacement of node 19 of the mass element in the dynamic model, z 19 It is the z-direction displacement of node 19 of the mass element in the dynamic model, θ x19 θ is the x-direction deflection angle of node 19 of the mass element in the dynamic model. y19 It is the y-direction deflection angle of node 19 of the mass element in the dynamic model, F. xL It is the excitation of the gas film force of the left radial gas bearing in the x direction, F yL It is the left radial gas bearing film force excitation in the y direction, M xL M is the left radial gas bearing film force bending moment in the x-direction. yL F is the left radial gas bearing film force bending moment in the y-direction. xR It is the excitation of the film force of the right radial gas bearing in the x-direction, F yRIt is the right radial gas bearing film force excitation in the y direction, M xR M is the right radial gas bearing film force bending moment in the x-direction, M. yR F is the right radial gas bearing film force bending moment in the y-direction. zT It is the film force excitation of the thrust gas bearing in the axial direction, M zx M is the bending moment of the thrust gas bearing film force in the x-direction. zy It is the bending moment of the thrust gas bearing film force in the y direction, F controlx It is an active sliding membrane control excitation in the x-direction, F controly It is an active sliding membrane control excitation in the y-direction, M controlx It is the active sliding membrane control bending moment in the x-direction, M controly It is an active sliding membrane control of bending moment in the y direction;

[0103] Step 5: Solve the dynamic equations based on the Newmark-β method to calculate the dynamic response of the step size node, and calculate the active sliding film control excitation, dynamic air film force excitation, unbalanced magnetic pull force excitation and mass eccentricity excitation by the changes caused by the displacement location;

[0104] Step 6: Based on the active sliding membrane control excitation, dynamic air film force excitation, unbalanced magnetic pull force excitation and mass eccentricity excitation, determine the instantaneous load, and couple it with the dynamic equation to solve the real-time state of the displacement field. If the displacement trajectory does not show a convergence trend, repeat steps 1 to 5 until the node changes are stable and converged, and obtain the main shaft displacement field response curve under sliding membrane control for dynamic behavior characteristic analysis.

[0105] On the other hand, this application proposes an electronic device, including: one or more processors, and a memory for storing instructions, which, when executed by the one or more processors, cause the one or more processors to perform the nonlinear vibration coupling analysis method for an air hydrostatic spindle under sliding diaphragm control.

[0106] Thirdly, this application proposes a computer-readable storage medium storing executable instructions that, when executed, cause a processor to perform the aforementioned nonlinear vibration coupling analysis method for an air hydrostatic spindle under sliding diaphragm control.

[0107] Fourthly, this application proposes a computer program product, including a computer program or instructions, which, when executed by a processor, implements the aforementioned nonlinear vibration coupling analysis method for an air hydrostatic spindle under sliding diaphragm control.

[0108] The beneficial effects of adopting the above technical solution are as follows:

[0109] This invention provides a method for nonlinear vibration coupling analysis of an air hydrostatic master shaft under sliding diaphragm control, which has the following advantages:

[0110] (1) This invention is based on the active sliding film control law. The sliding film control is designed with a specific sliding surface so that once the system state is reached, it is insensitive to external disturbances and changes in internal parameters. This can effectively overcome the vibration problem caused by the inherent uncertainty of air static pressure spindle due to nonlinearity of air film force and fluctuation of air supply pressure, and significantly improve the operating stability and accuracy of spindle under complex working conditions.

[0111] (2) The coupling modeling path proposed in this invention achieves accurate description and coordinated suppression of complex vibration phenomena. By establishing a nonlinear vibration coupling model, the energy transfer and modal interaction between multiple degrees of freedom such as radial and axial directions within the spindle system can be deeply revealed. The sliding membrane controller designed based on this model can perform multi-input multi-output coordinated control, suppressing coupled vibrations from a holistic rather than local perspective, thereby comprehensively improving the dynamic performance of the system.

[0112] (3) This invention establishes a multi-excitation vibration model including active control for an air static pressure spindle. By using multi-directional vibration displacement as correlated variables, the interaction relationship between various loads is derived to study transient dynamic characteristics. An active control closed-loop model is constructed, and real-time vibration suppression is achieved through error compensation. By comparing the nonlinear behavior of multiple parameters under active and non-active control conditions, a reference basis is provided for the selection and adjustment of operation control parameters. Attached Figure Description

[0113] Figure 1 This is a schematic diagram of active sliding film control for an air static pressure spindle provided in an embodiment of the present invention;

[0114] Figure 2 This is a schematic diagram of the dynamic change of the radial air film on the air static pressure spindle provided in an embodiment of the present invention;

[0115] Where (a) is the front view, (b) is the xoy section after translation at AA, and (c) is the xoy section after flipping at AA;

[0116] Figure 3 This is a schematic diagram of the dynamic change of the axial air film on the air static pressure spindle provided in an embodiment of the present invention;

[0117] Where (a) - front view, (b) - xoy section, (c) - xoz section;

[0118] Figure 4 This is a schematic diagram illustrating the analysis of the unbalanced magnetic pull of the built-in induction motor provided in an embodiment of the present invention;

[0119] (a) - kinematic analysis diagram, (b) - instantaneous air gap distribution diagram between stator and rotor, (c) - magnetic stress on rotor surface and its projection diagram;

[0120] Figure 5 A schematic diagram of unbalanced mass eccentric excitation analysis considering spindle vibration displacement provided in an embodiment of the present invention;

[0121] Figure 6 A schematic diagram of the dynamic model of the air static pressure spindle system under active sliding film control provided in an embodiment of the present invention. Detailed Implementation

[0122] The specific implementation methods of this application will be further described in detail below with reference to the accompanying drawings and embodiments.

[0123] Example 1:

[0124] On the one hand, this invention provides a nonlinear vibration coupling analysis method for an air-static spindle under sliding diaphragm control. Based on the analysis of these multi-load excitations, a multi-degree-of-freedom nonlinear dynamic model of the spindle system is established using the finite element principle. A closed-loop iterative algorithm is used to achieve coupled solution; the displacement response of each node of the air-static spindle under active control is calculated; including the following steps:

[0125] Step 1: Construct an air static pressure spindle system and obtain real-time active control excitation under sliding diaphragm control based on the position error signal and the closed-loop feedback principle;

[0126] Based on air static pressure spindle closed-loop control Figure 1 As shown, a custom-designed tool holder is fixed to the tool holder. This tool holder integrates a pair of deep groove ball bearings, four piezoelectric actuators, and a bolt-fixed tool holder housing. The inner ring of the deep groove ball bearing rotates at high speed with the tool, while the outer ring contacts the piezoelectric actuators to transmit real-time control force, effectively suppressing spindle vibration displacement. Two sets of radially arranged displacement sensors detect the dynamic amplitude signal of the spindle. The acquisition card converts the analog signal into digital data and transmits it to the computer. The computer calculates the compensation voltage in real time through an active sliding diaphragm control program. The calculated compensation voltage digital signal is converted back into an analog signal and sent to the piezoelectric controller. After signal amplification, it is applied to the four piezoelectric actuators, achieving closed-loop instantaneous active compensation control.

[0127] Specifically: Based on the core idea of ​​sliding membrane control, the sliding membrane surface s is calculated using the following formula: ;

[0128] In the formula, λ is the speed error, λ is the slope of the sliding surface, and e is the position error;

[0129] Calculate the control function u of the active slicker control program: ;

[0130] In the formula, b is the control gain, and f is the nonlinear function of the air static pressure spindle system. It's a speed error. λ is the second derivative of the desired position, K is the set gain, s is the sliding surface, λ is the sliding surface parameter, u is the control gain, and sgn is the switching control law.

[0131] Considering the four mutually perpendicular radial directions, calculate the control gain u experienced by the piezoelectric actuator in each of the four directions. x1 u x2 u y1 and u y2 :

[0132] ;

[0133] In the formula, u x1 It is the active control gain in the positive x-axis direction; u x2 It is the active control gain in the negative x-axis direction; u y1 It is the active control gain in the positive y-axis direction; u y2 It is the active control gain in the negative y-axis direction, k x1 It is the stiffness coefficient in the positive x-axis direction, k x2 It is the stiffness coefficient in the negative x-axis direction, k y1 It is the stiffness coefficient in the positive y-axis direction, k y2 It is the stiffness coefficient in the negative y-axis direction, c x1 It is the damping coefficient in the positive x-axis direction, c x2 It is the damping coefficient in the negative x-axis direction, c y1 It is the damping coefficient in the positive y-axis direction, c y2 It is the damping coefficient in the negative y-axis direction. It is the second derivative of the desired position in the positive x-axis direction. It is the second derivative of the desired position in the negative x-axis direction. It is the second derivative of the desired position in the positive y-axis direction. It is the second derivative of the desired position in the negative y-axis direction. It is the positive velocity error along the x-axis. It is the negative velocity error along the x-axis. It is the positive y-axis velocity error. It is the negative velocity error along the y-axis. It is the positive x-axis sliding surface. It is the negative x-axis sliding surface. It is the positive y-axis sliding surface. It is the negative y-axis sliding surface.

[0134] Based on the determination of the control function, the active control excitation for real-time compensation of the sliding membrane control is calculated:

[0135] ;

[0136] In the formula, F controlx1 It is the active control force in the positive x-axis direction; F controlx2 It is the active control force in the negative x-axis direction; F controly1 It is the active control force in the positive y-axis direction; F controly2 It is the active control force in the negative y-axis direction; M controlx It is the active sliding membrane control bending moment in the x-direction, M controly It is an active sliding membrane control of bending moment in the y-direction; L T It is the axial length of the piezoelectric actuator from the center of the spindle system; d 33 It is the piezoelectric constant; K s It refers to the stiffness of piezoelectric ceramics.

[0137] Step 2: Based on the initial air film conditions, considering the changes in air film under the influence of vibration, derive the dynamic air film force excitation of the bearing;

[0138] Step 2.1: Calculate the instantaneous radial air film force excitation experienced by the spindle during operation;

[0139] according to Figure 2 As can be seen from the embodiments of the present invention provided in (a), (b), and (c), in the radial direction, two internal air bearings provide support to prevent core wear during high-speed rotation of the spindle and to enhance stability. Considering the changes in the internal air film under vibration, the instantaneous radial air film thickness h is calculated. J :

[0140] ;

[0141] In the formula, Δx is the vibration displacement in the radial x-direction; Δy is the vibration displacement in the radial y-direction; Δφ x0 It is the deflection angle in the radial x-direction; Δφ y0 It is the deflection angle in the radial y-direction; It is the initial eccentricity deflection angle; c0 is the film gap of the radial bearing; e0 is the initial eccentricity of the radial bearing; φ x0 It is the initial deflection state in the radial x-direction; φ y0 It represents the initial deflection state in the radial y-direction; θ is the angular position coordinate; z is the axial position coordinate; z m is the axial position coordinate of the midpoint of the main shaft; h0 is the air film thickness considering vibration at the initial position;

[0142] Substituting the obtained gas film thickness into the two-dimensional finite difference equation, the gas film pressure at each point is calculated:

[0143] ;

[0144] In the formula, the dimensionless parameters are defined as follows:

[0145] ;

[0146] In the formula, p a It is atmospheric pressure; c a It is the radial initial film gap; L is the bearing length; p s ρ is the gas supply pressure; μ is the gas viscosity; u1 is the operating speed; p is the actual pressure at each grid point; z is the axial position of each grid point; h is the actual gas film thickness at each grid point.

[0147] Based on the principle of flow conservation, an equilibrium region is defined at the orifice location, and Newton's iterative method is applied within this region to solve for the air pressure at each point. The termination condition for the calculation is set as follows:

[0148] ;

[0149] In the formula, The pressure at node (i, j) after N iterations represents the air film pressure; i and j represent the node coordinates; N is the number of iterations.

[0150] Based on the determined pressure conditions, calculate the instantaneous radial film force excitation experienced by the spindle during operation:

[0151] ;

[0152] In the formula, It is the film force in the radial x-direction; It is the film force in the radial y-direction; It is the bending moment of the air film force in the radial x direction; It is the bending moment of the air film force in the radial y direction; It is the dimensionless air film pressure at each node;

[0153] Step 2.2: Calculate the instantaneous axial film excitation force experienced by the spindle during operation;

[0154] according to Figure 3 As can be seen from the embodiments of the present invention provided in (a), (b), and (c), in the axial direction, due to the use of a bidirectional closed thrust structure inside the air static pressure main shaft, the calculated air film thickness h on both sides of the axial direction is... T1 and h T2 :

[0155] ;

[0156] In the formula, Δz is the vibration displacement in the axial z-direction; Δφ x It is the deflection angle of the thrust bearing in the x-direction; Δφ y It is the deflection angle of the thrust bearing in the y-direction; c T1 and c T2 It is a bidirectional closed thrust film gap, φzx0 and φ zy0 It represents the initial deflection state of the thrust bearing, and r is the radius of the circle containing the differential grid point.

[0157] Substituting the obtained film thickness into the two-dimensional finite difference equation expressed in polar coordinates, the film pressure at each point is calculated:

[0158] ;

[0159] In the formula, the dimensionless parameter is defined as:

[0160] ;

[0161] In the formula, p a It is atmospheric pressure; c T R2 is the initial axial film gap; R2 is the thrust disc radius; p s It is the supply pressure; μ is the gas viscosity; w is the operating angular velocity; w s is the rated angular velocity; p is the actual pressure at each grid point; r is the radius of the circle containing the differential grid point; h is the actual air film thickness at each grid point.

[0162] Considering the boundary conditions of the inner and outer rings, the air film within the bidirectional closed thrust bearing is divided into finite difference grid points. Flow balance elements are set based on each throttling orifice, and the pressure distribution at each point is finally determined using Newton's iteration method, with the termination condition being:

[0163] ;

[0164] In the formula, The gas film pressure at node (r, θ) after N iterations represents the gas film pressure; r and θ represent the polar coordinates of the node; N is the number of iterations.

[0165] Based on the determined pressure conditions, calculate the instantaneous axial film force excitation experienced by the spindle during operation:

[0166] ;

[0167] In the formula, It is the dimensionless gas film pressure at each node. Dimensionless polar radius.

[0168] Step 3: Considering the relative positional relationship between the stator and rotor and the change in the center of mass deflection angle, derive the unbalanced magnetic pull excitation and mass eccentricity excitation under the influence of nonlinear vibration;

[0169] Step 3.1: Calculate the unbalanced magnetic pull excitation experienced by the spindle during operation;

[0170] Based on Kirchhoff's voltage and current laws, the stator and rotor voltages, electromotive force, and core winding excitation are given by the following formula:

[0171] ;

[0172] In the formula, , ,and These are the excitation current, stator current, and rotor phase current, respectively, α Fe It is the stator iron loss angle, φ1 and φ2 are the power factor angles of the stator and rotor, respectively, I em I s and The root mean square values ​​of the excitation current, stator current, and rotor current, respectively. It is the power supply voltage applied to each phase of the stator. It is a single-phase induced electromotive force in the stator. It is the equivalent rotor induced electromotive force, Z s It is the stator leakage reactance, Z em It is the magnetizing impedance. The equivalent impedance of the rotor is expressed as:

[0173] ;

[0174] In the formula, R s and X s These are the stator single-phase resistance and reactance, and These are the rotor equivalent resistance and reactance, R. em and X em These are the motor's excitation resistance and excitation reactance, respectively. R le is the rotor equivalent additional resistance, j is the phase coordinate, and s is the slip of the induction motor.

[0175] Calculate the slip rate of the induction motor through speed analysis: ;

[0176] In the formula, n is the spindle speed, n s =60f cur / p is the synchronous rotational speed of the rotating magnetic field, f cur p and p are the stator current frequency and the number of pole pairs, respectively, and the additional resistance R is... le The power dissipated above represents the mechanical power P generated by the electric spindle. e .

[0177] Considering the actual mechanical output power in the electric spindle, calculate the additional resistance R in the equivalent circuit. le : ;

[0178] An equivalent circuit for the induction motor was established based on the root mean square values ​​of voltage and current. Given the motor parameters and power supply voltage, the excitation current, stator current, and rotor phase current were calculated.

[0179] ;

[0180] Furthermore, the rotor equivalent current is obtained. With the output mechanical power P of the electric spindle e The relationship between them:

[0181] ;

[0182] Solving the above equation yields the current. Then the stator current is obtained. Excitation current and motor slip s;

[0183] Based on the unity of time and space, the rotor, stator, and air gap magnetic field B m The spatial phase relationship of the fundamental magnetomotive force is the same as the temporal phase relationship of the stator current, rotor current, and excitation current. Calculate the three-phase AC current corresponding to the built-in motor:

[0184] ;

[0185] In the formula, ω cur =2πf cur It is the angular velocity of the rotating magnetic field, θ. cur =pα is the electrical angle, p is the pole pair number, and α is the spatial angle. It is the equivalent phase current, i s It is the stator phase current, i em It is the excitation phase current, t is the operating time, and α is the excitation phase current. Fe It is the stator iron damage angle. It is the stator power factor angle. It is the rotor power factor angle. It is the stator space phase angle. This is the stator spatial phase angle. The fundamental magnetomotive force F of the rotor under symmetrical load is calculated using the following formula. rf The fundamental magnetomotive force F of the stator sf The fundamental magnetomotive force F of the air gap m :

[0186] ;

[0187] In the formula, m1=3 is the number of stator phases, N1 is the number of turns of the series coil per pole per phase of the stator, and k w1 It is the stator winding distribution factor. It is the equivalent phase current.

[0188] according to Figure 4 As can be seen from the embodiments of the present invention provided in (a), (b), and (c), in the air gap magnetic field between the rotor and the stator, the normal magnetic flux density plays a dominant role, and its expression is:

[0189] ;

[0190] In the formula, δ n It is the normal magnetic gap length on the rotor surface; It is the normal magnetic flux density on the outer surface of the rotor; It is the air gap permeability per unit area of ​​the rotor's outer surface; μ0 is the vacuum permeability; φ3 is the spatial azimuth angle; f m It is magnetic potential; ω cur It is the angular velocity of the rotating magnetic field; θ cur It is the electric field angle; α is the rotor power factor angle; k is the air gap azimuth angle; c It is the stator winding distribution coefficient.

[0191] When the magnetic material on the rotor comes into contact with the magnetic field in the air gap, a magnetic force is generated. Calculation of conventional magnetic force... :

[0192] ;

[0193] In the formula, It is the magnetic flux density perpendicular to the rotor surface; It is the magnetic pull perpendicular to the rotor surface; α Fe It is the stator iron loss angle; It is the spatial phase of σ.

[0194] 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:

[0195] ;

[0196] In the formula, It is the average 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... si It is the rotor outer diameter, R ro θ is the stator inner diameter, r is the instantaneous linear distance between the rotor cross-section considering vibration and the rotor cross-section not considering vibration, Δx is the vibration displacement in the x-direction, Δy is the vibration displacement in the y-direction, and θ is the vibration displacement in the y-direction. x It is the vibration deflection angle in the x-direction, θ y θ is the vibration deflection angle in the y direction, θ is the dynamic eccentricity azimuth angle, and l is the effective axial length of the rotor.

[0197] Calculate the instantaneous magnetic circuit length considering the vibration response effects in five degrees of freedom. :

[0198] ;

[0199] In the formula, D si =2R si and D ro =2R ro These are the inner diameter of the stator and the outer diameter of the rotor of the built-in motor; (R−R) ro This indicates the change in air gap length caused by the change in the cross-sectional shape of the built-in motor rotor; It is the rotor radial deflection angle. It refers to the axial relative position.

[0200] Based on geometric coordinate relationships, calculate the horizontal, vertical, and axial projected components of the magnetic pull force acting at any point on the rotor surface, as well as the transverse magnetic moment:

[0201] ;

[0202] In the formula, It is the horizontal projection component of the unbalanced magnetic pull. It is the vertical projection component of the unbalanced magnetic pull. It is the axial projection component of the unbalanced magnetic pull. It is the transverse magnetic moment component of the unbalanced magnetic pull.

[0203] Furthermore, by integrating the horizontal, vertical, and axial projection components of the magnetic pull force in the rotor surface normal direction, the projections of the unbalanced magnetic pull force excitation in the X, Y, and Z axes and the components of its additional transverse electromagnetic torque in the X and Y axes are calculated:

[0204] ;

[0205] In the formula, It is an unbalanced magnetic pull effect horizontal excitation. It is an unbalanced magnetic pull effect with vertical excitation. It is an axial excitation caused by unbalanced magnetic pull effect. It is the transverse magnetic moment in the x-direction due to the unbalanced magnetic pull effect. It is the transverse magnetic moment in the y-direction due to the unbalanced magnetic pull effect.

[0206] Step 3.2: Calculate the unbalanced mass eccentric excitation experienced by the spindle during operation;

[0207] according to Figure 5As can be seen from the provided embodiments of the present invention, due to the existence of angular deflection, the centrifugal force caused by the eccentric mass projected along the radial and axial directions of the stator results in an unbalanced mass excitation. Specifically, the projections of the unbalanced mass excitation in the X, Y, and Z axis directions and the components of the additional lateral moment in the X and Y axis directions are calculated:

[0208] ;

[0209] In the formula, m is the mass of the spindle system, w is the angular velocity of the spindle, v is the spatial tilt of the centrifugal force relative to the stator radial direction, and e is the spindle eccentricity. It is the initial eccentricity angle in the x-direction. It is the initial eccentricity angle in the y-direction. It is a quality eccentricity level excitation, It is a mass eccentric vertical excitation. It is axial excitation due to mass eccentricity. It is the eccentric bending moment of mass in the x-direction. It is the eccentric bending moment of mass in the y direction;

[0210] Step 4: Due to the coupling relationship between multiple physics fields, construct a multi-degree-of-freedom dynamic model under complex excitation influence;

[0211] according to Figure 6 As can be seen from the provided embodiments of the present invention, considering the coupling effect of multiple excitation loads, the specific dynamic equations of the multi-degree-of-freedom air hydrostatic master shaft under sliding diaphragm control are as follows:

[0212] ;

[0213] In the formula M s G s and C s These are the mass matrix, gyroscope matrix, and damping matrix of the unit node, respectively, F air It is based on the radial and axial air film force excitation calculated by dynamic air film thickness, F ump It is the unbalanced magnetic excitation generated by the built-in motor under the condition of air gap change, F e It is an unbalanced mass eccentric excitation under the action of instantaneous deflection angle, F control It is an active sliding membrane control excitation determined based on vibration signals, where q is the vibration displacement of each element node. It is the first derivative of the vibration displacement of each element node. It is the second derivative of the vibration displacement of each unit node.

[0214] Based on the nonlinear finite element dynamics model, the displacement and excitation distribution of each node are shown in the following matrix:

[0215] ;

[0216] ;

[0217] ;

[0218] ;

[0219] ;

[0220] In the formula, x1 is the x-direction displacement of node 1 of the mass element in the dynamic model, y1 is the y-direction displacement of node 1 of the mass element in the dynamic model, z1 is the z-direction displacement of node 1 of the mass element in the dynamic model, and θ x1 θ is the x-direction deflection angle of node 1 of the mass element in the dynamic model. y1 It is the y-direction deflection angle of node 1 of the mass element in the dynamic model, x 19 It is the x-direction displacement of node 19 of the mass element in the dynamic model, y 19 It is the y-direction displacement of node 19 of the mass element in the dynamic model, z 19 It is the z-direction displacement of node 19 of the mass element in the dynamic model, θ x19 θ is the x-direction deflection angle of node 19 of the mass element in the dynamic model. y19 It is the y-direction deflection angle of node 19 of the mass element in the dynamic model, F. xL It is the excitation of the gas film force of the left radial gas bearing in the x direction, F yL It is the left radial gas bearing film force excitation in the y direction, M xL M is the left radial gas bearing film force bending moment in the x-direction. yL F is the left radial gas bearing film force bending moment in the y-direction. xR It is the excitation of the film force of the right radial gas bearing in the x-direction, F yR It is the right radial gas bearing film force excitation in the y direction, M xR M is the right radial gas bearing film force bending moment in the x-direction, M. yR F is the right radial gas bearing film force bending moment in the y-direction. zT It is the film force excitation of the thrust gas bearing in the axial direction, M zx M is the bending moment of the thrust gas bearing film force in the x-direction. zy It is the bending moment of the thrust gas bearing film force in the y direction, F controlx It is an active sliding membrane control excitation in the x-direction, F controly It is an active sliding membrane control excitation in the y-direction, M controlx It is the active sliding membrane control bending moment in the x-direction, M controly It is an active sliding membrane control of bending moment in the y direction;

[0221] Step 5: Solve the dynamic equations based on the Newmark-β method to calculate the dynamic response of the step size node, and calculate the active sliding film control excitation, dynamic air film force excitation, unbalanced magnetic pull force excitation and mass eccentricity excitation by the changes caused by the displacement location;

[0222] Step 6: Based on the active sliding membrane control excitation, dynamic air film force excitation, unbalanced magnetic pull force excitation and mass eccentricity excitation, determine the instantaneous load, and couple it with the dynamic equation to solve the real-time state of the displacement field. If the displacement trajectory does not show a convergence trend, repeat steps 1 to 5 until the node changes are stable and converged, and obtain the main shaft displacement field response curve under sliding membrane control for dynamic behavior characteristic analysis.

[0223] This invention breaks through the limitations of traditional linearized, single-physics modeling by constructing a high-fidelity system dynamics model through nonlinear coupling calculations of film force, unbalanced magnetic pull, and eccentric mass excitation. The calculation of film force accurately characterizes the time-varying nonlinearity of bearing stiffness and the generation mechanism of self-excited vibration of the air hammer; the introduction of unbalanced magnetic pull reveals the bidirectional coupling effect between the electromagnetic field and mechanical motion; the eccentric mass, as a classic external synchronous excitation, interacts with the aforementioned two internal excitations, potentially generating complex combined frequency vibrations. The synergistic calculation of these three factors enables the model to reproduce, for the first time, the nonlinear phenomena such as jumping and bifurcation that occur in the spindle under high-speed, variable-speed, and load-changing conditions, laying a physical foundation for accurately predicting the stability boundary and dynamic accuracy limit of the system.

[0224] Based on this, this invention provides an indispensable and precise platform for the design and verification of sliding mode control. Since sliding mode control is highly robust to system model uncertainties and external disturbances, a precise coupled model encompassing multi-source excitations is crucial for its successful application. This invention not only allows controller design to pre-consider and compensate for the nonlinear time-varying nature of film force, the state-dependent disturbances of magnetic pull, and eccentric force excitation, but also to evaluate the inverse effects of the control action itself on film stability and the electromagnetic field, achieving true "mechanical-electrical-pneumatic-control" synergistic optimization. Therefore, this invention suppresses multi-source coupled vibration through active control, providing theoretical reference and simulation prediction for the active vibration suppression coupled nonlinear analysis of precision machining of air hydrostatic spindles.

[0225] Example 2:

[0226] This embodiment proposes an electronic device, including: one or more processors, and a memory, wherein the memory is used to store instructions, and when the instructions are executed by the one or more processors, the one or more processors execute the nonlinear vibration coupling analysis method of an air hydrostatic spindle under sliding diaphragm control.

[0227] The electronic device can be a mobile phone, computer, or tablet computer, etc., and includes a memory and a processor. The memory stores a computer program, which, when executed by the processor, implements a nonlinear vibration coupling analysis method for an air hydrostatic spindle under sliding diaphragm control as described in the embodiment. It is understood that the electronic device may also include an input / output (I / O) interface and communication components.

[0228] The processor is used to execute all or part of the steps in the nonlinear vibration coupling analysis method for an air hydrostatic spindle under sliding diaphragm control as described in the above embodiments. The memory is used to store various types of data, which may include, for example, instructions for any application or method in an electronic device, as well as application-related data.

[0229] The processor can be implemented as an Application Specific Integrated Circuit (ASIC), Digital Signal Processor (DSP), Programmable Logic Device (PLD), Field Programmable Gate Array (FPGA), controller, microcontroller, microprocessor, or other electronic components, and is used to execute the nonlinear vibration coupling analysis method for air hydrostatic spindle under sliding diaphragm control described in the above embodiments.

[0230] Example 3:

[0231] This embodiment proposes a computer-readable storage medium that stores executable instructions. When these instructions are executed, if they are implemented as software functional units and sold or used as independent products, they can be stored in a computer-readable storage medium.

[0232] The computer software product is stored in a storage medium and includes several instructions to cause a computer device (which may be a personal computer, a server, or a network device, etc.) to execute all or part of the steps of the nonlinear vibration coupling analysis method for an air hydrostatic spindle under sliding diaphragm control as described in the various embodiments of this application.

[0233] The aforementioned storage media include: flash memory, hard disk, multimedia card, card-type memory (e.g., SD (Secure Digital Memory Card) or DX (Memory Data Register, MDR) memory, random access memory (RAM), static random access memory (SRAM), read-only memory (ROM), electrically erasable programmable read-only memory (EEPROM), programmable read-only memory (PROM), magnetic memory, disk, optical disk, server, APP (Application) application store, and other media capable of storing program verification codes. These media store computer programs, which, when executed by a processor, can implement the various steps of the aforementioned nonlinear vibration coupling analysis method for an air hydrostatic spindle under sliding diaphragm control.

[0234] Example 4:

[0235] This embodiment proposes a computer program product, including a computer program or instructions, which, when executed by a processor, implements the aforementioned nonlinear vibration coupling analysis method for an air hydrostatic spindle under sliding diaphragm control.

[0236] Based on this understanding, the technical solution of this application, in essence, or the part that contributes to the prior art, or part of the technical solution, can be embodied in the form of a computer program product.

[0237] The various embodiments in this application are described in a progressive manner. The same or similar parts between the various embodiments can be referred to each other. Each embodiment focuses on describing the differences from other embodiments.

[0238] The scope of protection of this application is not limited to the embodiments described above. Obviously, those skilled in the art can make various modifications and variations to this disclosure without departing from the scope and spirit of this disclosure. If such modifications and variations fall within the scope of the methods disclosed herein and their equivalents, then the intent of this disclosure also includes such modifications and variations.

Claims

1. A method for coupled analysis of nonlinear vibration of an air hydrostatic principal shaft under sliding diaphragm control, characterized in that, Includes the following steps: Step 1: Construct an air static pressure spindle system and obtain real-time active control excitation under sliding diaphragm control based on the position error signal and the closed-loop feedback principle; Step 2: Based on the initial air film conditions, considering the changes in air film under the influence of vibration, derive the dynamic air film force excitation of the bearing; Step 3: Considering the relative positional relationship between the stator and rotor and the change in the center of mass deflection angle, derive the unbalanced magnetic pull excitation and mass eccentricity excitation under the influence of nonlinear vibration; Step 4: Due to the coupling relationship between multiple physics fields, construct a multi-degree-of-freedom dynamic model under complex excitation influence; Step 5: Solve the dynamic equations based on the Newmark-β method to calculate the dynamic response of the step size node, and calculate the active sliding film control excitation, dynamic air film force excitation, unbalanced magnetic pull force excitation and mass eccentricity excitation by the changes caused by the displacement location; Step 6: Based on the active sliding membrane control excitation, dynamic air film force excitation, unbalanced magnetic pull force excitation and mass eccentricity excitation, determine the instantaneous load, and couple it with the dynamic equation to solve the real-time state of the displacement field. If the displacement trajectory does not show a convergence trend, repeat steps 1 to 5 until the node changes are stable and converged, and obtain the main shaft displacement field response curve under sliding membrane control for dynamic behavior characteristic analysis.

2. The nonlinear vibration coupling analysis method for an air hydrostatic master shaft under sliding diaphragm control according to claim 1, characterized in that, Step 1 involves customizing a dedicated tool holder fixed to the tool shank based on the closed-loop control of the air-static spindle. This dedicated tool holder integrates a pair of deep groove ball bearings, four piezoelectric actuators, and a bolt-fixed tool holder housing. The inner ring of the deep groove ball bearing rotates at high speed with the tool, while the outer ring contacts the piezoelectric actuators to transmit real-time control force. Two sets of radially arranged displacement sensors detect the dynamic amplitude signal of the spindle. The acquisition card converts the analog signal into digital data and transmits it to the computer. The computer calculates the compensation voltage in real time through an active sliding diaphragm control program. The calculated compensation voltage digital signal is converted into an analog signal and sent to the piezoelectric controller. After signal amplification, the signal is applied to the four piezoelectric actuators to achieve closed-loop instantaneous active compensation control.

3. The nonlinear vibration coupling analysis method for an air hydrostatic master shaft under sliding diaphragm control according to claim 2, characterized in that, Step 1 specifically involves: Based on the core concept of sliding membrane control, the sliding membrane surface s is calculated using the following formula: ; In the formula, λ is the speed error, λ is the slope of the sliding surface, and e is the position error; Calculate the control function u of the active slicker control program: ; In the formula, b is the control gain, and f is the nonlinear function of the air static pressure spindle system. It's a speed error. λ is the second derivative of the desired position, K is the set gain, s is the sliding surface, λ is the sliding surface parameter, u is the control gain, and sgn is the switching control law. Considering the four mutually perpendicular radial directions, calculate the control gain u experienced by the piezoelectric actuator in each of the four directions. x1 u x2 u y1 and u y2 : ； In the formula, u x1 It is the active control gain in the positive x-axis direction; u x2 It is the active control gain in the negative x-axis direction; u y1 It is the active control gain in the positive y-axis direction; u y2 It is the active control gain in the negative y-axis direction, k x1 It is the stiffness coefficient in the positive x-axis direction, k x2 It is the stiffness coefficient in the negative x-axis direction, k y1 It is the stiffness coefficient in the positive y-axis direction, k y2 It is the stiffness coefficient in the negative y-axis direction, c x1 It is the damping coefficient in the positive x-axis direction, c x2 It is the damping coefficient in the negative x-axis direction, c y1 It is the damping coefficient in the positive y-axis direction, c y2 It is the damping coefficient in the negative y-axis direction. It is the second derivative of the desired position in the positive x-axis direction. It is the second derivative of the desired position in the negative x-axis direction. It is the second derivative of the desired position in the positive y-axis direction. It is the second derivative of the desired position in the negative y-axis direction. It is the positive velocity error along the x-axis. It is the negative velocity error along the x-axis. It is the positive y-axis velocity error. It is the negative velocity error along the y-axis. It is the positive x-axis sliding surface. It is the negative x-axis sliding surface. It is the positive y-axis sliding surface. It is the negative y-axis sliding surface; Based on the determination of the control function, the active control excitation for real-time compensation of the sliding membrane control is calculated: ； In the formula, F controlx1 It is the active control force in the positive x-axis direction; F controlx2 It is the active control force in the negative x-axis direction; F controly1 It is the active control force in the positive y-axis direction; F controly2 It is the active control force in the negative y-axis direction; M controlx It is the active sliding membrane control bending moment in the x-direction, M controly It is an active sliding membrane control of bending moment in the y-direction; L T It is the axial length of the piezoelectric actuator from the center of the spindle system; d 33 It is the piezoelectric constant; K s It refers to the stiffness of piezoelectric ceramics.

4. The nonlinear vibration coupling analysis method for an air hydrostatic master shaft under sliding diaphragm control according to claim 1, characterized in that, Step 2 specifically includes: Step 2.1: Calculate the instantaneous radial air film force excitation experienced by the spindle during operation; Calculate the radial instantaneous air film thickness h J : ； In the formula, Δx is the vibration displacement in the radial x-direction; Δy is the vibration displacement in the radial y-direction; Δφ x0 It is the deflection angle in the radial x-direction; Δφ y0 It is the deflection angle in the radial y-direction; It is the initial eccentricity deflection angle; c0 is the film gap of the radial bearing; e0 is the initial eccentricity of the radial bearing; φ x0 It is the initial deflection state in the radial x-direction; φ y0 It represents the initial deflection state in the radial y-direction; θ is the angular position coordinate; z is the axial position coordinate; z m is the axial position coordinate of the midpoint of the main shaft; h0 is the air film thickness considering vibration at the initial position; Substituting the obtained gas film thickness into the two-dimensional finite difference equation, the gas film pressure at each point is calculated: ； In the formula, the dimensionless parameters are defined as follows: ； In the formula, p a It is atmospheric pressure; c a It is the radial initial film gap; L is the bearing length; p s ρ is the gas supply pressure; μ is the gas viscosity; u1 is the operating speed; p is the actual pressure at each grid point; z is the axial position of each grid point; h is the actual gas film thickness at each grid point. Based on the principle of flow conservation, an equilibrium region is defined at the location of the throttling orifice, and the Newton-Raphson iteration method is applied within this region to solve for the air pressure at each point; the termination condition for the calculation is set as follows: ； In the formula, The pressure at node (i, j) after N iterations represents the air film pressure; i and j represent the node coordinates; N is the number of iterations. Based on the determined pressure conditions, calculate the instantaneous radial film force excitation experienced by the spindle during operation: ； In the formula, It is the film force in the radial x-direction; It is the film force in the radial y-direction; It is the bending moment of the air film force in the radial x direction; It is the bending moment of the air film force in the radial y direction; It is the dimensionless air film pressure at each node; Step 2.2: Calculate the instantaneous axial film excitation force experienced by the spindle during operation; Calculate the thickness h of the air film on both sides of the axis. T1 and h T2 : ； In the formula, Δz is the vibration displacement in the axial z-direction; Δφ x It is the deflection angle of the thrust bearing in the x-direction; Δφ y It is the deflection angle of the thrust bearing in the y-direction; c T1 and c T2 It is a bidirectional closed thrust film gap, φ zx0 and φ zy0 It represents the initial deflection state of the thrust bearing, where r is the radius of the circle containing the differential grid point; Substituting the obtained film thickness into the two-dimensional finite difference equation expressed in polar coordinates, the film pressure at each point is calculated: ； In the formula, the dimensionless parameter is defined as: ； In the formula, p a It is atmospheric pressure; c T R2 is the initial axial film gap; R2 is the thrust disc radius; p s It is the supply pressure; μ is the gas viscosity; w is the operating angular velocity; w s is the rated angular velocity; p is the actual pressure at each grid point; r is the radius of the circle containing the differential grid point; h is the actual air film thickness at each grid point; Considering the boundary conditions of the inner and outer rings, the air film inside the bidirectional closed thrust bearing is divided into finite difference grid points; a flow balance element is set based on each throttling orifice, and the pressure distribution at each point is finally determined by Newton's iteration method, with the termination condition being: ； In the formula, The gas film pressure at node (r, θ) after N iterations represents the gas film pressure; r and θ represent the polar coordinates of the node; N is the number of iterations. Based on the determined pressure conditions, calculate the instantaneous axial film force excitation experienced by the spindle during operation: ； In the formula, It is the dimensionless gas film pressure at each node. Dimensionless polar radius.

5. The nonlinear vibration coupling analysis method for an air hydrostatic master shaft under sliding diaphragm control according to claim 1, characterized in that, Step 3 specifically includes: Step 3.1: Calculate the unbalanced magnetic pull excitation experienced by the spindle during operation; Step 3.2: Calculate the unbalanced mass eccentric excitation experienced by the spindle during operation; The centrifugal force caused by the eccentric mass is projected along the radial and axial directions of the stator, resulting in unbalanced mass excitation. Specifically, the projections of the unbalanced mass excitation onto the X, Y, and Z axes, and the components of the additional lateral moment in the X and Y axes, are calculated: ； In the formula, m is the mass of the spindle system, w is the angular velocity of the spindle, v is the spatial tilt of the centrifugal force relative to the stator radial direction, and e is the spindle eccentricity. It is the initial eccentricity angle in the x-direction. It is the initial eccentricity angle in the y-direction. It is a quality eccentricity level excitation, It is a mass eccentric vertical excitation. It is axial excitation due to mass eccentricity. It is the eccentric bending moment of mass in the x-direction. It is the eccentric bending moment of mass in the y direction.

6. The nonlinear vibration coupling analysis method for an air hydrostatic master shaft under sliding diaphragm control according to claim 5, characterized in that, Step 3.1 specifically involves: Based on Kirchhoff's voltage and current laws, the stator and rotor voltages, electromotive force, and core winding excitation are given by the following formulas: ； In the formula, , ,and These are the excitation current, stator current, and rotor phase current, respectively, α Fe It is the stator iron loss angle, φ1 and φ2 are the power factor angles of the stator and rotor, respectively, I em I s and The root mean square values ​​of the excitation current, stator current, and rotor current, respectively. It is the power supply voltage applied to each phase of the stator. It is a single-phase induced electromotive force in the stator. It is the equivalent rotor induced electromotive force, Z s It is the stator leakage reactance, Z em It is the magnetizing impedance. The equivalent impedance of the rotor is expressed as: ； In the formula, R s and X s These are the stator single-phase resistance and reactance, and These are the rotor equivalent resistance and reactance, R. em and X em These are the motor's excitation resistance and excitation reactance; R le is the rotor equivalent additional resistance, j is the phase coordinate, and s is the slip of the induction motor; Calculate the slip rate of the induction motor through speed analysis: ; In the formula, n is the spindle speed, n s =60f cur / p is the synchronous rotational speed of the rotating magnetic field, f cur p and p are the stator current frequency and the number of pole pairs, respectively, and the additional resistance R is... le The power dissipated above represents the mechanical power P generated by the electric spindle. e ; Considering the actual mechanical output power in the electric spindle, calculate the additional resistance R in the equivalent circuit. le : ; An equivalent circuit for the induction motor was established based on the root mean square values ​​of voltage and current. Given the motor parameters and power supply voltage, the excitation current, stator current, and rotor phase current were calculated. ； Furthermore, the rotor equivalent current is obtained. With the output mechanical power P of the electric spindle e The relationship between them: ； Solving the above equation yields the current. Then the stator current is obtained. Excitation current and motor slip s; Calculate the three-phase AC power required for the built-in motor: ； In the formula, ω cur =2πf cur It is the angular velocity of the rotating magnetic field, θ. cur =pα is the electrical angle, p is the pole pair number, and α is the spatial angle. It is the equivalent phase current, i s It is the stator phase current, i em It is the excitation phase current, t is the operating time, and α is the excitation phase current. Fe It is the stator iron damage angle. It is the stator power factor angle. It is the rotor power factor angle. It is the stator space phase angle. It is the stator spatial phase angle; the fundamental magnetomotive force F of the rotor under symmetrical load is calculated by the following formula. rf The fundamental magnetomotive force F of the stator sf The fundamental magnetomotive force F of the air gap m : ； In the formula, m1=3 is the number of stator phases, N1 is the number of turns of the series coil per pole per phase of the stator, and k w1 It is the stator winding distribution factor. It is the equivalent phase current; In the air gap magnetic field between the rotor and stator, the expression for the normal magnetic flux density is: ； In the formula, δ n It is the normal magnetic gap length on the rotor surface; It is the normal magnetic flux density on the outer surface of the rotor; It is the air gap permeability per unit area of ​​the rotor's outer surface; μ0 is the vacuum permeability; φ3 is the spatial azimuth angle; f m It is magnetic potential; ω cur It is the angular velocity of the rotating magnetic field; θ cur It is the electric field angle; α is the rotor power factor angle; k is the air gap azimuth angle; c It is the stator winding distribution factor; When the magnetic material on the rotor comes into contact with the magnetic field in the air gap, a magnetic force is generated. Calculation of conventional magnetic force... : ； In the formula, It is the magnetic flux density perpendicular to the rotor surface; It is the magnetic pull perpendicular to the rotor surface; α Fe It is the stator iron loss angle; It is the spatial phase of σ; 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: ； In the formula, It is the average 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... si It is the rotor outer diameter, R ro θ is the stator inner diameter, r is the instantaneous linear distance between the rotor cross-section considering vibration and the rotor cross-section not considering vibration, Δx is the vibration displacement in the x-direction, Δy is the vibration displacement in the y-direction, and θ is the vibration displacement in the y-direction. x It is the vibration deflection angle in the x-direction, θ y θ is the vibration deflection angle in the y direction, θ is the dynamic eccentricity azimuth angle, and l is the effective axial length of the rotor. Calculate the instantaneous magnetic circuit length considering the vibration response effects in five degrees of freedom. : ； In the formula, D si =2R si and D ro =2R ro These are the inner diameter of the stator and the outer diameter of the rotor of the built-in motor; (R−R) ro This indicates the change in air gap length caused by the change in the cross-sectional shape of the built-in motor rotor; It is the rotor radial deflection angle. It refers to the axial relative position; Based on geometric coordinate relationships, calculate the horizontal, vertical, and axial projected components of the magnetic pull force acting at any point on the rotor surface, as well as the transverse magnetic moment: ； In the formula, It is the horizontal projection component of the unbalanced magnetic pull. It is the vertical projection component of the unbalanced magnetic pull. It is the axial projection component of the unbalanced magnetic pull. It is the transverse magnetic moment component of the unbalanced magnetic pull; Furthermore, by integrating the horizontal, vertical, and axial projection components of the magnetic pull force in the rotor surface normal direction, the projections of the unbalanced magnetic pull force excitation in the X, Y, and Z axes and the components of its additional transverse electromagnetic torque in the X and Y axes are calculated: ； In the formula, It is an unbalanced magnetic pull effect horizontal excitation. It is an unbalanced magnetic pull effect with vertical excitation. It is an axial excitation caused by unbalanced magnetic pull effect. It is the transverse magnetic moment in the x-direction due to the unbalanced magnetic pull effect. It is the transverse magnetic moment in the y-direction due to the unbalanced magnetic pull effect.

7. The nonlinear vibration coupling analysis method for an air hydrostatic master shaft under sliding diaphragm control according to claim 6, characterized in that, Step 4 specifically involves: Considering the coupling effect of multiple excitation loads, the multi-degree-of-freedom dynamic equations of the air hydrostatic master shaft under sliding film control are as follows: ； In the formula M s G s and C s These are the mass matrix, gyroscope matrix, and damping matrix of the unit node, respectively, F air It is based on the radial and axial air film force excitation calculated by dynamic air film thickness, F ump It is the unbalanced magnetic excitation generated by the built-in motor under the condition of air gap change, F e It is an unbalanced mass eccentric excitation under the action of instantaneous deflection angle, F control It is an active sliding membrane control excitation determined based on vibration signals, where q is the vibration displacement of each element node. It is the first derivative of the vibration displacement of each element node. It is the second derivative of the vibration displacement of each unit node; Based on the nonlinear finite element dynamics model, the displacement and excitation distribution of each node are shown in the following matrix: ； ； ； ； ； In the formula, x1 is the x-direction displacement of node 1 of the mass element in the dynamic model, y1 is the y-direction displacement of node 1 of the mass element in the dynamic model, z1 is the z-direction displacement of node 1 of the mass element in the dynamic model, and θ x1 θ is the x-direction deflection angle of node 1 of the mass element in the dynamic model. y1 It is the y-direction deflection angle of node 1 of the mass element in the dynamic model, x 19 It is the x-direction displacement of node 19 of the mass element in the dynamic model, y 19 It is the y-direction displacement of node 19 of the mass element in the dynamic model, z 19 It is the z-direction displacement of node 19 of the mass element in the dynamic model, θ x19 θ is the x-direction deflection angle of node 19 of the mass element in the dynamic model. y19 It is the y-direction deflection angle of node 19 of the mass element in the dynamic model, F. xL It is the excitation of the gas film force of the left radial gas bearing in the x direction, F yL It is the left radial gas bearing film force excitation in the y direction, M xL M is the left radial gas bearing film force bending moment in the x-direction. yL F is the left radial gas bearing film force bending moment in the y-direction. xR It is the excitation of the film force of the right radial gas bearing in the x-direction, F yR It is the right radial gas bearing film force excitation in the y direction, M xR M is the right radial gas bearing film force bending moment in the x-direction, M. yR F is the right radial gas bearing film force bending moment in the y-direction. zT It is the film force excitation of the thrust gas bearing in the axial direction, M zx M is the bending moment of the thrust gas bearing film force in the x-direction. zy It is the bending moment of the thrust gas bearing film force in the y direction, F controlx It is an active sliding membrane control excitation in the x-direction, F controly It is an active sliding membrane control excitation in the y-direction, M controlx It is the active sliding membrane control bending moment in the x-direction, M controly It is an active sliding membrane control of bending moment in the y direction.

Images