Stability Analysis of Hydrodynamic Mechanical Seals in Multifrequency Excitation

: The dynamic characteristics of the complex relationship among the sealing system, excitation, and response have a considerable impact on the operational reliability of hydrodynamic mechanical seals, which is a critical issue in the field of sealing theory and technology. Scholars at home and abroad have established dynamic models and calculated the displacement responses of dynamic and static rings in the time domain based on the force on these rings so that the response results can be used for system stability analysis. Neither are the excitation characteristics of cavitation load extracted, nor are the distance response and system leakage rate of the dynamic and static rings analyzed under coupled cavitation and random excitation. In this study, under different operating conditions of the hydrodynamic mechanical seal system, the liquid film evaporation load and seismic load are applied to study the frequency domain response of the distance between the dynamic and static rings and the system leakage rate. The following conclusions have been obtained: Assuming that the chamber pressure is 0.5 MPa and the spring specific pressure is 0.055 MPa, during stable operation, the distance between the moving and stationary rings at 1500 rpm~3000 rpm speeds is 1.12 μm~3.05 μm. For a specific spring pressure of 0.055 MPa, medium pressures of 0.2 MPa~1.0 MPa, and spindle speeds of 1500 rpm~3000 rpm, the excitation force is 30 N, and the frequency is 30 Hz, And the seismic load is assumed to be sinusoidal, the excitation force is 6 N, the fundamental frequency is 120 Hz, and the system leak rate is in 0.1 mL/min~1.3 mL/min. Under multi-frequency excitation coupling, the distance between the dynamic and static rings will decrease as the pressure of the medium in the sealing cavity increases, and this will increase with the increase in the rotating speed. The leakage rate of the system will increase with the increase in the rotating speed and the pressure of the medium, and the test value is largely consistent with the theoretical value.


Introduction
Fluid dynamic pressure-type liquid film seals have the characteristics of low leakage, little wear on the end faces of dynamic and static rings, and a long service life [1,2]. These are suitable for advanced sealing applications that are subjected to high-speed, high-temperature, and high-pressure conditions. Seals are widely applied and studied in aviation, aerospace, petrochemical, shipbuilding, and other industries. However, the operation of a hydrodynamic liquid film seal will be affected by fluid film cavitation, external impact load, and other excitations, and under the coupling effect of multiple-frequency excitations, the distance response of the dynamic and static rings will show steps, hysteresis, limit cycles, bifurcations, abrupt changes, and chaos, resulting in excessive system leakage [3].
To mitigate the aforementioned disadvantages of hydrodynamic liquid film seals, scholars at home and abroad have studied the system's leakage mechanism from the simplified excitation characteristics, film thickness response, and stability mechanisms. Xu [4] simplified the resultant force generated by the liquid film cavitation, the axial disturbance, and the transient operating conditions and expressed them as the axial impact excitation on the moving ring. The study found that the smaller the response of the moving ring, the longer the time taken to return to the stable state after being excited by external disturbances. Pang [5] considered the different pressure loads in the cavity resulting from the different spindle eccentricities as pulsating excitations on the dynamic ring and studied the influence of the seal excitation force on system stability. Salant [6] set the axial resultant force on the moving ring as the impact excitation due to the influence of the end groove and end-face deformation on the body film pressure and calculated the distance between the moving and stationary rings and the leakage of the system. Varney [7] studied the influence of the deviation in installation coaxiality on the system. Considering the impact excitation of resultant forces in all directions, he established a three-degree-of-freedom structural dynamic model of the dynamic and static rings, analyzed the distance response of the dynamic and static rings, and determined that the deviation in coaxiality during their installation would cause the dynamic and static rings to collide, thus causing the vibration of the sealing system. Falaleev [8] calculated the stiffness response of a liquid film using the finite volume method based on the influence of the liquid film cavitation and the thermal deformation of the seal ring and analyzed the reliability of the system operation using the pressure of the liquid film as the impact excitation. Xu [9] analyzed the influence of parameters such as the angle and height of the spiral groove on the liquid film pressure and deduced through simulation that the stability of hydrodynamic liquid film seals can be significantly improved by selecting an appropriate shaft angular velocity and height of the spiral groove. He [10] designed a new spiral groove structure that effectively improves the hydrodynamic performance of the liquid seal, resulting in a significant increase in dynamic stability. Wang [11] and Zhu [12] changed the influence of the cavitation excitation on the stability by changing the parameters of the seal end face. Su [13] increased the opening force of the seal interface by adding a conical diffuser on the end face of the dynamic ring, and this sealing system can quickly enter stable operation during the startup.
Zhang [14] studied the pressure pulsation of the liquid film flow field on the end face of mechanical seals. The results show that the pressure pulsation of the T-groove liquid film exists at different rotational speeds. Guo [15], Tian [16], Yang [17] et al. analyzed the root cause that affects the operation of the mechanical seal. However, none of the excitation studies involved the response of external excitation to the film's thickness. In conclusion, in the research and application of hydrodynamic liquid film seals, the cavitation load of the liquid film directly affects the sealing performance and will affect the hydrodynamic liquid film seals. Therefore, the excitation characteristics of the load need to be extracted to analyze the distance response of the dynamic and static rings, and the system leakage rate when the system is subjected to multiple coupled excitations, such as fluid film cavitation and seismic load.
The perturbation method was adopted to study the transient stiffness coefficient and the transient damping coefficient of the fluid-gas two-phase dynamic pressure seal based on Reynolds equation. The vibration characteristics of the system was analyzed, the excitation characteristic parameters of the fluid film cavitation load was extracted and the relationship between the distance of the dynamic and static rings and the rotational speed was obtained on the basis of previous work. In a further, the opening force was studied using a simulation software, and the clearance between the dynamic and static rings on steady operating and at different rotational speed was acquired. On this basis, the Volterra series is used to establish a mathematical model to calculate the distance response and the leakage rate of the dynamic and static rings of the hydrodynamic liquid film seal system under the coupling excitation of cavitation and seismic loads. This provides a theoretical basis for the research, development, use, and maintenance of sealing equipment for marine and aeronautical applications.

Physical Model
To analyze the relationship that connects the distance between the dynamic and static rings of hydrodynamic mechanical seals and the system leakage rate-as well as the operating conditions and physical parameters of the sealing components, considering the mass m of the moving ring, the fluid film stiffness Kq between the moving ring and the stationary ring Msp, fluid film damping Cq, spring stiffness K, and damping C of the Oring-a physical model of the excitation and film thickness response of the fluid dynamic pressure-type liquid film seal system was established, as shown in Figure 1. The cavitation excitation was set as x1, the random excitation was set as x2, and the film thickness when the mechanical seal system was under stable operation was △. The structure of the dynamic ring end face of the hydrodynamic mechanical seal involving liquid film evaporation is shown in Figure 2. The end face parameters of the movable ring are presented in Table 1. The working principle of the seal is that when the main shaft of the seal rotates clockwise, the lubricant forms a liquid wedge and exerts a dynamic pressure on the end face of the dynamic and static rings. It extrudes between the end faces of the dynamic and static rings, forming a liquid film on the end face to provide lubrication and cooling to the sealing end face. The end face of the dynamic ring contains spiral grooves, sealing dams, and sealing weirs.

Extraction of Excitation Characteristics of Liquid Film Vaporization Load
In the research and application of hydrodynamic mechanical seals, the phenomenon of liquid film vaporization directly affects the sealing performance, causing the seal leakage to exceed the specifications. Preventing such leakage requires an in-depth analysis of its mechanism. Therefore, based on the liquid-gas two-phase Reynolds equation, this study uses the perturbation method to study the stiffness coefficient and damping coefficient of the liquid-gas two-phase dynamic pressure seal and extracts the excitation characteristics of the liquid film evaporation load to analyze the response of the seal dynamic ring and the leakage of the system.

Control Equation of Liquid-Gas Two-Phase Membrane
As the fluid film on the end face of the hydrodynamic mechanical seal produces a liquid-gas two-phase fluid owing to the effect of the vaporization load of the liquid film, it is assumed that the volume fraction of the gas film in the liquid on the seal end face is c, the gas phase in the two-phase fluid on the end face is small, the gas is evenly distributed in the liquid, and the interaction between the gases is ignored: where μm is the equivalent viscosi t y of the liquid-gas two-phase fluid, μliq is the viscosity of the liquid, and c is the volume fraction of the gas film in the liquid of the seal face. The equivalent density of liquid-gas two-phase fluid is as follows: where ρm is the equivalent density of the liquid-gas two-phase fluid, ρliq is the liquid density, and ρgas is the gas phase density.
Equation (3) is obtained by introducing Equations (1) and (2) into the unsteady twodimensional cylindrical Reynolds equation for a compressible fluid with constant temperature and viscosity: where pm is the liquid-gas two-phase film pressure, v is the spindle speed of the sealing system, r is the radius, h is the fluid film thickness, and θ is the angle.
In the system model shown in Figure 1, to analyze the dynamic characteristics of the two phases, it is assumed that the steady-state liquid film thickness for the moving ring and the stationary ring is h0, and p0 is the steady-state position pressure of the fluid film. When the sealing system is subject to a small disturbance, the moving ring experiences axial vibration at the steady-state position, the vibration displacement is y(t), and the perturbation motion of the moving ring is defined as a simple harmonic motion, as shown in Equation (4): where ω is the frequency that the dynamic loop is perturbed.
The derivative of the above equation can be used to obtain the perturbation speed of the moving ring, as shown in Equation (5): Therefore, the transient thickness of the fluid film of the liquid-gas two-phase dynamic pressure seal is given with Equation (6): The transient pressure expression of the fluid film is given with Equation (7): Defining Formula (8) and (9): By making use of the equality between the real and imaginary parts of the complex numbers, the steady-state fluid film pressure p0 of the seal and the dynamic perturbation fluid film pressure pr and pi can be obtained. The specific derivation process is to substitute Equations (5)-(7) into Equation (3), ignoring the higher order terms, and obtaining the liquid-gas two-phase dynamic Reynolds equation, as shown in Equation (10): After separation from Equation (10), the steady-state and dynamic Reynolds equations of the two-phase fluid can be obtained. The steady-state equation is indicated as Equation (11), and the dynamic equation is indicated as Equation (12).
Equation (12) is derived from y(t) and simplified to (13): The complex variable in Equation (8) is introduced into Equation (13) to obtain Equation (14): After Equation (14) is separated, the pressure equations of the real and imaginary parts of the two-phase flow can be obtained, as shown in Equations (15) and (16), respectively: where v is the spindle speed of the sealing system, ω is the perturbation frequency of the moving ring, r is the radius, h0 is the thickness of the fluid film, θ is the angle, p0 is the steady-state seal fluid membrane pressure, py is the dynamic seal fluid membrane pressure, pr is the real part of the dynamic perturbation fluid membrane pressure, and pi is the imaginary part of the dynamic perturbation fluid membrane pressure. The opening force Equation (17) The dynamic stiffness and dynamic damping coefficient of the liquid film on the end face are given with Equation (18) and Equation (19), respectively.

Extraction of Excitation Characteristics of Load
Assuming that the excitation signal received by the system during the operation of the mechanical seal has finite energy, the axial nonlinear vibration equation of the moving ring can be simplified using the Duffing equation, as shown in Equation (20): where m is the mass of the moving ring, C is axial damping, K is axial stiffness, and b is nonlinear stiffness. When considering the influence of the liquid film evaporation load on the seal face, Equation (21) is obtained from the system balance: According to Equations (20) and (21), the excitation of the liquid film vaporization is given with Equation (22): After Equation (5) is introduced into Equation (22), it is known that Equation (22) results in an imaginary expression; therefore, the evaporation load on the liquid film can be defined as harmonic excitation, defined as Equation (23)

Dynamic and Static Ring Distance Responses at Different Speeds
It was assumed that the film's thickness represents the distance between the dynamic ring and the static ring. A three-dimensional modeling software UG was used to build the model, and the model was imported into the professional simulation software ICEM for grid division and dynamic simulation analysis in ANSYS. CFD calculations were performed on the opening force for different film thicknesses, as shown in Figure 3; at a film thickness of 1.4 μm, the total number of the grid units was 2,522,232. When the hydrodynamic mechanical seal is in operation, the opening force formed by the pressure of the fluid blocking the end face of the dynamic and static rings will change with the change in the spindle speed, forming a different balance with the closing force of the spring. Thus, the distance response of the dynamic and static rings will vary with the change in the spindle speed. When the pressure in the cavity was 0.5 MPa, the specific pressure of the spring was 0.055 MPa, the mass of the dynamic ring was 0.5 kg, the damping was 400 N·s/m, and the spring stiffness was 10,000 N/m; the relationship between the distance, speed, and opening force of the dynamic and static rings obtained with numerical simulation using Fluent software is shown in Figure 4. The results of this paper are consistent with reference [18]. Based on the above data, the relationship between the distance h0 between the movable and stationary rings at each speed and the opening force was fitted, as shown in Table 2. Table 2. Relationship between end-face distance and opening force at the same speed.

(25)
2500 Assuming that the chamber pressure is 0.5 MPa and the spring specific pressure is 0.055 MPa, the opening force of the data listed in Table 1 The distance between the moving and stationary rings at different speeds was calculated by introducing Equation (24)

Relationship between Excitation and Frequency Domain Response
When a non-contact mechanical seal operates and the system is subjected to external excitation x(t), such as the cavitation excitation of the fluid film or external impact, the relationship between the output response y(t) of the moving ring and the input excitation x(t) is the multiple convolution integral known as the Volterra series as shown in Equation (29): where hn (τ1,...,τn) is the nth-order Volterra time-varying kernel function.

Dynamic Loop Response in the Case of Multi-Frequency Coupling
We substitute Equation (32)  According to Equation (33), the output response spectrum Y(jω) is necessary to calculate the Volterra frequency domain kernels of each order. To facilitate the calculation, it was assumed that the seal operated stably. Before the excitation disturbance, the axial disturbance displacement and initial velocity of the moving ring were zero. Due to the collision between the moving ring and the stationary ring, and other reasons, the stiffness is assumed to be highly nonlinear, the mass of the moving ring is 1 kg, and the dynamic Equation is designated as Equation (34): Ky t by t x t dt dt yt y t t m dt We assume that Equation (29) is the solution of Equation (34). In Equation (29), we repeatedly use the derivative formula of the integral with parameter variables to find dy(t)/dt, d 2 y(t)/dt, y 3 (t), which can be substituted with Equation (34), compare the coefficients on both sides of the equation, obtain a set of kernels that determine each order  hn(τ1,...,τn), and then perform a Fourier transform on the linear differential equation to obtain the first fifth-order Volterra frequency domain kernel [19], as shown in Equations (35)-(39): )   3  1  2  3  4  4  5  5  1  5  2  1  2  3  1  2  3 3 ( , , ) ( ) ( ) ( , , )

System Leakage Rate
When the average system is running, we assume that the balance position of the dynamic ring and the end clearance of static rings are  , and the motion amplitude of the dynamic ring motion is given with Equation (40): The displacement of the moving ring is given with Equation (41): where p1 and p2 represent the pressure at the inner and outer diameter sides of the seal ring, respectively, and η indicates the viscosity of the sealing medium.

Numerical Example of Distance Response Analysis between Dynamic and Static Loops and System Leakage Rate under Multi-Frequency Excitation Coupling
In this section, based on the different speeds of the system, the parameters listed in Table 4 are selected for frequency domain analysis to obtain the distance response and leakage rate of the dynamic and static rings during the vibration balance. The sealing medium is water, and the viscosity of the sealing medium η = 0.9358 × 10 −3 Pa·s The pressure in the sealing chamber is 0.5 MPa, specific pressure of the spring is 0.055 MPa, excitation force is 30 N, and the frequency is 30−100 Hz. The seismic load is assumed to be sinusoidal; the excitation force is 6 N; the fundamental frequency is 10 Hz−120 Hz; and the harmonic quantity is taken as the cavitation excitation. The output spectrum of the free-moving end of the moving ring is shown in Figure 5. When the spindle speed is 1500 rpm, the vibration output spectrum of the dynamic ring was obtained as shown in Figure 6a. It can be seen from the figure that the displacement response of the active ring reaches 1.1 μm, and as excitation frequency range increased, the dynamic and static rings collided.
When the spindle speed is 2000 rpm, the vibration output spectrum of the dynamic ring was obtained as shown in Figure 6b. It can be seen from the figure that when the displacement response of the moving ring reaches 2.0 μm, the moving and stationary rings collide, and when the speed of the main shaft is 1500 rpm, the excitation frequency range of the moving and stationary rings with collision becomes narrower as the clearance distance between the moving and stationary rings becomes larger with the increase in the speed of the main shaft.
For a spindle speed of 2500 rpm, the vibration output spectrum of the dynamic ring was obtained as shown in Figure 6c. As shown in the figure, when the displacement response of the dynamic ring reaches 2.6 μm, the dynamic and static rings collide, and the excitation frequency band of the dynamic and static rings with collision becomes narrower than that when the spindle speed is 1500 rpm and 2000 rpm.
When the spindle speed is 3000 rpm, the vibration output spectrum of the dynamic ring was obtained as shown in Figure 6d. It can be seen from the figure that when the displacement response of the active ring reaches 3.1 μm, in the narrow range of the excitation frequency, the dynamic and static rings collide. This is because when the excitation force reaches a certain value, the displacement response of the dynamic ring remains unchanged, and the clearance distance between the dynamic and static rings increases with increasing speed. The clearance distance between the moving and stationary rings is incorporated into the leakage rate Formula (42) to obtain the leakage rate at different speeds, as shown in Table 5. 3.14 0.89

Result Analysis
For a specific spring pressure of 0.055 MPa; medium pressures of 0.2 MPa, 0.4 MPa, 0.6 MPa, 0.8 MPa, and 1.0 MPa; and spindle speeds of 1500 rpm, 2000 rpm, 2500 rpm, and 3000 rpm, the relationship between the system leakage rate and multi-frequency excitation is obtained via Equation (42), as shown in Figure 7. During system operation, eddy current sensors are used to measure the distance between the dynamic and static rings. After running for 1000 s, we measured the leakage rate using a beaker and divided it by the operating time to obtain the average leakage rate. The measured experimental values are shown in Figure 7. In the case of multi-frequency excitation coupling, the leakage rate of the system will increase with the increase in rotating speed and with the increase in medium pressure. The experimental value is consistent with the theoretical value as the increase in rotating speed is equivalent to the increase in system excitation, and the increase in medium pressure results in an increase in the fluid flow thrust. In summary, the theoretical values of the clearance distance and leakage rate of the dynamic and static rings of the system are consistent with the experimental values and are of the same order of magnitude, but the theoretical values are slightly smaller than the experimental values. According to the calculation, the error between the experimental value and the theoretical value is 15%. This is owing to the influence of the installation error of the moving and stationary rings, machining error of the moving and stationary rings, fluctuation of pressure in the cavity, and axial displacement of the moving ring caused by the fluctuation of spindle speed during the experiment. The experimental results of system leakage rate under multi-frequency excitation are given in Table A1.

Conclusions
1. In this study, based on the two-phase Reynolds equation, the transient stiffness coefficient of the fluid film and the transient damping coefficient of the fluid-gas twophase dynamic pressure seal were studied using the perturbation method, based on which the vibration characteristics of the system were analyzed, excitation characteristic parameters of the fluid film cavitation load were extracted, and the Volterra series was used to establish a mathematical model to calculate the distance response and leakage rate of the dynamic and static rings of the hydrodynamic liquid film seal system under the coupling excitation of cavitation and seismic loads. This provides a theoretical basis for the research, development, use, and maintenance of sealing equipment for marine and aeronautical applications.
2. The relationship between the opening force and rotating speed was numerically simulated using Fluent software (version 2022), and the disengagement speed under different cavity pressures was calculated. The distance response of the dynamic and static rings increased with the rotating speed under a certain cavity pressure. 3. With an increase in the rotating speed, under the same excitation conditions, the distance between the dynamic and static rings and the leakage rate of the system increased when the vibration was balanced. 4. To predict the clearance distance between the dynamic and static rings and system leakage rate, it was necessary to comprehensively examine the spindle speed, cavity pressure, spring specific pressure, and system load. Moreover, the results of this paper can be applied to the study of non-contact mechanical seals. Funding: This work was supported by the "Study on load characteristic parameter identification and theoretical basis of preventing leakage and collision of hydrodynamic liquid film seals" (KX202106101) and the National Natural Science Foundation of China's "theoretical basis research on zero leakage and long-term safe operation of diffuser self pumped hydrodynamic mechanical seal for nuclear main pump" (52075268).
Institutional Review Board Statement: Not applicable.

Informed Consent Statement: Not applicable.
Data Availability Statement: Data-sharing is not applicable to this article.

Conflicts of Interest:
The authors declare no conflict of interest. -Outer radius of sealing movable ring, mm ri -Inner radius of sealing movable ring, mm y(t) -Relative displacement of moving ring to balance position, mm μm -Equivalent viscosity of liquid gas two-phase fluid, Pa·s μliq -Viscosity of liquid, Pa·s ρm -Equivalent density of liquid gas two-phase fluid, g/cm³ ρliq -Liquid density, g/cm³ ω -Excitation frequency, HZ v -Spindle speed, rpm η -Viscosity of sealing medium water, Pa·s Appendix A