UNSTEADY BODY FORCE MODEL FOR ROTATING STALL IN AXIAL COMPRESSOR WITH VARIOUS INLET CONDITIONS

To analyze the dynamic stall of multistage axial compressors, a three-dimensional unsteady numerical model is established based on the body force model. For a two-stage axial compressor with a clean inlet, the calculated maximum steady static pressure rise coeﬃcient is only 0.1% diﬀerent compared with the experimental data. The characteristic frequency of the dynamic stall evolution basically agrees with the experimental results, which proves the eﬀectiveness of the model. For the compressor with a combined radial-circumferential total pressure inlet distortion, the predictions preliminarily verify the ability of the model for qualitative description of the ﬂow instability process with the complex inlet distortion.


Introduction
Aerodynamic instability characteristics of areoengine compressors due to the inverse pressure gradient of the internal viscous gas directly affect the operational reliability of engines (Day, 2016).The pursuit for a design of a highly loaded compressor has heightened the problem of flow stability.Meanwhile, with gradual application of new concepts such as boundary layer ingesting (BLI) jet engines (Uranga et al., 2014) in future aero propulsion systems, the compressor must continuously operate with a complex inlet distortion, which has a significant effect on the compressor stability.It is imperative to establish efficient and re-liable models to predict compressor instability in various inlet conditions.
The aerodynamic instability of the compressor is mainly manifested in two unsteady phenomena: rotating stall and surge (Day, 2016).Compared with surge, which is an axial oscillating unsteady flow pattern caused by the interaction between the compressor and downstream components with the cavity effect, rotating stall features more three-dimensionality.In the research area of axial compressors, the rotating stall has been the focus of attention.
One of the most direct approaches for predicting the dynamic stall process of the axial compressor is the three-dimensional unsteady Reynolds-averaged Navier-Stokes (URANS) method, which treats the flow stability as an initial boundary value problem and simulates the whole dynamic stall process (Huang et al., 2019a,b;Wang et al., 2020;Zhang and Vahdati, 2020).It is undeniable that the URANS simulation can effectively describe the multi-scale flow interaction in turbomachinery and play a positive role in the understanding of the instability mechanism.The method and other higher-order computational fluid dynamics (CFD) methods will be central to the forecasting technique of the flow stability of the compressor in the future.However, as for the current computing level, the full-annulus unsteady CFD simulation of multistage compressors is time-consuming, which leads makes it difficult to be popularized as a conventional technology in engineering practice.Therefore, the development of an efficient stability analysis model is still of great practical significance to improve the fidelity of the existing compressor design system.
An alternative approach of the flow stability analysis is to focus on the modeling of the main scale flow in the compressor within the framework of the time marching technology of CFD.Central to the modeling is to describe the blade effect in reasonable reduction of dimensionality to lower the computational cost.According to the ideology, the actuator disk model (Hu et al., 1999) and body force model (Gong et al., 1999;Longley, 2007;Righi et al., 2018) have been established.In the actuator disk model, the blade row is treated as a discontinuity and the flow process inside the blade passage is not solved.While in the body force model, the boundary effect of the blade is modeled as a distributed force source terms to effectively simulate main flow characteristics in the blade passage.Gong et al. (1999) developed a three-dimensional incompressible model for analyzing the dynamic stall process of multistage axial compressors based on the body force model.Longley (2007) investigated the modeling method for the blade force in the reverse-flow condition and established a three-dimensional unsteady compressible model to analyze the rotating stall and surge of a multistage axial compressor.Righi et al. (2018) proposed a highly robust three-dimensional unsteady throughflow model for the rotating stall and surge in an axial compressor based on the improved Gong's blade force model.Recently, Rosa Taddei (2021) developed a meanline body force approach to simulate the two-dimensional flow behavior during the rotating stall in a transonic compressor.
The aforementioned models focus on simulation of the stall dynamic process of a compressor within uniform inflow.Numerous numerical studies (Cao et al., 2017;Guo and Hu, 2018; Mao and Dang, 2020) have examined the effectiveness of the body force model for simulating steady transfer characteristics of a circumferential large-scale distortion in the axial compressor.There are rare comprehensive investigations on the simulation of the dynamic stall process for multistage compressors with an inlet distortion through the body force approach.The construction method of the blade force plays a significant role in this kind of model, especially, a comprehensive and effective correlation between blade force source terms and local aerodynamic and geometric parameters under inlet distortions remains to be addressed elaborately.
This study aims to develop a three-dimensional unsteady numerical model to analyze the dynamic stall process of a multistage axial compressor in various inlet conditions.The distributed blade force source terms are established by concepts of the deviation angle and loss coefficient in through flow theory, which can be suitable for stability evaluation of the compressor in the throughflow design stage.The model is implemented by an in-house program and capability to quantitatively describe the dynamic stall process of a compressor with a clean inlet is evaluated on a two-stage axial compressor at Nanjing University of Aeronautics and Astronautics (NUAA).Besides, the potential of the model for efficient simulation of dynamic stall processes with complex inlet distortions is demonstrated in detail.

Theoretical methods
The flow in the compressor is governed by a three-dimensional unsteady compressible equation in the absolute cylindrical coordinate system.The boundary effect of blades is modeled by distributed force source terms and a blockage coefficient of the blade profile, instead of a traditional body-fitted mesh approach.The main governing equation of the model is as follows In the non-bladed region, the blade blockage coefficient b is equal to The term e * and h * denotes the total energy and the total enthalpy per unit mass of the gas, respectively.The term τ is the turbulent Reynolds stress and q is the turbulent heat flux.The specific expressions are as follows 3) The assumption of a calorically perfect gas is employed for the model, and the ideal gas law p = ρR g T is satisfied.The value of c p and γ are assumed to be constant.The calculation methods for the eddy viscosity coefficient µ t and the turbulent heat conductivity coefficient k t are presented in Section 2.2.
In the bladed region, the terms b, S b and S F are expressed as where N is the number of blades.The terms θ p and θ s are the circumferential coordinates of the pressure and suction sides of the blades, respectively.The term F represents the instantaneous blade force per unit volume of the flow and Ω is the angular velocity of the rotor, as shown in Fig. 1. (2.5)

Unsteady blade force model
The modeling of the unsteady blade force F plays a central role in the developed model.A method of the blade force construction with a clear physical meaning is employed in the present work.A quasi-steady blade force is constructed directly from the equation of steady circumferential momentum.Subsequently, an instantaneous blade force is modeled by a first--order lagging response function.This provides an effective approach to consider the effect of the radial flow.Meanwhile, the robustness of the model is enhanced because there is no need to introduce the artificial iteration coefficient.
To begin this process, the blade force is decomposed into two parts: the turning force F T perpendicular to the relative velocity u rel and the loss force F L parallel but opposite to u rel by referring to Marble (1964), as shown as where The value of F L is determined by the local loss coefficient where u m = √ u 2 + w 2 is the instantaneous meridional velocity and u rel is the magnitude of relative velocity.The term ∆s denotes the entropy increase from the inlet to the outlet of the blade passage and ∆m represents the meridional length of the blade passage.For the rotor passage, the terms Ma and T * are the relative Mach number and the relative total temperature, respectively.For the stator passage, the two terms signify the absolute values.Subscripts 1 and 2 indicate the inlet and outlet of the blade passage, respectively.
The total circumferential force F θ is obtained by the steady circumferential momentum equation The imposed circumferential velocity v * in the bladed region is determined as follows where κ is the blade metal angle and δ is the deviation between the actual flow direction and blade metal angle due to the flow viscous effect and the centrifugal force of the curved blade passage.
The radial turning force F T,r is correlated to the three-dimensional blade mean camber surface.Hence, the end bend effect is reflected in this model to a certain extent, as shown in the equation where n r and n θ are the radial and circumferential components of the normal vector of the blade mean camber surface.
Up to now, every component of the quasi-steady blade force can be determined when the local loss coefficient ̟ and the deviation angle δ are obtained.When the gas moves forward normally, according to the deviation angle and loss model in the through flow theory of the compressor (Johnsen and Bullock, 1965;Creveling and Carmody, 1968), the local deviation angle and the loss coefficient are calculated by the incidence, Mach number, and typical blade geometry at different circumferential and radial positions.To modify and broaden the incidence characteristics of different elementary cascades, numerical blowing experiments can be employed by the three-dimensional RANS method.
When the compressor is throttled to a low mass flow rate near stall, the axial velocity in the flowpath may be negative.However, it should be noted that under the local circumferentially averaged theory, the reverse flow is not necessarily presented.In the reversed flow region, similar to the analysis of Longley (2007) and Righi et al. (2018), the loss force is no longer calculated.The turning force is limited to be perpendicular to the local mean camber surface of the blade instead of the relative flow direction.The benefit of this is that the robustness is guaranteed and the high loss is naturally introduced by adjusting the gas flow direction.The outlet flow angle should be specified at the blade inlet in the reverse flow.In the current numerical study, it can be found that the dynamic stall process obtained by the model basically agrees with the experimental data when the outlet direction of the flow is consistent with the blade stagger angle.
More importantly, in the dynamic stall process of the compressor, it suffers a delayed dynamic response of the local blade force to the upstream aerodynamic disturbance.To effectively describe the formation of the rotating stall, this dynamic effect should be quantified in the model.
Because the circumferential turning force is dominant, this term is modified with the unsteady response by a first-order lagging response function, as shown in the equation where F ss T,θ is the quasi-steady circumferential turning force, and λ = c m /u m represents the lagging response time.The term c m denotes the meridional blade chord length.

Numerical solution
A cell-centered finite-volume method is utilized to discretize Eq. (2.1).The time terms are solved using an explicit four-stage Runge-Kutta scheme.In the unsteady calculation, the global minimum time step according to the CFL number is selected as a unified time step of the computational domain.The inviscid flux of the space term is calculated by the low-diffusion flux-splitting scheme (LDFSS(2)) proposed by Edwards (1997).
The discretization of the viscous flux in Eq. (2.1) is much easier due to its elliptic characteristic compared with that of the inviscid flux.As can be found in Eqs.(2.3), the viscous flux at the interface can be determined only by obtaining the eddy viscosity coefficient µ t , the turbulent heat conductivity coefficient k t , and the partial derivatives of the velocity and static temperature.The model focuses on the energy exchange between different streamlines caused by the turbulent diffusion in the main flow instead of the small-scale viscous flow near the wall.The aerodynamic loss in the flow field is mainly loaded through the loss force term.Hence, the eddy viscosity coefficient µ t is calculated by a simple semi-empirical method proposed by Gallimore and Cumpsty (1986).Namely, the criterion number µ t /(ρuL ref ) is determined in advance in the global flow field.Subsequently, the local µ t is carried out according to the local instantaneous flow field.The parameter L ref is the compressor characteristic length.The turbulent Prandtl number in the global flow field is given as 0.9 and combined with the µ t to calculate k t .

Boundary conditions
To ensure the robustness of the model, the inlet and outlet boundaries of the computational domain are treated as open boundaries.Thus, the gas is allowed to flow in and out freely.When the gas normally flows into the computational domain, the total pressure, total temperature, and flow direction distribution are given at the domain inlet.The total pressure specified at the inlet would be treated as the static pressure where the flow exits the domain inlet.The outlet pressure is calculated by the classical throttling model based on the quadratic parabolic analytic function, as shown in where p e is the static pressure at the domain outlet.The parameters p * 0 and ρ 0 are the average total pressure and density of the inflow, respectively.The term U m is the rim speed at the midspan of the compressor inlet and K t is the throttle coefficient.The term φ = u e /U m denotes the flow coefficient at the compressor outlet.If the reverse flow is present at the outlet, p e will be regarded as the total pressure.
In the condition of the clean inlet, some circumferential non-uniform artificial initial small disturbances are introduced to induce the stall initiation process.Working at the steady condition, a small disturbance will disappear quickly without any influence on the compressor flow field.While the compressor is throttled to the stall point, the small disturbance will expand nonlinearly and eventually evolve into the mature rotating stall cell.In this model, a sinusoidal total pressure distortion of amplitude 0.1% is applied at the inlet boundary.The disturbance is canceled after 0.1 of rotor revolutions.

Numerical results and discussion
In this Section, a two-stage axial compressor at the laboratory located in NUAA is employed to quantitatively evaluate the ability of the model to effectively simulate the dynamic process of the rotating stall in a multistage compressor.A steady and dynamic experiment has been performed on the compressor within the clean inlet to collect abundant data to validate the developed model.Meanwhile, to further demonstrate the potential of the model for efficiently simulating the dynamic stall process with a complex inlet distortion, the model is applied to predict the stall flow behavior of the two-stage axial compressor with the combined radial-circumferential total pressure inlet distortion.
Typical aerodynamic and geometric parameters of the compressor are listed in Table 1.

Model validation for dynamic stall process under clean inlet
The dynamic stall process of the compressor has been measured at 800 rpm. Figure 2 presents a schematic of the test rig.The steady pressure parameters of the compressor inlet are measured by four static pressure taps on the casing wall and six six-hole total pressure combs (Fig. 2, plane I and plane II).The steady pressure parameters of the compressor outlet are measured by 16 static pressure taps on the casing and hub and eight total pressure combs (Fig. 2, plane III and plane IV).To investigate the dynamic process of the compressor from the steady-state to rotating stall, six single-point differential-pressure probes are distributed along the circumference of plane II to monitor the perturbation of the total pressure at 85% span of the compressor inlet.The pressure-sensing holes of the probes are facing the flow direction.High-frequency response Kulite sensors are selected as the pressure sensor.It should be noted that the model is based on the circumferentially averaged theory of the local flowpath and the governing equation is solved in the cylindrical coordinate system.Therefore, the solution is theoretically not sensitive to quantity of the circumferential mesh.This has also been proved by Mao and Dang (2020).Particularly, in calculating the steady performance of the compressor with the clean inlet, simply one grid is set in the circumferential direction.Figure 3b plots the computed characteristic curves of the steady static pressure rise coefficient based on two sets of meridional grids.The static pressure rise coefficient ψ st is defined as where p 2 represents the average static pressure at the last-stage stator outlet of the compressor.It can be found that characteristics obtained by two different sets of grids are nearly the same.The coarse mesh is selected for the later full-annulus unsteady computation to further shorten the time cost of simulating the dynamic stall process.The detailed computational mesh is shown in Fig. 4. The compressor is throttled by adjusting the throttle coefficient until the rotating stall occurs.Figure 5a gives a comparison between the calculated and experimental results of the pressure rise characteristics of the compressor at different throttle openings.It can be seen that the overall characteristics obtained by the model are in line with the experimental results.The comparison shows only about a 0.1% difference in the maximum static pressure rise coefficient between the calculated and experimental results.
Figure 5b further illustrates the time history of the flow coefficient during the compressor stall obtained by the model.As can be found in Fig. 5b, the throttle coefficient increases from 2.89 to 2.92 near the stall point and the dynamic stall process begins.After about 30 rotor revolutions, a fully developed rotating stall is formed.Subsequently, the average overall characteristic point of the compressor remains constant.
To further reveal formation of the stall cell, axial velocity contours at the midspan of the compressor in different rotor revolutions are displayed in Fig. 6.A slight closing of the downstream valve at the stable boundary point has a great influence on the development of internal flow.The generation of the stall cell is at the second-stage stator and it gradually expands to the whole compressor inlet.A stable single cell, about a third of the circumferential length of the compressor, is finally formed.Fig. 6.Formation process of the stall cell in the compressor obtained by the model Figure 7 shows the total pressure traces at different circumferential locations at the 85% span of the compressor inlet during the stall obtained by the experiment and calculation.As can be observed, the main characteristics of the dynamic stall evolution obtained by calculation basically agree with the experimental results.The stall precursor of the compressor is a typical large--scale modal wave disturbance.The experimental results show that the propagation frequency of the final stall cell is 42% of the rotor speed (5.6 Hz).While the prediction of the propagation frequency is 38% of the rotor speed (5.1 Hz).Overall, the characteristic frequency of the dynamic stall evolution obtained by the model is similar to that from the experimental data.

Numerical prediction of dynamic stall process for complex distortion
To exhibit the potential of the model for quantitative description of the dynamic stall process with a complex inlet distortion, the dynamic stall process of the compressor with a combined radial-circumferential total pressure inlet distortion is simulated by the developed model.Its capability to describe the aerodynamic instability process of the multistage compressor with the distorted inlet is temporarily investigated by numerical simulations at this stage.Besides, the rationality of calculation results is discussed based on the previous qualitative consensus.Figure 8 plots the contour of the total pressure distortion at the computational domain inlet.The low-pressure distorted region occupies a 50% span near the tip in the radial direction and 120 • in the circumferential direction.The absolute total pressure in the distorted region decreases by 0.5%.Fig. 8. Distribution of total pressure at the inlet A large-scale circumferential nonuniform disturbance exists at the compressor inlet with the distorted inflow, so the artificial small disturbance is no longer introduced in the model.Figure 9 provides the calculated time-averaged pressure rise coefficient of the compressor with the inlet distortion.It can be found that the distortion results in a significant decrease of the peak pressure rise coefficient and makes the stall condition in advance.The partially enlarged drawing near the stall stage is further displayed in Fig. 9b.As can be seen, the calculation results show that the compressor gradually transits from the full stable state (point A) to the small--scale oscillation state (points C, D, E) during the throttling process with the inlet distortion.Point B is the critical stable state when the compressor enters the small-scale oscillation region.In this region, the average pressure rise of the compressor no longer increases with the decrease of the flow coefficient and reaches the peak near the critical point B. With further throttling, the oscillation amplitude of the flow parameter increases gradually until the full deep stall (point F ) appears.
The transfer characteristics of the combined radial-circumferential distortion in the two-stage compressor are analyzed by taking the full stable state (point A) as an example.Figure 10 plots the dimensionless axial velocity distribution at the rotor inlet of each stage.Figure 11 presents the total pressure distribution at the stator outlet of each stage.Due to the low-pressure distorted region at the tip, the axial velocity near the side of the rotor away from the distorted region (240 • ) is smaller than that at other circumferential positions of the same blade height while the corresponding incidence is larger and the blade aerodynamic load increases.Therefore, the maximum total pressure is obtained at the root of the stator outlet.Correspondingly, the axial velocity near the rotor entering the distorted region (120 • ) is larger, the corresponding incidence is smaller, and the blade aerodynamic load decreases.The minimum total pressure is at the tip of the stator outlet.Moreover, as the distorted region transfers downstream, the extreme regions of flow parameters are significantly shifted along the rotating direction due to rotor rotation.The evolution of the internal disturbance in the compressor under a small-scale oscillation is analyzed taking point E as an example.Figure 12 displays the mass flux traces at different circumferential locations near the tip of each stage outlet.
Before the deep stall, some local flow disturbances are generated in the region with the large incidence on the side of the rotor leaving the distorted region.The flow disturbances are suppressed and even disappeared in the region with the small incidence on the side of the rotor entering the distorted region.Therefore, the circumferential propagation of the local flow disturbance presents great periodicity due to the circumferential inlet distortion.The stability of the whole compressor depends on the balance between the disturbance decay in the Figure 13 displays the total pressure traces near the tip of the compressor inlet at different circumferential locations during the stall.Figure 14 shows the time-frequency characteristics of the total pressure near the tip in the non-distorted region (0 • ) and distorted region (180 • ).The amplitude-frequency characteristic of the dynamic pressure signal is calculated by the fast Fourier transform (FFT).It is found that the frequency of the stall cell predicted by the model with the combined distortion is 4.7 Hz (35% of the rotor speed), which is close to the value (38% of the rotor speed) with the clean inlet.This phenomenon suggests that the frequency of the final stall cell is not significantly influenced by the distorted inlet.Its frequency is mainly determined by the inherent parameters of the compression system.The conclusion is also supported by the early experimental results of Fortin and Moffatt (1990).Besides, it can be seen in Fig. 14 that the amplitude of the stall cell in the distorted region is stronger than that in the undistorted region.Despite the lack of experimental data, the predicted results are consistent with the basic understanding of the early classical experimental research, and the ability of the model to quantitatively describe the flow instability process with the complex inlet distortion can be preliminarily verified.

Conclusions
This study develops a three-dimensional unsteady numerical model based on the body force model to analyze the dynamic stall process of a multistage axial compressor in various inlet conditions.The boundary effect of blades is transformed into distributed force source terms in the model.A quasi-steady blade force is constructed directly from the steady circumferential momentum equation combined with the deviation angle and the loss model in the through flow theory.Subsequently, an instantaneous blade force is modeled by a first-order lagging response function to quantify the dynamic response of the blade force to the upstream disturbance during the dynamic stall of the compressor.
According to the abundant experimental data of a two-stage axial compressor with a clean inlet, the capability of the developed model to quantitatively describe the dynamic stall process in a multistage compressor is evaluated.The calculated main performances and characteristic frequency of the dynamic stall evolution are in line with the experimental results, which proves the model has great effectiveness in predicting the compressor performance and the dynamic stall process with a low computational cost.
The dynamic stall process of the compressor with a combined radial-circumferential total pressure inlet distortion is also simulated by the developed model, which reveals that the model has potential for describing the dynamic stall process with the complex inlet distortion.The calculation shows a large non-uniformity of the rotor inlet incidence due to the inlet distortion in the critical stall condition, which has a considerable influence on the development of the local flow disturbance before the stall.The flow disturbance formed in the region of the high incidence decays sharply in the region of the low incidence.Once the balance between the disturbance decay in the non-distorted region and the disturbance growth in the distorted region is disturbed, a full rotating stall occurs.Besides, the frequency of the stall cell predicted by the model with the combined distortion is close to that with the clean inlet.The above-predicted results are consistent with the basic understanding of the early classical experimental research.The ability of the model for qualitative description of the flow instability process with the complex inlet distortion can be preliminarily verified.
Further experimental investigation on the inlet distortion will be promoted to improve simulation of the model.Besides, the time marching scheme for unsteady simulations can be improved by introducing acceleration techniques when the model is used for compressors with more stages and computational grids.Meanwhile, the proposed model focuses on the compressor with the long-scale modal wave stall initiation.The description of the short-scale disturbance stall initiation process of the high-speed compressor is not involved and remains to be further improved in the future work.

Fig. 1 .
Fig. 1.Schematic of aerodynamic and geometric parameters in the bladed region According to the circumferentially averaged theory of the local flowpath, the value of F v is 0 in the bladed region.Then the circumferential inviscid flux vector F is transformed into F = Ωr[ρ, ρu, ρv, ρw, ρe * ] T (2.5)

Fig. 2 .
Fig. 2. Schematic of the test rig and measurement arrangement The dimensionless process is carried out for the collected dynamic total pressure by p * = p * (t) − p * p * (3.1)where p * represents the average value of the sampling period.To focus on the main low-order frequency signals in the dynamic stall, low-pass filtering is utilized to eliminate the interference of high-order frequency signals in the flow field of the compressor.

Fig. 3 .
Fig. 3. (a) Schematic of the computational domain for the compressor in the model.(b) Meridional mesh independency study Figure 3a displays the computational domain of the two-stage compressor in the unsteady body force model.The dependency of the solution for the developed model on the mesh density has been assessed based on the steady performance of the compressor.The cell number along the radial, circumferential, and streamwise directions are expressed as I, J and K, respectively.

Fig. 4 .Fig. 5 .
Fig. 4. Three-dimensional computational mesh of the compressor in the model

Fig. 7 .
Fig. 7. Total pressure traces at 85% span of the compressor inlet at different circumferential locations during stall: (a) experimental result, (b) numerical results

Fig. 9 .
Fig. 9. Comparison of the pressure rise coefficient between the clean inlet and inlet distortion: (a) time-averaged pressure rise coefficient of the compressor, (b) partially enlarged drawing

Fig. 10 .
Fig. 10.Distribution of axial velocity at the rotor inlet of each stage at point A: (a) at the R 1 inlet, (b) at the R 2 inlet

Fig. 11 .Fig. 12 .
Fig. 11.Distribution of total pressure at each stage outlet at point A: (a) at the S 1 outlet, (b) at the S 2 outlet

Table 1 .
Design parameters of NUAA two-stage compressor