CC BY
16
DOI: 10.17516/1997-1397-2021-14-2-213-223 УДК 517.958:531.32
Numerical and Analytical Modeling of Centrifugal Pump
Viktor A. Ivanov*
Siberian Federal University Krasnoyarsk, Russian Federation
Nikolay V. Erkaev1"
Institute of Computational Modeling SB RAS Krasnoyarsk, Russian Federation
Received 10.08.2020, received in revised form 10.12.2020, accepted 21.01.2021 Abstract. A mathematical model of non-stationary rotation modes of the rotor of a centrifugal pump is constructed. It is based on preliminary calculation in ANSYS Fluent package with subsequent analytical calculation taking into account the specified parameters.
Keywords: mathematical modeling, hydrostatic bearing, centrifugal pump.
Citation: V.A. Ivanov, N.V. Erkaev, Numerical and Analytical Modeling of Centrifugal Pump, J. Sib. Fed. Univ. Math. Phys., 2021, 14(2), 213-223. DOI: 10.17516/1997-1397-2021-14-2-213-223.
Introduction
Low-flow centrifugal pumps are widely used in technology. They are especially often used in various cooling systems. The design and calculation of centrifugal pumps are fairly well studied and described in the literature [1-9]. When calculating such a pump and the hydrostatic bearing included in it (Fig. 1), simplified formulas are often used, containing rather vague empirical coefficients of nozzle flow and oil flow through the bearing ends. In principle, self-consistent determination of these coefficients is possible with a full simulation of hydrostatic bearings in the ANSYS Fluent package. However, in this case, a complete calculation would require a powerful computational resources. In this case, the best option seems to be the limited use of a computing package such as ANSYS Fluent in conjunction with an analytical model, which allows one to simulate the device operation with minimal computational costs. We call further this approach as "hybrid modeling". This work is devoted to hybrid simulation of the nonstationary dynamics of the centrifugal pump rotor during its rotational acceleration from the initial steady position to a final stationary rotational regime.
A schematic diagram of a centrifugal pump is shown in Fig. 1 with the following designations: 1 — impeller; 2 — hydraulic bearing; 3 — a pair of thrust bearings; 4 — inlet pipe; 5 — outlet branch pipe; 6 — auxiliary pump; 7 — the rotor of the electric motor; 8 — electric motor stator; 9 — motor housing; 10 — housing of one of the pumps of the electric pump unit; 11 — dashed lines indicate the direction of movement of the working fluid; 12 — throttle space; 13 — throttle; 14 — pocket; 15 — a system of holes in the shaft, providing a return of the working fluid back to the entrance to the impeller.
[email protected] https://orcid.org/0000-0001-8993-6400 © Siberian Federal University. All rights reserved
Fig. 1. Schematic diagram of a centrifugal pump
1. Mathematical model
Forces and moments acting on the rotor of a low-flow centrifugal pump with two symmetrical hydrostatic bearings are shown in Fig. 2, where mg is the weight of the rotor (in the figure this
force is divided by two, since it is distributed over two bearings), N is the normal reaction of the support, Rc is the static reaction of the hydraulic bearing, Pr is the radial force, taking into account the fluid pressure on the impeller blades (its direction depends on the design features of the pump), Mrk is impeller moment, Mvt is viscous friction moment, Mst is dry friction moment, Md is engine torque, Ri and R2 are the radii of the shaft and bearing housing, d is
the displacement of the center of mass of the rotor relative to the axis, W is the reaction of the liquid layer to the acceleration of the rotor axis. The rotor dynamics equations are as follows:
i'dty = 2Rcy(y,u) - 2PrM - mg + 2N(y) - Wy(y, dy)
= Md(u) - 2Mrk(u) - 2Mst(y) - 2Mvt(y, u),
k dt
where I is the moment of inertia of the rotor about its axis, u is the angular speed of the rotor, t is the time since the pump was started, A = R2 — R\ is the average clearance. The values Rc, Pr, Mrk, Mst, Mvt are indicated with a factor of 2, since the pump has two bearings.
Oscillations of the rotor along the y axis are caused by two force reactions of the liquid layer, designated Rc and W. The expression of the force Rc, depending on the displacement of the rotor, is given in [10], in which the load was assumed to be applied strictly along the vertical axis y. The expression for the force W caused by the hydrodynamic reaction of the liquid layer to the acceleration of the rotor is obtained in Appendix 1. In case of strictly vertical loading, the component of the layer reaction can be written in the following form:
RCy — l • d •
Pn • (Vc • ndl/4)2 Pn • {^c • nd?c/4)2
(^c • ndC/4)2 + (2^f (y))2 (Vc • ndC/4)2 + (2^f (-y))2
(2)
where l is the bearing width, d is the rotor diameter, ¡ic is the nozzle flow rate, ^ is the flow rate through the bearing ends, dc nozzle diameter, Pn is the pump outlet pressure, f (x) is the function given in [10] and characterizing the area of the gap sector
f (y) — l • A •
t/4
sin2(t)+(A +cos(4>))
1/2
dfi.
The pressure at the pump outlet is found from the following equation [11]:
Pn = pu2Rj n,
where n is the pump efficiency, p is the density of the working fluid.
When using formula (2), it is especially difficult to determine the coefficients /i,c and ^. In typical calculations, these coefficients are selected on the basis of tables obtained empirically for a certain range of diameters, as a result of which the calculation error increases. In addition, this formula assumes the average pressure, which creates a force directed strictly along the coordinate axes, which is an additional source of error. In reality, one must take into account the distribution of the pressure along the surface, at each point of which the elementary force vector is directed along the normal. Therefore, when finding the resulting force, it is necessary to integrate this vector along the surface. To eliminate these errors, a series of calculations was made for a hydrostatic bearing with different rotor eccentricities in the ANSYS Fluent package in order to determine the flow rates through the nozzle and the bearing ends, as well as the pressure correction function. Taking into account these coefficients, the modified formula (2) takes the following form:
Rcy — l • d •
(Vc(y) • ndC/A)2 x(y) __
(Vc(y) • nd2/4) + (2vt(y)f (y))
Pn ^ {^(-y) • nd2/4)2 X( J . (3)
(Vc(-y) • nd2c/4)2 + (2vt(-y)f (-y))2
2
0
Pn
Here x is a correction factor that takes into account the pressure distribution. Next, we pay a particular attention on a way to obtain these coefficients, depending on the rotor displacement, using the ANSYS Fluent package. To do this, in the software package, we simulate the bearing lubricant layer for a specific eccentricity. At the same time, inlets of the lubricating layer through the nozzles are created in the model. A model of a lubricating layer of finite thickness is shown in Fig. 3. Here, the inlets of the lubricating layer are labeled A, B, C, D.
Fig. 3. Finite element model of the lubricating layer
To obtain a uniformly ordered mesh, it is important to split the lubricant layer model into five elements, i.e. a layer around the shaft and four lubricant inlets through the nozzles. We divide the layer around the shaft into 100 elements along the corner, 10 elements along the layer thickness and 10 elements along the length. In addition, we divide the grease inlet through the nozzle into 20 elements in angle, 5 elements in thickness and 5 elements in length. After building the model, we set the load and boundary conditions. Here, the pump pressure Pn is set at the nozzle inlet ends, atmospheric pressure is set at one bearing end, and at the other end zero liquid flow rate is assumed. Since the shaft rotates in a stationary housing, we fix the outer boundary between the housing and the lubricating layer. In this case, the rotation speed is set at the inner interface between the shaft and the lubricating layer. The resulting pressure distribution is shown in Fig. 4. According to Fig. 4, the liquid entering through the nozzles is distributed along the rotation of the shaft. In this case, the maximum pressure peak is observed at the liquid outlet from the lower nozzle. We also observe increased pressure along the end with the cuff, since there is nowhere for the liquid to exit.
Based on the pressure difference from below (in the vicinity of the point of minimum thickness of the lubricating layer) and from above (in the vicinity of the point of maximum layer thickness), we determine the lifting force Rc. Knowing the oil flow rate through the ends, we determine the functions ^c(y) and ^t(y) based on the balance of the fluid flow in the pump [12]:
* • dl 2(Pn - P) ^ 2 f 2P Vc---:--\ - = Q = • J\ -,
4 V P VP
where Q is the flow rate through the nozzle/ends. The found functions are shown in Fig. 5.
The found functions are used in the formula (3) to determine the lift of the rotor. Fig. 6 shows plots of the Rc functions corresponding to the analytical (2), hybrid (3) and numerical (based on the ANSYS package) models for comparison.
Fig. 4. Pressure distribution in the lubricating layer
-0.8 -0.6 -0.4 -0.2 0
y!A
0.2
-----------1
----------2
— 3
0.4 0.6 0.8 1
1
0.75 0.5 * 0.25 0
Fig. 5. Functions and ^t, X, depending on the shaft displacement. Curves 1 and 2 correspond to the flow rate coefficients ^c and ^ through the nozzles and bearing ends, respectively; curve 3 shows correction factor for pressure x
Fig. 6. Lifting force in the bearing. Curves 1, 2 correspond to the formulas (2) and (3), curve 3 is obtained from the ANSYS Fluent calculations.
Fig. 6 shows that at low eccentricities all curves are fairly close to each other. However, at eccentricities greater than 0.3, the curves obtained on the basis of the formula (3) and direct ANSYS calculation go up significantly parallel to each other. These curves have a particularly steep increase for the eccentricities larger than 0.6, which is associated with the transition from hydrostatic to hydrodynamic operating regime of the bearing. In this case, the curve corresponding to the formula (2) has a weak growth and remains almost parallel to the abscissa axis.
Next, we consider the additional force W caused by a dynamic reaction of the liquid layer
on the rotor acceleration. This force is linearly dependent on the acceleration vector:
( d?x d2y\
W = Ka' a = (, dè) ■ (4)
Components of the matrix coefficient K are derived in Appendix 1.
The expression for the radial force Pr acting on the impeller is given in [13]:
Pr =0.1 • p ■ D3 ■ b2 ■ n ■ u2,
where D2 is the diameter of the impeller outlet, b2 is the width of the impeller at the outlet. The normal support reaction is:
¡Pr + mg - R y = -A
N = ^ 2 c y
\0, y> -A.
The dry friction moment is determined by the following expression:
( (2Pr + mg - Rc)dftr .
Mst = \ -4-' y = -A
0, y > -A,
where ftr is the dry friction coefficient.
The moment of viscous friction is determined from the following expression [10]:
n
O 1 f dS
Mvt = 0.25ld>w— 1
(^^ + cos(^))2
The torque generated by the engine is given by formula [14]:
Md = J - J^,
where J and Ji are dynamic and static coefficients of the torque-mechanical characteristics of the electric motor.
The moment on the impeller is determined by the expression:
unuRi „ n
Mrk = —— + pQRiu, s
where s is the axial clearance between the impeller and the casing, u is the dynamic viscosity of the oil, Qp is the flow (fluid leakage from the impeller), determined by formula [15]:
Qp = UpnDisuR2^2n,
where uP is the flow coefficient in the front axial clearance between the impeller and the pump casing [15], Di is the impeller diameter at the inlet.
2. Calculation results
For brevity, we use the following notation:
Fy = 2Rcy - 2Pr - mg + 2N.
Substituting all the determined force expressions into the equation (1) and taking into account the equation (3) and zero initial conditions, we write the system of equations in the Cauchy normal form:
dVy_ dt
Fy (m - Kn) + FxK21 (m - K11) (m - K22) - K12K21'
dy = v dt Vy'
dw = Md - 2Mrk - 2Mst - 2Mvt
dt I '
Vy (0) = 0,
(5)
y(0) = -A,
w(0) = 0.
Here y(0) = -A, since the origin of coordinates is chosen in the center of the body, and the rotor at the initial moment is in the lowest position being in contact with the body. For comparison, we carry out two variants of calculation with different expressions of the lifting force Rc, determined by the formulas (2) and (3). Solving the system of differential equations (5) by the fourth order Runge-Kutta method, we obtain the angular velocity, the displacement and speeds of the rotation axis along the coordinate axes as functions of time (Figs. 7-9). The calculations were performed for the following input parameter: p = 1000 kg/m3, D2 = 0.06 m, D1 = 0.02 m, b2 = 0.003 m, n = 0.5, R2 = 0.03 m, l = 0.02 m, d = 0.01 m, ^ = 0.01 Pa • s, A = 10~4 m, ftr = 0.1 m, J = 0.9 H • m, Ji = 0.0015 H • m • s, m = 0.2 kg, I = 10~5 H • m • s2, dc = 0.001 m, ¡ip = 0.2, flow rates through nozzles and ends for analytical calculation using the formula (2): ¡ic = 0.3, fit = 0.2 [15].
600 450
xn
«300 3
150
°0 0.1 0.2 0.3 0.4 0.5 0.6 0.7 0.8 0.9 1 600
450
l/i1
^,300 s
150
°0 0.1 0.2 0.3 0.4 0.5 0.6 0.7 0.8 0.9 1
f[s]
Fig. 7. The angular velocity in dependence on time
From Fig. 7 it can be seen that the rotor rotation reaches the operating angular speed rather quickly, in less than one second. At the same time, the angular velocity behaviours obtained on the basis of the expressions (2) and (3) are almost identical.
a)
r
f[s]
b)
r
10
5
lo
-5
"100 123456789 10 10
5
lo
-5
"100 123456789 10
t[ s]
Fig. 8. The speed of the rotor displacement along the y axis: (a) and (b) correspond to formulas (2) and (3), respectively
Fig. 8 shows significant fluctuations in the displacement velocity of the rotor along the y axis. In this case, the calculation using the equation (2) shows a smaller initial amplitude of oscillations (Fig. 8a), as well as a longer and smoother damping when passing to a stationary value, compared to the calculation when using modified formulas (3) and (Fig. 8b).
o
-0.25
-0.75
_10 1 2 3 4 5 6 7 8 9 10 0
-0.25
-0.75
"'o 1 2 3 4 5 6 7 8 9 10
t[s]
Fig. 9. Rotor displacement along the y axis: (a) and (b) correspond to (2) and (3), respectively
The results obtained on the basis of the equations (2) show gradual damping of the vertical oscillations (Fig. 9a). At the same time, the calculation using the modified equation (3) gives a faster decay of the oscillation amplitude along the y axis, as shown in Fig. 9b. Over time, the rotor asymptotically reaches the equilibrium rotation mode without oscillations.
a)
/ Vi u
<[s]
b)
, A Aaaa
fll/P /vw
1 "
Conclusion.
A mathematical model of the dynamics of a rotor in a low-flow centrifugal pump is built, taking into account the dynamic reaction of the liquid layer. Comparative characteristics of the "basic" calculation based on a purely analytical model and a "hybrid" calculation using the ANSYS Fluent package are given. Simulation of rotor dynamics with refined coefficients obtained on the basis of a single calculation in the ANSYS Fluent package shows a gradual, smooth transition from oscillations to a stationary operating mode. In this case, the calculation according to the approximate basic formulas shows a slow extinction of a non-stationary mode of operation, with an exit to constant small fluctuations. The resulting model can be used for preliminary calculation of pump operation modes, for optimal selection of design parameters in order to improve the dynamic properties of the rotor. At the same time, the proposed calculation method provides a sufficiently high accuracy of the results obtained without large computational costs.
Appendix 1. Dynamic layer response
Since the centrifugal pumps of spacecraft use low-viscosity fluids, we apply the Euler equation for the motion of an ideal fluid:
pd~L + pv w + VP = 0. (6)
Considering the effects associated with acceleration, we assume the velocity is small and neglect the convective term pVVV. In this case, the equation (6) takes the following form:
dV
Pit + VP = °. (7)
The equation (7) is supplemented with the continuity equation:
V- V = 0.
From the Euler and continuity equations, we obtain the equation for the pressure in cylindrical coordinates:
1 e_ ( BP \ dp dp
r dr \ dr J r2 84>2 dy2 On the surface of the rotor at r = R1 we have the boundary condition:
dVr dP
p~t + dr = 0. (9)
On a fixed surface with r = R2, the following condition is set:
£ = 0.
r
If the rotor has acceleration a in an arbitrary direction, then its radial component has the form:
dVr
= ax sin(^) - ay cos(^). (10)
Substituting the equation (10) into the equation (9), we find the pressure gradient along the normal on the rotor surface:
dP
— = p (ay cos(4) - ax sin(^j) .
Integrating equation (8) along r and using the boundary conditions on the surfaces, we obtain the averaged equation:
1 d2Ph d2Ph . . .
+ = P (ay COs(^) - ax sin(4))' (11)
where P is the mean pressure,
h = A — x sin(4) + y cos(4). For long bearings, we neglect the second term in equation (11):
d 2P h
d^2 = pR2 (ay cos(4) — ax sin(4)). (12)
By integrating equation (12) twice, we find the additional pressure associated with the acceleration of the rotor:
- pR2 (ay cos(^) — ax sin(4))
P = h . (13) Next, integrating equation (13), we determine the components of the additional force acting on the rotor due to its acceleration:
l>2n
Wx = PR sin(^)d^ = Knax + Kuay,
j 0
Wy = PR cos(^)d^ = K2iax + K22ay,
0
where
K = _ pR r2n sin2 4 d
11 A J0 1 — f cos(4) + f cos 4 ^ K = _ pR [2n cos2 4 4
22 a J 1 — f cos(4) + f cos 4 ^ K K pR3 i2n sin(4)cos(4) d4
k12 = k21 =--T— Z-x-, y-7 d4.
A Jo 1 — f cos(4) + f cos 4
References
[1] E.V.Morozova, Optimization of a centrifugal pump by criteria dimensions and cavitation characteristics, Gidravlika, 2(2018), 74-80 (in Russian).
[2] A.I.Petrov, N.Yu.Isaev, Hydrodynamic modeling of work centrifugal pump in the negative flow zone, Gidravlika, 3(2017), 91-101 (in Russian).
[3] E.S.Menon, Liquid Pipeline Hydraulics, NY. Basel: SYSTEK Technologies, Inc. Marcel Dekker, Inc, 2004.
[4] B.G.Kuharenko, Limiting cycle of self-oscillations of the rotor blades of a centrifugal compressor, Mechanical engineering and reliability problems, 2(2010), 10-14.
[5] D.S.Dolgin, A.E.Lebedev, I.S.Gudanov, Mathematical model of formation cavitation bubbles in a centrifugal pump, Engineering vestnik of Don's, 3(2019), 7-13 (in Russian).
[6] V.A.Pyhliy, A.A.Juravlev, E.V.Glushkova, Semi-analytical solution method boundary and initial problems in relation to the calculation of the rotor blades of centrifugal open type pumps, Power plants and technologies, 1(2019), 34-38.
[7] P.V.Pisarev, D.A.Ermakov, Yu.S.Kirova, O.A.Kashin, Numerical calculation of vibrating centrifugal pump blades loaded with fluid flow, Nauchno tekhnicheskij vestnik Povolzh'ya, 1(2018), 90-93 (in Russian).
[8] P.V.Pisarev, Numerical analysis of oscillatory processes in two-stage centrifugal pumps of various shapes, Mathematical modeling in natural sciences, 1(2017), 254-257.
[9] A.A.Lomakin, Centrifugal and axial pumps, Mashinostroenie, 1966 (in Russian).
[10] G.K.Borovin, A.I.Petrov, A.A.Protopopov, N.Yu.Isaev, Dynamics of rotors of oil-flow centrifugal pumps with hydrostatic bearings and driven by DC motors, Preprinty IPM im. M. V. Keldysha, (2016), no. 142 (in Russian).
[11] V.A.Voskresensiy, V.I.Diyakov, A.Z.Zile, Calculation and design of fluid friction supports: Handbook, Mashinostroenie, 1983 (in Russian).
[12] S.A.Chernavskiy, Plain bearings Mashgiz, 1963 (in Russian).
[13] V.M.Cherkasskiy, Pumps, fans, compressors, Energiya, 1977 (in Russian).
[14] V.O.Lomakin, A A.V.rtemov, A.I.Petrov, Determination of the influence of the main geometric parameters of the pump outlet NM 10000-210 on its characteristics, Science and education: electronic scientific and technical publication, (2012), no. 8, 71-83 (in Russian). DOI: 10.7463/0812.0445666
[15] G.K.Borovin, A.A.Protopopov, Calculation of the optimal axial clearance of a half-open impeller of a centrifugal low-flow pump of the thermal control system of a spacecraft, Preprinty IPMim. M.V. Keldysha, (2013), no. 086 (in Russian).
Численно-аналитическое моделирование работы центробежного насоса
Виктор А. Иванов
Сибирский федеральный университет Красноярск, Российская Федерация
Николай В. Еркаев
Институт вычислительного моделирования СО РАН Красноярск, Российская Федерация
Аннотация. Построена математическая модель нестационарных режимов вращения ротора центробежного насоса на основе предварительного расчета в ANSYS Fluent и последующего аналитического расчета с использованием уточненных параметров.
Ключевые слова: математическое моделирование, гидростатический подшипник, центробежный насос.
