Investigation on the unsteady pressure pulsations and related vortical structures in a molten salt pump

The high‐temperature molten salt pump (MSP) is the core equipment of the thermal storage system for concentrating solar power plants, which circulates the molten salt of the thermal storage medium. This study employs numerical simulation of a two‐stage MSP with a storage tank to find the correlation between pressure pulsations and the evolutions of vortical structures based on renormalization group k‐ε turbulence model. It is shown that rotor/stator interaction is the most significant excitation source of pressure pulsations with the prominent excitations being the diffuser passing frequency (fDPF) in the impeller and the blade passing frequency (fBPF). Furthermore, Q‐criterion is employed to capture the periodic shedding vortexes with a frequency of about 3fR at the suction surface and the trailing edge in diffusers. Impinging vortex and shedding vortex with the frequency 6fR on the leading edge of the diffuser is captured. The vortex shedding characteristics are significantly different under various conditions, and the separation vortexes appear at the leading edge of pressure surface in the diffusers at higher flow rates. Finally, it is obvious that the excited components and their magnitudes in the pressure spectrum are closely associated with the unsteady vortical structures within the model pump.

the bottom of the tank. Consequently, the length of the pump shaft is being more than 10 meters. Compared with conventional centrifugal pump with volute, the flow structures of the MSP are different and complex due to the suction type and the multistage. In Sandia National Laboratories (SNL) National Solar Thermal Test Facility (NSTTF), Barth and Smith et al. 6,7 operated the test loops to demonstrate the performance, bearing reliability, service life of MSPs and valves, respectively, for the commercial CSP plants in early 1990s. The combined effects of hydromechanics, structural mechanics, rotor dynamics and thermodynamics should be considered during the design, and it leads to a systematic complicated engineering for the design of the MSP. Especially in terms of the vibration properties of the structure itself, which are mainly caused by the internal unsteady flow structures.
Pressure pulsation is induced by the unsteady flow inside the centrifugal pump, and it may cause intense structural vibration by alternating stresses on the inner surface of the parts at extreme and abnormal condition. 8 The unsteady interaction between the rotating impeller and stationary diffuser is usually defined as the rotor/stator interaction (RSI), which is the major excitation source of pressure pulsations and attracts the researchers' attention. [9][10][11] To reduce the intensity of RSI, several optimization measures, including adding to splitters in the channel of impeller, modification of trailing edge shape, changing the gap between rotor and stator and so on, have been put forward validated by experiments. [12][13][14][15] As for pressure pulsations, many research are conducted to establish its relation with the complex internal flow fields. Liu et al. 16 proposed an optimize performance of multistage pump by theoretical prediction based on Oseen vortex. Tan 17 investigated the tip leakage vortex in a mixed flow pump by numerical simulation. An Yu et al. 18 investigated the unsteady natural cavitation and ventilated cavitation with the large eddy simulation method. Through the investigations, it is found that the shed vortex also develops at blade trailing of pumps, which is confirmed by a lot of scholars. Gonzalez and Santolaria 19 and Zhang et al. 20 used the numerical method to simulate the flow fields of the centrifugal pumps, and they focused on the trailing vortex structures at near the impeller outlet and the tongue. results discussed the relationship between vortex movement and the pressure pulsation. With particle image velocimetry (PIV) flow visualization measurement techniques, Keller 21 successfully captured the collision of wake vortex shedding, motion, development and dissipation at high flow rate. Later the more comprehensive research at various flow rates of the pump were carried out by Zhang et al. 22 The above phenomena may be more complicated for the pumps with diffuser due to the RSI interaction. Wang and Tsukamoto 23 founded two-dimensional (2D) vortexes method to investigate pressure pulsations both in the impeller and diffuser regions. By through numerical results, Shibata et al. 24 considered that the vortex at the trailing edge of the diffuser vanes blocked the flow and induced a separate flow structure at the leading edge, which finally resulted in the rotating stall under low flow rate. Thus, vortex-induced unstable factor would have a significant impact on the pressure pulsation energy of the model pump, which should be concerned and analyzed during the pump design.
At present, relevant investigation between flow structures and pressure pulsation in MSP is rare. Besides, some investigations do not take the effect of the storage tank into consideration, which is often replaced by the straight pipe. In fact, the existing of the salt storage tank would affect the inlet flow structures and also the pressure pulsations of the MSP. Thus, it is necessary to consider the effect caused by the tank. To ensure the safe and stable operation of the MSP, investigations on the unsteady flow fields and pressure pulsations are fundamental and essential.
Based on the Reynolds-Averaged Navier Stokes (RANS) approach for numerical calculation of MSP under different operation flow rates, pressure pulsation and vortex evolution will be discussed in this paper. The investigation on this study can provide a comprehensive flow fields and the induced pressure pulsations of the MSP, which will be helpful for the design and safe operation of the pump. Nomenclatures of professional terms is presented in Table 1.

| GEOMETRY AND MAIN PARAMETERS
The model pump investigated in the paper is used in a currently building 100 MW CSP station. The total length of the MSP is 19 m, and the submerged depth is up to 16 m. The MSP is closed to tank bottom, possibly generating unstable vortex strip attached to the bottom and sidewall. 25 Thus a molten salt tank with diameter of 16 m in reality is introduced to consider the real flow conditions for the pump. Since, the disturbance flow at the pump inlet does not affect the whole region of the tank, the reasonable size of the molten salt tank is reduced to 6 m after optimization to decrease the calculation amount. 26 The length of outlet pipe is shortened to 3 m, and the depth of the liquid surface is setting as the normal operating condition. The MSP computational domain contains impellers, diffusers, suction chamber, and molten salt tank, and its assembly model is shown in Figure 1. The impellers and diffusers are immersed in the fluid, and the gap between fluid domain and hydraulic components is the wall thickness of component. The axial projection drawing of impeller and diffuser is shown in Figure 2. The main design parameters are presented in Table 2.

| Mesh generation
Mesh quality largely determines the numerical simulation accuracy. The computational mesh for geometry domain is created by ANSYS-ICEM, which is shown in Figure 3. Wherein the meshes of impellers, diffusers, salt storage tank, bottom suction chamber and outlet pipe are generated with structural multiblock, respectively. The boundary layers are encrypted on impellers and diffusers blade wall, and local y+ value is less than 40, which is satisfied with computational requirements of renormalization group (RNG) k-ε turbulence model. Considering the mesh sensitive analysis, the total number of meshes is determined to be approximately 12 million, with 2.2 million mesh for a single impeller and 1.8 million mesh for a single diffuser. The mesh sensitive analysis is shown in Figure 4.

| Unsteady calculation setup
The continuity equation and the 3D incompressible unsteady Reynolds-averaged Navier-Stokes equations are solved by the CFD code Fluent 19.0. The viscosity of the liquid ternary nitrate delivered by MSP at working temperature 570°C is similar to clean water at room temperature, and liquid ternary nitrate has a stable physical properties at working temperature. There is evidence that the change of viscosity and density have little effect on the hydraulics characteristics and flow pattern of MSP. 27 Considering the convenience of calculation and convenience of comparison with prototype water experiments in the future, clean water is selected as medium. The turbulence flow is simulated by the RNG k-ε turbulence model with standard wall function near wall. 28 All the physical walls of the MSP are set as nonslip walls. The SIMPLEC algorithm disposes the coupling between velocity and pressure, and the second-order upwind approach is used to discretize the velocity, turbulent kinetic energy, and turbulent energy dissipation rate. A velocity-inlet is specified at the inlet pipe which is given by the flow rate and outlet boundary is set as outflow. Turbulent intensity is defined as the ratio of the root mean square of speed fluctuations to the average speed which below 1% is called low turbulence intensity, and above 10% is called high turbulence. The turbulent intensity is set to 5%, a median value, by theoretical calculation. For small undulating free surface and lager volume fluid motions, the rigid-lid hypothesis method are often used to avoid the problem of free surface which assumes that the free surface is a movable rigid cap, and the boundary condition meet the non-penetration condition. The storage tank liquid surface decline 0.003% of total height per second which satisfies the premise of rigid-lid hypothesis method when MSPs are normally operated, so the rigid-lid hypothesis method is used for the tank fluid free liquid surface to simplify the calculation. 29 The steady results are set as the initial boundary for the transient simulation. The time step Δt is set to 1.1495 × 10 −4 s, corresponding to a rotating angle change of Δϕ = 1°. per time step in the case of the rotating speed of 1450 rpm. The maximum number of iterations for each time step is set to 40 to guarantee all maximum residuals below 2 × 10 −6 . Besides, to obtain sufficient and accurate numerical results, nearly 17 impeller revolutions are calculated. 30 Local pressure changes have significant effects on internal flow of the MSP, hence several monitoring points were selected to analyze the unsteady flow structures and pressure pulsations. As a result of rotational symmetry structure of impellers and diffusers, only a single flow channel is selected to analysis. It's noted that all points are arranged on middle stream surface in impellers and diffusers as shown in Figure 5. In Figure 5A, IP1, IP2, and IP3 are, respectively, located on the leading edge, middle position and trailing edge on the pressure surface (IPS) of the Impeller Ⅰ. IS1, IS2, and IS3 are placed on its suction surface (ISS), and IM1, IM2, and IM3 are set in the middle region of the channel. Other monitoring points are defined in impeller Ⅱ and diffusers as shown in Figure 5B-D.

| Pressure pulsation
The static pressure coefficient C p (as seen in Table 1) is introduced to treat the pressure value of internal monitoring points, so that the calculated results can be used to discuss the pressure pulsation characteristics of the pump with the actual working medium (molten salt). As shown in Figure 7, the time domain signals of the monitoring points on pressure surface of impellers and diffusers during the 0-0.2 s, that is about five impeller rotation periods under the design working condition. Due to the periodic rotor-stator interaction between the impeller and the diffuser, the time-domain signal of pressure pulsation presents the periodic characteristic. There are significant differences among different monitoring points. Pressure pulsation signals in the first-stage are displayed regularly, especially at IP3 with eight waves and DP1 with six waves in single impeller rotating period, respectively. Compared with the first stage, the time domain signals in the secondary impeller seem to be F I G U R E 6 Performance curves of the model pump. In addition, there are large-scale waves of five periods during 0-0.2 s at monitoring points in impeller owing to the rotational points synchronously with the impellers. Fast Fourier transform (FFT) is used to convert the time-domain signal into frequency domain signal, and the spectral resolution is about 1.7 Hz. The fundamental frequency of pressure pulsation is determined by the impeller rotate speed n, number of impeller blade Z i and diffuser blade Z d . The fundamental frequency is f BPF = nZ i /60 in the diffuser, and f DPF = nZ d /60 in the impeller. 31 Figure 8 shows the pressure spectra of the two-stage model pump under 0.8, 1.0, and 1.2Q d flow rates, respectively. Due to the rotation of the impeller, the monitoring points on pressure surface are characterized by the shaft frequency f R = n/60 and its high harmonics, especially at 4f R , the corresponding amplitude is  prominent. When the position is closer to the trailing edge, the amplitude at f DPF = 8f R caused by RSI is greater. Compared with 1.0Q d operating condition, the flow structures at 0.8Q d are more complex leading to some low and more prominent frequency signals in spectrum. However, the pressure pulsation amplitude of typical frequency Xf DPF (X = 1,2, … N) at 1.2Q d increases significantly, and it is indicated that the RSI effect is stronger due to in the intense flow at the impeller outlet. Compared with the first-stage pressure spectra, a large number of low-frequency signals occur in the secondary impeller. At high flow rate, the phenomenon is alleviated.
Affected by the unsteady flow discharged from the impeller interacting with the diffuser, flow filed in the diffuser shows pulsation characteristics. So it is necessary to analyze pressure pulsations in the diffuser. Figure 8 shows pressure spectra of the points in the diffuser. As can be seen from the figure, the dominant frequencies in secondary diffuser are similar to the first one, and there are more low-frequency signals in diffuser Ⅱ. The pressure spectra show significant difference under different operating conditions, especially the prominent low-frequency components of 0-3f R under 0.8Q d condition. In fact, the leading edge of diffusers is the most intense RSI region, therefore the amplitude at fundamental frequency f BPF and 2f BPF gets much higher than the other locations in the diffusers. The effect of RSI near the tailing edge of diffuser is weakened, which leads to the amplitude at f BPF to be reduced. The pressure spectrum is dominated by the low-frequency signals.
To evaluate the whole components in the pressure spectra when the pump operates under different working conditions, the root-mean-square value of the pressure pulsation coefficient is introduced, as shown in the following formula: where A P1 is the pressure amplitude at the beginning of the frequency band, A PN is the amplitude at the end of the band, and N is the number of components occurring in the frequency band. The RMS value represents entire energy of pressure pulsation of the single point. As shown in Figures 8 and 9, the components in pressure pulsation are concentrated within 20f R , so the pressure pulsation amplitude in the 0-20f R frequency band is treated by the RMS method. The RMS values for the different points in the pump are shown in Figure 10. On the whole, the RMS value gradually decreases with the increase of flow rate, and the RMS difference among different monitoring points in the diffusers is greater than that in the impellers. As can be seen from Figure 10A, the change trends of RMS values of the impellers are consistent along the blade surface at different flow rates. the RMS value of secondary impeller is slightly greater than the first impeller, especially on the blade leading edge, the RMS value is much higher than the other positions (the monitoring points is exactly placed in high vorticity region as shown in Figure 13A).
The pressure pulsation in the diffusers at double sides is complicated with flow rates and monitoring points location as shown in Figure 10B. The RMS value at 0.8Q d is much higher than the other operating conditions, which further demonstrates that the flow structures at the low flow rates are more complex resulting in intense pressure pulsation. The RMS value changes slightly on the pressure surface, only enlarges at leading edge due to intense RSI effect directly. But there are large differences on suction surface under different operating conditions. Compared to the other positions on suction surface, the RMS value at the center of the blade is the largest under 0.8Q d operating condition, but increases significantly at trailing edge at the 1.0 and 1.2Q d operating conditions. As shown in Figure 10C,D, the RMS values on DM monitoring points are the lowest at 1.2Q d flow rate, and totally larger in the second stage than the first, which basically the same as blade surface of impellers and diffusers.
From Figure 10, it is evident that RMS values of the second impeller and diffuser are larger than the first stage. It is considered that there are complex flow structures originated from impeller Ⅰ and diffuser Ⅰ, directly affecting the flow in the secondary. Then local pressure distribution in the secondary pulsates dramatically under RSI effect. It means that the turbulent flow from the upstream stage will increase the pulsation characteristics in the second stage. In addition, pressure pulsation characteristics at the center of diffusers channel are consistent with the trend on the suction surface, which indicates that the flow condition of middle passage is closely related to the suction surface. Thus, it's inferred that there are unidentified flow structures affecting pressure pulsations except the RSI effect, especially in diffusers.

| Flow structure in suction chamber
Research show that there are multiple kinds of vortices due to the close distance between the pump inlet and the bottom of tank, which can be divided into free surfaces vortices, submerged vortices et al according to originating locations. 32 Uniform suction condition would be affected, which will deteriorate the pump efficiency, stability, cavitation, even intense noise and vibration due to the suction vortex. 33 In the MSP, suction chamber was placed 0.35 m above tank bottom, 0.5 m offset sidewall. Thus, it needs to verify the installation rationalization of pump by numerical simulation method to prevent the large-scale unstable vortices or vortex strip generating at the pump inlet. The pressure distributions on different sections in the suction chamber are shown in Figure 11A.
The pressure value gets lower from inlet to impeller Ⅰ. It seems that streamlines are almost smooth at the design flow rate, and the velocity is increased by impeller rotating effect. In Figure 11B, a long and narrow vortex strip occurs near sidewall of suction chamber with a very small z-vorticity legend，which has hardly any influence on the inflow conditions. Hence, it is considered that the install position on tank is reasonable, and obvious unstable flow structures were not generated.

| Pressure distributions in impeller and diffuser
The fluid will gain energy through impellers and change flow direction within the diffusers, which converts a part of kinetic energy into static pressure to reduce the flow losses. It is believed that the flow state in the impellers and diffusers basically determines the performance and unsteady pressure pulsation characteristics of the MSP. Figure 12 presents the static pressure contours and flow streamlines within the impeller and diffuser under the design operating point. It is shown from streamlines that when the fluid particle flows through the diffuser Ⅰ, its circular velocity is not completely eliminated. According to the streamlines inside impeller Ⅱ, it can be seen that the inflow conditions are relatively unsmooth and irregular. Both the above reasons account for the reduction of the real head of impeller Ⅱ. Besides the added value of the secondary static pressure is lower than that of the first stage in pressure contours. Thus, it indicates that the secondary impeller has a weak power capacity. The pressure value increases along the moving direction due to impellers working, and there is a slight increase from inlet to outlet in diffuser because of its section enlargement. In addition, it is quite clear that pressure values on pressure surface are higher than the suction surface both in the impeller and diffuser. Moreover, from streamlines, we find that the flow within the impeller is smooth, but a large scale back flow structure appears on the outlet of the diffuser because of excessive curvature in the latter half of diffusers blade.

| Vortex structures in the model pump
Vorticity is defined as the curl of fluid velocity vector, which is generally used to measure the strength and direction of vortices. 34 It can be seen from Figure 13A that, under the design condition, the strong vorticity region in the impeller concentrates on the leading edge of blade and near suction surface, and there is wake structures attached to the tail section. The large-scale shedding vortexes are generated near the trailing edge of the diffuser suction surface. It even moves to the secondary impeller with complete dissipation, which directly leads to the deterioration of the secondary inlet flow condition. High vorticity area the suction side of the diffusers vane is defined as α, and the trailing edge vortex is defined as β. So the typical vortex regions of the diffuser could be divided into α1, β1 and α2, and β2 for the two stage model pump.
Deviation from design conditions may induce more complex unsteady flow structures. Figure 13B,C describes the vorticity contours at 0.8 and 1.2Q d , and it is found that there is a significant difference between the off-design condition and the design condition. At 0.8Q d operating condition, the high vorticity zone near the suction surface of the impeller blade is observed, which is due to the change of flow rate and the change of the flow attack angle entering the blade. finally, it would form the detachment region on the suction surface of the blade. A vortex appears at the SS leading edge of diffuser, which is defined as vortex γ. It is inferred that the vortex γ is generated by the wake vortex discharged from upstream impeller. Compared with the design condition, the location of α structure moves forward to the middle of diffuser channel at 0.8Q d , and then to the trailing edge of the diffuser at 1.2Q d . At the same time, the vortex χ is captured at 1.2Q d , which divide into multiple branches to downstream.

| Vortex evolution and induced pressure pulsation
RSI effect is the major excitation source of pressure pulsations in a diffuser pump. 10 So it is necessary to analyze the vortex evolution. The Q criterion is based on the invariants of the velocity gradient tensor, which is widely used to study the evolution of the complex vortical structures. By comparison with vorticity criterion above, it is found that Q criterion can simultaneously capture vortical structures in the impeller and the shedding vortexes on the suction surface of diffuser. It can also identify positive and negative vortexes, which can more accurately show the periodic shedding of Karmen vortexes. To explore the transient evolution process of the vortex, Q criterion is adopted to resolve the transient flow fields. Figure 14 presents the transient evolution of vortical structures by Q criterion, and the periodic vortex shedding is captured in diffuser Ⅰ. Here, t = t + 0°is defined as original moment, and t = t + x°is the moment of impeller rotates x°. As shown in these pictures, positive vortex initially leaves in the α1 region at t = t + 0°moment, then it enlarges and moves toward the downstream of the diffuser. meanwhile, the negative vortex structure occurs at t = t + 20°moment. When at t = t + 120°moment, the first positive vortex begins to break and dissipate, and the second leaves the blade surface exactly. so, positive and negative vortical structures are alternately generated, which expand and shed from the diffuser surface, and finally dissipate when moving to the downstream. The vortex shedding process is repeated over time with a period of about 1/3T, and the corresponding vortex shedding frequency is 3f R . The vortex in the area β1 is generated at the trailing edge, and the evolution process is more complex probably due to the evolution of the upstream α1 vortexes, resulting in an uncertain vortex shedding period of 1/3-1/2T. Vortical structures in the diffuser Ⅱ is shown in Figure 15, and it is found that the v evolutions of α2, and β2 are similar to the diffuser Ⅰ, but their shapes are stretched and deformed apparently during the movements. The vortex scale is increased, which affects the unsteady pressure pulsation characteristics. Besides it is found that the vortex shedding periods of α2, and β2 equal to 1/3T approximately from the unsteady vortical contours, and it is identical with diffuser Ⅰ.
The pressure spectra at diffuser trailing edge regions are shown in Figure 16. It is observed that the amplitudes of low-frequency components at the trailing edge region are much larger than that at f BPF , which are caused by the weakened influence of RSI and the generated vortices in this region. Additionally, the spectra show that the pressure pulsation intensity in diffuser shows a descending order in different measuring region, namely suction surface, central region, and pressure surface, which is related to the vortical distribution on the suction surface. when the upstream shedding vortexes αstrike with trailing edge vortexes βat the diffuser outlet, some induced components concentrate at 2-3f R . Due to vortexes αmove downstream with flow, the 3f R component at the DS3 is obvious. from comparison, the characteristics of low-frequency components of the secondary diffuser are similar to those of the first stage.
To comprehensively and accurately analyze the internal unsteady flow structures in the model pump, pressure pulsation intensity C p * is introduced and calculated by the nondimensional pressure coefficient where N is the sampling number during one period, that is, the number of time steps in one period. The comparison of pressure pulsation intensity between two diffusers is revealed in Figure 17. It is noted that both the size of the spot and its color scale represents pressure pulsation intensity, that is the larger the scatter, the greater the intensity. It can be observed  that areas of high intensity are concentrated on SS trailing edge and PS leading edge, which is consistent with the typical vortex shedding regions. Therefore, it is believed that the pressure pulsation is induced by the vortices shedding. Besides, the pressure pulsation intensity in the diffuser Ⅱ is larger than the diffuser Ⅰ slightly due to more complex flow conditions. According to vortexes evolution process, it can be demonstrated that flow structures are basically consistent with each other in the diffuser. Thus, the study of flow structures under off-design operating conditions is just concentrated in the first diffuser. Figure 18 presents the transient evolution of vortical structures at 0.8Q d , and obvious periodic vortex shedding is captured in diffuser Ⅰ. From t = t + 0°to t = t + 120°moments, vortex α exactly complete a vortex shedding period, and the vortex pattern is not the same with the design working condition. Thus, the vortex shedding period of structure α is also 1/3T corresponding to 3f R frequency of pressure pulsation at diffuser trailing edge as shown in Figure 20A. Meanwhile, the evolution of vortex street β attached to trailing edge is almost irregular, and the period maybe equal to 1/3-1/2T. It is noticed that the vortex structure γ is striking the leading edge of diffuser, which is shed from impeller. It can be seen in detail that the γ structure impacts the leading edge of diffuser at t = T + 0°, and the secondary approaches at T + 40°, and attacks the forepart again at t = t + 60°, in turn, there is once again collision by the third one at t = t + 60°. Thus, the interaction period of vortex γ is 1/6T, which corresponds to the typical frequency 6f R of pressure pulsation at the leading point DP1 of diffuser in Figure 20C. In addition to 0.8Q d operating condition, it can be seen that the frequency f BPF is prominent in other conditions in Figure 20C. So there may be impacting vortical structures at different flow rates.
Q criterion contours at 1.2Q d in diffuser I are shown in Figure 19. As seen in the figures, vortex α on the suction surface moves downstream to the diffuser outlet compared with other flow rates, and its evolution period is 1/3T. Besides, vortex β at the trailing edge gets weakened due to larger flow condition, and its shedding period is also 1/3T. Thus, it means that both vortical structures α and β have the same vortex shedding frequency approximately, which corresponds to 3f R of pressure pulsation in the trailing region points of diffuser in Figure 20B. Clearly the vortex χ occurs from the leading edge, then divides into multiple branches to  downstream, mainly including positive vortexes near the pressure surface and negative vortexes in the middle of the flow channel. It is found that the period of vortex χ is 1/6T, which is connected with the typical frequency of 6f R in Figure 20B.
Pressure pulsation intensities under off-design operating conditions are presented in Figure 21. There are significant differences among design and off-design conditions. It can be observed that the locations of high intensity are concentrated in the middle regions and trailing edge region near SS in diffusers at 0.8Q d . According to the above vortical structures evolution analysis, the significant bluish areas are consistent with shedding regions of vortexes α and β. Thus, it is noted that the periodic flow phenomena due to separation, vortex generation and shedding have significant effects on pressure pulsations in diffusers under 0.8Q d . From comparison with intensities at the design and low flow rates, the values are weakened at 1.2Q d , especially in the trailing region near the SS, which is closely related to reduced energy of vortexes. It is considered that the vortex χ has a great effect on formation of high intensity region at the leading edge of the PS except for the rotorstator interaction effect. From the above analysis, it is concluded that transient vortical structures indeed play an important role on internal pressure pulsation of the MSP.

| CONCLUSIONS
In the present study, the complex flow structures and the induced pressure pulsations of the MSP are analyzed by the numerical method. The storage tank is introduced during the calculation to take the inlet flow condition into consideration. Some conclusions are summarized as follows.
The most prominent discrete frequencies in the pressure spectrum are the diffuser blade passing frequency (f DPF ), the impeller blade passing frequency (f BPF ) and their high harmonics due to RSI. The amplitude of f DPF and f BPF on the pressure surface of impellers and diffusers are the largest in the core interaction region. The pressure pulsation RMS value at the corresponding position of the two-stage impellers and diffusers has a similar trend, but the secondary RMS value is slightly larger than the first stage.
Periodic vortex shedding exists at the trailing edge of diffusers under different operating conditions, corresponding to the pressure pulsation characteristic frequency of 2~3f R at that position. At the same time, vortex shedding of frequency (3f R ) is observed on the suction surface of diffusers. Compared with the design condition, the initial positions of the shedding vortex move forward to the center of the suction surface at 0.8Q d and to the trailing edge at 1.2Q d . There is a periodic hitting vortex evolution process of 6f R in the flow transition region between the impeller and diffuser, which is much obvious at 0.8Q d . Under 1.2Q d , a separating vortex of 6f R is captured at the leading edge of pressure surface in the diffuser, then it dissipates to the downstream, which could be attributed to excessive flow angle.
Through the analysis of the pressure pulsation intensity, it is clarified that the regions of high-pressure pulsations are consistent with that of vortexes shedding, developing and moving. Thus, the evolutions of typical vortical structures are crucial for the pressure pulsations of the MSP.
In this current research, the unsteady vortical structures and the corresponding effects on pressure pulsations in the MSP are revealed. We expect that the results will be useful for the subsequent optimization of the MSPs and other similar pumps. Furthermore, the structural stability, fluid dynamic induced excitation for the MSP with long-shaft system can be considered in future research.