Construction of Solid-Liquid Two-Phase Flow and Wear Rate Prediction Model in Multiphase Pump Based on Mixture Model-Discrete Phase Model Combination Method

Abstract

Blade wear is the critical problem in the operation of multiphase pump. This paper presents a numerical study of the multiphase flow of multiphase pump. The trajectory of particles in the pump is calculated by the discrete phase model. Then, the simulation results are compared with the model test results of the pump to verify the correctness of the simulation method. The results show that the particles in the impeller domain are mainly near the hub, and the particles in the diffuser domain form a agglomerated area in the middle of the flow channel. The average wear rate of the impeller is more affected by the particle size than that of the diffuser. The maximum wear rate of blade surface increases first and then decreases with the increase of particle size. According to the wear data under different particle sizes, the regression model between particle size and wear rate is fitted to predict the wear of mixed transport pump in actual operation. The research results have important reference value for the prediction of the wear performance of the multiphase pump.

1. Introduction

In remote desert and offshore oilfields, volumetric pumps are currently the primary method for extraction [1]. The fluids produced in these fields often contain a high concentration of particles [2], which can lead to pump wear and significantly reduce its service life [3]. As a result, there is a growing trend towards the use of multiphase pumps with impeller structures in oilfields. These pumps effectively mitigate wear compared to traditional volumetric multiphase pumps. However, further research is needed to understand the solid-fluid two-phase flow and wear characteristics within the impeller of multiphase pumps.

The investigation of two-phase flow, specifically solid-liquid interactions, has consistently been a focal point within the realm of fluid mechanics. Scholars conduct experiments on their research subjects to obtain physical results of solid-liquid flow for analysis [4,5,6]. However, for the complex solid-liquid two-phase flow test, the design process and the experiment process need to consume huge manpower and material resources. In this regard, computational fluid dynamics is more applied to the research process. There are two prevalent approaches for calculating two-phase flow: the Euler-Euler method and the Euler-Lagrangian method. The Euler-Euler method treats the motion of solid particles as fluid flow and is suitable for high particle concentrations [7,8]. Scholarly research has revealed that the Euler-Euler method is capable of simulating specific two-phase flow scenarios involving solid and liquid [9,10], but it can also lead to significant errors in other research conditions [11,12,13]. The Euler-Lagrangian method achieves the tracking by integrating the particle motion trajectory in the whole system [14]. Scholars have used this method to study the particle motion and found that the initial velocity of the particle [15], forces on particles [16], the size and shape of the particle [17], the vortex in the flow field and the secondary flow [18] all have important effects on the particle motion.

Whether the consideration of particle velocity, particle force and particle size is only one aspect of solid-liquid two-phase flow calculation, the calculation of particle trajectory is also based on the calculation of flow field. When oil and water flow together, it is divided into oil-in-water fluid and water-in-oil fluid [19]. Oil and water two phase medium in the emulsified form, belong to the micro mutual mixing, need additional catalyst to separate [20,21]. When the two substances are fully mixed, it is not necessary to consider the interface between the media, and it is suitable to use the mixture model to solve the problem [22].

In the solid-liquid two-phase flow, the surface of the flow passage in the multiphase pump will experience significant erosion due to the impact of solid particles carried by the fluid [23]. Scholars have carried out a large number of wear tests, quantitative research and qualitative analysis of the cause [24,25] and location of wear [26,27]. The above publication of a large number of wear test data provides a basis for the simulation of wear phenomena in computational fluid dynamics. Researchers began to analyze the simulation calculation of wear, including: particle shape [28,29,30], particle concentration [31,32], particle size [33,34] and other effects [35,36] on wear. It is evident that there are well-developed computational methods for fluid flow, particle trajectories, and wear in the field of solid-liquid two-phase flow, which are highly accurate.

Based on the research methods employed in the aforementioned literature for two-phase flow studies and taking into account the specific engineering conditions addressed in this paper, the mixed model and discrete phase model (DPM) have been selected to analyze the flow behavior of oil, water, and particles. The trajectory of particle inside the multiphase pump and the distribution of particle in two-dimensional and three-dimensional space are revealed. Based on the Oka model, the wear distribution on the surface of impeller and diffuser and the average wear rate and maximum wear rate of impeller and diffuser in the pump with particle diameter are explored. Finally, the polynomial mathematical regression prediction model of average wear rate and maximum wear rate is constructed with particle diameter as variable.

2. Physical Model of Multiphase Pump

Figure 1 is a vane-type multiphase pump independently designed by Xihua University, which has been applied to Fengcheng Oilfield in China. Since the pump is a multi-stage type, each stage has the same structure. When the multiphase transport pump is directly selected for calculation, the flow of each stage is similar. In order to save the computational cost, the first stage of the pump is selected as the research object in this paper. The main parameters of impeller and diffuser are shown in

Table 1. Parameters of multiphase pump.

Design Parameter	Value
n	2980 (r/min)
Q	200 (m3/h)
Hydraulic Head	30 (m)
Impeller diameter	234 (mm)
Impeller blades number	3
Diffuser blades number	11

. The design flow is Q = 200 m³/h, and the speed is n = 2980 r/min. Using Unigraphics NX 10.0 software, the inlet and outlet are extended to ensure the full development of the inlet and outlet flow of the booster unit. The whole fluid domain model is simplified into four parts: inlet section, impeller, diffuser and outlet section. The specific three-dimensional calculation domain is shown in Figure 1. In the figure, the impeller and diffuser are made of stainless steel (04Cr13Ni5Mo) with stable mechanical and chemical properties and a density of 7.7 g/cm³.

3. Grid Division and Boundary Condition Setting

3.1. Grid Division and Independence Verification

Mesh quality has an important impact on the accuracy and efficiency of CFD calculations. Mesh generation is the most critical step in the preprocessing of CFD 16.0 software. In this paper, the Euler-Lagrange method is used to solve the solid-liquid two-phase flow problem in a multiphase pump. This method requires high computer performance and memory, long relative computing time, and low efficiency. In order to obtain good results, structured grids are drawn for each structure in the fluid domain. At the same time, the grids of blade surface is encrypted and a 30-layer grid is drawn at the blade tip gap. The specific fluid domain grid is shown in Figure 2.

Therefore, in order to improve the calculation efficiency, before the simulation of the solid-liquid two-phase flow field, the hydraulic head under the condition of oil-water medium is first solved to test the grid independence. A total of 6 sets of fluid domain models with different number of grids were drawn. After grid independence verification, the number of grids and hydraulic head of the computational domain are shown in

Table 2. Grid independence verification.

Inlet Section	Impeller	Diffuser	Outlet Section	Sum Grid Numbers	Hydraulic Head
0.11 million	1.88 million	1.23 million	0.21 million	3.43 million	35.59
0.19 million	2.17 million	1.38 million	0.30 million	4.04 million	35.62
0.30 million	2.36 million	1.49 million	0.41 million	4.56 million	35.65
0.39 million	2.64 million	1.55 million	0.49 million	5.04 million	35.72
0.42 million	2.74 million	1.67 million	0.60 million	5.43 million	35.53
0.51 million	2.93 million	1.72 million	0.71 million	5.87 million	35.61

. The numerical simulation results show that the number of grids has little effect on the hydraulic head of the multiphase pump. At the same time, in order to ensure the suitability of y+, the fourth group is finally determined as the computational grid of this paper.

3.2. Boundary Condition Setting

The liquid medium in the calculation domain is mainly liquid water with a density of 998.2 kg/m³ and crude oil with a density of 882.6 kg/m³. The proportion of oil is 10%, which is evenly distributed in the water in the shape of droplets with a diameter of 0.01 mm. The turbulence model is computed using the SST k-ω model. The wall of the flow channel is set to a no-slip wall, the inlet boundary condition is set to the velocity inlet, the outlet boundary condition is set to the pressure outlet, and the Impeller rotation speed is 2980 r/min. The intersection of the dynamic and static calculation domains is used as the frozen rotor model.

In this paper, Fluent is used for simulation calculation and the pressure-based solver is selected. The velocity field is first obtained from the momentum equation, and then modified from the pressure equation to make the velocity field satisfy the continuity condition. Taking into account the interaction between the two phases, the continuous phase and the discrete phase are alternately coupled for solution.

The particle incident source selects the surface incident source. The incident surface is the annular boundary of the inlet section, and the particles are uniformly perpendicular to the flow field from the normal direction of each grid unit on the incident surface. It is assumed that the particle velocity has been fully accelerated when entering the calculation domain. The incident velocity of the particle is set to be consistent with the inlet velocity. The inlet and outlet surfaces of the multiphase pump are designated as ‘Escape’, and the wall of the flow com-ponent is designated as ‘Reflect’.

The material of the particle is quartz sand with a density of 2650 kg/m³. Figure 3 shows part of the gravel remaining on the impeller surface when the pump is stopped after the actual operation of the multiphase pump. The thickness of the pump impeller blade is 10~12 mm. It can be seen that the diameter distribution of the gravel shown in the figure is wide, which is consistent with the gravel data file provided by the pump station. According to the sand content of the crude oil medium, the particle diameter (d_p) range is selected from 0.1 mm to 5 mm (0.1 mm, 0.2 mm, 0.3 mm, 0.4 mm, 0.5 mm, 1.0 mm, 1.5 mm, 2 mm, 2.5 mm, 3 mm, 3.5 mm, 4 mm, 4.5 mm, 5 mm), the particle concentration is 2%, and The shape of the particles is assumed to be spherical.

Among the spatial discretization schemes selected in this paper, the gradient discretization scheme uses Green-Gauss Cell-Based, the pressure phase discretization scheme is PRESTO!, and the density term discretization scheme uses First Order Upwind, the momentum equation, turbulent kinetic energy, and turbulent dissipation rate schemes are First Order Upwind. The boundary condition design and calculation conditions used in this paper are also summarized in

Table 3. Calculate Settings and operating conditions.

Item	Number/Content
Solver type	Presure-based
Inlet condition	2.54 m/s
Outlet condition	0.3 Mpa
Location of injection	Inlet
Particle velocity	2.54 m/s
particle diameter	0.1 mm, 0.2 mm, 0.3 mm, 0.4 mm, 0.5 mm, 1.0 mm, 1.5 mm, 2 mm, 2.5 mm, 3 mm, 3.5 mm, 4 mm, 4.5 mm, 5 mm
Inlet and Outlet	Escape
All wall	Reflect
Presure-Velocity Coupling Scheme	SIMPLE
Gradient	Green-Gauss Cell-Based
Pressure	PRESTO!
Density	First Order Upwind
Momentum	First Order Upwind
Turbulent kinetic energy	First Order Upwind
Turbulent dissipation rate	First Order Upwind

. Taking the stable flow field of oil and water as the initial flow field, the solid-liquid two-phase calculation is started.

In addition, the particle-wall collision model proposed by Pagalthivarthi et al. [37] is also selected. The normal and tangential restitution coefficients of particle after impacting the surface are polynomial functions of the impact angle. The specific formulas are as follows:

$$e_n=0.993-0.0307\gamma +0.000475{\gamma }^2-0.00000261{\gamma }^3$$

(1)

$$e_t=0.998-0.029\gamma +0.000643{\gamma }^2-0.00000356{\gamma }^3$$

(2)

In the formula: e_n is the normal restitution coefficient; e_t is the tangential recovery coefficient; γ is the wall impact angle.

4. Mathematical Model

4.1. Mixture Model

In this paper, the mixed model of oil and water two medium flow is used for calculation. When solving the oil-water mixture in the mixture model, the velocity slip and volume force are considered, and the heterogeneous model is used to solve the problem.

The expression of the oil-water two-phase continuity equation is as follows:

$${\displaystyle \frac{\partial {\rho }_m}{\partial t}}+\nabla \cdot \left({\rho }_m{u}_m\right)=0$$

(3)

In the formula: ρ_m is the mixing density; u_m is the average mass velocity.

$${\rho }_m={\alpha }_{oil}{\rho }_{oil}+{\alpha }_{water}{\rho }_{water}$$

(4)

$$u_m={\displaystyle \frac{{\alpha }_{oil}{\rho }_{oil}{u}_{oil}+{\alpha }_{water}{\rho }_{water}{u}_{water}}{{\rho }_m}}$$

(5)

In the formula: α_oil and α_water are the volume fractions of crude oil and water, respectively; ρ_oil and ρ_water are the density of crude oil and water, respectively; u_oil and u_water are the average mass velocity of crude oil and water, respectively.

The momentum equation is obtained by summing the dynamic equations of all phases. The oil-water two-phase momentum equation is as follows:

$$\begin{array}{l}{\displaystyle \frac{\partial {\rho }_m{u}_m}{\partial t}}+\nabla \cdot \left({\rho }_m{u}_m{u}_m\right)=-\nabla p+\nabla \cdot \left[{\mu }_m\left(\nabla u_m+\nabla u_m^T\right)\right]+{\rho }_{m}g\\ +F-\nabla \cdot \left({\alpha }_{oil}{\rho }_{oil}{u}_{dr,oil}{u}_{dr,oil}+{\alpha }_{water}{\rho }_{water}{u}_{dr,water}{u}_{dr,water}\right)\end{array}$$

(6)

In the formula: g is the acceleration of gravity; F is the volume force; ${\mu }_m$ is hybrid viscosity; ${u_m}^T$ is the transpose of the velocity matrix.

$${\mu }_m={\alpha }_{oil}{\mu }_{oil}+{\alpha }_{water}{\mu }_{water}$$

(7)

u_dr,oil and u_dr,water is the respective drift velocity of oil-water two phases:

$$u_{dr,oil}=u_{oil}-u_m$$

(8)

$$u_{dr,water}=u_{water}-u_m$$

(9)

4.2. Turbulent Model

The SST k-ω model is one of the k-ω models that can be used at present, including Wilcox k-ω model and Baseline k-ω model. The SST k-ω model uses the Wilcox k-ω model near the wall and the standard k-ε model outside the near-wall region. At the same time, the advantages of the two models in the calculation of the free shear layer are concentrated. In addition, the eddy viscosity coefficient is corrected, and the transmission effect of turbulent shear stress is considered, and then the starting point of flow separation and the size of the separation zone from the smooth surface are accurately calculated.

SST k-ω model mainly solves the convective transport equation of turbulent kinetic energy k and its specific dissipation rate ω. The expression of SST k-ω model is as follows:

$${\displaystyle \frac{\partial }{\partial t}}(\rho k)+{\displaystyle \frac{\partial \left(\rho v_{i}k\right)}{\partial x_i}}={\displaystyle \frac{\partial }{\partial x_j}}\left[\left(\mu +{\displaystyle \frac{{\mu }_t}{{\sigma }_{\kappa }}}\right){\displaystyle \frac{\partial k}{\partial x_j}}\right]+G_{\kappa }-\rho {\beta }^{*}k\omega$$

(10)

$${\displaystyle \frac{\partial }{\partial t}}\left(\rho \omega \right)+{\displaystyle \frac{\partial }{\partial x_j}}\left(\rho v_j\omega \right)={\displaystyle \frac{\partial }{\partial x_j}}\left[\left(\mu +{\displaystyle \frac{{\mu }_t}{{\sigma }_{\omega }}}\right){\displaystyle \frac{\partial \omega }{\partial x_j}}\right]+G_{\omega }-\rho \beta {\omega }^2+D_{\omega }$$

(11)

In the formula: t is time, ρ is fluid density, k is turbulent kinetic energy; v is the time average speed_, v_i and v_j are the time average speed in i and j directions, respectively. x_i and x_j are the displacements in Cartesian coordinates in i and j directions, respectively. β is the coefficient of thermal expansion, β* is the coefficient of thermal expansion due to the compressibility of the fluid, ω is the specific dissipation rate; μ is the dynamic viscosity of the fluid, μ_t is the turbulent viscosity; σ_k is the Prandtl number corresponding to the turbulent kinetic energy, and G_k is the production term of the turbulent kinetic energy caused by the average velocity gradient. σ_k is the turbulent Prandtl number with a specific dissipation rate, G_ω is the generation term of a specific dissipation rate, and D_ω is the orthogonal divergence term.

A new dissipation source term is added to the turbulent dissipation rate ω equation in the SST k-ω model to solve the phenomenon of excessive turbulence in the local area.

$$D_{\omega }={\displaystyle \frac{2\rho (1-F_1)}{\omega {\sigma }_{\omega 2}}}\left({\displaystyle \frac{\partial k\partial \omega }{\partial x\partial x}}+{\displaystyle \frac{\partial k\partial \omega }{\partial y\partial y}}+{\displaystyle \frac{\partial k\partial \omega }{\partial z\partial z}}\right)$$

(12)

In the formula: F₁ is the distance from the wall.

4.3. DPM Model

According to the force of flow field on particles, the DPM model solves the trajectory of particles through integration. The integral formula is as follows:

$$m_p{\displaystyle \frac{d{u}_p}{dt}}=F_D+F_B+F_{VM}+F_P+F_R+F_M+F_S$$

(13)

In the formula: the subscript p represents the parameters of the particles; m_p is particle mass; u_p is the particle velocity; F_D is drag force; F_B is buoyancy caused by gravity; F_VM is virtual mass force; F_p is the pressure gradient force; F_R is the Coriolis force and centrifugal force in the rotating system; F_M is Magnus lift; F_S is the Saffman lift.

Due to the inertia of the particle, the V_slip between the particle and the fluid is calculated as follows:

$$v_{slip}=v_f-v_p$$

(14)

In the formula, v_f is the fluid velocity and v_p is the particle velocity.

The drag force on the particle can be calculated by v_slip:

$$F_D=m_p{\displaystyle \frac{v_{slip}}{{\tau }_r}}$$

(15)

In the formula: the subscript f represents the parameters of the fluid; u is the speed; τ_r is the particle relaxation time.

$${\tau }_r={\displaystyle \frac{{\rho }_p{d}_p^2}{18{\mu }_f}}{\displaystyle \frac{24}{C_D{\mathrm{Re}}_p}}$$

(16)

$${\mathrm{Re}}_P={\displaystyle \frac{{\rho }_f{d}_p\left|v_{slip}\right|}{{\mu }_f}}$$

(17)

In the formula: ρ is the density; d_p is the particle diameter; μ_f is fluid dynamic viscosity; Re_p is the particle Reynolds number; C_D is the drag coefficient.

In this paper, the particles are assumed to be spherical, and the calculation formula of the drag coefficient proposed by Morsi and Alexander [38] is adopted. The drag coefficient is:

$$C_D=a_1+{\displaystyle \frac{a_2}{{\mathrm{Re}}_p}}+{\displaystyle \frac{a_3}{{\mathrm{Re}}_p^2}}$$

(18)

In the formula: a₁, a₂, and a₃ are empirical constants, using Morsi’s recommended values.

When particles are in the fluid, its own gravity creates a certain buoyancy:

$$F_B=m_p\left(1-{\displaystyle \frac{{\rho }_f}{{\rho }_p}}\right)g$$

(19)

In the formula: g is the acceleration of gravity.

The virtual mass force is generated when the velocity between the particle and the fluid is inconsistent:

$$F_{VM}=C_{VM}{m}_f\left({\displaystyle \frac{d{u}_f}{dt}}-{\displaystyle \frac{d{u}_p}{dt}}\right)$$

(20)

In the formula: C_VM is the virtual mass coefficient, taking 0.5.

There are pressure gradients in the flow field, resulting in the pressure gradient force on the particles:

$$F_P=m_p{\displaystyle \frac{{\rho }_f}{{\rho }_p}}{u}_p\nabla u_f$$

(21)

In rotating machinery, there must be centrifugal force on the particles. The calculation formula is:

$$F_R=m_p\left(-2\omega \times u_p-\omega \times \omega \times r_p\right)$$

(22)

In the formula: ω is the relative angular velocity between particles and fluid; r_p is the set vector from the particle to the coordinate origin.

When a velocity gradient exists in the flow field, particles are rotated around their center of mass and Magnus lift is generated. In cases where the velocity gradient is not aligned with the direction of particle motion, Saffman lift is experienced by particles. The calculation formulas of the two forces are as follows:

$$F_M={\displaystyle \frac{1}{2}}{A}_p{C}_M{\rho }_f{\displaystyle \frac{\left|v_{slip}\right|}{\left|\omega \right|}}\left[v_{slip}\times \omega \right]$$

(23)

$$F_S=K_{S}4r_p^2{\rho }_f{\left|{\nu }_f{\displaystyle \frac{\partial V_{fi}}{\partial x_i}}\right|}^{{\displaystyle \frac{1}{2}}}\left(V_{fi}-V_{pi}\right)\mathrm{sgn}\left({\displaystyle \frac{\partial V_{fi}}{\partial x_i}}\right)$$

(24)

In the formula: A_p is the projection area of particles; C_M is the rotational lift coefficient; K_S is the empirical coefficient; sgn is the sign function.

4.4. Wear Model

In this paper, the Oka [39] wear model is used to predict the wear of the surface of the multiphase pump. The wear rate in the Oka model is expressed as follows:

$$E=E_{90}{\left({\displaystyle \frac{V}{V_{ref}}}\right)}^{k_2}{\left({\displaystyle \frac{d}{d_{ref}}}\right)}^{k_3}f\left(\gamma \right)$$

(25)

In the formula: E₉₀ is the impact wear rate of particle at 90°; V is the impact velocity of particle; V_ref is the reference speed; d and d_ref are particle diameter and particle reference diameter; k₂ and k₃ are model velocity constants and model diameter constants. f(γ) is the impact angle function.

The calculation formula of the impact angle function is as follows:

$$f\left(\gamma \right)={\left(\sin \gamma \right)}^{n_1}{\left(1+H_V\left(1-\sin \gamma \right)\right)}^{n_2}$$

(26)

In the formula: H_V is the Vickers hardness of the wall material; n₁ and n₂ are model constants.

Due to the different values of each coefficient in the Oka model, it is suitable for different wear prediction situations. For the wear of steel by quartz sand, the values of each coefficient are shown in

Table 4. Wear formula parameter table.

C	Sand-Steel	C	Sand-Steel
$E_{90}$	6.154e−4	$k_2$	2.35
$H_V$	1.8 (GPa)	$k_3$	0.19
$n_1$	0.8	$V_{ref}$	326 ($\mathsf{\mu }\mathrm{m})$
$n_2$	1.3	$d_{ref}$	104 (m/s)

.

5. Experimental Study on External Characteristics of Multiphase Flow Pump

Before the simulation calculation, experimental tests were conducted on the mixer pump. The test device is a 4-stage multiphase pump, and the test medium is pure water. In the test, the pressure gage is used to measure the inlet and outlet pressure of the pump, the flow meter is used to measure the pump flow, the Tachometer is used to measure the rotational speed, and the Torquemeter is used to measure the shaft power. The specific instrument accuracy is shown in

Table 5. Test instrument precision table.

Instrument	Range	Precision	Unit
Inlet pressure gage	−0.1~0.25	±0.02%	MPa
Outlet pressure gage	0~4	±0.06%	MPa
Flow meter	0~380	±0.2%	m3/h
Tachometer	20~33,000	±0.03%	rpm
Torquemeter	0~800	±0.2%	N·m

. The external characteristics test system of pure water is shown in Figure 4, which includes the field layout of the test system and the schematic diagram of the test system. The external characteristic results of the test are compared with the corresponding simulation results, and the results are shown in Figure 5. It can be seen from Figure 5 that the simulation results of the external characteristics of the multiphase pump under the condition of 0.7 times the rated flow rate to 1.1 times the rated flow rate are similar to the test results, and the overall error is small near the rated flow rate. The large gap between the simulation results and the experimental data on the middle point of Figure 4 May originate from the instability of the calculation. Therefore, the numerical research method selected in this paper has high reliability and can be used for further research.

After pure water testing in the test bed, the pump was applied to the field site. When the pump is operating in the field, the delivery medium is a mixture of oil, water and particles. The flow rate and head were recorded in the test operation condition, and the accuracy of the calculation was further verified by combining the simulation calculation method used in this paper. The comparison between the test run results and the simulation calculation is shown in Figure 6. It can be seen that when the three media of oil, water and particles are transported, the gap between the test results and the simulation results of pump head, efficiency and power is small with the rise of speed and the increase of flow rate. At low speed and low flow, the test efficiency value of the pump and the simulation are larger, which is related to the pump working in the bias condition. As the speed and flow gradually approach the design value, the error value gradually decreases.

In order to validate the OKA model, the calculated results of the trials conducted by Romo et al. [40] were used, and the results are shown in Figure 7. The normalized erosion rate R_norm is used to characterize the wear degree of the surface, and the calculation formula is given in Equation (27). It is evident that with an increasing impact angle, the normalized erosion rate exhibits a rising trend until it reaches its peak, followed by a decline and stabilization. Through quantitative analysis of both experimental and simulated results, the predicted prediction error falls within the range of 0.56% to 23.33%. Although the error is larger when the impact Angle is 45 degrees, the influence on the result is relatively small because the magnitude of this value is 10⁻¹⁰.

R_norm = Q/(Vs × A)

(27)

where, Q is the volume flow rate (m³/s), Vs is the particle impact velocity, and A is the apparent erosion area of the sample.

6. Results Analysis

6.1. Solid-Liquid Two-Phase Movement in Multiphase Pump

Taking the 0.4 mm particle calculation condition as an example, the trajectory of particle and fluid and the distribution of particle in the solid-liquid two-phase condition are analyzed. Figure 8 shows the solid-liquid two-phase trajectory inside the multiphase pump.

It can be seen from the left figure that the liquid phase flow inside the impeller domain is more in line with the shape of the blade. The liquid phase velocity is the maximum at the inlet of the impeller domain, and the liquid phase velocity decreases gradually along the flow direction. The flow inside the diffuser domain is more disordered, and a vortex is generated in the middle of the flow channel of the diffuser domain. The flow channel vortex squeezes the flow channel area of the liquid phase, resulting in a slight increase in the liquid phase velocity at the outlet of the diffuser. In the right diagram, before entering the impeller domain of the multiphase pump, the particle flow rate is low and moves linearly along the inlet section. When the particle enters the impeller domain, its motion direction deflects and accelerates rapidly, as shown in the red box in the figure. At the same time, the particle will impact the head of the impeller blade, and the particle will separate along both sides of the blade head and continue to move with the fluid. In the black box in the figure, there is a part of low-speed particle inside the impeller domain. These low-speed particles are affected by centrifugal force on the blade surface and move radially along the blade. When they move to the tip clearance, they will follow the fluid and leak to the adjacent flow channel. After the particle enters the interior of the diffuser domain, the sand grains are more obviously affected by the liquid phase. In the middle of the diffuser region, due to the influence of the liquid phase vortex, the particle accumulates in groups. Comparing the liquid streamline and the particle trajectory, it is found that in most areas, the similarity of the two trajectories is high, which further indicates that the particle movement is closely related to the kinetic energy change of the fluid.

Figure 9 shows the phenomenon of particle accumulation in the flow channel of the multiphase pump and the distribution of particle mass concentration in different blade heights. The blade height = 0.1, 0.5 and 0.9 represent the hub, the middle of the blade height and the vicinity of the rim, respectively. It can be seen from Figure 9 that the particles in the impeller domain are mainly distributed, and the particle concentration at the rim is significantly smaller. It can be seen from the solid line black frame in Figure 9a,b that some of the particle in the impeller flow channel is distributed in strips. After the particle enters the diffuser domain, the concentration of particle appears at the hub, middle of the blade height and near the rim of the diffuser. The particle moves from the hub of the diffuser inlet along the red line in Figure 9a, and finally moves to the diffuser outlet near the rim. In Figure 9a and the sand concentration distribution cloud map at 0.5 times the leaf height, it can be seen that the sand formed a agglomerated area in the middle of the flow channel due to the influence of the vortex in the diffuser flow channel. Simultaneously, the spatial distribution of particles within the diffusion channel exhibits significant heterogeneity.

6.2. Effect of Particle Diameter on Wear of Multiphase Pump

In order to explore the wear distribution area inside the multiphase pump and the influence of particle diameter on the wear distribution, Figure 10, Figure 11, Figure 12 and Figure 13 show the wear distribution of the blade surface under different particle diameters. Figure 10 and Figure 11 show the effect of particle diameter on the wear of the pressure surface and the suction surface of the impeller, respectively. It can be seen from Figure 10 and Figure 11 that the particle diameter has a great influence on the wear distribution of the impeller. The increase of particle diameter makes the wear position of the pressure surface develop from the front to the back, and the wear area almost occupies the whole pressure surface of the blade when the particle diameter is large. At the same time, the particle near the pressure surface causes wear along the radial distribution of the blade due to the movement of the rotating machinery and the influence of centrifugal force. The suction surface is mainly the first half of the wear is more serious, but there is a region between the blade head and the suction surface wear position without wear, as shown in the black box in Figure 11. In the red box of Figure 11, it can be seen that the wear of the impeller blade head is very serious. Even if the wear on the blade surface is small, the blade head still belongs to the part where the wear occurs. Combined with the dotted red box in Figure 8, it can be seen that the wear of the blade head comes from the vertical impact of the particle relative to the head surface.

Figure 12 and Figure 13 are the wear distribution nephograms of the diffuser surface under different particle diameters of particle. It can be seen from Figure 12 and Figure 13 that the increase of diameter will lead to the aggravation of wear. The wear area of the pressure surface of the diffuser is mainly distributed at the edge of the diffuser, and the closer to the diffuser outlet, the more serious the wear is. With the increase of particle diameter, the wear of the pressure surface of the diffuser is first evenly distributed along the circumference, and when the particle diameter is large, it is concentrated on a few diffusers. On the suction surface of the diffuser, under different particle diameters, the wear occurs on the blade surface near the inlet of the diffuser. The wear degree of the suction surface of the diffuser increases with the increase of the particle diameter. It can be seen from the solid black box in Figure 8 that the inlet of the suction surface is impacted by the particle, and the surface wear of the pressure surface is formed by the impact of the particle entrained by the fluid in the red box of the solid line in Figure 8. Because the particle is affected by the vortex in the flow channel of the diffuser, the particle is not easy to move to the suction surface of the blade, so there is no wear phenomenon in the back section of the suction surface of the diffuser when the particle diameter is large.

Figure 14 shows the average wear rate (wear rate after averaging the blade area) and maximum wear rate of impeller and diffuser surface with the increase of particle diameter. It can be seen from Figure 14a that the average wear rate of the impeller surface is linearly distributed before the particle diameter of 3.5 mm, and the average wear rate after 3.5 mm changes little due to the influence of the gravity of the particle itself.

The average wear rate of the diffuser surface is smaller than that of the impeller, and the average wear rate is linearly related to the particle diameter in the range, as shown in Figure 14a. In Figure 14b, the maximum wear rate of the diffuser surface is much larger than that of the impeller surface. Only when the particle diameter is very small, the maximum wear rate of the impeller and the diffuser is close. The maximum wear rate of the impeller reaches the peak when the particle diameter is 3.5 mm, and the maximum wear rate of the diffuser surface is weakened by the particle diameter when the diameter is greater than 2.5 mm.

6.3. Wear Rate Prediction Mathematical Model of Multiphase Pump

It can be seen from the previous analysis that the particle diameter is an important factor affecting the wear in the pump. In order to predict the wear of equipment through the sediment concentration of crude oil extraction in actual production engineering, a polynomial regression mathematical model is established to predict the change rule of wear rate for particle diameter and wear rate. By introducing a polynomial regression model, we can set

$$y=B_n{x}^1+B_{n-1}{x}^2+B_{n-2}{x}^3+...+B_3{x}^{n-2}+B_2{x}^{n-1}+B_1{x}^n+B_0$$

(28)

Here, x, y represent different variables respectively. B₀ is the truncation constant of the model, and B₁ to B_n are the undetermined coefficients. Let (x_n, Y_n) be a set of samples, and (x_n, y_n) be the observed value of the sample. It can be seen that the n-order polynomial has 1 to n unknown fitting coefficients, and the purpose of fitting is to find the best n + 1 fitting coefficients.

Bring n sets of data (x_n, y_n) into Formula (27), then we have

$$y_i=B_n{x_i}^1+B_{n-1}{x_i}^2+B_{n-2}{x_i}^3+...+B_3{x_i}^{n-2}+B_2{x_i}^{n-1}+B_1{x_i}^n+B_{0,i}$$

(29)

where: i = 1, 2, 3..., n.

According to the above analysis, the set of fitting data is determined, as shown in

Table 6. Effect of particle diameter on wear rate of multiphase pump.

Diameter (mm)	Impeller	Diffuser	Impeller	Diffuser
	Ave Wear (kg × m−2 × s−1)		Max Wear (kg × m−2 × s−1)
0.1	3.00999 × 10−7	3.16854 × 10−8	0.00134	6.53395 × 10−4
0.4	1.60537 × 10−5	9.86224 × 10−7	0.03253	0.01113
1	3.68153 × 10−5	1.29888 × 10−5	0.04283	0.3055
1.5	5.54406 × 10−5	1.69755 × 10−5	0.02518	0.56944
2	7.96815 × 10−5	2.25272 × 10−5	0.03583	0.55437
2.5	1.23332 × 10−4	3.02421 × 10−5	0.07758	0.8168
3	1.39623 × 10−4	3.43214 × 10−5	0.08231	0.77447
3.5	1.5 × 10−4	3.538 × 10−5	0.23696	0.8765
4	1.79798 × 10−4	4.10062 × 10−5	0.1489	0.8959
4.5	1.6575 × 10−4	7.49759 × 10−5	0.13014	0.88
5	1.57479 × 10−4	7.97509 × 10−5	0.04195	0.75

below:

The data in the above table are processed to calculate the polynomial regression equation between particle diameter and wear rate. In order to avoid under-fitting and over-fitting, the R² of each order fitting equation is used as a reference to determine the best fitting order. The change of R² value of each polynomial regression equation with the order is shown in Figure 15. Considering the accuracy and simplicity of the fitting equation, the R² value should be as large as possible. At the same time, considering the actual situation, the fitted regression equation should satisfy that the average wear rate and the maximum wear rate are not zero under any particle diameter, and the wear should not occur when the particle diameter is 0.

Under the previous constraints, the polynomial regression model is determined and a series of statistical tests are performed to diagnose its predictive ability for sample data. The requirement of a good regression equation is that under the condition that the R² value is as large as possible, the F of the regression equation is less than F_0.05, and the F_0.05 value can be determined by the degree of freedom and the number of independent variables of the regression equation. Under the above requirements, the statistical test of the regression equation when the particle diameter changes is shown in

Table 7. Statistical test of regression equation when particle diameter changes.

Model	Degrees	Variables	R2	F		F0.05	SSE
${\mathrm{Ave}\mathrm{wear}}_{\mathrm{Impeller}}$	7	3	0.99648	1040.47	>	4.347	3.9072 × 10−10
${\mathrm{Ave}\mathrm{wear}}_{\mathrm{Diffuser}}$	4	6	0.98527	123.631	>	6.163	1.2033 × 10−10
${\mathrm{Max}\mathrm{wear}}_{\mathrm{Impeller}}$	6	4	0.87426	20.1197	>	4.534	0.00917
${\mathrm{Max}\mathrm{wear}}_{\mathrm{Diffuser}}$	6	4	0.98884	244.591	>	4.534	0.03481

.

It can be seen from Table 7 that the goodness of fit R² of the three regression equations is greater than 0.95, indicating that the regression equation can explain more than 95% of the variable changes. Only the prediction model of the maximum wear rate of the impeller has a goodness of fit R² greater than 0.87, and the overall goodness of fit is good. Moreover, the F of the four regression equations is less than F_0.05, indicating that the univariate polynomial regression equation obtained by polynomial regression analysis has a high degree of fitting, which can reflect the effect relationship between particle diameter and average wear rate. It is proved that the functional relationship between particle diameter and average wear rate is significant. At the same time, the residual sum of squares of the polynomial regression equation is also small and conforms to the normal distribution.

Under the above fitting requirements, the polynomial regression equation between particle diameter and wear rate is determined as follows:

$${\mathrm{Ave}\mathrm{wear}}_{\mathrm{Impeller}}=1.70496\times {10}^{-5}\times x^1+2.01786\times {10}^{-5}\times x^2-3.45993\times {10}^{-6}\times x^3$$

(30)

$$\begin{array}{c}{\mathrm{Ave}\mathrm{wear}}_{\mathrm{Diffuser}}=1.97868\times {10}^{-5}\times x^1-4.61762\times {10}^{-5}\times x^2+5.7232\times {10}^{-5}\times x^3\\ -2.76024\times {10}^{-5}\times x^4+5.6916\times {10}^{-6}\times x^5-4.19368\times {10}^{-7}\times x^6\end{array}$$

(31)

$${\mathrm{Max}\mathrm{wear}}_{\mathrm{Impeller}}=0.09651\times x^1-0.10638\times x^2+0.04625\times x^3-0.00572\times x^4$$

(32)

$${\mathrm{Max}\mathrm{wear}}_{\mathrm{Diffuser}}=0.20915\times x^1+0.13011\times x^2-0.04609\times x^3+0.00357\times x^4$$

(33)

The prediction curves of the average wear rate and the maximum wear rate with the change of particle diameter are drawn, as shown in Figure 16. It can be seen from Figure 16a that with the increase of particle diameter, the average wear rate of multiphase pump impeller increases sharply and then changes slowly, and the average wear rate of diffuser increases slowly with the increase of particle diameter. In Figure 16b, the maximum wear rate of impeller and diffuser surface increases first and then decreases with the increase of particle diameter, but the influence of particle diameter on the maximum wear rate of diffuser surface is much greater than that of impeller surface. The polynomial regression equations of average wear rate and maximum wear rate have good fitting degree at each sample point.

7. Conclusions

In this paper, the external characteristics and solid-liquid two-phase flow of multiphase pump are studied by means of model test and numerical calculation. The trajectory of particle in the booster unit, the distribution of particle in two-dimensional and three-dimensional space, the wear of impeller and diffuser surface, the average wear rate and maximum wear rate of impeller and diffuser in the pump with particle diameter are analyzed. Finally, the mathematical prediction models of average wear rate and maximum wear rate are constructed with particle diameter as variable. Through this study, the following conclusions are obtained:

Before entering the impeller domain of the multiphase pump, the particle velocity is low and moves linearly along the inlet section. When the particle enters the impeller domain, its movement direction deflects and accelerates faster. At the same time, the particle will impact the head of the impeller blade. The particles are separated, and continue to move. The particles are mainly distributed at the wheel rim and the particle concentration is obviously smaller. After entering the diffuser domain, the particles are concentrated and distributed near the flange near the outlet. Due to the influence of the vortex in the diffuser flow channel, the particle forms a agglomerated area in the middle of the flow channel.

The increase of particle diameter makes the wear position of the impeller pressure surface develop from the front to the back, and the wear area almost occupies the whole blade pressure surface when the particle diameter is large, while the suction surface is mainly the first half of the wear is more serious. The wear area of the pressure surface of the diffuser is mainly distributed at the edge of the guide impeller, and the closer to the outlet of the diffuser, the more serious the wear. On the suction surface of the diffuser, under different particle diameters, the wear occurs on the blade surface near the inlet of the diffuser, and the wear degree gradually increases with the increase of particle diameter. Overall, the average wear of the diffuser surface is less than that of the impeller, and the maximum wear rate of the diffuser surface is much larger than that of the impeller surface. Only when the particle diameter is very small, the maximum wear rate of the impeller and the diffuser is close.

The sand particle diameter was selected as the independent variable, and the polynomial regression models of the average wear rate and the maximum wear rate of the impeller and the diffuser were established respectively. The test value F of each equation is calculated, and the critical value F_0.05 is found through the degrees of freedom of the equation. The F of the four regression equations is less than F_0.05, indicating that the univariate polynomial regression equation obtained by polynomial regression analysis has a high degree of fitting, which can reflect the relationship between particle diameter and average wear rate. It is proved that the functional relationship between particle diameter and average wear rate is significant. With the increase of particle diameter, the average wear of the impeller of the multiphase pump first rises sharply and then changes slowly, while the average wear of the diffuser increases slowly with the increase of particle diameter. The maximum wear rate of impeller and diffuser surface increases first and then decreases with the increase of particle diameter, but the influence of particle diameter on the maximum wear rate of diffuser surface is much greater than that of impeller surface.

In summary, through the analysis of particle movement trajectory, particle distribution, blade surface wear and the establishment of wear rate prediction formula, it can explain the wear phenomenon of multiphase pump in engineering practice, and provide reference for the screening of pump particle size and the optimization of pump impeller. At the same time, the influence of multiple factors combined on blade wear can be explored in subsequent studies.