The performance of electric motors, especially in applications such as automotive, aerospace, and industrial systems, is heavily influenced by their electromagnetic vibration and noise. These aspects not only affect the operational efficiency but also the comfort and reliability of motor-driven systems. With the rapid development of electric vehicles and advanced industrial automation, there is an increasing demand for motors that are efficient, compact, and silent. Therefore, the accurate prediction and control of their electromagnetic vibration and noise have become paramount. A significant body of research has been dedicated to understanding and mitigating the sources of electromagnetic vibration and noise in permanent magnet synchronous motors (PMSMs) and other types of electric motors.

Studies by Ballo et al. [1] and Xing et al. [2, 3] have focused on developing simplified analytical models and numerical prediction models to forecast the noise and vibration in PMSMs at the design stage. The influence of electromagnetic forces on motor vibration [4, 5] has also been a central theme. Strategies to mitigate vibrations and noise have been explored [6, 7]. Experimental studies by Torregrossa et al. [8] and Zhao et al. [9] have validated theoretical models and numerical simulations, providing insights into the effectiveness of various vibration reduction methods. The detection and analysis of oscillations in rail vehicle systems has been significantly advanced [10–12], particularly focusing on pantograph control, which not only improved the understanding of signal processing in this context but also provided practical solutions for real-time applications, enhancing the safety and efficiency of rail transportation. A multidisciplinary approach by Chai et al. [13] and Wu et al. [14], combining electromagnetics, structural mechanics, and acoustics, has been employed to provide a comprehensive understanding of motor behavior. The dynamic behavior of the rotor and the acoustic performance of the entire motor system have been examined [15, 16]. Optimization studies by Mendizabal et al. [17] have provided guidelines for designing low-noise motors. Certain studies have focused on specific aspects such as the effects of laminations [18], axial forces [19], and the application of amorphous alloys in stators [20], offering specialized insights into motor design.

The finite element method (FEM) has been extensively employed by Mao et al. [21] and Wang et al. [22] to predict acoustics, fracture mechanics, electromagnetics, and vibrations. Additionally, [23, 24] investigated the natural frequencies of motor components. However, there are several problems with FEM when modeling infinite domains. The boundary element method (BEM) has been used to tackle potential problems because it offers good accuracy and easy mesh construction. Particularly for exterior acoustic problems, the Sommerfeld radiation condition at infinity is rapidly satisfied [25]. The boundary integral problem has been quantitatively solved using the Galerkin approach for BEM implementation [26, 27].

In order to directly resolve the equation system, the conventional boundary element method (CBEM) produces a dense and non-symmetric coefficient matrix that takes a long time to compute. The fast multipole method (FMM) [28–30], the fast direct solver [31, 32], and the adaptive cross approximation approach [33] are only a few of the methods that have been employed to expedite the resolution of the integral issue. Architects and designers are increasingly considering changing the structural geometry to reduce noise. There is much potential for radiated noise reduction with this structural-acoustic optimization [34–36]. FEM and BEM may be employed with some computer-aided engineering (CAE) software. However, contemporary CAE requires that the models produced by CAD software be converted into simulation-ready models as part of the preprocessing phase. The transfer of geometric model data by the CAE results in geometry errors. One proposed approach to this problem [12, 37] is to combine BEM with geometric modeling and numerical simulation using isogeometric analysis (IGA) [38–40]. IGABEM has been employed to tackle an extensive variety of problems, including elastic mechanics [41], potential problems [42–46], heat transfer problems [47], wave propagation [48–53], fracture mechanics [54], electromagnetics [55–60], and structural optimization [61–66].

In this study, the acoustic analysis under electromagnetic vibration is solved using the fast multipole boundary element method (FMBEM). With regard to the advantages of FMBEM over CBEM, please refer to Chen et al. [67].

A cylindrical shell with radius a and infinite length is the subject of the investigation. Consider a region of length l in the model where the cylindrical shell vibrates and outside of which there are no displacements of the cylindrical shell. The acoustic vibration coupling is not taken into account while analyzing the sound field. The region outside the cylindrical shell when r > a is the acoustical analysis domain. The model is presented in Figure 1.

Assuming a known radial displacement w in Equation 1: w z , θ , t = W m γ m z cos n θ e − j ω t , (1) where γ m ( z ) is the modal function in the axial direction, cos ( n θ ) is the modal function in the circumferential direction, W m is the amplitude, and j is the imaginary unit. After that, the time term e − j ω t is omitted from the calculation in the frequency domain.

Equation 1 is derived for time t in order to obtain the velocity expression, as shown in Equation 2. w ̇ z , θ = − j ω W m γ m z cos n θ . (2)

The sound field in the domain must follow Equation 3 after being excited by the displacement in Equation 1. p r , z , θ = P r γ m z cos n θ . (3)

A Fourier transform of Equations 1, 3 in the z direction yields Equations 4, 5. w ̃ k z , θ , t = W m γ ̃ m k z cos n θ , (4) p ̃ r , k z , θ = P r γ ̃ m k z cos n θ , (5) in which k z = ω / c z is the wave number, and we have Equation 6. γ ̃ m k z = ∫ − ∞ ∞ γ m z e − j k z z d z . (6)

In the cylindrical coordinate system, the acoustic Helmholtz equation has the form shown in Equation 7. ∂ 2 p ∂ r 2 + 1 r ∂ p ∂ r + 1 r 2 ∂ 2 p ∂ θ 2 + ∂ 2 p ∂ z 2 + k f 2 p = 0 , (7) where the wave number k f = ω / c f .

A Fourier transform of Equation 7 in the z direction yields Equation 8. ∂ 2 p ̃ ∂ r 2 + 1 r ∂ p ̃ ∂ r + k f 2 − k z 2 − n r 2 p ̃ = 0 . (8)

Substituting Equation 5 into Equation 8 yields Equation 9. ∂ 2 P r ∂ r 2 + 1 r ∂ P r ∂ r + k f 2 − k z 2 − n r 2 P r γ ̃ m k z cos n θ = 0 . (9)

We need Equation 10 to make Equation 9 constant. ∂ 2 P r ∂ r 2 + 1 r ∂ P r ∂ r + k f 2 − k z 2 − n r 2 P r = 0 . (10)

The homogeneous equation shown in Equation 10 is a Bessel equation with the solution in Equation 11. P r = A J n k r r + B Y n k r r , k r = k f 2 − k z 2 1 / 2 . (11)

When r → ∞ , the radiation field tends to be in the form of a plane wave, as shown in Equation 12. P r → A e j k r r − j ω t = A cos k r r + j ⁡ sin k r r e − j ω t . (12)

Furthermore, based on the approximation of two classes of Bessel functions at infinity in Equation 13, lim k , r → ∞ J n k r r → 2 π k r r cos k r r − 2 n + 1 4 π , lim k , r → ∞ Y n k r r → 2 π k r r sin k r r − 2 n + 1 4 π . (13)

B = j A is needed to satisfy the infinity condition. Then, we have Equation 14. P r = A J n k r r + j Y n k r r = A H n 1 k r r , (14) where H n ( 1 ) denotes the Hankel functions of the n -th order first kind.

The continuity conditions shown in Equation 15 must be satisfied at the interface: − 1 j ω ρ f ∂ p ∂ r r = a = w ̇ | r = a . (15)

Performing a Fourier inverse transform of Equation 15, we have Equation 16. ∂ p ̃ ∂ r = − ω 2 ρ f w ̃ . (16)

Substituting Equations 4, 5 into Equation 16 yields Equation 17. A = − ω 2 ρ f W m k r H n 1 ′ k r a , (17)

and then, we have Equation 18. p ̃ = − ω 2 ρ f H n 1 k r r k r H n 1 ′ k r a γ ̃ m k z W m ⁡ cos n θ . (18)

Performing a Fourier inverse transform of Equation 18, we finally obtain an analytical solution for the sound pressure at any point in the exterior acoustic domain for the cylindrical shell model, as shown in Equation 19. p r , z , θ = − ω 2 ρ f 2 π ∫ − ∞ ∞ H n 1 k r r k r H n 1 ′ k r a γ ̃ m k z e j k z z d k z W m ⁡ cos n θ . (19)

Note that r > a .

Consider the Helmholtz governing equation in time-harmonic acoustic analysis, as shown in Equation 20. ∇ 2 p x + k 2 p x = 0 , ∀ x ∈ Ω , (20) where ∇ 2 is the Laplace operator, p ( x ) is the sound pressure at field point x , k = ω / c is the wave number, ω is the angular frequency of the incoming wave, c is the sound speed in the domain Ω , and Ω is the domain for acoustic analysis. In the frequency domain, the time-dependent component e − j ω t can be excluded from the computation.

Applying Green’s second theorem to the Helmholtz equation yields the following integral equation, as shown in Equation 21. p x + ∫ Γ ∂ G x , y ∂ n y p y d Γ y = ∫ Γ G x , y q y d Γ y + p i n c x , y ∈ Γ , (21) where Γ = ∂ Ω is the border of the domain Ω , y is the sound source point on the boundary Γ , x is the field point affected by the sound source y , G ( x , y ) is the Green’s function, n ( y ) is the outward unit normal vector at point y , ∂ ( ) / ∂ n = ∇ ( ) ⋅ n is the exterior derivative, p inc is the sound pressure of the incident wave, and q ( y ) is the outward flux at point y . q ( y ) satisfies the Neumann boundary condition, as shown in Equation 22. q y = ∂ p y ∂ n y = j ρ ω v f y , (22) where ρ is the density of the medium in which the sound wave propagates, and v f ( y ) is the normal velocity at point y .

To obtain the unknown sound pressure on the boundary Γ , the field point x may converge to the boundary, and thus, the Kirchhoff–Helmholtz conventional boundary integral equation (CBIE) is obtained from Equation 21, as shown in Equation 23. c x p x + ∫ Γ ∂ G x , y ∂ n y p y d Γ y = ∫ Γ G x , y q y d Γ y + p i n c x , x , y ∈ Γ . (23)

The normal derivative boundary integral equation (NDBIE) of Equation 23 is given by Equation 24. c x q x + ∫ Γ ∂ 2 G x , y ∂ n y ∂ n x p y d Γ y = ∫ Γ ∂ G x , y ∂ n x q y d Γ y + ∂ p inc x ∂ n x , x , y ∈ Γ . (24)

In Equations 23, 24, c ( x ) = 1/2 when the boundary around point x is smooth.

Applying only CBIE or NDBIE leads to non-uniqueness of the solution to the exterior sound field analysis, which can be solved by linearly combining CBIE and NDBIE. This is called the Burton–Miller method [68], as shown in Equation 25. c x p x + α c x q x + ∫ Γ ∂ G x , y ∂ n y p y d Γ y + α ∫ Γ ∂ 2 G x , y ∂ n y ∂ n x p y d Γ y = ∫ Γ G x , y q y d Γ y + α ∫ Γ ∂ G x , y ∂ n x q y d Γ y + p i n c x + α ∂ p inc x ∂ n x , x , y ∈ Γ , (25)

in which α is the coupling coefficient, and we have Equation 26. α = j , k < 1 , − j k , k ≥ 1 . (26)

In Equation 25, for 2-D acoustic analysis, we have Equation 27. G x , y = j 4 H 0 1 k r , ∂ G x , y ∂ n x = − j k 4 H 1 1 k r ∂ r ∂ n x , ∂ 2 G x , y ∂ n y ∂ n x = j k 2 4 H 2 1 k r ∂ r ∂ n y ∂ r ∂ n x + j k 4 r H 1 1 k r n i y n i x , (27)

and for the 3-D acoustic problem, we have Equation 28. G x , y = e j k r 4 π r , ∂ G x , y ∂ n x = − e j k r 4 π r 2 1 − j k r ∂ r ∂ n x , ∂ 2 G x , y ∂ n y ∂ n x = e j k r 4 π r 3 3 − 3 j k r − k 2 r 2 ∂ r ∂ n y ∂ r ∂ n x + 1 − j k r n i y n i x , (28) where H n ( 1 ) is the Hankel function of n -th order first kind, r is the distance between point x and y , r = | x − y | , and n i is the component of the coordinate.

In Equation 25, discretizing the boundary Γ yields the linear equations shown in Equation 29. H p = G q + p i , (29) where p i is the nodal pressure from p inc .

The sound pressure at point x could be obtained by solving Equation 29.

In this part, we will calculate the sound pressure at some certain location using the analytical equation (Equation 19) and numerical analysis (using BEM), respectively. Both are implemented using our in-house Fortran code. The algorithm is crafted using the Fortran 90 programming language and compiled with the combination of Visual Studio 2022 and Intel^®oneAPI 2022 toolkit. It is executed on a PC with an Intel(R) Xeon(R) Bronze 3204 CPU @ 1.90 GHz and complemented by 128.0 GB of RAM. In this part, the coordinates of the points are described using the cylindrical coordinate system in Equation 19.

Consider the cylindrical shell model shown in Figure 1. The cylindrical shell of radius a = 0.5 m and length l = 0.4 m is closed at both ends for the three-dimensional (3-D) numerical analysis. Consider Equation 1. Let the vibration of the cylindrical shell be uniformly distributed in the axial direction; that is, let γ m = 1. Let the amplitude of the vibration W m = 2.24 × 1 0 − 8 m, the frequency f = ω / 2 π = [100, 1,000] Hz, and the order of the circumferential vibration order be n = 0.

The cylindrical shell model is meshed with quadrilateral elements. Figures 2A, B gives the meshing scheme with element number of 2112 and 8320, separately. Figure 3 provides the numerical solutions of the sound pressure at ( 2 / 2 , 0.0, π / 4 ) — (0.5, 0.5, 0.0) in the cartesian coordinates for the meshing schemes in Figure 2, where the center of the cylindrical shell is at (0.0, 0.0, 0.0). The comparison in Figure 3 is crucial, as it tests the sensitivity of the numerical solution to the mesh density. A fine mesh provides a more accurate representation of the continuous problem but at a higher computational cost. Conversely, a coarse mesh reduces computational demand but may lead to inaccuracies. The two mesh schemes have 2,112 and 8,320 elements, respectively. The agreement between these two schemes suggests that the coarser mesh (2,112 elements) is sufficiently refined for the problem at hand, capturing the essential features of the sound pressure distribution without the excessive computational cost associated with a very fine mesh. Good agreement between the two mesh schemes indicates that the numerical solution has converged to a value that accurately represents the sound pressure at the specified point. This convergence is a critical aspect of numerical methods, ensuring that the solution is reliable and independent of the discretization. The consistency of results from different mesh densities also serves as an indirect validation of the FMBEM approach. It demonstrates that the method is capable of providing accurate results even with relatively coarse meshes, which is particularly beneficial for complex geometries and larger models where fine meshing becomes impractical. Based on this, we will select the element size configuration shown in Figure 2A for the cylindrical shell model with radius a = 0.5 m and length l ≥ 0.4 m in the subsequent analysis.

Consider the model presented in Figure 1. The cylindrical shell model for analytical analysis is of infinite length, in which only a segment of length l vibrates. This infinite-length cylindrical shell is considered as the boundary of the exterior acoustic domain, meaning that the acoustic analysis does not take into account the interior region of the cylindrical shell. However, during the 3-D numerical analysis, the boundary of the exterior acoustic domain is not a cylindrical shell of infinite length but a shell of length l with closed ends, which causes the space outside the l -length segment in the infinitely long cylindrical shell to be involved in the acoustic analysis. Thus, computing errors are introduced. For the cylindrical shell model with a radius of 0.5 m and a length of 0.4 m, the analytical and 3-D numerical solutions for the sound pressure at ( 2 / 2 , 0.0, π / 4 ) are provided in Figure 4. The two curves do not fit well, as we have predicted. This discrepancy is expected and can be attributed to several factors. The analytical model assumes an infinite length for the cylindrical shell, which simplifies the boundary conditions but deviates from the actual finite length used in the numerical model. The numerical model includes the effects of the shell’s finite length, which introduces additional complexities not fully captured by the analytical model’s assumptions. The discrepancy highlights the limitations of analytical models when they are based on simplified assumptions that may not fully represent the physical system’s complexity. The comparison underscores the importance of validating analytical models with numerical simulations, especially when the models are used for design and optimization in engineering applications. The results suggest that the model configuration, particularly the length of the cylindrical shell, significantly affects the accuracy of the sound pressure predictions. This insight is crucial for selecting appropriate model parameters in both theoretical and numerical studies. The choice between analytical and numerical methods should be guided by the specific requirements of the problem, including the need for accuracy, computational resources, and the complexity of the model.

To reduce the computational errors resulting from the axial external space of the cylindrical shell, we attempt to lengthen the shell model with closed ends in the 3-D numerical analysis. For models of various lengths, the numerical results are compared with the analytical solution to determine an appropriate model configuration. Figure 5 gives the sound pressure at ( 2 / 2 , 0.0, π / 4 ) for cylindrical shell models of different lengths. The circumferential vibration order n = 0. Figures 5A–F show that the 3-D numerical result agrees well with the analytical solution when the length l ≥ 4 m. The analytical model assumes an infinite length for the cylindrical shell, which simplifies the boundary conditions but deviates from the actual finite length used in the numerical model. The numerical model includes the effects of the shell’s finite length, which introduces additional complexities not fully captured by the analytical model’s assumptions. The discrepancy between the analytical and numerical solutions for shorter shells highlights the limitations of analytical models when they are based on simplified assumptions that may not fully represent the physical system’s complexity. The comparison underscores the importance of validating analytical models with numerical simulations, especially when the models are used for design and optimization in engineering applications.

The data in Figures 5C–F, 6 give the relative error of the 3-D numerical solution to the analytical result. Similar characteristics are displayed by the four curves, namely, a rising relative error between the analytical and numerical solutions with increasing frequency, indicating that higher frequencies are more challenging to model accurately. This could be due to the higher spatial and temporal resolution required for capturing the dynamics of higher-frequency waves. For the cylindrical shell model of length l = 4 m, the relative error range is (−5%, 1%); for the model of l = 5 m, it is (−4%, 1%); and for the models of l = 9 m and l = 10 m, it is about (−3%, −1%). The relative error decreases, and the computational accuracy generally improves as the model lengthens. The present shell model has a radius of 0.5 m. The model length is 10 times the radius when l = 5 m and 20 times the radius when l = 10 m. The results suggest that the length of the model plays a critical role in the accuracy of numerical simulations. The analytical model assumes an infinite length for simplicity, but the numerical model must account for the finite length, which introduces additional complexities. The findings indicate that a model length of at least 10 times the radius is necessary for accurate uncoupling analysis. Although increasing the model length improves accuracy, it also increases the computational cost in terms of memory usage and processing time. Therefore, a balance must be struck between accuracy and computational efficiency, which is achieved when the model length is approximately 10 times the radius. In our future work, acceleration algorithms may be employed to speed up the calculation.

Consider the cylindrical shell model shown in Figure 1. The cylindrical shell is closed at both ends for the 3-D numerical analysis. The radius a = 0.5 m. The length l = 5 m or l = 10 m for the 3-D numerical analysis. Consider Equation 1. Let γ m = 1, the amplitude of the vibration W m = 2.24 × 1 0 − 8 m, and the frequency f = ω / 2 π = [100, 1,000] Hz. The analytical and 3-D numerical solutions for the sound pressure at ( 2 / 2 , 0.0, π / 4 ) for different circumferential vibration orders n are given in Figure 7. Figures 7A–F show that the 3-D numerical result agrees well with the analytical solution when the circumferential vibration order n = 1, 3, or 4. The data in Figures 7A, 8 provide, as an example, the relative error between the 3-D numerical result and the analytical solution.

Figure 7 shows that for n = 2, 6, and 10, the analytical solution goes to zero, and the numerical solution is less than 3.0E−12, which might be thought of as tending toward zero in numerical analysis. Note that the sound pressure at ( 2 / 2 , 0.0, π / 4 ) is calculated in Figure 7. Consider Equation 19. When θ = π / 4 and n = 2, 6, or 10, cos ( n θ ) = 0, meaning that the sound pressure p ( r , z , θ ) = 0 in Equation 19, which coincides with Figures 7B, E, F.

The agreement between analytical and numerical results for certain modes validates the theoretical model’s assumptions and the mathematical formulations that describe the sound pressure distribution. This is crucial for the reliability of predictions made using these models in practical applications. The results highlight the importance of considering different vibration modes in the analysis. The theoretical model predicts different behaviors for different modes, and the numerical simulations confirm these predictions, emphasizing the need for accurate modeling of vibration modes in the design and analysis of electric motors. The findings provide valuable insights for the design and optimization of electric motors, particularly in controlling the vibration modes that contribute to noise. Understanding which modes contribute to the sound pressure can guide the design of motor components to minimize noise emissions. The choice between analytical and numerical methods should be guided by the specific requirements of the problem, including the need for accuracy, computational resources, and the complexity of the model. The agreement between the two methods for certain modes suggests that, under these conditions, either method could be reliably used.

For the closed cylindrical shell model with a length of 5 m, Figures 8A–H provide the numerical sound pressure on the structure surface for various circumferential vibration orders, in which the frequency f = 500 Hz. The black line in the figure represents the section for θ = π / 4 . In Equation 19, when θ = 0, cos ( n θ ) = 1, and the sound pressure p ( r , z , θ ) reaches its maximum. This is in accordance with Figure 8I. Furthermore, in Equation 19, when n = 0, cos ( n θ ) = 1, indicating that the sound pressure p ( r , z , θ ) is constant along the cylindrical shell’s circumferential direction. This finding aligns with Figure 8A. Still in Equation 19, when θ = π / 4 and n = 2, 6, or 10, cos ( n θ ) = 0, indicating that the sound pressure p ( r , z , θ ) = 0. This conclusion is consistent with Figures 8C, F, H. Remaining in Equation 19, when θ = π / 4 and n = 4 or 8, abs ( cos ( n θ ) ) = 1, indicating that the sound pressure p ( r , z , θ ) is maximized. This is in accordance with Figures 8E, G. These results demonstrate good agreement between the results of our 3-D numerical analysis and the analytical solution. The agreement between the numerical results and the analytical predictions validates the analytical model’s accuracy in predicting the sound pressure distribution on the cylindrical shell under different vibration modes. This is essential for the model’s reliability in practical applications. The results highlight the importance of considering different vibration modes in the analysis. Each mode affects the sound pressure distribution, which is crucial for noise control and motor design.

This study concludes that the 3-D numerical analysis using FMBEM aligns well with the analytical solution for the sound pressure in the exterior acoustic domain of a cylindrical shell model, particularly when the model length is sufficient relative to its radius. The investigation reveals that computational accuracy improves as the model lengthens, suggesting a model length of at least 10 times the radius for accurate uncoupling analysis. The findings provide valuable insights for the design of low-noise motors and contribute to the understanding of electromagnetic vibration and noise control in electric motors.