Numerical Investigation of Gas Binding Dynamics in Centrifugal Pumps Using LBM–LES

1 Introduction

Centrifugal pumps are widely used in petroleum, chemical, electric power, metallurgy, and water supply [1,2,3,4]. As the classic rotating machines, continuous high-speed operation of centrifugal pumps can lead to various faults in harsh environments. Now the common faults mainly include mechanical faults [5,6], electrical faults [7], and hydraulic faults [8,9]. The mechanical faults include a stuck shaft fault, a seal fault, bearing damage, a motor stator phase short circuit, blade damage, cavitation, and blockage.

At present, scholars use different methods for different faults. Liu et al. [10] employed a combination of local linear embedding and wavelet transform to analyze the pump vibration signals, achieving high accuracy in identifying bearing faults. Panda et al. [11] investigated the changes in vibration signals at the shell and base of the pump under cavitation and clogging conditions, and obtained various statistical indices. And they trained and experimented with these indicators using the support vector machine (SVM) method to enhance the accuracy of diagnosing faults. Li et al. [12] constructed a dataset for pump cavitation identification based on vibration signals and used a particle swarm optimization support vector machine algorithm to achieve the identification of cavitation states. The method exhibited high recognition. Gu et al. and Hamomd et al. [13,14] investigated the pump disc wear faults based on motor rotational speed fluctuations. They conducted order-ratio analysis of velocity signals, which effectively compensated for the non-smooth and non-linear features of conventional monitoring signals. Using artificial intelligence algorithms, Sunal et al. and Wang et al. [15,16] identified various faults such as bearing faults, impeller breakage, pump cavitation, and shaft misalignment. The results demonstrated that this method achieved superior fault identification compared to traditional methods.

In recent years, although research on monitoring and diagnosing various faults of centrifugal pumps tends to be mature, research on monitoring and diagnosing the GBF, a classical fault, remains relatively limited. During the operation of the centrifugal pump, unstable internal flow often leads to significant performance degradation and severe energy loss. Dehghan et al. [17] systematically studied the numerical and experimental optimization of the pump volute, emphasizing the key mechanisms of energy loss and strategies for improving hydraulic efficiency at the optimal efficiency point. Moreover, the two-phase flow phenomenon caused by phase change seriously disrupts the internal flow field. Dehghan et al. [18] proposed a new physical criterion for the determination of initial cavitation, opening up new methods in the field of energy conservation. The phenomenon of the gas binding will lead to the reduction of the pump’s flow, head and efficiency, thus reducing the performance of the pump, which even impacts the pump’s life and stable operation [19,20,21,22]. Therefore, it is of great theoretical value and engineering significance to study the GBF in the centrifugal pump. However, the research of the GBF is mainly carried out from two aspects: the experiments and the numerical simulation at present.

Firstly, in terms of the experiments, Mousmoulis et al. [23] and Cui et al. [24] explored the phenomenon of the GBF in the pump by high-speed photographic technology. Long et al. [25] captured the transient flow structure in a jet pump using a high-speed camera, revealing the evolution of the cavitation-induced gas plugging phenomenon. Although the cavitation evolution to GBF can be directly observed by high-speed photography, accurate information on the internal flow field cannot be obtained because of the presence of gas bubbles. Coupled with the special experimental requirements and high costs, it is currently impractical for widespread application. In terms of numerical simulations, Shi et al. [26] and Liu et al. [27] investigated the cavitation features of a gas craft fuel pump by means of numerical simulation, revealing the mechanisms behind the occurrence and development of cavitation bubbles that lead to clogging in the pump. Li et al. [28] studied the distribution features of bubbles in the pump under different β based on the Eulerian-Eulerian method and the bubble number density equation, and its results showed that the bubble diameter reached the maximum value at the highest of β, and thus bubbles were prone to aggregation under high β conditions. Refined numerical simulation can enhance the accuracy of the numerical simulation results by LBM-LES [29]. Guo et al. [30] used a combined LBM-LES method to explore various vortex structures and their developments in an open-type pump inlet, and verified the method’s feasibility and stability by experiments. Pessoa et al. [31] used the same method to study the complex flow patterns in the intake of a pump, and the results demonstrated that the method can simulate free surface vortices accurately. Barrios and Prado’s pioneering experimental work [32,33], through visualization techniques, clearly captured and defined typical flow patterns within pumps (such as bubble flow and gas-bubble flow) for the first time, revealing the correlation between flow pattern evolution and operating conditions. The systematic review by Lei [34] and colleagues adopts a broader perspective to examine and underscore the critical importance and challenges of understanding transient gas-liquid two-phase flow patterns throughout the entire flow channel under demanding conditions such as those in nuclear main pumps. These authoritative reviews and key studies collectively form the cornerstone of knowledge in this field, pointing toward a future direction of deepening our understanding of the mechanisms underlying manifold evolution dynamics through the integration of advanced observations and numerical simulations.

In summary, scholars predominantly focused on constant rotational speed conditions at present [35,36]. They emphasized flow fields but neglected the mechanical responses of the shaft system and electromagnetic motor field influences. Actually, the rotational speed of the pump is also affected by mechanical responses of the shaft system and electromagnetic from the motor. Therefore, it is necessary to carry out hydro-mechanical-electrical multi-physics field numerical simulation studies, considering the mechanical properties of the rotor shaft system and the electromagnetic field influence of the motor. At present, research on hydro-mechanical-electrical multi-physics field numerical simulation focuses on the flow features of liquids in fluid-mechanical systems [37], pressure and energy conversion processes [38], synergistic interactions between mechanical equipment and electrical systems [39], and intelligent decision support.

This paper aims to conduct multi-physics field simulations and experiments under the condition of centrifugal pump air binding failure, define the turning point of the air binding failure state, and focus on analyzing the pressure rise performance and changes in motor characteristic parameters of the centrifugal pump under different gas contents. As follows specifically, this paper takes the G-2.2 centrifugal pump as the research object, adopts the LBM-LES method to conduct multi-physics field numerical simulation of the pump, which considers the influence of the motor’s electromagnetic field and the mechanical properties of the rotor shaft system. The paper investigates GBF development processes under different β conditions and analyzes pump torque, motor rotational speed, current, and electromagnetic torque time domain features under different GBFs. It provides reliable theoretical guidance for diagnosing gas cluster clogging faults under gas-containing incoming flows.

2 Research Methods

2.1 The Main Parameters and Three-Dimensional Model of the Pump

Based on the main parameters of the G-2.2 centrifugal pump in Table 1, the paper selects all the flow channels from the inlet to the outlet of the pump as the calculation domain, and completes the three-dimensional model using UG software at a 1:1 scale. The model includes inlet pipe, impeller, volute and outlet pipe, as detailed in Fig. 1.

Figure 1: Water body three-dimensional modeling schematic.

2.2 Numerical Simulation Methods

2.2.1 Lattice Boltzmann Equation

Lattice Boltzmann equation (LBE) can be viewed as a special “discretization” form of continuous Boltzmann equation. The Eq. (1) is the evolution equation of the distribution function fi with discrete time for LBE, which adopt cubic lattice with 19 particle discrete velocity directions (D3Q19 model).

In Eq. (1), fi is the distribution function for particles with discrete velocity ei at position x and time t along the i-th direction of velocity ei, ∆t is the time step, fi(eq) is the corresponding local equilibrium distribution, and τ is the single relaxation time.

Discrete velocity ei of the D3Q19 model are as follows.

In Eq. (2), C = ∆x/∆t, ∆x is the lattice spacing.

In the paper, the equilibrium distribution function fi(eq) is denoted as:

In Eq. (3), Cs2 = C2/3 is the sound speed of system, and C is the velocity of lattice. And the weighting coefficient ωi, macroscopic velocity u, density ρ and pressure p are given as follows.

2.2.2 LBM-LES Method

The Lattice Boltzmann Method (LBM) [40,41] is a fluid simulation technique based on mesoscopic dynamics theory. Unlike traditional computational fluid dynamics (CFD), which directly solves the macroscopic Navier-Stokes equations, LBM reconstructs macroscopic fluid motion by tracking the evolution and collisions of particle distribution functions in discrete velocity space. This method possesses inherent advantages in handling complex multiphase flow interface dynamics and intricate geometric boundaries, as it automatically captures interface changes without requiring complex interface reconstruction algorithms. However, when simulating high Reynolds number turbulence, the LBM faces challenges of numerical instability if no turbulence model is introduced. To this end, this paper combines LBM with Large Eddy Simulation (LES) [42,43], known as the LBM-LES method [30,34,44], to simulate the influence of small eddies on large eddies through a sublattice scale model, thereby effectively capturing the transient characteristics of turbulence.

The LBM-LES method mainly introduces the idea of large eddy simulation into the LBM. As Re increases, the viscosity decreases and the maximum velocity gradient converges, while the overly dense mesh causes compressibility errors. Under this condition, the Smagorinsky model is introduced in LBM. And the equivalent turbulent viscosity ν in this model is defined as follows.

In Eq. (5), νo is the molecular viscosity and νt is the sub-grid vortex viscosity coefficient, which is calculated as follows.

In Eq. (6), Cx is Smagorinsky’s constant. | S ¯ | is the modulus of strain rate tensor and the ∆ is the filter width. The | S ¯ | is calculated as follows.

The introduction of LES changes the flow field, which is directly reflected in the relaxation time. In the LBM-LES method, the relaxation time is replaced by the equivalent relaxation time. The relationship between the equivalent viscosity coefficient ν and the equivalent relaxation time τe is as follows.

The equivalent relaxation time τe is replaced by the relaxation time τ in the evolution Eq. (1) of the LBM single relaxation model. And the evolution equation of the LBM-LES coupled method is obtained.

In this paper, the LBM framework employs two core models: the BGK model and the Shan-Chen two-phase model. The BGK collision operator describes the transfer of mass and momentum, effectively handling the exchange process between two fluids. The Shan-Chen pseudopotential model is adopted to maintain force balance and solve the two-phase interface problem. The Shan-Chen model effectively assesses the inter-particle interaction forces between the gas and liquid phases. Through this pseudo potential interaction, macroscopic phase separation and interface surface tension naturally occur without the need to track the phase boundary. This eliminates the reliance on empirical models based on sub-grid bubble dynamics, significantly improving the physical accuracy of predicting gas accumulation and channel blockage phenomena.

2.2.3 Verification of Resolution Irrelevance

In terms of conventional CFD simulations, meshing and irrelevance verification are the key factors to improve the numerical simulation quality and effectively control the computational costs. In order to balance the accuracy of the numerical simulation with the computational efficiency, the effect of gradually increasing the resolution on head of the pump in the numerical simulation results is analyzed.

The boundary conditions are as follows: the speed of the pump is 2900 r/min, and the pressure is the standard atmospheric pressure at the inlet. The flow rate is 15 m3/h at the outlet, and the medium is water and the temperature is 25°C. As shown in Fig. 2, the head curve is essentially smooth when the resolution is greater than 0.0005. Therefore, the global resolution is determined to be 0.0005, and the local resolution at the tongue of the spacer and at the inlet of the impeller channel is refined to 0.0001.

Figure 2: Resolution independence verification.

In order to further verify the reliability and independence of the large eddy simulation, this paper conducts a comparative analysis of three sub-grid scale turbulence models. The centrifugal pump gas phase and vorticity distributions of the three turbulence models at an air content of 10% are analyzed. In Fig. 3, the results show that the WALE model outperforms the other models in accurately depicting the fine-scale vortex structure near the blade tip of the impeller and accurately capturing the dynamic aggregation of the gas phase towards the impeller hub. It is also more effective in simulating the fluid dynamics behavior of gas-liquid two-phase flow.

Figure 3: Distribution of gas phase and vortex of three turbulence model centrifugal pumps at 10% gas content under low flow condition.

2.2.4 Performance Parameters of the Pump for Numerical Simulation

The β is a parameter used to describe the proportion of gas in a fluid. It is defined as the ratio of the gas phase flow to the mixture flow. Its mathematical expression is expressed as follows.

In Eq. (11), Qg, Ql denote the flow of the gas and liquid phases, respectively. Qgl is the flow of the mixture, m3/s.

The pressure rise, Pc, is sensitive to β, and its definition is the pressure difference between inlet and outlet.

In Eq. (12), Pout, Pin are pressure at inlet and outlet, respectively, kPa.

The most direct indicators of the pump performance are head and efficiency. The mixture’s head HTp can be calculated as follows.

In Eq. (13), ν o u t − l and ν i n − l are velocity of liquid phases at the inlet and outlet, respectively, m/s. ν o u t − g and ν i n − g are velocity of gas phases at the inlet and outlet, respectively, m/s.

On the other hand, when the pump conveys gas-liquid mixture, the efficiency of the pump can be calculated as follows.

In Eq. (14), M T p is the gas liquid two-phase torque, N·m. ω is the angular speed of rotation, rad/s.

2.3 Hydro-Mechanical-Electrical Multi-Physics Field Simulation and Monitoring Methods

2.3.1 Numerical Simulation Technology

The performance of the pump at different β is affected by a variety of factors, of which the hydraulic excitation force is just one. During the power transmission process, there is a dynamic interaction among the water torque, the mechanical shaft and the motor. To have a better understanding of these complex dynamic features, especially the changes of torque, rotational speed, and the current of motor, this paper adopts the hydro-mechanical-electrical multi-physical field numerical simulation method to calculate the operating parameters of the pump unit under the GBF. The technology is as illustrated in Fig. 4.

Figure 4: Technology of hydro-mechanical-electrical multi-physical field numerical simulation.

As water flows into the pump, it applies a force to the vanes and pump body, generating a hydraulic torque. This water torque is transmitted to the motor through the shaft system. Mechanical properties such as stiffness and torsional deflection of the rotor shaft system affect the efficiency and stability of water torque transmission. When the water torque is transferred to the motor, the electromagnetic resistance torque changes with the water torque load, thereby affecting the stability of the motor stator current and rotational speed. This paper adopts a hydro-mechanical-electrical multi-physical field numerical simulation approach, which helps accurately monitor operating parameters of the pump unit at different β and provides an analytical basis for GBF.

This paper adopts a one-way sequential coupling strategy of hydro-mechanical-electrical to simulate the fault propagation process. The fluid pressure moment extracted from the LBM-LES fluid domain serves as the transient external excitation for the mechanical domain, while the dynamic response of the rotor shaft system is controlled by the multi-degree-of-freedom torsional vibration equation, as shown in Eq. (15).

In Eq. (15), [ J ] , [ C ] and [ K ] are the inertia, damping, and torsional stiffness matrices of the shaft system, respectively. θ ¨ , θ ˙ and θ represent the angular displacement, angular velocity, and angular acceleration vectors. T L is the load torque, T e is the electromagnetic resistance torque, and T f represents the mechanical friction losses.

In this paper, according to the results of hydro-mechanical-electrical multi-physical field numerical simulation of the pump unit, a method of monitoring the GBF of the pump unit based on the feature parameters is proposed, which is shown as the flow chart in Fig. 5.

Figure 5: Flow chart for monitoring and diagnosing the GBF of the pump unit.

2.3.2 Motor Numerical Model Construction

The paper builds the motor model using SCILAB. Firstly, the asynchronous motor and its open-loop constant voltage-to-frequency ratio control system is built using the motor module as well as the circuit control module in SCILAB. Secondly, the load torque of the shaft system is taken as the input of the motor system, the current, rotational speed and electromagnetic torque are the outputs of the system, and the calculation step is set to 0.00001 s.

The electromagnetic conversion inside the motor is mathematically described using the d-q synchronous rotating coordinate system. The control equation for the coupling between the electromagnetic torque and the mechanical shaft is shown in Eq. (16).

In Eq. (16), p n is the number of pole pairs, L m and L r are the mutual and rotor inductances, ψ d r and ψ q r are the d-axis and q-axis rotor flux linkages, i d s and i q s are the d-axis and q-axis stator currents, respectively.

Since the LBM-LES fluid solver and the SCILAB electrical solver operate on entirely different time scales, data exchange between the physical domains is managed via discrete time-series arrays. A cubic spline interpolation algorithm is employed at the physical interfaces to resample and synchronize the transient torque data from the fluid domain into the highly refined 0.00001 s electrical time step, ensuring smooth and accurate data transition across the domains, the SCILAB simulation model is shown in Fig. 6.

Figure 6: Model of rotational speed open-loop constant voltage-frequency ratio control.

2.4 Experimental Methods

2.4.1 Experiment Platform Construction and Its Data Acquisition System

The pump used in this experiment is a G-2.2 single-stage and single-suction centrifugal pump. The schematic diagram of the experiment platform is shown in Fig. 7 and the physical drawing is shown in Fig. 8a.

Figure 7: Schematic diagram of the experiment platform.

The open cycle experiment platform mainly includes components such as 1-water storage tank, 2-gas liquid two-phase mixing device, 3-inlet valves, 4-inlet pressure sensors, 5-centrifugal pump, 6-outlet pressure sensors, 7-gas flowmeters, 8-vortex pump, 9-turbine flowmeters, 10-outlet valves. In order to accurately simulate the operating of the pump unit under different β, a micro-bubble generator is used in the experiment as shown in Fig. 8b. According to Stokes’ Law, an ordinary-sized bubble will rise quickly to the surface in water and then burst, and its existence will be short. The rise velocity of micro-bubbles is significantly reduced. Therefore, micro-bubbles with higher stability are chosen for better observation of bubbles.

The data acquisition system of the experiment includes operating parameter acquisition module of the pump and motor. The acquisition module of the pump adopts the TP Automation pump experiment system, which is specially used to measure pump performance parameters such as flow, inlet and outlet pressures, head, and so on. The motor operating parameters are taken using a LabVIEW-based NI system, which consists of a variety of sensors, data transmission lines, NI acquisition boards, and LabVIEW acquisition programs. Instruments used in the system include current sensors, rotary encoder couplings, turbine flow meters, and additional inlet and outlet pressure sensors, with data transfer between these devices through the NI-USB6211 acquisition board, the current sensor employs an industrial-grade closed-loop current sensor for current signal acquisition, utilizing a multi-point zero-flux technology system with an accuracy level of 0.02. Inlet and outlet pressures are measured using capacitive absolute pressure transmitters, with an accuracy level of 0.075. Meanwhile, the LabVIEW is responsible for managing the data acquisition process, performing data processing and signal conversion work to ensure data integrity and accuracy.

Figure 8: Experiment platform and microbubble generation device: (a) Physical drawing of experiment platform, (b) Microbubble generator and micro-nano bubbles.

The specific experimental procedures are as follows:

• (1)Calibrate all measuring devices.
• (2)Install the measuring equipment and data acquisition lines, ensure that the test circuit meets the testing conditions, and adjust each measuring device to the optimal state.
• (3)Start the centrifugal pump, check whether the test system and data acquisition system are operating normally, and observe whether the waveform of the collected signals is stable.
• (4)By adjusting the gas rotor flowmeter connected to the high-pressure vortex pump, change the gas content of the inlet at the water intake, record the pump inlet and outlet pressures, flow rates, and gas content rates under different gas contents.
• (5)After determining the appropriate gas content range, record the motor output current, voltage, and rotational speed data under each working condition.
• (6)Gradually increase the gas content at the pump inlet, record the motor output current, voltage, and rotational speed data at each stage.
• (7)After each test, open the valve on the pressure measurement tube for venting and wait for the residual bubbles in the water storage tank to overflow.
• (8)Repeat the test, adjust the valve to conduct data collection under different working conditions, and record and analyze the data.

Conclude the test, summarize the variation patterns of the motor output current, voltage, and rotational speed under different gas content conditions, so as to facilitate subsequent analysis and research.

2.4.2 Uncertainty Analyses of Experimental Data

In order to ensure the accuracy of the test results, the experimental platform and data were analyzed through uncertainty indicators. Experimental uncertainties include Type A uncertainty, Type B uncertainty, and combined uncertainty. They are respectively determined by the equations as follows:

In Eq. (17), xi is the measured value, x ¯ is the average value, tn−1 is the confidence coefficient, taking 1.96, and n is the number of groups, fm is the accuracy of the sensor, fc is the accuracy of the acquisition device, k is the inclusion factor, taking 1.96, Er is the class A uncertainty, Es is the class B uncertainty, and E is the integrated uncertainty.

Under the same operating conditions, the flow rate, head, voltage, speed, and current signals were repeatedly collected 20 times to calculate the class A uncertainty. Class B uncertainty is obtained from the parameters of each sensor. Finally, the integrated uncertainty of the system is calculated from the class A and class B uncertainties. The uncertainty calculation results of the test data are shown in Table 2.

During the signal acquisition process, sensors play a crucial role, and their accuracy is often affected by interference. The actual signals are often disturbed by various noises, and the signals themselves will exhibit certain fluctuations. To ensure the accuracy and reliability of the signal sampling results, the most suitable signal time scale should be selected, and the signal quality can be evaluated and verified through transient uncertainty analysis. The types of signals collected in the experiment include voltage, current, and rotational speed signals. The influence of sampling frequency and sampling duration on transient uncertainty is analyzed. The transient uncertainty is as shown in Eq. (18):

In Eq. (18), x p i represents the peak value of the signal collected for the i-th time, and x p ¯ represents the average value of the effective values of the signals collected for n times. The calculation of the transient uncertainty of various signals is shown in Fig. 9.

Figure 9: Transient uncertainty values for various types of signals: (a) Transient uncertainty value of current signal, (b) Voltage signal transient uncertainty value, (c) Transient uncertainty value of rotational speed signal.

As shown in Fig. 9, the transient uncertainty values of the current and voltage signals start to decrease gradually after the sampling time exceeds two seconds. The transient uncertainty value of the rotational speed signal reaches its minimum between 1.5 and 2 s.

2.4.3 Data Verification

Under the condition of no gas, the external characteristic curve obtained from the centrifugal pump experiment and numerical simulation is shown in Fig. 10a, in which Q/Qd is the ratio of the actual flow rate to the rated flow rate. From the comparison of the external characteristic curves, it can be seen that the performance curves obtained from the experiment and the simulation are in a better match, and the error between them is small. In most cases, the head and efficiency obtained by numerical calculation are slightly higher than those obtained by experiment, because the friction loss and leakage loss of the disc are not considered in the numerical simulation calculation. In order to verify whether the flow field, torque and motor signals under the two-phase GBF condition have been reliably predicted, Fig. 10b presents the comparison of the centrifugal pump head and the variation curves of different gas contents between the experiment and the numerical simulation. Both decrease with the increase of gas content, and the change patterns are similar, the maximum relative error between the simulation results and the experimental results is 6.5%, which fully validates the accuracy of the numerical model under the gas-liquid two-phase condition.

Figure 10: Single-phase and multi-phase tests and simulation verification: (a) External characteristic curve, (b) Curve of head as a function of β.

3 Results and Discussion

3.1 Gas Phase Distribution and Vortex Distribution Law under the GBF

Fig. 11 shows the distribution of the gas phase within the impeller at β values of 0%, 3%, 6%, 9%, 12%, and 15%, under the operating condition of 0.25 Qd. When the β is about 3%, bubbles are mainly concentrated at suction surface and the impeller outlet. When the β is about 6%, bubbles gather near the impeller hub, and some bubbles move towards the impeller outlet, where gas-liquid separation occurs. With the β increasing, bubbles in the impeller hub continue to gather and become large gas clusters [45]. This is because when the gas-containing liquid enters the pump, the rapid rotation of the impeller generates a low-pressure region on the suction surfaces, where micro-bubbles can easily nucleate and aggregate. And the gas clusters are not easily dispersed by fluid shear. In addition, the forces on large gas clusters may be in an equilibrium, thus the clusters can aggregate together. When the volume of gas clusters increases sufficiently to block the flow of fluid, the phenomenon of flow channel blockage occurs. Gas clusters clogging the flow channel will lead to a decrease in flow, pressure, and efficiency. This situation will seriously affect the normal operation of the pump, known as the “gas binding” phenomenon. The analysis of gas phase volume fraction distribution at different β can show that on the one hand, due to the kinetic differences between the liquid and gas, gas will mainly aggregate in the inlet and blade suction surface. On the other hand, the larger the β is, the easier it is to produce “gas binding” phenomenon.

Figure 11: Cloud diagram of gas phase volume fraction distribution in the pump at 0.25 Q_d.

To further explore the generation of GBF and the flow features of the fluid, the vortex and velocity distributions in the pump are compared under different the β when the flow is 0.25 Qd, as shown in Fig. 12. It can be seen that under the small flow condition, with the increase of the β, the vortex is gradually formed at the impeller outlet and spreads to the pressure surface. At the same time, the flow velocity and the current increase accordingly. It is because the flow channels get narrow due to the increased concentration and diffusion range of the gas phase, which leads to a faster flow of liquid phase. In the region where the gas phase aggregates, the vortices are mainly concentrated at the vortex core location. As the concentration of the gas phase rises to a certain level, bubbles will form clusters and occupy the flow channel, resulting in a reduced liquid flow in that region. In the gas-liquid two-phase separation region, a secondary flow phenomenon similar to single-phase flow occurs due to the low concentration of the gas phase.

Figure 12: Cloud diagram of the pump vortex and velocity distribution at 0.25 Q_d: (a) Vortex distribution cloud map, (b) Velocity distribution cloud map.

3.2 Parameter Monitoring and Shaft System Response Analysis under the GBF Conditions

3.2.1 Variation Laws of Performance Parameters at Different β

Fig. 13 and Fig. 14 show when the flow is ranging from 4 m3/h to 10 m3/h (Rated flow Qd is 15 m3/h), and the rated speed is respectively nd, 0.8nd, 0.6nd, variation curves of head and efficiency with different β. It can be seen that under different flow, speed, head and efficiency will gradually decrease with the increase of the β. This phenomenon primarily stems from two key physical mechanisms: First, the presence of gas directly reduces the average density of the gas-liquid mixture. Since pressure rise theoretically scales proportionally with fluid density, a lower-density gas-liquid mixture naturally achieves a reduced pressure rise when the impeller performs the same amount of work on the fluid. Second, and more significantly, the energy dissipation mechanism triggered by gas plays a crucial role: the deformation, oscillation, coalescence, and rupture of bubbles under high-speed shear forces within the impeller, coupled with the reduction in effective flow area and flow separation caused by uneven gas distribution within the flow path, substantially increase turbulent dissipation and impact losses. This converts a greater portion of the shaft power into useless thermal energy rather than fluid pressure energy. Fig. 15 shows the average rate of change curves of head and efficiency with flow and speed. The average rate of change of head and efficiency, the slope of the straight line that is used to fit the curve, are increased as flow elevated and speed. This indicates that at large flow and high speed, the head and efficiency are affected by the β to a greater extent.

Figure 13: Curve of pump operating parameters with β under different flow: (a) H T p curve under different flow conditions, (b) η T p curve under different flow conditions.

Figure 14: Curve of pump operating parameters with β under different speed: (a) H T p curve under different speed conditions, (b) η T p curve under different speed conditions.

Figure 15: The average rate of change of head and efficiency: (a) Curve under different flow conditions, (b) Curve under different speed conditions.

3.2.2 Variation of Impeller Torque at Different β

At high rotational speeds, the frictional torque and centrifugal torque cannot be ignored. To ensure the physical rigor of dynamic coupling, the total load torque T L acting on the motor shaft is decomposed mathematically into a constant average component T L 0 and a periodic fluctuating component T L a s c , as shown in Eq. (19).

In Eq. (19), T L 0 and T e 0 represent the steady-state mean values, T L a s c and T e a s c denote the periodic dynamic fluctuations induced by the pump’s internal flow instabilities.

At the stable operating condition point, these components remain relatively stable. For the GBF research aimed at extracting dynamic characteristics from the fluctuating signals, the constant baseline shift caused by the friction component does not affect the research results, and thus is incorporated into the constant baseline load T L 0 . Among several torques, water torque occupies a larger proportion, so this paper effectively replaces the role of friction and centrifugal torque and focuses on the water torque. The pressure, velocity of fluid, and torque on both the pressure and suction surfaces of the blades throughout one-fifth of a cycle are chosen for analysis, as shown in Fig. 16.

It can be seen that, the velocity of fluid, pressure and torque at the pressure surface of the blades follow the same change trend. And the change trend of velocity of fluid and torque on the suction surface is opposite, while it is the same for pressure and torque. The faster the rotational speed of the pump to promote the fluid is greater, that is, the greater the resulting hydraulic torque. At the same time, the greater centrifugal force fluid is subjected to, the faster the fluid will flow. Due to the Bernoulli effect, the pressure is small in the place where the velocity of fluid is fast, which leads to the opposite trend of the velocity and pressure change at the suction surface. Gas-liquid two-phase flow leads to unsteady distributions of pressure and velocity fields within the impeller, and the increase of the β generates vortex turbulence, which further affects the pressure and velocity distributions. Dynamic changes in the pressure and velocity fields affect the force on blades and ultimately the hydraulic torque.

Figure 16: Time-domain plot of pressure velocity torque variation at various locations of blades: (a) Time-domain diagram of the pressure face, (b) Time-domain diagram of the suction surface.

To save computing resources and more accurately reflect the situation, the fluctuation of the sampled signal is relatively small. The sampling rate of the current and voltage signals is set to 20,000 Hz, and the sampling time is uniformly set to 2 s. This is aimed at optimizing the computing efficiency while retaining the ability to accurately describe the characteristics of the signals. After the signals stabilize, 0.5 cycles are selected for analysis. The Hanning window is used as the window function to reduce spectral leakage to a certain extent.

Fig. 17 shows the time-domain and frequency-domain characteristics of the impeller’s hydraulic torque under varying gas content ratios β. The presence of the gas phase not only alters the fluid’s equivalent density and viscosity, but more critically, introduces compressibility and phase-phase interactions. These factors profoundly influence the system’s flow state and energy transfer processes, ultimately modulating the impeller’s hydraulic torque response. As β increases, the hydraulic torque of the impeller exhibits a gradually rising trend. The primary mechanism involves gas phase aggregation reducing the effective flow area within the passages, increasing local flow velocities. Concurrently, the dynamic behavior of bubbles—such as deformation, oscillation, and coalescence—significantly amplifies the pulsation intensity of the pressure field. The combined effect of these factors elevates both the time-domain mean torque and the amplitude of torque fluctuations.

Figure 17: Hydrodynamic torque time domain and frequency domain plots at different β: (a) Time domain diagram at different β, (b) Frequency domain diagram at different β.

In the frequency domain, a distinct energy redistribution phenomenon can be observed: as the β value increases, the amplitude at the axial frequency (impeller rotational frequency) gradually decays, while the amplitude of low-frequency harmonic components significantly increases. This phenomenon originates from the complex energy transfer and dissipation mechanisms within the two-phase flow system. The primary mechanism behind the low-frequency harmonic amplification is that, as β increases, gas accumulations within the impeller—forming gas pockets or large bubbles—exhibit dynamic behaviors (such as shedding, coalescing, and collapsing) that generate pronounced low-frequency pressure pulsations. This phenomenon aligns closely with Hallez’s [46] classic conclusion that “gas phase introduction excites strong low-frequency excitation.” These low-frequency pulsations, dominated by gas phase aggregation dynamics, act upon the blades, directly causing a significant rise in low-frequency harmonic amplitudes within the hydraulic torque spectrum. Concurrently, the reduction in shaft frequency amplitude can be attributed to energy redistribution across frequency bands. The presence of the gas-liquid interface substantially enhances the system’s effective damping and nonlinearity. Processes such as bubble deformation and oscillation dissipate and scatter energy within the flow field. Concurrently, increased β values induce unstable flow phenomena—such as intensified vortex turbulence—generating turbulent kinetic energy across a broader frequency spectrum. Collectively, these effects shift system energy away from the regular impeller rotational frequency (shaft frequency) toward broader low-frequency bands. This “energy redistribution” effect aligns with Monte’s [47] findings that in multiphase flow systems, “the dominant frequency energy decays due to flow instability, while energy spreads toward lower frequencies and broader bandwidths.” Consequently, in the frequency domain, this manifests as a reduction in axial frequency amplitude.

3.2.3 Shaft System Response Features of the Pump Unit under Different β

Pressure pulse in the pump is closely related to changes in the Pc. When the Pc change is small, it indicates that the pressure difference between the inlet and outlet remains relatively stable, and the pressure pulse is small at this time. This is because the small fluctuation of the Pc indicates that the transport of fluid is more stable. During the flow of liquid, there is no obvious vortex, vortex and separation phenomenon, and no violent oscillation and disturbance. Conversely, when the Pc varies significantly, there are more intense vortex flow, swirling and separation during fluid transport, resulting in sharp fluctuations in fluid pressure. Therefore, the pressure pulse will be larger. In addition, the fluctuation of the Pc causes a non-uniform distribution of fluid flow and pressure within the pump, and this non-uniform distribution further exacerbates pressure pulse.

When pressure fluctuates due to pump faults, the torque on the impeller shaft changes. Changes of torque are transmitted through the rotor shaft system to the motor. Changes in the motor torque load are reflected in the electromagnetic torque signal, the current signal and rotational speed signal of motor through the electromagnetic coupling between the stator and rotor. During the transmission of torque, the dynamic response of the shaft system has an effect on the torque signal. Therefore, when a non-stationary analysis of the torque of the pump is performed, the dynamic response of the shaft system needs to be analyzed.

According to the structural drawings of the pump, the shaft system dynamics model is constructed on the Workbench soft, and the transient dynamics is analyzed by the finite element method. The hydraulic torque calculated by the LBM-LES method is used as an input to the shaft system model, as shown in Fig. 18a. The rotor shaft system of the experiment platform is rigid and its material chosen is #45 steel, the torque response of the terminal output of the shaft system in Workbench is shown in Fig. 18b. It can be seen that the input and output have a linear correlation tendency under the rigid rotor axis system.

Figure 18: Pump hydraulic torque input and shaft end torque output time domain plot: (a) Hydraulic torque input Time domain diagram, (b) Shaft end torque output Time domain diagram.

3.3 Analysis of Motor Operation Parameters at Different β

The output signals of the SCILAB three-phase asynchronous AC motors show a high degree of similarity at different β and need to analyze time-frequency domain feature indexes. The commonly used feature parameters in time-frequency domain are shown in Table 3, where x(i) is the analogue signal value and P(f) is the power spectrum of the signal [48,49,50].

For dynamic feature extraction and early fault detection, the specific effects of gas volume fraction on dynamic signals are strictly separated, and any changes in motor current, speed, and electromagnetic torque are attributed to the development of two-phase flow anomalies. The features indexes are extracted from the current, rotational speed and electromagnetic torque signals of the pump unit at 4 m3/h and speed nd condition, and the feature indexes that are sensitive to the change of β are selected as shown in Fig. 19.

Figure 19: Feature indexes of the output signal of motor at different β: (a) RMSA of electromagnetic torque at different β, (b) Impulse factor of electromagnetic torque at different β, (c) Waveform factor of electromagnetic torque at different β, (d) RMSA of rotational speed at different β, (e) Kurtosis of rotational speed at different β, (f) Frequcenty variance of current at different β, (g) RMS of current at different β, (h) THD of current at different β, (i) SINAD of current at different β.

It can be seen that the RMSA of electromagnetic torque, kurtosis factor of rotational speed, frequency variance of the current, SINAD of the current, and RMS of the current decrease gradually with the increase of the β. And the pulse factor and waveform factor of electromagnetic torque, RMSA of rotational speed, and THD of the current increase gradually with the increase of the β.

3.4 Experiment Verification of Signal Variation Law under GBF Conditions

Fig. 20 shows the Pc with the different β. It can be seen that the Pc calculated by numerical simulation has the same variation rule as that in the experiment, and the error is small. This shows the consistency between experiment and simulation. The Pc decreases with the increase of the β. Combined with the results of the analysis in Fig. 10, it can be seen that bubbles will block the flow channels causing a reduction in the Pc under GBF conditions, so you can judge whether the GBF occurs by the change in the Pc.

Figure 20: The P_c with different β changes law schematic diagram.

When the GBF does not occur, bubbles do not gather at the impeller hub, and the gas-liquid separation phenomenon is not yet obvious. At the same time, bubble flow or cohesion bubble flow state occurs in the impeller, and the Pc does not descend rapidly. With the increase of the β, bubbles begin to gather and move to the impeller hub, blocking flow channels and then the Pc decreases rapidly. And the stage is the early state of GBF. As the β increases further, the variation of the Pc tends to slow down and close to the lowest point, a large number of bubbles gathered in the impeller hub, at this time the pump for the serious GBF.

Fig. 21 shows the feature indexes extracted from the current and the rotational speed signals in the experiment and the numerical simulation. It can be seen that with the increase of the β, the RMSA of the rotational speed signal, SINAD and THD of the current show a certain law of change. This is because when the β increases, the presence of gas decreases the density of the gas-liquid mixture and enhances its compressibility. As a result, the inertial load on the pump diminishes, leading to an increase in the motor rotational speed. And there are two key points in curves, which can be used to judge the GBF. When the β is less than 6%, the RMSA of the rotational speed signal increases steadily from about This phenomenon primarily 70 r/min to about 3010 r/min. When the β is 6%, the feature index curve of the RMSA of the rotational speed begins to rise rapidly, and this is the first turning point. Combined with Fig. 20, it can be seen that when the β is 6%, there is a transition from the non-gas binding stage to the early stage of gas binding. The RMSA of the rotational speed tends to slow down when the β reaches about 12%, the second turning point, which is the stage of progression from the early stage of gas binding to severe gas binding. The variation of the motor rotational speed under the multi-physics field numerical simulation is monotonic in the range of about 2950.8–2999.4 r/min and the variation law of the RMS of the motor rotational speed is similar to that under experiment. Secondly, comparing the THD of the stator current signals under experiment and numerical simulation, respectively, it can be found that the variation law is similar to that of the RMSA of rotational speed, and the turning point of the GBF state occurs more obviously in the current signals.

Figure 21: Variation curves of feature indexes of the pump unit with increasing β, (a) RMSA of rotational speed with increasing β, (b) THD of current with increasing β.

The analysis of feature index curves, which extracted from motor operating paraments, along with the Pc shows that the turning point of the curve variation are used as the basis for judging the GBF. And the operating state of the pump unit is divided into three stages. No gas binding stage which is before the first turning point of change of the curve, and the β reaches about 6%. Gas-binding initial stage which is between the first and second turning points of change of the curve, and the β ranges from about 6 to 12%. Gas-binding severe stage which is after the second turning point of change of the curve, and the β reaches about 12%.

4 Conclusions

In this paper, a multi-physics field numerical simulation of a centrifugal pump unit under GBF is carried out based on the LBM-LES method, taking into account the effects of the electromagnetic field of the motor and the dynamic response of the shaft system. This paper has obtained the gas-phase distribution laws under different β conditions, and at the same time, the feature extraction and analysis of the motor signals under different β are realized. The main conclusions are as follows.

• 1.When the β is around 3%, bubbles are mainly concentrated on the suction surface of the blades and at the impeller outlet. When the β is about 6%, the bubbles converge at the impeller hub, and some bubbles move to the impeller outlet. Meanwhile, the gas-liquid separation phenomenon occurs. As the β increases further, bubbles continue to accumulate at the impeller hub and form large gas clusters that clog the flow channels. When the pump flow is small, with the β increasing, vortices gradually form at the impeller outlet and spread to the pressure surfaces of the impeller.
• 2.As the β increases, the shaft frequency amplitude of the hydraulic torque gradually decreases and the low frequency harmonic amplitude gradually increases. The feature indexes vary monotonically with the β at small flow, and the feature indexes show large fluctuations at large flow. As the β increases, the head, efficiency and the Pc decrease. And the head and efficiency are affected by the β to a greater extent under the large flow and high-speed conditions.
• 3.Based on the β, the Pc and feature index can be divided into three stages. In the first stage, the β is less than 6%, and the Pc and feature indexes change steadily, which is no gas binding stage. In the second stage, the β is at 6%–12%, and the Pc and feature index decrease rapidly, which is the stage when gas binding occurs. In the third stage, the β exceed 12%, and the Pc and feature index decline slows down, which is the severe gas binding stage.

However, this paper only deals with a single GBF in the case of an incoming gas-containing stream, and the actual situation is often more complex and may involve cavitation as well as gas-liquid phase transitions. Therefore, the mechanisms and laws when multiple faults are concurrent can be further explored in the future.