Numerical analysis of aerodynamic damping for a centrifugal impeller

ABSTRACT In order to realize the numerical analysis of aerodynamic damping for the centrifugal impellers, a numerical code integrating the flow-structure data transfer, grid deformation and flow simulation under the moving grid system is first developed. To decrease the huge consumptions of CPU time and memory space, the compactly supported radial basis function is adopted to carry out the data transfer and grid deformation, and the alternating digital tree technique is utilized to calculate the wall distances. By the test cases of an oscillating cascade and a centrifugal impeller, the correctness of the code for the flow simulations with moving boundaries and the flow fields in centrifugal impellers is validated. Then taking a centrifugal impeller as the research object, the calculations of aerodynamic damping characteristics were carried out. The results show that, the modal aerodynamic damping ratio has no relationship with the vibrating amplitude under small vibrations. For the disk-dominant vibration mode, the modal aerodynamic damping ratio decreases as the operating point shifts toward the stall point. The aerodynamic damping caused by the backward traveling wave vibration is larger than that by the forward traveling wave vibration.


Introduction
For modern industrial centrifugal compressors, the forced vibration of the impellers induced by gas excitations is an important factor causing the fatigue failures of impeller blades. In order to predict the forced response characteristics of centrifugal impellers, all mechanical parameters of the impeller component should be obtained, but the damping characteristic is an important data hard to be obtained. For unshrouded centrifugal impellers, the mechanical damping can be neglected because the blades and the wheel disk are integrally machined. In addition, some studies have shown that the blade material damping is small, and has little effect on the dynamic response. Thus the aerodynamic damping becomes the dominant factor affecting the damping characteristics of the impeller. The aerodynamic damping is used to describe the energy transfer effect between the impeller and the surrounding gas, and its value is affected by many flow parameters. However, presently the literature data on the aerodynamic damping measurement are limited, and existing studies are mainly focused on the axial blades. An earlier work by Crawley (1983) realized the measurements of the blade aerodynamic damping of a transonic compressor rotor. The blades were excited by the upstream inlet disturbance device. The experiments obtained the aerodynamic damping under different interblade phase angles. Kielb and Abhari (2003) measured the aerodynamic damping and the structural damping of the blades in a gas turbine. By subtracting the damping measured under the idle condition from the total damping under a certain operation condition, the corresponding aerodynamic damping of turbine blades can be obtained. Kammerer and Abhari (2009) measured the damping data of a high speed centrifugal impeller under different operating conditions, and the emphasis was put on the influences of inlet pressure and operating point on the modal damping data of the impeller. By using the similar experimental method, Zemp, Abhari, and Ribi (2011) studied the aerodynamic damping and the strain response of the blades in a centrifugal impeller. The results showed that when the inlet pressure was 1 bar, the aerodynamic damping for blade mode 1 was about three times that of the material damping. Seeley, Wakelam, Zhang, Hofer, and Ren (2017) built a linear turbine cascade. The blade was excited about its pitch axis by an electromagnetic actuation system over a range of frequencies and amplitudes. Experimental results indicated that the aerodynamic damping was independent of displacement amplitude within the tested range.
With the CFD techniques, there were also some studies adopted the numerical simulations to obtain the aerodynamic damping of turbomachinery blades. Grüber and Carstens (2001) carried out a parametric study investigating the influence of viscous effects on the damping behavior of vibrating compressor blades. An important result is that viscous effects may cause a significant change of aerodynamic damping. This behavior is demonstrated by two cases in which an Euler calculation predicts a damped vibration while a Navier-Stokes calculation leads to an excited vibration. Parthasarathy (2011) numerically calculated the aerodynamic damping of a transonic fan blade, and investigated the influence of inter-blade phase angles on the damping characteristics. Srivastava, Bakhle, and Keith (2003) calculated the blade aerodynamic damping by adopting the energy method, and analyzed the influence of back pressure on the aerodynamic damping, as well as the influence of the shock positions on the aerodynamic damping. Seinturier, Lombard, and Dumas (2004) calculated the aerodynamic damping of the second bending mode of a compressor blade respectively with the Euler equations and the linearized north-south (N-S) equations, and the obtained damping was utilized in the calculations of gas-excited forced responses. Robin et al. (2013) adopted the Fourier transformation method to solve the periodically varying flow field of NASA rotor 67, and only two blade passages were used in the calculation. The aerodynamic damping data corresponding to different inter-blade phase angles agreed well with the results by the full passage simulations. Sadeghi and Liu (2001) studied the flutter characteristics of a mistuning compressor cascade by the full-channel flow simulation. The results show that the mistuning of inter-blade phase angle exerts little effect on the aerodynamic damping, while the frequency mistuning has a significant influence on the aerodynamic damping. An appropriate frequency mistuning can improve the aeroelastic stability. Im and Zha (2013) investigated the flutter mechanism of a transonic fan rotor with a forward traveling wave vibration by a fully coupled fluidstructure interaction method. From the calculated blade displacement responses, the aerodynamic damping can be obtained. The results show that as the operating point shifts toward the stall, the aerodynamic damping gradually varies from positive to negative. Peruzzi et al. (2016) employed a non-linear uncoupled method for the flutter stability analysis to estimate the aerodynamic damping of a high-pressure steam turbine blade row. The obtained results are useful to give a better insight into the aeroelastic response of this type of blades.
The earlier computational studies of aerodynamic damping are all focused on axial flow blades, while few studies are carried out for the numerical analyses of the aerodynamic damping characteristics of centrifugal impellers. Due to the structural characteristic of rotor disks, the blades usually exhibit traveling wave vibrations. Especially for the high and short structures like axial flow blades, it can be assumed that all the blades are fixed at the roots, each blade vibrates with the same frequency and amplitude, and the vibrations between two adjacent blades only differ by a phase angle. This approach has been widely utilized in the calculations of aerodynamic damping for axial flow blades. But for the centrifugal impeller structure, the blade-disk vibration must be analyzed as a whole because the vibrations of the blades and the wheel disk are closely coupled. The vibrations of different nodal diameters are different not only in the interblade phase angle, but also in the frequency. Therefore, the analysis of aerodynamic damping for the centrifugal impellers is more difficult than that for the axial flow blades.
The calculation of aerodynamic damping is a kind of aeroelastic problem which involves multiple numerical modules, such as fluid-structure data transfer, grid deformation and flow simulation with moving boundaries. However, how to integrate these modules into one code is tough work. So far, besides the earlier mentioned studies on the aerodynamic damping, there are also many aeroelastic studies contriving to realize this work (Alder, 2015;Fairuz, Abdullah, Yusoff, & Abdullah, 2013;Im & Zha, 2012;Keye, 2011;Sadeghi, Yang, Liu, & Tsai, 2003;Zhang, Ding, Ji, & Zhang, 2016;Zhang, Jiang, & Ye, 2007). But presently few commercial software can achieve a seamless integration of multiple modules. Furthermore, to calculate the aerodynamic damping the numerical modeling of the impeller will involve multiple blade passages, so the grid scale is usually very large. Presently most commercial software lack efficient solutions for the large-scale grid deformation and wall distance calculation, which leads to the huge consuming of CPU time and memory space, and hence works against the numerical investigation of aerodynamic damping for centrifugal impellers.
A numerical code integrating data transfer, grid deformation and flow analysis is first developed in this work to realize the numerical calculations of aerodynamic damping for the centrifugal impellers. By using the compactly supported radial basis function and the alternating digital tree technique, the consumptions of computational resources for grid deformation and distance calculation are greatly reduced. Taking a centrifugal impeller as the research object, through the calculation of its modal aerodynamic damping ratios, the influences of vibration amplitude and operating condition on the aerodynamic damping are analyzed, and then the distributions of aero power density on the impeller surface are discussed. Finally, the aerodynamic damping characteristics of the impeller under the traveling wave vibrations are studied.

Governing equations
To solve the flow field with the finite volume method, the compressible Reynolds-averaged N-S equations are converted from the relatively rotating system to the arbitrary curvilinear coordinate system (ξ , η, ζ ).
(1) whereQ denotes the conservative quantities, andF v , G v ,Ĥ v are viscous fluxes. The inviscid fluxesF,Ĝ,Ĥ and the source termŜ can be expressed aŝ where J denotes the determinant of Jacobian matrix. ω is the angular velocity around the z axis. p is the pressure, ρis the density, and E is the total energy. U, V, W are the covariant velocities, which can be expressed as where u, v, w denote the absolute velocity components, and w is the z-directional velocity. x t ,y t ,z t denote the moving velocity components of cell faces.
In the code the inviscid terms can be discretized using the second-order central scheme and the upwind scheme. The viscous terms are discretized using the central scheme. Considering that current PCs and workstations have widely adopted the multi-core and multithread processors, the shared memory parallel compilation technology of OPENMP is introduced in the programing. By adding a small amount of pseudo code to the program, the single machine multi-thread parallel computation has been achieved, which significantly reduces the cost of time.

Inter-block data exchange for multi-block mesh
Because of the low efficiency of space filling for unstructured grids, the program is developed based on the multiblock structured grids in order to adapt to the complex geometry of the impeller passage. For the multiblock grid topology, the inter-block interface is a nonphysical boundary, on which the boundary conditions are imposed by data exchange between the two adjacent blocks. At present, only the full matching connection is realized in the code. As shown in Figure 1, for the cell faces on the inter-block interface, to construct the spatial discretization scheme the cell-center flow data on both sides of the interface should be obtained. In this case, the two-layer cells on either side of the inter-block interface can be regarded as the dummy cells of the opposite block. Thus, the inter-block data exchange can use the flow variable values on the real cells.

Discretization of viscous terms
The viscous terms are discretized by the central difference. The second derivative of flow variable on the cell center is calculated as the difference of the first derivatives on the cell faces. Based on this principle, the i-directional difference of viscous flux can be written as For a multi-block structured grid with the topology shown in Figure 2, when applying the Gauss theorem on the offset control volume to compute the cell-face velocity derivatives and temperature gradients, the interblock data exchange of geometric information will lead to ambiguity. Usually the tip clearance grid at the blade leading edge possesses this type of topology. Under this circumstance, the Jacobi transform method can be used to calculate the derivatives on the cell faces. Take the idirection cell face as an example, the first derivative can be calculated as where φ denotes the physical quantities such as velocity and temperature. The value of φ on the half point can be interpolated by the values on the surrounding cell centers. For instance

Temporal discretization
The unsteady flow calculation is achieved by adopting the dual-time method (Jameson, 1991). A pseudo-time iteration process is introduced at each physical time step. When the flow field is marched to a steady state in the pseudo-time domain, the flow solution at each physical time step is obtained. The iteration can be expressed as where τ and t respectively denote the pseudo time and the real time. The superscripts m and n respectively denote the pseudo time step and the real-time step. Time marching in the pseudo time domain adopts the LU-SGS implicit algorithm. For the flow simulations with moving boundaries, in order to avoid the errors of conservative variables introduced by the deformation of control volumes, the geometric conversation law (Sitaraman & Baeder, 2006) must be satisfied.
where 1/J equals the cell volume, and represents the summation over cell faces. The physical meaning of the geometric conservation law is that the volume increment caused by the grid deformation equals the net volume swept by the moving cell faces.

Grid deformation
The aerodynamic damping is caused by the impeller vibration, but the calculation of the impeller vibration is based on the finite element mesh. In order to obtain the aerodynamic damping, the flow field surrounding the deforming impeller must be solved. Therefore the deformation data should be transferred from the structural surface to the aerodynamic surface. The study in this work is focused on the aerodynamic damping caused by the modal vibration of the impeller. To this end, first the nodal displacements of the structural surface under a certain mode are transferred to the aerodynamic surface points, and then the deformation of the aerodynamic mesh is carried out. Flow-structure data transfer and grid deformation are generally considered to be of a different nature. However, with the help of radial basis function, the two problems can be easily handled in the same way.
The interpolation problem based on the radial basis function (Allen & Rendall, 2007;Buhmann, 2003) can be expressed as follows where s(x)is the function of point coordinate x, representing the amount of deformation here. ϕ is the radial basis function, and x i denotes the known structural surface nodes. The p(x) is a polynomial function, which can be expressed as the following linear polynomial.
The displacement equations of N structural nodes can be written as To determine the interpolation coefficient α i ,γ 0 ,γ x ,γ y and γ z , the following additional equations should be supplemented.
where q(x) denotes the polynomials with the degree less than or equal to that of p(x). However, the radial basis function method has the problem of large consumptions of computational time and memory space. In order to make the matrix present a sparse and banded distribution and thereby beneficial to solving the large-scale data interpolation problem, Wendland (1995) defines the compactly supported radial basis function as follows where k is a positive integer. g(r) is a polynomial function of r, and r = (||x − x i ||/R), where R is the compactly supported radius. A reasonable value of R can ensure the interpolation accuracy and meanwhile obtain a good matrix sparsity. For the aerodynamic grid deformation of the impeller, besides the impeller surface, the constraint of the casing surface should also be taken into account. After the displacements of structural nodes are transferred to the aerodynamic grid points, the attenuation control of grid deformation should be implemented as follows where d r and d s respectively represent the distances from the grid point to the impeller surface and the casing surface, and σ represents the attenuation rate.

Fast calculation of wall distance
For the analysis of unsteady flows with moving boundaries, after the grid deformation at each time step the minimum distance from the center of each cell to the wall should be calculated. The wall distances will be used in the turbulence calculation. At present, the wall distance is usually calculated by the direct method, which calculates the distance of each cell to all the wall points. Then through point-by-point comparison, the minimum value is taken as the wall distance of that cell. For largescale meshes, the amount of computation for the direct method is very large, which greatly reduces the computational efficiency of unsteady flows. In this work, a kind of alternating digital tree searching technique (Boger, 2001;Bonet & Peraire, 1991) is adopted to realize the fast calculation of wall distance. Take the two-dimensional space as an example, the specific implementation steps are as follows: Establish a root node A for the first wall point, which points to a rectangular region R containing the flow field. Then the region is equally divided into two parts along the x coordinate, one of which named R l is the left half of R, and the other named R r is the right half of R. For the second wall point, find the subregion (R l or R r ) which containing that point, and establish a child node B pointing to that subregion. Then the region is equally divided into two subregions along the y coordinate. For each of the followup wall points, traverse the current digital tree, search the minimum subregion containing that point, and establish a child node pointing to that subregion. Again the region is divided alternately along the x or y coordinate. Repeat these procedures until the digital tree containing all the wall points is built. The building process can be shown in Figure 3.
Assume the minimum distance from the grid point x 1 to the wall points is d 1 , then the minimum distance d 2 from an arbitrary grid point x 2 to the wall points must satisfy d 2 ≤ d 1 + d 1-2 , where d 1-2 is the distance from x 1 to x 2 . Let r = d 1 + d 1-2 , then the closest wall point to x 2 must be located in the sphere with the center of x 2 and the radius of r. For the convenience of calculation, the search range can be expanded from the sphere to its outer tangent box Q.
After the search range is determined, the earlier built digital tree is traversed. The specific steps are as follows: If the corresponding wall point of current node is located in the rectangular box Q to be searched, then the distance from x 2 to that wall point is calculated. If the left child node is not null, and the left subregion R l of current node intersects with the box Q, then the left subtree is traversed. If the right child node is not null, and the right subregion R r of current node intersects with the box Q, then the right subtree is traversed.

Oscillating cascade
The cascade adopts the NACA 0012 airfoil. The inlet Mach number is 0.593, the Reynolds number is 3.21 × 10 6 , and the inlet airflow angle is 0°. The pitch chord ratio is 1. The blade oscillates around its leading edge with a torsional amplitude of 2°and an inter-blade phase angle of 0°. The oscillation frequency is 64 Hz. The inlet and outlet sections both adopt the H-type grid, while the blade section adopts the C-type grid. The grid sizes of each block are respectively 49 × 13 × 9, 57 × 37 × 9 and 129(around the blade) × 29(blade-to-blade) × 9(span).
In the calculation, the grid points on the blade surface are set to vibrate with the blade, and with the increase of distance from the blade surface the grid deformation attenuates gradually. At the inlet, outlet and periodic boundaries the grid deformation decays to zero. Usually there exists a phase lag between the flow responses and the oscillation of the cascade. The calculated first-order unsteady pressure coefficients are shown in Figure 4, where the abscissa represents the relative chord. The real part represents the variation amplitude multiplied by the cosine of the phase, while the imaginary part represents the variation amplitude multiplied by the sine of the phase. It is obvious that the calculated curves agree well with the data by Huff (1987), and the results of the real part are almost the same. The results verify the correctness of the code to simulate the flow field with moving boundaries.

Subsonic radial impeller
The test case targets to simulate the unsteady flow of a subsonic radial impeller. The rotational speed is 22,790 rpm. The impeller has 16 blades with a diameter of 175 mm and an outlet width of 10.5 mm. As shown in Figure 5, the mesh consists of four H-type blocks and the total grid number is 550420. The flow field is calculated with an inlet total pressure of 101.3 kPa, a total temperature of 293.15 K and an average back pressure of 130 kPa. Figure 6 shows the calculated surface pressure distributions at the blade mid-span. It can be seen that the  surface pressure curves obtained by this code are in good agreement with those by the commercial software. At the front of the blade the pressure on the pressure side is lower than that on the suction side, which indicates that the impeller intake is in the state of negative attack angle.
By exerting a total pressure pulse at the upstream inlet of the impeller, and meanwhile keeping other boundary conditions unchanged, the flow responses to inlet pressure pulse were calculated. The pressure pulse appears as a sudden drop. The variation amplitude is 5% of the initial total pressure, and the pulse width is 0.001 s. Defining the non-dimensional time t* as the physical time divided by the pulse width. The pulse occurs at t* = 0.1, reaches a peak at t* = 0.6, and ends at t* = 1.1. The time step size is 0.0001 s. Figures 7(a) and (b) show the response curves of mass flow rate and aerodynamic torque. It is obvious that the response curves of mass flow rate obtained from this code agree well with those from the software. There exists no distinct difference at any time instant. While for the results of aerodynamic torque, there exists a small

Description of computational model
In view of the complicate operating conditions of centrifugal compressors and the high costs of hardware construction, experiment design and equipment debugging, the calculation and analysis of aerodynamic damping for centrifugal impellers are carried out by numerical approaches in this work. As mentioned earlier, there are few computational studies on the aerodynamic damping  characteristics of centrifugal impellers, and the difficulty lies in the influence of blade-disk coupling vibration on the flow response.
The research object is a semi-open radial impeller of a single-stage centrifugal compressor with a wheel diameter of 280 mm, an outlet width of 16.8 mm and a tip clearance of 0.8 mm. The blade number is 16. For the convenience of multi-passage grid deformation, the full model of the impeller flow field is built in this work. The aerodynamic mesh consists of 80 H-type blocks with a total of 8,722,384 grid points and a minimum angle of 8.6°between two cell edges. The mesh is shown in Figure 8. It is obvious that the numerical modeling takes into account the leakage flow between the wheel disk and the casing. There are many studies existing which only modeled the blade passages when analyzing the aerodynamic forces of the impellers. The flow model without consideration of the leakage flow on the backside of the wheel disk will cause a significant deviation in the calculation of aerodynamic damping. The following discussed results are calculated by the developed code. The convection terms are discretized by the central scheme, and the turbulence calculation adopts the S-A model. The inlet boundary is applied with a total pressure of 101.3 kPa and a total temperature of 293.15 K. The rotational speed is 14,250 rpm. The leakage outflow through the backside gap of the wheel disk is set to 0.003 kg/s. By regulating the downstream back pressure, different operating points can be simulated. The unsteady flow field is initialized by the steady flow, and the number of iteration steps in the pseudo time domain is 30.
The modal aerodynamic damping ratios under two impeller modes are calculated. Before the flow calculations, the modal analysis of the impeller is carried out by the software Abaqus. The structural mesh has a total of 173,968 nodes, as shown in Figure 9. The 10-node tetrahedron element is adopted here. The material density is 2700 kg/m 3 , elastic modulus is 7.2 × 10 10 Pa, and Poisson's ratio is 0.33. An angular velocity of 1492.3 radian/s is applied on the impeller. The axial restraint is exerted on the shaft shoulder position, and an  initial interference is exerted on the shaft hole. Figure 10 shows two mass orthonormalized modal shapes of the extracted modes. From the modal shape contours it can be seen that the two modes both exhibit the characteristics of 3-nodal diameter vibration, which means that in the case of traveling wave vibration, the phase difference between any two neighboring sectors is 3π/8. For a centrifugal impeller, the vibration modes of different nodal diameters have different frequencies and wavenumbers. Furthermore, the vibration coupling between the blades and the wheel disk is very strong, so the impeller vibration must be analyzed as a whole, which is distinctly different from the vibration of axial blades.

Modal deformation of aerodynamic mesh
To calculate the modal aerodynamic damping ratio of the impeller, the modal deformation should be transferred from the structural mesh to the aerodynamic mesh. For the radial basis function method, in order to determine the interpolation coefficients in Equation (12), a dense coefficient matrix of 173,972 × 173,972 should be constructed, but the memory requirement cannot be met with PCs. By utilizing the compactly supported radial basis function method, the number of non-zero elements of the coefficient matrix is reduced to 17,051,442. Figure 11 shows the modal shape of mode 6 expressed by the aerodynamic surface mesh. The contour is almost identical to that in Figure 10(a), which reflects a good interpolation accuracy.
By adopting the compactly supported radial basis function, the deformation of aerodynamic mesh is realized. Taking mode 6 as an example, the maximum deformation of the blade tip leading edge is set to 1 mm. The deformed mesh is shown in Figure 12. The dark shadow entity is the undeformed impeller, the green transparent entity is the deformed impeller, and the green lines are the edge curves of the deformed tip clearance mesh. In the gap region from the blade tip to the casing, the grid deformation varies from the maximum to zero. Although there exists some grid distortion due to the deformation, the minimum angle of the deformed mesh is 8.1°, which is not much different from that of the original mesh. This result indicates that the conformality of the grid deformation algorithm is fairly good. When applying Equation (12) to calculate the displacement of each aerodynamic grid point, the Euclidean distance x−x i from that point to the structural surface node x i should be calculated. For the computational model in this work, the Euclidean distance needs to be calculated 8,722,384 × 173,968 times. If the direct method is adopted, the amount of calculation will be very huge. Whereas for the compactly supported radial basis function method, the Euclidean distance is only calculated within the neighborhood R of each aerodynamic grid point. To this end the alternating digital tree technique described in section 2.6 can also be used here to realize the fast searching of structural nodes. The calculation results show that, on a dual-core and four-thread PC, if the direct method is utilized, the computational time for one grid deformation is 44 min 10 s, and if the alternating digital tree algorithm is utilized, the computational time reduces to only 3 min 20 s.

Results of aerodynamic damping
For the simulation of the flow field with moving boundaries, the grid deformation, the wall distance calculation and the solution of flow equations are in turn performed at each physical time step. To complete the distance calculation for the grid scale of 8,722,384 points, the direct method will be very time consuming. On a dual-core and four-thread PC, the elapsed time is 17 min 10 s. If the alternating digital tree algorithm is used, then the elapsed time is only 36 s. The iterative solution of flow equations at each physical time step needs 12 min 45 s, and thus the fast calculation of wall distance saves more than half of the time cost, which is very beneficial to realize the aerodynamic damping calculation of centrifugal impellers in this work.
First the modal aerodynamic damping ratio is calculated for mode 6. Each vibration cycle is divided into 48 time steps, and the maximum deformation of the blade tip is set to 1 mm. Figure 13 shows the timevarying curves of generalized aerodynamic force at three operating points with different back pressures. Here the generalized force is defined as where f and are respectively the aerodynamic load vector and the modal shape vector, and S denotes the aerodynamic surface. Usually there exist significant differences in the flow fields of different operating points, but as shown in Figure 13, except for the two cycles after the start of calculation the time-varying generalized aerodynamic forces at different operating points are fairly consistent with each other. The generalized force is a scalar, which represents the work done by the force along the modal displacement. The calculated results show that the abilities of the aerodynamic force to resist the blade deformation are not quite different for the studied operating points. The aerodynamic damping work in one period is defined as where q(t) is the generalized vibration displacement on the modal shape, andf (t)q(t) is the aerodynamic power. Figure 14 shows the aerodynamic power curves of the impeller under different back pressure conditions. As with the generalized aerodynamic force, after the calculation of several cycles, the overall differences between the aerodynamic power curves of different operating points are small. Different from the time varying of generalized aerodynamic force, the aerodynamic power in the wide time domain is almost always negative, which indicates that the aerodynamic force always hinders the vibration of the impeller. Combined with the definition of aerodynamic power, it can be deduced that there exists phase difference between the generalized aerodynamic force and the generalized vibration velocity. Namely the flow responses do not synchronize with, but lag behind the impeller vibration. The numerical simulations are also carried out for the cases when the maximum deformations of the blade tip are 0.5 mm and 0.2 mm, and the obtained results of aerodynamic force and power show similar characteristics. By introducing the definition of viscous damping, the modal aerodynamic damping ratio can be calculated by where denotes the modal frequency, and A max denotes the generalized vibration amplitude. Figure 15 shows the trend of modal aerodynamic damping ratio of mode 6 with the variation of mass flow rate at different operating points. The mass flow rates corresponding to different back pressures are shown in Table 1. When the back pressure is 137.8 kPa, the operating point is near stall. As shown in Figure 15, the aerodynamic damping ratio does not increase or decrease monotonically with the variation of mass flow rate. The values are all positive and differ little. The relative difference between the maximum and minimum values is less than 2.1%. Furthermore, the differences of aerodynamic damping ratio are very small under different vibration amplitudes. According to the linearized small-perturbation theory, the unsteady aerodynamic force caused by the structural vibration is  approximately a linear function of the vibration amplitude, while the aerodynamic work is approximately proportional to the square of vibration amplitude. According to the definition of Equation (19), the aerodynamic damping ratio should be independent of the vibration amplitude. The calculated results are in good agreement with this conclusion. Taking the operating point with a back pressure of 137 kPa as an example, when the maximum vibration amplitude of the blade tip is 1 mm, the aerodynamic work per cycle is −0.0349 J, and when the maximum vibration amplitude of the blade tip is 0.5 mm, the aerodynamic work per cycle is −0.00873 J. In addition, the results without consideration of the leakage flow are also illustrated in Figure 15. It is obvious that the leakage flow has little effect on the aerodynamic damping ratio. This is not difficult to understand, because mode 6 is dominated by the blade vibration, and the vibration of the wheel disk is weak. Therefore the aerodynamic work done by the leakage flow is small. Figure 16 shows the distribution of average aerodynamic power density for mode 6 during the 1/4 period from the maximum deformation to zero deformation. The blade surfaces seen in Figure 16 are all suction surfaces. The aerodynamic power density on the pressure side has a similar distribution but with an opposite sign. Without consideration of the sign of aerodynamic power density, its distribution characteristics will be fairly similar to the modal shape of the impeller, which indicates  that the aerodynamic power density is mainly affected by the local vibration velocity.
Then the modal aerodynamic damping ratio and the aerodynamic power density for mode 23 are calculated. Each vibration cycle is divided into 52 time steps, and the maximum deformations of the wheel disk are respectively set to 0.1 and 0.2 mm. Figure 17 shows the trend of modal aerodynamic damping ratio of mode 23 with the variation of mass flow rate at different operating points. Just like with mode 6, the influence of vibrating amplitude on the aerodynamic damping is also weak. Furthermore, it is obvious that the aerodynamic damping ratio reduces gradually as the mass flow rate decreases, which indicates that as the operating point shifts toward the stall, the aeroelastic stability of the impeller becomes worse.
The difference between the maximum and minimum values of aerodynamic damping ratio is 11.65%. Compared with mode 6, the aerodynamic damping ratio of mode 23 is significantly affected by the change of operating points. Evidently this mode is dominated by the vibration of the wheel disk, which is different from the vibration pattern of mode 6. The aerodynamic work done by the leakage flow on the backside of the wheel disk inevitably affects the aerodynamic damping characteristics of the impeller. However, the geometric detail of the leakage passage is usually complex, and furthermore the leakage flow exerts little influence on the main stream, so most numerical analyses do not include the leakage model. The results in Figure 17 show that the aerodynamic damping ratio reduces a lot without consideration of the leakage flow on the backside of the wheel disk, which implies the necessity to include the leakage model for the disk-dominant mode. Figure 18 shows the distribution of average aerodynamic power density for mode 23 during the 1/4 period from the maximum deformation to zero deformation. It is obvious that the contour also approximately reflects the characteristics of modal shape. Furthermore, because the wheel disk vibration is mainly axial, the aerodynamic power density of the pressure side and the suction side near the blade exit is close to zero, and the aerodynamic damping ratio of mode 23 mainly reflects the ability of aerodynamic force to resist the axial deformation of the impeller.
The earlier discussed is around the modal aerodynamic damping ratios of two impeller modes under the standing wave vibrations. In what follows, the aerodynamic responses and damping characteristics of the   impeller under the forward and backward traveling wave vibrations will be calculated and analyzed. In reality, the traveling wave vibrations of centrifugal impellers may be excited by a variety of excitation sources. Observed in the rotating coordinate system, the forward traveling wave propagates along the rotating direction, while the backward traveling wave propagates back the rotating direction. Figure 19 shows the instantaneous axial deformation of the impeller under the backward traveling wave vibration of mode 23. It is obvious that when the traveling wave propagates a circle, the flow boundaries of the impeller vary three cycles.
The calculations are carried out for both modes at the operating point with a back pressure of 135 kPa. The maximum deformation of the blade tip is set to 1 mm for mode 6, and the maximum deformation of the wheel disk is set to 0.2 mm for mode 23. After the unsteady flow simulation begins from a steady flow, the flow field undergoes a period of transient response and gradually converges to a quasi steady state. The forward and backward traveling    wave vibrations lead to different flow responses. Figure 20 shows the convergence curves of the aerodynamic thrust of the impeller during the calculation process. It can be seen that the time histories of aerodynamic thrust for the forward and backward traveling wave vibrations are different. For the disk-dominant mode, the aerodynamic thrusts converge to different values, and the aerodynamic thrust caused by the backward traveling wave vibration is greater. But for the blade-dominant mode, the final difference of aerodynamic thrust is little. Figure 21 shows the time histories of the aerodynamic power of the impeller during the calculation process. It can be seen that the aerodynamic power under the backward traveling wave vibration is larger for both modes, which indicates that the backward traveling wave vibration has larger aerodynamic damping. According to the definition of Equation (19), for the forward traveling wave vibration the aerodynamic damping ratios of mode 6 and mode 23 are respectively 0.00247 and 0.00408, while for the backward traveling wave vibration the values are respectively 0.00697 and 0.00653.
For each mode, the distribution of aerodynamic power density at any time under the traveling wave vibration is similar to that in Figure 16 or Figure 18, which shows the characteristics of modal shape. However, by the averaging of the aerodynamic power density in a vibration cycle, the average power density will present a cyclic symmetric distribution. Figure 22 shows the distribution of average aerodynamic power density of mode 23 under the backward traveling wave vibration. It is obvious that the distribution of aerodynamic power density for each sector is the same. The average aerodynamic power density reflects the aerodynamic accumulative work per unit area in a cycle. It is shown in Figure 22 that the average aerodynamic power density is negative in most areas of the impeller surface, and the average aerodynamic power density on the backside of the wheel disk is negative everywhere, which indicates that the leakage flow exerts a positive damping effect on the impeller vibration. By setting the same range of contour legend, the average power density under the forward traveling wave vibration is shown in Figure 23. In comparison with the backward traveling wave vibration, the absolute value of average power density is smaller at the same position, which indicates smaller accumulative work. Furthermore, it can be seen from the front views that the distribution of average power density on the windward side is more axisymmetric for the backward traveling wave. This phenomenon indicates that the flow disturbance along the circumference is more homogeneous for the backward traveling wave. Figure 24 shows the distributions of average aerodynamic power density of mode 6 under the backward and forward traveling wave vibrations. As with the results in Figures 22 and 23, the average power density presents a cyclic symmetric distribution. Under the forward traveling wave vibration, there exists some area on the blade suction side where the average power density is positive. While under the backward traveling wave vibration, the average power density is negative almost everywhere on the blade surface.

Conclusions
By adopting the compactly supported radial basis function and the alternating digital tree technique, the consumptions of CPU time and memory space for grid deformation and wall distance are reduced greatly. By using the developed numerical code, efficient computations of aerodynamic damping of centrifugal impellers are first realized in this work. The analysis of aerodynamic damping characteristics of a radial impeller indicates the following.
The adoption of the compactly supported radial basis function greatly reduces the memory requirement of the classic radial basis function method, and in the meantime possesses good numerical accuracy. By applying the alternating digital tree technique, the time cost of the compactly supported radial basis function method is also greatly reduced. The actual tests show that the alternating digital tree technique can improve the computational efficiencies of grid deformation and wall distance by at least one order of magnitude.
Under the studied blade-dominant mode and diskdominant mode the modal aerodynamic damping ratios are both positive. The modal aerodynamic damping ratio of the impeller is little affected by the vibration amplitude, and the calculated results show that the aerodynamic work is proportional to the square of vibration amplitude under the studied vibration conditions. The modeling of the leakage flow has an obvious effect on the aerodynamic damping for the disk-dominant mode.
For the blade-dominant mode, the ability of the aerodynamic force to resist the blade deformation varies little at the studied operating points, so the variation of aerodynamic damping ratio is small. While for the diskdominant mode, the aerodynamic damping ratio of the impeller shows a gradual decline trend as the operating point shifts toward the stall.
In the case of traveling wave vibration, the aerodynamic force still performs negative work on the whole impeller. The forward and backward traveling wave vibrations lead to different flow responses and aerodynamic dampings, and the aerodynamic damping caused by the backward traveling wave vibration is larger.
The aerodynamic damping characteristics of centrifugal impellers may be affected by many factors from aerodynamic and structural aspects. Only the operating points on one speed line is considered in this work. The variation trends of aerodynamic damping characteristics under different speed lines can be studied in future work. Another problem requiring investigation is to what extent will the intake parameters such as pressure and temperature affect aerodynamic damping. The intake pressure is thought to be an important factor affecting aerodynamic damping. Considering from structural aspects, the investigated impeller in this work is unshrouded. Whether the aerodynamic damping characteristic of a shrouded impeller is different from that of an unshrouded impeller is also worth studying. Besides the earlier mentioned physical issues, there are still algorithmic difficulties to be solved in the next work. The unsteady CFD calculation is very time-consuming for the structure of the centrifugal impeller, and therefore more efficient flow simulation solutions should be sought.

Disclosure statement
No potential conflict of interest was reported by the authors Funding This work is supported by the National Natural Science Foundation for Young Scholars of China [grant no. 51405130].