Static aeroelasticity analysis of a rotor blade using a Gauss-Seidel fluid-structure interaction method

The present study introduces a Gauss-Seidel fluid-structure interaction (FSI) method including the flow solver, structural statics solver and a fast data transfer technique, for the research of structural deformation and flow field variation of rotor blades under the combined influence of steady aerodynamic and centrifugal forces. The FSI method is illustrated and validated by the static aeroelasticity analysis of a transonic compressor rotor blade, NASA Rotor 37. An improved local interpolation with data reduction (LIWDR) algorithm is introduced for fast data transfer on the fluid-solid interface of blade. The results of FSI calculation of NASA Rotor 37 show that when compared with the radial basis function (RBF) based interpolation algorithm, LIWDR meets the interpolation accuracy requirements, while the calculation cost can be greatly improved. The data transmission time is only about 1% of that of RBF. Moreover, the iteration step of steady flow computation within one single FSI has little impact on the converged aerodynamic and structural results. The aerodynamic load-caused deformation accounts for nearly 50% of the total. The effects of blade deformation on the variations of aerodynamic performance are given, demonstrating that when static aeroelasticity is taken into account, the choke mass flow rate increases and the peak adiabatic efficiency slightly decreases. The impact mechanisms on performance variations are presented in detail.


Introduction
For modern high-performance and high-load fan/compressor rotor blades, especially for the ones with high aspect ratios, blade deformations resulted by the compound actions of high aerodynamic load and centrifugal forces are common.The flow in the passage of the deformed blade subsequently deviates from the design, thus greatly affecting both the aerodynamic and structural performance.The geometric deformation of the blade under steady aerodynamic load and centrifugal force is defined as a static aeroelasticity problem [1].Predicting blade deformation and aerodynamic performance accurately under different operating conditions is of great significance for blade aerodynamic design.
Fluid-structure interaction (FSI) is a typical multidisciplinary analysis (MDA) problem.A widely used strategy for solving MDA is the nonlinear block Gauss-Seidel (NLBGS) approach [2], in which the fluid and structure fields are solved sequentially with mutual iteration until convergence.Compared with weak coupling methods, NLBGS is more in line with the physical mechanism and has been extensively studied in external flows [3][4][5][6][7][8].As for the internal flow, some developments in the numerical modeling of FSI are described as follows.Doi [9] developed a FSI calculation platform based on an unsteady three-dimensional Reynolds-averaged Navier-Stokes (RANS) solver and a structural solver, by which the aeroelastic response of Rotor 67 under various operating conditions was successfully predicted.Dowell [10] reduced the order of unsteady aerodynamic load and successfully reduced the computational cost of FSI.Liu [11][12][13] presented a computational method for flutter prediction of turbomachinery cascades, and the computational results for a cascade were presented and compared with those obtained by the conventional energy method and with experimental and numerical data by other researches.Carstens [14] analyzed the flutter behavior of turbomachinery blading in the time domain and introduced the Newmark method into the structural analysis for time discretization.Rzadkowski [15,16] simulated the aeroelastic behavior of an oscillating turbine blade row based on a full 3D moving grid.Zheng [17][18][19] examined the deformation of transonic fan blades under aerodynamic load and how it affected aerodynamic performance using the two-way FSI.Zhou [20] predicted and analyzed the aeroelasticity of transonic fan rotor blades.
FSI involves the transmission of interface data between fluid field and structure field.The commonly used data transfer methods are generally divided into global interpolation and local interpolation.For the global interpolation method, radial basis functions (RBFs) [21,22] are often used, which allow the fluid field and structure field to adopt any grid form and are widely used to solve the data transfer problem in FSI.However, RBF is computationally expensive and has poor interpolation precision on the interface.Thus it greatly depends on the design of the interpolation matrix [21,23].The local interpolation method [24,25] calculates the corresponding unknown data of the nodes/elements on the fluid-solid interface from some known nodes/elements on the interface.However, local interpolation in the three-dimensional coordinate system involves a large number of search procedures and repeatedly solving the mapping equations.Moreover, due to the different grid resolution of fluid and solid fields, usually it is unable to detect the master element on the interface, leading to an increased transmission error.
Currently, commercial software is often utilized to calculate the two-way FSI problems by solving unsteady governing flow equations and structural dynamics equations [26][27][28].Such FSI studies based on unsteady flow computations require a large number of computational resources and are difficult to be applied in multidisciplinary blade design.The present study focuses on the static aeroelasticity problem of rotor blades and investigates the impact mechanisms of blade deformation on the variations of aerodynamic performance.The organization of the paper is as follows.First, the static aeroelasticity analysis platform based on Gauss-Seidel fluid-structure interaction (GS-FSI) method is introduced, in which the RANS flow solver, structural statics solver and a fast data transfer technique are included.Second, a local interpolation with data reduction (LIWRD) method is proposed for fast data transfer on the fluid-solid interface.The improvements of data transfer efficiency and accuracy of LIWDR are verified and validated through the FSI calculation on a transonic compressor rotor blade, NASA Rotor 37. Finally, the improved GS-FSI method is used in the static aeroelasticity analysis of Rotor 37. The results are given in detail and the influence of blade deformation on aerodynamic performance variations is illustrated.

Gauss-Seidel fluid-structure interaction method
In this section, the Gauss-Seidel method is introduced to investigate the fluid-structure interaction phenomena of rotor blades.GS-FSI allows a simple combination of the existing flow solver and structural statics solver.Meanwhile, the technique of data transfer on the fluid-solid interface is also required.Figure 1 presents the procedure of GS-FSI method.As shown in the figure, the fluid equations are firstly solved on the original blade, and then the aerodynamic load is transferred to the solid field.By solving the structural statics equations, a new deformed blade can be obtained, which is then used in the fluid field calculation again.The sequential steps are repeated until FSI convergence is reached.In the study, the FSI calculation converges when where ε CFD and ε CSM represent errors of the normalized flow and structural solutions, respectively between K and K − 1 FSI steps, where the solutions at the first FSI step are regarded as the references; U CFD and U CSM,max are the aerodynamic parameters and the maximum blade deformation. (1)

Flow solver
The transonic rotor blade NASA Rotor 37 [29,30] will be used in the following static aeroelasticity analysis.NASA Rotor 37 was designed and tested originally by Reid and Moore [29], and the specifications are given in Table 1.Flow simulation of this rotor blade is carried out using an in-house program Turbo90 [31,32], which solves the steady-state N-S equations and the Spalart-Allmaras (SA) turbulence equation [33].
The JST scheme is used for spatial discretization.A fifth-order Runge-Kutta method is used for time stepping.In addition, multi-grid method, local time step, implicit residual smoothing are used to accelerate the convergence.
A single-block H-type grid is generated for the computational fluid dynamics (CFD) calculation of Rotor 37, and grid-independent study is performed sequentially in the pitchwise (node numbers: 33/41/49/57/65) and spanwise (node numbers: 37/49/61/73/85) directions (Fig. 2).The total pressure ratio π and mass flow rate MFR of all the grids are shown in Fig. 3. Finally, the grid in the spanwise grid-independent study with 161, 57 and 49 nodes in the streamwise, pitchwise and spanwise directions, respectively, is used in the following study, which is shown in Fig. 4.
Figure 5 shows the operation characteristics of π and adiabatic efficiency η at 100% design speed.Generally, the computed total pressure ratio is higher than the experimental one, while the computed adiabatic efficiency is higher.Although obvious deviations between CFD and experimental results can be found, the CFD results are quite similar to those of others [34][35][36].

Structural statics solver
An in-house finite element analysis (FEA) program is used to solve the structural statics equations.The material of the blade is titanium-alloy, with a density of 4500 kg/m 3 , a Young's modulus of 116 GPa, and a Poisson's ratio of 0.32.The grid-independent study of computational structural mechanics (CSM) is performed and the results are given in Fig. 6.Five gradually refined grids (element numbers: 9000/13000/19000/29000/38000/ 42000) are used to calculate the maximum blade deformations in the three coordinate directions, which are named by dx, dy and dz, respectively in the figure.Notice that only the centrifugal force is considered.The fourth grid with a total element number of 38000 is ultimately selected for further research, which is shown in Fig. 7.
In addition, ANSYS Workbench is used for FEA on the blade under centrifugal force (CF) only and both centrifugal force and aerodynamic load (CF+AL), separately.The results are compared with those obtained by the in-house program.The relative deviations of the maximum deformations obtained by these two solvers are shown in Table 2.When both centrifugal force and aerodynamic load are considered, the relative deviation is about 2.4%, indicating that the FEA program used in this study has a high degree of reliability.
3 Data transfer method on fluid-solid interface

Principles of LIWDR
GS-FSI firstly transfers the aerodynamic load from the fluid field to the solid field, and then feeds back the displacement from the solid field to the fluid field.Taking the displacement transfer as an example, the global interpolation method usually requires the construction of an interpolation matrix H fs (N f × N s ): where g s is the known data on the CSM grid; g f is the data needing to be determined on the CFD grid; N f and N s are the grid numbers of the fluid field and the solid field, respec- tively on the fluid-solid interface.H fs can be calculated using RBF [21,22].
where x and x i denote the coordinates of the unknown nodes and the known nodes, respectively.φ is the chosen RBF, α i is the coefficient to be solved, and p(x) is the pol- ynomial function.As previously stated, the accuracy of data transfer is significantly impacted by the construction of H fs .Besides, the size of H fs is large, resulting in high computation cost of interpolation.
For general three-dimensional local interpolation, the nodes on the CFD/CSM grid must be projected onto the plane where the elements on the CSM/CFD grid are (2)  placed, with the foot of the perpendicular as the mapping point.When the mapping point is inside the element, the element is called the master element.By using the data on the master element, the data on the unknown node can be interpolated.The search for master elements is generally based on the alternating digital tree algorithm [37].
For each element on the fluid-solid interface, its geometry can be described using shape functions through a transformation.As shown in Fig. 8, considering the twodimensional four-node quadrilateral isoparametric transformation, given the local coordinates of a point S within the element and the global coordinates of the four boundary nodes of the element, the global coordinates of point S can be calculated by Eq. ( 4).
(4)  where R s is the data at point S, and R i is the data on the i-th boundary node.
Three-dimensional local interpolation involves a large number of search procedures and repeatedly solving the mapping equations.Meanwhile, it may be unable to detect the master element on the interface since the fluid and structural fields have different grid resolutions.In such situations, the nearest-neighbor interpolation method is usually employed, but will lead to an increased transmission error [38].Therefore, an improved fluid-solid data transfer technique based on two-dimensional local interpolation method is proposed.The procedure of this method is depicted in Fig. 9. (

1) Coordinate Transformation
In this step, the coordinates in the three-dimensional X-Y-Z Cartesian system are firstly transferred to the x-rθ-r cylinder system.Then the CFD and CSM blades are unripped at the trailing edge (TE) and stretched; and meanwhile, the coordinates of x and r are normalized by the blade chord at hub and the blade height at TE.It is necessary to ensure that both blades are unripped at the same TE.In reality, after determining the TE of the CFD blade, a B-spline function is used to interpolate the TE of the CSM blade.
(5)  (2) Search of Master Element Possible master element search: For each node K of either the CFD or CSM grid on the fluid-solid interface, find out the possible master element of the K-th node.The global coordinates of the element boundary nodes satisfy: Exact master element search: Then the two-dimensional inverse isoparametric transformation is used to calculate the local coordinates of node K on the possible master element.When the local coordinates (ζ k , η k ) of node K on a particular element are all within the range of [−1, 1] , the found element is then regarded as the master element of node K.
The mapping relationship from any quadrilateral to the master element is called the inverse isoparametric transformation, as shown in Fig. 8. Taking a four-node quadrilateral as an example, the expression for the inverse isoparametric transformation can be derived based on Eq. ( 4).
where Equations ( 7) and ( 8) can be solved by the Newton iterative method.
(3) Data Interpolation Based on the data of the four boundary nodes of the master element, such as load and displacements, we can interpolate the data on node K using Eq. ( 5).
The introduced data transfer method uses spatial coordinate transformation to reduce the data size and is the so-called local interpolation with data reduction (LIWDR).(6)

Fig. 8 Isoparametric transformation and inverse isoparametric transformation of two-dimensional four-node quadrilateral
Compared with the three-dimensional local interpolation method with a serious risk of not being able to detect the master element on the interface, LIWDR can find the master element fast for each node on the interface.Furthermore, by LIWDR, the computational cost of searching for the master element can be decreased.To verify the efficiency and accuracy of data transfer by LIWDR, the traditional RBF global interpolation method without any improved algorithm will also be used in the following subsection.

Verifications and validations
Rotor 37 is a transonic rotor blade and the distribution of static pressure on the blade near the shock wave is strongly nonlinear.In the present study, using the global interpolation method, the interpolated pressure distributions on the CSM grid deviate far away from the CFD ones.Therefore, only the interpolated pressure distributions using LIWDR are given to verify the data transfer accuracy of LIWDR in this section.LIWDR is used to interpolate the pressure from the CFD grid to the CSM grid at the operation condition near peak efficiency (NP) of Rotor 37. Figure 10 compares the pressure distribution contours on the fluid-solid interface.It can be found that the pressure distributions on both the CSM grid and the CFD grid are quite close, even near the shock wave, demonstrating that the data transfer accuracy of LIWDR is acceptable.
Besides pressure transfer from CFD to CSM, the performance including efficiency and accuracy of displacement transfer from CSM to CFD by LIWDR also requires to be verified.Figure 11 shows two different CSM grids.The one in Fig. 11(a) is used for FEA and the one in Fig. 11(b) is the "background grid" uniformly extracted from Fig. 11(a).First, the displacements of the nodes in Fig. 11(a) are calculated by FEA and used as the reference solutions.And meanwhile, the displacements of the nodes in Fig. 11(b) can be determined regarding as the "background solutions".Then, the Fig. 9 Local interpolation with data reduction method displacements of the nodes in Fig. 11(a) are interpolated using RBF and LIWDR methods separately from the "background solutions".By comparing with the reference solutions, the interpolation error of the two methods can be calculated.
Table 3 shows the computer time for a single data transfer cycle.Compared with RBF, LIWDR has superior data transfer efficiency because the expensive matrix inversion is not necessary for LIWDR.The relative interpolation error δ is defined as: where D is the interpolated displacement and L is the reference solution.In this study, δ is calculated only at the nodes with relative displacements higher than 1% of the maximum.Table 4 and Table 5 illustrate the maximum and average δ in each dimen- sion, respectively obtained by LIWDR and RBF, where δx , δy and δz represent the rela- tive interpolation errors in the corresponding directions.According to Table 4, it can be observed that LIWDR has a significantly lower maximum relative error compared to RBF.Besides, although the average relative errors obtained by LIWDR and RBF are small, the ones obtained by LIWDR are lower than those of RBF in the Y and Z directions.
The maximum deformation of a rotor blade typically occurs at the blade tip.It is necessary to ensure that the interpolation accuracy at the blade tip meets the requirements.Figure 12 shows the comparisons of data transfer error distributions at the blade tip.Notice that the blade is unripped at TE and stretched in the streamwise direction.As shown in Fig. 12, the interpolation error of LIWDR in each direction is significantly lower than that of RBF.
The displacement is transferred from the CSM grid to the CFD grid by LIWDR at NP condition of Rotor 37. Taking the X direction displacement transfer as an example, Fig. 13 compares the displacement contours on the two-dimensional fluid-solid interface, where x represents the displacement in the X direction, and R/R 0 represents the normalized blade height.The displacement contour on the fluid interface is determined by LIWDR from the one on the solid interface.It can be seen that the LIWDR data transfer method introduced in this paper can meet the accuracy requirements of fluidsolid interface interpolation.
Figure 14 shows the original blade, the deformed blade of CSM, and the deformed blade of CFD from displacement transfer, where Fig. 14

Static aeroelasticity analysis
In this section, the GS-FSI method is used to predict the static aeroelasticity performance of Rotor 37.As shown in Fig. 1, the CFD and CSM governing equations are sequentially solved.After certain steps of RANS iterations, the aerodynamic loads are transferred to the solid field.Then the structural statics equations are solved, the solutions of which, i.e., blade deformations are fed back to the fluid field.Until both the flow and solid fields converge, blade deformation and aerodynamic performance in the equilibrium state can be obtained, regarding as the converged FSI solutions.
To investigate the effects of the RANS iteration step on the converged FSI solutions, the structural statics equations are solved after 100, 400, and 800 RANS iteration steps, which are referred to as F 100 , F 400 and F 800 , respectively in the following.The convergence histories of maximum displacements in three dimensionalities, dx, dy and dz and the aerodynamic performance parameters including total pressure ratio, π and adiabatic efficiency, η are shown in Fig. 15.It can be found that as the RANS  iteration step increases, the necessary FSI iteration cycle decreases.After 10 iterations, F 100 reaches an equilibrium state, while F 400 and F 800 reach an equilibrium state after 5 iterations.Table 6 gives the converged FSI solutions with different RANS iteration steps.It can be seen from the table that the RANS iteration step within a single FSI has little impact on the converged FSI solutions.Since the time required for solving the structural statics equations is far less than that for solving RANS equations, the aerodynamic load is transferred to the solid field when the flow field is still not fully converged after only 100 RANS iteration steps, by which the computational cost of the FSI calculation can thus be significantly reduced.
Table 7 shows the maximum deformations considering only the centrifugal force and static aeroelasticity at NP condition and near stall (NS) condition, separately.It can be seen that the deformation of Rotor 37 due to static aeroelasticity is significantly increased compared with that considering only the centrifugal force.As the total pressure ratio of Rotor 37 increases from NP condition to NS condition, the deformation contributed by the aerodynamic load increases.Ignoring the effects of aerodynamic load on blade deformation will reduce the accuracy of the static aeroelasticity calculation.Figure 16 shows the spanwise distributions of π and η at the outlet under NP condi- tion.The total pressure ratio increment of the deformed blade, identified by NP-FSI in the figures, is obvious and it is intensified along the span over 30% of the blade height, while the adiabatic efficiency is slightly decreased on most spans.
Figures 17 and 18 show the blade profiles and the corresponding pressure distributions on the blade surface at 30%, 70%, and 90% of the blade height, respectively.It can be found that the deformation of the FSI blade increases in the spanwise direction.Evident variations of the blade profile on 70% and 90% spans can be found, while the profile of the FSI blade on the 30% span is almost a duplicate of the original profile.Generally, such deformations illustrated in Fig. 17 result in an increased incidence flow angle, which consequently increases the blade loading.In Fig. 18, the loading increases for the FSI blade, especially on the front portions.Moreover, due to the most increased incidence angle, the loading increment on the 90% span is higher than those on 30% and 70% spans.Moreover, as shown in Fig. 18, on the suction side of Rotor 37 there exists a shock wave on the middle portion of the blade in the axial direction.Compared with the original blade, the shock wave of the FSI blade slightly moves upstream and it becomes stronger.In such situations, the total pressure ratio and flow loss increase and the adiabatic efficiency decreases, which have already been illustrated in Fig. 16.
Finally, the effects of static aeroelasticity on the overall performance of Rotor 37 are investigated.Figure 19 presents the operation characteristics of the original and FSI blades.It is obvious that the choke mass flow rate of the FSI blade significantly increases (0.604%).That is because the deformation, as shown in Fig. 17, decreases the stagger angle of the original blade and meanwhile increases the throughflow area.Furthermore, the total pressure ratio and adiabatic efficiency exhibit a slight increase and decrease, respectively, in the whole operation range.Moreover, as the operation condition moves from NP to NS, more total pressure ratio increments and adiabatic efficiency decrements are resulted by blade deformation.

Conclusions
This paper presents a fast and accurate data transfer method, LIWDR, between fluid and solid interfaces of rotor blades.The method is integrated in the Gauss-Seidel sequential fluid-structure interaction calculation platform and then applied   to the static aeroelasticity study of a transonic fan rotor blade.The accuracy and efficiency of LIWDR in data transfer are verified and validated by comparing with the RBF method.The effectiveness of the improved GS-FSI method is illustrated through the static aeroelasticity study of NASA Rotor 37. The effects of blade deformation on aerodynamic performance are analyzed as well.The main conclusions are as follows: (1) The computer time needed by LIWDR is only about 1% of that of the global interpolation method, RBF, while the relative errors of displacement transfer are close for these two methods.For the GS-FSI calculation, the RANS iteration step within a single FSI has little impact on the converged FSI solutions including the blade deformation and aerodynamic performance parameters.

Fig. 10
Fig. 10 Pressure contours on: a pressure side; b suction side.Left: pressure on the CSM grid; right: pressure on the CFD grid

Fig. 12
Fig. 12 Comparison of data transfer errors at the blade tip: a X direction; b Y direction; c Z direction

Fig. 13
Fig. 13 Deformation of blade in the X direction: a CSM grid; b CFD grid

Fig. 15
Fig. 15 Convergence histories of FSI calculation: a deformation; b aerodynamic performance

Table 1
Specifications of NASA Rotor 37

Table 2
Comparisons of deformation calculated by FEA program and ANSYS Workbench

Table 3
Computer time of data transfer by LIWDR and RBF

Table 4
Maximum relative error in each direction

Table 5
Average relative error in each direction

Table 6
Converged FSI solutions

Table 7
Maximum deformation under different operation conditions