Engineering Application of an Identification Method to Shock- Induced Vortex Stability in the Transonic Axial Fan Rotor

In the present paper, the steady RANS (Reynolds-Averaged Navier-Stokes) simulations based on our independently developed CFD (Computational Fluid Dynamics) solver NUAA-Turbo 2.0, are carried out to investigate the shock wave/tip leakage vortex (SW/TLV) interaction in two representative transonic axial fan rotors, NASA Rotor 67 and NASA Rotor 37. The intent of this study is mainly to verify if an identification method derived from relevant theories is suitable for shock-induced vortex stability in the real engineering environment. As the additional findings, a universal tip vortex model is established and the characteristics of vortex breakdown or not are also summarized under different load levels. To ensure the prediction accuracy of all numerical methods selected in this research, detailed comparisons are made between computational and experimental results before flow analysis. The excellent agreement between the both indicates that the current code is capable of capturing the dominant secondary flow structures and aerodynamic phenomenon, especially the vortex system in tip region and SW/TLV interaction. It is found that three vortical structures such as tip leakage vortex (TLV), shock-induced vortex (SIV), tip separation vortex (TSV) in addition the tip leakage vortex-induced vortex (TLV-IV, which only occurs when the TLV strength increases to a certain extent) frequently exist near the blade tip and then abstracted as a tip vortex model. A stable TLV after passing through the passage shock is commonly characterized by tight rolling-up, slow deceleration and slight expansion. Conversely, the vortex behaves in a breakdown state. The final verification results show that the above two vortex states can be satisfactorily detected by the theoretical discriminant introduced in this work.


Introduction
The secondary flow structures with their related aerodynamic phenomena distributed in the tip region, such as various vortices and their interactions with passage shock, have significant influences on further improvement in total pressure ratio, adiabatic efficiency and stable operating range [1,2] for a highly-loaded transonic axial compressor. For the aspect of establishing vortex model, some efforts have been made in terms of subsonic compressor rotor. Based on the results of their detailed measurements, Inoue and Kuroumaru [3] firstly proposed a vortex model outlining the basic structure of vortex system in a axial-flow rotating blade row, which includes a leakage vortex, scraping vortex, horseshoe vortex, trailing vortex, etc. Yu et al. [4] also utilized a more sophisticated measurement technique called SPIV (stereoscopic particle image velocimetry) to construct the three-dimensional complex flow structures near the blade tip, they are clearly depicted with schematic model incorporating the tip leakage vortex accompanied by its induced vortex, corner vortex, etc. Recently, based on the numerical investigations in a mixed flow pump, Liu and Tan [5,6] divided the TLV into two main types. One is called the primary tip leakage vortex (PTLV), which presents a continuous tubular pattern; the other is the secondary tip leakage vortex (STLV), which distributes along the suction surface with a discontinuous strip pattern. The unstable PTLV was further classified as two parts [7], oscillating PTLV-A and shedding PTLV-B. However, they can only provide a limited guidance on understanding of typical vortices in a transonic axial compressor rotor due to lack of shock effect. In addition, few researchers has payed attention to studying its vortex model especially near the blade tip, and thus there is still no a relatively complete and acceptable tip vortex model involving the vortex types and distributions. Wisler [8] attempted to get insight into specific vortex model for the transonic case, they suggested that two main vortices existed in the tip passage, namely the TLV originating from the leading edge and the tip corner vortex (TCV) generated by interaction of trailing edge separation and radial migration flow. Within the scope of literatures known by the author, a relatively perfect vortex system in the whole passage is numerically investigated by Wang et al. [9], their noteworthy achievement is the discovery of a transversal vortex structure induced by passage shock in the tip region (called shock-induced vortex). In order to be more convincing, a general tip vortex model in present work is obtained through extracting the same vortical structures from two typical transonic fan rotors mentioned earlier.
As known, the SW/TLV interaction is an inherent aerodynamic phenomenon in a transonic axial fan rotor, the extent of its influence on the overall performance depends on the interaction strength and whether the vortex bursts. Therefore, it is desirable to accurately identify the flow state of TLV after passing through the passage shock. Since the vortex breakdown was initially observed in engineering practice, many researchers [10][11][12][13][14][15] had widely performed the experimental basic studies with regard to the normal (or oblique) shock wave/supersonic streamwise vortex interaction. The purpose of these basic investigations is to remove other irrelevant factors for the parametric studies. Their findings contain the properties of vortical structure both upstream and downstream of the undisturbed shock front, vortex breakdown types (bubble and spiral types) and schematic breakdown model, etc. In the author's opinion, the most important aspect is that Delery et al. [15] presented the shock-induced vortex breakdown limit curve based on the extensively experimental measurements for the normal shock wave/supersonic vortex (NSV/SV) interaction. A key vortex intensity parameter τ max = ðV θ /V ∞ Þ max (V θ and V ∞ represent the tangential and freestream velocity within the vortex diameter, respectively), which is the maximum swirl ratio of the vortex, is introduced to determine the vortex stability limit for a given freestream Mach number. However, due to the spatial restriction for the experimental measuring devices in the compressor rotor tip region, the satisfactory resolution of the flow structures that can be used to deeply understand the physical mechanism of vortex breakdown is unavailable. Meanwhile, the more general and reliable criteria for the shock-induced vortex stability in a highlyloaded transonic axial compressor is still missing. Yamada et al. [16] attempted to apply the above-mentioned vortex breakdown limit curve for NSV/SV interaction to detecting the TLV breakdown in the NASA Rotor 37 at near stall condition. The same intension can be seen in the numerical investigation on the influence of SW/TLV interaction in a stage environment performed by Zhang Zhuoxun et al. [17]. Considering that the application of vortex breakdown limit curve to a realistic transonic axial compressor can only give a relatively qualitative discrimination, it is necessary to combine with other parameters those enable to visualize the flow state in the vortex core for mutual verification. The value of relative Mach number within the disturbed TLV core just downstream of the interaction location was decreased to almost zero, which indicated the possible appearance of backflow zone and its consequent vortex breakdown, as discussed in [16]. As described in the literatures [16][17][18], the distribution of the absolute vorticity normalized by rotor angular velocity shown that it was concentrated in the vortex core when the TLV is stable and it became dispersed in a breakdown state. As a more quantitative parameter, the normalized helicity was capable of assessing the nature of vortex that it implied the occurrence of vortex breakdown when its sign changed downstream of the passage shock in comparison with the upstream one, this detection method was also commonly utilized in many researches [16][17][18][19]. From the above aspect of judging the shock-induced vortex breakdown in a comprehensive manner, a strong deceleration of vortex flow, a drastic expansion of vortex diameter, the dispersed absolute vorticity in the vortex core and the change in the sign of normalized helicity are all regarded as the distinctive features of vortex breakdown.
The main object of the present work is to extract a widely acceptable tip vortex model and to clarify some distinct characteristics of vortex breakdown or not. The most significant aspect is that a theoretical criterion is applied to the shockinduced vortex stability as well as the verification of its effectiveness for the first time. The research findings in this paper are all based on the steady RANS numerical investigations of two well-known typical transonic axial fan rotors, NASA Rotor 67 and NASA Rotor 37, which have different load levels.

Research Object and Numerical
Analysis Method 2.1. About NUAA-Turbo 2.0 Code. NUAA-Turbo 2.0 code is a CFD analysis software specially developed for the highprecision numerical simulation of complex flow in the turbomachinery. Based on the long-term efforts made by our research group to devote to the deep understanding on the internal flow mechanism of impeller components in the aeroengine, the function modules of the code are regularly enriched and optimized. Currently, this code not only meets the research needs of our group but also has been popularized and used by some aeroengine institutes in China. For the preprocessing module, it can support the grid data files generated by a commercial software AutoGrid IGG and obtain any desired grid partition strategies for various mesh topologies to prepare for subsequent parallel computation based on MPI (Message Passing Interface) technique. With regard to the core solver, various numerical methods (including the numerical schemes with different accuracy and the boundary conditions, etc.) can be provided for the flow fields having different natures. To visually present and analyse the predicted results, a variety of flow variables in the form of profiles and contours enable to be outputted into the Tecplot software for post-processing. More functions will be expanded and improved in the future.  [21]. A schematic representation of meridional flow-path for this fan rotor is shown in Figure 1. In this figure, AP and RP represent the axial and radial positions, respectively. Two aerodynamic survey positions both upstream (labeled as "AERO STATION 1") and downstream (labeled as "AERO STATION 2") of the rotor and several laser anemometer measuring locations along the spanwise are clearly identified. Specific measurements, detailed data processing method and final experiment results, which contains the overall performance and the radial distributions of total pressure, total temperature and flow angle, are also completely documented in [21].

NASA Rotor 37.
A transonic axial compressor rotor with a low-aspect-ratio design, NASA Rotor 37, is usually used as the CFD code assessment case since the ASME/IGTI 39th International Gas Turbine Conference [22]. It consists of 36 blades, which has an aspect ratio of 1.19, a tip solidity of 1.288 and an inlet hub-tip radius ratio of 0.70. At the design wheel speed of 17188.7 rpm (tip speed of 454.14 m/s), it leads to the inlet Mach number of 1.13 at the hub and 1.48 at the tip. The design pressure ratio is 2.106 at a mass flow of 20.19 kg/s with a choked mass flow of 20.93 kg/s. The running tip clearance is 0.356 mm (0.5% span). The detailed specification of its basic design parameters are listed in Table 2. Figure 2 shows a meridional view of the blade passage to describe the aerodynamic and laser anemometer survey locations. The radial distributions of static and total pressure, total temperature and flow angle are measured by aerodynamic probes at two axial positions, designated as Stn 1 and Stn 2. The flow field is also surveyed at several blade-toblade planes for revealing the flow structures within the blade passage. Details of experimental investigations on Rotor 37 performed by Suder can be seen in [23].

Computational Grid
2.3.1. NASA Rotor 67. In this study, the structured H-type mesh topology is employed to generate the computational domain for a single blade passage of NASA Rotor 67. All grids are automatically generated by IGG AutoGrid5 in the commercial software NUMECA. The total cells are about 1,200,000 with 120 nodes in the streamwise direction (65 nodes on the blade), 90 nodes in the spanwise direction, and 65 nodes in the pitchwise direction. Particularly, 17 nodes are reserved in the radial direction below the casing wall as the tip clearance for its periodic boundary condition treatment. Employing this simpler H-grid topology with the periodic gap approximation for this case is in view of the following considerations: 1) The H-grid solver was developed earlier than the "butterfly" grid solver and is computationally inexpensive; 2) It is easier to present the distributions of parameters PanBie and Delta on the pitchwise planes by the post-processing software Tecplot; 3) Most importantly, our previous computation experience found that the prediction results of Rotor 67 have almost the same flow structures (especially in the tip region) between these two tip gap treatments. The simulation accuracy for this "periodic gap" is also guaranteed by its code validation (seen in section 3.1.1). A mesh clustering has been done near the all solid wall regions to provide enough resolution for solving the viscous boundary layer flow. The first grid spacing at all solid walls is set to 1 × 10 -3 mm, which can give y + (or z + ) <2. To accelerate the convergence rate, a single H-type grid is divided into 8 subgrid blocks having approximately equal mesh size for load balance in the MPI-based parallel computation. The   2.3.2. NASA Rotor 37. Unlike Rotor 67, the structured 4HOtype mesh topology is adopted to discretize the computational domain of the main flow region. The whole singlepassage grids are also generated by a same interactive grid generator as the Rotor 67. The total number of grid cells, with the tip clearance region, are about 1,860,000. Especially for the O-type grid around the blade, it consists of 298 nodes around the blade, 35 nodes in the direction perpendicular to the blade surface and 100 nodes in the spanwise direction. The fully gridded tip clearance, called "butterfly" mesh, is utilized with 17 nodes in radial direction to accurately capture the tip leakage flow field. The grid distributions near all solid walls are fine enough to directly model the viscous turbulent flow in the boundary layer without wall function method. The value of y + (or z + ) of the first node away from the wall can also obtain sufficient near wall grid resolution like Rotor 67, namely y + (or z + ) <2. Total 20 subgrid blocks, which make up the whole rotor passage fluid domain incorporating the tip clearance gap, are given for parallel computation based on the MPI technique. The detailed informations about grid distribution are presented in Figure 4. In the present studies, the steady, three-dimensional, compressible RANS equations in a conservative form are solved in terms of primitive variables to carry out the following relevant numerical investigations. A cell-centered finite volume method is used to discretize all derivative terms in space by numerically integrating the above governing equations over the closed boundaries of a control volume. The temporal derivative term is discretized by employing a first-order accurate Euler scheme with implicit treatment of convective fluxes and explicit treatment of diffusive fluxes. The Roe's approximate Riemann solver having high-resolution of discontinuity is applied to evaluating the inviscid fluxes, which is extended to third-order accurate by the WENO reconstruction. To approximately determine the viscous fluxes, all spatial derivative terms within them are numerically calculated in a central differencing manner. The nonlinear system of equations achieved through above all discretization are linearized by the iterative Newton's method, then they are solved using the lower-upper Symmetric Gauss-Seidel (LU/SGS) relaxation iteration at each pseudo-time step. A more suitable turbulence model for the adverse pressure gradient flow, so-called Shear-Stress Transport (SST) turbulence model without any modification, is utilized to estimate the eddy viscosity as well as for the final closure of RANS equations.
With respect to the boundary conditions set in the both cases, "ghost" cell technique [19] is introduced to deal with all related boundary conditions. At the inlet boundary, constant uniform distributions of stagnation quantities (total temperature and pressure) and on inlet swirl are assumed and specified by using the characteristic-based boundary condition. The periodic boundary conditions are applied to the axisymmetrical re-entry boundaries for a single-passage computation domain. Adiabatic and no-slip boundary conditions are imposed on all smooth solid walls. For the outflow boundary, the radial equilibrium condition with the back pressure prescribed at the radial position of the casing is utilized to obtain the pressure distribution at exit. The setting of the boundary conditions adopted for the turbulence quantities are as follows: the turbulence intensity of 1% and turbulent eddy viscosity of 100 μ l are specified at the fully turbulent inflow boundary; at the exit boundary, the zero gradients are assumed for all turbulent quantities.

A Steady-State Analysis Method of the Vortex Stability.
This section introduces a steady-state analysis method of the vortex stability, which is obtained from a series of theoretical derivations early reported by Zhang Hanxin and Deng [22].       International Journal of Aerospace Engineering It provides a more accurate and deep understanding for the behavior of streamline near the vortex core on the cross section perpendicular to the vortex axis. A brief description of this method is as follows: In order to simplify the research issue, assuming that the flow field has the steady-state and incompressible natures and is analyzed in the Cartesian coordinate system. Meanwhile, it is prescribed that the vortex axis along the z coordinate axis is straight, the other two coordinate axes x and y are located on the cross section of a slender vortex. The vortex core is denoted as a point "o". By analyzing the characteristics of the streamline equation on a crossflow plane normal to the vortex axis, the conditions for the existence of vortex movement and for its stability are mathematically proved. To further reveal the exact implications of these theoretical criteria, the steady incompressible NS equations are combined to derive the more specific discriminants which can clearly reflect their physical influence factors. Resultantly, two main discriminants form a complete methodology for the judgement of vortex stability, defined as follows: where the velocity components u, v, w correspond to the directions along the coordinate axes x, y, z, respectively. The parameters p, ρ and υ separately represent the static pressure, density and kinematic viscosity coefficient. The character string "PanBie" is a simplified symbol for the expression on the existence of swirl flow, numbered expression (1). The expression (2) labeled as "Delta" is responsible for determining whether the vortex is stable. In the light of  7 International Journal of Aerospace Engineering their research conclusions [22], the swirl flow movement can be detected when PanBie >0. Under this premise, the streamline near the vortex core displays a stable spiral for the case of Delta <0, otherwise the opposite.

Results and Discussions
3.1. Code Validation. In order to verify the prediction accuracy of the simulation code used in the present numerical investigations, detailed comparisons between the predicted and measured results are made at two operating conditions, namely the near-peak efficiency and near stall conditions. They are labeled as "NPE" and "NS", respectively. Particularly, so-called "numerical stall" operating point of the steady simulation is determined by two criteria: one is that the variation of inlet mass flow rate is greater than 0.3% of the average of the maximum and minimum value in 1000 iterative steps (It usually presents continuous and sharp decrease in   International Journal of Aerospace Engineering it); the other is that the difference between inlet and exit mass flow rate exceeds 2.5% of the inlet one. As the abscissa of the overall performance curves, the mass flow is normalized by their corresponding choking mass flow (34.72 kg/sec for Rotor 67 and 20.83 kg/sec for Rotor 37). The current grid resolutions for both simulation cases are sufficient to achieve the grid independence solutions, which can be con-firmed from the evidence that the same mesh topologies with less total number of grids used in the literature [24] can obtain the almost unchanged performance results. The above statements are all suitable for both Rotor 67 and Rotor 37, furthermore, their available detailed experiment data have been reported by Strazisar et al. [21] and Suder [23], respectively. Next, some validations and discussions  10 International Journal of Aerospace Engineering will be conducted with respect to above two representative simulation cases.
3.1.1. NASA Rotor 67. The overall characteristic curves at 100% speed obtained from the simulation and experiment are shown in Figure 5. The both results have a good agreement within most mass flow range, although the computational total pressure ratio and adiabatic efficiency are a little higher than experimental ones in the range of 97%~98% mass flow. In addition, the numerical stall point is located at a little lower mass flow compared to the experiment one. With regard to the fact that RANS simulation results have a certain degree of discrepancies with the experimental data, it largely attributes to the inability of capturing some unsteady effects such as SW/TLV and shock/suction surface boundary layer interactions. Figures 6 and 7 present the radial distributions of measured and computed total pressure ratio, total temperature ratio, adiabatic efficiency and absolute flow angle at nearpeak efficiency and near stall conditions, respectively. Their comparisons show that the predicted results can accurately capture the experimental quantitative trends. The only apparent disagreement is a slight underestimation of the absolute flow angle near the tip region, which is probably due to a periodic boundary treatment of the tip clearance effect. It is well-known that a periodic boundary ignores the real flow in the tip clearance region, which generally brings about the changes in mass, momentum and energy across the blade tip. Neglecting the blockage introduced by the vena contracta leads to a smaller value of absolute flow angle near tip region.
The relative Mach number distributions on the blade-toblade planes at 30%, 70% and 90% spans from hub for two operating conditions, are given in Figures 8 and 9. Closer examinations of these comparisons indicate that the flow structures from the simulations have a strong resemblance to those of experiments, as is the location and its change of the shock wave system. Overall, the satisfactory accuracy of the current code well supports the subsequent numerical researches for Rotor 67.
3.1.2. NASA Rotor 37. The overall performance curves at design speed, which are numerically and experimentally obtained, are displayed in Figure 10. By careful observations and comparisons of their results, the predicted trends are nearly consistent with the measured values. In terms of the visible discrepancies between them, the computational results are slightly underestimated in comparison with the  11 International Journal of Aerospace Engineering experimental data in the whole stable flow range. Although the prediction accuracy gradually improves with the operating point marching toward the stall condition, it indicates in general that the shock system is not well reproduced, which promotes us to further improve the solver performance in predicting it under a higher loading level.
The simulation assessments pertaining to the spanwise distributions of total pressure ratio, total temperature ratio and adiabatic efficiency are presented in Figures 11 and 12, which are, respectively, for near-peak efficiency and near stall conditions. It can be found that the simulation results of all quantities qualitatively follow the same trends as the experimental ones at both operating points. Especially, their adiabatic efficiency have a quantitatively perfect consistency under two loading levels. However, comparing with the measured results, a little higher values predicted by the current code are quantitatively examined in the total pressure ratio and total temperature ratio profiles. This is most likely due to neglecting the shock-induced unsteady behavior of tip leakage vortex, which commonly brings about greater blockage, entropy increase and premature triggering of flow instability.
Furthermore, in order to test the ability of the code to predict the more complex intra-blade flow structures, the rel-ative Mach number distributions on the blade-to-blade planes at 70% and 95% spans are given for the comparison purpose, as shown in Figures 13 and 14. The prediction can successfully capture the passage shock location as well as its strong interaction with the tip leakage vortex. Nevertheless, the radial influence of this aerodynamic interaction is somewhat underestimated by the numerical solutions. This can also be confirmed from a locally weaker bulged-forward shock shape in the simulation. As stated, the predicted accuracy of the present numerical methods strengthen the confidence in investigating the complex aerodynamic phenomena near the tip region.

Blade Tip Vortex Model.
As known, the complex secondary flow or vortex structures in the tip region of a highly-loaded transonic axial compressor, which always have a significant influence on the blockage, flow loss and stability. So, establishing a generally accepted blade tip vortex model is beneficial for deeper understanding of the internal flow mechanism and then reorganizing them to obtain some positive effects. For getting a more convincing result, two representative transonic axial fan rotors having different loading levels, the Rotor 67 and Rotor 37, are selected to analyze and summarize the dominant vortex structures near the tip

12
International Journal of Aerospace Engineering region at near-peak efficiency and near stall conditions. All vortical structures are identified employing the widely used Q-criterion [25]. Additionally, several abbreviations in Figures 15 and 16, which are not mentioned above, are interpreted as follows: PS-pressure surface、SS-suction sur-face、LE-leading edge、TE-trailing edge.
Firstly, Figure 15 clearly shows the large-scale vortex structures which govern the flow field in the tip region at two operating conditions for Rotor 67. All vortices are colored with the relative Mach number to display their flow states. At near-peak efficiency condition, two distinct tip leakage vortices TLV1 and TLV2 are formed from the front part of the blade chord. Both of them occur to merge at a certain axial position and then develop downstream together. One vortex called "SIV" is originated from the location of the passage shock/suction surface interaction in the tip region. Due to the strong adverse pressure gradient created by the shock effect, the tip leakage flow is decelerated and accumulated in the shock region (reflecting the shockinduced role), which is "rubbed" with the casing boundary layer to roll up a vortical structure. It develops toward the pitchwise direction with an angle with the suction surface, which is larger than those of TLV1 and TLV2. Similarly, the SIV meets with the above mixed tip leakage vortex of TLV1 and TLV2, again merging with it into the bigger one. On the suction side near the trailing edge, the tip leakage and separation flows are rolled together under the shear of casing wall to form a vortex named TSV. Increasing the back pressure to near stall operating point, except for the apparent changes in the size and spatial relative location, the types of the vortical structures in the blade tip vortex system nearly maintain constant. As the shock front is pushed upstream, the origin of TLV2 gradually replaces that of TLV1 for becoming a stronger single tip leakage vortex. As a result, more weaker and slender vortices adjacent to this single TLV are produced by the strong interaction with the casing boundary layer, that is why they are called TLV-IV. Responding to the forward movement of the shock front, the TSV covers a larger axial range along the suction surface.
Most interestingly, the tip vortex system for Rotor 37 has a strong resemblance to that of Rotor 67 although it has a higher loading level, as indicated in Figure 16. At both operating conditions, four dominant vortex structures possessing their own generation mechanisms, including the TLV、TLV-IV、SIV and TSV, can also be observed in the tip region. With respect to some discrepancies in comparison (1) (2) 13 International Journal of Aerospace Engineering with Rotor 67, mainly reflected in the following aspects: 1) only one TLV accompanied with its TLV-IV are generated at near-peak efficiency condition; 2) the vortex core undergoes a dramatical expansion after crossing the strong passage shock at near stall condition (It indicates the appearance of vortex breakdown, which is confirmed not only by   14 International Journal of Aerospace Engineering some previous research results but also by our subsequent studies). Consequently, based on the comparative analyses of the blade tip vortex structures for above two typical transonic axial fan rotors, a tip vortex model containing the four types of vortices those frequently coexist is proposed in Figure 17. Due to limited research object, further validation of its generality is needed based on more extensive highlyloaded transonic compressor rotors in the future.

Characteristics of Shock-Induced Vortex Stability.
In this section, several qualitative discriminant methods of the shock-induced vortex stability, frequently employed in a highly-loaded transonic axial fan rotor, are combined to examine the flow pattern of the TLV from different dimensions. The characteristics of the TLV after passing through the passage shock, which can also be revealed from these flow variables and physically visual images. Generally, the relative Mach number、axial velocity、absolute vorticity、normalized helicity、3D tip leakage streamline and the visualization technique of vortex structure, are widely used in assessing the nature of shock-disturbed TLV by many scholars [26,27]. For convenience of describing the flow field, all flow parameters in the following figures are presented in a dimensionless form. They are defined as follows: where V ! , ω ! , u ∧ , U r and L r are the relative velocity vec-tor、absolute vorticity vector、axial velocity with dimen-sion、reference values of velocity and length, respectively. Besides, V、u and H denote the dimensionless values of the magnitude of absolute vorticity、nondimensional axial velocity and the normalized helicity. Next, the pre-and postshock developments and flow properties of the tip leakage vortex for both Rotor 67 and Rotor 37 at two typical operating conditions, which are detailedly analyzed and summarized by above identifiable flow variables and 3D tip leakage streamlines. Figures 18 and 19 display the vortex structures colored with relative Mach number and normalized helicity as well as the distributions of absolute vorticity and axial velocity on the four crossflow planes for Rotor 67 at both operating conditions, to reveal the characteristics of shock-induced vortex stability. It is found from Figure 18(a) that the incoming supersonic TLV 1 and TLV 2 undergo the deceleration after they cross the passage shock, but the shock is not strong enough to cause the flow stagnation. The tight rolling-up 3D tip leakage streamlines can also indicate that the TLV always keeps stable from leading to trailing edges. This point well corresponds to the concentrated vorticity near the vortex core shown on the cross-sectional planes. Figure 18(b) depicts all vortices colored with normalized helicity, the value of it near the vortex core is almost 1 upstream of the shock but only slightly decreased downstream of it. It implies that the stability of the TLV still keeps constant under the current blade load level. This can also be confirmed by the fact that the TLV passes through the areas with positive axial velocity on all four crossflow planes. As the operating point is marched to the near stall condition, the vortex diameter becomes larger although no recirculation flow in the postshock TLV. Seen from the distributions of absolute vorticity on the plane 2 located immediately downstream of the shock, the shock-disturbed vortex core still behaves the same concentrated vorticity as that on the plane 1 upstream of the shock. These flow states reflecting the vortex stability are consistent with the tight rolling-up TLV colored with relative Mach number, as shown in Figure 19(a). As observed in Figure 19(b), the high positive values of the normalized helicity along the TLV can still be remained even after interacting with the shock, accompanying with no reversed flow detected from the distributions of axial velocity in the vortex region. It further demonstrates that the stability of the TLV has not suffer from apparent changes in nature.
Unlike Rotor 67, the Rotor 37 possesses a so higher blade load level that the possibility of vortex breakdown is increased. Respectively, Figures 20 and 21 present the more strong SW/TLV interactions at near-peak efficiency and near stall conditions, which are depicted with the same flow variables as those used in Rotor 67. At near-peak efficiency condition, although a sudden deceleration and the moderate expansion happen in the TLV due to such aerodynamic interaction, the tip leakage flow is still capable of rolling up a complete and stable vortical structure with concentrated vorticity at the vortex center, as indicated in Figure 20(a). The high positive value of normalized helicity along with the distribution of axial velocity in the vortex region further strengthen this judgment, as seen in Figure 20(b). However, the local stagnation (almost 0 value of relative Mach number in Figure 21(a) and slightly negative axial velocity in Figure 21(b)) flow near the vortex core region on the plane 2, which is positioned just downstream of the shock, takes place at near stall condition, as presented in Figure 21. In addition, the drastic expansion of the vortex core also gives rise to a essential change in the vorticity distribution at the vortex center. It means that such vorticity distribution is characterized by a "ring-type" distribution, namely low in the center and high around (without concentrated vorticity at the vortex center), as described on the plane 2 in Figure 21(a). Correspondingly, the value of normalized helicity drops from almost 1 of pre-shock undisturbed TLV to partly negative value downstream of the shock, this can be explained the occurrence of an abrupt change in the nature of vortical structure, as validated in Figure 21(b). Unlike some experimentally parametric investigations on the normal (or oblique) shock/streamwise vortex interaction [28], the shock wave and vortex intensity will usually strengthen simultaneously as the operating condition advances to the

16
International Journal of Aerospace Engineering near stall point in a transonic axial compressor. It is the double-enhancement effect that deteriorates the tip leakage flow and even triggers its breakdown. Although the shock-induced vortex stability is evaluated based on the steady-state calculation, it does not affect the assessment of an identification method introduced in the current study, mainly based on the following reasons: a) The focus of this paper is mostly on the physical mechanism of SW/TLV interaction rather than its unsteady flow behavior and process. In other words, the physical mechanism of this aerodynamic interaction remain constant at every instant, but what changes is just its flow state (stable or unstable); b) In the steady solutions of two computation cases, the key features of two TLV flow patterns mentioned above enable to be observed, which has almost consistency with those shown in some transient process of previous time-   17 International Journal of Aerospace Engineering accurate studies [16]; c) The reliability evaluation of the detection technique emphasized in the next section is done by comparing with these traditional methods.
As stated, some distinct characteristics of the shockinduced vortex stability, especially the vortex breakdown, can be obtained from the numerical studies of above two typical transonic axial fan rotors having different load levels. It is clearly realized that the vortex breakdown, which generally emerges just downstream of the shock, is always accompanied with a sharp deceleration to almost stagnation (nearly 0 value of relative Mach number or slightly negative axial velocity) at the center of a bubble-type structure formed by the drastic vortex core expansion, the dispersed absolute vorticity distribution in the vortex core (like a "ring-type" distribution) and the occurrence of partly opposite helicity sign to that of upstream undisturbed vortex. These comprehensive flow phenomena provide a more multidimensional consideration for judging whether the vortex is bursting or not.

Validation of a Discriminant Method for Shock-Induced
Vortex Stability. Although the previous section presents the some qualitative descriptions and distinctive features on the shock-induced vortex stability, a theoretical discriminant which can more detailedly reflect the influence factors of vortex stability has hardly been seen in the open published literatures. This section will focus on introducing a vortex stability criterion while developing along the vortex axis, early proposed by Zhang Hanxin and Deng [22]. Then, in the author's knowledge, it is the first time to apply this discriminant to the SW/TLV interaction in the transonic axial fan rotor. It can be known from the previous introduction of this complete methodology given by equations (1) and (2) that under a big premise of the vortex existence identified by equation (1), its stability dominantly depends on the sign of equation (2). The value of Delta is determined by three influence aspects, seen from the right side of equation (2), they represent the axial velocity、pressure gradient and viscous dissipation. The NASA Rotor 67 and NASA Rotor 37

18
International Journal of Aerospace Engineering are also taken as the verification objects to check the feasibility of these detection parameters in comparison with the analysis results of previous section. For NASA Rotor 67, in order to have a better visualization display for the tip vortices and their relative locations, six and four pitchwise sections are selected for near-peak efficiency and near stall conditions, as shown in Figures 22(a) and 23(a), respectively. At near-peak efficiency condition, all six pitchwise planes are sequentially marked as 1-1, 2-2, ......, 6-6 from suction to pressure side, namely they correspond to 30%, 44%, 48%, 50%, 58% and 61% pitch, respectively. It is observed that high positive values of PanBie are concentrated in several zones where the tip vortices exist, as seen in Figures 22(b)-22(g). This means that the parameter  PanBie is capable of capturing the existence of vortical structures, and it can still remain the concentrated and high positive values after crossing the passage shock. The judgment based on the PanBie that the TLV still keeps stable just downstream of the shock can also be confirmed by the following parameter Delta. As indicated in Figure 24, the value of Delta is always negative within the vortex diameter range even after passing through the shock. All these characteristic phenomena demonstrate that the TLV has no obvious instability at near-peak efficiency point. As the operating point  With regard to NASA Rotor 37, four and six crosssectional planes are chosen for examining the ability of parameter PanBie to capture the TLV at near-peak efficiency (in Figure 26) and near stall (in Figure 27) conditions, respectively. In the following analyses, only the distributions of PanBie in the TLV region are needed to pay attention to, but not other areas. When the fan rotor runs at the nearpeak efficiency point, as shown in Figure 26, it is clearly found that the value of PanBie is not only positive but also highly concentrated near the vortex core before and after the passage shock. The only differences are that the magnitude and distribution scope of its value gradually become smaller with the downstream evolvement, due to the shock effect and viscous dissipation. But the nature of TLV stability is unchanged. As the mass flow rate is throttled to the near stall point, the PanBie distribution inside the TLV experiences the three stages from leading to trailing edges: high concentration、incomplete breakdown and slight concentration stages. The first stage is defined as the highly concentrated distribution of PanBie in the vortex region, which manifests that the vortex is quite stable, as shown on the

21
International Journal of Aerospace Engineering plane 1 (upstream of shock) in Figure 27(b). Just downstream of the shock, from plane 2 to plane 4, a "ring-type" distribution of PanBie, which is similar to that of vorticity on plane 2 in Figure 21(a), appears with a great change in the nature of vortex structure. This phenomenon means the occurrence of local instability in the vortex core region, but not the complete breakdown. It can also be evident from the gradual recovery of PanBie to weak concentration distribution, as presented on planes 5 and 6 further downstream of shock. Nevertheless, from the observation of Delta distributions along the TLV trajectory, they have strong resemblance to each other at both operating conditions. As displayed in Figure 28, the parameter Delta in the TLV has a higher negative value upstream of the shock compared to the downstream one, which sharply increases to almost zero after crossing the shock. Although there is no obviously positive value of Delta in the expansion region of TLV, this just can illustrate that the post-shock vortex behaves the weak stability not the complete breakdown.
According to the above series of validations about the applications of parameters PanBie and Delta to detecting the shock-induced vortex stability, they give the basically same conclusions as the identification methods utilized in previous section but have a more explicit physical implications. Notablely, on the one hand, these two detection parameters are yielded through a series of theoretical derivations which leads to possessing a more physical connotation, three main influencing factors on the right side of their equations can be further used to evaluate which factor dominates the vortex stability (from this point of view, they are superior to the common parameters like the helicity); on the other hand, parameter PanBie performs better than Delta from the results of their effectiveness validations, the incomplete burst of TLV based on the steady-state simulation is the most likely reason for the poor performance of Delta. In the future research, further verifications of their detection abilities in an unsteady flow environment will be conducted, especially the behavior of parameter Delta when the vortex is completely broken.

Conclusions
In this paper, the steady RANS simulations based on NUAA-Turbo 2.0 solver, are performed to investigate the application of an identification method proposed by Zhang Hanxin and Deng [22] to analyzing the shock-induced vortex stability in two well-known transonic axial fan rotors, NASA Rotor 67 and NASA Rotor 37. Before that, a relatively general blade tip vortex model are proposed and some distinctive characteristics of the shock-induced vortex stability are also comprehensively summarized as a reference for the methods verified in this study. The findings are summarized as follows: (1) A general blade tip vortex model in a transonic axial fan rotor consists of four main vortex structures: TLV、TLV-IV (occurs when the tip leakage vortex/casing boundary layer interaction reaches to a certain strength)、SIV and TSV, which dominantly govern the flow field in the tip region (2) Based on several flow variables (relative Mach number、axial velocity、absolute vorticity and normalized helicity) and 3D tip leakage streamline, frequently used by many researchers to assess the post-shock flow state of TLV, a vortex instability or breakdown is normally characterized by the following features: a sudden deceleration to almost stagnation (nearly 0 value of relative Mach number or slightly negative axial velocity) at the vortex center, the drastic vortex core expansion accompanied with dispersed absolute vorticity distribution and the occurrence of locally opposite helicity sign relative to that of pre-shock undisturbed vortex (3) The parameter PanBie not only enables to detect the existence of vortical structure but also has the ability of identifying the occurrence of shock-induced vortex instability, which behaves a "ring-type" distribution of PanBie in the vortex core region.

22
International Journal of Aerospace Engineering not show obvious positive value due to no complete vortex breakdown takes place even at near stall condition, it can clearly reveals the development process of vortex stability strength along the vortex trajectory and indicates that the TLV presents strongly weak stability after crossing the passage shock Nomenclature RANS: Reynolds-averaged Navier-Stokes SW/TLV: Shock wave/tip leakage vortex SIV: Shock-induced vortex TSV: Tip separation vortex TLV-IV: Tip leakage vortex-induced vortex NSW/SV: Normal shock wave/supersonic vortex NPE: Near-peak efficiency NS: Near stall PS: Pressure surface SS: Suction surface LE: Leading edge TE: Trailing edge x: x-direction coordinate y: y-direction coordinate z: z-direction coordinate u: Velocity component in the x-axis direction v: Velocity component in the y-axis direction w: Velocity component in the x-axis direction p: Static pressure ρ: Density υ: Kinematic viscosity coefficient μ l : Molecular viscosity m: Mass flow rate m_choke: Choking mass flow rate y + : Nondimensional distance of the first node away from the hub or casing wall z + : Nondimensional distance of the first node away from the blade surface τ Max : Maximum swirl ratio V θ : Tangential velocity V ∞ : Freestream velocity V ! : Relative velocity vector ω ! : Absolute vorticity vector u ∧ : Dimensional axial velocity U r : Reference value of velocity L r : Reference value of length V: Dimensionless value of the magnitude of absolute vorticity H: Normalized helicity PanBie: Parameter representing the existence of vortex Delta: Parameter representing the vortex stability.

Data Availability
The data used to support the findings of this study are available from the corresponding author upon request.