A Mathematical Model of the Hydrodynamic Pressure Acting on a Turbine Runner Blade under the Rotor-Stator Interaction Based on the Quasi-Three-Dimensional Finite Element Method

In this paper, a turbine runner blade under hydraulic excitation was taken as the object of study and the mathematical hydrodynamic pressure model of a turbine runner blade under the rotor-stator interaction was established using quasi-three-dimensional ﬁnite element model. The regularity of the distribution of the hydrodynamic pressure was then investigated. First, the quasi-three-dimensional ﬁnite element model of a runner blade was established using the quasi-three-dimensional theory. Next, according to the pressure distribution from the numerical simulation, the mathematical mean pressure model of a node was established using the relative pressure diﬀerence method; the pressure ﬂuctuation was then established taking the rotor-stator interaction into consideration. Then, based on the mean pressure and the pressure ﬂuctuation of a node, the mathematical hydrodynamic pressure model for any position of a runner blade was obtained. Finally, the feasibility of the mathematical model was veriﬁed by an example analysis, and the internal relationship between the hydrodynamic pressure and its hydraulic parameters and the structural parameters of the blade was revealed, which provides a theoretical basis for further study of the dynamic characteristics of turbine runner blades under the rotor-stator interaction.


Introduction
Runner blades and guide vanes are the significant flow components of a turbine, and the relative rotation of the runner blades and the guide vanes results in the rotor-stator interaction phenomenon; this also creates the flow field to produce periodic oscillations and causes abnormal pressure pulses [1]. ese abnormal pressure pulses can cause the runner blades to vibrate intensely, and this can cause fatigue cracks to develop in the runner blade, which can seriously endanger the operation of the unit [2]. erefore, it is necessary to study the hydrodynamic pressure acting on a runner blade under the rotorstator interaction in the process of studying the dynamic characteristics of a runner blade.
Currently, the hydrodynamic pressure acting on turbine runner blades under the rotor-stator interaction has mainly been studied through the use of numerical simulation and experimental testing. When numerical simulation has been used to study the hydrodynamic pressure of runner blades under the rotor-stator interaction, the pressure fluctuation under different working conditions was mainly studied through the use of finite element software such as CFD (Computational Fluid Dynamics). Wang et al. [3] used CFD to investigate three-dimensional unsteady turbulent flow in the numerical simulation of a turbine. ey obtained a simulation of the pressure fluctuation under the rotor-stator interaction and used the double-amplitude (peak value) of the mixing frequency to calculate the variation of the pressure fluctuations acting on the runner blade. Anup et al. [4] used CFD to carry out a transient numerical analysis of the rotor-stator interaction in a Francis turbine at 100% of the optimum load and 72% part load. In their work, a periodical behavior of the pressure distribution was observed acting on the runner blades. When experimental testing has been used to study the hydrodynamic pressure acting on runner blades under the rotor-stator interaction, the influence of the rotor-stator interaction on the amplitude and frequency of the hydrodynamic pressure acting on the runner blades under different working conditions was studied. Tanaka et al. [5] calculated the transient hydrodynamic pressure of turbine runners with a very high head under the rotor-stator interaction based on both theoretical and experimental investigations. Rivedi et al. [6] studied the generation of unsteady pressure fluctuations in hydraulic turbines. e varying pressure amplitude of turbine blades with time was measured experimentally and simulated numerically. Chalghoum et al. [7] analyzed the characteristics of both the time domain and frequency domain of a pressure pulse under different coupling conditions. e Fast Fourier transform was used to obtain the spectra of the pressure pulse, and they compared the pump's performance curves that were obtained numerically with those obtained experimentally. Fu et al. [8] analyzed the cavitation characteristics of the runner blade pressure and optimized the operation of existing turbines for better fish passage conditions. Liu et al. [9] investigated the tip leakage vortex evolution mechanism and associated unsteady flow characteristics of a mixed flow pump under the design flow rate and proposed a theoretical prediction model of the primary tip leakage vortex trajectory. Han et al. [10] used ANSYS CFX to investigate the complicated structure of tip leakage vortex in an evolution period for a mixed flow pump as turbine at pump mode, obtained its dynamic and kinematic characteristics. Liu et al. [11] established the theoretical model of energy performance prediction for centrifugal pump as turbine under pump and turbine modes based on detailed analysis of losses and proposed a flowrate-based iteration method to determine the best efficiency point under turbine mode. e variations in the hydrodynamic pressure were mainly studied under different working conditions through the use of numerical simulation and experimental testing. However, a mathematical model of the hydrodynamic pressure containing both hydraulic parameters and structural parameters was not obtained, so it was difficult to analyze the hydrodynamic pressure of the transient process of runner blades [12]. Nevertheless, the mathematical model can reflect the internal relationship between the hydrodynamic pressure, the hydraulic parameters, and the structural parameters, which can analyze the hydrodynamic pressure of the transient process of runner blades. e transient process analysis is the basic work of studying the dynamic characteristics of runner blades. However, studies based on this aspect of the research have rarely been published. e quasi-three-dimensional finite element method can be used not only to design Francis runner blades, but also to effectively analyze the flow on the runner blades [13]. A turbine runner blade was taken as the object of study in this paper. According to the characteristics of the quasi-threedimensional design of turbine runner blades, the quasithree-dimensional finite element model of the blade was established, and the mathematical hydrodynamic pressure model of a turbine runner blade under the rotor-stator interaction was then obtained. e regularity of the distribution of the hydrodynamic pressure was taken as an example to be studied.

The Quasi-Three-Dimensional Finite Element Model of a Turbine Runner Blade
e global coordinate system of a turbine runner blade was established using the cylindrical coordinate system, as shown in Figure 1. e origin was located at the intersection of the lower surface of the upper crown and the axis of the runner blade; the three variables were R, θ, and Z, where R denotes the radius of the blade, θ denotes the rotation angle at an instantaneous moment, and Z denotes the depth. According to the quasi-three-dimensional design method of a turbine runner blade [14], it was assumed that the flow of water on the runner blade was incompressible, inviscid, and steady. Based on the flow streamline, the blade can be divided into (m + 1)(n + 1) units using m + 1 radial lines and n + 1 axial lines. e node at the intersection of the (i+1)th axial line C i0 C im ⌢ and the (j+1)th radial line C 0j C nj ⌢ was C ij , and its coordinates were (r ij , θ ij , z ij ). e region between the radial lines C 0j C nj ⌢ , C 0(j+1) C n(j+1) ⌢ and the axial lines was defined as unitJ ij ; the nodes at this unit were C ij , C (i+1)j , C (i+1)(j+1) , C i(j+1) . e hydrodynamic pressure at (r, θ, z) in unit under the rotor-stator interaction can be divided into two terms: the mean pressure p J ij (r, θ, z) and the pressure fluctuation p J ij ′ (r, θ, z) that are caused by the rotor-stator interaction [7]; that is, e pressure p J ij (r, θ, z) can be calculated from where p ij , p (i+1)j , p (i+1)(j+1) , and p i(j+1) are the pressures of C ij , C (i+1)j , C (i+1)(j+1) , and C i(j+1) and N 1 , N 2 , N 3 , and N 4 are the shape functions of the pressures at C ij , C (i+1)j , C (i+1)(j+1) , and C i(j+1) . According to the quasi-three-dimensional design method of a turbine runner blade [14], and based on the regularity of the mean pressures, the shape functions of the pressures can be obtained by stream curve of the turbine runner blade. Based on the different surface modeling techniques of the runner blade, the shape functions can be obtained by many stream curve parameterization methods [15], including Ferguson parameterization method, Coons parameterization method, Bezier parameterization method, and B spline curve method.

e Flow Velocity Based on the Quasi-ree-Dimensional
Finite Element Model. According to Bernoulli equation, the hydrodynamic pressure at any position is related to the water head and the flow velocity at that point, so it is necessary to calculate the flow velocity of the runner blade of quasi-threeelement finite element model. Since the relative pressure difference method is adopted, the reference point needs to be set. According to the characteristics of the movement of the turbine runner system, the reference point is usually set at the central point of the inlet section of the volute [16], as shown in Figure 2. e pressure of the reference point is related to the average flow velocity at the central point of the inlet section of the volute v c , which can be expressed as [17] v c � α where α is the velocity coefficient at the inlet section of the volute; for a metal spiral case, α � 0.7 − 0.8. e velocity moment of the outlet edge can be determined by empirical formula, and the velocity moment of the inlet edge can be obtained by fundamental equation of the hydraulic excitation. e variation of velocity moment from the inlet edge to the outlet edge can be determined by exponential function [18]. Based on the quasi-three-dimensional method of a turbine runner blade and the velocity moment distribution regulation in the flow line, the velocity moment of node C ij in Figure 1 v ij r ij can be expressed as follows [18]: where κ is a constant that is determined by the blade profile of the turbine [19]; the pressure side constant is different from the suction side. λ ij � L ij /l i ; L ij is the length of the arc C 0j C ij ⌢ , and l i is the length of the arc C 0j C nj ⌢ . v 0j r 0j , v nj r nj represent the velocity moment of the inlet edge node C 0j and the outlet edge node C nj in the streamline C 0j C nj ⌢ , respectively. According to equation (4), the absolute velocity of node C ij in Figure 1 where v 0j , v nj represent the absolute velocity of the inlet edge node C 0j and the outlet edge node C nj in the streamline C 0j C nj ⌢ , respectively, which can be expressed as [18] v where z a , z a ′ are the Z-axis coordinates of the inlet water edge and the outlet water edge on the central streamline of the runner blade, respectively; h, h ′ are the heights at the inlet water edge and the outlet water edge on the central streamline of the runner blade respectively; and v a , v a ′ are the circumferential components of the absolute velocity at the inlet water edge and the outlet water edge on the central streamline of the runner blade, respectively; these can be determined according to the flow rate, the volume efficiency, the guide vane height, the outlet flow angle of the guide vane, the cross-sectional area at the outlet edge of the runner blade, and other parameters [20]. Substituting equations (6) and (7) into (5), and rearranging allow the following to be obtained: where Λ ij , Λ ιj ′ are the quantities related to the turbine structural parameters and position parameters, and

e Mean Pressure at a Node of the Turbine Runner Blade.
Using the relative pressure difference method, the mathematical expression of the mean pressure of the reference point and the relative pressure difference between the node and the reference point were established, and then the mathematical mean pressure model at a node of the runner blade was then obtained. e quasi-three-dimensional finite element method can be used not only to design Francis runner blades, but also to effectively analyze the flow on the runner blades [13]. According to the characteristics of the quasi-three-dimensional finite element model of the runner blade, the finite element was divided based on the stream curve function Figure 1: e global coordinate system of a turbine runner blade.
Mathematical Problems in Engineering 3 [14]. So, Bernoulli's equation can be used when the hydrodynamic pressure distribution regulation is studied. According to Bernoulli equation, the mean pressure of the reference point in Figure 2 can be expressed as where ρ is the water density, g is the acceleration due to gravity, and H is the water head of the turbine. Substituting equation (3) into (10), and rearranging, the following can be obtained: According to Bernoulli's equation, and because a runner blade is close to the reference point, the energy loss between a runner blade and the reference point can be ignored; the relative pressure difference Δp ij between the node C ij and the reference point will be where z ij , z c denote the Z-axis coordinates of node C ij and the reference point, respectively, andv ij represents the absolute velocity of node C ij . Substituting equations (3) and (8) into equation (12)  results in e mean pressure of the node C ij is the sum of the pressure of the reference point and the pressure difference between the node C ij and the reference point; that is, where From equation (14), it can be seen that the mean pressure of the node of a runner blade is not only related to its positional parameters, such as the Z-axis coordinates, the radius, the length of the arc, etc., but also related to the hydraulic parameters, such as the flow rate, the water head, the velocity, etc., as well as the structural parameters, such as the height of the guide vane, the outlet flow angle of the guide vane, the outlet flow angle of the runner blade, etc.

e Pressure Fluctuation at a Node of the Turbine Runner
Blade under the Rotor-Stator Interaction. Under the rotorstator interaction, the absolute velocity of the runner blade contains a periodic pressure fluctuation term [6], the frequency of which is the multiple passing frequency of a runner blade. e pressure fluctuation p ij ′ of the runner blade node C ij under the rotor-stator interaction can be expressed as where K is a natural number, K � 1, 2, . . . G; because the higher the frequency order of the pressure fluctuation, the smaller the amplitude of the pressure fluctuation, G � 2 is a desirable value in the actual runner blade system. Z r denotes the number of runner blades, and k ij denotes the variation coefficient of the amplitude of the pressure fluctuation at node C ij , which can be determined based on the profile of the turbine blade [3]. φ K is the phase of the pressure fluctuation when there are K multiples of the passing frequency, which can be determined from the initial position of the guide vane and the runner blade [6].B K is the amplitude of the pressure fluctuation when there are K multiples of the passing frequency, which can be calculated from the average velocity of the reference point [14] and can be expressed as The reference point The volute The guide blade The turbine runner where δ K is the pressure fluctuation coefficient when there are K multiples of the passing frequency, which can be determined from the empirical value of the experimental data and the simulation data [20]. Substituting equation (16) into equation (15), and rearranging, the following can be obtained: It can be seen from equation (17) that the pressure fluctuation under the rotor-stator interaction has a periodic function, which is related to the number of runner blades, the water head, and the speed of rotation.

e Mathematical Hydrodynamic Pressure Model of a Turbine Runner Blade.
e hydrodynamic pressure at node C ij of a turbine runner blade under the rotor-stator interaction is the sum of the mean pressure and the pressure fluctuation under the rotor-stator interaction; that is, Substituting equations (15) and (17) into (18), and rearranging, the following can be obtained: From equation (19), the hydrodynamic pressure at the nodes C (i+1)j , C (i+1)(j+1) , C i(j+1) on the element J ij can be calculated. Substituting equations (19) into (2), the hydrodynamic pressure mathematical model for position (r, θ, z) on the element J ij at t second can be written as As can be seen from equation (20), the hydrodynamic pressure is a function of both time and space. e pressure is not only related to the position parameters, such as the Zaxis coordinate, the radius, the length of the arc, etc., but also related to the hydraulic parameters, such as flow rate, water head, speed of rotation, etc.; and it is also related to the structural parameters, such as the height of the guide vane, the outlet flow angle of the guide vane, the outlet flow angle of the runner blade, etc. It reflects the relationship between hydrodynamic pressure and structural parameters, hydraulic parameters. Transient process analysis can be carried out by changing hydraulic parameters.

Example.
In this paper, the feasibility of the mathematical model has been verified using an example. e internal relationship between the hydrodynamic pressure acting on runner blades, the hydraulic parameters, and the structural parameters under the rotor-stator interaction has been revealed. e regularity of the distribution of the hydrodynamic pressures at different positions on the runner blades under the rotor-stator interaction has been studied, and the transient process of the hydrodynamic pressure has been analyzed.
According to the quasi-three-dimensional finite element model of the Francis-99 turbine runner blade characteristic, the shape functions of the pressures on the streamline direction can be expressed by cubic B spline curve [21]. Due to the different design parameters of the turbine runner blade, the shape of cubic B spline curve is also different, resulting in the difference of the shape functions values in equation (2); the general form of cubic B spline curve function can be expressed as [22] where ϖ is the shape change parameter of cubic B spline curve; the value is related to the quasi-three-dimensional design parameters of the runner blade. When ϖ � 0, it is cubic uniform B spline curve. According to the characteristics of a Francis-99 turbine runner blade from the literature [23], ϖ � 0.5. W is the correlation coefficient between hydrodynamic pressure and position, and Mathematical Problems in Engineering where δ ′ is the comprehensive variation coefficient along the axial and radial directions and can be obtained on the basis of test results or simulation results and 2a, 2b are the length and width of the element, respectively. Substituting equation (21) into (20) and rearranging, the hydrodynamic pressure mathematical model for position (r, θ, z) can be obtained: e pressure surface of a Francis-99 turbine runner blade from the literature [23] was selected as the object of this study, some parameters of which can be obtained through the website [24]. e parameters are as follows: the rated water head was H � 11.91 m, the rated flow rate was Q � 0.203 m 3 /s, and the rated speed of rotation was 335.4 r/min. e velocity coefficient of the volute inlet section was α � 0.75, and the Z-axis coordinate of the central section of the volute inlet section was z c � 0.0216 m. e height of the guide vane was b 0 � 0.0224 m, and the outlet flow angle of the guide vane was α 0 � 13.2°. e diameter of the runner blade was D 0 � 0.349 m. e Z-axis coordinate of the inlet water edge on the central streamline of the runner blade was z a � 0.0488 m, and the Z-axis coordinate of the outlet water edge on the central streamline was z a ′ � 0.192 m; the height at the inlet water edge on the central streamline of the runner blade was h � 0.0596 m; the height at the outlet water edge on the central streamline was h ′ � 0.0576m. e number of runner blades was Z r � 30; there were 15 long leaves and short leaves among them. Other parameters can be obtained through experimental data and simulation data of previous similar research: the pressure side of the blade profile constant was κ � 0.56, and the suction side one was κ ′ � 0.89 [18]. e rated working condition was calculated as the working condition. e first two orders of the pressure fluctuation coefficient under the rotor-stator interaction were δ 1 � 1.2 × 10 − 2 , δ 2 � 2.3 × 10 − 3 , the first two phases of which wereφ 1 � π/6 and φ 2 � 5π/13 [6,23]. e variation coefficient of the amplitude of the pressure fluctuation was k ij � (1 − 0.2L ij /l i )(1 + 0.2z ij /z mj ) [25,26]. e comprehensive variation coefficient was δ ′ � 0.1 [23]. e elements quantity can influence calculation results, and node VL01 was selected to study the relationship between elements quantity and the hydrodynamic pressures simulation results. e specific results are shown in Table 1.
It can be seen from Table 1 that the higher the elements quantity divided, the smaller the error of the hydrodynamic pressure simulation results. When the elements quantity exceeded 11, the error of the hydrodynamic pressure simulation results of node VL01 was less than 5.6%, which was within the allowable range.
e No. 4 way of elements quantity was selected; the runner blade was then divided into 32 elements by 5 radial lines and 9 axial lines: the 1st radial line was located in the upper crown; the 5th radial line was located in the lower ring; the 1st axial line was located in the inlet edge; the 9th radial line was located in the outlet edge. A structural sketch of the long blade is shown in Figure 3. e nodes in the inlet edge were P00, P01, P02, P03, and P04; the nodes in the outlet edge were P80, P81, P82, P83, and P84.

Verification of the Mathematical Model.
According to equation (19), the hydrodynamic pressures at nodes VL01, P42, and P71 were calculated. e time domain simulation curves of the hydrodynamic pressures and the experimental data from the literature [23] are shown in Figures 4-6.
As shown in Figures 4 and 6, the error in the hydrodynamic pressure of node VL01 was 4.78%. However, at nodes P42 and P71, the pressure fluctuation in the experimental data was caused not only by the rotor-stator interaction, but also by draft tube back-flow, so the errors were relatively larger. e errors in the hydrodynamic pressure of nodes P42 and P71 were 9.17% and 8.78%, respectively. From the data comparison, the suitability of equation (19) in calculating the hydrodynamic pressure of nodes can be proved, and the feasibility of the mathematical model can be illustrated.

e Regularity of the Distribution of Hydrodynamic Pressures at Different Positions.
In order to research the dynamic characteristics of runner blades under the rotorstator interaction, the regularity of the distribution of the hydrodynamic pressures at different positions has been studied. e selected nodes were P01, P02, and P03 at the inlet edge; P81, P82, and P83 at the outlet edge; and P02, P22, P42, P62, and P82 at the quasi-three-dimensional streamline.
According to equation (19), the hydrodynamic pressures of nodes P01, P02, and P03 at the inlet edge and nodes P81, P82, and P83 at the outlet edge were then calculated. e mean pressures are shown in Figure 7 and the amplitudes of the pressure fluctuations are shown in Figure 8. It can be seen from Figures 7 and 8 that the mean pressures gradually decreased in the axial direction, and the amplitudes of the pressure fluctuations caused by the rotor-stator interaction gradually increased in the axial direction.
According to equation (19), the hydrodynamic pressures of nodes P02, P22, P42, P62, and P82 at the quasi-threedimensional streamline were then calculated. e mean pressures are shown in Figure 9 and the amplitudes of the pressure fluctuations are shown in Figure 10. It can be seen from Figures 9 and 10 that the mean pressures and the pressure fluctuations caused by the rotor-stator interaction gradually decreased along the streamline.
Combining the regularity of the mean pressure and the pressure fluctuations caused by the rotor-stator interaction in the axial direction and the streamline direction, it can be known that the maximum mean pressure occurred near the upper ring at the inlet edge (1.772 × 10 5 Pa). e maximum pressure fluctuation caused by the rotor-stator interaction occurred near the lower ring at the inlet edge (1.12 × 10 3 Pa), which also provided a reference for analyzing the dynamic characteristics of runner blades.
According to equation (23), the hydrodynamic pressure diagram of the runner blade pressure side and the suction side at t � 8 s under the rotor-stator interaction could be obtained as shown in Figures 11 and 12, respectively. It can be seen from Figures 11 and 12 that the instantaneous hydrodynamic pressure gradually decreased in the axial  Figure 3: A structural sketch of the runner blade.    direction and the instantaneous hydrodynamic pressure gradually decreased in the streamline direction. Compared with the diagram from the CFD simulation of the hydrodynamic pressures of the runner blade under the rotor-stator interaction from the literature [23], the maximum error occurred near the lower ring of the outlet edge, which was 9.92%. As the pressure fluctuations in the experimental data were caused not only by the rotor-stator interaction, but also by draft tube back-flow, the errors were relatively larger. erefore, the feasibility of the mathematical model has been illustrated again.

e Hydrodynamic Pressure Analysis of the Transient
Process. In order to study the dynamic characteristics of the runner blades under the rotor-stator interaction, the regularity of the hydrodynamic pressure in the transient process has been studied. For comparison with the experimental data in the literature [23], the selected nodes were P03 at the inlet edge, P42 near the middle, and P71 near the outlet edge. Four working conditions were selected as follows: part load 1, the head was H � 12.29 mand the flow rate was Q � 0.071 m 3 /s; part load 2, the head was H � 12.00 m and the flow rate was Q � 0.163 m 3 /s; best efficiency point (BEP), the head was H � 11.91 m and the flow rate was Q � 0.203 m 3 /s; high load, the head was H � 11.84 m and the flow rate was Q � 0.221 m 3 /s. According to equation (19), the hydrodynamic pressures at node P03 were then calculated. e mean pressures are shown in Figure 13 and the amplitudes of the pressure fluctuations are shown in Figure 14. As can be seen from Figures 13 and 14, with the increase of the working conditions, the mean pressure at node P03 decreased, while the pressure fluctuations at node P03, which were caused by the rotor-stator interaction, increased.
According to equation (19), the hydrodynamic pressures at node P42 were then calculated. e mean pressures are shown in Figure 15 and the amplitudes of the pressure fluctuations are shown in Figure 16. As can be seen from Figure 15, with the increase in the working conditions, the mean pressures at node P42 increased at the start and then decreased later. e calculation of equation (19) shows that the mean pressure of P42 reached its maximum at part load of H � 12.00 m and Q � 0.169 m 3 /s. As can be seen from Figure 16, with the increase in the working conditions, the pressure fluctuations at node P42 decreased at the start and then increased later. e pressure fluctuations at node P42 reached their minimum at part load of H � 12.00 m and Q � 0.169 m 3 /s. According to equation (19), the hydrodynamic pressures at node P71 were calculated. e mean pressures are shown in Figure 17, and the amplitudes of the pressure fluctuations are shown in Figure 18. As can be seen from Figure 17, with the increase of working conditions, the mean pressures of node P71 increased at the start and then decreased later. e calculation of equation (19) shows that the mean pressure of P71 reached the maximum at part load of H � 12.00 m and Q � 0.186 m 3 /s. As can be seen from Figure 18, with the    increase in the working conditions, the pressure fluctuations at node P71 decreased at the start and then increased later. e pressure fluctuations at node P71 reached their minimum at part load of H � 12.00 m and Q � 0.186m 3 /s. In this paper, the regularity of the distribution of the hydrodynamic pressures at different positions of the runner blades under the rotor-stator interaction has been studied, and the transient process of hydrodynamic pressure has been analyzed.

Conclusion
According to the characteristics of the quasi-three-dimensional design of turbine runner blades, a mathematical model of the hydrodynamic pressure acting on the blades under the rotorstator interaction has been established. e regularity of the hydrodynamic pressure in the transient process has also been revealed. e research has shown the following: (1) Under the same working condition, the pressure fluctuation caused by the rotor-stator interaction increased gradually from the upper crown to the lower ring, and the maximum value of the pressure fluctuation appeared at the position near the outlet water edge of the lower ring.