Computational Approaches in Industrial Centrifugal Pumps

Abstract

The growing energy demand is expected to be met with increased oil and gas production. Hence, there is a need to design high-performance industrial centrifugal pumps. Recent improvements in CFD are considered as a valuable research tool to investigate the flow inside the pump and its influence on the performance of the centrifugal pump. The scope of the chapter is to emphasize the use of CFD and theoretical analysis for design and to show the prospect of improving the efficiency of a centrifugal pump. The chapter discusses the computational approaches to the CAD modeling and CFD simulation of the industrial centrifugal pumps, and the strategies and methodologies adopted. The chapter would be relevant and useful to both the pump designers, manufacturers, and industrial users.

1. Introduction

Centrifugal pump has wide application areas as it is being used in the sectors like agriculture, oil and natural gas, public water supply, sanitation, domestic, household utilities, petroleum refining, petrochemicals, mining, etc. That is why the centrifugal pump has a significant effect on the nation’s economy. The present chapter describes the steps in detail that include the design of industrial centrifugal pumps, construction of geometrical models of various components of centrifugal pumps, grid generation, governing equations, boundary conditions, discretization, CFD solver settings, and the importance of CFD in improving pump design and performance.

Despite wide applications, the design of the pump and its performance prediction using a conventional approach involves trial and error which is a time-consuming task. CFD is a robust tool for prediction of performance of the centrifugal pumps that can be used to analyze and design the pump for performance and cavitation prediction [1], and condition monitoring [2] with a better understanding of flow physics. CFD analysis of any problem of engineering interest consists of three major steps, namely, pre-processing, simulation, and post-processing.

2. Pre-processing for centrifugal pump

Pre-processing consists to define the domain of interest (i.e., creation of geometry), then dividing it into a finite number of elements known as computational grid or mesh generation of the flow domain, and formulation of the problem.

2.1 Geometrical modeling of a centrifugal pump

Modeling is a process to create the complex engineering model of an industrial machine or device that can be shaped into the desired model by assembling the different parts (called solid volumes) and thus the final solid model is a virtual replica of the actual product that can be functioned (for example, rotated) like a real product. Therefore, such a complex geometrical machine or device can be modeled (constructed) by any computer-aided design (CAD) software, such as CATIA, Solid works, ProE, AutoCAD, and also with the help of non-CAD software, like BladeGen, Ansys Design Modeler (DM), etc. Modeling of an industrial centrifugal pump is always a tedious task for engineers as it deals with complex and intricate shapes associated with the impeller and volute geometries [3, 4].

A two-dimensional geometry of a single-suction, single-stage industrial centrifugal pump having a double-volute casing and impeller with three blades is generated as per the drawings made available with the help of CAD software- CATIA. The two-dimensional geometry of the impeller and double volute of an industrial pump is then converted into three-dimensional geometry for further processing. The major dimensions of the centrifugal pump are furnished in Figure 1 shows the impeller and volute casing as separate entities, assembled and cut section of the centrifugal pump.

Figure 1.

Impeller and casing assembly with cut-section of the centrifugal pump.

The three-dimensional pump geometry constructed in CATIA is imported to Ansys Design Modeler (DM) through Solid works using appropriate format to avoid loss of surface in the process and to ensure the accurate geometry. Next, the flow domain is extracted from the imported geometry after using capping option available in Ansys DM. The extracted flow domain requires cleanup before the meshing. A flow chart explaining the steps involved in the complete process is shown in Figure 2.

Figure 2.

Steps in the import of the geometry and extraction of the flow domain.

Wireframe of the flow domain after cleanup is shown in Figure 3 while parts of the final flow domain are shown in Figure 4 which is used for further processing.

Figure 3.

Wireframe of the flow domain of the centrifugal pump after cleanup.

Figure 4.

Flow domain of the centrifugal pump.

A pump impeller consists of an array of blades (or vanes) arranged in a certain fashion and covered with the shrouds on either side. In order to study the effect of blade geometry on the performance of the centrifugal pump, pump impellers are constructed with six different blade outlet angles. The baseline impeller consists of three backward curved-twisted blades with a 16° outlet angle. The blade thickness is varied from 6 mm at the leading edge to 3 mm at the trailing edge with a wrap angle of 278°. The construction of the impeller blade with different blade outlet angle in CATIA software is a daunting and time-consuming task. Hence, many researchers [5, 6, 7, 8] relied on BladeGen tool available in the Ansys Blade Modular software in modeling the blade geometry of turbomachines. A two-dimensional meridional view of the impeller blade is shown in Figure 5, which is referred to as baseline profile. The top layer of the profile is termed as the shroud layer, while the bottom layer is known as the hub layer. The liquid flows around the blade profile radially and discharges at the outlet. Three-dimensional profile of a single blade of impeller is shown in Figure 6a, while the geometry of three blades of the impeller is shown in Figure 6b. The geometry of the blades thus created is imported from BladeGen to Ansys-Workbench. The extracted flow domain of a three bladed impeller is shown in Figure 7a.

Figure 5.

Two-dimensional master profile of the impeller blade.

Figure 6.

Generation of blade profile of the centrifugal pump.

Figure 7.

Flow domain of and meshing of three bladed impeller.

2.2 Generation of computational grid

Grid generation (or meshing) is the next step required to perform after extraction of flow domain. It is the process of dividing the given flow domain into a number of small domains or parts called control volumes. It is seen from the literature survey that most of the researchers [8, 9, 10, 11, 12, 13, 14, 15, 16, 17, 18, 19] used unstructured tetrahedral elements for grid generation of the centrifugal pumps. Some of the researchers [20, 21] used both structured and unstructured element, while others [22] used combined unstructured hexahedral and tetrahedral elements for centrifugal pump. Zhang et al. [23] used tetrahedral grids for meshing double-volute casing while impeller was meshed with hexahedral grids.

Due to complex profile of the impeller and the volute casing, unstructured tetrahedral grid is used in the present study for all flow domains of the centrifugal pump using grid generation software available in Ansys-Workbench. Patch independent method has been used, which performs meshing from surface to volume, to maintain the quality of the grid to a great extent. The accuracy of the computational results depends upon the grid density and quality. Therefore, a high-density and high-quality grid is desirable to capture the complex flow phenomena, like flow separation and recirculation in the diffusing passages accurately [24]. However, high grid density increases computational time and memory. Therefore, a tradeoff between the computational cost and the grid density is required.

Grid independence test is performed to ensure the independency of solution with the grid size. This test is performed by refining the grid using patch independent algorithm. The number of elements are gradually increased from 1.6 million to 5.2 million in the present study to carry out the grid independence test. The results of the grid independence test are presented in Table 1. The results show that no significant change in the computed value of the head is observed from 4.2 to 5.2 million grid size. Therefore, in order to save computational time, the grid with around 4.16 million elements is chosen for further simulation. The corresponding average element size is 3.86 mm and the convergence time is 181 minutes.

Grid size	Average element size (mm)	CPU time (min.)	No of iteration	Head (m)
1,628,749	5.28	52	245	215.72
2,162,341	4.9	67	347	214.91
2,833,632	4.39	84	458	213.62
3,682,491	4.02	135	667	213.08
4,164,202	3.86	181	804	210.55
4,725,137	3.7	224	976	210.46
5,152,315	3.56	258	1014	210.56

Table 1.

Grid independence study.

Apart from size of element, the quality of a computational grid is also determined by the shape of the individual cells in the computational flow domain. The skewness of the element describes the shape of the element. CFD results are greatly affected by skewness and are therefore considered as one of the main measures for grid quality. Ansys-CFX provides a general guideline to determine the grid quality based on the value of skewness as reproduced in Table 2 for ready reference.

Value of Skewness	0.0–0.25	0.25–0.50	0.50–0.75	0.75–0.90	0.90–0.99	1.00
Grid quality	Excellent	Good	Acceptable	Poor	Bad	Degenerate

Table 2.

Range of skewness and cell quality.

It is not always easy to reduce the skewness of all the elements in the acceptable limit as indicated by Table 1 in the complex flow domain as in the case of centrifugal pump. The literature [17, 23, 25, 26] also reports maximum skewness of centrifugal pumps’ grid between 0.82 to 0.95, indicating poor to bad quality of the grid as per Table 1. However, it is not only the high skewness of the elements but their number and location are also important. It is found that if the highly skewed elements are few and are away from the region of interest, they do not bring a significant effect on the computational results. Histogram in Figure 8 shows the frequency distribution of the skewness of the elements corresponding to 4.2 million elements in the present study. The histogram shows that insignificant number of elements are present above 0.5 skewness and almost negligible above 0.75 skewness.

Figure 8.

Histogram of number of elements versus skewness.

The maximum skewness reported in the present study is 0.852 (affecting only 29 elements) and their locations are shown in Figure 9. It is observed in the present study that their presence does not influence the computational results. The final grid corresponding to 4.2 million elements is thus selected for further simulation as shown in Figures 10 and 11.

Figure 9.

Most skewed elements in the computational domain of the centrifugal pump.

Figure 10.

Grid in the different flow domains of the centrifugal pump.

Figure 11.

(a) Grid in the impeller and casing of centrifugal pump.

2.3 Governing equations

Various governing equations used in the present study to simulate the performance of the centrifugal pump are discussed in this section.

2.3.1 Continuity equation

$$\frac{\partial \rho }{\partial t}+\nabla •\left(\rho \mathit{u}\right)=0$$E1

2.3.2 Momentum equations

$$\frac{\partial \left(\rho \mathit{u}\right)}{\partial t}+\nabla •\left(\rho \mathit{u}\otimes \mathit{u}\right)=-\nabla p+\nabla •\mathit{\tau }+S_M$$E2

where τ is the stress tensor related to the strain rate by.

$$\mathit{\tau }={\mu }_{\mathit{eff}}\left(\nabla \mathit{u}+{\left(\nabla \mathit{u}\right)}^T-\frac{2}{3}\delta \nabla •\mathit{u}\right)$$E3

where S_M is the sum of body forces. For the study of flow in a rotating frame of reference having constant angular velocity ω, additional sources of momentum are required to include the effect of the Coriolis force and that of the centrifugal force:

$$S_{M,\mathit{rot}}=S_{\mathit{Cor}}+S_{\mathit{cfg}}$$E4

where,

$$S_{\mathit{Cor}}=-2\mathit{\rho \omega }\times \mathit{u}$$E5

$$S_{\mathit{cfg}}=-\mathit{\rho \omega }\times \left(\omega \times r\right)$$E6

where r is the location vector and u is the relative frame velocity.

μ_eff is the effective viscosity accounting the turbulence through turbulence model described in next section.

For the simulation of steady case, time derivative term will be zero in Eqs. (1) and (2).

2.3.3 Turbulence model

Various turbulence models are available to capture the turbulence behavior of the flow. Thus, the selection of appropriate turbulence model for CFD analysis of the centrifugal pump is a difficult task. The literature review reveals that most researchers used two-equation turbulence models in turbo machinery application as they offer a good compromise between computational accuracy and computational effort. A number of researchers [2, 5, 6, 9, 12, 20, 21, 22, 23, 24, 27, 28] used the standard k-ε turbulence model. Shahin et al. [27] and Jafarzadeh et al. [20] have used renormalized group (RNG) k-ε turbulence model. Siddique et al. [19], Deshmukh et al. [8], Paul et al. [28] and Bhattacharyya et al. [29, 30, 31, 32, 33, 34, 35] used the shear stress transport (SST) turbulence model. Alemi et al. [25] used the standard k-ε, low-Re k-ω, and SST turbulence models. Thus, it is evident that the standard k-ε turbulence model is widely used as compared to any other turbulence model for the simulation of turbomachines and other simulations involving complex geometries. Few researchers used the RNG k-ε model as an alternative to the standard k-ε model. Ansys CFX-solver theory guide suggests that the k-ω based SST model is more accurate for the prediction of the onset and the magnitude of flow separation under adverse pressure gradients by the inclusion of transport effects into the formulation of the eddy-viscosity.

Thus, four well-known two-equation turbulence models, namely standard k-ε, RNG k-ε, Wilcox k-ω, and SST turbulence models, are tested in the present study to select the one for the present study. Keeping all other conditions same, the results for these four turbulence models are obtained and compared with experimental results in Table 3.

Inlet boundary condition	Outlet boundary condition	Turbulence model	Head (m)	Output power (kW)	Input power (kW)	Efficiency (%)
Total pressure (Pa)	Mass flow rate (kg/s)	Standard k-ε	210.55	138.84	176.27	78.76
		RNG k-ε	210.49	138.80	176.15	78.79
		k-ω	203.53	137.19	174.00	78.84
		SST k-ω	204.62	134.93	172.57	78.18
Experimental results			206.54	136.65	185.11	73.82

Table 3.

Comparison of turbulence models for Q = 242.8 m³/h at 2933 rpm.

It is found that predicted result by k-ω model is in more agreement with the experimental results compared to the other turbulence models. But the computation time is high as compared to others. Since standard k-ε turbulence model has taken minimum computation time and its results are also in reasonable agreement with the experimental ones, standard k-ε turbulence model is selected for further simulation in the present study. The standard k-ε model used in the present study is mentioned in the following equations. Since the standard k-ε model is based on the eddy viscosity concept,

$${\mu }_{\mathit{eff}}=\mu +{\mu }_t$$E7

where μ_t is the turbulence viscosity. The k-ε model assumes that the turbulence viscosity (μ_t) is linked to the turbulence kinetic energy (k) and dissipation rate (ε) via the relation.

$${\mu }_t=C_{\mu }\rho \frac{k^2}{\epsilon }$$E8

$$\frac{\partial \left(\mathit{\rho k}\right)}{\partial t}+\frac{\partial }{\partial x_j}\left(\rho U_{j}k\right)=\frac{\partial }{\partial x_j}\left[\left(\mu +\frac{{\mu }_t}{{\sigma }_k}\right)\frac{\partial k}{\partial x_j}\right]+P_k-\mathit{\rho \epsilon }+P_{\mathit{kb}}$$E9

$$\frac{\partial \left(\mathit{\rho \epsilon }\right)}{\partial t}+\frac{\partial }{\partial x_j}\left(\rho U_j\epsilon \right)=\frac{\partial }{\partial x_j}\left[\left(\mu +\frac{{\mu }_t}{{\sigma }_{\epsilon }}\right)\frac{\partial \epsilon }{\partial x_j}\right]+\frac{\epsilon }{k}\left(C_{\epsilon 1}{P}_k-C_{\epsilon 2}\mathit{\rho \epsilon }+C_{\epsilon 1}{P}_{\mathit{\epsilon b}}\right)$$E10

where model constants are $C_{\mu }=0.09,{\mathrm{C}}_{\epsilon 1}=1.44,{\mathrm{C}}_{\epsilon 2}=1.92,{\sigma }_{\mathrm{k}}=1.0\text{and}{\sigma }_{\epsilon }=1.3$.

P_kb and P_εb represent the influence of the buoyancy forces associated with k and ε, respectively.

P_k is the turbulence production term due to viscous force, which is modeled for incompressible flow as.

$$P_k={\mu }_t\left(\frac{\partial U_i}{\partial x_j}+\frac{\partial U_j}{\partial x_i}\right)\frac{\partial U_i}{\partial x_j}-\frac{2}{3}\frac{\partial U_k}{\partial x_k}\left(3{\mu }_t\frac{\partial U_k}{\partial x_k}+\mathit{\rho k}\right)$$E11

Scalable wall function is used to capture near wall flow physics in k-ε turbulence model. Judicious choice of y⁺ value affects the quality of computed results as well as computation time [36]. This approach limits the value of y⁺ used in the logarithmic formulation as follows:

$$y^{+}=\mathit{max}\left(y,*,,,h_s^{+},/,2,,,11.06\right)$$E12

This approach also incorporates the effect of wall surface roughness as dimensionless sand-grain roughness, h_s⁺.

2.3.4 Turbulence intensity

The turbulence intensity is defined as the ratio of root mean square of fluctuating velocity (u′) to the mean flow velocity (u_avg.)_. For internal flow, it is found that turbulence intensity less than 1% represents low turbulence and turbulence intensity greater than 10% represents high turbulence. Thus, the allowable range of turbulence intensity is from 1–10% corresponding to very low to very high levels of turbulence in the internal flow. The turbulence intensity in the fully developed duct flow can be calculated from the following formula:

$$I=\frac{u^{\text{'}}}{u_{\mathit{avg}.}}=0.16{\left({Re}_{D_H}\right)}^{-1/8}$$E13

Nominal turbulence intensity ranges from 1–5% depending up on the specific application. In absence of experimental data, the default turbulence intensity value of 3.7% is considered good estimate for nominal turbulence through a circular inlet.

It is found that many researchers [25, 37, 38, 39, 40, 41] used 5% turbulence intensity at the inlet boundary conditions for turbomachines, while Shahin et al. [27] used 2% turbulence intensity, Jafarzadeh et al. [20] used 1% turbulence intensity, Murmu et al. [42] and Alam et al. [43] investigated varied turbulence intensity at the inlet boundary condition.

Three standard turbulence intensity, 1%, 5%, and 10%, are tested in the present study for CFD simulation of flow through the centrifugal pump. The results are presented in Table 4 for comparison. The results show almost negligible variation in the value of head, output power, input power, and efficiency of centrifugal pump. Hence turbulence intensity equal to 5% is finally chosen for further study.

Turbulence intensity (%)	Head (m)	Output power (kW)	Input power (kW)	Efficiency (%)
1.0	210.5	138.83	176.26	78.76
5.0	210.5	138.83	176.26	78.76
10.0	210.5	138.84	176.27	78.77

Table 4.

Comparison of turbulence intensity for Q = 242.8 m³/h at 2933 rpm.

2.3.5 Boundary conditions

While solving the Navier–Stokes equations and continuity equation, an appropriate boundary conditions need to be applied. Boundary conditions are a very important set of the properties and conditions applied on the surfaces of computational domains. These are required to define the flow field completely. The solid and fluid regions are generally represented by domains or cell zones while internal surfaces and boundaries are represented by the face zones. The data corresponding to boundary conditions are then assigned to these face zones. Each variable equation requires the meaningful values at the boundaries at the domain in order to generate the values in entire domain. Boundary conditions used by various researchers working in computational study of centrifugal pump are summarized in Table 5.

Authors	Year	Inlet	Outlet
Blanco-Marigorta et al. [26]	2000	Flow rate	Static pressure.
Stickland et al. [9]	2000	Flow rate	Static pressure.
Gonzalez et al. [10]	2002	Total pressure	A variable static pressure proportional to the kinetic energy at the outlet for their unsteady pump.
Majidi [37]	2005	Mass flow rate	For all variables (with exception of pressure) a zero-gradient condition.
Gonzalez and Santolaria [13]	2006	Total pressure	A pressure drop proportional to the kinetic energy at the outlet.
Cheah et al. [14]	2007	Mass flow rate	Mass flow is being specified.
Cheah et al. [38]	2008	Total pressure	Mass flow rate.
Spence and Teixeira [44]	2009	Mass flow rate	Static pressure.
Feng et al. [45]	2010	Total pressure	Mass flow rate.
Perez et al. [16]	2010	Relative Pressure of 0 Pa	Mass flow rate.
Houlin et al. [17]	2010	Velocity	Outflow.
Li [46]	2011	Velocity	Mass flow rate
Jafarzadeh et al.[20]	2011	Volume flow rate	Pressure.
Hedi et al. [47]	2012	Total pressure of 1 atm	Mass flow rate.
Kumar et al. [22]	2013	Maas flow rate	Static pressure.
Bellary and Samad [5]	2013	Pressure	Mass flow rate.
Lei et al. [48]	2014	Velocity	Pressure.
Kim et al. [40]	2014	Pressure of 0 atm	Mass flow rate.
Li [41]	2014	Absolute velocity	Zero Static pressure.
Bellary and Samad [49]	2015	Pressure	Mass flow rate.
Alemi et al. [25]	2015	Mass flow rate	Pressure.
Siddique et al. [19]	2017	Static pressure	Mass flow rate.
Deshmukh et al. [8]	2017	Static pressure of 1.01 kPa	Mass flow rate.

Table 5.

Boundary conditions used by researchers for centrifugal pump.

After analyzing all the boundary conditions, it is found that the computational results of centrifugal pump corresponding to the boundary condition of total pressure at inlet and mass flow rate matches most closely with the experimental measurement. The inlet total pressure is known from the measured data corresponding to different discharge. The turbulence kinetic energy and turbulence dissipation rate at the inlet totally depends on the upstream history of the flow.

The inlet pipe and casing are considered as stationary reference frames while impeller is considered as rotating reference frame. These frames are coupled through frozen-rotor interface model. There are two interfaces, one connecting pipe and impeller and other connecting impeller and casing. Walls associated with impeller are considered as rotating walls while those associated with pipe and casing are considered as stationary walls. No-slip boundary conditions are applied over the impeller blades and casing walls. A list of boundary conditions used for the CFD simulation of the pump in the study is given in Table 6.

Parameters	Description
Inlet	Total Pressure
Outlet	Mass flow rate
Blades	Rotating wall
Pipe	Stationary reference frame
Impeller	Rotating reference frame
Casing	Stationary reference frame
Interface-1 (Between pipe outlet and impeller inlet)	Fluid–Fluid
Interface-2 (Between impeller outlet and casing inlet)	Fluid–Fluid

Table 6.

Boundary conditions used.

2.4 CFD solver settings

The governing equations along with boundary conditions are discretized using second order upwinding scheme. The discretized equations result in a set of algebraic equations. Ansys-CFX, an element based finite volume method (EbFVM) solver with a cell vertex formulation is used for the present study to solve the algebraic equations. Table 7 shows various solver settings used during the CFD simulation.

2.4.1 Convergence criteria

Ansys-CFX uses a multigrid (MG) accelerated incomplete lower upper (ILU) factorization technique for solving the discretized system of equations. This technique is iterative in nature. The measure of convergence of iterative solution is the residual which quantifies the error in the computational solution. In a CFD analysis, the residual quantifies the local imbalance of each conservative control volume equation. Thus, after specifying all the boundary conditions, numerical schemes, convergence criteria are specified to solve the problem.

Researchers used various convergence criteria for numerical simulation of centrifugal pumps. Feng et al. [45] taken convergence criteria 10⁻³. Majidi [37], Cheah et al. [14], Perez et al. [16], Hedi et al. [47], Li [46], Alemi et al. [25], Siddique et al. [19] used 10⁻⁴ residual, while Gonzalez et al. [10] used a residual of 10⁻⁵ for mass and momentum and 10⁻⁴ for k-ε turbulence model. Gonzalez et al. [9], Houlin et al. [17], Bellary and Samad [5, 28], Deshmukh et al. [8] have considered the residual at 10⁻⁵.

The computational effort and computational accuracy are always competing parameters. Thus, the accuracy of the computation in terms of root mean square (RMS) of residual up to10⁻⁶ is calculated with the corresponding computational effort for two different turbulence models, k-ε and k-ω, at a flow rate of 242.8 m³/h and rotation speed of impeller equal to 2933 rpm. Keeping computation time in consideration, the convergence criteria for further run of flow simulation is set to 10⁻⁵ for residuals of mass and momentum equations.

2.4.2 Steady state flow simulation

Moving reference frame (MRF) allows the computational analysis of cases involving domains that are rotating relative to one another. Since MRF is based on the general grid interface (GGI) technology, this feature allows rotor/stator interaction in the investigation of turbomachines in Ansys-CFX. In GGI, a control surface approach is used to connect across the interface or periodic condition. This intersection algorithm allows the complete freedom to change the physical distribution and the grid topology across the interface. This permits the use of most appropriate meshing style for each component involved in the analysis. Many researchers [14, 15, 18, 20, 27, 39, 41, 45, 50] carried out computational investigation of centrifugal pump to assess its performance using MRF approach in Ansys-Fluent and Ansys-CFX.

The steady state CFD simulation is carried out using MRF approach considering the volute and pipe sections as stationary frame while impeller flow field as rotating frame. In this case, the relative position between volute casing and impeller flow field is considered fixed in time and space. The grids of pipe section, impeller, and volute that are generated separately with the help of tetrahedral elements are connected by means of a “frozen-rotor” interface. The interface treatment is fully implicit that preserves the flow field variation across the interface and does not adversely affect convergence of the overall solution.

2.4.3 Unsteady flow simulation

Transient/unsteady behavior can be caused by changing the boundary conditions of the flow or can be related to the flow where steady state condition is never reached, even when all other aspects of the flow conditions are not changing. For the unsteady simulation of the present problem, the surface fluxes at each side of the interface are first computed at the start of each time step at the current relative position. The dissimilar meshes at the pipe, impeller, and volute interfaces are connected by GGI connections to permit transient flow interaction between a stator and rotor passages across the sliding (transient rotor-stator/frame change) interface. This approach allows the accounting of all interaction effects between the components that are in the transient relative motion to each other across the interface. The interface position is updated at each time step as the relative position of the grid’s changes on each side of the interface. To determine the real time information, the transient simulations require appropriate time interval at which the Ansys-CFX solver calculates the flow field.

A few researchers [10, 22, 26, 27, 37, 42] also studied the transient CFD simulation. Three different time step sizes (3°, 6°, 9°) are compared for Q = 280 m³/h at 2933 rpm in the present study to choose the appropriate time step size which provides the necessary time resolution as shown in Table 8. From the study, it is found that 6° impeller rotation per step is enough to reduce the maximum residuals to 10⁻⁵. Thus, impeller rotation step size of 6° is chosen for further study corresponds to 60 steps/rotation and 360 steps to complete 6 full rotations of the impeller corresponding to the pump running time equal to 0.12274 s.

Solver	Coupled
Analysis method	Multiple reference frame (MRF) or Frozen rotor method
	Steady/Unsteady state analysis
Reference Pressure	1 atm

Table 7.

Solver setting options used in Ansys-CFX.

Impeller rotation/step	Total no of time steps/6 revolutions	Time taken/revolution (ms)	Time per step size (ms)	Total time taken for 6 revolution (ms)
9°	240	20.46	0.51	122.74
6°	360	20.46	0.34	122.74
3°	720	20.46	0.17	122.74

Table 8.

Comparison of impellerrelative positions at Q = 280 m³/h at 2933 rpm.

2.4.4 Cavitation simulation using Rayleigh–Plesset equation

In the present study, the steady state numerical simulation with cavitation is carried out with a multiple frames of reference (MRF) approach using the Rayleigh Plesset equation to investigate the cavitation in the centrifugal pump. The Rayleigh–Plesset equation provides the basis for the rate equation controlling vapor generation and condensation. The Rayleigh–Plesset equation that describes the growth of a gas bubble in a liquid is derived from a mechanical balance, assuming no thermal barriers to bubble growth. The equation is given as:

$$R_B\frac{d^2{R}_B}{{\mathit{dt}}^2}+\frac{3}{2}{\left(\frac{{\mathit{dR}}_B}{\mathit{dt}}\right)}^2+\frac{2\sigma }{{\rho }_f{R}_B}=\frac{p_v-p}{{\rho }_f}$$E14

where R_B is the bubble radius, p_v is the pressure in the bubble (assumed to be equal to the vapor pressure at the temperature of the liquid), p is the liquid pressure surrounding the bubble, ρ_f is the density of liquid, and σ is the coefficient of surface tension between the liquid and vapor. Neglecting the surface tension and the second order terms (which is acceptable for low oscillation frequencies), the equation reduces to:

$$\frac{{\mathit{dR}}_B}{\mathit{dt}}=\sqrt{\frac{2}{3}\frac{p_v-p}{{\rho }_f}}$$E15

The rate of change of bubble volume is:

$$\frac{{\mathit{dV}}_B}{\mathit{dt}}=\frac{d}{\mathit{dt}}\left(\frac{4}{3}\pi R_B^3\right)=4\pi R_B^2\sqrt{\frac{2}{3}\frac{p_v-p}{{\rho }_f}}$$E16

thus, the rate of change of bubble mass is:

$$\frac{{\mathit{dm}}_B}{\mathit{dt}}={\rho }_g\frac{{\mathit{dV}}_B}{\mathit{dt}}=4\pi R_B^2{\rho }_g\sqrt{\frac{2}{3}\frac{p_v-p}{{\rho }_f}}$$E17

The volume fraction r_g may be expressed as:

$$r_g=V_B{N}_B=\frac{4}{3}\pi R_B^{3}N{}_B$$E18

where N_B is number of bubbles per unit volume,

and the total inter-phase mass transfer rate per unit volume is:

$${\dot{m}}_{\mathit{fg}}=N_B\frac{{\mathit{dm}}_B}{\mathit{dt}}=\frac{3r_g{\rho }_g}{R_B}\sqrt{\frac{2}{3}\frac{p_v-p}{{\rho }_f}}$$E19

This expression has been derived assuming bubble growth due to vaporization. The expression can be generalized to include condensation as follows:

$${\dot{m}}_{\mathit{fg}}=F\frac{3r_g{\rho }_g}{R_B}\sqrt{\frac{2}{3}\frac{p_v-p}{{\rho }_f}}sgn\left(p_v-p\right)$$E20

where F is an empirical factor that is designed to account for rates of condensation and vaporization. Vaporization is usually much faster than condensation. Thus, F has different values for vaporization and condensation. Despite the fact that Eq. (20) has been generalized for vaporization and condensation, it requires further modification in the case of vaporization. Vaporization is initiated most commonly at the nucleation sites for non-condensable gases. Thus, the bubble radius R_B will be replaced by the nucleation site radius R_nuc for the purpose of modeling of vaporization. The nucleation site density decreases as the vapor volume fraction increases as there is less liquid available. For the case of vaporization, r_g in Eq. (20) is replaced by $r_{\mathit{nuc}}\left(1-r_g\right)$to give mass transfer rate as:

$${\dot{m}}_{\mathit{fg}}=F\frac{3r_{\mathit{nuc}}\left(1-r_g\right){\rho }_g}{R_{\mathit{nuc}}}\sqrt{\frac{2}{3}\frac{\left|p_v-p\right|}{{\rho }_f}}sgn\left(p_v-p\right)$$E21

where r_nuc is the nucleation site volume fraction. Eq. (21) is maintained in the case of condensation.

To obtain an inter-phase mass transfer rate, the values of bubble concentration and radius are required. The Rayleigh–Plesset cavitation model uses the following default values for the model parameters as implemented in Ansys-CFX:

$$R_{\mathit{nuc}}=1\mathit{\mu m};r_{\mathit{nuc}}=5E-4;F_{\mathit{vap}}=50;F_{\mathit{\text{cond}}}=0.01$$

2.4.5 Acoustics simulation

Noise in the pump is produced primarily by the unsteady pressure field. To predict mid- to far-field noise, the methods based on Lighthill’s acoustic analogy offer sustainable alternative to the direct method. The acoustic analogy decouples the generation of sound from its propagation. Thus, enabling the separation of the flow solution process from the acoustics analysis. Therefore, in this approach, the near-field flow is obtained from the appropriate flow governing equations such as unsteady RANS equations. Then the sound is predicted with the help of analytically derived integral solutions of wave equations using determined flow.

2.4.6 The Ffowcs-Williams and Hawkings Equation for Acoustic Study

Ansys-Fluent offers an acoustic study method based on the Ffowcs-Williams and Hawkings (FW-H) formulation. The Ffowcs-Williams and Hawkings (FW-H) equation is essentially an inhomogeneous wave equation that is derived by manipulating the Navier–Stokes equations and the continuity equation. The FW-H formulation uses the most general form of Light hill’s acoustic analogy neglecting the interaction between fluid and solid. This formulation predicts sound generated by equivalent acoustic sources. ANSYS-Fluent uses a time domain integral formulation to compute time histories of sound pressure, or acoustic signals, at prescribed receiver locations. The FW-H equation can be written as:

$$\begin{array}{l}\frac{1}{a_0^2}\frac{{\partial }^2{p}^{\text{'}}}{\partial t^2}-{\nabla }^2{p}^{\text{'}}=\frac{{\partial }^2}{{\partial }_i{\partial }_j}\left\{T_{\mathit{ij}}H\left(f\right)\right\}-\frac{\partial }{\partial x_i}\left\{\left[p_{\mathit{ij}}{n}_j+\rho u_i\left(u_n-v_n\right)\right]\delta \left(f\right)\right\}\\ +\frac{\partial }{\partial t}\left\{\left[{\rho }_0{v}_n+\rho \left(u_n-v_n\right)\right]\delta \left(f\right)\right\}\end{array}$$E22

where

$u_i$= fluid velocity component in the$X_i$direction

$u_n$= fluid velocity component normal to the surface$f=0$

$v_i$= surface velocity components in the $X_i$direction

$v_n$= surface velocity component normal to the surface

$\delta \left(f\right)=$ Dirac delta function

$H\left(f\right)=$ Heaviside function

$p^{\text{'}}=p-p_0$ is the sound pressure at the far field.

In order to facilitate the use of generalized function theory and the free space Green function to obtain the solution, a mathematical surface represented by f = 0 is introduced to ‘embed’ the exterior flow problem (f > 0) in an unbounded space. The surface (f = 0) corresponds to the source (emission) surface and can be made to coincide with a body (impermeable) surface or off the body permeable surface.

n_j is the unit normal vector pointing toward the exterior region (f > 0),

a₀ is the far-field speed of sound,

T_ij is the Lighthill stress tensor, defined as

$$T_{\mathit{ij}}=\rho u_i{u}_j+P_{\mathit{ij}}-a_0^2\left(\rho -{\rho }_0\right){\delta }_{\mathit{ij}}$$E23

For a Stokesi an fluid, the compressive stress tensor, P_ij, is given by

$$P_{\mathit{ij}}=p{\delta }_{\mathit{ij}}-\mu \left[\frac{\partial u_i}{\partial x_j}+\frac{\partial u_j}{\partial x_i}-\frac{2}{3}\frac{\partial u_k}{\partial x_k}{\delta }_{\mathit{ij}}\right]$$E24

The subscript 0 denotes the free-stream properties.

2.4.7 Proudman’s formula

Proudman [51] derived a formula for acoustic power generation by isotropic turbulence without mean flow using Light hill’s acoustic analogy. For a given turbulence field, the Proudman’s formula yields an approximate measure of the local contribution to total acoustic power per unit volume. Proudman’s original derivation neglected the retarded time difference. Later, the formula was re-derived by accounting for the retarded time difference. Both derivations provide acoustic power due to the unit volume of isotropic turbulence in (W/m³) as

$$P_A={\mathit{\alpha \rho }}_0\left(\frac{u^3}{\mathrm{\ell }}\right)\frac{u^5}{a_0^5}$$E25

where u is the turbulence velocity,

$\mathrm{\ell }$₌ turbulence length scale, a₀₌ speed of sound, α = model constant.

Eq. (25) can be rewritten in terms of k and ε as.

$$P_A={\alpha }_{\epsilon }{\rho }_0\epsilon M_t^5$$E26

where

$$M_t=\frac{\sqrt{2k}}{a_0}$$E27

${\alpha }_{\epsilon }$ is the rescaled constant, set to 0.1

In the present study, the FW-H model is used to measure acoustics in terms of the level of noise radiation induced by the inner flow field in the centrifugal pump encompassing the impellers and volute. In FW-H acoustic model, it is essentially required that all the receivers are located far away (at least at a distance of 1 m) from the primary sources of sound. Thus, for monitoring, 24 receivers, indicated by P1-P24, are arranged in circular surface at an angular gap of 15° is shown in Figure 12.

Figure 12.

Arrangement of sound receivers outside the centrifugal pump.

Receiver P7 is placed in line with the volute tongue and receiver P19 is placed in line with the volute splitter. To extract the relevant acoustic spectra, a Fourier transformation is applied connecting the acoustic signals found at these receivers. The sound pressure level (SPL) in decibel (dB) is calculated by:

$$\mathit{SPL}=20{log}_{10}\left(\frac{P_e}{P_{\mathrm{Ref}}}\right)$$E28

where P_Ref is the reference sound pressure = 2 × 10⁻⁵ Pa for air,

P_e is the effective acoustic pressure defined by:

$$P_e=\sqrt{\frac{1}{T}}\underset{0}{\overset{T}{\int }}p{\text{'}}^2\mathit{dt}$$E29

It is necessary to derive a temporal intensity profile to obtain the intensity of noise radiation, involving a superposition of acoustic pressures at each Fourier frequency. In this regard, the total sound pressure level(TSPL) is introduced and expressed as:

$$\mathit{\text{TSPL}}=10lg\sum _{i=1}^n{10}^{{\mathit{SPL}}_i/10}$$E30

Detailed analysis of cavitation detection in a centrifugal pump impeller using acoustic generation can be found in Jaiswal et al. [52].

3. Conclusions

Intensive computation coupled with geometric modeling helps the engineers and industrial users of the centrifugal pumps to investigate its performance characteristics [53] and cavitation prediction while dealing with multiple working fluids [54] and geometric modification of blades and its number [55, 56]. The present chapter discussed the computational strategies that encompasses the CAD modeling, computational grid generation, steps involved in the CFD solvers, and the various issues associated with it.

Further studies in this field may include computational analysis on a double-suction, double-volute centrifugal pump impeller [57] to investigate the effects of double-suction on its hydraulic performance and cavitation characteristics. Besides, optimization of impeller blades using surrogate modeling [58] can also lead to design of an improved centrifugal pump impeller.