Time-Marching Throughflow Analysis of Centrifugal Compressors with Boundary Conditions Based on Newton’s Method

Featured Application

This method can be used to quickly predict the overall performance as well as the meridional flow details of the centrifugal compressors. After being verified, it can provide an efficient and robust analysis tool for the design system of different types of centrifugal compressors.

Abstract

The meridional distribution of the flow parameters inside the centrifugal compressor is of great importance to its overall performance, as well as its matching performance under a thermal cycle. A time-marching throughflow method for the off-design performance analysis of the centrifugal compressor is described. The method is based on the strictly conservative throughflow-governing equations, and an improved method of boundary-condition enforcement is developed based on Newton’s method to achieve a robust and fast throughflow simulation. An inviscid blade force model was adopted to obtain the flow deflection inside the blade passage. Empirical loss models were integrated into the throughflow model to simulate the viscous force effects in the real three-dimensional flow. Two test cases are presented to validate the throughflow method by comparisons with the experimental data or CFD results, including the NASA low-speed centrifugal compressor (LSCC) and the Allison high-performance centrifugal compressor (HPCC). The simulation indicated that the developed enforcement method for the inlet and outlet boundary conditions significantly improves the computational robustness. For both the LSCC and HPCC cases, reasonable flow-parameter distribution was obtained and accurate overall characteristics were also predicted under the off-design conditions. The results indicated that the developed time-marching throughflow method is effective and efficient for the performance analysis of centrifugal compressors.

1. Introduction

Turbomachinery design is among the most challenging tasks in modern industrial design mainly due to the highly complex, non-linear physical environment in which the turbomachines operate. This task demands that the most advanced methods and tools are employed in order to gain the necessary understanding of flow phenomena and to exploit the flow physics to achieve optimum turbomachine performance. To cope with this, the hierarchical design system of turbomachinery has been developed consisting of the low-order models and high-fidelity models, such as the one-dimensional (1D) mean-line method, the two-dimensional (2D) throughflow method, and the three-dimensional (3D) method. Among them, the throughflow method has become the backbone of turbomachinery design. Benefitting from its good balance between computation efficiency and accuracy, the throughflow method has been commonly adopted in the preliminary design and optimization stage both for axial [1,2,3,4,5,6] and radial turbomachinery [7,8,9,10]. As has been reported, up to 80% or more of the design parameters may be fixed during the throughflow modeling phase [11]. Apparently, an improved throughflow method is still of great interest to the industrial community.

Since the concept of the S2 stream surface was introduced by Wu in the 1950s [12], several kinds of throughflow methods have been developed based on different mathematic approaches. These approaches can be mainly categorized into three kinds: the stream function method (SFM) first introduced by Marsh [13], the streamline curvilinear method (SCM) introduced by Smith and Novak [14,15] and the time-marching-based method pioneered by Spurr [16]. In the first two methods, usually called the classical throughflow flow methods, the original set of governing equations are reduced to a single differential equation: the stream function equation for the SFM and the radial equilibrium equation for the SCM. Though they are different numerical approaches in principle, it has been concluded by Davis that there is little to choose between them [17]. As the common drawbacks, the simplifications make these approaches conceptually unable to predict choking conditions or capture shock waves [18,19]. This causes striking difficulties in the throughflow analysis of the turbomachinery, especially when the flow inside the blade passage becomes transonic. Although several strategies were implemented [6,8,20] to enhance their performances in transonic flow analysis, the classical throughflow approaches are inherently weak in these situations.

To overcome the major drawbacks of the two classical throughflow approaches mentioned above, the time-marching throughflow method has been proposed, in which the complete set of Euler/Navier–Stokes equations are discretized and solved by advancing in pseudo time on the S2 surface. Benefitting from the progress in the theory of computational fluid dynamics (CFD), the time-marching throughflow method can automatically predict the choking conditions, and it has the shock-capturing feature. This method provides an efficient and robust analysis tool for the subsonic, transonic, and supersonic flow and has been becoming more and more popular in the modern turbomachinery design and analysis.

The first time-marching throughflow method based on the Euler equations was presented by Spurr [16]. Then, systematic works were carried out by Baralon [21,22,23] and some key points were studied, including the modeling of inviscid blade force and blade blockage, and the analysis of the properties of shock waves. Gu and Anderson [24,25] have been dedicated to developing a time-marching throughflow analysis tool, which incorporates the effects of viscous loss, flow deviation, air bleeding and injection, and throat choking. Topp [26] attempted to construct a turbomachinery analysis and design system (TADS) at Allison Engine Company by coupling the time-marching throughflow and blade-to-blade models, which is similar to Spurr’s work. Further developments led to the Navier–Stokes-based throughflow model and progress in modeling the effects of 3D loss in Simon and Léonard’s work [27,28]. The potentials of the time-marching throughflow method were also investigated in the design tasks of both the axial compressor and turbine based on the solution of the inverse problem [29,30,31] and its coupling with advanced optimization methods [32,33,34]. Detailed works about axial turbine simulation were performed by Ricci et al. [35], aiming to incorporate the real gas capabilities and 3D flow features into the time-marching throughflow models. Yang et al. [36] raised a novel inviscid blade force model and proved its superiority in computing robustness compared to the widely used ordinary-differential-equation-based one.

Although the aforementioned works have shown the considerable superiorities of the time-marching throughflow method compared to the classical ones, its applications are still limited within the scope of axial turbomachinery. To the author’s best knowledge, the time-marching throughflow analysis of the centrifugal compressor has not been implemented in the open literature. This may be attributed to the very strong adverse pressure gradient and complex viscous effects commonly existing inside the centrifugal compressor. The former can result in a very low Mach number flow and even backflow, which significantly increases the stiffness of the throughflow-governing equations; meanwhile, the latter negatively affects the accuracy of the performance prediction for the centrifugal compressor.

In the present paper, a viscous time-marching throughflow method for the centrifugal compressor simulation is presented. The method has the capability to predict both the overall performance at off-design points and the corresponding meridional flow details. A normally averaged throughflow model, rather than the circumferentially averaged one, was elaborately selected, preserving the rigorous conservation of flow parameters in transonic flow analysis. To overcome the convergence problem occurring in the throughflow simulation of the centrifugal compressors with a high pressure ratio, characteristic boundary conditions (CBCs) based on Newton’s iteration method were developed and can greatly improve the computation robustness. To simulate the viscous effects, complete loss correlations deduced from the open literature were adopted, accounting for both the internal losses and parasitic losses inside the centrifugal compressor. Owing to all the efforts above, throughflow solutions of different kinds of centrifugal compressors can be efficiently obtained using the time-marching method.

This paper is organized as follows. In Section 2, the time-marching throughflow flow model of the centrifugal compressor is elucidated exhaustively in terms of the governing equations, the inviscid blade force model, the viscous force model, and the empirical models. Section 3 then presents the specific numerical method utilized in the solver, especially introducing the developed non-reflective boundary condition (NRBC). The developed throughflow solver is validated against experimental and numerical data in Section 4, taking the well-known NASA low-speed centrifugal compressor (LSCC) and the impeller of the high-performance centrifugal compressor (HPCC) as test cases. At last, all the conclusions obtained in this paper are summarized in Section 5.

2. Throughflow Model

2.1. Governing Equations

Governing equations for the time-marching throughflow model can be derived from the 3D viscous equations by taking a simple axisymmetric assumption or a circumferential average. For body-fitted curvilinear coordinates $\left(\xi ,\eta ,\zeta \right)$ satisfying Equation (1),

$${\xi }_{\phi }={\eta }_{\phi }\equiv 0$$

(1)

their governing equations can be written as Equation (2), which represents the throughflow model adopted in most studies, both for the non-bladed duct and the bladed region.

$$\frac{\partial }{\partial t}\left(\frac{rU}{J}\right)+\frac{\partial }{\partial \xi }\left[\frac{r}{J}\left(F{\xi }_z+G{\xi }_r\right)\right]+\frac{\partial }{\partial \eta }\left[\frac{r}{J}\left(F{\eta }_z+G{\eta }_r\right)\right]=\frac{Q}{J}+\frac{r}{J}{f}_v+f_b+Q_b$$

(2)

However, this throughflow model may encounter the computing difficulty in the transonic flow simulation when shock waves exist inside the blade passage, for which the details have been presented in [37]. To overcome this problem, quasi-orthogonal body-fitted coordinates satisfying Equation (3)

$$\{\begin{array}{c}\xi \cdot \zeta =0\\ \eta \cdot \zeta =0\end{array}$$

(3)

were adopted in the present study for the bladed region. This resulted in the normally averaged governing equations of throughflow in the blade passage, as in Equation (4),

$$\frac{\partial }{\partial t}\left(\frac{rU}{J}\right)+\frac{\partial }{\partial \xi }\left[\frac{r}{J}\left(F{\xi }_z+G{\xi }_r+H\frac{1}{r}{\xi }_{\phi }\right)\right]+\frac{\partial }{\partial \eta }\left[\frac{r}{J}\left(F{\eta }_z+G{\eta }_r+H\frac{1}{r}{\eta }_{\phi }\right)\right]=\frac{Q}{J}+\frac{r}{J}{f}_v+f_b+Q_b$$

(4)

where ${\xi }_z$, ${\xi }_r$, ${\xi }_{\phi }$, ${\eta }_z$, ${\eta }_r$, ${\eta }_{\phi }$, ${\zeta }_z$, ${\zeta }_r$, ${\zeta }_{\phi }$ are the metrics including the blockage effects of blade thickness. $\xi \left({\xi }_z,{\xi }_r,{\xi }_{\phi }/r\right)$,$\eta \left({\eta }_z,{\eta }_r,{\eta }_{\phi }/r\right)$, and $\zeta \left({\zeta }_z,{\zeta }_r,{\zeta }_{\phi }/r\right)$ correspond to vectors of three orientations, respectively, and $J$ is the metric Jacobian. The conservative variables $U$ are defined as

$$U={\left[\rho ,\rho u,\rho v,\rho w,e\right]}^T$$

and the flux vectors $F$, $G$ and $H$ are

$$F={\left[\rho u,\rho u^2+p,\rho uv,\rho uw,(e+p)u\right]}^T$$

$$G={\left[\rho v,\rho uv,\rho v^2+p,\rho vw,(e+p)v\right]}^T$$

$$H={\left[\rho w,\rho uw,\rho vw,\rho w+p,(e+p)w\right]}^T$$

where $\rho$ is density, $p$ is pressure, $e$ is specific energy, and $\left(u,v,w\right)$ are the components of the relative velocity vector $W$. In the right-hand side (RHS) of Equation (4), the term $Q$ contains contributions from the cylindrical-coordinate form of the governing equations, and $Q_b$ accounts for the variation of the area in the blade passage. Both terms are defined, respectively, as

$$Q={\left[0,0,p+\rho {(w+\omega r)}^2,\rho v(w+2\omega r),{\omega }^2{r}^2\rho v\right]}^T$$

$$Q_b={\left[0,-p\frac{\partial }{\partial \zeta }\left(\frac{r}{J}{\zeta }_z\right),-p\frac{\partial }{\partial \zeta }\left(\frac{r}{J}{\zeta }_r\right),-p\frac{\partial }{\partial \zeta }\left(\frac{{\zeta }_{\phi }}{J}\right),0\right]}^T$$

The inviscid blade force term

$$f_b={\left[0,-\frac{r{\zeta }_z}{J}\left(\frac{\partial p}{\partial \zeta }\right),-\frac{r{\zeta }_r}{J}\left(\frac{\partial p}{\partial \zeta }\right),-\frac{{\zeta }_{\phi }}{J}\left(\frac{\partial p}{\partial \zeta }\right),0\right]}^T$$

guarantees a tangential condition between the flow and the mean flow surface inside the blade passage and vanishes in the non-bladed duct.

The viscous force term $f_v={\left[0,f_{vz},f_{vr},f_{v\phi },0\right]}^T$ is introduced to simulate viscous losses. Both the inviscid blade force term $f_b$ and viscous force term $f_v$ need to be modeled to close the system in Equation (4), as will be introduced in the following sections.

2.2. Inviscid Blade Force Model

To simulate the flow deflection effects made by a blade, many attempts have been made, which are generally based on two kinds of ideas, i.e., the relaxed method and the direct method. The former introduces a fictitious evolution equation for $f_b$, then relaxes to the steady state [21,26,38], while the latter usually directly expresses the circumferential component of $f_b$ as the gradient of the circulation [29,34]. However, both of these kinds of methods have their shortcomings. In terms of the relaxed methods, the relaxation factor needs to be manually selected, of which an unsuitable value may significantly degrade the convergence and accuracy of the throughflow calculation. As for the direct method, projection of the circumferential component of the $f_b$ to the other two orientations may cause severe computing instabilities in the case of a large blade angle. To ensure the stability of throughflow simulations for the centrifugal compressor, an inviscid blade force model developed by Yang et al. [36] was adopted. This model derives the modulus of the $f_b$ directly from the discretized normal momentum on the mean flow surface, which has proven to be of great benefit to the solver robustness and can be expressed as

$$f_b=R\cdot \zeta /\left|\zeta \right|$$

(5)

where $R$ is the residual vector term in the RHS of the discretized momentum equations. In the solving process, $f_b$ is calculated at each iteration in conjunction with a slip condition constraint

$$W^{*}=W-\left(W\cdot \zeta \right)\zeta /{\left|\zeta \right|}^2$$

(6)

$W^{*}$ represents an updated relative velocity.

2.3. Viscous Force Model

In the present study, the viscous effects were introduced by reconstructing a distribution of the entropy field. In general, this was achieved by adopting a viscous force in the opposite direction of the flow velocity

$$f_v=-f_v\frac{W}{\left|W\right|}$$

(7)

It is natural, then, to find a way to relate the modulus of the viscous force model $f_v$ with the entropy distribution. This problem has been addressed by Bosman and Marsh [39], in whose work a derivation of Crocco’s equation, incorporating the conservation of total enthalpy under an absolute frame and rothalpy under a relative frame, gives the specific expression as

$$f_v=-\frac{\rho TW\mathbf{·}\nabla s}{{\left|W\right|}^2}W$$

(8)

Notice that only the relative velocity was adopted here because it naturally degrades into the absolute velocity in an absolute frame.

2.4. Empirical Models

To reconstruct a physically reasonable entropy field in the S2 stream surface for the centrifugal compressor off-design performance analysis, the complete empirical correlations available in the open literature were adapted to the present throughflow method. Viscous effects from both the impeller and vaneless diffuser were considered. Moreover, the slip models were employed in the present throughflow model to simulate the deviation between the real flow angle and the blade angle. Referring to Oh’s work [40] and the authors’ previous work [41], a list of the empirical models already implemented in the present throughflow solver is given in Table 1. These selected empirical models were validated in the aforementioned works and showed satisfying performance in terms of compatibility and thus the accuracy of the analysis results.

Table 1. Empirical correlations list in centrifugal compressor simulation.

Categories of Models	Correlations
Slip factor	Qiu et al. [42]
Incidence loss	Aungier [43]
Blade loading loss	Coppage and Dallenbach [44]
Skin friction loss	Jansen [45]
Clearance loss	Jansen [45]
Mixing loss	Johnston&Dean [46]
Disc friction loss	Daily and Nece [47]
Recirculation loss	Yang et al. [48], Oh et al. [40]
Vaneless diffuser loss	Coppage and Dallenbach [44]

Table 1. Empirical correlations list in centrifugal compressor simulation.

Categories of Models	Correlations
Slip factor	Qiu et al. [42]
Incidence loss	Aungier [43]
Blade loading loss	Coppage and Dallenbach [44]
Skin friction loss	Jansen [45]
Clearance loss	Jansen [45]
Mixing loss	Johnston&Dean [46]
Disc friction loss	Daily and Nece [47]
Recirculation loss	Yang et al. [48], Oh et al. [40]
Vaneless diffuser loss	Coppage and Dallenbach [44]

In general, most of these empirical correlations were established by correlating the losses to the flow parameters at the blade leading and trailing edges. Therefore, special treatment is required when they are applied to the throughflow simulation, which means a redistribution of the evaluated losses along the spanwise and streamwise directions is necessary.

2.4.1. Spanwise Loss Distribution

In terms of the loss distribution in the throughflow simulation of the axial compressor, some studies have been implemented, in which the results seem to be basically matched with the 3D simulation or the real flows. However, the loss distribution in the throughflow simulation of the centrifugal compressor has not been noticed enough in the open literature. Additionally, difficulties are encountered in the application of the empirical correlations along the spanwise direction, since almost all the empirical correlations for the centrifugal compressors were originally established and verified in a 1D perspective, as concluded by Oh et al. [40] and Zhang et al. [10]. To deal with this problem, the spanwise-averaged/meanline parameters were employed to evaluate the “total” flow loss in the present solver. This “total” loss was then distributed with a prescribed form along the spanwise direction according to the specific type of loss. A uniform loss distribution along the spanwise is advised by Casey [8], and was adopted in the distribution of incidence loss, blade loading loss, skin friction loss, and mixing loss. Then, considering the source of loss production, a simple locally linear distribution was chosen for the clearance loss, disc friction loss, and recirculation loss. As for the vaneless diffuser loss, a parabolic distribution was employed to approximate the viscous effects from real endwall effects. The conditions prescribed for the linear and parabolic function are as follows:

For the linear distribution, once the distribution range is specified:

• The minimum value of the loss at one end is set as zero;
• The integral of the loss along the span is equal to the estimated “total” loss by a 1D approach.

For the parabolic distribution:

• The spanwise location of the minimum loss value needs to be manually inputted and is usually chosen at mid-span;
• The ratio of the minimum loss value to the averaged one needs to be manually inputted. It takes the value of zero for the vaneless diffuser loss distribution in this study;
• The integral of the loss along the span is equal to the estimated “total” loss by a 1D approach.

Both the linear and parabolic distribution forms along the span are depicted in Figure 1. The values that have to be manually entered in these distributions, such as the spanwise location and value of the minimum vaneless loss with a parabolic distribution, are specified with the view that they can represent some basic features of the meridional field of a compressor. In this way, the throughflow solver can be adopted as an independent analysis tool in the preliminary design stage with moderate accuracy in terms of meridional flow field prediction. The accuracy of these distributions can be further improved by the results from the 3D viscous computations and experimental measurements.

2.4.2. Streamwise Loss Distribution

In addition to the spanwise loss distributions, the loss distributions along the streamwise direction have to be prescribed. To simulate the entropy field as real as possible, either a linear distribution or a cosinoidal distribution referring to Pacciani’s work [49], which is simple but still realistic, can be selected for each kind of loss in the present throughflow solver.

2.5. Blade Leading and Trailing Edge Treatment

As pointed out by many researchers [21,50], an inviscid blade force discontinuity happens at the blade leading edge in the time-marching throughflow calculation if no special treatment is implemented under an off-design performance analysis for the turbomachinery. Accompanied by this discontinuity, unphysical entropy is produced and is directly proportional to the mismatch degree between the real inflow angle and the inlet metal angle. This kind of unphysical entropy can pollute the predicted flow fields and make an accurate off-design performance analysis impossible. Meanwhile, the deviation angle, which also originates from the mismatch between the real flow surface and the blade camber surface, needs to be simulated at the trailing edge. Casey [8] proposed that in the SCM calculation, the actual value of the swirl at the blade inlet and the deviation angle at the outlet should be iteratively updated, which is essentially a flow-regularization procedure. A similar but slightly different flow-regularization approach proposed by Yang et. [36] was employed in the present time-marching throughflow solver. This approach eliminates the inviscid blade force discontinuity at the blade leading edge by smoothing the flow field. It also provides a realistic outflow angle by incorporating the well-accessed deviation or slip factor models. When it was applied to the present throughflow model, prescribed portions at both the leading and trailing edges of the blade camber surface were locally modified to regularize the flow in that region. A cubic spline function was adopted in the prescribed portions so as to ensure that the modified camber surface regions remained C1 continuous.

3. Numerical Methods

The present throughflow code, S2CFD, was programmed with the Fortran 95 language and supports multi-block computation by adopting a zonal method and thus relieving the difficulty of the grid generation in simulations of complex configurations. In this code, the governing equations of Equation (4) were solved by a numerical procedure based on a Godunov-type finite volume method [51]. Solutions were obtained by marching at the pseudo-time steps until the steady state was reached. Second-order accuracy was obtained by using MUSCL reconstruction in conjunction with a limiter function. The governing equations system can be advanced in pseudo-time either by using an explicit four-stage Runge–Kutta scheme, or an implicit LU-SGS scheme. The local time-stepping method was adopted both in the explicit and implicit schemes. Residual smoothing was applied to the explicit scheme to speed up convergence to the steady state.

Improved Boundary Condition Enforcement

The numerical enforcement of the boundary condition plays an important role in the CFD solver development and can significantly affect the computational efficiency and robustness of the marching process in a steady simulation. A commonly used characteristic boundary condition (CBC) for perfect gas involves the so-called “Riemann-invariant” at the inlet and adopts linearized characteristic equations at the outlet, of which the details and an example can be found in reference [52]. However, in practice, the writers found that difficulties in convergence were encountered when the Riemann-invariant-based CBCs were enforced in the time-marching throughflow simulation of the centrifugal compressor. Especially for the high-pressure-ratio cases, the computation can blow up even near the peak efficiency operating point. To ensure a stable computation, a method of imposing the CBCs based on Newton’s iteration method was developed for the throughflow method with the help of the local one-dimensional inviscid (LODI) relations. The developed CBCs refer to the methodology in [53], but were adapted for the throughflow equations. Usually, four values need to be specified at the inlet boundary for the subsonic inflow in the throughflow simulation, therefore a residual vector can be defined as

$$R=\left[\begin{array}{c}{R}_1\\ R_2\\ R_3\\ R_4\end{array}\right]=\left[\begin{array}{c}{h}_t-h_{t0}\\ s-s_0\\ v-u\tan \sigma \\ w-u\tan \alpha \end{array}\right]=\left[\begin{array}{c}dh\\ ds\\ v-u\tan \sigma \\ w-u\tan \alpha \end{array}\right]$$

(9)

where $h_{t0}$ and $s_0$ are the specified total enthalpy and entropy, respectively, and $h_t$ and $s$ are the corresponding calculated values. $\sigma$ and $\alpha$ denote the radial and tangential flow angle. The residual vector $R$ is dependent on the primitive variables at the boundary and should vanish after the computational convergence. Using Newton’s method:

$$R+\frac{\partial R}{\partial U_p}d{U}_p=0 \Rightarrow \frac{\partial R}{\partial U_p}d{U}_p=-R$$

(10)

$U_p$ represents the primitive variable $U_p={\left[\begin{array}{c}\begin{array}{ccccc}p& u& v& w& T\end{array}\end{array}\right]}^T$. According to the characteristic theory, another characteristic variable propagating outside the inlet is added here as:

$$d{C}_5=d{u}_n-\frac{dp}{\rho c}=0$$

(11)

where $u_n$ is the local velocity normal to the inlet boundary. Equation (11) represents a left-propagating acoustic wave. Combining Equations (9) and (11), a system can be presented as:

$$A_{in}\cdot d{U}_p=b_{in}$$

(12)

with

$$A_{in}=\left[\begin{array}{ccccc}0& u& v& w& c_p\\ -\frac{R_g}{p}& 0& 0& 0& \frac{c_p}{T}\\ 0& -\tan \sigma & 1& 0& 0\\ 0& -\tan \alpha & 0& 1& 0\\ 1& -\rho c{n}_z& -\rho c{n}_r& -\rho c{n}_{\phi }& 0\end{array}\right],b_{in}=\left[\begin{array}{c}-R_1\\ -R_2\\ -R_3\\ -R_4\\ 0\end{array}\right]$$

The analytical solution of Equation (12) can be obtained as

$$dp=\frac{\rho dH-\frac{pds}{R}-\frac{\rho u_s{u}_n}{\cos \delta }+\rho \left(u^2+v^2+w^2\right)}{1+\frac{u_s}{c\cdot \cos \delta }}$$

(13)

$$d\rho =\frac{\rho }{\gamma }\frac{dp}{p}-\frac{\rho ds}{c_p}$$

(14)

$$d{u}_n=\frac{dp}{\rho c}$$

(15)

$s\left(\frac{1}{A},\frac{\tan \sigma }{A},\frac{\tan \alpha }{A}\right)$ denotes the unit vector of the prescribed flow direction at the inlet boundary. Where $A=\sqrt{{\tan }^2\sigma +{\tan }^2\alpha +1}$, $u_s=W\cdot s$, $\cos \delta =n\cdot s$.

For the subsonic outlet boundary, one characteristic variable propagates into the computational domain and the other four propagate outside. In the time-marching throughflow, the static pressure is usually imposed at the outlet boundary, and thus can be chosen as the inward propagating characteristic

$$R_5=p-p_b$$

(16)

where $p_b$ denotes the specified backpressure. In terms of the outgoing characteristic variables, four correlations, which represent entropy wave, right-propagating acoustic wave and vorticity waves, respectively, are added and can be expressed as

$$dT-\frac{dp}{\rho c_p}=0$$

(17)

$$d{u}_n+\frac{dp}{\rho c}=0$$

(18)

$$d{u}_{t1}=0$$

(19)

$$d{u}_{t2}=0$$

(20)

where $u_{t1}$ and $u_{t2}$ denote the other two velocity components, which are orthogonal to the normal velocity $u_n$. Similar to that at the inlet, a system can be obtained by combing correlations (17)–(21)

$$A_{out}\cdot d{U}_p=b_{out}$$

(21)

$$A_{out}=\left[\begin{array}{ccccc}1& 0& 0& 0& 0\\ \left(\gamma -1\right)& 0& 0& 0& \frac{-\gamma p}{T}\\ 1& \rho c{n}_z& \rho c{n}_r& \rho c{n}_{\phi }& 0\\ 0& t_{1z}& t_{1r}& t_{1\phi }& 0\\ 0& t_{2z}& t_{2r}& t_{1\phi }& 0\end{array}\right],b_{out}=\left[\begin{array}{c}-R_5\\ 0\\ 0\\ 0\\ 0\end{array}\right]$$

It is not difficult to obtain the analytical solution for Equation (21) as

$$dp=-R_5$$

(22)

$$dT=\frac{dp}{\rho c_p}$$

(23)

$$d{u}_n=-\frac{dp}{\rho c}$$

(24)

which is exactly the same as the subsonic outlet BC correlations shown in [52]. Therefore, the CBCs based on Newton’s method developed in this study actually change the inlet boundary treatment for the throughflow simulation. It calculates the variations of the flow parameters at the inlet boundary, rather than directly calculating their absolute values as conducted in the Riemann-invariant-based method, which is considered to be more robust. Velocities both at the inlet and outlet boundaries can be computed as

$$du=d{u}_n{n}_z,dv=d{u}_n{n}_r,dw=d{u}_n{n}_{\phi }$$

(25)

By using Equations (13)–(15) and (22)–(25), the primitive variables at both the inlet and outlet boundaries can be updated during the throughflow iterations. In the following sections, a comparison will show the superiority of this kind of boundary-condition enforcement compared to the Riemann-invariant-based method.

4. Validation and Applications

In this section, two test cases are presented to validate the robustness and accuracy of the developed time-marching throughflow method in the prediction of the off-design performances and the flow details. The two test cases include the low-speed centrifugal compressor (LSCC) and the high-performance centrifugal compressor (HPCC).

4.1. Low-Speed Centrifugal Compressor (LSCC)

The LSCC is a low-speed centrifugal compressor with a semi-open and backswept impeller that was designed by the NASA Lewis research center, and it has a design performance of 30 kg/s in the mass flow rate and 1.166 in terms of the total pressure ratio. As a validation for the developed throughflow method, the numerical simulation of the LSCC case is compared with the experimental data by Hathaway et al. [54]. The computations were performed on a 365 × 40 grid, which is sufficient for the Euler-based time-marching throughflow computation, as can be seen in the grid-independence validation shown in Figure 2b. For clarity, a coarser grid is depicted in Figure 2a instead of the actual grid.

With the grid, a design-point throughflow simulation costs less than 2 min on an R5-5600H processor to achieve convergence. The convergence histories enforcing different CBCs, i.e., the Riemann-invariant-based method and the developed CBCs based on Newton’s method in this study, are compared in Figure 3. It can be observed that, except for a slight saving of two iterations, the CBCs based on Newton’s method and the Riemann-invariant obtained almost the same convergence curves. This indicates that there exists little difference between these two kinds of CBCs when they are adopted into the computation of the low-speed centrifugal compressor.

The computed overall performances at the design mass flow rate are compared with the experimental data in

Table 2. Computed and experimental LSCC performance on design point.

	Computation	Experiment	Error (%)
Total pressure ratio	1.159	1.141	1.8
Isentropic efficiency	0.919	0.922	−0.3

. The results are in good agreement, with an error of 1.8% for the total pressure ratio and −0.3% for the isentropic efficiency. This shows the adequacy of the present throughflow method of design-point simulation in predicting the overall performance of the LSCC.

A further investigation of the off-design performance of the LSCC was implemented. It was reported by Oh et al. [40] that the recirculation flow model has a great influence on the off-design predictions. Therefore, both the recirculation loss models developed by Oh et al. [40] and Yang et al. [48] were validated in this case. Figure 4 shows the comparison between different settings of the recirculation loss model in terms of both the total pressure ratio and the isentropic efficiency. In Figure 4a, it seems that the recirculation loss has little effect on the total pressure ratio. This is reasonable because as one kind of parasitic loss, the recirculation loss only leads to a change in enthalpy and a corresponding change in efficiency. It can also be seen that the discrepancy between the throughflow results and the experimental data increased with the mass flow rate. Nonetheless, the maximum error was less than 2.5% within all operating ranges. Unlike the total pressure ratio, the off-design performance of the isentropic efficiency was greatly influenced by the recirculation loss model setting, as shown in Figure 4b. When the mass flow rate was larger than the design value, Yang’s recirculation model was in good agreement with the experiment, whereas it seems that no significant loss was produced by Oh’s model and thus resulted in higher isentropic efficiency. For the small mass flow rate, both Yang’s and Oh’s models showed dramatically increased discrepancy compared to the experiment as the mass flow rate decreased. Nevertheless, the trend of the rapid efficiency decrement was predicted by both models. The comparison shown in Figure 4b indicates that a more accurate recirculation flow loss model needs to be developed.

Table 3. Contributions of different losses for the LSCC case.

Mass Flow Rate (kg/s)	Incidence Loss (%)	Blade Loading Loss (%)	Skin Friction Loss (%)	Clearance Loss (%)	Mixing Loss (%)	Disc Friction Loss (%)	Recirculation Loss (%)	Vaneless Diffuser Loss (%)
22.0	7.965	12.389	13.843	2.212	6.005	3.603	49.621	4.362
30.0	3.209	19.652	33.422	5.080	26.738	6.417	0.000	5.481
40.0	0.765	7.957	27.085	3.060	34.124	3.366	21.653	1.989

further reports the detailed contributions of each loss to the total loss in the throughflow simulation at different operating conditions when Yang’s recirculation model was adopted.

The meridional flow field of the LSCC at its design mass flow rate is detailed in Figure 5. As can be seen from Figure 5a, the relative Mach number inside the blade passage remained lower than 0.3, except near the shroud at the impeller inlet, and achieved the lowest value at the impeller exit. This kind of low-speed flow shows a slightly incompressible characteristic. Figure 5b shows the static pressure distribution along the flow passage. Though an overall diffusion process was achieved, the static pressure did not exhibit a monotonous lift. Instead, the static pressure drop was observed at the vaneless diffuser outlet. This was mainly due to the convergent passage area resulting from the dramatic decrease in the diffuser width at the passage outlet, though the radius increased. Figure 5c presents the total temperature rise due to the work done by the impeller on the airflow. It is noticeable that a uniform distribution of the total temperature along the blade span was observed at the impeller exit. The predicted entropy distribution is shown in Figure 5d. It can be seen that along the streamwise direction, the entropy increased continuously. Meanwhile, in the spanwise direction, the entropy was predicted to be of low value at midspan and at a higher value near the endwall. Specifically, the entropy product at the shroud was higher than that at the hub, which can be accounted for by the significant tip leakage loss. Overall, the reasonable flow details presented in Figure 5 demonstrate the adequacy of the developed throughflow method on the meridional flow simulation in the low-speed centrifugal compressors.

The spanwise distribution of the flow parameters at the impeller outlet is always the key point that the designers are concerned with. Therefore, a comparison of the spanwise distribution of total temperature ratio at the design point between the throughflow results and experimental data was implemented and is shown in Figure 6. It can be seen that, except for the nearby shroud region, the calculated total temperature ratio was slightly lower than the experiment below the 25% spanwise location. At the other positions, the calculated results remained higher than the experimental values. This accounts for the higher total pressure prediction at the design point discussed above, as the isentropic efficiency was accurately predicted. Moreover, it can also be observed that the experimental distribution of the total temperature ratio had a distinct change below 40% spanwise. As a comparison, the throughflow result had a nearly flat distribution curve. This was mainly due to a uniformly distributed slip factor along the whole blade span in the throughflow simulation in this study. A more realistic distributing method of the slip factor suiting the throughflow calculation is still expected.

4.2. High-Performance Centrifugal Compressor (HPCC)

The HPCC is an Allison Engine Company design with a design mass flow rate of 4.39 kg/s, a design pressure ratio of 4, and a design revolution speed of 21,789 RPM. This centrifugal compressor originates from a geometric scaling up to the DDA’s 404-III compressor, making use of the backward impeller and the truncated-impeller splitter. Detailed experimental measurements on its flow field and off-design performances have been determined by Skoch and Prahst et al. [55]. Compared with the LSCC case, the HPCC case apparently has a much higher design aerodynamic loading and can operate at the choke conditions if the boundary conditions are satisfied. This is feasible for the validation of the choking-prediction capability for the developed throughflow model in the centrifugal flow compressor simulation. Another challenge in the throughflow simulation of the HPCC is how to simulate the effects of the splitter blades, which do not exist in the LSCC case. In a throughflow model, the impeller blade is usually represented by its camber surface adding the corresponding blockage distribution. In this study, within the meridional region of the splitter, the new camber surface was generated by a weight average between the cambers of the impeller and splitter. For simplicity, an arithmetic average was adopted in this case. Meanwhile, the blockage factor corresponding to the thickness of the splitter was added to the original blockage factor caused by the impeller. Additionally, the losses caused by the splitter were considered by adopting the concept of “equivalent blade number”, i.e., the number of splitters was converted to the equivalent number of the impeller blades according to their chord lengths. This equivalent number was then added to the impeller blade number as the input of the empirical loss models. The generated mesh for the HPCC case had 161 × 20 points and is presented with corresponding blockage distribution in Figure 7.

A typical computation for one operating condition can be converged within 1 min for this case. Convergence histories of the throughflow simulation of the HPCC at a typical operating point are shown in Figure 8 to present a comparative study of the different methods of CBC enforcement. It can be observed that with the same initialization and computational set, the developed CBCs based on Newton’s method achieved convergence after approximately 3500 iterations, whereas the Riemann-invariant-based CBCs blew up at the beginning of the computation after about 40 iterations. The reason why the Riemann–invariant-based CBCs presented in [51] experience such difficulty in convergence is that this CBC method explicitly utilizes the velocity of the internal domain. When the initial value of the internal domain is not close to the real inlet boundary value, convergence difficulty occurs. In contrast, the CBCs based on Newton’s method developed in this paper implicitly utilize the velocity at the inlet boundary, which preserves the computational accuracy at the inlet boundary. This result proves the superiority of the developed CBCs based on Newton’s method in terms of computational robustness, at least for the throughflow simulation of centrifugal compressors with a high pressure ratio.

The predicted characteristics of the total pressure ratio and the isentropic efficiency by the developed throughflow model were compared with the experimental data and 3D-CFD results at different speed lines as shown in Figure 9a,b, respectively. Although overestimated at the design speed, the total pressure ratio characteristic computed by the throughflow solver matched well with the experimental data and the 3D-CFD results in the whole operating range. At the design mass flow rate, the throughflow results showed relative errors of 8.3% and 4.9% compared to the experimental data and 3D-CFD results, respectively. This should be attributed to the underestimated slip effects as the isentropic efficiency predicted by the throughflow solver was in good agreement with the experiment at the design mass flow rate. At the off-design speed, the throughflow results showed even better accuracy compared with the experimental data and the 3D-CFD results. Regarding the isentropic efficiency characteristic, the throughflow and 3D-CFD results showed good agreement at all the different speed lines. At the design speed, both calculations matched well with the experimental data, whereas they underestimated the efficiency at lower speeds. It can also be seen from Figure 9 that the choke mass flow rate was accurately predicted at the design speed by the throughflow solver with a value of 5.07 kg/s, which was almost the same as the experimental data and 3D-CFD results. For the near-stall conditions, it was observed that the throughflow calculation could not reach the stall boundary shown in the experimental characteristics, especially for the 60%~80% speed. This problem may have been caused by the setup of the loss correlations or the boundary conditions and will be further studied in later work.

Table 4. Contributions of different losses for the HPCC case at 100% speed.

Mass Flow Rate (kg/s)	Incidence Loss (%)	Blade Loading Loss (%)	Skin Friction Loss (%)	Clearance Loss (%)	Mixing Loss (%)	Disc Friction Loss (%)	Recirculation Loss (%)	Vaneless Diffuser Loss (%)
3.60	2.210	8.287	22.541	2.541	4.751	7.348	33.646	18.674
4.39	0.313	11.346	26.448	3.756	9.624	8.842	18.936	20.736
5.07	0.090	11.622	32.973	4.595	16.216	9.009	6.757	18.739

further reports the detailed contributions of each loss to the total loss in the throughflow simulation of the HPCC at different operating conditions.

The meridional flow fields obtained by the throughflow and 3D-CFD simulations are compared at the design mass flow rate in Figure 10 in terms of the relative Mach number. In the comparison, the 3D solution was circumferentially averaged. It can obviously be seen that both calculations predicted similar subsonic flow fields with a maximum relative Mach number of about 0.8 at the leading edge of the blade tip. The blockage effects resulting from the blade thickness of the splitter were successfully captured by the throughflow solver as it was predicted by the 3D solution, which locally exhibited a positive gradient of the relative Mach number at the leading edge of the splitter. In addition, the low Mach number region at the impeller exit near the shroud was predicted by both different calculation methods despite the slight discrepancy in its strength. This kind of low-energy flow means that the blade leakage effects were successfully simulated by the present throughflow solver.

5. Conclusions and Discussion

A time-marching throughflow method for the centrifugal compressor simulation was presented. To improve the robustness of the throughflow calculation, a method of enforcing the characteristic boundary conditions based on Newton’s iteration method was developed and introduced. The flow deflection inside the impeller passage was achieved by the well-validated inviscid blade force model. The viscous effects were simulated by the classic viscous force model coupling with the complete empirical loss and slip model for the centrifugal compressor. This throughflow solver was finally employed to analyze the performance and flow fields of the LSCC and HPCC cases as validations. The results were compared with the experimental data and 3D calculations. The main conclusions are summarized below.

• Compared to the classical Riemann-invariant-based method, the developed enforcement method for the boundary conditions at the inlet and outlet of computation models significantly improved the robustness of the time-marching throughflow simulation of centrifugal compressors. For the high-load test case, i.e., the HPCC case, the latter obtained converged throughflow results while the former encountered blowing up during the computation.
• Concerning the computational accuracy, the developed throughflow method accurately predicted the overall performance of the centrifugal compressor at design conditions, and also obtained correct off-design characteristics that were close to the experimental or the 3D-CFD results. Reasonable results were also obtained with respect to the spanwise distribution of the flow parameters and even the whole meridional flow field. Particularly, the throughflow method precisely predicted the choke mass flow rate automatically and precisely for the high-load centrifugal compressor.
• The simulation of the LSCC and HPCC cases using the developed throughflow method only cost less than two minutes, which is approximately two orders of magnitude lower than the 3D-CFD method.

In conclusion, the throughflow method developed in this study provides an accurate and efficient tool for the performance analysis of centrifugal compressors. Although some kinds of empirical models, such as the recirculation and the slip models, are expected to be further improved, the prediction results demonstrated that this throughflow method is of great value to the preliminary design of centrifugal compressors.