Rapid calculation method for electromagnetic vibration of permanent magnet synchronous motor

Abstract

The invention discloses a rapid calculation method for electromagnetic vibration of a permanent magnet synchronous motor, and belongs to the technical field of electromagnetic field calculation. According to the method, a magnetic potential and magnetic conductance method is innovatively combined with an improved sub-domain method, the problem of complexity of electromagnetic field distribution of a nonlinear material is effectively solved through the magnetic potential and magnetic conductance method, and the capture precision of electromagnetic field changes is remarkably improved by dividing a computational domain into a plurality of sub-domains through the improved sub-domain method. According to the method, the influence of the tangential electromagnetic force wave and the harmonic current is fully considered, and high-precision calculation of the electromagnetic force can be realized. On the basis, a stator core shell equivalent theoretical model is established, the inherent frequency of the stator is calculated by adopting an energy method, and the accuracy of the model is verified through a modal experiment. Furthermore, by determining the frequency response function of the stator in the low-order mode, the acceleration spectrum and electromagnetic force wave distribution of the prototype stator in each mode are obtained based on the linear superposition principle. Compared with a traditional electromagnetic-structure coupling finite element model, the method has the advantages that the calculation efficiency is remarkably improved while the calculation precision is guaranteed, the calculation time is shortened by 59.3%, and efficient and reliable technical support is provided for optimization design of electromagnetic equipment.

Description

TECHNICAL FIELD

[0001] The present application relates to the technical field of permanent magnet synchronous motor (PMSM) electromagnetic vibration calculation, in particular to a method for quickly calculating electromagnetic vibration of PMSM stator tooth surface by comprehensively considering the influence of harmonic current and tangential electromagnetic force. BACKGROUND

[0002] With the continuous expansion of the application field of permanent magnet synchronous motor (PMSM) and the increasing requirements of electric drive system on performance and environment, the vibration and noise problem of the motor has attracted widespread attention from domestic and foreign scholars. In the application scene with strict requirements on vibration and noise, the vibration and noise performance has become a key technical index for measuring the comprehensive performance of PMSM. Electromagnetic vibration, as an important index for evaluating the performance of the motor, is generated by the electromagnetic force acting on the stator core. Therefore, the fast and accurate calculation of electromagnetic force and stator natural frequency is the basis for predicting and reducing electromagnetic vibration.

[0003] When studying the electromagnetic vibration and noise of PMSM, the electromagnetic force in the air gap needs to be analyzed based on the magnetic field model, and the magnetic field modeling and solving are the basis for electromagnetic force analysis. Common magnetic field modeling methods include analytical method, numerical method and semi-analytical method, among which the magnetic potential-permeability method and sub-domain model method in the analytical method are widely used.

[0004] Although the equivalent magnetic circuit method can consider factors such as core saturation, slot effect, leakage, etc., and has high calculation accuracy, the accuracy depends on the number of nodes. For example, in 1985, Rahman M.A. proposed an equivalent magnetic circuit method for solving the magnetic field of embedded PMSM, which only considered the conductivity of permanent magnet and air gap, and the accuracy was insufficient. Hemeida used the magnetic equivalent circuit (MEC) to accurately solve the fine parts of complex geometry in a simple way, and verified the overall quantity of the motor by the finite element model results. Forstner proposed a motor modeling framework based on MEC model, which used network theory to systematize and automate the modeling task, and could be widely applied to various electromagnetic actuators.

[0005] The sub-domain modeling technology was proposed by Gysen B.L.J, and later used by L.J Wu and Z.Q Zhu for electromagnetic field analysis of surface-mounted PMSM. This method is simple and fast, but ignores the core saturation and is not suitable for high saturation conditions. L.J. Wu and Z.Q. Zhu proposed an improved analytical sub-domain model, which considered the tooth tip part of the slot effect, and divided the whole field into four simple sub-domains for solving. Z.Li, X.Huang, etc. proposed a hybrid sub-domain model (HSM), which considered the slot effect, rotor asymmetry and ferromagnetic nonlinearity, and was used to calculate the magnetic field distribution of internal permanent magnet motor. Kyung-Hun Shin carried out electromagnetic analysis on PMSM based on the sub-domain method, and considered the shaft leakage and end winding leakage inductance for experimental verification.

[0006] In the calculation of stator natural frequency, the existing research often uses the finite element method (FEM) to establish a three-dimensional model of the stator core for calculation. Although the calculation is accurate, the modeling process is complex and not suitable for rapid calculation and parameter adjustment. Existing analytical methods are mostly based on single-loop electromechanical analogy, often ignoring the stiffness of the stator teeth and windings, resulting in reduced accuracy. Based on three-dimensional elastic theory to calculate the natural frequency of cylindrical shell, the axial, circumferential and coupled vibration characteristics can be considered comprehensively, but the modeling process is complex and time-consuming due to multi-physical field coupling. Some analytical methods do not consider the influence of the casing, which differs from the actual operation of the motor.

[0007] In summary, the existing theoretical methods have the shortcomings of ignoring the influence of the stator core and the casing on the actual operating conditions, and the complexity and time-consuming of finite element modeling. Therefore, it is of great theoretical significance and engineering application value to develop a method for calculating the electromagnetic vibration of the stator teeth of a permanent magnet synchronous motor that can balance calculation accuracy and efficiency and consider the influence of various factors. SUMMARY

[0008] The purpose of the present application is to provide a method for quickly calculating the electromagnetic vibration of the stator teeth of a permanent magnet synchronous motor, to solve the problems of insufficient accuracy in dealing with the complexity of electromagnetic field distribution of nonlinear materials, considering the influence of tangential electromagnetic force wave and harmonic current, and low efficiency in calculating the natural frequency and electromagnetic vibration of the stator.

[0009] The technical solution of the present application is as follows:

[0010] The electromagnetic vibration fast calculation method of the permanent magnet synchronous motor comprises the following steps:

[0011] 1) Calculate the electromagnetic force by combining the magnetic potential permeability method with the improved subdomain method;

[0012] 2) Establish an equivalent theoretical model of the stator core and the casing, and calculate the natural frequency of the stator using the energy method;

[0013] 3) Introduce the electromagnetic force into the structure model in the form of node force, and calculate the electromagnetic vibration using the modal superposition method;

[0014] 4) Verify the calculation results of the stator natural frequency to ensure the accuracy of the calculation results.

[0015] Further, the specific process of step 1 is as follows:

[0016] S1.1 Calculation of air gap magnetic motive force of permanent magnet

[0017] The slotless motor model is divided into two subdomains: PM subdomain I and gap subdomain II as shown in Figure 1 F1(θ) and B o(θ) are the air-gap MMF and flux density distributions respectively, δ is the air-gap length, h pm F1(θ) is calculated as:

[0018]

[0019] F1(θ) = F 1r (θ) + jF 1t (θ) (2)

[0020] B o (θ) = B or (θ) + B ot (θ) (3)

[0021] where F 1r (θ), F 1t (θ), B or (θ) and B ot (θ) are the radial and tangential components of F1(θ) and B o (θ), j is an imaginary operator; F 1r (θ) and F 1t (θ) have the same harmonic content, but their amplitudes and phases are different;

[0022] The air-gap flux density B o (θ) of the slotless machine is obtained by using the subdomain method, i.e. by dividing the above-mentioned subdomains, combining the boundary conditions and Maxwell's equations to solve the electromagnetic field distribution;

[0023] The PM air-gap MMF F1(θ, t) is expressed as:

[0024]

[0025] where ω r is the angular speed, p is the pole number, F 1r_μ and F 1t_μ are the μth harmonic amplitudes of the radial F 1r (θ, t) and tangential F 1t (θ, t) components of F1(θ, t), which are obtained by Fourier decomposition of F 1r (θ) and F 1t (θ);

[0026] S1.2 Calculation of air-gap permeability:

[0027] The basis for dividing the subdomains is:

[0028] a) the relative permeability of the fixed and rotor ferromagnetic regions is infinite;

[0029] b) the PM has linearized characteristics, and its relative permeability μ r = 1;

[0030] c) ignoring the end effect of motor shaft;

[0031] d) the stator slot of motor is equivalent to a slot;

[0032] The model of slotted motor is divided into four sub-domains: sub-domain I is permanent magnet PM, sub-domain II is air gap, sub-domain III is slot opening, and sub-domain IV is slot body, as shown in Figure 2 By applying surface current on both sides of the stator tooth, the saturation effect in the tooth is effectively simulated by using iterative formula, and the surface current generated by the magnetic flux is opposite to the main magnetic flux, and it is generally considered that the saturation on both sides of the slot opening is in a vacuum state, considering the surface current on both sides of the opening, the vector potential A z3i in sub-domain III is represented as:

[0033]

[0034] In the formula, A z3i : the vector potential of the i-th slot in sub-domain III (axial component, unit: Wb / m); A 3i0 : the constant component of the vector potential, reflecting the direct current or average effect of the magnetic field; A 3im : the m-th spatial harmonic amplitude of the vector potential, representing the periodic fluctuation of the magnetic field; F n : the spatial harmonic number, defined as F n =nπ / b oa , related to the slot width angle; θ: circumferential coordinate (radian); θ zi : the central position angle of the i-th slot (radian), used to locate the harmonic spatial phase; R s : the inner radius of the stator; R t : the top radius of the stator slot, b oa : the slot width angle, Q 3i and D 3i are the coefficients to be determined;

[0035] The amplitude of the current density constant component j 3io and the n-th harmonic component J 3im on both sides of the i-th slot is represented as:

[0036]

[0037] Where S c is the area of each surface current;

[0038] The vector potential A z4i in sub-domain IV can be directly derived;

[0039] By solving the coupled equations defined by the boundary conditions of each sub-domain, the Q 3i , D3i and the value of the pending coefficients, by analyzing the air-gap flux density along the circumference B δt (θ) distribution within a specified tooth pitch range, the air-gap permeability Λ(θ) within the tooth pitch can be calculated.

[0040]

[0041] where B δt_r (θ) and B δt_t (θ) are the radial and tangential components of Bδt(θ), which can be obtained by the subdomain method. Λ r (θ), Λ t (θ) are the radial and tangential components of Λ(θ), which can be obtained as follows.

[0042]

[0043] The air-gap permeability Λ(θ) can also be expressed as:

[0044]

[0045] where Λ r_0 is the constant component of Λ r (θ), Λ r_n and Λ t_n are the n-th harmonic amplitudes of Λ r (θ) and Λ t (θ), respectively. These amplitudes are obtained from the Fourier decomposition of Λ r (θ) and Λ t (θ).

[0046] S1.3 Calculation of the armature reaction magnetomotive force

[0047] When magnetic saturation is considered, the permeability varies nonlinearly. From equation (1), it can be seen that the flux generated by the permanent magnet is nonlinear with respect to the current. To simplify the analysis in this section, we consider the saturation effect of magnetic saturation with respect to the basic permanent flux (ψ(i d ,i q )). Thus, the relationship between the d-axis flux and the current is given by:

[0048]

[0049] where: ψ d : d-axis flux linkage, which is composed of the flux linkage generated by the d-axis current id and the inductance Ld, superimposed with the harmonic components;

[0050] ψ q : q-axis flux linkage, which is composed of the flux linkage generated by the q-axis current iq and the inductance Lq, superimposed with the harmonic components;

[0051] Ld ,L q : d-axis and q-axis inductance, representing the linear relationship between current and flux linkage;

[0052] ψ1: fundamental flux linkage component, generated by permanent magnet or main magnetic field;

[0053] ψ3, ψ5, ψ7, ψ9: 3rd, 5th, 7th, 9th harmonic flux linkage components, caused by magnetic field spatial harmonics;

[0054] cos6θ, cos12θ, sin6θ, sin12θ: spatial harmonic functions, θ is the mechanical angle, the coefficients 6 and 12 correspond to the combined effect of pole pair number and slot number (such as 6-pole motor or slot harmonics);

[0055]

[0056] where: ψ dk : d-axis k-th harmonic no-load flux linkage, representing the k-th harmonic component in the d-axis direction of the permanent magnet magnetic field, such as 6th, 12th harmonic, etc.

[0057] ψ qk : q-axis k-th harmonic no-load flux linkage, representing the k-th harmonic component in the q-axis direction of the permanent magnet magnetic field, which acts together with the d-axis harmonic component to produce torque ripple;

[0058] T em : electromagnetic torque, is the core physical quantity of energy conversion of motor, representing the torque generated by the interaction of magnetic field and current;

[0059] transpose or equivalent expression of no-load flux linkage fundamental component (ψ i-1 ), representing the coupling effect of the fundamental flux linkage generated by the permanent magnet and the current;

[0060] p: number of motor pole pairs;

[0061] As can be seen from the expression, ψ 2i-1 represents the no-load flux linkage of harmonic component (2i-1). Where R is the phase resistance, L d and L q are phase inductances, both are equal in surface-mounted PMSM (L d = L q ). Parameters ω r , J, T load and D fcorrespond to the mechanical angular velocity, the moment of inertia, the load torque and the damping coefficient, respectively. By establishing the analytical model of the surface-mounted permanent magnet synchronous motor vector control system based on MATLAB space vector pulse width modulation (SVPWM), the inverter output voltage u d and u q are derived. The fourth-order Runge-Kutta method is used to solve the coupled equations of flux, voltage, electromagnetic torque and motion dynamics, so as to determine the phase current i d and i q .

[0062] The calculation model of the armature magnetic motive force is shown in Figure 3 , which only includes the air gap sub-domain II. The center line of each slot corresponds to a surface current J i uniformly distributed along the inner surface of the stator core, which is represented as a line current in a two-dimensional field. The tangential component of the magnetic field intensity H θ at the inner and outer boundaries of the air gap sub-domain II satisfies specific boundary conditions.

[0063]

[0064] where I i is the phase current of the i-th slot, S c is the single-slot surface current distribution area. When the current on both sides of the slot (I 1i = I 2i = I i / 2) is half of the phase current divided by the slot area.

[0065]

[0066] where n is the harmonic order, and sin(nπ / 2)) reflects the non-zero characteristics of odd harmonics (n = 1, 3, 5, …), and the harmonic component amplitude is proportional to the phase current (I i ) and inversely proportional to the harmonic order n.

[0067] ω el = lω e = (6n+1)ω e (n = 0, 1, 2, …) (21)

[0068] where ω e is the fundamental angular frequency (corresponding to the synchronous speed of the motor).

[0069] (l = 6n ± 1) represents the low-order harmonic order, such as (n = 0) (l = 1) (fundamental frequency), (n = 1) (l = 5, 7…) (5th, 7th harmonic), corresponding to the low-frequency component (F L (θ, t)) of the armature magnetic motive force.

[0070] ω ep= upω r (22)

[0071] where u is the harmonic number of the permanent magnet MMF (u = 1, 3, 5,...), p is the number of poles of the machine, (ω r ) is the mechanical angular speed of the rotor.

[0072] represents the spatial harmonic angular frequency generated by the permanent magnet rotation, such as the fundamental frequency (u = 1) (pω r )

[0073] The armature reaction MMF is expressed as

[0074]

[0075] The low and high frequency harmonics generated by the armature MMF are denoted as F L (θ, t) and F H (θ, t), respectively, where β l and β h represent their phase angles. Here, l = 6n ± 1, where n = 0, 1, 2,..., ω h corresponds to the electrical angular frequency of f h . The radial and tangential components of F Lr (θ, t) and F Hr (θ, t) are denoted as F Lr_v , F hr_v , F lt_v , and F ht_v , representing the amplitude of the vth harmonic. The F L (θ, t) and F H (θ, t) of the armature MMF can also be analyzed by referring to the process of the no-load air-gap MMF.

[0076] S1.4 Calculation of electromagnetic force density

[0077] The radial and tangential components of the no-load air-gap flux density and the armature air-gap flux density can be obtained from:

[0078] B r_no_load (θ, t) = F 1r (θ, t) Λ r (θ) + F 1t (θ, t) Λ t (θ) (24)

[0079] B t_no_load = F 1t (θ, t) Λ r (θ) - F 1r (θ, t) Λ t (θ) (25)

[0080]

[0081] Therefore, the air-gap flux density on the load is:

[0082] B load (θ, t) = B r_load (θ, t) + jB t_load (θ, t) (28)

[0083] In the formula:

[0084] B r_load (θ, t) = B r_load (θ, t) + B r_arm (θ, t) (29)

[0085] B t_load (θ, t) = B t_load (θ, t) + B t_arm (θ, t)

[0086] The calculation of the electromagnetic force includes the calculation of the radial and tangential electromagnetic force density, which is specifically realized by the following formula:

[0087]

[0088] In the formula, P r (θ, t) and P t (θ, t) are the radial and tangential electromagnetic force densities, respectively, B r (θ, t) and B t (θ, t) are the radial and tangential components of the air-gap flux density, respectively, and μ0 is the vacuum permeability.

[0089]

[0090] In the formula, F r (t) and F t (t) are the radial and tangential components of the concentrated force, respectively, R s is the length and inner radius of the stator core, θ zi is the position of the i-th tooth, and θ st is the tooth surface width angle of the stator.

[0091] Further, the step 2 of the stator natural frequency calculation process based on the energy method is specifically as follows:

[0092] 2.1 Establishment of equivalent model:

[0093] The stator core and the shell of the permanent magnet synchronous motor are equivalent to a cylindrical shell structure with inner and outer ribs, as shown in Figure 12 . Let the average radius of the cylindrical shell be R, the average thickness be h, the axial length be L, and the density be p cs, the elastic modulus is E, and the Poisson's ratio is υ. The axial, circumferential, and radial displacements of any point P on the surface of the cylindrical shell are u, v, and w, respectively. The eccentricity (e arn ) of the equivalent tooth structure is considered, and the displacement expression is:

[0094]

[0095] where, and ψ are the axial and circumferential rotation angles, representing the bending deformation of the shell; r: the radius of any point on the cylindrical shell; w cs : the radial displacement of any point P on the surface of the cylindrical shell; u cs : the axial displacement of any point P on the surface of the cylindrical shell, where u is the axial basic displacement, R is the average radius of the cylindrical shell, is the axial rotation angle, representing the bending deformation of the shell, represents the additional axial displacement due to bending deformation; v cs : the circumferential displacement of any point P on the surface of the cylindrical shell, where v is the circumferential basic displacement, ψ is the circumferential rotation angle, representing the bending deformation of the shell, represents the additional circumferential displacement due to bending deformation;

[0096] 2.2 Potential energy and kinetic energy analysis

[0097] According to the three-dimensional elasticity theory, the potential energy of the cylindrical shell includes membrane strain energy and bending strain energy, and the expression is:

[0098]

[0099] where δ and ε are the normal stress and strain, τ and γ are the shear stress and strain, and after expansion through geometric equations and constitutive relations, they can be separated into membrane force and bending moment contributions.

[0100] The kinetic energy considers the inertial effect of the shell, and the expression is:

[0101]

[0102] After neglecting high-order small quantities, it is simplified to the contribution of mid-surface motion kinetic energy and moment of inertia, and the mid-surface refers to an ideal cylindrical surface with an average radius R, which is located at the middle position of the shell thickness, i.e., the distance from the inner and outer surfaces of the shell is half of the shell thickness h:

[0103]

[0104] 2.3 Variable separation and mode assumption

[0105] Assuming that the stator vibration is a simple harmonic mode, the variable separation method is adopted, and the displacement form is:

[0106]

[0107] where n is the circumferential mode order, ω is the natural angular frequency, U(z), V(z), W(z) are the axial mode shape functions.

[0108] Considering the axial symmetry of the stator, the axial mode shape is assumed to be a sine or cosine function, such as (U(z) = A cos(kz), (k = mπ / L) (m is the axial mode order).

[0109] 2.4 Derivation of frequency equation

[0110] Substitute the displacement assumption into the potential energy and kinetic energy expressions, and use Hamilton's principle (δ∫(K E -P E )dt = 0) to obtain a homogeneous equation system for amplitudes A, (V n ), (W n ):

[0111]

[0112] where (K ij ) and (M ij ) are the stiffness and mass matrix elements, respectively, obtained by integrating the strain energy and inertial terms in the potential and kinetic energies, V n : amplitude of the axial mode shape function V(z) corresponding to the circumferential displacement v at the circumferential mode order n; W n : amplitude of the axial mode shape function W(z) corresponding to the radial displacement w at the circumferential mode order n;

[0113] To have a non-zero solution for the equation system, the determinant of the coefficient matrix must be zero, i.e.:

[0114] det(K - ω 2 M) = 0

[0115] After expansion, the modal frequency equation is obtained: ω 4 - ξ1ω 2 + ξ2 = 0

[0116] D = Eh 3 / [12(1 - v 2 )] is the bending stiffness of the cylindrical shell, (K = Eh) is the membrane stiffness.

[0117] 2.5 Solution of natural frequencies

[0118] Solving the frequency equation gives: For low-order circumferential modes (n = 0, 1, 2,...), ignoring axial high-order vibrations (m = 0), the simplified equation gives the natural frequencies of radial and circumferential coupled vibrations:

[0119]

[0120] where, is the bending stiffness, K = Eh is the membrane stiffness, embodying the elastic property of the stator structure; p cs is the density, h is the thickness, R is the average radius, reflecting the mass distribution; n is the circumferential modal order, determining the vibration mode (e.g. (n = 0) is rigid body vibration, (n = 1) is first order bending vibration).

[0121] Further, the step 3 is specifically as follows:

[0122] For the mass M n , stiffness K n , damping coefficient C n of the mechanical system, under the action of the periodic alternating force F m (t), the electromagnetic vibration process can be expressed as:

[0123]

[0124] The solution of the system vibration acceleration is:

[0125]

[0126] In the formula:

[0127]

[0128] ξ n is the n-order circumferential modal corresponding damping ratio, which can be obtained by an empirical formula:

[0129]

[0130] F m (t) = 2πR i lp m (45)

[0131] In the formula, R i and l are the inner radius and axial length of the equivalent ring of the stator respectively; PM is the amplitude of the electromagnetic force density.

[0132] Further, the step 2) further comprises:

[0133] 2.6 equivalent tooth structure correction

[0134] 2.6.1 potential energy and kinetic energy analysis of equivalent tooth structure:

[0135] The potential energy and kinetic energy of the equivalent teeth attached to the inner surface of the cylindrical shell need to be added to the overall model. Let the equivalent tooth elastic modulus be (E ar ), the cross-sectional area be (A ar), the inertia moments around x, y axes are I xar , (I yar ), then the additional potential energy is:

[0136]

[0137] where N ar : represents the number of equivalent teeth (or the total number of teeth), used for summing the additional potential energy of all equivalent teeth; u arn : represents the displacement component of the nth equivalent tooth in the axial (z direction), reflecting the deformation of the tooth structure in the axial direction; w arn : represents the displacement component of the nth equivalent tooth in the radial (r direction), reflecting the deformation of the tooth structure in the radial direction; is the partial derivative symbol with respect to the axial coordinate z, representing the rate of change of the corresponding displacement component along the axial direction;

[0138] Additional kinetic energy:

[0139]

[0140] P ar is the equivalent tooth density; are the axial and radial velocities, respectively;

[0141] 2.6.2 Modification of stiffness matrix K:

[0142] The additional potential energy corresponds to the increment of the stiffness matrix, which needs to consider the axial stiffness and bending stiffness:

[0143] Axial stiffness modification:

[0144]

[0145] where φ n (z) is the axial mode function (such as (cos(kz)) or (sin(kz)), k = m\π / L.

[0146] Bending stiffness modification

[0147]

[0148] ψ n (z): radial mode function (similar to axial mode, take low order mode).

[0149] Total stiffness matrix modification

[0150]

[0151] 2.6.3 Modification of mass matrix M:

[0152] The additional kinetic energy corresponds to the increment of the mass matrix ΔM, considering the mass distribution of the equivalent teeth:

[0153]

[0154] Mode function n (z) is consistent with the axial mode, reflecting the distribution of mass along the axial direction.

[0155] Total mass matrix correction

[0156] M 修正 = M 原 + ΔM ar

[0157] 2.6.4 Frequency equation after considering the tooth structure:

[0158] Substitute the corrected stiffness matrix and mass matrix into the eigenvalue equation to solve the natural frequency:

[0159] det(K 修正 - ω 2 M 修正 )

[0160] After expansion, the modified modal frequency equation is obtained, and the natural frequency f n considering the influence of the tooth structure is obtained by solving the equation:

[0161]

[0162] By modifying the stiffness matrix K and the mass matrix M, the influence of the tooth structure on the natural frequency is considered.

[0163] The design idea of the application is as follows:

[0164] Electromagnetic force calculation steps:

[0165] Determine the basic parameters of the motor, such as the number of poles, the number of slots, the diameter of the stator and rotor, the thickness of the permanent magnet, etc. These parameters are the basis for subsequent calculations.

[0166] Using the magnetic potential permeability method and the improved subdomain method, the above air gap magnetic motive force, air gap permeability, armature magnetic motive force, air gap flux density and electromagnetic force density calculation method is used to calculate in turn. In the calculation process, the parameter values are reasonably determined according to the actual situation, and accurate calculation is carried out according to the corresponding formula.

[0167] Stator natural frequency calculation steps

[0168] According to the structural characteristics of the motor stator core and the shell, determine the relevant parameters of the equivalent cylindrical shell with inner and outer ribs, such as average radius, average thickness, average axial length, etc.

[0169] According to the three-dimensional elasticity theory, accurately substitute the parameters into the potential energy and kinetic energy expressions of the cylindrical shell, equivalent teeth and stator yoke.

[0170] The stator natural frequency is solved by using energy method, and is verified by modal hammer test.

[0171] Electromagnetic vibration calculation step

[0172] According to the structure and material characteristics of the motor, the parameters such as mass Mn, stiffness Kn and damping coefficient Cn of the mechanical system are determined.

[0173] The calculated electromagnetic force component is substituted into the vibration equation to calculate the acceleration response function spectrum of the low-order mode of the stator with shell.

[0174] Through the above detailed implementation steps, the method of the application can accurately calculate the electromagnetic vibration of the stator tooth surface of the permanent magnet synchronous motor, provide reliable data support for the design optimization of the motor, and effectively improve the performance and quality of the motor.

[0175] The beneficial effects of the application are as follows:

[0176] Precision improvement: the method proposed by the application effectively solves the complexity problem of electromagnetic field distribution in nonlinear materials, accurately simulates the influence of tangential electromagnetic force wave and harmonic current, and greatly improves the accuracy of the calculation results.

[0177] Efficient calculation: based on energy method, the stator natural frequency is quickly calculated by establishing reasonable equivalent theoretical models of stator core and shell, realizing efficient and accurate acquisition of anisotropic stator core and isotropic shell natural frequency.

[0178] Verification effectiveness: compared with traditional electromagnetic-structure coupled finite element analysis, the method of the application not only effectively verifies its accuracy, but also reduces the total calculation time by 59.3% when calculating electromagnetic vibration related physical indicators, showing high efficiency in engineering application and providing strong technical support for design and optimization of permanent magnet synchronous motor. BRIEF DESCRIPTION OF DRAWINGS

[0179] Figure 1 Table-mounted permanent magnet synchronous motor slotless structure schematic diagram;

[0180] Figure 2 Table-mounted PMSM slot structure diagram;

[0181] Figure 3 Armature magnetic motive force calculation schematic diagram;

[0182] Figure 4 Permanent magnet synchronous motor electromagnetic model schematic diagram;

[0183] Figure 5 Permanent magnet synchronous motor no-load air gap magnetic motive force schematic diagram;

[0184] Figure 6 Permanent magnet synchronous motor air gap permeability schematic diagram;

[0185] Figure 7 Harmonic current schematic diagram for phase a;

[0186] Figure 8 Air gap flux density schematic diagram for 6-pole 36-slot motor under rated load; wherein (a) is a radial component diagram; (b) is a tangential component diagram;

[0187] Figure 9 Time-space distribution diagram of electromagnetic force; wherein (a) is a radial electromagnetic force density diagram; (b) is a tangential electromagnetic force density diagram;

[0188] Figure 10 Air gap electromagnetic force magnetic density two-dimensional Fourier decomposition comparison diagram; wherein (a) is a radial electromagnetic force density diagram under rated load, (b) is a tangential electromagnetic force density diagram;

[0189] Figure 11 Stator three-dimensional finite element model and equivalent model schematic diagram;

[0190] Figure 12 Ribbed cylindrical shell structure schematic diagram;

[0191] Figure 13 Stator object and model schematic diagram;

[0192] Figure 14 Modal test platform schematic diagram based on hammering method;

[0193] Figure 15 Data acquisition coherence curve diagram;

[0194] Figure 16 Analytical method, finite element and experimental modal comparison diagram of the stator;

[0195] Figure 17 Analytical method, finite element modal comparison diagram of the stator with shell;

[0196] Figure 18 Stator electromagnetic force schematic diagram; wherein (a) is a stator induced electromagnetic force diagram; (b) is a stator tooth surface electromagnetic force equivalent diagram;

[0197] Figure 19 Low-order modal frequency response function diagram;

[0198] Figure 20 Acceleration spectrum of the shell stator prototype. DETAILED DESCRIPTION

[0199] The technical solutions in the embodiments of the present application will be clearly and completely described below with reference to the drawings in the embodiments of the present application. Obviously, the described embodiments are only part of the embodiments of the present application, rather than all the embodiments of the present application.

[0200] The present application combines the magnetic potential magnetic permeability method with the improved sub-domain method, and proposes a new method for calculating electromagnetic force. The magnetic potential magnetic permeability method effectively solves the complexity of electromagnetic field distribution, especially when dealing with nonlinear materials. The improved sub-domain method divides the calculation domain into multiple sub-domains, improving the accuracy of capturing electromagnetic field changes. By combining the two methods, the method considers the influence of tangential electromagnetic force wave and harmonic current, significantly improving the accuracy of the results. On this basis, a reasonable equivalent theoretical model of the stator core shell is proposed, and the natural frequency of the stator is calculated by using the energy method. The accuracy of the model is verified through modal experiments. In addition, the frequency response function of the stator under low-order modes is determined, and the acceleration spectrum and electromagnetic force wave of the prototype stator under each mode are obtained based on the linear superposition principle. Compared with the established electromagnetic-structure coupled finite element model, the calculation time of this method is shortened by 59.3%.

[0201] Specifically includes the following steps:

[0202] S1. Calculation of electromagnetic force

[0203] According to the Maxwell stress tensor method, the electromagnetic force density in the motor air gap is:

[0204]

[0205]

[0206] In the formula, P r (θ,t) and P t (θ,t) are the radial and tangential electromagnetic force densities, respectively, B r (θ,t) and B t (θ,t) are the radial and tangential components of the air gap flux density, respectively, and μ0 is the vacuum permeability.

[0207]

[0208] In the formula, F r (t) and F t (t) are the radial and tangential components of the concentrated force, respectively, R s is the inner radius of the stator, θ zi is the position of the i th tooth, and θ stis the stator tooth face width angle.

[0209] S1.1 Calculation of air gap magnetomotive force under no-load conditions

[0210] The slotless structure of permanent magnet synchronous motor is as follows Figure 1 As shown. The slotless motor model can be divided into two subdomains: PM subdomain I and gap subdomain II. F1(θ) and B o (θ) are the distributions of no-load air gap magnetomotive force and air gap magnetic flux density, δ is the air gap length, h pm is the PM thickness. F1(θ) can be calculated as:

[0211]

[0212] F1(θ)=F 1r (θ)+jF 1t (θ) (6)

[0213] B o (θ)=B or (θ)+B ot (θ) (7)

[0214] Among them F 1r (θ), F 1t (θ), B or (θ) and B ot (θ) is F1(θ) and B o The radial and tangential components of (θ). J is an imaginary operator. F 1r (θ) and F 1t (θ) have the same harmonic content, but their amplitude and phase are different. The air gap flux density B can be obtained quickly and accurately using the subdomain method. o (θ). In addition, the magnetomotive force F1(θ,t) can be expressed as:

[0215]

[0216] Where ω r is the angular velocity, p is the number of poles, F 1r_μ and F 1t_μ is the radial component F1(θ,t) 1r (θ, t) and the tangential component F 1t The μth harmonic amplitude of (θ, t) is calculated by F 1r (θ) and F 1t The Fourier decomposition of (θ) is obtained.

[0217] S1.2 Calculation of air gap permeability

[0218] To simplify the analytical calculation model, the following assumptions are made:

[0219] 1) The relative permeability of the fixed and rotor ferromagnetic regions is infinite;

[0220] 2) The permanent magnets have linearized characteristics with a relative permeability ur = 1 ;

[0221] 3) The axial end effects of the machine are neglected;

[0222] 4) The stator slots of the machine are equivalent to slotted.

[0223] The slot structure of a permanent magnet synchronous machine is shown in Figure 2 where Q1 is the number of slots and (-π / Q1, π / Q1) is the pitch of the teeth within the arc range of the permanent magnet. The slotted machine model can be divided into four sub-domains: sub-domain I is the permanent magnet (PM), sub-domain II is the air gap, sub-domain III is the slot opening, and sub-domain IV is the slot body. The intrinsic characteristics of the machine become particularly evident in the stator teeth. By applying surface currents on both sides of the stator teeth, the calculation can be performed using an iterative method. The magnetic flux generated by the surface currents is opposite to the main flux, effectively simulating the saturation effect within the teeth. It is generally considered that the saturation on both sides of the slot opening is in a vacuum state. Considering the surface currents on both sides of the opening, the vector potential A z3i in sub-domain III can be expressed as.

[0224]

[0225] where R t is the radius of the stator slot top, b oa is the slot width angle, Q 3i and D 3i are the coefficients to be determined, the amplitude of the constant component j 3io and the nth harmonic component J 3in of the current density on both sides of the i-th slot can be expressed as:

[0226]

[0227] where S c is the area of each surface current.

[0228] The vector potential A z4i in the slot sub-domain IV can be directly derived. By solving the coupled equations defined by the boundary conditions of each sub-domain, the values of Q 3i , D 3i and the undetermined coefficients of the other sub-domains can be determined. By analyzing the distribution of the air gap flux density B δt (θ) along the circumference within the specified tooth pitch range, the air gap permeability Λ(θ) within the tooth pitch can be calculated.

[0229]

[0230]

[0231] where B δt_r (θ) and B δt_t (θ) are the radial and tangential components of Bδt(θ) and can be obtained by the subdomain method. r (θ),Λ t (θ) are the radial and tangential components of Λ(θ) and can be obtained from equations (15), (16), (17).

[0232]

[0233] The air-gap permeability Λ(θ) can also be expressed as:

[0234]

[0235] where Λ r_0 is the constant component of Λ r (θ), Λ r_n and Λ t_n are the n-th harmonic amplitudes of Λ r (θ) and Λ t (θ), respectively. These amplitudes are obtained from the Fourier decomposition of Λ r (θ) and Λ t (θ).

[0236] S1.3 Calculation of the armature reaction magnetomotive force

[0237] The permeability is nonlinear when magnetic saturation is considered. From equation (1), the flux generated by the permanent magnet is nonlinear with respect to the current. To simplify the analysis in this section, we consider the saturation effect of the magnetic saturation with respect to the fundamental permanent flux (ψ(i d ,i q )). Thus, the relationship between the dq-axis flux and the current is given by equations (18) and (19):

[0238]

[0239]

[0240] From the expression, ψ 2i-1 represents the open-circuit flux linkage of the (2i-1)-th harmonic component. Here, R is the phase resistance, L d and L q are the phase inductances, which are equal in surface-mounted PMSMs (L d = L q ). The parameters ω r , J, T load and D fcorrespond to the mechanical angular velocity, the moment of inertia, the load torque and the damping coefficient, respectively. By establishing the analytical model of the surface-mounted permanent magnet synchronous motor vector control system based on MATLAB space vector pulse width modulation (SVPWM), the inverter output voltage u d and u q are derived. The fourth-order Runge-Kutta method is used to solve the coupled equations of flux linkage, voltage, electromagnetic torque and motion dynamics, so as to determine the phase current i d and i q .

[0241] The calculation model of the armature magnetic motive force is shown in Figure 3 , which only includes the air gap sub-domain II. The center line of each slot corresponds to a surface current Q1 uniformly distributed along the inner surface of the stator core, which is represented as a line current in a two-dimensional field. The magnetic field intensity H θ satisfies specific boundary conditions at the tangential components of the inner and outer boundaries of the air gap sub-domain II.

[0242] The armature reaction magnetic motive force is represented as

[0243]

[0244] The low-frequency and high-frequency harmonics generated by the armature magnetic motive force are denoted as F L (θ, t) and F H (θ, t), respectively, where β l and β h represent their phase angles. Here, l = 6n ± 1, where n = 0, 1, 2, …, and ω h corresponds to the electrical angular frequency of f h . The radial components F Lr (θ, t) and F Hr (θ, t) are denoted as F Lr_v , F hr_v , F lt_v , and F ht_v , representing the amplitude of the vth harmonic. By referring to the process of the no-load air gap magnetic motive force, F L (θ, t) and F H (θ, t) of the armature magnetic motive force can also be analyzed.

[0245] S1.4 Calculation of electromagnetic force density

[0246] The radial and tangential components of the no-load air gap flux density and the armature air gap flux density can be obtained from:

[0247] B r_no_load (θ, t) = F 1r (θ, t) Λ r (θ) + F 1t (θ, t) Λ t (θ) (23)

[0248] B t_no_load = F 1t (θ, t) Λ r (θ) - F 1r (θ, t) Λ t (θ) (24)

[0249]

[0250] Therefore, the air-gap flux density on the load is:

[0251] B load (θ, t) = B r_load (θ, t) + jB t_load (θ, t) (27)

[0252] In the formula:

[0253] B r_load (θ, t) = B r_load (θ, t) + B r_arm (θ, t) (28)

[0254] B t_load (θ, t) = B t_load (θ, t) + B t_arm (θ, t) (29)

[0255] Electromagnetic force density calculation

[0256] Ignoring the interaction between the magnetic fields generated by the harmonic currents, the main components of the electromagnetic force density of the inverter motor are ranked as shown in Table 2.

[0257] Table 1 Motor main parameters

[0258]

[0259] The air-gap flux density and electromagnetic force density of a 6-pole 36-slot surface-mounted permanent magnet synchronous motor are calculated. The motor model and main structural parameters are shown in Figure 4 The no-load air-gap magnetomotive force and air-gap permeability of the 6-pole 36-slot surface-mounted PMSM calculated by the amplitude modulation method and the finite element method are shown in Figure 5 and Figure 6 The a-phase harmonic current and its harmonic components when the motor is running stably at rated load when the inverter switching frequency is 5 kHz are shown in Figure 7 The fundamental frequency amplitude and low harmonic amplitude of the prototype test current are slightly larger than the amplitude modulation calculation value. The radial and tangential components of the air-gap flux density under rated load are further calculated, as shown in Figure 8 Obviously, considering the influence of saturation, the calculation result is more accurate.

[0260] S2. Calculation of stator natural frequency

[0261] Table II. Main low order error analysis of the radial electromagnetic force density at no-load

[0262]

[0263] Table III. Main low order error analysis of the tangential electromagnetic force density at no-load

[0264]

[0265] The spatial and temporal distribution of the radial and tangential electromagnetic force densities in a permanent magnet synchronous machine is shown in Figure 9 . The frequency characteristics and spatial harmonics obtained by two-dimensional Fourier decomposition of the finite element analysis (FEA) and analytical method are compared as shown in Figure 10 . This comparison verifies the accuracy of the analysis and calculation in this paper. The stator core and frame can be modeled as a cylindrical shell with inner and outer ribs as shown in Figure 11 . As shown in Figure 12 , for any point P on the surface of the ribbed cylindrical shell, let u, v, w represent the displacement of point P in the axial (z), circumferential (θ), and radial (r) directions, respectively. The parameters r, h, L represent the average radius, average thickness, and average axial length, respectively, and e arn represents the eccentricity of the i-th frame. The displacement of any point on the cylinder can be expressed as shown in equation (30).

[0266]

[0267] In the equation and ψ / R represent the axial and circumferential directions.

[0268] r: radius of any point on the cylindrical shell.

[0269] w cs : radial displacement of any point P on the surface of the cylindrical shell.

[0270] u cs : axial displacement of any point P on the surface of the cylindrical shell, where u is the axial fundamental displacement, R is the average radius of the cylindrical shell, is the axial rotation angle, representing the bending deformation of the shell, represents the additional axial displacement due to bending deformation.

[0271] v cs : circumferential displacement of any point P on the surface of the cylindrical shell, where v is the circumferential fundamental displacement, (ψ) is the circumferential rotation angle, representing the bending deformation of the shell, represents the additional circumferential displacement due to bending deformation.

[0272] 2.2 Potential and kinetic energy analysis

[0273] According to the three-dimensional elastic theory, the potential energy of the cylindrical shell contains membrane strain energy and bending strain energy, expressed as:

[0274]

[0275] where δ and ε are normal stress and strain, τ and γ are shear stress and strain, which can be separated into membrane force and bending moment contributions after expanding through geometric equations and constitutive relations;

[0276] The kinetic energy considers the inertia effect of the shell, expressed as:

[0277]

[0278] Neglecting high-order small quantities, it is simplified to the contribution of mid-surface motion kinetic energy and rotational inertia:

[0279]

[0280] 2.3 Variable separation and mode assumption:

[0281] Assuming that the stator vibration is a simple harmonic mode, the variable separation method is used, and the displacement form is set as:

[0282]

[0283] where n is the circumferential mode order, ω is the natural angular frequency, U(z), V(z), and W(z) are the axial mode functions;

[0284] Considering the axial symmetry of the stator, it is assumed that the axial mode is a sine or cosine function;

[0285] 2.4 Frequency equation derivation:

[0286] Substitute the displacement assumption into the potential energy and kinetic energy expressions, and use the Hamilton principle to obtain the homogeneous equation group about the amplitudes A, V n , and W n :

[0287] V n : Amplitude of the axial mode function V(z) corresponding to the circumferential displacement v at the circumferential mode order n.

[0288] W n : Amplitude of the axial mode function W(z) corresponding to the radial displacement w at the circumferential mode order n.

[0289]

[0290] where K ij and M ijare the stiffness and mass matrix elements, respectively, obtained by integrating the strain energy and kinetic energy in potential and kinetic terms;

[0291] For non-zero solution of the equation set, the determinant of the coefficient matrix needs to be zero, i.e.:

[0292] det(K-ω 2 M)=0

[0293] After expansion, the modal frequency equation is obtained: 4 -ξ1ω 2 +ξ2=0

[0294] D=Eh 3 / [12(1-ν 2 )]

[0295] K is the bending stiffness of the cylindrical shell, and K = Eh is the membrane stiffness;

[0296] 2.5 Solution of natural frequencies:

[0297] The frequency equation is solved to obtain: For low-order circumferential modes (n = 0, 1, 2, …), the axial high-order vibration (m = 0) is ignored, and after simplification, the natural frequencies of the radial and circumferential coupled vibration are obtained:

[0298]

[0299] 2.6 Modification of equivalent tooth structure:

[0300] 2.6.1 Potential energy and kinetic energy analysis of equivalent tooth structure:

[0301] The stator teeth are considered to be equivalent to elastic beams attached to the inner surface of the cylindrical shell, and their potential energy and kinetic energy need to be added to the overall model. Let the equivalent tooth elastic modulus be E ar , the cross-sectional area be A ar , and the moments of inertia around the x and y axes be I xar and I yar , respectively. Then the additional potential energy is:

[0302]

[0303] where N ar : represents the number of equivalent teeth (or the total number of teeth), and is used for summation to calculate the additional potential energy of all equivalent teeth; u arn : represents the displacement component of the nth equivalent tooth in the axial (z-direction), reflecting the deformation of the tooth structure in the axial direction; w arn : represents the displacement component of the nth equivalent tooth in the radial (r-direction), reflecting the deformation of the tooth structure in the radial direction; is the partial derivative symbol with respect to the axial coordinate z, indicating the rate of change of the corresponding displacement component along the axial direction; E ar is the equivalent tooth elastic modulus; A ar is the equivalent tooth cross-sectional area; I xar is the equivalent tooth moment of inertia about the x-axis; L is the stator axial length;

[0304] Additional kinetic energy:

[0305]

[0306] P ar is the equivalent tooth density; are the axial and radial velocities, respectively;

[0307] 2.6.2 Modification of the stiffness matrix K:

[0308] The increment of the stiffness matrix corresponding to the additional potential energy, which needs to consider both the axial stiffness and the bending stiffness:

[0309] Axial stiffness modification:

[0310]

[0311] where φ n (z) is the axial mode shape function;

[0312] Bending stiffness modification

[0313]

[0314] ψ n (z): radial mode shape function;

[0315] Total stiffness matrix modification

[0316]

[0317] 2.6.3 Modification of the mass matrix M:

[0318] The increment of the mass matrix ΔM corresponding to the additional kinetic energy, considering the mass distribution of the equivalent tooth:

[0319]

[0320] The mode shape function φ n (z) is consistent with the axial mode shape, reflecting the mass distribution along the axial direction;

[0321] Total mass matrix modification

[0322] M 修正 = M 原 + ΔM ar (46)

[0323] 2.6.4 Frequency equation after considering the tooth structure:

[0324] Substitute the modified stiffness matrix and mass matrix into the eigenvalue equation to solve the natural frequency:

[0325] det(K 修正 -ω 2 M 修正 )(47)

[0326] After expansion, the modified modal frequency equation is obtained, and the natural frequency f n considering the influence of the tooth structure is obtained by solving the equation:

[0327]

[0328] To verify the accuracy of the energy method calculation results, a modal hammering test was conducted on the stator of a 6-pole 36-slot permanent magnet synchronous motor. The stator and test device are shown in Figure 13 , and the modal test arrangement is shown in Figure 14 . The stator is suspended by nylon rope, and the response signal generated by hammering is transmitted to the signal acquisition and analysis instrument through the accelerometer. During data acquisition, the coherence curve is shown in Figure 15 . The coherence function can be calculated by the spectrum of the excitation force and the response acceleration. Coherence refers to the consistency of the frequency response function between multiple averages, and in engineering practice, if the value is greater than 0.75, it is generally considered to be good. From Figure 15 , it can be seen that the coherence of this modal test is relatively high. From Figure 16 , it can be seen that the equivalent method used in this study can accurately determine the natural frequency corresponding to the low-order circumferential mode of the stator, with an error of no more than 10%.

[0329] Table 4 Anisotropic material parameters of stator core and casing

[0330]

[0331] Where:

[0332] 1. Anisotropic elastic modulus (E ar )

[0333] The stator core exhibits different elastic moduli in different directions:

[0334] Radial (x, y direction) elastic modulus: (E ar (x) = E ar (y) = 198 GPa)

[0335] Axial (z direction) elastic modulus: (E ar (z) = 155 GPa)

[0336] 2. Anisotropic shear modulus (G ar )

[0337] The shear modulus varies with the direction of the shear plane:

[0338] Shear modulus along the yz, xz plane: (G ar (yz) = G ar (xz) = 73.2 GPa)

[0339] Shear modulus along the xy plane: (G ar (xy) = 81.6 GPa)

[0340] 3. Anisotropic Poisson's ratio (v ar )

[0341] Poisson's ratio reflects the ratio of strain in different directions:

[0342] Poisson's ratio along the yz, xz plane: (v ar (yz) = v ar (xz) = 5.4)

[0343] Poisson's ratio along the xy plane: (v ar (xy) = 4.6)

[0344] The material properties of the stator core exhibit significant orthotropic anisotropy, showing significant anisotropy in the axial and transverse directions. The anisotropic material parameters of each part of the stator are shown in Table 4.

[0345] S3. Electromagnetic vibration calculation

[0346] For a mechanical system with mass M n , stiffness K n , and damping coefficient C n , under the action of a periodic alternating force with amplitude F m (t), its vibration process can be represented as:

[0347]

[0348] Then the solution of the system vibration acceleration is:

[0349]

[0350] In the formula:

[0351]

[0352] ξ n is the damping ratio corresponding to the n-th circumferential mode, which can be obtained from the empirical formula:

[0353]

[0354] F m (t)=2πR i lp m (53)

[0355] Where R i and l are the inner radius and axial length of the stator equivalent ring respectively; PM is the amplitude of the electromagnetic force density.

[0356] Substituting the electromagnetic force component obtained each time into formula (38), the low-order modal acceleration response function spectrum of the shell stator can be obtained, as follows: Figure 19 As shown. Assuming n = 1, the acceleration spectrum of the prototype stator with shell under the action of low-order electromagnetic force waves can be obtained. From the results, it can be seen that the vibration of the motor is similar to the vibration of a multi-degree-of-freedom rigid body, and the modes of each order are basically decoupled. Therefore, the total acceleration response spectrum can be obtained by superimposing the acceleration spectrum of the prototype stator and each mode. Assume that the radial electromagnetic force applied to the stator core acts entirely on the center line of the stator teeth and is evenly distributed along the axial direction, as shown in Figure 18 As shown, the electromagnetic force acting on each stator tooth can be calculated by multiplying the electromagnetic force density at the tooth center by the surface area of ​​the stator tooth. Figure 20 The vibration spectrum of a stator with a housing is obtained using the finite element method after applying the electromagnetic force density at the center of the stator teeth to the inner surface. Although this method exhibits some discrepancies compared to the results obtained using the finite element method, it still accurately reflects the stator's vibration characteristics and the distribution of the vibration acceleration spectrum. Therefore, the linear superposition method of vibration responses can quickly predict the motor's vibration response, significantly improving the efficiency of motor vibration prediction.

[0357] The frequency response function of the stator low-order mode is as follows: Figure 17 The results show that the vibration of the stator has a large decoupling between different modes, and the contribution of the zero-order mode of the stator is significant and cannot be ignored.

[0358] Table 5 shows the calculation time for various physical parameters related to electromagnetic vibration in the prototype. The finite element software was run using a 3.0 GHz CPU and 64 GB of memory. Compared to the multiphysics coupled finite element method, the proposed additive method reduces the total calculation time by over 59.3%. In particular, the superiority of the proposed analytical method further highlights its ability to avoid tedious modeling work during the motor optimization design phase.

[0359] Table 5 Calculation time of various physical indicators

[0360]

Claims

1. A fast calculation method for electromagnetic vibration of a permanent magnet synchronous motor, characterized in that: The following steps are involved: 1) Calculate the electromagnetic force by combining the magnetic potential permeance method with the improved subdomain method; 2) Establish an equivalent theoretical model of the stator core and housing, and use the energy method to calculate the natural frequency of the stator; 3) Introduce electromagnetic force into the structural model in the form of nodal force and use the modal superposition method to calculate electromagnetic vibration; 4) Verify the calculation results of the stator natural frequency to ensure the accuracy of the calculation results.

2. The method for rapid calculation of electromagnetic vibration of a permanent magnet synchronous motor according to claim 1, characterized in that: The specific process of step 1 is as follows: S1.1 Calculation of permanent magnet air gap magnetomotive force: The slotless motor model is divided into two subdomains: PM subdomain I and gap subdomain II; F1(θ) and B o (θ) are the distributions of no-load air gap magnetomotive force and air gap magnetic flux density, δ is the air gap length, h pm is the PM thickness; F1(θ) is calculated as: F1(θ)=F 1r (θ)+jF 1t (i) (2) B o (θ)=B or (i)+B ot (i) (3) Among them F 1r (θ), F 1t (θ), B or (θ) and B ot (θ) is F1(θ) and B o The radial and tangential components of (θ), j is an imaginary operator; F 1r (θ) and F 1t (θ) have the same harmonic content, but their amplitudes and phases are different, μ0 is the magnetic permeability of vacuum; The permanent magnet air gap magnetomotive force F1(θ,t) is expressed as: Where ω r is the angular velocity, p is the number of poles, F 1r_μ and F 1t_μ is the radial component F1(θ,t) 1r (θ, t) and the tangential component F 1t The μth harmonic amplitude of (θ, t) is calculated by F 1r (θ) and F 1t The Fourier decomposition of (θ) yields; S1.2 Calculation of air gap permeability: The basis for dividing subdomains: a) The relative permeability of the fixed and rotor ferromagnetic regions is infinite; b) Permanent magnets have linear characteristics, and their relative permeability μ r =1; c) Ignore the motor axial end effect; d) The stator slots of the motor are equivalent to slots; According to the above division basis, the slotted motor model is divided into four subdomains: subdomain I is the permanent magnet PM, subdomain II is the air gap, subdomain III is the slot, and subdomain IV is the slot body; by applying surface current on both sides of the stator teeth, the saturation effect in the teeth is simulated, and the saturation on both sides of the slot is considered to be in a vacuum state. Considering the surface current on both sides of the opening, the vector potential A in subdomain III is z3i Expressed as: Where A z3i is the vector potential of the ith slot in subdomain III; A 3i0 A is the constant component of the vector potential, reflecting the DC or average effect of the magnetic field; 3im is the mth spatial harmonic amplitude of the vector potential, representing the periodic fluctuation of the magnetic field; F n is the spatial harmonic order, and F is defined n =nπ / b oa , which is related to the slot width angle; θ is the circumferential coordinate; θ zi is the center position angle of the i-th slot, used to locate the harmonic spatial phase; R s is the inner radius of the stator; R t is the stator slot top radius, b oa is the slot width angle, Q 3i and D 3i is the coefficient to be determined; The current density constant component j on both sides of the i-th slot 3io and the nth harmonic component J 3im The amplitude is expressed as: Among them S c is the area per surface current, I 1i =I 2i =I i / 2, where I i is the phase current of the i-th slot; By analyzing the air gap magnetic flux density along the circumference B within the specified tooth pitch range δt (θ) distribution, calculate the air gap permeability Λ(θ) within the tooth spacing; Among them, B δt_r (θ) and B δt_t (θ) are the radial and tangential components of Bδt(θ); Λ r (θ),Λ t (θ) is the radial and tangential components of Λ(θ), which are obtained from the following equations: The air gap permeability Λ(θ) is expressed as: Where, Λ r_0 Λ r The constant component of (θ), Λ r_n and Λ t_n are Λ r (θ) and Λ t The amplitude of the nth harmonic of (θ); S1.3 Calculation of armature reaction magnetomotive force: When considering magnetic saturation, the magnetic permeability changes nonlinearly, and the magnetic flux generated by the permanent magnet is nonlinear with respect to the current. In order to simplify the analysis, the saturation effect of magnetic saturation relative to the basic permanent magnetic flux (ψ(i d ,i q )); therefore, the relationship between the dq axial flux and the current is as follows: Where: d is the d-axis flux, which is composed of the flux generated by the d-axis current id and the inductance Ld and the superposition of the harmonic components; ψ q is the q-axis flux, which is composed of the flux generated by the q-axis current iq and the inductance Lq and the superposition of the harmonic components; L d ,L q are the d-axis and q-axis inductances, characterizing the linear relationship between current and flux; ψ1 is the fundamental flux component, generated by the permanent magnet or main magnetic field; ψ3, ψ5, ψ7, ψ9 are the 3rd, 5th, 7th, and 9th harmonic flux components, caused by the spatial harmonics of the magnetic field; cos6θ, cos12θ, sin6θ, sin12θ are spatial harmonic functions, θ is the mechanical angle, and the coefficients 6 and 12 correspond to the combined effect of the number of pole pairs and the number of slots; ψ dk is the kth harmonic no-load flux linkage component of the d-axis, indicating the kth harmonic component in the d-axis direction of the permanent magnet magnetic field; ψ qk is the kth harmonic no-load flux linkage component of the q-axis, indicating the kth harmonic component in the q-axis direction of the permanent magnet magnetic field, which works together with the d-axis harmonic component to produce torque pulsation; T em Electromagnetic torque is the core physical quantity of motor energy conversion, which represents the torque generated by the interaction between magnetic field and current; No-load flux fundamental component ψ i-1 The transpose or equivalent expression of , which represents the coupling effect between the fundamental magnetic flux generated by the permanent magnet and the current; p: number of motor pole pairs; ψ 2i-1 represents the no-load flux linkage of the harmonic component (2i-1); R is the phase resistance; parameter ω r , J, T load and D f Corresponding to mechanical angular velocity, moment of inertia, load torque and damping coefficient respectively; The inverter output voltage u is derived by the analytical model of the surface mounted permanent magnet synchronous motor vector control system based on MATLAB space vector pulse width modulation SVPWM established above. d and u q The fourth-order Runge-Kutta method is used to solve the flux, voltage, electromagnetic torque and motion dynamics coupling equations to determine the phase current i d and i q ; The calculation model of armature magnetomotive force only includes the air gap subdomain II; the center line of each slot corresponds to a surface current J uniformly distributed along the inner surface of the stator core. i , expressed as a line current in a two-dimensional field; the magnetic field intensity H θ The tangential components at the inner and outer boundaries of the air gap subdomain II meet specific boundary conditions; Where, I i is the phase current of the i-th slot, S c is the current distribution area on the surface of a single slot; when the current on both sides of the slot I 1i =I 2i =I i / 2, the constant component is half of the phase current divided by the slot area; Where n is the harmonic order, sin(nπ / 2)) is non-zero when reflecting odd harmonic characteristics (n=1,3,5,…), and the amplitude of the harmonic component is related to the phase current (I i ) is proportional to the harmonic order n and inversely proportional to the harmonic order n; oh el =lω e =(6n+1)ω e (n=0,1,2,…) (21) Where ω e is the fundamental electrical angular frequency, corresponding to the synchronous speed of the motor; (l=6n±1) represents the low harmonic order, corresponding to the low-frequency component of the armature magnetomotive force (F L (θ,t)); oh ep =upω r (22) Where u is the harmonic order of the permanent magnet magnetomotive force (u = 1, 3, 5, ...), p is the number of motor poles, ω r is the rotor mechanical angular velocity, reflecting the electrical angular frequency of the spatial harmonics generated by the rotation of the permanent magnet; The armature reaction magnetomotive force calculation is expressed as: The low-frequency and high-frequency harmonics generated by the armature magnetomotive force are denoted as F L (θ,t) and F H (θ,t), where β l and β h represents their phase angle; here, l = 6n ± 1, where n = 0, 1, 2, ..., ω h Corresponding to f h The electrical angular frequency; radial component F Lr (θ,t) and F Hr (θ,t) is denoted as F Lr_v 、F hr_v 、F lt_v and F ht_v , represents the amplitude of the vth harmonic; S1.4 Calculation of electromagnetic force density: The radial and tangential components of the no-load air gap flux density and the armature air gap flux density can be expressed as: B r_no_load (θ,t)=F 1r (θ,t)Λ r (θ)+F 1t (θ,t)Λ t (i) (24) B t_no_load =F 1t (θ,t)Λ r (θ)-F 1r (θ,t)Λ t (i) (25) Therefore, the air gap flux density at the load is: B load (θ,t)=B r_load (θ,t)+jB t_load (θ,t) (28) Where: B r_load (θ,t)=B r_load (θ,t)+B r_arm (θ,t) (29) B t_load (θ,t)=B t_load (θ,t)+B t_arm (θ,t) The calculation of the electromagnetic force includes the calculation of radial and tangential electromagnetic force densities, which is specifically achieved through the following formula: Where, P r (θ,t) and P t (θ, t) are the radial and tangential electromagnetic force densities, respectively, B r (θ,t) and B t (θ, t) are the radial and tangential components of the air gap flux density, μ0 is the vacuum permeability; Where F r (t) and F t (t) are the radial and tangential components of the concentrated force, R s is the inner radius of the stator, θ zi is the position of tooth i, θ st is the stator tooth face width angle.

3. The method for rapid calculation of electromagnetic vibration of a permanent magnet synchronous motor according to claim 1, characterized in that: The step 2 is specifically as follows: 2.1 Equivalent model establishment: The stator core and casing of the permanent magnet synchronous motor are equivalent to a cylindrical shell structure with internal and external ribs. Assume that the average radius of the cylindrical shell is R, the average thickness is h, the axial length is L, and the density is ρ cs , the elastic modulus is E, the Poisson's ratio is υ; the axial, circumferential and radial displacements of any point P on the cylindrical shell surface are u, v and w respectively, considering the eccentricity e of the equivalent tooth structure arn , the displacement expression is: in, and ψ are the axial and circumferential rotation angles, representing the bending deformation of the shell; r is the radius of any point on the cylindrical shell; w cs is the radial displacement of any point P on the cylindrical shell surface; u cs is the axial position of any point P on the surface of the cylindrical shell, where u is the axial basic displacement, R is the average radius of the cylindrical shell, is the axial rotation angle, representing the bending deformation of the shell, Indicates the additional axial displacement caused by bending deformation; v cs is the circumferential displacement of any point P on the surface of the cylindrical shell, where v is the circumferential basic displacement and ψ is the circumferential rotation angle, which represents the bending deformation of the shell. Indicates the additional circumferential displacement caused by bending deformation; 2.2 Potential energy and kinetic energy analysis: According to the three-dimensional elastic theory, the potential energy of a cylindrical shell includes membrane strain energy and bending strain energy, and is expressed as: Among them, δ and ε are normal stress and strain, τ and γ are shear stress and strain, which can be separated into membrane force and bending moment contributions after expansion through geometric equations and constitutive relations; The kinetic energy takes into account the inertia effect of the shell and the expression is: After neglecting high-order small quantities, it is simplified to the contribution of mid-surface motion kinetic energy and moment of inertia: 2.3 Variable separation and mode shape hypothesis: Assuming that the stator vibration is a simple harmonic mode, the variable separation method is used and the displacement form is: Where n is the circumferential modal order, ω is the natural angular frequency, U(z), V(z), and W(z) are the axial vibration mode functions; Considering the axial symmetry of the stator, it is assumed that the axial vibration mode is a sine or cosine function; 2.4 Derivation of frequency equation: Substituting the displacement assumption into the potential energy and kinetic energy expressions, and using Hamilton's principle, we can obtain the following about the amplitude A and V: n 、W n Homogeneous system of equations: Among them, K ij and M ij are the stiffness and mass matrix elements, respectively, obtained by integrating the strain energy and inertia terms in the potential energy and kinetic energy, V n is the amplitude of the axial vibration mode function V(z) corresponding to the circumferential displacement v at the circumferential modal order n; W n is the amplitude of the axial mode function W(z) corresponding to the radial displacement w at the circumferential mode order n; In order for the system of equations to have non-zero solutions, the determinant of the coefficient matrix must be zero, that is: det(K-ω 2 M)=0 After expansion, the modal frequency equation is obtained: ω 4 -ξ1ω 2 +ξ2=0 D=Eh 3 / [12(1-ν 2 )] is the bending stiffness of the cylindrical shell, K = Eh is the membrane stiffness; 2.5 Solve the natural frequency: Solving the frequency equation yields: For low-order circumferential modes (n = 0, 1, 2, ...), ignoring high-order axial vibration (m = 0), the natural frequency of the radial and circumferential coupled vibration is obtained after simplification: in, is the bending stiffness, K=Eh is the membrane stiffness, reflecting the elastic characteristics of the stator structure; ρ cs is the density, h is the thickness, R is the average radius, which reflects the mass distribution; n is the circumferential modal order.

4. The method for rapid calculation of electromagnetic vibration of a permanent magnet synchronous motor according to claim 1, characterized in that: The step 3 is specifically as follows: For mass M n , stiffness K n , damping coefficient C n The mechanical system is subjected to an amplitude F m (t), the electromagnetic vibration process can be expressed as: Then the solution of the system vibration acceleration is: Where: ξ n is the damping ratio corresponding to the n-order circumferential mode, which can be obtained from the empirical formula: F m (t)=2πR i lp m (45) Where R i and l are the inner radius and axial length of the stator equivalent ring respectively; PM is the amplitude of the electromagnetic force density.

5. The method for rapid calculation of electromagnetic vibration of a permanent magnet synchronous motor according to claim 3, characterized in that: The step 2) further comprises: 2.6 Equivalent tooth structure correction: 2.6.1 Potential and kinetic energy analysis of equivalent tooth structure: Considering that the stator teeth are equivalent to elastic beams attached to the inner surface of the cylindrical shell, their potential energy and kinetic energy need to be superimposed on the overall model; let the elastic modulus of the equivalent tooth be E ar , the cross-sectional area is A ar , the moment of inertia around the x and y axes is I xar , I yar , then the additional potential energy is: in, N ar Indicates the number of equivalent teeth, used to sum up the additional potential energy of all equivalent teeth; u arn represents the axial displacement component of the nth equivalent tooth, reflecting the axial deformation of the tooth structure; w arn It represents the displacement component of the nth equivalent tooth in the radial direction, reflecting the deformation of the tooth structure in the radial direction; is the sign of the partial derivative with respect to the axial coordinate z, indicating the rate of change of the corresponding displacement component along the axial direction; L is the axial length of the stator; Additional kinetic energy: P ar is the equivalent tooth density; are the axial and radial velocities, respectively; 2.6.2 Modification of stiffness matrix K: The additional potential energy corresponds to the increment of the stiffness matrix, and the axial stiffness and bending stiffness must be considered: Axial stiffness correction: where φ n is the axial vibration mode function, φ m Radial mode function; Bending stiffness correction ψ n (z): radial mode function, ψ m : axial vibration mode function; Total stiffness matrix correction 2.6.3 Modification of mass matrix M: The additional kinetic energy corresponds to the increment ΔM of the mass matrix, considering the mass distribution of the equivalent tooth: Mode function φ n (z) is consistent with the axial vibration mode, reflecting the distribution of mass along the axial direction; Total mass matrix correction M 修正 =M 原 +ΔM ar 2.6.4 Frequency equation after taking into account the tooth structure: Substitute the modified stiffness matrix and mass matrix into the eigenvalue equation to solve the natural frequency: det(K 修正 -ω 2 M 修正 ) After expansion, the modified modal frequency equation is obtained. By solving this equation, the natural frequency f considering the influence of the tooth structure is obtained. n : K i is the effective stiffness matrix of order i, M eff is the effective mass matrix; By modifying the stiffness matrix K and the mass matrix M, the influence of the tooth structure on the natural frequency is taken into account. That is, in the calculation of the natural frequency of the permanent magnet synchronous motor stator, the influence of the tooth structure must be taken into account by modifying the stiffness matrix K and the mass matrix M.

Images