The Temporal-Spatial Features of Pressure Pulsation in the Diffusers of a Large-Scale Vaned-Voluted Centrifugal Pump

: A large-scale, vaned-voluted centrifugal pump can be applied as the key component in water-transfer projects. Pressure pulsation will be an important factor in affecting the operation stability. This paper researches the propagation and spatial distribution law of blade passing frequency (BPF) and its harmonics on the design condition by numerical simulation. Experimental and numerical monitoring is conducted for pressure pulsation on four discrete points in the vaneless region, which shows that the BPF is dominant. The pulsation tracking network (PTN) is applied to research propagation law and spatial distribution law. It provides a reference for frequency domain information and visualization vaned diffuser. The amplitude of BPF and its harmonics decays rapidly in the vaneless region. BPF and BPF’s harmonics inﬂuence each other. BPF has local enhancement in the vaneless region when its harmonics attenuate. In the vaned diffuser, the pulsation amplitude of BPF attenuates rapidly, but the local high-pressure pulsation amplitude can be found on the vane blade concave side because of obstruction and accumulation of the vaned diffuser. In the volute, the pulsation amplitude of BPF is low with the decelerating attenuation. This study provides an effective method for understanding the pressure pulsation law in turbomachinery and other engineering ﬂow cases.


Introduction
Large scale centrifugal pumps suit both large flow rate and high head requirements. When applying a large-scale centrifugal pump in liquid delivering, the operation stability is very important due to the contained strong flow energy and its pulsations. In centrifugal pumps, flow is usually fully developed and turbulent, and the flow pulsation performs directly as the pressure and velocity pulsation. It will cause the alternately changing force on the pump body and induce noise, vibration, and damage [1][2][3].
In actual pump and pumping station cases, the pressure pulsation is easy to measure by using pressure pulsation sensors. Especially in stators, sensors can be stably installed by touching the local pulsing flow field. Thus, many researchers tried to study the pressure pulsations at typical positions. By measuring the pressure pulsation in the suction chamber, it was found that this pressure is mainly composed of rotation frequency and other low frequencies [4,5]. Some researchers monitored the volute and studied the influence of cut-water and volute area on pressure pulsations under different working conditions [6][7][8].
The vaned diffuser and vaneless region of pump turbine were respectively measured and studied [9][10][11]. Researchers [12,13] could also measure the pressure pulsation caused by rotating stall phenomenon by arranging sensors at the impeller outlet. Low frequency bands could be well acquired and recognized. Among them, the blade-excited frequencies are always found dominant because the impeller pumps liquid periodically. Liquid obtains energy from impeller blades and flows downstream. The turbulent flow pulses, propagates, and attenuates not only streamwise but also in the secondary directions, which makes a more complex flow. Unfortunately, researchers cannot measure the pressure pulsation everywhere in the flow field. Several discrete monitors can only show the local amplitude of the blade-excited frequencies. The propagation and attenuation details, especially the spatial distribution of temporal states, are still unknown.
Computational fluid dynamics (CFD) is a method that is used widely to analyze the flow pattern and regime with a computer. By using the correct solution method and choosing an appropriate and accurate turbulence model, accurate results matching the experimental results can be obtained. For the single-stage centrifugal pump, most of the research focuses on the rotor-stator interaction and rotation stall. There will be a high amplitude of pressure pulsation with the condition deviating from the design point [1]. Due to the interaction between the vaned diffuser and volute, blade passage frequency and its high-order harmonics were dominant downstream of the impeller [14][15][16][17]. The amplitude of blade passing frequency was affected by the position of cut-water and the distance between the impeller trailing edge and cut-water [15,16]. However, in some flow conditions, due to rotating stall and other reasons, various forms of vortex structures appeared downstream and produced a lower dominant frequency [16,17], which was also affected by volute cut-water [13]. Cavitation can also cause pressure pulsation. Researchers have focused on the relationship between cavitation and pressure pulsation. It was found that cavitation has a strong influence on the amplitude of pressure pulsation under partialload conditions [18]. When cavitation occurred, the dominant frequency changed from blade passing frequency (BPF) to a lower frequency [19]. In axial-flow pump cases, the maximum pressure pulsation may be upstream to the impeller inlet, and in vaned diffuser channels it will mainly be with the BPF. The pressure pulsation amplitude distributes non-uniformly along a spanwise direction. When the axial-flow pump operated off-design, pressure pulsation amplitude became higher [20,21]. The vortex structure formed by the leading-edge separating flow in a vaned diffuser could cause strong pressure pulsation, which is up to 3~4 times higher than that in the vaneless region [22]. For double-suction centrifugal pumps, pressure pulsation in a semi-spiral suction chamber, impeller, and volute was found dominated by impeller frequency and lower frequency bands, while volute was being dominated by BPF [10]. During the start-up process, pressure pulsation was dominated by impeller frequency, low frequency broadband, and BPF and its harmonics [23]. The impeller blade distributions of a double-suction centrifugal pump also have a significant impact on the amplitude of pressure pulsation. Research shows that an alternative distribution of the blade is helpful to reduce the amplitude of pressure pulsation [24,25]. For pump turbines, the pressure pulsation is mainly affected by a rotorstator interaction in the vaneless region. The frequency of pressure pulsation is also mainly the BPF. The BPF attenuates rapidly both upstream (to the draft tube) and downstream (to volute) under the influence of a vaned diffuser opening [26][27][28]. In other fields, CFD is also widely used to analyze pressure pulsation. Farzaneh-Gord [29] et al. researched a snubber between the turbine flow meter and the reciprocating compressors by CFD to decrease flow pulsation effects. Lai [30] et al. analyzed the Francis turbine pressure pulsation by CFD.
The method to weaken pressure pulsation is studied widely. Wu [31] et al. weakened the pressure pulsation by modifying the blade pressure side for uniform energy distribution. Lin [32] et al. used a Sinusoidal Tubercle Trailing-Edge to weaken the pressure pulsation caused by rotor-stator interaction. Posa [33] found that the wide rotor-stator clearance weakened the pressure pulsation in the vaneless region. Wang [34] et al. used alternate loading methods to weaken pressure pulsation. Zeng [35] et al. weakened the pressure pulsation by the positive lean mode. Tanasa [36][37][38][39] studied the flow-feedback technique, pulsating jet technique, adjustable diaphragm method, and axial water jet method to weaken the pressure pulsation of the vortex rope.
Generally, in these numerical predictions for the engineering, the internal flow at the pump's design-load is proven accurate with the unsteady Reynolds-averaged Navier-Stokes (URANS) method [14]. By using the URANS method, the temporal characteristics of pressure pulsation can be well captured at typical points. The BPF is already found dominant and important. However, understanding the spatial distribution of pressure pulsation and other temporal characteristics i4s still difficult. There are great challenges such as the computational cost, arrangement law of monitoring points, setup of preprocessing, resolution of visualization, and algorithms of post-processing.
Considering the disadvantages above, the pulsation tracking network based on the program that is compiled is used. The research object is a large scale vaned-voluted (with both a vaned diffuser and a volute as a diffuser) centrifugal pump that is used in the water transfer project. The accuracy of the calculation is verified by experiment. The blade-excited frequencies are deeply focused towards the downstream region of the impeller, which covers vanes, volute, and cut-water. The pulsation amplitude, phase, dominant frequencies, propagation, and attenuation are analyzed in detail. The generation, propagation, and attenuation of blade-excited frequencies are numerically investigated. It gives a completely new view of the turbulent flow inside the pump and provides new ideas for solving the operation stability problems in pump and pumping station engineering cases.

Parameters and Condition
In this research study, a large-scale vaned-voluted centrifugal pump unit is studied in model scale. The pump model consists of four parts that are inflow section, impeller, vaned diffuser, and volute. The vaned diffuser is designed for guiding water and supporting the weight of the shaft system and motor. The volute is designed for flow guidance and energy recovery. Figure 1 shows the schematic map of the pump model with the indication of the components. Coordinate X-Y-Z is defined for analysis, where X is the inflow and outflow direction, and Z is the axial direction. The design parameters and geometric parameters are listed in Table 1.  In Table 1, the specific speed n q is an important similarity parameter for comparison among similar pump units. It can be defined as: where Q d is the design flow rate (m 3 /s), H d is the design head (m), and n d is the design rotational speed (r/min). The dimensionless flow rate coefficient ϕ, the dimensionless head coefficient ψ, the efficiency of the pump η, the dimensionless coefficient of pressure C p , and the pressure pulsation dimensionless coefficient of A p are defined as: where Q is the flow rate (m 3 /s), H is the head (m), ω r is the rotation angular speed of the impeller (rad/s), R iout is the impeller outlet radius (mm), p is pressure (Pa), p ∞ is the reference pressure at the inflow inlet (Pa), ∆p is the difference between instantaneous and average pressure value (Pa), ρ is the fluid density (m 3 /kg), η h is the hydraulic efficiency, η v is the volume efficiency, and η m is the mechanical efficiency. In this case, the flow rate coefficient and head coefficient are ϕ d = 0.092 and ψ d = 0.225.

Model Test Rig, Apparatus, and Method
An experimental study was conducted on the test rig shown in Figure 2. The guideline of the test acceptance procedure is IEC-60193-1999 [40]. Physical quantities including flow rate, head, efficiency, and pressure pulsations were measured. The sampling frequency f sap of the transient pressure sensor is 2000 Hz, and the total acquisition time t acq is 12 s. Four pressure pulsation sensors were installed at the top of the vaneless region, as shown in Figure 3. They are denoted as VS + X, VS + Y, VS − X, and VS − Y based on the X-Y-Z coordinate. Due to the inconvenience of installation, the four sensors are deviated from X and Y axes. The deviation angles can be also found in Figure 3. The total uncertainty E can be decomposed into the systematic uncertainty E s and random uncertainty E r as: According to the efficiency measurement, it is with an uncertainty of less than ±0.25%. In detail, the apparatus with their types and uncertainties in the test are listed in Table 2.

Governing Equations and Turbulence Model
In this calculation, based on the URANS equation and CFD, the three-dimensional incompressible flow is numerically simulated. Because the fluid is incompressible, the continuity equation and momentum equation under the URANS method are as follows: where u i is the average velocity, u i is the fluctuating velocity, p is the pressure, t is the time, ρ is the density, µ is the dynamic viscosity, and x i is the spatial coordinate. According to these two equations, the pressure field and the velocity field can be solved out. However, the Reynolds average equations are not closed. Therefore, the turbulence model should be chosen to make the equation closed. The k-ε model [41] can simulate the complex flow, especially in the free shear layer. However, this model is hard to reflect the turbulence anisotropy characteristics. It is inaccurate to simulate the separated flow. RNG k-ε model [42] and Realizable k-ε model [43] can be better to deal with separated flow and swirling flow. There are strict requirements for the first layer of mesh at the boundary layer for these models. Wilcox k-ω model [44] improves the simulation of a low Reynolds number flow near the boundary layer, which can improve the applicability of the model and predict the flow separation phenomenon with reverse pressure gradient accurately. The Wilcox k-ω model is very sensitive for the ω value of free flow. The BSL k-ω model [45] combines Wilcox k-ω model and k-ε model by blending functions to solve the problem. However, the BSL k-ω model does not consider the transmission effect of turbulent shear stress. The SST k-ω model [45] remedies this defect by eddy viscosity coefficient correction. The SST k-ω model retains the advantages of the k-ω model in the boundary layer and k-ε model in the free shear layer. Many research results show that the SST k-ω model has higher accuracy than other RANS models [46][47][48]. The SST k-ω model is adopted as the turbulence model to close the Reynolds-averaged Navier-Stokes (RANS) equations by establishing the relationship between Reynold's stress and eddy viscosity. The kinetic energy k equation and specific dissipation rate ω equation are as follows: where P k is the production term of turbulent kinetic energy, P ω is the generating term of the specific dissipation rate. Parameters β, β * , σ k , and σ ω are constants of the turbulence model, µ is the dynamic viscosity, and F 1 is the blending function. The turbulent eddy viscosity µ t is defined as follows: where a 1 is a constant coefficient, S is an invariant measure of the strain rate, and F 2 is a second auxiliary function within 0~1. In this paper, fast Fourier transform (FFT) [49] is used for time-frequency analysis.

Pulsation Tracking Network (PTN) with Fourier Transform
Pulsation tracking can be applied to analyzing and visualizing the unknown frequencies. It clarifies the spatial distribution of temporal states in a high resolution. In this case, the pressure pulsation is mainly focused. Hence, the time-frequency conversion is necessary, and the fast Fourier transform is considered as a well-known way in engineering cases. The Fourier series (cosine signal form) of trigonometric function form of X T (t) is as follows: where c n are the amplitude component, ω f is the angular frequency, ϕ n is the initial phase angle of each sinusoidal function, C 0 is a constant term that represents a deviation from zero, and n is the integer from 1 to infinity. In this way, a complex period signal can be divided into several cosine (or sine) signals with different frequencies. Among them, each frequency component is mainly composed of three important elements that are amplitude c n , frequency ω f , and phase ϕ n , which determine the shape of the signal. Among them, c n determines the amplitude and the proportion of frequency division in the whole signal, ω f determines its fundamental nature, and ϕ n determines the state at the initial time or the specific temporal state. In the traditional time-frequency analysis, the amplitude c n and frequency nω are mainly concerned. By using this innovative PTN method, all three of the above important elements of time-frequency transformation can be specifically focused on. Well-ordered monitoring points can be arranged with a spatial high resolution. The arrangement is based on both the domain geometry and the CFD mesh node distribution. By arranging large amounts of well-ordered monitoring points, the pressure pulsation in the diffuser from the vaneless region to the volute can be detected. It can be effective to find out the source, propagation law, and attenuation law of pressure pulsation.
Spatial interpolation realizes the ability of continuous surface modeling based on discrete sampling points and estimates the variables in the unsampled area at the same time. It is a powerful method to analyze and deal with the spatial distribution and change trend. As illustrated in Figure 4, the natural neighbor interpolation is an interpolation method based on spatial auto correlation. It is widely used in the research field with the formula as follows [50]: where g(a) is the result of interpolation at the point A to be interpolated, w i (a) is the weight of the interpolation point A about the sample point P i (i = 1, . . . , n), and f i is the value at the sample point P i . The weight is determined by the following formula: where s i is the area of the polygon where the sample points participate in the interpolation, and s(a) is the area of the polygon where point A is to be interpolated. Thus, the variable size of the required point can be obtained, which is convenient to observe the spatial distribution of pressure fluctuation or other temporal states.

Setup of PTN in Diffusers
In this study, PTN is set in the vaneless region, vaned diffuser, and volute on the Z = 0 plane, which is also the mid-span of vaned diffuser as shown in Figure 5. The distribution of points is based on the vane and volute geometry and the CFD mesh distribution. Considering the influence of the vane blade, PTN is relatively dense in the vaned diffuser and coarse in the volute. The minimum space between each two PTN points are not smaller than the local CFD meshes space. The maximum space between each two PTN points depends on the local resolution. There are 12,030 monitoring points in the vaneless region and vaned diffuser, and there are 6441 monitoring points in the volute. Therefore, in total, 18,471 monitoring points are arranged as indicated. The monitoring points are distributed based on the mesh node distribution. PTN monitoring can be executed based on transient CFD simulation within a total of 10 impeller revolutions. According to the pressure pulsation data on each point, the spatial distribution and variation can be analyzed by applying interpolations.

Flow Domain Modeling
The computational domain of the pump model is shown in Figure 6. It consists of an inflow section, impeller, vaned diffuser, and volute. These domains are discretized by finite volume method with high quality hexahedral elements. The boundary layer meshing of all domains is considered. Complex geometric regions are also locally refined.

Setup of CFD
The CFD simulation was based on ANSYS CFX. The impeller was defined as a rotating domain by 1200 r/min and others were defined as stationary under the multiple reference frame. Environment pressure was set as 101,325 Pa. The fluid medium was water at 25 • C, which did not involve energy exchange. A mass flow inlet boundary was set at the inlet interface of inflow section, and a static pressure outlet boundary was set at the volute outlet interface. All the walls were set as non-slip wall boundaries. General grid interface (GGI) was used as a grid connection between the two domains. The convergence criterion was the root mean square residual (RMS) of the momentum equation and the continuity equation less than 5 × 10 −5 . The steady CFD result was set as the initial condition. Then, the following transient simulation simulated 10 impeller revolutions. The time of the impeller rotating a circle was divided to 720-time steps. And the per time step was divided to 10 iterations. The total time was 0.5 s.

Flow Domain Meshing and Checking
The grid convergence analysis was conducted using the Grid Convergence Index (GCI) method based on the Richardson extrapolation method [51] for the precision of the simulation. Three mesh schemes of fine mesh (M 1 ), medium mesh (M 2 ), and coarse mesh (M 3 ) were determined in the mesh refinement ratio, which is about 1.3. The elements number of M 1 , M 2 , and M 3 are 3473556, 1466672, and 554149. Because the research object is the characteristics of pressure pulsation in the vaneless region, the head and the efficiency of the pump and the pressure at the points VS + X Z=0 , VS + Y Z=0 , VS − X Z=0 , and VS − Y Z=0 (as shown in Figure 7) are taken as the parameters to compare as shown in Table 3. In this table, φ 1~φ3 are evaluated variables, p is the observation accuracy order, φ 21 ext and φ 32 ext are extrapolated values, e 21 ext and e 32 ext are the relative error of the extrapolated values, and GCI 21 fine and GCI 32 fine are the convergence index of the grid. For the different key parameters between M1 and M2, the relative errors of extrapolated values are between 0.06% and 6.2%. The grid convergence index is between 0.07% and 8.3%. Obviously, the M1 and M2 mesh schemes meet the requirement of convergence criteria, and the accuracy of the simulation can be well guaranteed.
The best mesh scheme has the higher prediction accuracy and better capture-ability of the flow details. However, it has a high computational cost, which should be considered. In this case, the fine mesh M 1 can be accepted to have an accurate simulation. The detailed information of mesh node and element number is shown in Table 4. The meshing result of M1 is shown in Figure 8 Table 3. Evaluation results of discretization error.

Experimental-Numerical Verification
For the accuracy of the numerical simulation results, validation is necessary to compare CFD and experimental data. The head versus the flow rate curve and the efficiency versus the flow rate curve are shown in Figure 9. Highly matched results of CFD and the experiment can be observed, especially around the design point that ϕ = ϕ d . The efficiency at ϕ d is higher than 90%, which shows a good hydraulic performance. The gap between the experimental and numerical values of efficiency at ϕ d is 0.1%, and the gap between the experimental and numerical values of the head at ϕ d is 5.2%. The comparisons of frequency-domain and time-domain of pressure fluctuation between the experiment and CFD results at ϕ d are shown in Figure 10. The frequency resolution is 4 Hz. FFT is used for the time-frequency analysis without the windowing and weighing function. Among them, the time-domain diagram is obtained by taking three impeller revolutions into consideration. Both the time-domain and frequency-domain data were found accurate especially in capturing the blade-excited frequencies. The one time, two times and three times of BPF were well simulated out. Therefore, the CFD simulation was accurate and reliable enough in the analyses below.

Pressure Contour and Velocity Vectors
To understand the pressure pulsation in the vaneless region, it is necessary to have an overview of the local flow pattern in the pump. Figure 11 shows the pressure contour and velocity vectors at ϕ d condition on the Z = 0 plane and impeller mid-span surface. Local flow regime is smooth without obvious vorticial structures. C p globally increases along the streamwise direction. Fluid gets both kinetic and pressure energy well from impeller. Then, the kinetic energy is well converted into pressure energy when flow passes through the diffusers. There are local high-pressure sites in front and local low-pressure sites on the back of the vaned diffuser leading-edge. At the leading-edge of the vaned diffuser, obvious flow separation can be found, which is the main reason why pressure suddenly varies. In detail, the flow striking causes local high pressure, and the flow separation causes local low pressure.

Pressure Pulsation at Scatter Points in Vaned Diffuser
To study the characteristics of pressure pulsation, four monitoring points were added in the vaneless region and the middle of the vaned diffuser passage along the flow (streamwise) direction on the Z = 0 plane. The distributions of added points Pa1~Pa4 are shown in Figure 3. Pa1 was located in the middle of the vaneless region. Pa2, Pa3, and Pa4 were orderly distributed in the vaned diffuser channel from leading edge to trailing edge. The time resolution was 6.944 × 10 −5 s and the frequency resolution was 4 Hz. The time-frequency characteristic diagram is shown in Figure 12.  As shown in Figure 12, the main amplitude in the vaneless region and between the vaned diffuser is a narrow band frequency of 140 Hz, which is 7 (impeller blade number) times that of the impeller frequency. Along the flow direction, the amplitude of the BPF gradually decreases, and the second and third harmonics of the BPF rapidly decrease and disappear, while other frequencies are basically unchanged. The amplitude of BPF decreases and is lower than the amplitude of impeller frequency. In the whole process, there is always a low frequency broadband. There is a low frequency band of about 60 Hz between the vaneless region and vaned diffuser. It disappears in the volute.
In the vaneless region, BPF with a high amplitude is dominant in affecting pressure pulsation. As shown in Figure 12, the amplitude of BPF attenuates rapidly along streamwise direction. This is mainly because the working condition is a design point working condition, and the flow pattern is smooth without an obviously vorticial flow or flow separation. Moreover, there is always a broadband between 0-20 Hz and the rotation frequency is 20 Hz.
In this paper, the main research object was the source, the propagation, and the attenuation law of the amplitude of BPF. Therefore, a tracking of the BPF is supplied in a perspective of innovation in the following sections.

Tracking and Visualizing of Pressure Pulsation in Diffusers 4.3.1. Main Frequency
Main frequency is an important quantity in frequency-domain analysis. It represents the dominant frequency in a certain signal. Figure 13 shows the main frequency distribution on Z = 0 plane in diffusers. In the vaneless region and vaned diffuser, the main frequency is BPF of 140 Hz and its harmonics. Along the streamwise direction in the vaned diffuser, main frequency gradually changes from BPF to the impeller frequency of 20 Hz. Based on the BPF part with green color as shown in Figure 13, there is an obvious outward diffusion phenomenon of BPF along the radial direction with the volute section area growing.

Absolute and Relative Amplitude of Typical Frequencies
The amplitude of pressure pulsation is an important parameter in this case. Based on Figures 10 and 12, there are typical frequencies that are the impeller frequency 20 Hz, BPF 140 Hz, and the second, third, and fourth harmonics of BPF 280 Hz, 420 Hz, and 560 Hz. The absolute amplitude of BPF and its harmonics is analyzed by comparing the A p contour on the Z = 0 plane in Figure 14.  The BPF and its harmonics strongly vary in diffusers as shown in Figure 14a-d. In the vaneless region, the amplitude of BPF and its harmonics gradually decreases. Among all the blade-excited frequencies, a higher frequency has a shorter propagation distance. The highest amplitude region of BPF is located in the middle of the vaneless region instead of close to the impeller blade's trailing edge. However, the highest amplitude region of the other BPF's harmonics is very close to the impeller blade's trailing-edge position. The amplitude in the vaned diffuser and the vaneless region is obviously higher than that behind the vaned diffuser.

Relative Amplitude of Typical Frequencies
Compared with the absolute amplitude, the relative amplitude will perform better in tracking the spatial propagation and distribution of a specific frequency. Thus, parameter P eqv is defined to show the relative amplitude: where p is the pressure at a certain point and a certain time, and p max and p min are the maximum and minimum pressures at this certain point. This parameter is the purity of main frequency, which eliminates the influence of its own strength. Larger P eqv is, higher than the proportion of frequency is. It is accurate in judging the pressure amplitude change of a certain frequency. Figure 15 shows the distribution of the purity P eqv on the Z = 0 plane. Figure 15a-d are the P eqv contours of BPF and its harmonics. The P eqv contour shows the blade-excited frequency series better in the vaneless region. High order harmonics of BPF dominate the region close to the impeller trailing edge. Low order harmonics of BPF dominate the middle of the vaneless region. BPF dominates the vaned diffuser channel. Outward diffusion patterns can be seen in the contour of BPF. The variation law of P eqv contour shows that the impeller blade generates complex high-order BPF harmonics. When flow deviates away from the impeller and is interrupted by the vane, high order harmonics interact on each other to enhance lower-order BPF characters. The high amplitude region in the vaneless region is wider where the volute section-area is larger along the spiral direction. It shows that the high amplitude range in the vaneless region is affected by the volute.

Detailed Distributions of Pressure Pulsation amplitudes
The A p contour and P eqv contour provide a good spatial visualization of pressure pulsation. However, more detailed information is required, especially for the local sudden change of pressure pulsation amplitudes. Thus, reference angle and positions are defined to have a better analysis on the Z = 0 plane. Considering the different radius position r, the dimensionless radius C R is defined as: Figure 16 shows the reference angle and positions. The circumferential angle θ is defined along the impeller rotation direction ω. Three streamwise lines L A , L B , and L C are set in a vaned diffuser channel that is dominated by 60 Hz frequency. Five ISO-C R circles are set within C R = 1.071~1.480. Five ISO-θ lines are set at every 5 degrees within θ = 312~332 degrees. After the θ = 332 degrees' line, three ISO-θ lines are set at every 1 degree within θ = 333~335 degrees as a refinement. Therefore, the absolute amplitude A p and relative amplitude P eqv of typical frequencies are analyzed in Figures 17-19.  The absolute amplitude A p and relative amplitude P eqv on the ISO-C R circles from θ = 0 degree to 360 degrees are shown in Figure 17. Figure 17a,b are about the BPF and its second harmonic frequency. The circumferential variation law of A p and P eqv are similar. In total, 15 (vaned diffuser blade number) periods can be seen at different C R positions. For example, for A p of BPF, the amplitude increases with the increasing volute section-area. The phenomenon becomes stronger when the C R is smaller. The circumferential fluctuation period can be divided into four parts along θ direction that are (a) 285-345 degrees, (b) 345-90 degrees, (c) 90-195 degrees, and (d) 195-285 degrees. In each part, the local peak and valley are at the same level. It shows that the volute not only affects amplitude distribution range, but it also affects the maximum value of amplitude, and the volute cut-water especially has great influence, which causes the sudden decrease of amplitude. The vaned diffuser geometry also affects the A p amplitude distribution that the peak corresponds to mid-channel when the valley is near the blade leading-edge. There is a special phenomenon at the three radial positions that are C R = 1.071, 1.173, and 1.276, which are upstream from the vaned diffuser. Their peak and valley positions strongly overlap in circumferential direction with the connection line (as indicated) perfectly towards the radial direction. It reveals the radial propagation of pressure pulsation in the vaneless region. At the radial position that C R = 1.378, the local peak is still high with the valley becoming very low. The maintaining high peaks are because of the accumulation phenomenon of pressure pulsation on the vaned diffuser blade concave side. At C R = 1.480, both the peaks and valleys strongly drop to a lower level. This is because C R = 1.480 circle is on the convex side of the blade, which is out of the accumulation region. In general, the maximum A p amplitude level decreases from about 0.008 to about 0.0015 along the radial direction. For BPF, the P eqv curves are similar to the A p curves but with an amplification that all the distribution laws are much clearer. The distribution laws A p and P eqv amplitudes of the second harmonics of BPF are similar to that of BPF. The positions of peaks and valleys still correspond to the vaned diffuser. However, it is not obviously affected by the volute. The A p and P eqv amplitude levels are stable along the circumferential direction. The amplitude of second harmonics of BPF has a rapid drop from C R = 1.071 to 1.173 as indicated. It drops by at least 0.002, which is approximately 0.3% of the pump head.  Because some specific pressure pulsation propagations are found along the radial direction, it is necessary to have detailed analyses along the radial reference lines. As the pressure wave can pass through objects, the radial distributions, especially across blades, can indicate special variations. To understand the detailed radial variation laws, the absolute amplitude A p on the ISO-θ lines from C R = 1.03 to 1.68 are shown in Figure 18.
To have more direct results of the attenuation of typical frequencies, the derivative of A p against C R is also plotted as dA p /dC R . Figure 18a,b are the radial distribution of A p and dA p /dC R of BPF and its second harmonic frequency. Their variation laws are very complex. For a better comparison, six different regions are divided as shown in Figure 18a,b. They are, respectively, named 1 impeller blade trailing-edge downstream region, 2 amplitude fluctuation region, 3 vaneless attenuation region, 4 pre-guide-vane fluctuating region, 5 inter-channel attenuation region, and 6 post-guide-vane attenuation region. The 1 impeller blade trailing-edge region is near the impeller blade trailing-edge. It is a complex region where both the A p of BPF and the second harmonic of BPF are strong, about 0.0055~0.0058. Different θ positions have different variation tendencies. For example, the θ = 327 degrees' site is near one of the peaks in Figure 17a,b. The A p of BPF and the second harmonic of BPF may rise locally. On the contrary, the θ = 335 degrees' site is near one of the valleys in Figure 17a,b. The A p of BPF and the second harmonic of BPF may decrease locally. In the figures dA p /dC R , both positive and negative values can be found in region 1 . The 2 amplitude fluctuation region, 3 vaneless attenuation region, and 4 pre-guide-vane fluctuating region are in the middle of the vaneless region. In the 2 amplitude fluctuation region, BPF suddenly changes from decreasing to increasing, but the second harmonic of BPF still keeps rA p idly decreasing. The dA p /dC R values of BPF are positive, and that of the second harmonic of BPF is negative. It may indicate the interaction and transmission of pressure pulsation amplitude from the second harmonic of BPF to BPF. Along the radial direction, the BPF's A p increases slower with the quicker attenuation of 2BPF's A p when the 2BPF's A p becomes very low. This is also why the maximum A p of BPF is in the middle of the vaneless region instead of near the impeller trailing edge. In the 3 vaneless attenuation region, the A p of both BPF and the second harmonic of BPF decreases. The values of dA p /dC R are negative for BPF and the second harmonic of BPF, which indicates an accelerated attenuation in this region. However, for the second harmonic of BPF, the attenuation becomes slower than that in the 2 region. The 4 pre-guide-vane fluctuating region is slightly similar to the 3 region because A p is also locally decreasing. However, it is different because dA p /dC R is fluctuating in the 4 region. It indicates a very complex local flow pulsation due to the influence of the vaned diffuser. The 5 inter-channel attenuation region is in the vaned diffuser channel near the blade concave surface. There are two different typical conditions. First, near one of the peaks in Figure 17a,b, the A p attenuates stably and flatly with a low level in the 5 region because of the prevention effect of the vaned diffuser blade. Second, near the leading edge of the vaned diffuser, the A p is rapidly rising and then suddenly decreasing down with a local peak. There is a spatial shifting phenomenon of local peaks that can be indicated on the C R -A p figure. With the decreasing of θ, for example, from 335 to 332 degrees, the propagation direction shifts from one blade leading edge to the upstream blade's concave surface. It indicates the deflection of the pressure wave under the influence of the vaned diffuser leading edge. The 6 post-guide-vane attenuation region is downstream to the vaned diffuser from the blade convex side to the volute. The A p in this region attenuates slower and slower. The A p of the second harmonic of BPF decreases to a very low level. The value of dA p /dC R is at a small negative level. It indicates the slow attenuation of pressure pulsation amplitude without the vaned diffuser's interference.
The variation of A p on the streamwise lines is also important for understanding the pressure pulsation propagation in the vaned diffuser channels. The streamwise length is normalized into L st = 0~1, as shown in Figure 19. As illustrated in Figure 16, L A is the streamwise line on the blade concave side, and L C is the streamwise line on the blade convex side. The A p on the two lines are similar as there is a local increasing and decreasing around the blade leading-edge. It shows the importance of the leading-edge geometry again in affecting pressure pulsation. As illustrated in Figure 16, L B is in the mid-channel, which is relatively far away from the leading edge. Therefore, A p continually decreases on L B from L st = 0 to L st = 1.
The vaned diffuser blade, especially the leading edge, has a strong effect on pressure pulsation. BPF and its harmonics are found interacting with each other on local strong amplitude variations. The attenuation and propagation laws of pressure pulsation amplitude are somehow indicated, which provides a new understanding of pulsation flow driven by rotor-stator interaction in the vaneless region, vane channels, and downstream diffusers.

Phase and Phase Difference
The initial phase is important to analyze phase and phase difference. In this study, it is defined at the volute cut-water. Figure 20 shows the distributions of phase difference of typical frequencies on the Z = 0 plane. As the sine and cosine functions are the basis functions in Fourier transform, the phase will change within -π to π. If the phase difference is π, which is called antiphase, all the temporal states, including displacement, velocity, and acceleration, are completely on the contrary. If the phase difference is 0 or 2π, all the temporal states are in the same status. Phase π and −π in Figure 20 is also the same. Due to the interpolation problem, mutations can be seen between adjacent π and −π, which are smooth transitions. Figure 20a-d shows the phase difference of BPF and its harmonics. In the vaneless region, the pattern of phase difference is strongly periodic. The number of periods relates to the impeller blade number. For example, 140 Hz frequency (BPF), 280 Hz frequency, 420 Hz frequency, and 560 Hz frequency have respectively 7, 14, 2,1 and 28 periods. Without the interference of vane blades, the pressure wave fluctuates along the circumferential direction repeating the same states. This means that the BPF and its harmonics in the vaneless region are only influenced by the impeller, and the impeller influences the BPF and its harmonics mainly in the circumferential direction. The BPF and its harmonics propagate along the radial direction. The vaned diffuser, especially the leading edge, interferes with the propagation of pressure wave and broke the circumferential phase fluctuation pattern. In a single vaned diffuser channel, the phase is almost the same. As pressure wave propagates along the radial direction, the circumferential phase fluctuation pattern is outward diffusing from the vaned diffuser to the volute with diffraction patterns. The lower the order the harmonic of BPF is, the further it diffuses. This means in the outlet of the vaned diffusers, as indicated in Figure 20a, the outward diffusion of BPF is obvious. On the contrary, as shown in Figure 20d, the circumferential phase fluctuation pattern of four times of BPF stops on the vaned diffuser leading edge. The volute cut-water is important in affecting the pressure wave propagation. As indicated in Figure 20a, the outward diffusion of BPF is fluctuating between −0.5π to −π. The two times and four times of BPF have also a phase change from 0 to −π there. In the volute straight section, the phase difference pattern is divided into three parts, as indicated by upper, middle, and lower. The upper and lower part have the same phase. The middle part has the phase with a difference of π with the upper and lower parts. This three-part pattern is related to the local merging flow in the volute straight section, as shown in Figure 11. Phase differences can be found among different frequencies. The spatial distribution of phase well indicates the trend and state of BPF and its harmonics' pressure variations. As shown in Figure 20a, the phase difference between adjacent vane channels is about π, which means that the pressure change between adjacent vanes channels is completely opposite. Similarly, at the outlet of the vaned diffusers, the diffraction fringes indicate that there is variable pressure change. However, it has little effect because the amplitude of BPF at the outlet is not strong. It is difficult to observe these phenomena by conventional methods.

Discussions
According to the PTN-based analyses above, there are some discussions for a better summary of the tracked pulsation characters. Figure 21 is drawn for detailed illustrations. The pulsation amplitude is found alternatively weak and strong in the vaneless region. As the blade-excited pressure wave propagates along the radial direction, this alternative amplitude pattern indicates the interference of the vaned diffuser on pressure pulsation amplitude. Based on the PTN analysis, it is always weaker for the amplitude of BPF near the vaned diffuser leading edge in the vaneless region due to the vaned diffuser blocking effect. It is always stronger for the amplitude of BPF near the vaned diffuser mid-channel.
Vaned diffuser geometry has two main influences that are deviation and accumulation on pressure pulsation propagation. As the vaned diffuser is designed for flow guidance, the radial flow from the impeller is guided into a more circumferential direction for downstream volute. Leading edge, which is the site nearest to the impeller, may deviate flow onto the concave side of the adjacent vane blade. This why a deviation can be found on the connection of pulsation amplitude peaks in the vaned diffuser channel. At the same time, the blade concave side prevents the deviated flow. The accumulation effect can be found on the vane blade concave side with local strong pressure pulsation amplitude. Similar effects can be found at the volute cut-water, which can be treated as a single vane. It also strongly affects the pressure pulsation. Therefore, the impeller is the main source of pulsating signal, and the vaned diffuser enhances pulsation amplitude in the local region, especially near the head of blade. Then, the pulsation amplitude attenuates quickly in the vaned diffuser channel. Overall, pressure pulsation amplitude generally attenuates in the vaned diffuser. In the vaneless region, the BPF series frequencies rapidly attenuate but with BPF local enhancements due to the influence on each other. High order harmonics of BPF attenuate earlier. Low order harmonics such as BPF and two times BPF may shortly self-enhance by receiving the energy from attenuated high-order harmonics. As the vaned diffuser geometry has the obstruction effect on the pressure wave, the pressure pulsation amplitude attenuates slower, which even increases near the vane blade, especially at the leading edge. However, there is a small region with low amplitude at the leading edge of blade. There is a high pressure and a low flow velocity region in the front of the blade leading edge because of the impact between flow and the leading edge. In this region, the pressure pulsation is stable and with low amplitude. The closer to vane blade concave side, the slower the attenuation is, but at the convex side, the amplitude attenuates rapidly. The amplitude attenuates to a lower level. After the pressure wave bypasses the vaned diffuser, the amplitude attenuation decelerates without strong preventions. Hence, the vaned diffuser made the amplitude of the pulsating signal attenuate overall.
The spatial distribution of phase difference is interesting. Phase change from π to 0 to −π indicates a complete waveform and sub-flow of a certain frequency. In the vaneless region, BPF and its harmonics have several complete waveform periods along the circumferential direction. It indicates the circumferential sub-flow in the vaneless region. However, there are some local incomplete phase changes. For example, the phase of BPF fluctuates between −0.5π to −π near the volute cut-water. The composition of the pattern is two or more sub-flow that have similar frequencies, but different directions, interfere with each other and mix, and then they result in a phenomenon such as wave interference.
Vaneless region, guide vaned diffuser, impeller blades and vaned diffuser blades have obvious influences on the propagation characteristics of BPF and its harmonics based on the above analysis. The vaneless region can be enlarged for the amplitude rapid attenuation when designing the pump, and close attention must be paid to the head and the blade angle at the medium of the vaned diffuser blade to avoid the enhancement of pressure pulsation amplitude.

Conclusions
In this study, the propagation and attenuation of pressure pulsation in the diffusers of large-scale vaned-voluted centrifugal pumps were numerically tracked and visualized. The experiment verifies the accuracy of the simulation. Conclusions can be drawn as follows: (1) When a large-scale vaned-voluted centrifugal pump is operating at design-load, impeller frequency and blade-excited frequency are the pressure pulsations dominant in diffusers. The impeller frequency of 20 Hz, the blade passing frequency (BPF) of 140 Hz, and its high order harmonics that are 280 Hz, 420 Hz, and 560 Hz can be obviously found. By analyzing the variation law of pressure pulsation in one vaned diffuser channel, attenuations of the BPF and BPF's harmonics are found along the streamwise direction. (2) By using the pulsation tracking network (PTN), it is possible to have a more detailed visualization of the source, propagation law, and attenuation law of pulsation frequencies. The BPF series are obviously related to the impeller and vaned diffuser. The BPF series is concentrated in the vaneless region and the inlet of the vaned diffuser channel. The significant influence region and amplitude of BPF have a positive correlation with the sectional area of the volute. (3) PTN results reveal the propagation and attenuation law of pressure pulsation frequencies in diffusers of the vaned-voluted pump. In the vaneless region, the pressure wave propagates along the radial direction. Vaned diffuser blades (and volute cut-water) have influence on the pressure pulsation. The accumulation effect and the flowprevention effect of the vaned diffuser will cause local enhancement of amplitude. Hence, the pulsation amplitude attenuation has mainly three stages. In the vaneless region, it is rapid. Near the vaned diffuser, it attenuates first slowly then fast. After the vaned diffuser, the attenuation decelerates. (4) The phase distribution can indicate sub-flow with certain frequencies. BPF and BPF's harmonics fluctuate from π to 0 to −π in the vaneless region along the circumferential direction. It shows the circumferential sub-flow driven by the impeller rotation. The outward diffusing pattern of BPF in the volute spiral section shows the sub-flow along the vaned diffuser in the streamwise direction. The fluctuating phase pattern between −0.5π and −π indicates local interference by several sub-flows, with some similar specific pulsation frequencies.
This study only analyzes the BPF and its harmonics at design condition. Considering the simplicity of the operation, the FFT method is selected. It has certain limitations. More reasonable distribution of the PTN points need to be studied, and better frequency domain analysis methods can be used so that more working conditions can be analyzed. The selection of the interpolation method also needs to be improved. This study can be combined with some typical physical phenomena, and further research will be developed in the future.