# A CFD-Based Throughflow Solver Using Cascade Potential Theory for Axial Compressor Flow Modeling

Florent Kravtsoff, Michaël Schvallinger, Xavier Ottavy

► To cite this version:

Florent Kravtsoff, Michaël Schvallinger, Xavier Ottavy. A CFD-Based Throughflow Solver Using Cascade Potential Theory for Axial Compressor Flow Modeling. ASME Turbo Expo 2023, Jun 2023, Boston, United States. 10.1115/GT2023-102875. hal-04305803

HAL Id: hal-04305803

https://hal.science/hal-04305803v1

Submitted on 27 Nov 2023

HAL is a multi-disciplinary open access
archive for the deposit and dissemination of sci-
entific research documents, whether they are pub-
lished or not. The documents may come from
teaching and research institutions in France or
abroad, or from public or private research centers.

L'archive ouverte pluridisciplinaire HAL, est
destinée au dépôt et à la diffusion de documents
scientifiques de niveau recherche, publiés ou non,
émanant des établissements d'enseignement et de
recherche français ou étrangers, des laboratoires
publics ou privés.

<figure><img src="image_1.png" /></figure>

Distributed under a Creative Commons CC BY 4.0 - Attribution - International License
# A CFD-BASED THROUGHFLOW SOLVER USING CASCADE POTENTIAL THEORY FOR AXIAL COMPRESSOR FLOW MODELING

Florent Kravtsoff¹*, Michël Schvallinger², Xavier Ottavy³*,

¹Univ Lyon, École Centrale de Lyon, CNRS, Univ Claude Bernard Lyon 1, INSA Lyon, LMFA, UMR5509, 69130, Écully, France

²Safran Aircraft Engines, 77550, Moissy-Cramayel, France

³Univ Lyon, CNRS, École Centrale de Lyon, INSA Lyon, Univ Claude Bernard Lyon 1, LMFA, UMR5509, 69130, Écully, France

## ABSTRACT

This paper describes the use of the Hess & Smith cascade potential method as an additional feature for throughflow modeling. This method is incorporated inside ASTEC, a CFD-based throughflow solver developed by Safran Tech providing local blade loading information for blade force and loss modeling. The test cases are all subsonic, composed of two compressors cascades and a low speed compressor. The two profile families for the cascades are a classical NACA 65 009 and a controlled diffusion cascade designed by the Whittle Laboratory. The compressor is the well studied CME2, a research single stage designed by Safran Aircraft Engines. From the potential calculation, boundary layers and wakes are calculated on the blades providing deviation angles and profile losses. The tip leakage losses are also estimated from blade loading at the tip of a rotor. The results are compared to MISES code predictions, RANS simulations, experimental data and classic Lieblein empirical correlations. This work aims to model compressors flow with geometry-based classic analytical methodologies for a minimal use of empiricism. The predictions are good on bi-dimensional cascades cases compared to several other methodologies and at mid-span for the multiple incidences considered compared to experimental cascade results. For the compressor case, the blade loading near mid-span of the rotor blades show good agreement compared to RANS results leading to an overall good capture of the performance characteristic trends near design point. There is still room for improvements near the end-walls, mostly due to the lack of hub boundary layer model in the actual methodology. These results still build a step on which the enrichment of the flow model can be pursued.

## NOMENCLATURE

### Roman letters

b Blockage factor [-]

$C_f$ Friction coefficient in boundary layer theory [-]

$C_E$ Head entrainment coef. in boundary layer theory [-]

$C_d$ Discharge coef. in vena contracta theory [-]

$C_D$ Profile drag coefficient [-]

$C_p$ Profile pressure coefficient [-]

$d$ Normalized wake distance [-]

$e_t$ Total specific energy [J kg⁻¹]

$\mathbf{f}_b$ Blade force volumic source term [N m⁻³]

$\mathbf{f}_{bi}$ Inviscid blade force volumic source term [N m⁻³]

$\mathbf{f}_{bv}$ Viscous blade force volumic source term [N m⁻³]

$\bar{\mathbf{I}}$ Identity tensor [-]

$h_t$ Total specific enthalpy [J kg⁻¹]

$H_{12}$ Shape factor in boundary layer theory [-]

$N$ Number of blades in a blade row [-]

$\mathbf{q}$ Heat fluxes [J m⁻² s⁻¹]

$Q$ Source intensity [m² s⁻¹]

$\Delta s$ Specific entropy increment [J K⁻¹ kg⁻¹]

$T$ Static temperature [K]

$R_e$ Reynolds number [-]

$U_e$ Free stream velocity [m s⁻¹]

$\mathbf{V}$ Velocity vector in absolute frame of reference [m s⁻¹]

$\mathbf{W}$ Velocity vector in relative frame of reference [m s⁻¹]

$w$ Normalized wake width [-]

$y^+$ Dimensionless wall distance [-]

### Greek letters

$\alpha$ Absolute angle [rad]

$\beta$ Blade relative angle [rad]

$\delta$ Boundary layer thickness [m]

$\delta_1$ Boundary layer displacement thickness [m]

$\delta_2$ Boundary layer momentum thickness [m]

$\varepsilon$ Blade circumferential thickness [m]

$\varepsilon_1$ Aerodynamic circumferential thickness [m]

$\rho$ Density [kg m⁻³]

$\sigma$ Lineic source density [m s⁻¹]

$\gamma$ Lineic vortex density [m s⁻¹]

$\bar{\tau}$ Constraint viscous stress tensor [Pa]

$\Omega$ Angular velocity vector [rad s⁻¹]

### Superscripts

$\cdot$ Tensor notation

$\cdot k$ Variable at solver k iteration

<footnote>*Corresponding author: florent.kravtsoff@ec-lyon.fr, xavier.ottavy@ec-lyon.fr</footnote>
• Perturbation variable

• * Complex conjugate quantity

### Subscripts

• *t* Total quantity

• *te* Trailing edge quantity

• *ps* Pressure side quantity

• *ss* Suction side quantity

• *w* Wake quantity

• *∞* Infinite upstream quantity

### Geometry

x, r, θ Cylindrical coordinates

z Complex number in 2D plane

c Blade chord [m]

c_l Blade camber-line [m]

h/H Reduced blade height [-]

g Tip gap [m]

m Meridional coordinate [m]

m' Normalized meridional coordinate [-]

s Blade pitch [m]

n Camber-line normal vector for inviscid blade force

n_p Profile normal vector

t_p Profile tangent vector

### Operators

∇· Gradient operator

div· Divergence operator

• *t* Transpose operator

# 1. INTRODUCTION

The use of throughflow methods in aeronautical compressors design is historically strongly established. However, the past few decades demonstrated a growing interest for high fidelity methods as the computational costs decreased. This led to two phenomena: on one hand, the research interest for throughflow modeling declined but on the other hand understanding of 3D complex features of compressor flows has strongly increased. In compressor design actual trends, a reduced number of stages is observed. This involves increased blade loading as already observed by Cumpsty [1] and the generation of secondary flows of greater importance. To control these undesirable flows, 3D complex shapes of blades have emerged thanks to better understanding of compressor flow topologies. Although throughflow methodologies were lacking popularity compared to more advanced high fidelity simulations they are still an essential design tool in the process of making a new engine [2]. In this context, a recent interest in CFD-based throughflow solvers is observed [3–5]. Thanks to increased CFD power, the computational time required for this type of simulations is now of the same order of the well-known streamline curvature method (SCM) during its advent [6]. CFD-based method does not suffer from some drawbacks imputed to SCM method such as the handling of recirculations or compressible effects (although the quality of shock capturing properties of CFD-based throughflow must be treated with caution [7]). Also, it is easier to consider a multi-fidelity conception approach when the throughflow code and the higher fidelity methods share the same CFD solver. Finally, keeping the throughflow models up-to-date is easier in CFD-based methodologies, especially if the throughflow routines are external to the CFD solver. If the recent

better understanding of 3D flow features is admitted, the reduction of this knowledge to reliable reduced order model adaptable to throughflow simulations is not yet systemically observed as noticed by Banjac et al. [8]. In addition to this, the legacy models used in the past show limited accuracy for modern designs [9]. This is why models must continuously be developed to improve new generation of throughflow codes suitable for modern compressor design methodologies.

To address this objective, this paper present the use of literature models showing a low level of empiricism and using the local blade geometry. They are incorporated in ASTEC (Aerodynamic Source Term for Efficient Computation), a CFD-based throughflow solver relying on elsA CFD solver. In details, the use of the well-known Hess & Smith cascade potential method [10] allows estimating the local blade pressure coefficients which are used to compute profile boundary layers using an integral boundary layer method. This leads to the estimation of profile drag losses. The deviation angles resulting from the development of the profile boundary layers are estimated using a simple coupling between the potential method and the boundary layer method by the modification of the profile geometry. In addition, wakes are also estimated as an axial evolution of the state of the boundary layer at the trailing edge. Finally, in the case of a rotor, tip leakage flow losses are computed from the pressure distribution at the blade tip. The results are confronted to MISES, RANS, experimental and empirical correlation data.

# 2. ASTEC: A CFD-BASED THROUGHFLOW SOLVER

ASTEC refers to a newly developed<fnref n="1" /> throughflow code which solves the circumferentially averaged Navier-Stokes equations. The circumferentially averaged flow is calculated with the state-of-the-art finite volumes CFD solver elsA [11] developed by ONERA and the volumic source terms arising from the average process in blades rows are calculated with ASTEC external modules. This non-invasive coupling between elsA and the external source terms occurs at each elsA iterations. Similar coupling exists in the literature such as the work done by Thollet et al. [12] to study air intake distortions. The resulting solver is fully parallelized. It uses one processor per blade zone and one processor per inter-blade zone making the calculation response time independent of the stage number. For commodity and simplicity ASTEC will refer to the ensemble elsA + source terms in the following.

## 2.1 General governing equations

CFD-based throughflow solvers intend to resolve the circumferentially-average flow. The rigorous establishment of the circumferentially averaged 3D-Navier-Stokes equations is well described by Simon [3]. It consists in the successive application of Reynolds ensemble-average, a time-average, the average-passage and finally the circumferential-average as detailed by Adamczyk [13]. Each successive average process will create additional terms. Those terms, Reynolds, unsteady, aperiodic and circumferential stresses, added to the stationary axisymmetric

<footnote n="1">¹ASTEC was recently developed by Safran Tech the Research & Technology center of Safran Group.</footnote>
flow field, allow the reconstruction of the original 3D-unsteady turbulent flow. The general governing throughflow equations written in conservative form following these average process read:

$$ \frac{\partial b\rho}{\partial t} + \text{div } b\rho\mathbf{V} = 0 \qquad (1) $$

$$ \frac{\partial b\rho\mathbf{V}}{\partial t} + \text{div } b(\rho\mathbf{V}\mathbf{V}^t - \bar{\tau} + p\bar{\mathbf{I}}) = b\mathbf{S}_{\rho\mathbf{V}} \qquad (2) $$

$$ \frac{\partial b\rho e_t}{\partial t} + \text{div } b(\rho h_t \mathbf{V} - \bar{\tau} \cdot \mathbf{V} + \mathbf{q}) = bS_{e_t} \qquad (3) $$

In those equations one can notice the blockage factor b appearing in the average process. Its value is less than one and equals one in the perfect case of no existent aerodynamic blockage outside a blade row. In the left hand-side, one can recognize the Navier-Stokes conservation equations omitting b. In the right hand-side two composite source terms appears, $S_{\rho\mathbf{V}}$ and $S_{e_t}$. They are composed of stresses arising from the average process and, in presence of blades, also blade forces, contributing in momentum conservation Eq. (2) and total energy conservation Eq. (3). For the sake of readability, although nothing denote the consecutive averaging operations in the above system, the variables must be understood as the multiple-average resultant version of themselves obtained at the end of the whole averaging process.<fnref n="2" />

The left-hand side can be solved with a CFD solver and the right-hand side must be estimated. In the special case of no further simplifications, the estimation of the right-hand side can be tedious since it implies the precise estimation of all the stresses and blade terms obtained in the average process. Such a task has been partly addressed by Thomas et al. [14] with the use of the circumferential stresses from 3D Navier-Stokes simulations, showing gain in results precision compared to the case without the extracted terms. This is relevant to high order throughflow modeling topic. For the work described in this paper, a particular importance has been placed in the time-efficiency of the simulation. For the authors, what made throughflow simulation so popular in the history of turbomachinery conception tools is their robustness and fast-response abilities in multistage cases. In order to respect this guiding principle, hypothesis are made to simplify the system of equations solved as explained in the following section.

## 2.2 ASTEC governing equations

In this work, the equations solved with ASTEC are obtained after four main hypotheses:

* The Euler form of Navier-Stokes equations are used;

* Heat fluxes with the exterior of the domain are neglected;

* All stresses obtained in the averaging process part of $S_{\rho\mathbf{V}}$ and $S_{e_t}$ are neglected;

* The blade force $\mathbf{f}_b$ are taken into account through two volumic source terms $\mathbf{f}_{bi}$ and $\mathbf{f}_{bv}$ that are the inviscid and the viscous part respectively.

The previous system immediately become:

$$ \frac{\partial b\rho}{\partial t} + \text{div } b\rho\mathbf{V} = 0 \qquad (4) $$

$$ \frac{\partial b\rho\mathbf{V}}{\partial t} + \text{div } b(\rho\mathbf{V}\mathbf{V}^t + p\bar{\mathbf{I}}) = b\mathbf{f}_b \qquad (5) $$

$$ \frac{\partial b\rho e_t}{\partial t} + \text{div } b\rho h_t \mathbf{V} = b\mathbf{f}_b \cdot \mathbf{\Omega} r \qquad (6) $$

With $r$ the local radius and $\mathbf{\Omega}$ the blade rotation velocity vector around the shaft axis. The final non-invasive ASTEC coupling system is obtained using property of the divergence operator and the non-dependence in time of the blockage factor $b$:

$$ \frac{\partial \rho}{\partial t} + \text{div } \rho\mathbf{V} = -\frac{1}{b}(\rho\mathbf{V} \cdot \nabla b) \qquad (7) $$

$$ \frac{\partial \rho \mathbf{V}}{\partial t} + \text{div} (\rho \mathbf{V} \mathbf{V}^t + p \bar{\mathbf{I}}) = \mathbf{f}_{bi} + \mathbf{f}_{bv} - \frac{1}{b} (\rho \mathbf{V} \cdot \nabla b) \mathbf{V} \qquad (8) $$

$$ \frac{\partial \rho e_t}{\partial t} + \text{div } \rho h_t \mathbf{V} = (\mathbf{f}_{bi} + \mathbf{f}_{bv}) \cdot \mathbf{\Omega} r - \frac{1}{b}(\rho h_t \mathbf{V} \cdot \nabla b) \qquad (9) $$

In this system one can identify the Euler equations integrated by elsA at the left hand-side and an additional blockage source term, identified by terms with $b$ in their expression calculated by ASTEC jointly with $\mathbf{f}_{bi}$ and $\mathbf{f}_{bv}$.

## 2.3 Baseline source term models

**Blockage model.** The blockage factor involve blade circumferential thickness $\varepsilon$ visible in Fig. 1 and blade number $N$.

$$ b = 1 - \frac{N\varepsilon}{2\pi r} \qquad (10) $$

This formulation is consistent with the circumferential averaging process. However, for thick or highly staggered blades the blockage factor may reach values that prevent the correct estimation of compressor choke massflow rate for high Mach number cases. In these cases the use of a normal to flow direction blockage factor may allow better predictions as studied by Baralon et al. [7]. The blockage factor may also be increased with aerodynamic blockage as it will be seen in Sect. 5.1. For this work, the formulation in Eq. (10) is used as the test cases presented in the results sections belong to subsonic flow domain.

**Blade force model.** As illustrated in Fig. 1 the blade force is decomposed into the two terms $\mathbf{f}_{bi}$ and $\mathbf{f}_{bv}$ to account for the two effects of blades on the flow: flow turning and viscous loss.

The inviscid term $\mathbf{f}_{bi}$ must not create entropy since it is the part of $\mathbf{f}_{bv}$. Thus the flow turning model must deliver a force that is normal to the mean relative flow field, summarized by the orthogonality condition:

$$ \mathbf{f}_{bi} \cdot \mathbf{W} = 0 \qquad (11) $$

In practice, the flow turning procedure is an iterative process and the orthogonality condition is satisfied when the calculation is converged. More precisely the procedure assume that the flow deflection is known a priori, for example it is possible to assume

<footnote n="2">²For more details see Simon [3] Appendix B. p.173.</footnote>
that the flow follows the camber line direction so an orthogonality condition between the mean relative flow and the normal vectors **n** to the camber line is resolved to obtain the inviscid blade force modulus:

$$ |\mathbf{df}_{\mathbf{bi}}|^k = -\mathcal{H}(W_x^k n_x + W_r^k n_r + W_\theta^k n_\theta) \quad (12) $$

Where $|\mathbf{df}_{\mathbf{bi}}|^k$ is the contribution obtained at iteration $k$ to total inviscid blade force modulus $|\mathbf{f}_{\mathbf{bi}}| = \sum_k |\mathbf{df}_{\mathbf{bi}}|^k$ during ASTEC convergence process. The minus sign is due to the outward orientation of the normal vectors and the constant $\mathcal{H}$ must be interpreted as a mean to adjust the rate of response to a violation of orthogonality condition as stated by Simon [3]. The modulus of the force is then projected on the normal vectors **n** and added to the right-hand side of ASTEC conservative equations.

The estimation of $\mathbf{f}_{\mathbf{bv}}$ is still an active research topic for throughflow simulation accuracy improvement. The overall precision of a throughflow model rely on the accurate estimation of the most relevant sources of loss and deviation and their correct distribution in the $(x, r)$ meridional plane. As a pure dissipative force $\mathbf{f}_{\mathbf{bv}}$ is oriented in the opposite direction of the relative mean flow and may be estimated from entropy generation sources as suggested by Denton [15] or by estimation of a drag force from boundary layer momentum thickness calculations.

Concerning the deviation angles arising from viscous flows such as profile boundary layers, or secondary flows, they are taken into account through a rotation of the normal vectors, used in $\mathbf{f}_{\mathbf{bi}}$ calculation, of the corresponding predicted deviation angle from viscous models. The calculation then converge to an updated prescribed flow angle taking into account viscous effects and so during the complete convergence procedure.

<figure><img src="image_2.png" /><figcaption>FIGURE 1: Blade thickness and blade force components</figcaption></figure>

## 2.4 Comments on the described methodology

About the first ASTEC hypothesis described in Sect. 2.2, the use of Euler equations in place of their viscous RANS version is motivated by two complementary considerations:

* Blade rows are fictitious in the simulation, so only axisymmetric end-wall boundary layers would be predicted by the turbulence model of a RANS solver;

* It is known that the turbulent boundary layers (TBL) show 3D features at the end-wall region of the blades. Moreover, the TBL at the hub of a multistage compressor does not behave as a regular TBL at end walls [16].

Running ASTEC with a viscous set of equations would require a refined mesh at end walls, to satisfy turbulence model usage with still questionable end-wall predictions once the flow passes blades rows. The gain in precision of the results by running viscous seems not worth the time invested in heavier simulation, that is why Euler mode is preferred by the authors as it also encourages the use/development of accurate secondary flow models. For the record, the use of viscous conservatives equations in ASTEC is described in the work of Budo et al. [17] applied to a realistic modern low pressure compressor.

Recent designs of axial compressors blade showed diversified profile shapes, increased solidity due to larger chords and reduced blade aspect ratios compared to the precedent generations. This brings difficulties to estimate the correct loss levels with classical semi-empirical models mostly based on cascade test cases such as the well-known Lieblein correlations [9, 18]. This fact has motivated the testing of a more locally driven approach based on local geometry as it will be described in the next section of this paper.

To conclude this section, in the following paper, ASTEC is used in analysis mode. Design mode is theoretically possible, but it has not been addressed in this work.

## 3. POTENTIAL FLOW APPROACH FOR INVISCID BLADE FORCE

As a first step toward detailed geometry driven throughflow modeling, the estimation of a realistic inviscid blade loading is incorporated in ASTEC. This is currently done by the use of the well-known singularity method also known as panel method in the literature [19].

### 3.1 Singularity method

In this work, a 2D singularity method is used on slices of the real 3D blade geometry along the blade span as illustrated in Fig. 2. It consists in the estimation of a steady, incompressible, irrotational non-viscous potential flow by the use of singularities, particular solution of Laplace equation, placed on the blade profile. If needed, blades slices are flattened in the well-known $(m' - \theta)$ stream surface coordinates before any potential calculation as described in [20]:

$$ m' = \int \frac{dm}{r} = \int \frac{\sqrt{dx^2 + dr^2}}{r} \quad (13) $$

In this conformal transformation, profiles angles are preserved and $\theta$ is unchanged. Any calculation is then performed on the flattened profile and the results are transposed to the original slice.

<figure><img src="image_3.png" /><figcaption>FIGURE 2: Blade slice definition</figcaption></figure>
Laplace equation for the incompressible flow potential $\phi$ reads:

$$ \Delta\phi = 0 \qquad (14) $$

Under the small perturbation hypothesis, the flow potential can be separated into the contribution of the onset flow potential $\phi_\infty$ and a perturbation potential $\tilde{\phi}$ induced by the presence of the profile:

$$ \phi = \phi_{\infty} + \tilde{\phi} \qquad (15) $$

The Laplace equation being linear, $\tilde{\phi}$ still satisfies it. It is then possible to use particular solutions of Laplace equation to find the perturbation potential. Those solutions are called the singularities. The commonly used singularities are sources, vortices or doublets. In addition, the potential flow is linked to the flow speed by the relation $\mathbf{V} = \nabla\phi$ so in terms of speed Eq. (15) is written:

$$ \mathbf{V} = \mathbf{V}_{\infty} + \tilde{\mathbf{V}} \qquad (16) $$

**Singularities definition.** The presence of the profile is materialized by an ensemble of singularities, sources and vortices distributed on the profile. The profile is thus discretized in $N_{\text{elem}}$, the so-called panels. Each element is given lineic source and vortex densities. For great simplification of relations, the complex formalism is used in the following.

<figure><img src="image_4.png" /><figcaption>FIGURE 3: Profile element definition</figcaption></figure>

The complex conjugate velocity induced at a point $M(z)$ of 2D complex space by a point source of intensity $Q$ located in $z_0$ reads:

$$ v_{\sigma}^{*}(z) = \frac{Q}{2\pi(z - z_0)} \qquad (17) $$

Thus, the complex conjugate velocity induced by an element $dz_i = |dz_i|e^{j\alpha_i}$ of the profile possessing a lineic source density is obtained by integration on the length of the element $|dz_i| = dz_i e^{-j\alpha_i}$ (see Fig. 3):

$$ v_{\sigma i}^{*}(z) = \int_{z_i}^{z_{i+1}} \frac{\sigma_i e^{-j\alpha_i}}{2\pi(z - z_i)} dz_i \qquad (18) $$

Where $\alpha_i$ if the orientation angle of the element in the complex plane, $e$ the complex exponential and $j$ is the complex number defined as $j^2 = -1$.

The complex conjugate velocity induced by the lineic vortex density $\gamma_i$ on $dz_i$ is obtained by an indirect rotation of $\pi/2$ of the complex conjugate velocity induced by a lineic source density.

$$ v_{\gamma_i}^{*}(z) = \int_{z_i}^{z_{i+1}} \frac{\gamma_i e^{-j(\alpha_i + \pi/2)}}{2\pi(z - z_i)} dz_i \qquad (19) $$

For the whole set of singularities the total complex conjugate perturbation speed at $M(z)$ is the added contributions of all singularities on the profile so that:

$$ \tilde{V}^{*}(z) = \sum_{N_{\text{elem}}} v_{\sigma i}^{*}(z) + v_{\gamma i}^{*}(z) \qquad (20) $$

**Neumann boundary condition.** A non penetration condition is then applied on the surface of the profile. In practice the condition is applied at the centers $c_j$ of the elements $dz_j, 1 \le j \le N_{\text{elem}}$ using real velocity vectors, giving $N_{\text{elem}}$ equations:

$$ \forall j, \quad \mathbf{V}(c_j) \cdot \mathbf{n}_{p,j} = \mathbf{V}_{\infty} \cdot \mathbf{n}_{p,j} + \tilde{\mathbf{V}}(c_j) \cdot \mathbf{n}_{p,j} = \mathbf{0} \qquad (21) $$

$$ \Leftrightarrow \sum_{i=1}^{N_{\text{elem}}} (\mathbf{v}_{\sigma i}(c_j) + \mathbf{v}_{\gamma i}(c_j)) \cdot \mathbf{n}_{p,j} = -\mathbf{V}_{\infty} \cdot \mathbf{n}_{p,j} \qquad (22) $$

Where $\mathbf{n}_{p,j}$ is the normal at the center of element $j$ located at $z_j$.

**Kutta-Joukowski condition.** A particular flow escaping condition is applied at the trailing edge to reflect the real flow behavior. It consists in preventing the flow to go from pressure side to suction side or vice-versa at the trailing edge and so, allows lift generation by the profile. In a pressure point of view it is equivalent to say that the pressure at the suction side equals the pressure at the pressure side at the trailing edge. Using Bernoulli theorem, this condition reads in terms of velocities:

$$ \mathbf{V}(z_{te,ps}) \cdot \mathbf{t}_{te,ps} = \mathbf{V}(z_{te,ss}) \cdot \mathbf{t}_{te,ss} \qquad (23) $$

$$ \Leftrightarrow (\mathbf{V}_{\infty} + \tilde{\mathbf{V}}(z_{te,ps})) \cdot \mathbf{t}_{te,ps} = (\mathbf{V}_{\infty} + \tilde{\mathbf{V}}(z_{te,ss})) \cdot \mathbf{t}_{te,ss} \qquad (24) $$

Where $\mathbf{t}_{te,ps}$ and $\mathbf{t}_{te,ss}$ are the tangent vectors at the trailing edge of the profile at pressure side, respectively suction side, oriented from leading edge to trailing edge. In practice this condition is applied at the center of the two trailing edge elements.

The application of Kutta-Joukowski condition therefore require a sharp trailing edge and an adapted refinement process at the profile trailing edge. In practice, if the analyzed profile presents round trailing edge, it is then cut and closed conserving profile chord. For the refinement process, distribution laws based on cosine function showed good results in the literature [19]. Also, leading edge refinement is important for a good accuracy of singularity method, so the use of symmetric cosine-based points distribution on suction side and pressure side fulfill both requirements and is used in this work.

**Solved linear system.** In the singularity method, the unknown of the problem are the singularities intensities $\sigma_i$ and $\gamma_i$ that simulate the profile influence on the onset potential flow. Equations (22) and (24) total $N_{\text{elem}} + 1$ equations, but there is a total of $2N_{\text{elem}}$ unknown singularities. To solve this system an assumption is made:

$$ \forall i, \quad \gamma_i = \gamma \qquad (25) $$
Physically this means that the circulation around the profile is homogeneously distributed. The number of unknown then passes from $2N_{elem}$ to $N_{elem} + 1$. The unknown vector reads:

$$ \mathbf{X}^l = [\sigma_1 \ \sigma_2 \ \dots \ \sigma_{N_{elem}} \ \gamma] \quad (26) $$

The system solved is of the form:

$$ \mathbf{AX} = \mathbf{B} \quad (27) $$

With the $N_{elem}$ first lines of A and B filled with Eq. (22) and their respective last line filled with the Kutta-Joukowski relation Eq. (24) using Eqs. (18), (19) and (25). This system is easily solved with linear algebra routines. Once X is known it is possible to calculate the flow velocity everywhere around the blade and on the blade profile using Eq. (16).

## 3.2 Cascade extension: Hess & Smith method

The previous detailing of the singularity method allows to easily understand its extension to cascade geometries done by Hess & Smith [10]. The principle of the extension is to consider a cascade of as an infinite repetition of same geometry spaced apart by a distance equals to the pitch s in cascade both directions as pictured in Fig. 4.

**Singularities expression modification.** In the considered cascade the influence of repeated profiles accumulate and the complex conjugate velocity induced by the entire cascade of lineic sources on corresponding elements on a point M(z) can be written using Eq. (18):

$$ v_{\sigma_i}^*(z) = \int_{z_i}^{z_{i+1}} \left[ \frac{\sigma_i e^{-j\alpha_i}}{2\pi(z-z_i)} + \frac{\sigma_i e^{-j\alpha_i}}{2\pi(z-z_i+js)} + \frac{\sigma_i e^{-j\alpha_i}}{2\pi(z-z_i-js)} + \frac{\sigma_i e^{-j\alpha_i}}{2\pi(z-z_i+2js)} + \frac{\sigma_i e^{-j\alpha_i}}{2\pi(z-z_i-2js)} + \dots \right] dz_i \quad (28) $$

The expression in the brackets regroups to an analytic expression as an alternate series decomposition of coth(z) on $\mathbb{C}$. With the same remark made to obtain Eq. (19), expressions of Eqs. (18) and (19) become in the case of a cascade problem:

$$ v_{\sigma_i}^*(z) = \int_{z_i}^{z_{i+1}} \frac{\sigma_i e^{-j\alpha_i}}{2s} \coth\left(\pi \frac{z-z_i}{s}\right) dz_i \quad (29) $$

$$ v_{\gamma_i}^*(z) = \int_{z_i}^{z_{i+1}} \frac{\gamma_i e^{-j(\alpha_i + \pi/2)}}{2s} \coth\left(\pi \frac{z-z_i}{s}\right) dz_i \quad (30) $$

The effect on the resolved system Eq. (27) is to modify the A matrix influence coefficients.

**Implicit boundary condition for cascades.** As stated by Lewis [21], singularity method developed in Sect. 3.1 is adapted for isolated profile with inlet flow angle $\alpha_{\infty}$ set from profile chord reference. In a cascade or turbomachinery application the relative inlet velocity vector $\mathbf{W}_1$ is known and $\mathbf{W}_{\infty}$ required for the singularity method can be deduced from blades circulation.

$$ \mathbf{W}_{\infty} = \frac{\mathbf{W}_1 + \mathbf{W}_2}{2} \quad (31) $$

$$ \mathbf{W}_{\infty} = \mathbf{W}_1 - \frac{\Gamma}{2s} \mathbf{y} \quad (32) $$

Where $\Gamma$ is the circulation generated by the cascade, counted positively clockwise and y the direct unit vector in the direction of the cascade as illustrated in Fig. 4. $\Gamma$ is linked to the lineic vortex distribution on profiles by integrating on the profile length ($\Gamma$ and the $\gamma_i$ having opposite sign convention):

$$ \Gamma = - \int_{\mathcal{C}} \gamma dl \quad (33) $$

Numerically, Eq. (33) using Eq. (25) become:

$$ \Gamma = - \sum_{N_{\text{elem}}} \gamma_i |dz_i| \quad (34) $$

$$ \Gamma = -\gamma \sum_{N_{\text{elem}}} |dz_i| \quad (35) $$

Then $\mathbf{W}_{\infty}$ is substituted to $\mathbf{V}_{\infty}$ using Eq. (32) and Eq. (35). The system Eq. (27) is thus modified accordingly: the last column of A associated with $\gamma$ unknown are updated as well as the last line of A associated with the Kutta-Joukowksi condition.

<figure><img src="image_5.png" /><figcaption>FIGURE 4: Cascade definition</figcaption></figure>

## 3.3 Prediction of inviscid blade force

The application of the above procedure make possible the estimation of potential flow relative velocity on the profile and in the ensemble of stream surfaces defined by $(x, r)$ slices. Thus, incompressible pressure coefficients $C_p = 1 - \mathbf{W}^2/\mathbf{W}_1^2$ on the blades are available and also, relative mean flow angles in blade passage since it is possible to analyze the potential velocity everywhere. In particular mean flow angles are used to set the rotation of blade normal vectors $\mathbf{n}$ in Eq. (12) for inviscid blade force source term $\mathbf{f}_{bi}$ and pressure gradients on the profiles can be used for development of viscous effects for $\mathbf{f}_{bv}$ source term as detailed in the next section.

## 4. VISCOUS FLOW MODELING

Using previously introduced potential methodology, particular identified loss sources are modeled in ASTEC: profile loss and deviation, wake losses and tip leakage flow losses.
## 4.1 Profile loss and deviation

For boundary layers growing on the surface of the blade, an integral boundary layer method for turbulent boundary layers on flat plate with pressure gradient is used. It consists in the space marching integration on pressure side and suction side of the Von Kármán boundary layer integral equation using Head turbulent entrainment equation [22]:

$$ \text{Von Kármán:} \quad \frac{d\delta_2}{dx} = \frac{C_f}{2} - (H_{12} + 2)\frac{\delta_2}{U_e}\frac{dU_e}{dx} \quad (36) $$

$$ \text{Head:} \quad \frac{d(\delta - \delta_1)}{dx} = C_E - \frac{\delta - \delta_1}{U_e}\frac{dU_e}{dx} \quad (37) $$

To close this system the empirical friction and entrainment coefficients of Ludwieg-Tillman [23] and Head [22] respectively, are provided:

$$ C_f = 0.24610^{-0.678H_{12}} Re_{\delta_2}^{-0.268} \quad (38) $$

$$ C_E = 0.0306(H^* - 3)^{-0.653} \quad (39) $$

$$ \text{where:} \quad H^* = \frac{\delta - \delta_1}{\delta_2} = 1.535(H_{12} - 0.7)^{-2.715} + 3.3 \quad (40) $$

Where $U_e$ is the local free stream velocity above the blade surface. It is thus provided by the potential method in this approach. $\delta_1$ and $\delta_2$ are respectively the displacement and momentum boundary layer thicknesses. $H_{12}$ is the boundary layer shape factor defined as the ratio of $\delta_1$ over $\delta_2$. $Re_{\delta_2}$ in Eq. (38) is the Reynolds number based on $\delta_2$.

**Stall.** The previous introduced boundary layer method does not allow stalled regime which appears when $H_{12} \ge 2.4$. A simple procedure is adopted to avoid the method to fail in this case. In a physical point of view, when stall appears on the suction side of a cascade profile for example, the profile is no longer able to provide deflection to the flow. The following assumption is thus made: in case of stall, the flow on the suction side will follow the streamline given by the profile enlarged with displacement thickness, which is then extrapolated keeping the same direction from the point of detected stall to the trailing edge. A stalled enlarged profile is thus obtained and $\delta_1$ in stalled region can be deduced geometrically (see Fig. 5). Similarly, $H_{12}$ is linearly extrapolated on the blade surface from the point of detected stall until the profile trailing edge allowing to deduce $\delta_2$ in stalled region. In practice, the stall threshold value of $H_{12}$ is fixed to 2.0 in order to perform extrapolation before the boundary layer predicted by the integral method enters exponential growth regime.

**Profile loss.** The profile losses are then calculated from the well-known Von Kármán relation between momentum thickness and profile drag at a station x:

$$ D(x) = (\rho U_e^2 \delta_2)(x) \quad (41) $$

Equation (41) is then discretized to locally update ASTEC f\_bv source term for each blade slices.

<figure><img src="image_6.png" /><figcaption>FIGURE 5: View of trailing edge region of a profile showing enlarged profile with $\delta_1$ and new camber line definition for deviation in a stalled case</figcaption></figure>

**Deviation.** The deviation induced by the growth of profile boundary layers is obtained by profile geometry modification in the potential calculation. From the estimation of profile displacement thicknesses on pressure and suction sides, an enlarged profile is built from the original one. A new camber line of the profile is then deduced from the mean line of the enlarged profile as shown in Fig. 5. The original profile thickness law is then applied to this new camber line and the potential calculation take into account the new geometry for the next ASTEC iterations until the geometry is modified again. The calculation of the mean angle given by the potential solution, obtained by averaging the velocity in the inter blade channel, is then used to rotate the blade normal vectors n in Eq. (12) for f\_bi source term. By doing so, a blade force including deviation from profile boundary layers is then calculated for each blade slices.

## 4.2 Wake evolution

Since the boundary layer characteristics are known until the trailing edge, it is possible to use a wake decay model. To do so the decay correlations of Lakshminarayana et al. [24] are used. They consist in the estimation of the decay of the velocity defect at the center of the wake after the trailing edge of the profile:

$$ \text{Near Wake:} \quad \frac{v_c}{U_e \sqrt{C_D}} = \frac{0.168}{(d + 0.12)^{0.22}} + \frac{0.353}{(d + 0.12)^{1.22}} \quad (42) $$

$$ \text{Far Wake:} \quad \frac{v_c}{U_e \sqrt{C_D}} = \frac{0.9792}{d + 0.688} \quad (43) $$

Where $U_e$ is the free stream flow velocity, $v_c$ is the velocity deficit compared to free stream velocity at the center of the wake, d is the chord normalized distance in the direction of the wake and $C_D$ is the profile drag coefficient based on cascade mean velocity ($W_1 + W_2$)/2 as in Fig. 4. The frontier between near and far wake is located at d = 0.4. The paper indicate that the correlations are not valid at the vicinity of the trailing edge for $0 < d < 0.06$. In this application the correlations are simply extrapolated for this domain.

In their work, Lakshminarayana et al. [24] correlated the velocity profile in the wake as a Gaussian function. The correlation
use two length scale $L_{ss}$ and $L_{ps}$, the distances at suction side and pressure side, respectively, where the speed deficit $v = v_c/2$. The expression of speed deficit as a function of the normalized distance normal to wake center line $\eta = l/L$ is:

$$ \frac{v}{v_c} = e^{-\ln 2 \eta^2} \qquad (44) $$

Assuming a total wake width as $2L_{ss} + 2L_{ps}$ the paper also defines the pitch normalized width of the wake $w = 2(L_{ss} + L_{ps})/s$:

$$ \text{Near Wake:} \qquad \frac{w}{\sqrt{C_D}} = 1.39(d + 0.55)^{1.22} \qquad (45) $$

$$ \text{Far Wake:} \qquad \frac{w}{\sqrt{C_D}} = 1.18(d + 0.688) \qquad (46) $$

By building an equivalent symmetric wake of characteristic length scale $L_{eq} = \frac{L_{ps}+L_{ss}}{2}$ it is possible to integrate Eq. (44) on wake width using Eqs. (45) and (46) and Eqs. (42) and (43) to obtain wake displacement thickness $\delta_{1,w}$. As shown in [25], $U_e$ decreases when the flow leaves the trailing edge of a blade row. To model this, the evolution of $U_{e,w}$ in the wake is obtained in a small iterative procedure involving basic mass conservation equation$^3$. The procedure converge to find $\delta_{1,w}$ and $U_{e,w}$ satisfying Eq. (47) in the entire wake. As the wake displacement thickness rapidly decrease (i.e., the aerodynamic blockage induced by the wake), $U_e$ rapidly decrease causing increase of wake momentum thickness $\delta_{2,w}$ according to Eq. (36) applied to wake with $C_f = 0$.

$$ U_{e,te}(s - \delta_{1,te}) = U_{e,w}(s - \delta_{1,w}) = \text{const.} \qquad (47) $$

To compute loss arising from the wake decay, the correlation for the shape factor of [25] is used:

$$ 1 - \frac{1}{H_{12,w}} = \left(1 - \frac{1}{H_{12,te}}\right) \left(40\frac{d}{c} + 1\right)^{-0.5} \qquad (48) $$

$$ \text{where:} \qquad H_{12,te} = \frac{H_{12,ss} \delta_{2,ss} + H_{12,ps} \delta_{2,ps}}{\delta_{2,ss} + \delta_{2,ps}} \bigg|_{te} \qquad (49) $$

Since $H_{12,w}$, $\delta_{1,w}$ and $U_{e,w}$ evolutions are known, $\delta_{2,w}$ can be estimated using Eq. (36) applied to wake and loss can be calculated in wakes using Eq. (41), allowing to set ASTEC $f_{bv}$ source term in wake regions for each blade slices.

### 4.3 Tip leakage flow

The tip leakage losses are incorporated in ASTEC with the well-known model of Denton [15] since the velocity on the blade tip profile is known from potential calculation. It consists first in the estimation of the leakage massflow rate crossing over the blade tip and then the estimation of entropy created by the mixing of leakage flow arriving at the other side of the blade. In case of a compressor with a leakage flow blowing from pressure side to suction side the specific entropy created $\Delta s$ for the whole rotor row is modeled with:

$$ T \Delta s = \frac{1}{m_m} \int_c V_{ss}^2 \left( 1 - \frac{V_{ps}}{V_{ss}} \right) dm \qquad (50) $$

Where $dm$ is the local leakage massflow rate estimated by the pressure difference across the blade and $m_m$ is the massflow rate of the entire blade channel. The integration of the local leakage massflow rate of Eq. (51) take place on the camber-line of the blade tip profile in the proposed implementation.

$$ dm = C_d g \sqrt{\rho 2 \Delta P} dc_l \qquad (51) $$

Where $C_d$ is the vena contracta discharge coefficient. In this work the theoretical value of Moore & Tilton [26] is used ($C_d \approx 0.84$). With an incompressible assumption using Bernoulli equation: $\Delta P = 0.5\rho(V_{ss}^2 - V_{ps}^2)$. The volumic viscous force, in opposite flow direction, obtained from an increase of specific entropy $\Delta s$ on a streamline reads [15]:

$$ f_v = \rho T \frac{\Delta s}{\Delta l} \qquad (52) $$

Where $l$ is the distance along the streamline. Equation (50) is then discretized to locally update $f_{bv}$ source term for rotor using Eqs. (51) and (52).

### 5. AERODYNAMIC BLOCKAGE COUPLING

As introduced in Sect. 4.1 the effect of the boundary layer on the blade to blade passage averaged flow angle is taken into account by modification of the shape of the profile in the potential method. Still, the influence of the profile boundary layer on the profile potential velocity distribution is not taken into account by this mean. The purpose of the following section is to model that effect. In the literature accurate techniques exist for viscid-inviscid interactions implying a strong coupling between the boundary layer and the potential solver [27]. However, applying this kind of techniques on multiple slices for multistage application seems not realistic in terms of computation time for a throughflow solver. A simple procedure is described to simulate the potential-boundary layer interaction in the following.

#### 5.1 Aerodynamic blockage

The boundary layers that develop on the pressure side and suction side of the profiles create additional blockage in the bladed region as well as secondary flows. This aerodynamic blockage is directly linked to the displacement thickness induced by viscous effects in the blade channel. For now, in ASTEC, aerodynamic blockage is produced by the profile boundary layer and wakes displacement thicknesses. The blockage factor $b$ in Eq. (10) is updated as follows:

$$ b = 1 - \frac{N(\varepsilon + \varepsilon_1)}{2\pi r} \qquad (53) $$

Where $\varepsilon_1$ is the additional thickness in the circumferential direction due to viscous effects displacement thicknesses.

#### 5.2 Viscid-inviscid coupling

The consequence in ASTEC of the use of Eq. (53) is an acceleration of the computed mean flow in presence of aerodynamic blockage compared to the inviscid case. In order to adjust potential velocities on profile according to the development of viscous effect, a comparison is made between the mean meridional velocity in bladed region given by ASTEC $V_m$ and the mean potential

<footnote>$^3$Neglecting changes in hub or casing radius.</footnote>
meridional velocity obtained in the blade to blade plane $V_m^p$ from each stream surfaces. A ratio coefficient, $V_m/V_m^p$, is obtained for each slice. This ratio is thus $\ge 1$. Then the hypothesis of a homogeneous increase of speed in the blade to blade plane thanks to aerodynamic blockage is made to modify the velocity on suction sides respectively pressure side in potential calculations. This coupling update velocities obtained by the potential procedure as a pre-treatment according to Eqs. (54) and (55) before the calculation of the profile boundary layer. As an example of interaction, starting with an already thick boundary layer, the coupling will reduce the potential pressure gradients for the next boundary layer calculation. Then, the following calculated boundary layer become thinner, reducing aerodynamic blockage and increasing pressure gradient for next boundary layer calculation and so goes on.

$$V_{ss}^{\text{updated}} = \frac{V_m}{V_m^p} V_{ss}^{\text{old}} \qquad (54)$$

$$V_{ps}^{\text{updated}} = \frac{V_m}{V_m^p} V_{ps}^{\text{old}} \qquad (55)$$

The procedure appears to converge in few iterations and the coupling is done simultaneously with the convergence of ASTEC.

# 6. RESULTS AND DISCUSSION

## 6.1 Linear cascades

The first validations where made on two cases of linear cascades: a NACA 65 009 and a controlled diffusion arc blade (CDA), both tested at LMFA by Zambonini [28] and Dawkins [29] respectively. They have been mounted on the same test rig delivering an inlet Mach number of 0.12. The shape of the blades are visible in Fig. 6 and more details on the test cases are available in the previously mentioned papers. It is important to note that the NACA cascade is equipped with trips at the leading edge allowing the development of a fully turbulent boundary layer whereas the other cascade has no trip. Thus, a boundary layer transition appear on the latter at $\approx 30\%$ of blade chord as shown in [29].

<figure><img src="image_7.png" /><figcaption>FIGURE 6: Blade passages of validation cascades</figcaption></figure>

**Bi-dimensional results.** First, results from ASTEC are compared to the well-known turbomachinery design tool MISES developed by M. Drela and H. Youngren.<fnref n="4" /> The latter uses an Euler solver coupled with viscous boundary layer calculation allowing boundary layer transition. The simulated case is a profile

from the NACA 65 009 cascade at inlet Mach number of 0.12. The Reynolds number is around $385 \times 10^3$. Here, to be comparable with the fully turbulent hypothesis of the boundary layer method of ASTEC, the transition is forced at the leading edge in MISES simulations.

The blade pressure coefficients $C_p$ for two incidences, $i = 0^\circ$ and $i = 4^\circ$ are compared on Fig. 7. The agreement between the two methods is excellent. At the trailing edge small differences still exists, due to the sharp trailing edge required by the Hess & Smith method. The results validate the implemented potential method and the viscous coupling with the integral boundary layer method of Sect. 5.2.

<figure><img src="image_8.png" /><figcaption>FIGURE 7: ASTEC vs MISES $C_p$ on NACA 65 009 profile for $i = 0^\circ$ and $i = 4^\circ$</figcaption></figure>

The blockage factor evolution in the blade and wake regions is also showed in Fig. 8. It includes aerodynamic blockage from the boundary layer calculation in both methodologies. In the bladed region, the agreement is matching, which means that both boundary layers calculated from ASTEC and MISES show very close displacement thickness evolution. However, in the wake region, the blockage in ASTEC shows a steeper decrease close to the trailing edge, revealing differences between used wake models of Lakshminarayana and MISES wake handling methods. Since the aerodynamic blockage decays more quickly in ASTEC to a lower level than in MISES, the loss level in the wake can be expected to be slightly higher in ASTEC.

<figure><img src="image_9.png" /><figcaption>FIGURE 8: ASTEC vs MISES blockage factor Eq. (10) on NACA 65 009 profile for $i = 0^\circ$ and $i = 4^\circ$</figcaption></figure>

The total pressure loss coefficient $\omega$ and the exit flow angle $\beta_2$ for this bi-dimensional case are shown in Table 1, again for $i = 0^\circ$ and $i = 4^\circ$. Data from ASTEC, MISES and bi-dimensional blade-to-blade steady RANS simulations are compared together with estimations from empirical Lieblein correlations, fitted by Aungier [30]. The RANS simulations are carried out with elsA using 2-equations $k-l$ turbulence model of Smith [31] with a Reynolds number around $385 \times 10^3$ on a $y^+ \approx 1$ mesh. The results exposed show an overall good agreement between ASTEC and other methods on a bi-dimensional case for both incidences.

<footnote n="4">⁴See A User’s Guide to MISES 2.63, http://web.mit.edu/drela/Public/web/mises/</footnote>
In details, the loss level and exit angle predicted by ASTEC agree particularly well with RANS and MISES results. Regarding the results of the Lieblein correlations, they tend to predict lower deviation angles and lower loss levels than other methods. These results validate the loss and deviation models, previously introduced in ASTEC, based on boundary layer integral method, wake decay and the coupling with the potential calculation.

<table><caption>TABLE 1: NACA 65 009 profile performances comparison $0.5c_{ax}$ after trailing edge for $i = 0°$ and $i = 4°$</caption><thead><tr><th></th><th colspan="2">Lieblein</th><th colspan="2">ASTEC</th><th colspan="2">MISES</th><th colspan="2">RANS</th></tr><tr><th><i>i</i></th><th>0°</th><th>4°</th><th>0°</th><th>4°</th><th>0°</th><th>4°</th><th>0°</th><th>4°</th></tr></thead><tbody><tr><th><i>β</i><sub>2</sub></th><td>37.8°</td><td>38.2°</td><td>38.9°</td><td>39.8°</td><td>38.7°</td><td>39.6°</td><td>38.6°</td><td>39.5°</td></tr><tr><th><i>ω</i></th><td>1.5%</td><td>2.1%</td><td>2.1%</td><td>2.4%</td><td>2.0%</td><td>2.3%</td><td>2.1%</td><td>2.4%</td></tr></tbody></table>

**Cp comparison with experiment at mid-span.** End-wall secondary flows are not yet modeled in this work, but it is necessary to take into account their induced blockage effect to validate against experimental data the presented methodology at mid-span. The choice is made to include the experimental end-wall aerodynamic blockage of the cascades in ASTEC. To do so, a span-wise distribution of pitch-wise displacement thickness is calculated from experimental data just after the trailing edge of the cascades. The resulting end-wall aerodynamic blockage contribution is linearly distributed in the axial direction from trailing edge to leading edge, with leading edge value set to 0, and then added to the local aerodynamic blockage due to the profile boundary layers. In order to determine axial distribution of $ε_1$ of Eq. (53) for each blade slices, the axial distribution of end-wall blockage is added to the boundary layer contribution to aerodynamic blockage estimated in the methodology.

<figure><img src="image_10.png" /><figcaption>FIGURE 9: Cp comparison on NACA 65 009 at mid-span at $i = 4°$ with versus without exp. end-wall blockage</figcaption></figure>

Figure 9 compares the pressure coefficient at mid-span estimated with ASTEC on NACA 65 009 case for 4° degrees of incidence with and without experimental end-wall aerodynamic blockage. Regarding those results, once the end-wall blockage is introduced, a good agreement is observed. Figure 10 summarizes the pressure coefficients for the CDA cascade at mid-span using experimental end wall blockage as well for multiple incidences also showing an excellent agreement. This again validates the simple inviscid/viscid coupling introduced in Sect. 5.2.

<figure><img src="image_11.png" /><figcaption>FIGURE 10: Cp comparison on CDA cascade at mid-span for multiple incidences</figcaption></figure>

**Performances comparison at mid-span.** The performances of the NACA 65 009 at mid-span are summarized in Table 3 for 4° of incidence where $β_2$ is the mean cascade exit angle and $ω$ the total pressure loss coefficient, both measured at mid-span. The comparison with experimental measurements are done in the wake at $x = 0.363 c_{ax}$ after the trailing edge. About the CDA cascade, the blade mean exit angle and total pressure loss coefficient are plotted against inlet flow incidence in Fig. 11. The plane of measurement is situated $0.2c$ after trailing edge in axial direction.

<table><caption>TABLE 2: NACA 65 009 mid-span performances at $i = 4°$ $0.363c_{ax}$ after trailing edge</caption><thead><tr><th></th><th>ASTEC</th><th>exp.</th></tr></thead><tbody><tr><th>β<sub>2</sub></th><td>38.9°</td><td>38.5°</td></tr><tr><th>ω</th><td>2.2%</td><td>2.3%</td></tr></tbody></table>

<table><caption>TABLE 3: NACA 65 009 mid-span performances at $i = 4°$ $0.363c_{ax}$ after trailing edge</caption><thead><tr><th></th><th>exp.</th><th>ASTEC w/ EW blockage</th><th>ASTEC w/o EW blockage</th></tr></thead><tbody><tr><th>β<sub>2</sub></th><td>38.5°</td><td>38.9°</td><td>39.8°</td></tr><tr><th>ω</th><td>2.3%</td><td>2.2%</td><td>2.4%</td></tr></tbody></table>

The agreement between the experiment and ASTEC on performance data is good on the NACA 65 009 test case. Although the pressure coefficients showed good matching in Fig. 10 for the second cascade case, the models predicted over deviation and over loss. This can be explained by the nature of the boundary layer developing on the CDA cascade. A transition exists on the suction side of the blade at ≈ 30% of chord and no transition modeling is implemented in ASTEC at the moment. Consequently, levels discrepancies on loss and deviation obviously appear. But the trends are good despite the observed offsets caused by the ignored laminar fraction of the boundary layer. However, in
the most aeronautical axial compressor applications, Reynolds numbers are high enough to make profile boundary layers transitioning close to the leading edge allowing the authors to say that the presented methodology is adapted for future applications.

<figure><img src="image_12.png" /><figcaption>FIGURE 11: CDA cascade mid-span performances against incidence measured 0.2c after trailing edge</figcaption></figure>

## 6.2 CME2 low speed compressor

The third test case is the research low speed compressor CME2. It is a single stage axial compressor with relative inlet mach number about 0.6 at blade tip at design point. More characteristics are summarized in [32]. For comparison purpose, single passage RANS steady simulations of the compressor are carried out with elsA using 2-equations k-l turbulence model of Smith [31] on a y⁺ ≈ 1 mesh. About ASTEC simulations, the mesh used is visible in Fig. 14. Each blade have approximately 60 points in the axial direction.

**CME2 characteristics.** The compressor total to total pressure coefficient and isotropic efficiency are plotted against standardized massflow rate for nominal speed in Fig. 12. The data are also compared to experimental results from Gourdain [32]. The design point is evidenced by the only available experimental point in the efficiency plot. However since the efficiency is difficult to measure in a low speed compressor test-rig, extra care must be taken for the interpretation of its position.

In terms of pressure ratio, the agreement is good near design point but deteriorates as the compressor is throttled. This can be the evidence of several phenomena occurring and not taken into account in the actual methodology yet: as the massflow rate decreases, the blade loading increases and the secondary flows intensify leading to high deviation angles at hub and casing and so reduced useful work. For secondary flows, only tip leakage loss are today taken into account. The differences can also be explained by the lack of additional losses coming from hub (and casing for the stator) boundary layers in the current version of the code. As Denton remarked [15], end-wall secondary flows remains the most difficult flow mechanism to precisely model. Nonetheless, the integral boundary layer method coupled with the potential cascade calculation allow capturing a correct level of flow turning near the design point in the present methodology.

<figure><img src="image_13.png" /><figcaption>FIGURE 12: CME2 characteristic comparison at 100Nn</figcaption></figure>

In terms of efficiency, the relative position of ASTEC and RANS results is easily explained by the same lack of end-wall boundary layer losses, but the trend is quite well predicted.

**Rotor pressure coefficients.** Rotor blade pressure coefficients from ASTEC compared to RANS results are visible on Fig. 13. They derive from the operation points circled on Fig. 12 for several reduced blade heights $h/H$.

<figure><img src="image_14.png" /><figcaption>FIGURE 13: CP comparison on rotor blade at 100Nn near design point for multiple reduced blade heights</figcaption></figure>

This time, there is no additional aerodynamic blockage coming from experimental values. Still, for $0.3 \le h/H \le 0.7$ a very good agreement is observed. Below $h/H \approx 0.3$ the influence of hub end-wall boundary layer is observed in the RANS results. Above $h/H \approx 0.7$ the influence of casing end-wall boundary layer and increasing compressible effects cause the observed discrepancies.
<figure><img src="image_15.png" /><figcaption>FIGURE 14: CME2 ASTEC mesh colored with static pressure at 100NN near design point, Approx. 6000 cells</figcaption></figure>

**Tip leakage model evaluation.** First a comparison is made between RANS and ASTEC at the operation points circled on Fig. 12 presenting same compressor massflow rate. Figure 15 shows the distribution of tip leakage massflow rate $dm$ of Eq. (50) normalized by the total leakage massflow rate against the chord-normalized axial position.

<figure><img src="image_16.png" /><figcaption>FIGURE 15: Tip leakage massflow distribution at 100NN near design point</figcaption></figure>

A very good agreement on the absolute value of total leakage massflow rate with an error less than 1% is observed and the overall trend along the blade tip is well captured by the model compared to the RANS simulation.

Then the evolution of tip leakage loss predicted by the model while throttling the compressor is shown on Fig. 16. The loss and the channel massflow rate $m_m$ are normalized by reference values extracted at the ASTEC circled operation point on Fig. 12. During the compressor throttle, the loss produced by the leakage flow shows an intense growth, up to close to tree time the loss produced at the reference point.

<figure><img src="image_17.png" /><figcaption>FIGURE 16: Tip leakage loss evolution with throttle at 100NN from the model</figcaption></figure>

## 7. CONCLUSION

The present paper benefits consist in the supply of 2.5D features to a CFD-based throughflow solver. The implementation of

the Hess & Smith method constitutes a base for further throughflow modeling with a reduction of empiricism. This method allows the estimation of the bi-dimensional potential flow field in blade-to-blade streamsurfaces. Pressure gradients on both suction and pressure sides of the blades are coupled with the calculation of profile boundary layers using an integral boundary layer method. Flow angles on streamsurfaces are deduced from the blade-to-blade averaged flow velocities and used to drive the throughflow blade force model. Wakes are taken into account with a decay model since boundary layer quantities are well-estimated at the trailing edge. The estimation of the tip leakage loss benefits from the blade tip pressure distribution which leads to an accurate prediction of the tip leakage mass flow rate. The loss induced by the hub secondary flows has not yet been implemented in the solver.

Bi-dimensional full turbulent results agree very well with data from empirical correlations of Lieblein, MISES and RANS simulations. Comparison with experiments (two cascades and a low speed single stage axial compressor) have shown encouraging results for pressure coefficient distributions on the blades far enough of the influence of hub secondary flows. Tip leakage loss is accurately predicted. The historical fast-response trait of throughflow solvers was also respected as a constraint in this work leading to a usual response time of 1 minute.

The main objective for the future of ASTEC is to pursue the modeling of secondary flows, especially the compressor hub boundary layers. This would bring radial distribution of losses and blockage, improved results on the axial compressor case and make the code ready for future realistic multistage applications. The current methodology developed involves subsonic potential flow theory, so a second objective is to extend this methodology to be suitable for higher Mach numbers. On the long term, ASTEC would be used in an analysis optimization process for design purposes.

## ACKNOWLEDGMENTS

The authors would like to thank Safran Aircraft Engines and the French Ministry of Higher Education and Research for the funding of this work. They are also thankful to Safran Tech for accepting to share ASTEC source code to conduct this research.

For the purpose of Open Access, a CC-BY public copyright licence has been applied by the authors to the present document and will be applied to all subsequent versions up to the Author Accepted Manuscript arising from this submission.

<figure><img src="image_18.png" /></figure>
REFERENCES

[1] Cumpsty, N. A. *Compressor Aerodynamics*. Longman Scientific & Technical; J. Wiley, Harlow, Essex, England: New York (1989).

[2] Horlock, J. H. and Denton, J. D. "A Review of Some Early Design Practice Using Computational Fluid Dynamics and a Current Perspective." *Journal of Turbomachinery* Vol. 127 No. 1 (2005): pp. 5-13. DOI 10.1115/1.1650379.

[3] Simon, J-F. "Contribution to Throughflow Modelling for Axial Flow Turbomachines." Ph.D. Thesis, University of Liege. 2007.

[4] Pacciani, R., Rubechini, F., Marconcini, M., Arnone, A., Cecchi, S. and Daccà, F. "A CFD-based Throughflow Method with an Explicit Body Force Model and an Adaptive Formulation for the S2 Streamsurface." *Proceedings of the Institution of Mechanical Engineers, Part A: Journal of Power and Energy* Vol. 230 No. 1 (2016): pp. 16-28. DOI 10.1177/0957650915607091.

[5] Ricci, M., Pacciani, R., Marconcini, M., Macelloni, P., Cecchi, S. and Bettini, C. "Computational Fluid Dynamics-Based Throughflow Analysis of Transonic Flows in Steam Turbines." *Journal of Turbomachinery* Vol. 141 No. 11 (2019): p. 111005. DOI 10.1115/1.4044871.

[6] Novak, R. A. "Streamline Curvature Computing Procedures for Fluid-Flow Problems." *Journal of Engineering for Power* Vol. 89 No. 4 (1967): pp. 478-490. DOI 10.1115/1.3616716.

[7] Baralon, S., Eriksson, L-E. and Häll, U. "Validation of a Throughflow Time-Marching Finite-Volume Solver for Transonic Compressors." *Turbomachinery*, Vol. 1: p. V001T01A015. 1998. American Society of Mechanical Engineers, Stockholm, Sweden. DOI 10.1115/98-GT-047.

[8] Banjac, M., Petrovic, M. V. and Wiedermann, A. "Secondary Flows, Endwall Effects, and Stall Detection in Axial Compressor Design." *Journal of Turbomachinery* (2015)DOI 10.1115/1.4028648.

[9] Manfredi, M. and Fontaneto, F. "Transonic Axial Compressors Loss Correlations: Part I — Analysis and Update of Loss Models." *Volume 2A: Turbomachinery*: p. V02AT32A029. 2020. American Society of Mechanical Engineers, Virtual, Online. DOI 10.1115/GT2020-14713.

[10] Hess, J.L. and Smith, A.M.O. "Calculation of Potential Flow about Arbitrary Bodies." *Progress in Aerospace Sciences* Vol. 8 (1967). DOI 10.1016/0376-0421(67)90003-6.

[11] Cambier, L., Heib, S. and Plot, S. "The Onera elsA CFD Software: Input from Research and Feedback from Industry." *Mechanics & Industry* Vol. 14 No. 3 (2013): pp. 159-174. DOI 10.1051/meca/2013056.

[12] Thollet, W., Dufour, G., Carbonneau, X. and Blanc, F. "Body-Force Modeling for Aerodynamic Analysis of Air Intake - Fan Interactions." *International Journal of Numerical Methods for Heat & Fluid Flow* Vol. 26 No. 7 (2016): pp. 2048-2065. DOI 10.1108/HFF-07-2015-0274.

[13] Adamczyk, J. J. "Aerodynamic Analysis of Multistage Turbomachinery Flows in Support of Aerodynamic Design." *Journal of Turbomachinery* Vol. 122 No. 2 (2000): pp. 189-217. DOI 10.1115/1.555439.

[14] Thomas, J. P. and Léonard, O. "Towards a High Order Throughflow: Part I—Investigating the Effectiveness of a Harmonic Reconstruction for 3D Flows." *Turbomachinery*, Vol. 7A, B, C: pp. 831-841. 2010. ASME, Glasgow, UK. DOI 10.1115/GT2010-22841.

[15] Denton, J. D. "Loss Mechanisms in Turbomachines." *ASME 1993 International Gas Turbine and Aeroengine Congress and Exposition*: p. V002T14A001. 1993. American Society of Mechanical Engineers, Cincinnati, Ohio, USA. DOI 10.1115/93-GT-435.

[16] Smith, L. H. "Casing Boundary Layers in Multistage Axial-Flow Compressors." Dzung, L. S. (ed.). *Flow Research on Blading*. 1970. Elsevier Pub. Co., Amsterdam.

[17] Budo, A., Terrapon, V. E., Arnst, M., Hillewaert, K., Mouriaux, S., Rodriguez, B. and Bartholet, J. "Application of a Viscous Throughflow Model to a Modern Axial Low-Pressure Compressor." *Turbomachinery* — *Axial Flow Fan and Compressor Aerodynamics*, Vol. 2A: p. V02AT31A038. 2021. American Society of Mechanical Engineers, Virtual, Online. DOI 10.1115/GT2021-59926.

[18] Lieblein, S. "Incidence and Deviation-Angle Correlations for Compressor Cascades." *Journal of Basic Engineering* Vol. 82 No. 3 (1960): pp. 575-584. DOI 10.1115/1.3662666.

[19] Katz, J. and Plotkin, A. *Low-Speed Aerodynamics*, 2nd ed. Cambridge University Press (2004). DOI 10.1017/CBO9780511810329.

[20] McFarland, E. R. "A Rapid Blade-to-Blade Solution for Use in Turbomachinery Design." *Journal of Engineering for Gas Turbines and Power* (1984)DOI 10.1115/1.3239574.

[21] Lewis, R. I. *Vortex Element Methods for Fluid Dynamic Analysis of Engineering Systems*, 1st ed. Cambridge University Press (1991). DOI 10.1017/CBO9780511529542.

[22] Head, M. R. "Entrainment in the Turbulent Boundary Layer." Technical Report ADA951554. Aeronautical Research Council. 1958.

[23] Ludweig, H; and Tillmann, W. "Investigations of The Wall-Shearing Stress in Turbulent Boundary Layers." Rapport Technique 1285. National Advisory Committee for Aeronautics. 1950.

[24] Lakshminarayana, B. and Davino, R. "Mean Velocity and Decay Characteristics of the Guidevane and Stator Blade Wake of an Axial Flow Compressor." *Journal of Engineering for Power* Vol. 102 No. 1 (1980): pp. 50-60. DOI 10.1115/1.3230231.

[25] Raj, R. and Lakshminarayana, B. "Characteristics of the Wake behind a Cascade of Airfoils." *Journal of Fluid Mechanics* Vol. 61 No. 4 (1973): pp. 707-730. DOI 10.1017/S002211207300090X.

[26] Moore, J. and Tilton, J. S. "Tip Leakage Flow in a Linear Turbine Cascade." *Journal of Turbomachinery* Vol. 110 No. 1 (1988): pp. 18-26. DOI 10.1115/1.3262162.

[27] Veldman, A. E. P. "Quasi-Simultaneous Viscous-Inviscid Interaction for Transonic Airfoil Flow (Invited)." *4th AIAA Theoretical Fluid Mechanics Meeting*. 2005. American Institute of Aeronautics and Astronautics, Toronto, Ontario, Canada. DOI 10.2514/6.2005-4801.

[28] Zambonini, Gherardo, Ottavy, Xavier and Kriegseis, Jochen. "Corner Separation Dynamics in a Linear Compressor Cascade." *Journal of Fluids Engineering* Vol. 139 No. 6 (2017): p. 061101. DOI 10.1115/1.4035876.

[29] Dawkins, I., Taylor, J., Ottavy, X. and Miller, R. "The Unsteady Topology of Corner Separations." *Journal of Turbomachinery* Vol. 144 No. 11 (2022): p. 111001. DOI 10.1115/1.4054676.

[30] Aungier, R. H. *Axial-Flow Compressors*. ASME Press (2003). DOI 10.1115/1.801926.

[31] Smith, B. "A near Wall Model for the k-1 Two Equation Turbulence Model." *Fluid Dynamics Conference*. 1994. American Institute of Aeronautics and Astronautics, Colorado Springs, CO, U.S.A. DOI 10.2514/6.1994-2386.

[32] Gourdain, N. "Prediction of the Unsteady Turbulent Flow in an Axial Compressor Stage. Part 1: Comparison of Unsteady RANS and LES with Experiments." *Computers & Fluids* Vol. 106 (2015): pp. 119-129. DOI 10.1016/j.compfluid.2014.09.052.