Vibration analysis of permanent magnet synchronous motor using coupled finite element analysis and optimized meshless method

Conventional finite element analysis (FEA) performed in electromagnetic-vibration coupling calculation of motor suffers from several significant drawbacks, such as large memory space, high computational cost and heavy reliance on mesh quality for accurate solution. With the traditional meshless method, special attention needs to be paid to correctly impose the boundary condition like FEA. Besides, the matrix is easily prone to be ill-conditioned due to introducing large amount of higher-order basis functions. We propose a novel multi-physical coupling method combining FEA and optimized meshless method. The proposed methodology is further evaluated on the vibration analysis of a 12-slot 10-pole permanent magnet synchronous motor (PMSM). Firstly, 2D stator electromagnetic force is simulated and derived based on the local Jacobian derivative method through FEA. The electromagnetic force spectrum is calculated using FFT analysis and further imported into commercial meshless structural simulation software SimSolid for stator harmonic response analysis. Correct force boundary condition and data mapping between meshless and FEA simulation interface are key to the accuracy of the proposed combined multi-physical modeling methodology. This is achieved by introducing new high-dimensional ramp function in the transition region between FEA and meshless domains, which are defined with the shape functions composed of the FEA and meshless method. This function satisfies the continuity and consistency of the displacement function and ensures the convergence of coupled FEA-meshless method. Subsequently, construction of basis function is key to the establishment of convention meshless discrete equation for the elastic problem of rotating machinery. This is designed by using moving least square theory in cylindrical coordinate system. A harmonic response with meshless method is analyzed by using the mode superposition method to obtain detailed mode shape data, acceleration and displacement distribution of stator. Finally, the tangential continuity and robustness are not well considered in the traditional simulation with FEA coupled meshless method. To mitigate this problem, we propose an optimized meshless method based on modified local basis functions to recalculate the harmonic response motion. Then, the coupling electromagnetic-vibration simulation results of traditional coupled FEA-meshless method, optimized coupled FEA-meshless method and complete FEA coupled method are compared. It is worth noting that the optimized method significantly improves accuracy, robustness and computational speed at the same time. In short, the proposed electromagnetic-structure coupling calculation method provides a novel alternative for the multi-physical coupling calculation of rotating machinery combining FEA and meshless simulation methods.

Abstract Conventional finite element analysis (FEA) performed in electromagnetic-vibration coupling calculation of motor suffers from several significant drawbacks, such as large memory space, high computational cost and heavy reliance on mesh quality for accurate solution. With the traditional meshless method, special attention needs to be paid to correctly impose the boundary condition like FEA. Besides, the matrix is easily prone to be ill-conditioned due to introducing large amount of higher-order basis functions. We propose a novel multi-physical coupling method combining FEA and optimized meshless method. The proposed methodology is further evaluated on the vibration analysis of a 12-slot 10-pole permanent magnet synchronous motor (PMSM). Firstly, 2D stator electromagnetic force is simulated and derived based on the local Jacobian derivative method through FEA. The electromagnetic force spectrum is calculated using FFT analysis and further imported into commercial meshless structural simulation software SimSolid for stator harmonic response analysis. Correct force boundary condition and data mapping between meshless and FEA simulation interface are key to the accuracy of the proposed combined multi-physical modeling methodology. This is achieved by introducing new high-dimensional ramp function in the transition region between FEA and meshless domains, which are defined with the shape functions composed of the FEA and meshless method. This function satisfies the continuity and consistency of the displacement function and ensures the convergence of coupled FEA-meshless method. Subsequently, construction of basis function is key to the establishment of convention meshless discrete equation for the elastic problem of rotating machinery. This is designed by using moving least square theory in cylindrical coordinate system. A harmonic response with meshless method is analyzed by using the mode superposition method to obtain detailed mode shape data, acceleration and displacement distribution of stator. Finally, the tangential continuity and robustness are not well considered in the traditional simulation with FEA coupled meshless method. To mitigate this problem, we propose an optimized meshless method based on modified local basis functions to recalculate the harmonic response motion. Then, the coupling electromagnetic-vibration simulation results of traditional coupled FEA-meshless method, optimized coupled FEA-meshless method and complete FEA coupled method are compared. It is worth noting that the optimized method significantly improves accuracy, robustness and computational speed at the same time. In short, the proposed electromagnetic-structure coupling calculation

Introduction
Permanent magnet synchronous motor has been widely used in the fields of electric vehicles and aerospace due to its higher efficiency, larger power density, smaller size, and robust structure [1]. However, due to the lightweight trend of motor stators and rotors, severe vibrations will occur during operation [2]. From the design and manufacture to maintenance of products, numerical simulation of electromagneticstructure coupling analysis plays a very important role in the preliminary design phase of industrial parts [3]. Because of its significant importance, numerous researchers have been focusing on improving the accuracy and efficiency of this procedure for decades. Furthermore, under certain extreme conditions, severe vibrations will result in significant noise, and even damage to the machine [4]. Therefore, the establishment of a fast and accurate multi-physics coupling simulation technology can not only effectively predict and evaluate risks of development, but also reduce costs of research. At present, conventional electromagnetic-structure coupling analysis of electrical equipment includes: experimental method, analytical method, and finite element method [5]. Among them, experimental method requires an actual prototype, which normally introduces high cost and long time due to iterative design; analytical method is widely used in the electromagnetic and vibration calculation [6]. For electromagnetic analysis, the conformal transformation method is employed to solve the two-dimensional magnetic field in the literature [7]. However, the tangential magnetic field component is ignored. In addition, the infinite depth slot model was adopted, leading to the neglect of influence between adjacent slots; hence, it is unable to know the detailed magnetic field distribution of motor [8]. To overcome this problem, the Schwarz-Christoffel mapping is introduced for magnetic field calculation in SPMSM [9]. However, SPMSM stator has many irregular slots, and shape is too complex to get a good boundary correspondence. Hence, most regions require large numbers of tangent polygons to guarantee the Riemann mapping theorem, and this results in complicated mapping functions. In addition, it is greatly worth noting that the improper selection of function parameters may lead to serious crowding defects. The lumped parameter model proposed by Sun to calculate the electromagnetic field is relatively clear in terms of physical concept and is straightforward to be directly applied. However, only regular and linear model can be solved [10]. The skin effect, proximity effect, cogging effect and magnetic nonlinearity are not considered. More importantly, the boundary conditions are too complicate for solving in homogeneous Laplace-Poisson equations of magnetic field due to many sub-domains divided, which makes the theoretical results deviate greatly from the actual situation [11]. For structural vibration analysis, simplified finite-length cylindrical shell model is usually selected to calculate the vibration [12]. However, the tangential electromagnetic force harmonics are ignored. For the fractional-slot concentrated-winding motor investigated in this paper, high-order radial and tangential electromagnetic forces may still lead to large vibration due to large number of rotor poles and slot area [13]. In addition, structure of PMSM exhibits rich, diverse and continuous innovation; hence, it is difficult to obtain an unified analytical model, and only approximate calculation can be realized. Genetic algorithm is employed as an effective tool used for design optimization of electric machines [14]. Although computational convergence becomes easier to be achieved, the multi-physical field coupling calculation of rotating machinery involves many variables. Meanwhile, its objective function has characteristics of multiple extreme points and shows extremely nonlinearity, which leads to the problems such as premature convergence, poor local search ability, low efficiency for solution. Finite element analysis is a robust and mature method. It has been widely used in electromagnetic-structure coupling analysis for electric motors due to its versatility and flexibility characteristics in dealing with boundary problem [15]. However, the finite element analysis is based on the regional variation principle and subdivision interpolation method. It is not easy to modify once the elements are divided and the mesh is generated. Hence, the poor adaptability will affect the accuracy of calculation. In addition, the finite element method involves many regions in meshing, all nodes of the domain and boundary need to be solved simultaneously, leading to many unknown parameters to be solved and large computational burden. Simplified model is a common method to accelerate the solution, such as ignoring bolts and nuts in assemblies [16]. However, in case of improper omission or simplification of some parts in the motor vibration analysis, large error may be generated due to inaccurate modeling of mass and stiffness. Especially for the modal superposition method used in this paper, the resolution of higher-order vibration modes is mostly limited by the quality of the structure meshes [17].The excessively small girds result in large amount of memory and computation time consumed in electromagnetic-vibration coupling analysis [18]. According to the study of Sandia National Laboratories, the creation and edit of analysis solid model, geometry decomposition, meshing and mesh manipulation account for 73% of the overall analysis time costs [19]. Therefore, it is urgent to develop a fast and high precision electromagnetic and structural coupling calculation method for rotating machine with complex shape without tedious mesh generation process.
In recent years, scatter point approximation based on Galerkin meshless method has gained much attention [20][21][22][23][24][25]. It provides accurate numerical calculation even under irregular distribution of discrete points, which is superior to the conventional mesh-based FEA. Thus, it has technical advantages on numerical PDE solution for complex geometrical boundary conditions and provides discrete points distribution insensitive calculation. In addition, it constructs the shape function of field variable approximation by a set of nodes rather than the elements in the FEA, which not only gets rid of the tedious meshes generation process, but also reduces the calculation error caused by the distortion of meshes. Hence, this method has been successfully adopted in solving fluid dynamics, optimization design, and structural vibration problems [26][27][28]. However, the moving least square method is adopted to approximate the trial functions. This method results in one serious problem that the unknown computation quantities are usually the approximate values rather than the real node values [29]. Therefore, the shape function does not satisfy the property of Kronecker function, and the essential boundary conditions cannot be imposed as easily and directly as FEA in dealing with multi-physical field coupling analysis of motor [30]. To overcome this problem, several techniques have been proposed such as Lagrange multiplier method and penalty function method [31][32][33][34]. Although the accuracy of calculation is improved, the unknown parameters increase and lead to an awkward equation structure. The matrix which has to be inverted possesses zeroes on its main diagonal and the coefficients of matrix are no longer positive definite [35]. Hence, it is not good for efficient solution. The modified Lagrange method was adopted by Lopes et al. in the literature [36]. Generally, this method guarantees the symmetrical and positive character of stiffness matrix. However, the accuracy of the solution decreases because the modified function no longer maintains the extreme value property of the original function at the stagnation point. The penalty function method was used to solve the elastodynamic equation by Wang, which needed to solve sequential unconstrained minimization problems [37]. When the number of iterations was large, the condition number of the matrix will be very large. Therefore, the calculation workload is very large, and the convergence speed is slow. In addition, the calculation accuracy is greatly affected by the value of penalty function parameters, leading to matrix to be ill-conditioned due to penalty factor improperly chosen [38]. To overcome disadvantages of meshless methods, the FEA combined with the meshless method is proposed due to more convenient in dealing with the boundary problem by FEA. This is achieved by introducing the new multiple dimensions ramp functions in the transition region between FEA and meshless domains, which are defined with the shape functions composed of FEA and EFG shape functions. Then, the consistency of displacement function is ensured, and the problem of information exchange for the electromagnetic-structural coupling calculation of domain nodes is solved.
Studies that focused on the meshless method coupled with FEA have been conducted by many researchers. In the literature [39], Lluch proposed a fully coupled model of cardiac electromechanics based on smoothed particle hydrodynamics (SPH), and then, hemodynamics and mechanical behavior of the heart were reproduced. In the literature [40], a twoway coupling method of the three-dimensional based on SPH method is presented for handling the fluidstructure interaction problems. However, it is generally known that the traditional SPH scheme for derivative has lower accuracy and numerical stability under the non-uniform particle distributions or near the computing region boundary. In some worse cases, the numerical solution not only deteriorates at the boundary condition, but also has a large error in phase and amplitude because of the worse consistency [41]. To overcome these drawbacks, the conventional SPH method is improved by adding correction function to the governing equation and extended to higher-order Taylor series. Although the stability is increased, the calculation efficiency will be reduced in three-dimensional calculation, and hence, it is not the optimum scheme for the calculation of large rotating machinery such as PMSM. Belytschko and Organ introduced an interface element named ramp function to solve coupling of FEA and meshless, which is defined with the shape functions composed of FEA and meshless shape functions [42]. This method has been used for enforcing the essential boundary conditions to satisfy the property of Kronecker function. Nevertheless, even this method satisfies continuity and consistency of the displacement function variable exactly, the continuity of derivative is ignored. Moreover, the approximate calculation of function and its continuity of derivatives is exactly the technical key of the meshless numerical method for solving partial differential equations. In the literature [43], Dolbow adopted a new ramp function, which realized the continuous transition of shape function and its derivative between FEA and meshless region. However, for the cantilever beam model selected in this article, only the continuity of derivatives in one dimension is considered, while the continuity of tangential direction is ignored. The multi-physical coupling analysis of rotating machinery may contain the first-order spatial derivatives such as gradient, divergence and curl, meanwhile the second-order spatial derivative such as Laplace or Poisson equations [44]. When the support region is small, ignoring the derivative continuity will easily lead to the addition of constraint condition to stress partial differential equation based on the variation principle and the Gaussian divergence theorem, resulting in volumetric locking phenomenon; therefore, the accuracy and convergence speed is affected [45]. Moreover, most meshless method calculations involved in the current papers are limited to the Cartesian coordinate system, and hence, the calculation of motor in this paper also needs to be extended to the cylindrical coordinate system. In the literature [46], a weak form approach of FEAmeshless coupling calculation is proposed, and this is achieved by the help of Lagrange multiplier to build domains between the FEA and meshless method region. It is necessary to perform the inverse operation on the matrix A (shown in Eq. 40) during the process of constructing shape function matrix by using large numbers of higher-order basis functions; hence, the matrix is easily prone to be ill-conditioned. To overcome this problem, the method that number of support points is much larger than basis vector terms adopt by Liu and Gu (N ! m) [47]. Unfortunately, there is no theoretical optimal value of N, which can only be determined through endless numerical experiment calculation, resulting in lower computational efficiency [48]. From the literature review, it is clear that the vibration analysis of three-dimensional PMSM solved by coupled FEA-meshless approach has not attempted so far.
In our study, the electromagnetic-vibration field coupling analysis of a 12-slot 10-pole PMSM was performed combining FEA and meshless method. The focus is on the data transmission of transition region between the electromagnetic field solved by FEA and the structural field solved by meshless method. The structure and organization of this paper are as follows. Firstly, two-dimensional electromagnetic force of stator is calculated by the local Jacobian derivative method through FEA. Afterward, the stator radial electromagnetic force of each order and frequency is summarized through FFT analysis. Secondly, the electromagnetic force was imported into a meshless software named SimSolid as excitation to calculate harmonic response analysis. The shape function of the meshless method in cylindrical coordinate system is constructed through moving least squares theory. This process includes the construction of basis functions in cylindrical coordinates, the selection of weight functions and the discussion of domain. Then, the meshless discrete equation for the three-dimensional elastic problem of rotating machinery is established by coordinate transformation, and the meshless harmonic response motion control equation based on the mode superposition method is derived. With this method, the harmonic response calculation can be realized without ascertaining the precise range and distribution of the nodes. In addition, it should be noted that the node value in the FEA region is real, while the node value of meshless method is virtual; hence, it cannot be matched directly [49]. This is achieved by introducing the transition region between FEA and meshless domains, which are defined with ramp functions constructed by higher-order isoparametric elements. It not only solves the boundary value problem using FEA-meshless method coupling calculation, but also ensures the continuity and consistency of displacement function. The detailed data of stator mode shape, acceleration and displacement distribution are obtained. Thirdly, in view of the problem that the matrix is prone to be ill-conditioned when Joldes adopted moving least squares to calculate the shape function matrix using higher-order basis functions, the optimized moving least squares method is used for correction [50]. The optimized ramp functions of radial and tangential direction are reconstructed by higher-order isoparametric elements, respectively. Based on the previous electromagnetic-structure coupling calculation, the meshless method is adopted to construct the modified basis function of equation. Then, new stiffness matrix is obtained, without introducing any additional nodes or degrees of freedom to recalculate the harmonic response motion in SimSolid software [51]. Finally, the simulation results of traditional coupled FEA-meshless method, optimized coupled FEA-meshless method and complete FEA coupled method are compared, which not only improves the coupling calculation speed of electromagnetic structure, but also verifies the accuracy of optimized scheme. Although there are local differences by comparing the calculation results, it can be found that the amplitude and order are generally consistent. The optimized coupled FEA-meshless method of electromagnetic structure investigated in this paper not only provides a new reference for multiphysical coupling calculation of rotating machinery, but also has important theoretical significance and value for reducing the simulation time of the motor. Hence, the computing potential is improved.

Permanent magnet synchronous motor model
In this paper, a 12-slot 10-pole permanent magnet synchronous motor shown in Fig. 1a is chosen for the study, and the parameters of motor are listed in Table 1. From inside to outside, the whole calculation region can be divided into five parts: rotor core, permanent magnets, air gap, stator core and calculation boundary condition. Figure 1b shows the meshes  of two-dimensional electromagnetic field FEA model. Figure 1c, d shows the three-dimensional structural field FEA model and meshes of motor, respectively. In the calculation of the two-dimensional transient electromagnetic field, the value of field usually changes with coordinate. Whereas the electromagnetic field components of conventional three-node triangular element are constant, and hence, it is difficult to be employed in the calculation of rapidly changing field [52]. In order to ensure the calculation accuracy, huge calculation data will generate due to large numbers of meshes employed, which is the barrier to improvement of computational efficiency. Therefore, a higher-order element such as six-node triangular is utilized due to its higher-order polynomial with relatively small number of elements.
For three-dimensional structural fields, due to the electromagnetic force excitation and the transition region are applied to the tooth tips. Therefore, the mesh selection of tooth tips appears to be particularly important for the calculation accuracy. In this paper, two types of elements are employed, named tetrahedral and hexahedral element, respectively [53,54]. The hexahedral element has obvious advantages such as higher distortion resistance, higher calculation efficiency, and smaller number of meshes [55]. However, the structure of rotating machinery proposed in this paper has more complex shape, and hence, the tetrahedral element is more suitable to realize the automatic meshing and self-adaptation for arbitrarily complex geometries [56]. Therefore, in order to make full use of the advantages of hexahedral and tetrahedral element, the method combining with two types of elements is employed to achieve the optimized of computational efficiency. In particular, the eight-node hexahedral element is chosen in the radial region of tooth tips set at the tooth tips of stator due to lots of curved surfaces. By contrast, the tennode tetrahedral element is more suitable for complex geometric shapes with narrow chamfers and slender faces, and hence, it is selected as the mesh for the tangential region of tooth tips. To satisfy the continuity and consistency of nodes between tetrahedral and hexahedral element region, the surface quadrilateral element of hexahedral element is transformed into triangular element. Thus, it can be directly connected to the triangular elements on the tetrahedral surface. For other parts of stator, the eight-node hexahedral element is employed. The total number of elements in electromagnetic field is 753924, while in structural field is 12888098. Therefore, meshes generation will take up a lot of calculation time. This is also the problem that the meshless method proposed in this paper aiming to solve.
The process of electromagnetic-vibration coupling analysis is shown in Fig. 2a. First, the radial and tangential electromagnetic force is calculated for temporal and spatial using FEA in two-dimensional electromagnetic field. Next, the radial and tangential electromagnetic force is converted into a frequencyresponse curve by applying the FFT. Then, the results were applied to three-dimensional FEA and meshless method structure field as excitation source to calculate harmonic response analysis. The transition region between FEA and meshless domain is shown in Fig. 2b. Finally, the calculated results were compared with the meshless method to verify the accuracy.
The motor material and their properties are shown in Table 2. All the mechanical material properties for all motor components were considered as homogeneous and isotropic [57].
3 Electromagnetic vibration coupling analysis based on finite element method

Electromagnetic field simulations
In this paper, a 12-slot 10-pole surface-mounted PMSM established by FEA is used as an example. The simplified two-dimensional model is shown in Fig. 1, and the parameters are shown in Table 1. To improve the analysis efficiency and reduce the analysis cost, following assumptions were used in the analysis [58,59]: (1) The unipolar permanent magnet is divided into three pieces circumferentially, and the permanent magnet is taken as uniformly magnetized; thus, the relative permeability of permanent magnet is equal to 1; (2) Due to the small magnetic flux intensity of motor core in external region, the external magnetic field of motor is ignored, and the   outer region of stator is assumed to be ideally magnetically insulated; (3) The motor core is assumed to be infinite length, and end winding flux leakage is ignored. In addition, the vibration and deformation effect on magnetic field are not considered; (4) The turn-to-turn insulation of slots is neglected, namely all windings of each slot are treated as excellent homogeneous conductors; (5) The electromagnetic field model is established without considering the hysteresis effect. Besides, the material of motor core is treated as infinitely permeable; (6) In the pre-processing of electromagnetic field simulation, the influence of harmonics introduced by the driver and control method is not considered. The ideal three-phase symmetrical alternating current is adopted as excitation source.
The magnetic vector interface and boundary conditions of transient electromagnetic field equation of motor are [60]: The finite element discretization equation of electromagnetic field is established by the weighted residual approach, which is shown as: ZZ Regarding the surface integral as the sum of unit integrals, Eq. (2) is discretized. The first term of Eq. (2) is: The second and third term of Eq. (2) is: where A u is the magnetic potential, l e is the magnetic permeability, u is the electric potential, x is the angular frequency, r e is the conductor conductivity, and Jz is the current source density. N E is the shape function in the two-dimensional electromagnetic field composed of six-node triangular plane unit [61], expressed in Eq. (6). The schematic diagram of node in area coordinate is shown in Fig. 3.
where L 1 -L 3 can be expressed in area coordinate, 4 is the area of a triangle, the parameter a, b and c can be expressed as Eq. (9).

Electromagnetic force analysis
The electromagnetic force of nodes is expressed as [62]: The electromagnetic force transformed via isoparametric method can be expressed as: where J E is the Jacobian determinant after the change of coordinate shown in Eq. (2). The parameters g, f and n are in correspondence to the radial, circumferential and axial directions, respectively.
The radial and tangential electromagnetic forces are, respectively [63]: 3.3 The basic equations of 3D elastoplasticity problems in FEA When the three-dimensional FEA is used to analyze the structural dynamics of the stator, the following assumptions are often adapted [64,65]: (1) Because the vibration is mainly caused by the electromagnetic force of stator, the influence of the rotor system is not considered; (2) The minor physical structures on motor that exert negligible influence on the structural strength, such as wedge holder near the slot opening, were neglected (3) The electromagnetic force is applied to the inner circular surface of stator as excitation source to calculate harmonic response analysis. In addition, it is considered that the electromagnetic 1 2 3 P Fig. 3 The schematic diagram of area coordinate force remains constant along the z-axis direction; (4) The outer surface region of stator is fully constrained, and the displacements of nodes in three directions are constrained.
The equilibrium equation of three-dimensional elastoplasticity problems of cylindrical coordinate system is [66]: The stress-strain relationship is: The geometric equation is: where e and u are the strain vector and the displacement vector, respectively, and L is the derivative operator matrix of cylindrical coordinate system: For any point x in the domain, the displacement vector using can be expressed as [67]: From Eqs. (16) and (18), we have: The displacement boundary condition is [49]: The force boundary condition is [68]: where r is the vector of stress rate; E and l are elastic modulus and Poisson's ratio, respectively; C u and C r are the displacement and force boundary, respectively, C is boundary of domain X and C = C u [ C r ; f q , f h and f z are the vector of body force in correspondence to the radial, circumferential and axial directions, respectively; t q , t h and t z are the prescribed traction rate on boundary C r in correspondence to radial, circumferential and axial directions, respectively; and n q , n h , and n z are the normal direction components of the boundary, respectively; The stress components for the displacement u q , u h and u z in cylindrical coordinate system are: According to weighted residual method, the corresponding Galerkin weak form of elastoplasticity problems is [69]: Substituting Eqs. (15) and (16) into Eq. (24): Substituting Eqs. (15) and (19) into Eq. (25): then, the final discrete equation in FEA form without damping matrix can be obtained as [49,70]: where M FEA is element stiffness matrix, K FEA is element stiffness matrix, and t t (x) is concentrated force that acts on the stress boundary.
where B s is the strain displacement transformation matrix of the cylindrical coordinate system in FEA. The element strain e is derived from the global coordinates. The converted Jacobian matrix J S of natural coordinates is introduced as: The vibration response can be expressed as a function of the mode and the excitation force [71]: where W n is the nth mass-normalized mode vector; and x and p n are, respectively, angular frequencies of the force.

Moving least squares approximation
The nodal displacement function u is approximated over a series of discrete nodes with the MLS approximation as follows [72]: where p(x) is the polynomials basis function for threedimensional problems [73]: It is represented in cylindrical coordinates as: The coefficient a(x) is determined by minimizing the weighted discrete norm, defined as: where u J ¼ u J ðxÞ ¼ uðx J Þ, w J ðxÞ is weight function, shown as [74]: The deviation of discrete scale in each dimension can be ignored for plasticity problem calculation of cylinder in meshless method, so that the discrete nodes are distributed more uniformly and leading to calculation relatively simple. Therefore, circular support domain is used to construct the weight function [74]: The minimum of J can be achieved by [75]: where matrices A(x), B(x) and P are [76]: . . .
Finally, the approximation function is provided as: The shape function constructed by meshless method is [77]: The final discrete equation in meshless method form without damping matrix can be obtained as: where M meshless is the equivalent mass matrix, K meshless is the stiffness matrix, B m is the strain matrix, and L is the derivative operator matrix shown in Eq. 17.

Coupling of finite element and meshless method
In terms of calculation of the electromagnetic-structure coupling shown in Fig. 2, the whole domain is divided into three sections, which is FEA domain X FEA , meshless domain X meshless , and transition domain X INF , respectively. The FEA domain has been calculated in the third chapter. In order to coordinate the displacement function, a ramp function is applied in the transition domain, leading to the displacement of nodes closed to FEA boundary consistent with FEA domain [78]. Similarly, the displacement of nodes closed to meshless method boundary is consistent with meshless method domain. Hence, the continuity and consistency of displacement function are satisfied, and the convergence of coupled FEA-meshless method is ensured. The displacement vector u INT (x) of any point in the transition domain is [79]: The shape function of the transition domain can be expressed as [79]: where R (x) is the ramp function and can be expressed as: The derivative of ramp function is: The new form of discrete equation is [80]: The whole domain is divided into three regions: FEA domain, meshless domain, and transition domain; hence, the stiffness matrix can be divided into three integral units:

Optimized meshless method
There are large numbers of matrix inverse calculation in the process of constructing the approximate function. In order to accelerate convergence, the usual approach is to increase the node density of field and reduce the scale of support domain. Therefore, the matrix A is easily prone to be ill-conditioned. This is improved by introducing an optimized meshless method named moving least squares core (MLSC) [48]. Compared with the traditional MLS, we propose a new local basis function (core basis) to redefine the approximate function. This is achieved by defining the local coordinate system n = x-x I . For the three-dimensional model, the core basis of new complete polynomial can be expressed as [81]: ðx À x I Þ; ðy À y I Þ; ðz À z I Þ ; ðx À x I Þ Á ðy À y I Þ; ðy À y I Þ Á ðz À z I Þ; ðz À z I Þ Á ðx À x I Þ; ðx À x I Þ 2 ; ðy À y I Þ 2 ; ðz À z I Þ 2 2 6 6 6 4 The polynomial approximation function of MLSC can be expressed as: For the scatter aggregate: x J f g N J¼1 ð8x J 2 X S ðx I ÞÞ, the new moving least squares weighted functional can be expressed as: when oJ=oa ¼ 0, the new coefficient is: The new shape function is: In order to ensure the continuity of the function, a new ramp function R(x) is constructed on the basis of Eq. (51), which is taken as the sum of FEA shape functions of nodes located at the meshless boundary region considering both tangential and radial region. The ramp function of tangential region is constituted by ten-node quadratic tetrahedral element, and the radial region is eight-node hexahedral element, expressed as: Þði ¼ 1; 2; :::7; 8Þ; 5 Calculation results

Electromagnetic field calculation result
The analytic prediction of electromagnetic-vibration coupling calculation in PMSM contains the electromagnetic force calculation, the modal calculation, and the harmonic response calculation. The electromagnetic force acting on the stator tooth tips is the main excitation source of vibration. Therefore, the first step for coupling analysis is to determine the force accurately. Due to the slotting effect and the irregular air gap of the proposed motor, it is difficult to formulate electromagnetic force by the analytical method. Therefore, commercial FEA software ANSYS Maxwell is adapted to calculate the transient electromagnetic field of PMSM. The initial degree of motor is set as 0°. Meanwhile, the termination time of solution is set as 2.67 ms, and the step size is set as 2.67 9 10 -3 ms. The magnetic flux field distribution of motor at 2.67 ms are shown in Fig. 4. As can be seen from Fig. 4a, the main magnetic field of motor shows a trend of symmetrical distribution, while there is a mutation at the junction of the tooth and slot, and the maximum value is 2.91 T. As can be seen from Fig. 4b, the flux distribution is unevenly distributed, mainly concentrated in the iron core. It can be seen that most of the magnetic lines at the air gap enter the stator teeth. The main flux that passes the air gap only makes up a tiny portion among all flux lines produced by armature currents excitation, while the majorities are the slot leakage flux and the tooth-top leakage flux [59]. Some of the flux lines exit from the tooth top and thus result in radial electromagnetic force on the stator tooth tips. Meanwhile, other flux lines exit from tooth edges and are more concentrated near the slot opening. This part of the flux results is the main cause of tangential electromagnetic force. Taking the annulus of tooth tip inside the stator as observation path, the temporal-spatial distribution of two-dimensional electromagnetic force density is shown in Fig. 5  the largest change is at the cogging junction. In terms of spatial distribution, both radial and tangential force waves are mainly divided into four cycles along the circumferential direction, and hence, the dominant order of space is 4th order. In terms of temporal distribution, both radial and tangential force waves change periodically with the rotation of the rotor. Figure 5c presents the radial and tangential electromagnetic force density at 2.67 ms, the maximum value of the radial and tangential electromagnetic forces are 8.14 9 10 5 N/m 2 and 3.13 9 10 5 N/m 2 , respectively. Although the amplitude of tangential electromagnetic force wave is smaller, only 38.5% of the radial component. However, the tangential electromagnetic force of stator will still mutate due to the larger fluctuation of positive and negative directions. In addition, frequency is higher under high-speed conditions, and hence, the vibration caused by tangential electromagnetic force cannot be ignored either.
Frequency analysis of Fig. 5c acting on tooth tips was performed using the FFT analysis, illustrated in Fig. 6. It can be seen that the harmonic content contained of the radial and tangential electromagnetic force waves is nearly the same, the amplitudes peak when the frequency is f = 2f 0 = 1500 Hz. Since the flux weakening control of PMSM could result in serious current harmonics distortion. Therefore, the harmonic components such as 4f 0 , 6f 0 and 8f 0  According to the literature [82], the time order can be expressed as: n 1 AE n 2 j j; l 1 AE l 2 j j ; n 1 AE l 2 j j , where n and l are integer. Due to there are only odd non-third harmonic components, l and n are both odd numbers, resulting in an even number of times regardless of the radial or tangential electromagnetic force.
The force density obtained by Maxwell stress method is not convenient to load onto the stator tooth tips for further vibration calculation. Instead, force on the element nodes is easily employed to represent the specific distribution of local force. Figure 7 is the twodimensional electromagnetic force of stator teeth mapped to the three-dimensional structure field at 2.267 ms. It can be seen that electromagnetic force is mainly distributed on the surface of stator and the maximum value is 7.89 N/m 3 . In addition, it can be clearly seen that the electromagnetic force on the air boundary of rotor is much greater than that inside the core. The reason is that without considering the internal force caused by the deformation of the medium, the change of permeability in the core is small, while the permeability changes sharply at the core-air boundary.

Modal analysis results
In case of electromagnetic force distribution and corresponding frequency confirmed, the distribution of natural frequency determines the modulus of structural vibration, and the corresponding mode shape determines the distribution of structural displacement. Therefore, accurate calculation of stator mode is the premise of harmonic response analysis. Table 3 presents the part former 20th partial natural frequencies of stator calculated by FEA and meshless method, respectively. Accordingly, Table 4 presents the partial mode shapes. It can be seen from two tables that the natural frequency distribution is dense due to stator structure model with abundant degrees of freedom (DOF). Hence, when the natural frequency changes in a small range, the mode shape will change obviously. It can be found that the natural frequency calculated by two meshless methods is nearly the same. It is also worth noting that the natural frequency calculated by the meshless method is greater than that by FEA generally, while the errors of 5th and 7th order are relatively greater, namely 5.9% and 6.3%, respectively. The main reason is that when the FEA method is applied in the modal calculation, the dynamic coupling between the secondary DOF and the main DOF is omitted. On the contrary, the two meshless method improves this deficiency, so the calculation results tend to be smaller, the detailed proof process can be referred to the literature [83]. In addition, it can be seen from Table 4 that both FEA and meshless method calculation results express the same vibration trend. For the lower-order mode shapes, the local variation is dominant. The specific description is expressed as follows. The first mode shapes are mainly the first-order vertical bending deformation of upper right part teeth, while the third mode shapes are mainly the first-order vertical bending deformation of upper left part teeth. It can be found that the deformation expresses larger in the fifth-order mode shape. For the higher-order mode shapes such as10th, the vibration is relatively obvious, which is a second-order bending mode. At this frequency, the modal deformation corresponding to the stator has become more complicated.

Harmonic response analysis
To confirm the influence of the tangential electromagnetic force on vibration, vibration simulation is carried out. Figure 8a, b expresses displacement and acceleration calculated by the traditional meshless method only considers radial electromagnetic force as excitation source, respectively. It is compared with the calculation results of the FEA shown in Fig. 9 which considers both radial and tangential electromagnetic Body Force Density(N/m 3 ) 7.89×10 10 7.01×10 10 6.14×10 10 5.26×10 10 4.38×10 10 3.51×10 10 2.63×10 10 1.75×10 10 8.76×10 9 7.21×10 -2 Fig. 7 Electromagnetic force distribution forces. It can be seen that when the radial and tangential electromagnetic force excitation is both taken into account, the displacement and acceleration of the stator surface are significantly greater than that considering radial force alone. The maximum amplitude of displacement is 4.06 9 10 -5 mm and 7.64 9 10 -5 mm, respectively, with the error of 47%. Meanwhile, the maximum amplitude of acceleration is 3.42 m/s 2 and 6.78 m/s 2 , respectively, with the error of 49.56%. To evaluate the accuracy of traditional meshless method, a comparison was made between the displacement and acceleration at which peak value occurs. Figure 10a, b expresses displacement and acceleration calculated by the traditional meshless method considering both radial and tangential electromagnetic forces, respectively. Compared with Fig. 9, the maximum amplitude of displacement is 7.09 9 10 -5 mm and 7.64 9 10 -5 mm, respectively, with the error of 7.2%. Meanwhile, the maximum amplitude of acceleration is 8.45 m/s 2 and 6.78 m/s 2 , respectively, with the error of 24.6%. Thus again proves that the tangential electromagnetic force should not be ignored in the process of harmonic response analysis. It is worth noting that the error of acceleration is too large to meet the requirement of calculation accuracy. Hence, the optimized of traditional meshless method should be proposed urgently.
In order to validate the simulation results against the vibration analysis of PMSM using coupled FEA and optimized meshless method. We assume that both traditional meshless method and optimized meshless method share the same structure model with complete FEA. In addition, the harmonic response analysis is calculated in the same condition as complete FEA coupling method. Figure 11a, b expresses the calculation results of the optimized meshless method. It can be found that the maximum amplitude is 7.67 9 10 -5 mm. Compared with the FEA, the error has been reduced to 0.4%. Accordingly, the maximum acceleration is 6.77 m/s 2 , and the error is 0.15%. It is worth noting that the calculation accuracy is greatly improved. It can also be found that the vibration distribution trends in Figs. 9 and 11 are nearly the same. Due to the greater external constraints of the stator, the stiffness of stator is large, and hence, the main deformation concentrates on the stator teeth and yoke.
In order to demonstrate the advantages of double meshless methods in calculating time. The time of each step of the electromagnetic-structure coupling spent on calculation is given in Table 5. It can be found 7.64×10 -5 7.00×10 -5 6.36×10 -5 5.73×10 -5 5.09×10 -5 4.45×10 -5 3.82×10 -5 3.18×10 -5 2.55×10 -5 1.91×10 -5 1.27×10 -5 6.36×10  that the FEA generates a large number of meshes, and the calculation time far exceeds the meshless method, especially at the harmonic response analysis step. The simulation time with complete FEA is about 43 times longer than the simulation time with the proposed method combining FEA and traditional meshless method. Even the calculation workload increases due to the reconstruction of the basis function of optimization meshless method, leading to certain increase in time. The simulation time with complete FEA is about 30.7 times longer than the simulation time with the proposed method combining FEA and optimized method. Hence, it can still save a lot of calculation time.

Conclusions
This paper proposes a new method for electromagnetic-structure coupling calculation of rotating machinery with complex shape and boundary conditions. It is achieved by combining FEA and meshless method. On this basis, the method of combining FEA and meshless method is optimized, and the following conclusions are drawn: (1) The electromagnetic force is calculated by the local virtual displacement method, and the radial and tangential electromagnetic forces are obtained, respectively. Then, the source of electromagnetic force of each order and frequency is summarized through FFT analysis. It is found that the influence of tangential electromagnetic force on vibration cannot be ignored. In addition, by comparing the harmonic spectrum of electromagnetic force waves of tooth tips, it can be seen that both radial and tangential   (2) In order to achieve the data mapping between FEA and meshless method, a ramp function composed by shape functions of FEA and the meshless method is constructed to enforce the essential boundary conditions as directly as FEA. Thus, the continuity of the function is ensured. The final calculation results show that the traditional meshless method can greatly reduce the calculation time, especially in the step of mesh division and harmonious response; (3) In order to mitigate the problem that the tangential continuity and robustness are not well considered in the vibration analysis using coupled FEA and tradition meshless method. This paper proposes a moving least squares core approximation method for constructing local basis functions, which constructs a new stiffness matrix and recalculates the harmonic response motion. By comparison with FEA, it has been confirmed that the optimized meshless method proposed can be used to calculate the vibration. This paper provides a novel alternative for the multi-physical coupling calculation of rotating machinery.