Vibration modelling in full-scale elastic theory of large generator stator end-winding and optimized complex boundary

Abstract

Comparative data for a high-precise finite element model of a 600 MW large turbo-generator stator end winding already exist. Then, the double-layered winding modelling and characteristic equation under the theory of composite materials are implemented in detail. In the model, the shell is used to model double-layered winding and the supporting structures are treated as ring stiffeners and stringer stiffeners. Based on the discrete element method, the equivalent model of stiffened and ribbed conical shell for end winding can be established. After that, the natural and forced vibration equations of the end winding are established by Rayleigh–Ritz method with improved Fourier series which is suitable for different complex elastic boundary conditions. The vibration modeling is innovatively extended to the optimization model of characteristic equation for stator-winding. The optimized semi-analytical solution can find a better spring stiffness configuration to simulate complex boundary conditions. The analytical solution of the modal parameters of the end winding is obtained by calculation; a new, complete derivation of the frequency response function is carefully presented; Rayleigh damping and potential energy function of excitation force are innovatively introduced, and the displacement response analysis in multi-dimension is established. The complex evolution law in frequency domain and space domain is consistent with that of finite element analysis. Then, the rationality and correctness of the established equivalent digital mechanism model are verified.

Introduction

The vibration and noise of motor have been concerned widely [1, 2]. With the continuous enlargement of single-machine capacity of large turbo-generator, the electromagnetic force generated by powerful electromagnetic field has a greater impact on the structure [3, 4]. If the structural design of the stator end winding is unreasonable, the vibration frequency and operating amplitude of the end winding will exceed the safety standards. When the amplitude of the stator end winding exceeds a certain limit, it will seriously affect the normal operation of the generator and even lead to accidents [5, 6]. The research shows that resonance may occur when the vibration frequency of the stator end winding is close to double the power frequency. It is clearly stipulated in the national standard that for turbo-generator over 200 MW, the natural frequency of the elliptical vibration mode of the stator end winding shall be moved farther from 95 to 115 Hz [7]. Besides, the standard also specifies that the maximum amplitude of the stator end winding under long-term operation in the working environment cannot exceed 250 μm. Therefore, the research of the vibration characteristics and the dynamic mechanism model of the stator end winding have great value for engineering application under the background of industrial big data and smart manufacturing.

At present, the research on the vibration characteristics of stator end winding is mainly divided into experimental modal analysis and computational modal analysis. In the research of experimental modal analysis, Liao Wei et al. [8] used hammer method to apply excitation forces at different positions, and the response signals of displacement were detected by displacement sensors, so as to obtain modal parameters such as natural frequency through analyzing these signals. Sewak et al. [9] investigated the vibration behavior of an end winding structure by means of modal analysis and obtained that the end winding structure has elliptical vibration mode.

Computational modal analysis is mainly divided into the establishment of finite element simulation model and digital mechanism model. In previous research of authors [10, 11][40], the fine finite element modeling based on the measured data in the early stage has paved the way for the digital mechanism modeling in this paper. Wang Yixuan et al. [12, 13] established a simplified finite element model of end winding based on the simplified model of main components such as stator bars and pressure plate. Drubel et al. [14] analyzed the dynamic characteristics of the stator end winding of the turbo-generator by establishing a simplified finite element model. In this model, the three-dimensional solid element is used to simulate the bar; the spring is used to simulate the bandage; the shell is used to simulate other components. Iga et al. [15] focused on modeling accurately the stator bars and bindings. The bindings have been modeled with beam elements and multipoint constraint equations to account for the bending, twisting, and shearing rigidities of the bindings, and a bar is modeled as a composite material comprising a conductor bar with equivalent elastic modulus and ground insulation.

To sum up, both experimental modal analysis and finite element modeling are based on specific experimental objects, which do not have the general applicability of digital mechanism model, and cannot understand the inherent laws of the object, such as the mathematical relationship between physical parameters and modal parameters. Therefore, the two methods have limited guiding significance for the design of large-scale structures in intelligent manufacturing and the dynamic optimization monitoring of big data.

Digital mechanism modeling can not only clearly reflect the inherent laws of the object, but also has strong universality. Because the structure of stator end winding itself is extremely complex, it is difficult to analyze the mathematical relationship of various components of the structure. Therefore, there is little research on digital mechanism modeling of stator end winding at present. According to other existing analogies of the actual structure studied, the stator end winding can be equivalent to the anisotropic truncated conical shell structure. In detailed modeling approach of composite material theory (Sect. 3.2), the stator end winding can be considered as a composite conical shell structure, which is used to model double-layer winding, and the supporting structures are treated as ring-stiffeners and stringer-stiffeners [40]. And to simulate the actual complex boundary condition, the new optimization model of characteristic equation for stator-winding has been proposed and solved by heuristic intelligent algorithms (Sect. 4.3), then new derivation of the frequency response function and Rayleigh damping and potential energy function of excitation force are innovatively constructed (Sect. 5.2). Moreover, the research on the vibration model of shell has some basic models. Liew [16] and Li Fengming [17] obtained the free vibration solutions of conical shells by using Rayleigh–Ritz method. Jin et al. [18] researched the free vibration analysis of truncated conical shells with general elastic boundary conditions by using an accurate modified Fourier-series solution. Kouchakzadeh et al. [19] investigated the natural frequencies and mode shapes of two joined cross-ply laminated conical shells based on the thin-walled shallow shell theory of Donnell type and Hamilton’s principle. Ma et al. [20] researched a free and forced vibration analysis of coupled conical-cylindrical shells with arbitrary boundary conditions by using a modified Fourier-Ritz method. In addition, some scholars have studied the vibration characteristics of the motor. Xing Zezhi et al. [21] regarded the stator core and enclosure as cylindrical shells with axial ribs attached inside and outside the shells respectively. They discreted the vibration displacement of the continuous system by Galerkin discrete and Rayleigh–Ritz method, and derived the characteristic equations of the stator under free vibration based on the three-dimensional shell theory, and obtained the natural frequencies of the stator with anisotropic material. Hu Shenglong et al. [22] regarded the assembly of the enclosure, stator and windings as a cylindrical shell and it is supported at both-end covers. They established a theoretical model of the equivalent cylindrical shell with five orthotropic material parameters, and obtained the natural frequencies and the modal shapes by considering the simple-support boundary condition and the orthotropic characteristics of the motor. Zhao Guodong et al. [23] studied the stator vibration characteristics of the motor shell structure under the boundary conditions. Artificial springs are introduced to simulate boundary conditions. Further considering the influence of teeth, winding, frame and cooling ribs, the teeth with complex shape are regarded as longitudinal ribs of infinite cross-section.

Firstly, as a contrastive study, the comparison data of accurate finite element modeling were retained [10, 11]. Secondly, considering the equivalent detailed modeling approach of composite material theory (Sect. 3.2), an anisotropic truncated conical shell model for stator end winding has been implemented [40]. The dynamic characteristics of the equivalent conical shell model are solved by Rayleigh–Ritz method with improved Fourier series, which is suitable for different complex elastic boundary conditions. Thirdly, the stator end winding is equivalent to a model of double-layered stiffened and ribbed conical shell by the discrete element method. The natural and forced vibration equations of the end winding, with better new intelligent optimization modelling compared with empirical formula for simulating constraints of actual complex boundary (Sect. 4.3), are established, and the modal parameters of the digital mechanism model are obtained. Finally, a new, complete derivation of the frequency response function is carefully presented, Rayleigh damping and potential energy function of excitation force are innovatively introduced, and the displacement response analysis in multi-dimension is established (Sect. 5.2).

The finite element simulation and modal analysis of stator end winding

The finite element model of stator end winding

The structure of stator end winding is similar to cantilever beam, that is, the straight-line segment of the stator bar is fixed in the iron core, and the involute parts depends on the supporting structures. The whole structure is mainly composed of two parts, one is basket double-layered winding that is formed from the upper and lower bars, spacer blocks and rings, and the other is outermost layer supporting structure that is composed of radial brace, L-shaped bracket, copper bracket and lead ring, as shown in Fig. 1. According to the actual connection relationship between components, the 3D finite element model of stator end winding can be established by ANSYS, as shown in Fig. 1.

[See PDF for image]

Fig. 1

Stator end winding and three-dimensional diagram

The modal analysis of stator end winding (FEM based on experimental data)

After measurement and calculation, the physical and mechanical properties of main components are shown in Table 1. The connection of each component in the finite element model can be set as shown in Table 2. Figure 1 shows the schematic diagram of the three-dimensional model of the stator-end winding (In particular, pay attention to the small-diameter end, large-diameter end). Specific details of finite element modeling are detailed in the Appendix F.

Table 1. Parameters for main components of physical and mechanical properties

Component name	Young’s modulus (MPa)	Poisson ratio	Density (kg/m3)
stator bar (xb)	Exb = 40,000	μxb = 0.300	ρxb = 4033
inter-layer ring (cjh)	Ecjh = 39,045	μcjh = 0.267	ρcjh = 1893
radial brace (si)	Esi = 24,826	μsi = 0.173	ρsi = 1842
support ring (ri)	Eri = 39,045	μri = 0.267	ρri = 1893

Table 2. Connections of the components

Contact body	Target body	Connections
Stator bar	Spacer block	Bonded and Frictional
Stator bar	Ring and support ring	Bonded
Stator bar	Box	Bonded
Radial brace	Ring	Bonded
Copper bracket	Pin	Translational
Radial brace	Pin	Revolute

Through the modal analysis by ANSYS Workbench, the mode shapes of each order of stator end winding and their corresponding natural frequencies can be obtained. The calculated results are compared with the measured results in reference [11], as shown in Fig. 2, and the relative error of the corresponding natural frequencies is within 2.5%.

[See PDF for image]

Fig. 2

Comparison of mode shapes determined by modal analysis of finite element and experimental data

The natural frequency studied is concentrated in the frequency band of 80 ~ 114 Hz, and the following settings are carried out in the harmonic response analysis of ANSYS Workbench:

1. Frequency Spacing select Linear;

2. Range Minimum is set to 75 Hz;

3. Range Maximum is set to 125 Hz;

4. Solution Intervals is set to 100, i.e., two neighboring frequencies are separated by 0.5 Hz.

Harmonic response analysis using ANSYS Workbench requires the setting of a Modal Frequency Range, which refers to the frequency range of the structure involved in the vibration modal analysis.

The overview of the equivalent conical-shell model of stator end winding

The overview of double-layered stiffened and ribbed conical-shell model

The stator end winding is mostly a non-continuous structure, which is composed of double-layered winding structure and supporting structures, and is in the shape of truncated cone. The structure of the stator end winding of a large turbine-generator is not a continuous medium component by nature. The test shows that its vibration characteristics reflect good integrity and have the vibration characteristics of a complete continuous medium [24]. And it is difficult to establish its vibration model in the forward direction. Although the vibration mode of stator end winding is regular and the components are closely connected in the actual structure, it can be regarded as an integrated structure. The equivalent of simple isotropic material or anisotropic material will make it difficult to characterize the physical properties of the main components. Therefore, considering it as composite material, it can fully reflect the influence of the physical parameters of the stator bars and rings on the structure.

The structure of stator-end winding is equivalent to a continuous shell, which can realize the forward digital modeling of the structure on the basis of retaining the vibration characteristics. In addition, in the existing research, some scholars equivalent the structure to continuous anisotropic cylindrical shell [21, 22–23] (stator core and enclosure) or partial preliminary study of isotropic conical shell [13] and so on. So, according to the actual structure, the main part of the binding double-layered winding, the stator bar, is equivalent to fiber. And the binding part, the interlayer ring, is equivalent to the matrix filled in the gap. The structure is equivalent to the continuous truncated laminated composite-conical shell. For the supporting structures, in the existing research, some scholars equivalent it to beam structure, thin plate structure and stiffened structure respectively [13]. Therefore, according to the characteristics, that the radial brace distributed plays the main supporting role, and the nasal-tip support-ring plays the secondary supporting role, adopts the stiffened and ribbed method to model it.

Differing from [21, 22], this paper aims at the stator-end winding of 600 MW large turbo-generator, which has stronger discontinuity and larger installed volume of the motor (In order to better intuitively understand the larger installed volume of the motor, you can refer to the size scale in Fig. 2, which undoubtedly increases the difficulty of mechanism modeling of vibration characteristics.). Equivalent conical shell model is solved by Rayleigh–Ritz method with improved Fourier series, which is suitable for different complex elastic boundary conditions (refer to Sect. 4.3). Double-layered stiffened and ribbed conical shell modeled by the discrete element method has higher solving precision (refer to Sects. 4.1 and 5.1). In conclusion, based on the structural characteristics of the stator end winding and ignoring the influence of other components, the whole structure can be equivalent to the stiffened and ribbed laminated composite conical-shell model, as shown in Fig. 3(a), and boundary constraints are shown in Fig. 3(b).

[See PDF for image]

Fig. 3

Stator end winding structure and equivalent structure

Support structure and boundary constraints: the linear section of the stator winding has a part firmly embedded in the core groove, and the other part extends out of the core groove. Therefore, the section of the protruding part adopts a fixed constraint (Fixed Support B in Fig. 3(b)). The end face of the copper bracket is firmly fixed with the copper shield, so the end face of the copper bracket is fixed (Fixed Support A in Fig. 3(b)). The radial support is connected with the copper bracket fixed on the copper shield through a sliding pin, which can slide along the chute of the copper bracket to form an elastic structure along the axis [11].

The mechanism modeling is to establish the mathematical relationship of physical parameters, and the accuracy of the mathematical relationship of these leading factors depends on reasonable parameter. Therefore, in order to ensure the reliability of the comparison between the analytical solution and the finite element results, the main parameters in the formula are consistent with the structural parameters of the finite element model (refer to Sect. 2).

A laminated composite material composed of two layers

There are 84 stator bars in the stator end winding, and the stator bars at the same layer are separated by spacer blocks. There are three interlayer rings between the upper and lower bar, and the lower bar is connected with the radial brace through three outer-layer rings. Thus, the whole structure of stator end winding is formed.

The stator end winding (double-layer basket structure) is composed of the upper and lower layers of winding. The single-layer is equivalent to the lamellar series composite material (See Appendix A for effectiveness). According to the theory of composite materials [25], the following equations are obtained for the elastic modulus E₁₍₂₎, Poisson's ratio, shear modulus and density of the single-layer winding in the direction of coordinates 1–2:

1

$$\left\{\begin{array}{c}{E}_1=K_1(E_{f1}{V}_f+E_{m1}{V}_m)\\ E_2=K_2{E}_{f2}{E}_{m2}\left./\phantom{K_2{E}_{f2}{E}_{m2}\left(E_{m2}{V}_f+E_{f2}{V}_m\right)}\right)\left(E_{m2}{V}_f+E_{f2}{V}_m\right)\\ G_{12}=G_{f12}{G}_{m12}\left./\phantom{G_{f12}{G}_{m12}(G_{m12}{V}_f+G_{f12}{V}_m)}\right)(G_{m12}{V}_f+G_{f12}{V}_m)\\ {\mu }_{12}={\mu }_{21}={\mu }_{f12}{V}_f+{\mu }_{m12}{V}_m\\ \tilde{\rho }={\rho }_f{V}_f+{\rho }_m{V}_m\\ \end{array}\right)$$where E_f1(2), E_m1(2), G_f12, G_m12, and μ₁₂₍₂₁₎ denote the modulus of elasticity, shear modulus, and Poisson's ratio of the fiber to the substrate in the direction of coordinates 1–2, respectively; V_f and V_m are the volume ratio of the shells occupied by the fiber and the substrate, respectively; and K₁ and K₂ are the correction coefficients for the 1 and 2 directions. In this paper, K₁ and K₂ are set to be 1. The single-layer is equivalent to the lamellar series composite material (See Appendix A.4 for more details).

The establishment and solution of vibration equation of conical shell based on Rayleigh–Ritz method

In this section, based on references [13, 18], considering the complex boundary conditions of stator end winding in the actual working environment, the Rayleigh–Ritz method with improved Fourier series is used to solve the dynamic characteristics of the equivalent conical shell model.

The establishment of double-layer equivalent conical shell model

The double-layer winding can be regarded as the superposition of single-layer winding (refer to Appendix A and Sect. 3.2), and constitutive equation is established by using the theory of laminated composites [40]. As shown in Fig. 3, the orthogonal curvilinear coordinate (x,θ,z) of the equivalent conical shell model is depicted. Specifically, x, θ and z are coordinates along the generatrix direction of the conical shell surface, the tangential direction of radius of the conical shell and the vertical direction of the conical shell surface. The displacements of a point on the middle surface are indicated by u, v and w along x, θ and z directions, respectively. R₁ and R₂ are the radius of the conical shell at its iron core and lead ring end, respectively. α is the semi-vertex cone angle of the conical shell. The generatrix length is L and the thickness of shell is H. For the double-layer winding structure, a stacked anisotropic conical shell model consisting of a superimposed combination of two layer-winding structures, the three-dimensional model is shown in Fig. 4. It is equivalent to a laminated composite continuously truncated conical shell structure, and the geometry is schematically shown in Fig. 5. The geometrical parameters are taken from the finite element model in subsection 2.1.

[See PDF for image]

Fig. 4

The three-dimensional model of double-layer winding structure

[See PDF for image]

Fig. 5

Geometric structure diagram of truncated conical shell

According to the Reissner’s shell theory [26, 27–28], the geometric equations of the mid-plane strain and curvature of the shell can be obtained, as shown in Eq. (2).

2

$$\left\{\begin{array}{c}{\epsilon }_x^0=\frac{\partial u}{\partial x},{\epsilon }_{\theta }^0=\frac{1}{R(x)}\frac{\partial v}{\partial \theta }+\frac{usin\alpha +wcos\alpha }{R(x)}\\ {\gamma }_{x\theta }^0=\frac{1}{R(x)}\frac{\partial u}{\partial \theta }+\frac{\partial v}{\partial x}-\frac{vsin\alpha }{R(x)}\\ k_x^0=-\frac{{\partial }^{2}w}{\partial x^2},k_{\theta }^0=-\frac{1}{R^2(x)}\frac{{\partial }^{2}w}{\partial {\theta }^2}+\frac{cos\alpha }{R^2(x)}\frac{\partial v}{\partial \theta }-\frac{sin\alpha }{R(x)}\frac{\partial w}{\partial x}\\ {\tau }_{x\theta }^0=2\left(-\frac{1}{R(x)}\frac{{\partial }^{2}w}{\partial x\partial \theta }+\frac{sin\alpha }{R^2(x)}\frac{\partial w}{\partial \theta }+\frac{cos\alpha }{R(x)}\frac{\partial v}{\partial x}-\frac{vsin\alpha cos\alpha }{R^2(x)}\right)\\ \end{array}\right)$$where R(x) = R₁ + xsinα.

For the double-layered winding structure, the constitutive equation of the truncated conical shell can be written as (See Appendix B for more details. Equation (3) is also equation Eq. (61)):

3

$$\left[\begin{array}{c}{N}_x\\ N_{\theta }\\ N_{x\theta }\\ M_x\\ M_{\theta }\\ M_{x\theta }\end{array}\right]=\left\{\begin{array}{cccccc}{A}_{11}& A_{12}& A_{16}& B_{11}& B_{12}& B_{16}\\ A_{12}& A_{22}& A_{26}& B_{12}& B_{22}& B_{26}\\ A_{16}& A_{26}& A_{66}& B_{16}& B_{26}& B_{66}\\ B_{11}& B_{12}& B_{16}& D_{11}& D_{12}& D_{16}\\ B_{12}& B_{22}& B_{26}& D_{12}& D_{22}& D_{26}\\ B_{16}& B_{26}& B_{66}& D_{16}& D_{26}& D_{66}\end{array}\right\}\left\{\begin{array}{c}{\epsilon }_x^0\\ {\epsilon }_{\theta }^0\\ \gamma _{_{x\theta }^{}}^0\\ k_x^0\\ k_{\theta }^0\\ {\tau }_{x\theta }^0\\ \end{array}\right\}$$

The establishment of vibration mode function

For the traditional solution method of conical shell, the discontinuity of displacement function on the shell boundary is often ignored, so it is difficult to obtain the accurate solution of conical shell model. In order to ensure the continuity of the displacement function under any complex boundary conditions, the displacement function can be expanded into Fourier series and an auxiliary function can be introduced, which can not only ensure the continuous value of the whole integral domain L, but also greatly improve the convergence speed. The geometrical parameters are taken from the finite element model in subsection 2.1.

According to the vibration characteristics of the conical shell along the generatrix direction of the conical shell surface and the tangential direction of radius of the conical shell [30], the vibration displacements u, v and w can be written as

4

$$\left\{\begin{array}{c}u=U_n(x)cos(n\theta )e^{j\omega t}\\ v=V_n(x)sin(n\theta )e^{j\omega t}\\ w=W_n(x)cos(n\theta )e^{j\omega t}\\ \end{array}\right)$$where the vibration mode function U_n(x), V_n(x) and W_n(x) can be expanded into Fourier series [31, 32] and can be written as

5

$$\left\{\begin{array}{c}{U}_n(x)=\sum _{m=0}^{\infty }{A}_{m}cos\frac{m\pi }{L}x+a_1{f}_1(x)+a_2{f}_2(x)\\ V_n(x)=\sum _{m=0}^{\infty }{B}_{m}cos\frac{m\pi }{L}x+b_1{f}_1(x)+b_2{f}_2(x)\\ W_n(x)=\sum _{m=0}^{\infty }{C}_{m}cos\frac{m\pi }{L}x+c_1{f}_1(x)+c_2{f}_2(x)+c_3{f}_3(x)+c_4{f}_4(x)\\ \end{array}\right)$$where A_m, B_m and C_m are the amplitudes of vibration mode functions in each direction, and f₁(x), f₂(x), f₃(x) and f₄(x) are auxiliary functions. Refer to Appendix C for more details of the auxiliary function.

The establishment of energy function of conical shell

The kinetic energy of the shell can be written as

6

$$T_{\epsilon }=\frac{1}{2}\rho H{\int }_0^L{\int }_0^{2\pi }\left[{\left(\frac{\partial u}{\partial t}\right)}^2+{\left(\frac{\partial v}{\partial t}\right)}^2+{\left(\frac{\partial w}{\partial t}\right)}^2\right]R(x)d\theta dx$$

According to the principle of virtual work, for the equivalent conical shell model, the strain energy of the shell can be written as

7

$$\begin{array}{c}{U}_{\epsilon }=\frac{1}{2}{\int }_0^L{\int }_0^{2\pi }\left[N_x{\epsilon }_x^0+N_{\theta }\epsilon _{_{\theta }^{}}^0+N_{x\theta }\gamma _{_{x\theta }^{}}^0\right]R(x)d\theta dx\\ +\frac{1}{2}{\int }_0^L{\int }_0^{2\pi }\left[M_{x}k_{_x^{}}^0+M_{\theta }k_{_{\theta }^{}}^0+M_{x\theta }\tau _{_{x\theta }^{}}^0\right]R(x)d\theta dx\\ \end{array}$$

The modelling of complex boundary: In addition to the strain energy of the shell, in order to equivalent the elastic boundary constraints (Consider the support structure and boundary constraints shown in Fig. 3(b) in Sect. 3.1.), the boundary potential energy should also be considered. Due to the discontinuity between bars and bars, it is not possible to obtain suitable results using any classical boundary conditions.

To establish an accurate vibration model of conical shell, two sets of springs are introduced at both ends of the conical shell. Specifically, one set has 4 springs. k_u0, k_v0, k_w0 and k_θ0 denote the springs of iron core end, and k_u1, k_v1, k_w1 and k_θ1 denote the springs of lead ring end. k_u, k_v and k_w represent the displacement-constraint spring respectively along the u, v and w directions, and k_θ represent the rotation-constraint spring along the rotational θ direction. By setting the stiffness factor of the spring to zero or infinity, the classic boundary as shown in Table 3 is obtained.

Table 3. Comparison of the relationship of parameters between classical boundary and elastic boundary (N/m)

Free	Simple	Clamped
ku0 = ku1 = 0	ku0 = ku1 = 0	ku0 = ku1 = 1013
kv0 = kv1 = 0	kv0 = kv1 = 1013	kv0 = kv1 = 1013
kw0 = kw1 = 0	kw0 = kw1 = 1013	kw0 = kw1 = 1013
kθ0 = kθ1 = 0	kθ0 = kθ1 = 0	kθ0 = kθ1 = 1013

Similarly, if the stiffness values of eight boundary constraint springs are set, elastic boundary conditions can be simulated. The stator end winding is bound by multiple stator bars, so the two ends are discontinuous. Moreover, the boundary constraint can only restrict the stator bar part, but not the gap part (that is, it belongs to the complex elastic boundary condition), as shown in Fig. 6.

[See PDF for image]

Fig. 6

Diagram of boundary constraints

For the same stator end winding model, when the displacement boundary conditions at both ends change, the characteristics and energy distribution of vibration will also change, and the natural frequency obtained will also be different. Therefore, by establishing an elastic boundary, a vibration model under arbitrarily complex boundary conditions can be established without changing the parameters of the structure itself (refer to Eqs. (8) and (9)). As shown in Fig. 6 the boundary constraint cannot act on the gap part. And the larger the gap area is, the smaller the constraint effect is. Therefore, elastic boundary is required for equivalent boundary constraints (When Clamped are used, the calculation results are too large. Refer to Fig. 14(a)–(b).).

According to the data analysis in Fig. 7, the natural frequency of the structure changes obviously and increases with the increase of stiffness in a certain stiffness interval (sensitive interval), and changes weakly outside this interval. Because the actual boundary conditions are not continuous, it is necessary to consider the influence of the actual constraint area and the gap area. As the larger the proportion of constrained area, the closer to the clamped support constraint, and then the natural frequency will also increase, which is related to the increase of natural frequency as the spring stiffness increases in the sensitive interval.

[See PDF for image]

Fig. 7

The curve of natural frequency and spring stiffness value of cone shell structure under elastic boundary

Therefore, two approaches for constructing complex elastic boundaries are proposed. One is to use the proportion of constrained area and the sensitive interval to construct the elastic boundary formula (forward calculation of virtual spring stiffness), and the other is an elastic boundary optimization based on heuristic algorithm (reverse construction of equivalent boundary). The first method may be difficult to meet engineering requirements accurately, but is used in qualitative analysis in early development. The second method uses the intelligent optimization algorithm to find the exact solution to meet the engineering requirements by refining the optimization objective function, so as to obtain a more perfect elastic boundary. Further, two research methods have been developed.

The empirical formula of complex boundary: The empirical formula for constructing the elastic boundary is determined by the following formula,

8

$$k_i={10}^{\beta },\beta =log{k}_i^{min}+\frac{S_x}{S_x+S_k}\left(log{k}_i^{max}-log{k}_i^{min}\right)$$where i takes the vibration displacement direction u, v, w and the rotational direction θ, $k_i^{}$ is the spring stiffness in a particular direction, $k_i^{min}$ is the lower limit of the sensitive interval, and $k_i^{max}$ is the upper limit of the sensitive interval (in Fig. 8, the sensitive interval of $k_{u0}^{}$ is [1e7, 1e12]); S_x and S_k denote the cross sectional area of the stator bar and the gap part on the constraint surface respectively.

[See PDF for image]

Fig. 8

The curve of natural frequency and spring stiffness value under elastic boundary

In order to establish the influence characteristics of the boundary conditions in all directions on the vibration of the laminated shell, the stiffness value of one of the springs is gradually increased from 0 to 10¹³ after setting the boundary constraint condition. Among them, the change curve of the natural frequency of the truncated conical shell structure with boundary conditions (Clamped and Clamped, C–C: change the spring stiffness of large-diameter end) is shown in Fig. 7, and the change curve of the natural frequency of the double-layer winding structure with boundary conditions (Clamped and Free, C-F: change the spring stiffness of small-diameter end) is shown in Fig. 8. The physical parameters are as follows: density ρ = 7800 kg/m³, Poisson's ratio μ = 0.3, half cone angle α = 30°, Young's modulus E = 210 GPa, small-diameter end radius R₁ = 10 m, large-diameter end radius R₂ = 20 m, and shell thickness H = 0.2 m.

For the conical shell structure, among the four sets of boundary spring constraints, Fig. 7 illustrates that the circumferential and radial springs exert the most important impact on the natural frequency for a given wave-number, whereas the influence of the coiled spring remains relatively minor when the circumferential wave-number is below 5. However, with increasing circumferential wave-number, the effect of the radial and coiled springs on the natural frequency gradually becomes more pronounced. When the circumferential wave-number reaches 10, the influence of the radial and coiled springs surpasses that of the circumferential and axial springs, respectively. Additionally, for stiffness values below 10⁵ N/m and beyond 10¹² N/m, the natural frequency remains nearly unchanged, indicating that a stiffness bigger than 10¹² N/m can be regarded as a fixed constraint. Consequently, a stiffness value of 10¹³ N/m is justifiably selected as a fixed constraint. Furthermore, Fig. 7 reveals that the structure exhibits a distinct sensitivity interval concerning the natural frequency for each spring.

It can be seen from Fig. 8 that for the double-layer winding structure, among the four sets of boundary spring constraints, the influence of the spring in u direction is greater than that of the spring in other directions at n = 2. When n = 3, the influence relationship is ordered from large to small as u, w, θ and v, respectively. Further, from Fig. 8, it can be found that the sensitive interval for the natural frequency of the structure is determined relative to the change of the spring stiffness in a certain direction, [1e7, 1e12] for the u-direction, [1e8, 1e11] for the v-direction, [1e7, 1e11] for the w-direction, and [1e6, 1e9] for the θ-direction (to calculate the equivalent boundary of the winding, interval ranges of four spring need to be brought into Eq. (8).

From the parametric simulation experiments, it can be seen that the four boundary springs will have a significant impact on the frequency in the specific range of 10⁶ ~ 10¹¹. Therefore, when using the elastic boundary-constraint spring, the elastic potential energy of the boundary-constraint spring stored at both ends of the conical shell model can be written as

9

$$\begin{array}{c}{U}_b=\frac{1}{2}{\int }_0^{2\pi }{\left[k_{u0}{u}^2+k_{v0}{v}^2+k_{w0}{w}^2+k_{\theta 0}{\left(\frac{\partial w}{\partial x}\right)}^2\right]}_{x=0}{R}_{1}d\theta \\ +\frac{1}{2}{\int }_0^{2\pi }{\left[k_{u1}{u}^2+k_{v1}{v}^2+k_{w1}{w}^2+k_{\theta 1}{\left(\frac{\partial w}{\partial x}\right)}^2\right]}_{x=L}{R}_{2}d\theta \\ \end{array}$$

The total kinetic energy of the shell can be written as

10

$$T_k=T_{\epsilon }$$

The total strain energy of the shell can be written as

11

$$U_k=U_{\epsilon }+U_b$$

The evolutionary computational optimization of complex boundary

Constructing the optimization model of complex elastic boundary based on heuristic algorithm is as follows. It is necessary to establish the optimization objective function according to Eqs. (10) and (11). The natural frequency error and mode shape error are considered at the same time, and the dynamic penalty function is used to construct the objective function. First, the Lagrange function of the vibration system for double-layer winding is as follows:

12

$$L=U_{\epsilon }+U_b-T_{\epsilon }$$

The characteristic equation of the system obtained by the derivation of Ritz method is as follows (refer to the derivation of Eq. (35)):

13

$$\begin{array}{c}\left(\left[{\mathbf{K}}_{\epsilon }\right]+\left[{\mathbf{K}}_{b0}\right]+\left[{\mathbf{K}}_{b1}\right]-{\omega }^2\left[\mathbf{M}\right]\right)\mathbf{X}=0\\ \mathbf{X}={[{\mathbf{u}}^T{,\mathbf{v}}^T{,\mathbf{w}}^T]}^T,\left\{\begin{array}{c}\mathbf{u}={\left[A_0,A_1,\dots ,A_m,a_1,a_2\right]}^T\\ \mathbf{v}={\left[B_0,B_1,\dots ,B_m,b_1,b_2\right]}^T\\ \mathbf{w}={\left[C_0,C_1,\dots ,C_m,c_1,c_2,c_3,c_4\right]}^T\\ \end{array}\right)\\ \end{array}$$where [K_ε] represents the stiffness matrix associated with strain energy, while [K_b0] and [K_b1] correspond to the stiffness matrices of the boundary spring's potential energy at the small-diameter and large-diameter ends, respectively. Given that the large-diameter end of the studied vibration model is a free boundary, it follows that [K_b1] = 0. Furthermore, [K_b0] is expanded and rewritten as the product of each spring's stiffness and its corresponding stiffness matrix. To enhance computational efficiency, stiffness variables are separately extracted, and the coefficient matrix associated with these variables is precomputed. This approach avoids redundant integral and differential operations, thereby significantly improving optimization performance.

14

$$\left[{\mathbf{K}}_{b0}\right]=k_{u0}\left[{\mathbf{K}}_{u0}\right]+k_{v0}\left[{\mathbf{K}}_{v0}\right]+k_{w0}\left[{\mathbf{K}}_{w0}\right]+k_{\theta 0}\left[{\mathbf{K}}_{\theta 0}\right]$$

Since the physical and geometric parameters of the same type winding structure have been determined, the natural frequency and vibration mode are only related to the spring stiffness. The optimized model has general applicability to a wide range of complex elastic boundaries. In particular, [K_b1] = 0, then the optimization variable is defined as:

15

$$\begin{array}{c}X={\left[x_1,x_2,x_3,x_4\right]}^T={\left[k_{u0},k_{v0},k_{w0},k_{\theta 0}\right]}^T\\ x_{imin}\le x_i\le x_{imax},i=1,2,...,4\\ \end{array}$$

Given that the primary focus of the winding structure is on the vibrations associated with the elliptical and triplet modes, the objective function can be formulated as:

16

$$f\left(\mathit{X}\right)=\left\{\begin{array}{c}\begin{array}{c}\end{array}\sum _{i=2}^3\left|{\omega }^{\left(1,i\right)}\left./\phantom{{\omega }^{\left(1,i\right)}\left(2\pi \overline{f}_{_{(1,i)}^{}}^{}\right)}\right)\left(2\pi \overline{f}_{_{(1,i)}^{}}^{}\right)-1\right|\text{if}\mathrm{\Delta }\le \epsilon \\ Z\sum _{i=2}^3\left|{\omega }^{\left(1,i\right)}\left./\phantom{{\omega }^{\left(1,i\right)}\left(2\pi \overline{f}_{_{(1,i)}^{}}^{}\right)}\right)\left(2\pi \overline{f}_{_{(1,i)}^{}}^{}\right)-1\right|\text{if}\mathrm{\Delta }>\epsilon \\ \mathrm{\Delta }=\left|{\omega }^{\left(1,i\right)}\left./\phantom{{\omega }^{\left(1,i\right)}\left(2\pi \overline{f}_{_{(1,i)}^{}}^{}\right)}\right)\left(2\pi \overline{f}_{_{(1,i)}^{}}^{}\right)-1\right|(i=2,3)\\ \end{array}\right)$$where ω^(1,i) denotes the characteristic root obtained by solving Eq. (16) for a given wave number, while $\overline{f}_{_{(1,i)}^{}}^{}$ represents the corresponding target frequency for that wave number, determined using the modal frequency from finite element analysis. The deviation threshold ε is set to 0.05. Moreover, Z is a predefined positive penalty factor, assigned a value of 1000.

To enhance the accuracy of the characteristic equation optimization model, four sets of boundary constraints are incorporated to regulate the vibration behavior at the small-diameter end, thereby improving its consistency with the finite element model. To effectively constrain the vibration at the small-diameter end, the vibration vectors of specific nodes are required to approach zero across all global nodes. Additionally, the maximum displacement at the small-diameter end is restricted within a predefined threshold. Based on these conditions, the optimization model is formulated as follows:

17

$$\left\{\begin{array}{c}minf\left(\mathit{X}\right)\text{, s}\text{.t}.x_{imin}\le x_i\le x_{imax},i=1,2,...,4\\ \left|{\varphi }_{\mathit{uj}}\left(x_j,{\theta }_j,\mathit{X}\right)\right|<\mathrm{\Delta }{S}_u,{\varphi }_{\mathit{uj}}\left(x_j,{\theta }_j,\mathit{X}\right)=\frac{u_j\left(x_j,{\theta }_j,\mathit{X}\right)}{max(u_{\mathit{all}},v_{\mathit{all}},w_{\mathit{all}})}\\ \left|{\varphi }_{\mathit{vj}}\left(x_j,{\theta }_j,\mathit{X}\right)\right|<\mathrm{\Delta }{S}_v,{\varphi }_{\mathit{vj}}\left(x_j,{\theta }_j,\mathit{X}\right)=\frac{v_j\left(x_j,{\theta }_j,\mathit{X}\right)}{max(u_{\mathit{all}},v_{\mathit{all}},w_{\mathit{all}})}\\ \left|{\varphi }_{\mathit{wj}}\left(x_j,{\theta }_j,\mathit{X}\right)\right|<\mathrm{\Delta }{S}_w,{\varphi }_{\mathit{wj}}\left(x_j,{\theta }_j,\mathit{X}\right)=\frac{w_j\left(x_j,{\theta }_j,\mathit{X}\right)}{max(u_{\mathit{all}},v_{\mathit{all}},w_{\mathit{all}})}\\ max(\sqrt{u_j^2+v_j^2+w_j^2})<\mathrm{\Delta }{S}_L\\ \end{array}\right)$$

In the stiffness configuration X, the vibration vector of node j in the u, v, and w direction is represented as ${\varphi }_{u{(v,w)}_j}(x_j,{\theta }_j,\mathit{X})$.The maximum absolute displacement across all global nodes (all nodes situated on the neutral plane of the conical shell) in each direction is given by $max(u_{\mathit{all}},v_{\mathit{all}},w_{\mathit{all}})$. Furthermore, the displacement of node j in u(v,w) direction is denoted by $u_{{(v,w)}_j}(x_j,{\theta }_j,X)$.

In the large-diameter region (x = L), the maximum node displacement is observed. For n = 2, due to symmetry, the maximum displacement in each direction occurs at (x,θ) = (L,0) or (L,pi/2). When n = 3, it appears at (x,θ) = (L,0) or (L,pi/3). Consequently, it is unnecessary to compute the displacement of every global node; instead, evaluating three nodes (x,θ) = (L,0), (L,pi/3) and (L,pi/2) suffices to determine $max(u_{\mathit{all}},v_{\mathit{all}},w_{\mathit{all}})$. This approach significantly enhances iteration efficiency. The parameters $\mathrm{\Delta }{S}_u,\mathrm{\Delta }{S}_v,\mathrm{\Delta }{S}_w$ and $\mathrm{\Delta }{S}_L$ govern the vibration mode. When the parameters are decreased, namely, when the vibration-mode vector is reduced (leading to an enhancement in the mode matching degree), the vibration mode at the minor end diminishes, resulting in an increased deviation in natural frequency. Conversely, as the mode matching degree declines, the deviation in natural frequency is correspondingly mitigated. Further details regarding the optimization process and parameter configurations of the optimization model are provided in Sect. 6.1.

Owing to the incorporation of inequality constraints, direct optimization through intelligent optimization algorithms is not feasible. To address this, the dynamic penalty function method can be employed to convert the imposed inequality constraints into penalty terms, thereby reformulating the constrained optimization problem into an unconstrained one [33]:

18

$$\begin{array}{c}\left\{\begin{array}{c}{g}_{1j}(\mathit{X})=\left|{\varphi }_{\mathit{uj}}\left(x_j,{\theta }_j,\mathit{X}\right)\right|-\mathrm{\Delta }{S}_u,g_{2j}(\mathit{X})=\left|{\varphi }_{\mathit{vj}}\left(x_j,{\theta }_j,\mathit{X}\right)\right|-\mathrm{\Delta }{S}_v\\ g_{3j}(\mathit{X})=\left|{\varphi }_{\mathit{wj}}\left(x_j,{\theta }_j,\mathit{X}\right)\right|-\mathrm{\Delta }{S}_w,g_4(\mathit{X})=max(\sqrt{u_j^2+v_j^2+w_j^2})-\mathrm{\Delta }{S}_L\\ \end{array}\right)\\ F(\mathit{X},\sigma )=f(\mathit{X})+\sigma \left(\sum _{j=1}^P\sum _{i=1}^{3}max\left(0,g_{\mathit{ij}}(\mathit{X})\right)+max\left(0,g_4(\mathit{X})\right)\right),\sigma =T\sqrt{T}\\ \end{array}$$g(X) represents the penalty term, P denotes the number of nodes, σ is the dynamic penalty factor, and T corresponds to the number of iterations. By dynamically adjusting the penalty factor, a balance between exploration and exploitation within the population can be achieved. Specifically, a smaller penalty factor in the early stages facilitates population exploration, whereas a progressively increasing penalty factor in the later stages enhances exploitation and accelerates convergence.

In conclusion, the formulated optimization model is expressed as follows:

19

$$\left\{\begin{array}{c}minF(X,\sigma )\hfill \\ \mathrm{s}.\mathrm{t}.x_{imin}\le x_i\le x_{imax},i=1,2,...,4\hfill \end{array}\right)$$

Three heuristic algorithms (LDIW-PSO [34], RIW-PSO [35], IA [36]) are used to optimize the model. The optimization results are shown in Sect. 6.1.

The establishment of energy function of stiffened and ribbed

The stator end winding is equivalent to the double-layer stiffened and ribbed conical shell model; the double-layer winding structure is mainly equivalent to the laminated shell, and the inner and outer support structures are equivalent to the ring ribs and ribs. Thus, the energy equation of the reinforced ribbed conical shell model includes two parts: a laminated shell and a rib strip.

As shown in Fig. 9, since the outer-layer ring is directly bound to the radial bracket and the rest of the external support structure is bound to the radial bracket, the radial bracket can be equivalent to a rib outside the laminated shell. In addition, the support ring is bound to the upper winding structure, so it can be equivalent to the ring rib inside the laminated shell.

[See PDF for image]

Fig. 9

Three-dimensional model diagram of winding and support structure

The equivalence of stiffeners and ribs are based on the discrete element method [37, 38], [40]. According to the geometric model of stiffened and ribbed conical shell as shown in Fig. 3, the sectional view of stiffeners and ribs are shown in Fig. 10. The number of stiffeners distributed on the whole shell is N_s, in which d_si, b_si and θ_i are the thickness, the width and the coordinate position of the θ direction of the i-th stiffener, respectively. Moreover, set the physical parameters of the i bar as follows: ρ_si, E_si, μ_si and G_si are density, elastic modulus, Poisson ratio and shear modulus of elasticity (G_si = E_si/(2(1 + μ_si))). Similarly, the number of ribs is N_r, and d_ri, b_ri and x_i are the thickness, the width and the coordinate position of the x direction of i-th rib, respectively. And the other physical parameters are ρ_ri, E_ri, μ_ri and G_ri. Among them, the physical parameters are taken from Table 1, and the geometrical parameters are taken from the finite element model in subsection 2.1.

[See PDF for image]

Fig. 10

Sectional view of stiffeners and ribs

The cross-sectional area of the i-th stiffener(rib) is A_s(r)i and A_s(r)i = d_s(r)ib_s(r)i. The total kinetic energy T_s of N_s stiffeners can be written as

20

$$T_s=\frac{1}{2}\sum _{i=1}^{N_s}{\rho }_{\mathit{si}}{A}_{\mathit{si}}{\int }_0^L\left[{\left(\frac{\partial u_{\mathit{si}}}{\partial t}\right)}^2+{\left(\frac{\partial v_{\mathit{si}}}{\partial t}\right)}^2+{\left(\frac{\partial w_{\mathit{si}}}{\partial t}\right)}^2\right]dx$$

The total strain energy U_s can be written as

21

$$U_s=\frac{1}{2}\sum _{i=1}^{N_s}{\int }_0^L\left[D_{\mathit{si}}{\left({\epsilon }_x^0\right)}^2-2S_{\mathit{si}}{\epsilon }_x^0{k}_x^0+K_{\mathit{si}}{\left(k_x^0\right)}^2+P_{\mathit{si}}{\tau }_s^2\right]dx$$

The total kinetic energy T_r of N_r ribs can be written as

22

$$T_r=\frac{1}{2}\sum _{i=1}^{N_r}{\rho }_{\mathit{ri}}{A}_{\mathit{ri}}{\int }_0^{2\pi }\left[{\left(\frac{\partial u_{\mathit{ri}}}{\partial t}\right)}^2+{\left(\frac{\partial v_{\mathit{ri}}}{\partial t}\right)}^2+{\left(\frac{\partial w_{\mathit{ri}}}{\partial t}\right)}^2\right]R_{\mathit{ri}}d\theta$$

The total strain energy U_r can be written as

23

$$U_r=\frac{1}{2}\sum _{i=1}^{N_r}{\int }_0^{2\pi }\left[D_{\mathit{ri}}{\left({\epsilon }_{\theta }^0\right)}^2-2S_{\mathit{ri}}{\epsilon }_{\theta }^0{k}_{\theta }^0+K_{\mathit{ri}}{\left(k_{\theta }^0\right)}^2+P_{\mathit{ri}}{\tau }_r^2\right]R_{\mathit{ri}}d\theta$$where symbol definition of energy function of stiffened and ribbed can be found from Appendix D.

The total kinetic energy of the stiffened and ribbed conical shell model can be written as

24

$$T=T_{\epsilon }+T_r+T_s$$

The total strain energy of the stiffened and ribbed conical shell model can be written as

25

$$U=U_{\epsilon }+U_b+U_r+U_s$$

The establishment of equation of forced vibration and frequency response function

The external force f(x,θ) is applied at any coordinate (x,θ) on the surface of the stiffened and ribbed conical shell model. This force can be expressed as the component f_u, f_v and f_w, respectively, along the generatrix direction of the conical shell surface, the tangential direction of radius of the conical shell and the vertical direction of the conical shell surface. The potential energy W caused by force can be written as

26

$$W={\int }_0^L{\int }_0^{2\pi }[f_u(x,\theta )u+f_v(x,\theta )v+f_w(x,\theta )w]R(x)d\theta dx$$

Point force is the typical loading case, which is frequently encountered in practice. When the external force is a point load and the form is a harmonic excitation with frequency ω, change Eq. (26) to the point force form as follows:

27

$$\begin{array}{c}f_{_{u(v,w)}^{}}^{}=F_{u(v,w)}\delta (x-x_f,\theta -{\theta }_f)e^{j\omega t}\\ W_{}^{}={\int }_0^L{\int }_0^{2\pi }(F_{u}u+F_{v}v+F_{w}w)R(x)\delta (x-x_f,\theta -{\theta }_f)e^{j\omega t}d\theta dx\\ \end{array}$$where x_f and θ_f are the positions of action of point force. δ is the Dirac function. F_u(v,w) is the magnitude of the amplitude of the component force in each direction. ω is the frequency of the harmonic excitation.

By combining all the energy variations described above, the entire Lagrangian energy function can be written as follows:

28

$$L=(U_{\epsilon }+U_r+U_s)+U_b-(T_{\epsilon }+T_r+T_s)+W$$

By Substituting Eqs. (6)–(7), (9), (20–23) and (26) together with the admissible functions defined in Eqs. (4)–(5) into Eq. (28), according to the solution procedure of the Ritz method and minimizing the entire Lagrangian energy with respect to the unknown coefficients,

29

$$\frac{\partial \mathit{L}}{\partial \mathbf{\Phi }}=\mathbf{0},\mathbf{\Phi }={[A_0,A_1,\dots ,A_m,a_1,a_2,B_0,B_1,\dots ,B_m,b_1,b_2,C_0,C_1,\dots ,C_m,c_1,c_2,c_3,c_4]}^T.$$

Then, in matrix format, they can be summarized as follows:

30

$$\left(\mathbf{K}-{\omega }^2\mathbf{M}\right)\mathbf{X}=\mathbf{F}$$where K and M are stiffness and mass matrices, respectively, and F is the external force matrix. X denotes the unknown coefficient vector, and it is written as follows:

31

$$\mathbf{X}={[{\mathbf{u}}^T{,\mathbf{v}}^T{,\mathbf{w}}^T]}^T,\left\{\begin{array}{c}\mathbf{u}={\left[A_0,A_1,\dots ,A_m,a_1,a_2\right]}^T\\ \mathbf{v}={\left[B_0,B_1,\dots ,B_m,b_1,b_2\right]}^T\\ \mathbf{w}={\left[C_0,C_1,\dots ,C_m,c_1,c_2,c_3,c_4\right]}^T\\ \end{array}\right)$$

Then, the Eq. (30) is partially expanded and derived. First, kinetic energy and potential energy can be expressed in terms of unknown coefficients as follows

32

$$\begin{array}{c}U=\frac{1}{2}{\mathbf{X}}^T\mathbf{KX},\mathbf{K}={\mathbf{K}}_{\epsilon }+{\mathbf{K}}_r+{\mathbf{K}}_s+{\mathbf{K}}_b,\mathbf{K}=\frac{{\partial }^{2}U}{\partial {\mathbf{X}}^2},{\mathbf{K}}_{\epsilon (r,s,b)}=\frac{{\partial }^2{U}_{\epsilon (r,s,b)}}{\partial {\mathbf{X}}^2}\\ T=\frac{{\omega }^2}{2}{\mathbf{X}}^T\mathbf{MX},\mathbf{M}={\mathbf{M}}_{\epsilon }+{\mathbf{M}}_r+{\mathbf{M}}_s,\mathbf{M}=\frac{{\partial }^{2}T}{\partial {\mathbf{X}}^2},{\mathbf{M}}_{\epsilon (r,s)}=\frac{{\partial }^2{T}_{\epsilon (r,s)}}{\partial {\mathbf{X}}^2}\\ W={\mathbf{F}}^T\mathbf{X},\mathbf{F}=\frac{\partial W}{\partial \mathbf{X}}\\ \end{array}$$where, K_ε, K_r, K_s, and K_b are stiffness matrix of strain energy of the shell, stiffness matrix of strain energy of the rib, stiffness matrix of strain energy of the stiffener and stiffness matrix of elastic potential energy of the spring, respectively. M_ε, M_r and M_s are mass matrix of kinetic energy of the shell, mass matrix of kinetic energy of the rib and mass matrix of kinetic energy of the stiffener, respectively. The matrix is further expanded as follows:

33

$$\begin{array}{c}{\mathbf{M}}_{\epsilon (r,s)}={\left[\begin{array}{ccc}{\mathbf{M}}_{\epsilon (r,s)}^{\mathit{uu}}& 0& 0\\ 0& {\mathbf{M}}_{\epsilon (r,s)}^{\mathit{vv}}& 0\\ 0& 0& {\mathbf{M}}_{\epsilon (r,s)}^{\mathit{ww}}\end{array}\right]}_{3m+11}\\ {\mathbf{M}}_{\epsilon (r,s)}^{\mathit{uu}}=\frac{{\partial }^2{T}_{\epsilon (r,s)}}{\partial {\mathbf{u}}^2},{\mathbf{M}}_{\epsilon (r,s)}^{\mathit{vv}}=\frac{{\partial }^2{T}_{\epsilon (r,s)}}{\partial {\mathbf{v}}^2},{\mathbf{M}}_{\epsilon (r,s)}^{\mathit{ww}}=\frac{{\partial }^2{T}_{\epsilon (r,s)}}{\partial {\mathbf{w}}^2}\\ \end{array}$$

34

$$\begin{array}{c}{K}_{\epsilon (r,s,b)}={\left[\begin{array}{ccc}{K}_{\epsilon (r,s,b)}^{\mathit{uu}}& K_{\epsilon (r,s,b)}^{\mathit{uv}}& K_{\epsilon (r,s,b)}^{\mathit{uw}}\\ K_{\epsilon (r,s,b)}^{\mathit{uvT}}& K_{\epsilon (r,s,b)}^{\mathit{vv}}& K_{\epsilon (r,s,b)}^{\mathit{vw}}\\ K_{\epsilon (r,s,b)}^{\mathit{uwT}}& K_{\epsilon (r,s,b)}^{\mathit{vwT}}& K_{\epsilon (r,s,b)}^{\mathit{ww}}\end{array}\right]}_{3m+11}\\ K_{\epsilon (r,s,b)}^{\mathit{uu}}=\frac{{\partial }^2{U}_{\epsilon (r,s,b)}}{\partial u^2},K_{\epsilon (r,s,b)}^{\mathit{vv}}=\frac{{\partial }^2{U}_{\epsilon (r,s,b)}}{\partial v^2},K_{\epsilon (r,s,b)}^{\mathit{ww}}=\frac{{\partial }^2{U}_{\epsilon (r,s,b)}}{\partial w^2}\\ K_{\epsilon (r,s,b)}^{\mathit{uv}}=\frac{{\partial }^2{U}_{\epsilon (r,s,b)}}{\partial u\partial v},K_{\epsilon (r,s,b)}^{\mathit{uw}}=\frac{{\partial }^2{U}_{\epsilon (r,s,b)}}{\partial u\partial w},K_{\epsilon (r,s,b)}^{\mathit{vw}}=\frac{{\partial }^2{U}_{\epsilon (r,s,b)}}{\partial v\partial w}\\ \end{array}$$

By substituting Eqs. (31-34) into Eq. (30), the detailed characteristic equation can be obtained as shown below:

35

$$\left\{\left[\begin{array}{ccc}\sum _{i=\epsilon ,r,s,b}{\mathbf{K}}_i^{\mathit{uu}}& \sum _{i=\epsilon ,r,s,b}{\mathbf{K}}_i^{\mathit{uv}}& \sum _{i=\epsilon ,r,s,b}{\mathbf{K}}_i^{\mathit{uw}}\\ \sum _{i=\epsilon ,r,s,b}{\mathbf{K}}_i^{\mathit{uvT}}& \sum _{i=\epsilon ,r,s,b}{\mathbf{K}}_i^{\mathit{vv}}& \sum _{i=\epsilon ,r,s,b}{\mathbf{K}}_i^{\mathit{vw}}\\ \sum _{i=\epsilon ,r,s,b}{\mathbf{K}}_i^{\mathit{uwT}}& \sum _{i=\epsilon ,r,s,b}{\mathbf{K}}_i^{\mathit{vwT}}& \sum _{i=\epsilon ,r,s,b}{\mathbf{K}}_i^{\mathit{ww}}\end{array}\right]-{\omega }^2\left[\begin{array}{ccc}\sum _{i=\epsilon ,r,s}{\mathbf{M}}_i^{\mathit{uu}}& \mathbf{0}& \mathbf{0}\\ \mathbf{0}& \sum _{i=\epsilon ,r,s}{\mathbf{M}}_i^{\mathit{vv}}& \mathbf{0}\\ \mathbf{0}& \mathbf{0}& \sum _{i=\epsilon ,r,s}{\mathbf{M}}_i^{\mathit{ww}}\end{array}\right]\right\}\left[\begin{array}{c}\mathbf{u}\\ \mathbf{v}\\ \mathbf{w}\end{array}\right]=\mathbf{F}$$

When the external force F is 0, Eq. (30) is

36

$$\left(\mathbf{K}-{\omega }^2\mathbf{M}\right)\mathbf{X}=0$$

By solving Eq. (36), the natural frequency ω and the coefficient X (A_m, a₁, a₂, B_m, b₁, b₂, C_m, c₁, c₂, c₃ and c₄) can be obtained with the corresponding wave number n. In addition, substituting the obtained coefficient into Eq. (5), the corresponding mode shape can be determined.

When the excitation frequency is ω_F and the forces in each direction are f_u, f_v and f_w. Meanwhile, Rayleigh damping is introduced to establish the damping matrix C, Eq. (30) is

37

$$\left(\mathbf{K}-{\omega }_F^2\mathbf{M}+j{\omega }_F\mathbf{C}\right)\mathbf{X}=\mathbf{F},\mathbf{C}=\alpha \mathbf{M}+\beta \mathbf{K}$$where the damping coefficients α and β are calculated by the following formula [39]:

38

$$\left[\begin{array}{c}\alpha \hfill \\ \beta \hfill \end{array}\right]=\frac{2{\omega }_m{\omega }_n}{{\omega }_n^2-{\omega }_m^2}\left[\begin{array}{cc}{\omega }_n& -{\omega }_m\\ -1/{\omega }_n& 1/{\omega }_m\end{array}\right]\left[\begin{array}{c}{\xi }_m\hfill \\ {\xi }_n\hfill \end{array}\right]$$where ω_m and ω_n are the natural frequencies corresponding to the m and n modes respectively. ζ_m and ζ_n are the damping ratios of the m and n modes respectively. The damping ratios of each order of winding structure in practical engineering range from 0.01 to 0.03 [10], so the damping ratio is set as 0.02 in this paper.

Further, Eq. (37) can be expressed as:

39

$$\mathbf{X}=\mathbf{H}\left({\omega }_F\right)\mathbf{F},\mathbf{H}\left({\omega }_F\right)={\left(\mathbf{K}-{\omega }_F^2\mathbf{M}+j{\omega }_F\mathbf{C}\right)}^{-1}$$

When the external force is a harmonic excitation of a point load, the external force matrix F can be described as follows:

40

$$\begin{array}{c}\mathbf{F}={\left[{\mathbf{F}}_{\mathbf{u}},{\mathbf{F}}_{\mathbf{v}},{\mathbf{F}}_{\mathbf{w}}\right]}^T,{\mathbf{F}}_{\mathbf{u}(\mathbf{v},\mathbf{w})}=\frac{\partial W}{\partial {\mathrm{\Theta }}_{\mathbf{u}(\mathbf{v},\mathbf{w})}}=F_{u(v,w)}{\psi }_{u(v,w)}{e}^{j\omega t}\\ {\psi }_{u(v,w)}={\int }_0^L{\int }_0^{2\pi }\frac{\partial {\vartheta }_{u(v,w)}}{\partial {\mathrm{\Theta }}_{\mathbf{u}(\mathbf{v},\mathbf{w})}}R(x)\delta (x-x_f,\theta -{\theta }_f)d\theta dx\\ \end{array}$$where, define abbreviations: X_u(v,w) = X_u, X_v or X_w; X_u(v,w) = X_u, X_v or X_w; let ${\mathrm{\Theta }}_{\mathbf{u}}=\mathbf{u},{\mathrm{\Theta }}_{\mathbf{v}}=\mathbf{v},{\mathrm{\Theta }}_{\mathbf{w}}=\mathbf{w};$${\vartheta }_u=u,{\vartheta }_v=v,{\vartheta }_w=w$, and the definition of $\mathbf{u},\mathbf{v},\mathbf{w}$ refers to Eq. (31), whereas the definition of $u,v,w$ refers to Eq. (4).

By solving Eq. (37), the coefficient X under external excitation can be obtained. In addition, substituting the obtained coefficient into Eq. (5), the steady state responses of each point in each direction under external force excitation can be calculated, as follows:

41

$$\text{by}_{_{}^{}}^{}\text{Eq}.(20),\left\{\begin{array}{c}u=U_n(x)cos(n\theta )e^{j{\omega }_{F}t}\\ v=V_n(x)cos(n\theta )e^{j{\omega }_{F}t}\\ w=W_n(x)cos(n\theta )e^{j{\omega }_{F}t}\\ \end{array}\right)$$$$\begin{array}{c}and,\left\{\begin{array}{c}{U}_n=\left[\begin{array}{cccc}{A}_{0}cos\frac{0\pi }{L}x& ...& A_{m}cos\frac{m\pi }{L}x& \begin{array}{cc}{f}_1(x)& f_2(x)\end{array}\end{array}\right]u\\ V_n=\left[\begin{array}{cccc}{B}_{0}cos\frac{0\pi }{L}x& ...& B_{m}cos\frac{m\pi }{L}x& \begin{array}{cc}{f}_1(x)& f_2(x)\end{array}\end{array}\right]v\\ W_n=\left[\begin{array}{cccc}{C}_{0}cos\frac{0\pi }{L}x& ...& C_{m}cos\frac{m\pi }{L}x& \begin{array}{cccc}{f}_1(x)& f_2(x)& f_3(x)& f_4(x)\end{array}\end{array}\right]w\\ \end{array}\right)\\ substituting\text{int}o,_{}^{}{[u^T,v^T,w^T]}^T=X=H\left({\omega }_F\right)F\\ \end{array}$$

The simulation results and comparative analysis

Natural frequency

According to the actual working state of the stator end winding, Clamped and Free (C-F) edges are applied at the iron core end and the lead ring end of the conical shell, as well as the corresponding physical and geometrical parameters are set. Based on these, the natural vibration equations of the double-layered winding structure (shell) and the whole structure are solved respectively. Specifically, the geometrical parameters are shown as follows: core end radius of cone shell model R₁ = 0.7352 m, half-cone Angle α = π/16, generatrix length L = 1.3445 m, thickness H = 0.1095 m, θ₁ = 28.748°, thickness of wire rod and interlayer ring H_xb = 0.05892 m, H_cjh = 0.0426 m, and width ratio L_m0 = 0.1258. The height and width of stiffener and rib are: d_si = 0.42066 m, b_si = 0.05 m, d_ri = 0.0409 m, and b_ri = 0.05967 m. The calculation results are compared with the solution results of the finite element model, as shown in Table 4.

Table 4. Comparison for different structures of calculation results of natural frequency (Hz)

Mode shape	Shell (laminated)			whole structure
	Analytical solution by empirical formula	Analytical solution by IA	Finite element solution	Analytical solution by IA	Finite element solution
Elliptic mode	112.93	118.96	118.86	86.17	84.47
Triangular mode	158.59	162.37	161.34	109.52	114.64

It can be seen that the analytical solution of each natural frequency of the digital mechanism model established in this paper is in good agreement with that of the finite element model. Therefore, the digital mechanism model of stator end winding established can reflect the actual vibration characteristics of stator end winding. The empirical formula and heuristic algorithm are used to comprehensively calculate the equivalent boundary, and the calculation results of corresponding natural frequency are shown in Table 4 (in Table 7, parameters with group number 1). Some optimization details and key points of the optimization model of complex boundary in Sect. 4.3 are shown below.

The optimization model in Eq. (19) is evaluated using the three aforementioned algorithms under identical parameter configurations. The tolerance threshold $\epsilon$ is fixed at 0.05, with constraint limits $\mathrm{\Delta }{S}_u,$$\mathrm{\Delta }{S}_v,$$\mathrm{\Delta }{S}_w,$$\mathrm{\Delta }{S}_L$ assigned values of 0.18, 0.01, 0.055, and 0.02, respectively. These limits require model-specific calibration. Initially, a uniform adjustment is implemented to ensure stable optimization results (with natural frequency deviations within tolerance). This is then followed by mode-shape comparisons, including u-direction constraint relaxation, for precise tuning. As evidenced by Fig. 8, the u-direction spring exerts the strongest influence on natural frequency, justifying its higher constraint limit. Setting this limit too low may force the small-diameter end's vibration mode to resemble a fixed support, artificially inflating the natural frequency. Thus, optimal limits must balance frequency deviation and modal-shape fidelity. The model undergoes iterative optimization via all three algorithms, with results depicted in Fig. 11(c) and algorithm parameters detailed in Table 5.

[See PDF for image]

Fig. 11

Comparison of free vibration frequencies in double-layer winding

Table 5. Parameter configuration

LDIW-PSO		RIW-PSO		IA
Learning factor c1	1.5	Learning factor c1	1.5	Variation probability pm	0.7
Learning factor c2	1.5	Learning factor c2	1.5	Excitation constant α	1
Inertia weight wmax	0.9	Mean weight μmax	0.8	Excitation constant β	1
Inertia weight wmin	0.2	Mean weight μmin	0.4	Similarity threshold δs	0.2
Number of iterations	200	Number of iterations	200	Number of iterations	200
Population size	20	Population size	20	Population size	20
Velocity vmax	2.6	Velocity vmax	2.6	Number of clones	10
Velocity vmin	− 2.6	Velocity vmin	−2.6
		variance σ	0.2

Analysis of Fig. 11(c) reveals that all three algorithms initially converge to high-quality solutions. However, RIW-PSO exhibits a tendency to become trapped in local optima with limited escape capability. While LDIW-PSO demonstrates gradual improvement in solution quality, it requires significantly more iterations than the immune algorithm (IA) to achieve comparable results. The IA demonstrates superior robustness and global search performance, as evidenced by its early convergence to optimal solutions in Fig. 11(c), with local convergence occurring only after 100 generations. Corresponding optimization results are presented in Table 6.

Table 6. Optimization result of equivalent stiffness of winding

	ku0	kv0	kw0	kθ0
IA spring stiffness unit: N/m	9.9183e8	1.8539e10	1.1455e9	2.7941e6
LDIW-PSO spring stiffness unit: N/m	9.4629e8	0.8802e10	1.2279e9	4.1779e6
RIW-PSO spring stiffness unit: N/m	7.0151e8	1.3698e9	2.1346e9	1.3553e7

Table 6 reveals that the solutions obtained by IA and LDIW-PSO are largely comparable, but IA achieves a higher fitness value. Figure 11(c) reveals RIW-PSO's tendency to converge to suboptimal solutions. Notably, optimization processes are prone to becoming trapped in local minima, making it challenging to escape. In this case, the stiffness parameters k_u0 and k_w0 found by RIW-PSO are similar to those obtained by IA and LDIW-PSO. However, significant discrepancies are observed in k_u0 and k_w0, indicating that RIW-PSO succumbed to local convergence in the later optimization stages. In conclusion, the stiffness values identified by IA are adopted as the optimal results. The spring stiffness can be calculated using the empirical formula Eq. (8), and the displacement in all directions can be obtained by substituting the stiffness into Eq. (13), and then the displacement field Fig. 13(a)-(b)) and modal parameters (Fig. 14(c)) can be calculated. Under the same steps, the displacement field (Fig. 13(c)–(d)) and modal parameters (Fig. 14(d)) corresponding to the optimized model can be calculated.

Furthermore, the free vibration frequencies of the fiber and matrix under different material parameters are calculated. Parameter settings are shown in Table 7. The accuracy of parametric modeling of the proposed mathematical model was verified by changing the parameters’ settings, and rationality of the equivalent composite structure adopted was demonstrated. The nonlinear relationship between natural frequency and material parameters can be characterized by this equivalent method. The data comparison is shown in Fig. 11.

Table 7. Parameter setting

Group number	Fiber (stator bar) parameters		Matrix (ring) parameters
	Exb (MPa)	ρxb (kg/m3)	Ecjh (MPa)	ρcjh (kg/m3)
1 (original)	40,000	4033	39,045	1893
2	40,000	4033	30,000	1893
3	40,000	4033	40,000	1893
4	40,000	4033	50,000	1893
5	40,000	4033	39,045	893
6	40,000	4033	39,045	2893
7	40,000	4033	39,045	3893
8	30,000	4033	39,045	1893
9	50,000	4033	39,045	1893
10	60,000	4033	39,045	1893
11	40,000	3033	39,045	1893
12	40,000	5033	39,045	1893
13	40,000	6033	39,045	1893
14	40,000	7033	39,045	1893
15	30,000	3033	30,000	893
16	30,000	7033	30,000	3893
17	60,000	7033	50,000	3893
18	60,000	3033	50,000	893
19	50,000	6033	35,000	2893
20	35,000	4500	50,000	2500

The data of three boundary configurations are compared. Figure 11(a) presents free constraint (F-F), that is, the finite element model is unconstrained, and the spring stiffness of the mathematical model is 0. Figure 11(b) and (d) show the actual boundary situation of the winding (C-F), that is, the small-diameter end of the finite element model has fixed constraints, while the large-diameter end has no constraints. The small-diameter end of the mathematical model is the elastic constraint, and the spring stiffness setting is calculated according to the empirical formula in Fig. 11(b), and the optimized model in Fig. 11(d), respectively, and the spring stiffness of the large-diameter end is 0.

In particular, the first boundary condition is unconstrained at both ends (F-F), so the boundary configuration does not require a mathematical model to establish an equivalent boundary. The modeling of equivalent boundary and equivalent mechanism will introduce some errors. Therefore, when the boundary conditions are set to free (the first boundary condition), the error that may be introduced by the equivalent boundary will be excluded, and the accuracy of the equivalent mechanism modeling can be judged. Therefore, this method is necessary and meaningful to verify the mechanism modeling, and can lay a solid foundation for verifying the equivalent boundary.

It can be seen from Fig. 11(a) that the method can accurately calculate the natural frequency of the structure without constraints at both ends when the settings of material parameter are changed, which proves the rationality of the equivalent method and the accuracy of the model. Therefore, the method can reflect the influence of material parameter characteristics on the natural frequency of the structure. Further, use of the empirical formula Eq. (8) can reflect the complex boundary constraints to a certain extent, which proves the rationality of the formula in Fig. 11(b); then, establishes an optimization model of spring in complex boundary and uses the heuristic algorithm in Fig. 11(d) to optimize the stiffness.

By comparing Fig. 11(b) and (d), it can be seen that the proposed optimization of complex boundary can obtain better solutions. Based on the data in Fig. 11(b) and (d), the relative error of natural frequency after optimization of stiffness of complex boundary is calculated, as shown in Fig. 12. The red dotted line is the accuracy requirement ${\delta }^E$ of engineering, and ${\delta }^E=0.05$ is set in this paper. Define the relative error calculated as ${\delta }_i=\left|f_i-f_i^{\ast }\right|\left./\phantom{\left|f_i-f_i^{\ast }\right|f_i^{\ast }}\right)f_i^{\ast }$ using material parameters of Group i, where $f_i$ is the natural frequency calculated by substituting parameters of Group i into the mathematical model, and $f_i^{\ast }$ is the natural frequency of modal analysis of finite element under parameters of group i.

[See PDF for image]

Fig. 12

Relative error of natural frequency before and after optimization

Further, the average value of the relative error under each wave-number before and after optimization is calculated: ${\delta }_n^{\mathit{average}}=1\left./\phantom{1N}\right)N{\sum }_{i=1}^N{\delta }_{\mathit{ni}}$, where ${\delta }_{\mathit{ni}}$(${\delta }_{\mathit{ni}}=\left|f_{\mathit{ni}}-f_{\mathit{ni}}^{\ast }\right|\left./\phantom{\left|f_{\mathit{ni}}-f_{\mathit{ni}}^{\ast }\right|f_{\mathit{ni}}^{\ast }}\right)f_{\mathit{ni}}^{\ast }$) is the relative error of group i with wave number n, and N is the number of groups, then ${\delta }_n^{\mathit{average}}$ is the average relative error at wave number n. Before optimization, ${\delta }_2^{\mathit{average}}$ = 0.0691, ${\delta }_3^{\mathit{average}}$ = 0.0293; after optimization, ${\delta }_2^{\mathit{average}}$ = 0.0350, ${\delta }_3^{\mathit{average}}$ = 0.0141. It can be seen that the optimization method can effectively reduce the average relative error and improve the accuracy of the model. When $\delta$ > ${\delta }^E$, the precision is considered to be low, and the number of groups with low precision in the data set is recorded as N^LE. Before optimization, N^LE = 13 for n = 2 and N^LE = 2 for n = 3; after optimization, N^LE = 5 for n = 2 and N^LE = 0 for n = 3. It can be seen that after optimization, the relative error of most groups is guaranteed to meet ${\delta }^E$. Therefore, it proves the effectiveness of the proposed method that the relative error can be effectively reduced by optimization. The boundary optimization model with heuristic algorithm can realize the equivalent complex boundary under the complex structure, and can obtain the solution set that meets the engineering requirements considering natural frequency and mode (refer to Sect. 6.2).

Mode shapes

By substituting the coefficient vector obtained by solving the characteristic equation into the vibration mode function Eq. (4), the mode shapes of each order of the double-layer winding structure and whole structure under the digital mechanism model can be obtained. The spring stiffness calculated by the empirical formula and the spring stiffness optimized by the optimization model are respectively brought into the characteristic equation to calculate the displacement field, as shown in Fig. 13.

[See PDF for image]

Fig. 13

Displacement field of winding structure

From the displacement field presented, it can be seen that the boundary calculated by IA algorithm is closer to the actual situation of the winding, that is, the displacement of the small end (x = 0) approaches 0. Therefore, the spring stiffness optimized by IA algorithm is taken as the equivalent complex boundary of the model. Take the elliptic and triangular modes and the mode shapes of double-layer winding structure are shown in Fig. 14, and the mode shape of the whole structure is shown in Fig. 15. Figure 14 presents finite element analysis and the analytical solution under Clamped (Fig. 14(a)–(b)), the analytical solution of elastic boundary obtained under empirical formula (Fig. 14(c)) and the analytical solution of elastic boundary obtained under optimization (Fig. 14(d)). The selection of the overall winding node of analysis is referred to Fig. 15, and the vibration mode vector is extracted from the finite element analysis and compared with the numerical calculation results of the mathematical model, as shown in Fig. 15.

[See PDF for image]

Fig. 14

Comparison of the mode shapes of double-layer winding structure

[See PDF for image]

Fig. 15

Comparison of the mode shapes of whole structure

As can be seen from Fig. 14, the mode shapes simulated by using the ANSYS is consistent with the mode shapes obtained by using the MATLAB at the same order of mode. It can be shown that for the double-layer winding structure, the continuous truncated conical shell model of the laminated composite material established in this paper can reflect its vibration characteristics, which proves the correctness and effectiveness of the equivalent model.

Figure 15 shows the comparison between the modal shapes of the overall winding structure established in this paper and the finite element method. It can be seen that the modal matching degree of the two is high, which further indicates that the proposed method can be applied to complex winding structures and calculate its modal parameters. The formula for calculating mode shapes in Fig. 15(f) and (g) is as follows:

42

$${\varphi }_{ix(y,z)}={\varphi }_{ix(y,z)}\left./\phantom{{\varphi }_{ix(y,z)}\sqrt{{\left(p{\varphi }_{12x}\right)}^2}}\right)\sqrt{{\left(p{\varphi }_{12x}\right)}^2}$$where, p is the scaling factor, 30 is set in this paper, and ${\varphi }_{ix(y,z)}$ is the mode shape vector of the i node in the x(y,z) direction respectively. Based on the displacement of the 12th node in the x direction, the vibration mode of the 36 nodes in three directions is calculated. See Tables in Appendix E for specific calculation results.

Therefore, the method of calculating the free vibration of the whole winding proposed can help researchers to study the winding structure in the early stage, and to explore the vibration characteristics of the whole winding structure, and greatly save the development time and improve the development efficiency under the early qualitative analysis.

Frequency response function

The stiffness matrix and mass matrix in Eq. (39) can be uniquely obtained through structural parameters determined. When the position (x,θ) of action of the fixed-point load is given, x_f and θ_f in Eq. (27) are uniquely determined. When the magnitude and direction of the fixed-point load f(x,θ) are given, F_u(v,w) in Eq. (27) is uniquely determined (the force is expressed in component form). Therefore, the external force matrix F can be calculated by Eq. (40).

By substituting the mass matrix M, stiffness matrix K, damping matrix C, excitation force matrix F and excitation frequency ω_F into Eq. (39), the coefficient matrix X of the structure under this excitation frequency can be calculated. By substituting X into Eq. (4), the displacement response of each node can be obtained. To analyze the response in the excitation frequency range ω_L Hz—ω_H Hz, the analysis step size ω_step Hz can be set. When the response point coordinates are determined, the displacement response of the response point in the excitation frequency range can be obtained by calculating the frequency calculated by step length in the frequency range.

In order to analyze the forced vibration of the structure, it is necessary to solve the frequency response function of the structure under excitation. Because the forced vibration characteristics of the stator end winding are mainly reflected in the shell, the excitation point Node 1 coordinates (x_f,θ_f) = (x_Node1,θ_Node1) are taken at the lead-ring end of the shell, and the excitation force f_x(y,z) (coordinate-system refers to Fig. 15(e), and conversion of excitation-force refers to Eq. (43)) is a steady state single-point excitation (frequency is ω_F and varies with the response frequency band), which is $f_z=F_z\delta (x-x_f,\theta -{\theta }_f)e^{j\omega t}$ (F_z = 100 N). Take the response points Node 4 and Node 20 as (x,θ) = (x_Node4,20,θ_Node4,20) an example, the displacement response in the x-direction(transformation of response coordinate refers to Eq. (44)) is shown in the Fig. 16. Taking the response frequency band (ω_L = 75, ω_H = 125 and ω_step = 0.5).

[See PDF for image]

Fig. 16

Comparison of displacement response, taking Node 4 and Node 20 as an example

The coordinate system transformation of excitation and response is as follows:

43

$$\left\{\begin{array}{c}{f}_u=\left(f_{x}cos{\theta }_f+f_{y}sin{\theta }_f\right)sin\alpha +f_{z}cos\alpha \\ f_v=f_{y}cos{\theta }_f-f_{x}sin{\theta }_f\\ f_w=\left(f_{x}cos{\theta }_f+f_{y}sin{\theta }_f\right)cos\alpha -f_{z}sin\alpha \\ \end{array}\right)$$

44

$$\left\{\begin{array}{c}x=\left(usin\alpha +wcos\alpha \right)cos\theta -vsin\theta \\ y=\left(usin\alpha +wcos\alpha \right)sin\theta +vcos\theta \\ z=ucos\alpha -wsin\alpha \\ \end{array}\right)$$

It can be seen from Fig. 16 that the frequency response function calculated by the digital mechanism model and obtained by finite element simulation will reach the local peak when the frequency of excitation force approaches the natural frequency. Further, the responses of all nodes were calculated and two-dimensional surface plots were drawn to express the variation of displacement response in frequency domain and spatial domain, as shown in Fig. 17

[See PDF for image]

Fig. 17

Displacement response surface of all Nodes (excitation point is Node 1, excitation force F_x = 100 N)

Figure 17 shows the displacement responses of 36 Nodes in each direction when the exciting force of Node 1 is F_x = 100 N. It includes the top view, which can observe the change rule of displacement response in frequency domain and spatial domain (node information). It is worth noting that the law of calculation is consistent with the law of finite element analysis, which proves that the complex law of displacement response of complex winding structures on multiple scales can be captured by the method. However, it is extremely difficult to establish an accurate mathematical model to calculate the response to ensure that the response accuracy meets the requirements. One of the reasons is that the winding structure is complex, including the coupling of various structures. The second reason is that the mathematical model contains equivalent errors, calculation errors and neglects the influence of some structures. Of course, the method presented is helpful in the qualitative analysis of the vibration of the stator end winding of the large turbogenerator. The vibration law can be analyzed by the method presented, and a more accurate response model can be established based on the further improvement of the method.

Conclusion

According to the modeling approach of composite material theory, the stator bar can be equivalent to fiber, and the interlayer ring can be equivalent to matrix, which was considered to fill the gap between the stator bars. Based on the model with flake tandem, it can be equivalent to the composite conical shell model.

The dynamic characteristics of the equivalent conical-shell model were solved by using the Rayleigh Ritz method with improved Fourier series, which is suitable for different complex elastic boundary conditions. By studying the effect of elastic boundary, the sensitive interval was found. The optimized semi-analytical solution of the new optimization model of characteristic equation for stator-winding can find a better spring stiffness to simulate the actual complex boundary condition fully. Considering the influence of support structure fully, the support frame and support ring were regarded as stiffener and rib.

The stator end winding was equivalent to a double-layered stiffened and ribbed conical shell considering actual complex boundary condition, and the natural and forced vibration equations were established. The modal parameters calculated by the vibration model of double-layer winding established agreed well with the finite element analysis. Furthermore, the calculated modal parameters and the frequency response under the harmonic response excitation were fully derived.

Therefore, it is proved that the displacement response calculated by the mechanism model can reflect the complex evolution in the frequency domain and the space domain to a certain extent. And it showed that the digital mechanism model established can be applied to the study of the dynamic characteristics of the stator end winding.

The research work can further introduce the mechanism modeling of electromagnetic force of stator end winding, and obtain the analytical solution of the whole structure model under the electromagnetic force of current load. Therefore, combined with the digital mechanism modeling, the comprehensive cross-research will play an inestimable role in the dynamic monitoring driven by industrial big data.

Replication of results

Tables 1 and 2 clearly show the physical parameters and connection relationship of the finite element model structure of the 600 MW stator end winding studied in this paper. The finite element analysis software is ANSYS Workbench 19.2. For more details, please refer to the previous research [10, 11]. The mechanism model is derived in detail in this paper. MATLAB R2017b programming software is used for programming calculation. The structural parameters of isotropic conical shell refer to Sect. 4.3 and the structural parameters of mechanism model refer to Sect. 6.1. For equivalence verification, refer to Table 7, Fig. 11(a), Fig. 15, and Appendix A. Objective function of core optimization model reference Eq. (18), optimization algorithm reference LDIW-PSO [34], RIW-PSO [35] and IA [36], and specific parameter configuration refer to Table 5. Refer to Appendix B for the calculation of specific vibration modes in Fig. 15.

Funding

The work was supported by National Natural Science Foundation of China (Grant No.51807019), the Science and Technology Research Program of Chongqing Municipal Education Commission (KJZD-K201900604), the Open Projects of State Key Laboratory for Strength and Vibration of Mechanical Structures (Grant No. SV2020-KF-15), Chongqing Special Postdoctoral Science Foundation (Grant No.XmT2018030).

Data availability

The datasets generated during and analysed during the current study are available in the figshare repository, https://doi.org/https://doi.org/10.6084/m9.figshare.25802767.

Declarations

Conflict of interest

On behalf of all authors, the corresponding author states that there is no Conflict of interest.

Ethical approval

This study did not involve animal studies and human studies.

Adriano Todorovic Fabro

Publisher's Note

Springer Nature remains neutral with regard to jurisdictional claims in published maps and institutional affiliations.

References

1. Islam, MS; Islam, R; Sebastian, T. Noise and vibration characteristics of permanent-magnet synchronous motors using electromagnetic and structural analyses. IEEE T Ind Appl; 2014; 50, 5 pp. 3214-3222. [DOI: https://dx.doi.org/10.1109/TIA.2014.2305767]

2. Dong, Q; Liu, X; Qi, G; Sun, C; Wang, Y. Analysis and evaluation of electromagnetic vibration and noise in permanent magnet synchronous motor with rotor step skewing. Sci China Technol Sci; 2019; 62, pp. 839-848. [DOI: https://dx.doi.org/10.1007/s11431-018-9458-5]

3. He, YL; Xu, MX; Zhang, W; Wang, XL; Lu, P; Gerada, C; Gerada, D. Impact of stator interturn short circuit position on end winding vibration in synchronous generators. IEEE T Energy Conver; 2021; 36, 2 pp. 713-724. [DOI: https://dx.doi.org/10.1109/TEC.2020.3021901]

4. Gupta, N; Saini, M; Kumar, A. Operational availability analysis of generators in steam turbine power plants. SN Appl Sci; 2020; 2, 779. [DOI: https://dx.doi.org/10.1007/s42452-020-2520-y]

5. Letal, J; Satmoko, B; Manik, N; Stone, G. Stator end-winding vibration in two-pole machines: avoiding generator failure. IEEE Ind Appl Mag; 2020; 26, 6 pp. 29-39. [DOI: https://dx.doi.org/10.1109/MIAS.2020.2982725]

6. Zhou, Y; Zhang, QL; Duan, JG. Dynamical characteristics and influencing factors of stator end-windings of a turbine generator analyzed via heterogeneous element fusion modeling. Energy Rep; 2021; 7, pp. 658-672. [DOI: https://dx.doi.org/10.1016/j.egyr.2021.09.204]

7. GB/T 20140–2016, Measurement method and evaluation criteria of dynamic characteristic and vibration on stator end windings of cylindrical synchronous generators.

8. Liao W, Gao Y, Zhang H (2009) Analytical investigation of dynamic characteristic for stator end-winding of turbo-generator of 600MW. International conference on measuring technology and mechatronics automation, Zhangjiajie, pp 19–21. https://doi.org/10.1109/ICMTMA.2009.178

9. Sewak, R. Experimental modal analysis of stator overhangs of a large turbogenerator. Energy Power Eng; 2011; 3, 3 pp. 221-226. [DOI: https://dx.doi.org/10.4236/epe.2011.33028]

10. Zhao, Y; Yan, B; Zeng, C; Huang, S; Chen, C; Deng, J. Optimal scheme for structural design of large turbogenerator stator end winding. IEEE T Energy Conver; 2016; 31, 4 pp. 1423-1432. [DOI: https://dx.doi.org/10.1109/TEC.2016.2597151]

11. Zhao, Y; Yan, B; Chen, C; Deng, J; Zhou, Q. Parametric study on dynamic characteristics of turbogenerator stator end winding. IEEE T Energy Conver; 2014; 29, 1 pp. 129-137. [DOI: https://dx.doi.org/10.1109/TEC.2013.2294334]

12. Wang, YX; Wang, Y; Lin, L. Virtual prototype and modal analysis of stator system of large turbo-generator. Appl Mech Mater; 2012; 190, pp. 232-236. [DOI: https://dx.doi.org/10.4028/www.scientific.net/AMM.190-191.232]

13. Wang YX, Wang Y, Qiu H, Liu X (2009) Dynamic design and simulation analysis of 1000MW large turbo-generator. International conference on mechatronics and automation, Changchun, pp 1650–1655. https://doi.org/10.1109/ICMA.2009.5246294

14. Drubel, O; Kulig, S; Senske, K. End winding deformations in different turbo generators during 3-phase short circuit and full load operation. Electr Eng; 2000; 82, 3–4 pp. 145-152. [DOI: https://dx.doi.org/10.1007/s002020050005]

15. Iga, Y; Yamamoto, Y; Takahashi, K. Finite element modelling of turbine generator stator end windings for vibration analysis. Iet Electr Power Appl; 2016; 10, 2 pp. 75-81. [DOI: https://dx.doi.org/10.1049/iet-epa.2015.0142]

16. Liew, KM; Ng, TY; Zhao, X. Free vibration analysis of conical shells via the element-free kp-Ritz method. J Sound Vib; 2005; 281, 3–5 pp. 627-645. [DOI: https://dx.doi.org/10.1016/j.jsv.2004.01.005]

17. Li, FM; Kishimoto, K; Huang, WH. The calculations of natural frequencies and forced vibration responses of conical shell using the Rayleigh-Ritz method. Mech Res Commun; 2009; 36, 5 pp. 595-602. [DOI: https://dx.doi.org/10.1016/j.mechrescom.2009.02.003]

18. Jin, G; Ma, X; Shi, S; Ye, T; Liu, Z. A modified Fourier series solution for vibration analysis of truncated conical shells with general boundary conditions. Appl Acoust; 2014; 85, pp. 82-96. [DOI: https://dx.doi.org/10.1016/j.apacoust.2014.04.007]

19. Kouchakzadeh, MA; Shakouri, M. Free vibration analysis of joined cross-ply laminated conical shells. Int J Mech Sci; 2014; 78, pp. 118-125. [DOI: https://dx.doi.org/10.1016/j.ijmecsci.2013.11.008]

20. Ma, X; Jin, G; Xiong, Y; Liu, Z. Free and forced vibration analysis of coupled conical-cylindrical shells with arbitrary boundary conditions. Int J Mech Sci; 2014; 88, pp. 122-137. [DOI: https://dx.doi.org/10.1016/j.ijmecsci.2014.08.002]

21. Xing, ZZ; Wang, XH; Zhao, WL; Sun, LL; Niu, N. Calculation method for natural frequencies of stator of permanent magnet synchronous motors based on three-dimensional elastic theory. IEEE T Energy Conver; 2021; 36, 2 pp. 755-766. [DOI: https://dx.doi.org/10.1109/TEC.2020.3030042]

22. Hu, SL; Zuo, SG; Wu, H; Liu, MT. An analytical method for calculating the natural frequencies of a motor considering orthotropic material parameters. IEEE T Ind Electron; 2019; 66, 10 pp. 7520-7528. [DOI: https://dx.doi.org/10.1109/TIE.2018.2883194]

23. Zhao, GD; Sun, LL; Niu, N; Wang, XH; Xing, ZZ. Analysis of vibration characteristics of stators of electrical machines under actual boundary. Mech Syst Signal Pr; 2023; 185, [DOI: https://dx.doi.org/10.1016/j.ymssp.2022.109778] 109778.

24. Kumar PGS, Reddy KM, Kumar KK (2019) Reinforcement of generator stator end winding structure by modal analysis approach. International conference on condition assessment techniques in electrical systems, Chennai, pp 1–5. https://doi.org/10.1109/CATCON47128.2019.CN0071

25. Jones, RM. Mechanics of Composite Materials; 1999; Boca Raton, CRC Press:

26. Leissa, AW; Nordgren, RP. Vibration of shells; 1973; New York, Acoustical Society of America:

27. Thambiratnam, DP; Zhuge, Y. Axisymmetric free vibration analysis of conical shells. Eng Struct; 1993; 15, 2 pp. 83-89. [DOI: https://dx.doi.org/10.1016/0141-0296(93)90002-L]

28. Lam, KY; Hua, L. Influence of boundary conditions on the frequency characteristics of a rotating truncated circular conical shell. J Sound Vib; 1999; 223, 2 pp. 171-195. [DOI: https://dx.doi.org/10.1006/jsvi.1998.1432]

29. Reddy, JN. Mechanics of laminated composite plates and shells: theory and analysis; 2003; 2 Boca Raton, CRC Press: [DOI: https://dx.doi.org/10.1201/b12409]

30. Tian, L; Ye, T; Jin, G. Vibration analysis of combined conical-cylindrical shells based on the dynamic stiffness method. Thin Wall Struct; 2021; 159, [DOI: https://dx.doi.org/10.1016/j.tws.2020.107260] 107260.

31. Li, W. Free vibrations of beams with general boundary conditions. J Sound Vib; 2000; 237, 4 pp. 709-725. [DOI: https://dx.doi.org/10.1006/jsvi.2000.3150]

32. Li, W. Vibration analysis of rectangular plates with general elastic boundary supports. J Sound Vib; 2004; 273, 3 pp. 619-635. [DOI: https://dx.doi.org/10.1016/S0022-460X(03)00562-5]

33. Joines JA, Houck CR (1994) On the use of non-stationary penalty functions to solve nonlinear constrained optimization problems with GA's. First IEEE Conference on Evolutionary Computation, Florida, pp 579–584. https://doi.org/10.1109/ICEC.1994.349995

34. Shi Y, Eberhart RC (1999) Empirical study of particle swarm optimization. Congress on Evolutionary Computation-CEC99, Washington DC, pp 1945–1950. https://doi.org/10.1109/CEC.1999.785511

35. Eberhart RC, Yuhui S (2001) Tracking and optimizing dynamic systems with particle swarms. Congress on Evolutionary Computation, Seoul, pp. 94–100. https://doi.org/10.1109/CEC.2001.934376

36. Wang, L; Pan, J; Jiao, L. The immune algorithm. Acta Electron Sin; 2000; 28, 7 pp. 96-101.

37. Lee, YS; Kim, YW. Effect of boundary conditions on natural frequencies for rotating composite cylindrical shells with orthogonal stiffeners. Adv Eng Softw; 1999; 30, 9–11 pp. 649-655. [DOI: https://dx.doi.org/10.1016/S0965-9978(98)00115-X]

38. Zhao, XL; Liew, KM; Ng, TY. Structures, Vibrations of rotating cross-ply laminated circular cylindrical shells with stringer and ring stiffeners. Int J Solids Struct; 2002; 39, 2 pp. 529-545. [DOI: https://dx.doi.org/10.1016/S0020-7683(01)00194-9]

39. Song, Z; Su, C. Computation of Rayleigh damping coefficients for the seismic analysis of a hydro-powerhouse. Shock Vib; 2017; 2017, 2046345. [DOI: https://dx.doi.org/10.1155/2017/2046345]

40. Ting, W; Hui, W; Yang, Z; Xinghong, L. Digital mechanism modeling and vibration characteristic analysis for stator end winding of large turbo-generator. J Vibr Shock; 2023; 42, 7 pp. 143-153.