# Etude de l'Écoulement Tridimensionnel dans un Rotor d'une Turbine Semi-axiale

Omar Hireche, Miloud Abidat, Leila Merahi et Abbes Azzi

U.S.T.O. Mohammed Boudiaf, BP 1505, Oran El M'Naouer

(Accepté le 22. 05. 99)

**Résume** — Cet article fait le point d'une étude pour la simulation numérique d'un écoulement tridimensionnel dans le rotor d'une turbine semi-axiale par la résolution numérique des équations d'Euler (en fluide parfait) exprimées en coordonnées cylindriques. La méthode de calcul qui se caractérise par l'utilisation d'une technique de résolution par stationnarisation associée à un schéma explicite à viscosité corrigée et une approche par volumes finis, présente de bonnes propriétés de convergence et de précision.

**Abstract** — This paper presents a study carried out for the numerical simulation of tridimensional flow in a mixed flow turbine rotor by solving of the Euler equations (inviscid flows) written in a cylindrical system of coordinates. This numerical method is characterized by the use of a time marching technique with an explicit scheme and a finite volume approach. The method presents good properties of convergence and accuracy.

**MOTS CLES** : Turbine semi- axiale – Tridimensionnel - Non visqueux - Schéma explicite - Méthode numérique

## 1. INTRODUCTION

La résolution des équations de Navier-Stokes est le moyen le plus représentatif du calcul d'un écoulement réel, mais cette alternative est restée inaccessible pendant plusieurs années, et ce pour deux raisons essentielles. La première est liée au nombre de mémoire et aux temps de calcul, qui étaient incompatibles dans le cadre d'une utilisation fréquente avec les ordinateurs mis à la disposition des chercheurs. La seconde, qui est plus importante, résulte du constat que tous les schémas de turbulence proposés, et qui sont indispensables à la résolution des équations de Navier-Stokes en régime turbulent, ne rendaient pas compte avec une précision suffisante de la réalité physique, même pour des applications relativement simples.

Face aux problèmes soulevés, de nombreuses approches ont été entreprises; elles font appel à une simplification du problème en considérant l'écoulement non visqueux et obéissant aux équations différentielles d'un fluide parfait (équation d'Euler). Bien sûr cette hypothèse qui a ses limites d'application, reste valable dans la majorité des écoulements inter-aubes à grand nombre de Reynolds.

La grande majorité des méthodes présentées dans la littérature [1, 2] pour la résolution des équations d'Euler, sont des méthodes explicites dont on pourra distinguer deux grandes catégories d'approche. La première famille d'approche qui est connue sous le nom de Lax-Fiedrichs [2], a fait l'objet de plusieurs études ayant pour but l'amélioration de la précision et la diminution du temps de calcul. La seconde famille d'approche, qui se base sur le fameux schéma de MacCormack du type prédicteur - correcteur [2], présente un grand intérêt pour la résolution des problèmes instationnaires du fait que ce schéma est précis au second ordre par rapport au temps et l'espace. La méthode de résolution présentée dans cet article a été largement utilisée dans le calcul des écoulements dans les turbomachines [3,4]. Cette méthode
a été récemment utilisée pour la résolution d'un écoulement bidimensionnel à travers une grille d'aubes plane pour les cas non visqueux [5] et visqueux [6]. Le choix du schéma à viscosité corrigée qui est une extension du schéma de Lax Friedrichs a été dicté par le fait qu'on ne s'intéresse qu'à la solution stationnaire.

## 2. EQUATIONS DE BASE

L'écoulement adiabatique et non visqueux d'un fluide compressible dans un rotor d'une turbomachine est décrit par les équations d'Euler exprimées dans un repère cylindrique (r,θ,z) et tournant autour de l'axe de la machine avec une vitesse angulaire ω. La vitesse relative $\vec{w}$ et la vitesse absolue $\vec{v}$ d'une particule fluide éloignée de l'axe de rotation d'une distance r, sont reliées entre elles par la relation suivante:

$$ \vec{V} = \vec{W} + \omega r \vec{i}_0 \quad (1) $$

Dans un système de coordonnées cylindriques, les équations d'Euler s'écritent sous la forme compacte conservative suivante:

$$ \frac{\partial \vec{\sigma}}{\partial t} + \frac{1}{r} \frac{\partial}{\partial r} \left[ \vec{r} \vec{f}(\vec{\sigma}) \right] + \frac{1}{r} \frac{\partial}{\partial \theta} \left[ \vec{g}(\vec{\sigma}) \right] + \frac{\partial}{\partial z} \left[ \vec{h}(\vec{\sigma}) \right] + \vec{b}(\vec{\sigma}) = \vec{0} \quad (2) $$

où :

$$ \vec{\sigma} = \begin{bmatrix} \rho \\ \rho w_r \\ \rho w_\theta \\ \rho w_z \\ \rho e \end{bmatrix}, $$

$$ \vec{f}(\vec{\sigma}) = \begin{bmatrix} \rho w_r \\ \rho w_r^2 + p \\ \rho w_r w_\theta \\ \rho w_r w_z \\ \rho w_r H \end{bmatrix}; $$

$$ \vec{g}(\vec{\sigma}) = \begin{bmatrix} \rho w_\theta \\ \rho w_\theta w_r \\ \rho w_\theta^2 + p \\ \rho w_\theta w_z \\ \rho w_\theta H \end{bmatrix}; $$

$$ \vec{h}(\vec{\sigma}) = \begin{bmatrix} \rho w_r \\ \rho w_r w_z \\ \rho w_z w_\theta \\ \rho w_z^2 + p \\ \rho w_z H \end{bmatrix}; $$

$$ \vec{b}(\vec{\sigma}) = \begin{bmatrix} 0 \\ -\rho w_\theta^2 / r - p / r - \rho r \omega^2 - 2\rho\omega w_\theta \\ \rho w_r w_\theta / r + 2\rho\omega w_r \\ 0 \\ 0 \end{bmatrix} $$

$$ \text{avec : } \begin{cases} H = e + \frac{p}{\rho} \\ p = (\gamma - 1) \left( \rho e - \frac{\rho W^2}{2} + \rho \omega^2 \frac{r^2}{2} \right) \\ e = c_V T + \frac{W^2}{2} - \frac{\omega^2 r^2}{2} \end{cases} $$
### 3. METHODE NUMERIQUE

Le domaine physique de l'écoulement est constitué par le canal compris entre deux aubes successives. Pour pouvoir définir correctement les conditions à l'entrée et à la sortie du domaine de l'écoulement, celui-ci est prolongé en amont et en aval pour éviter les perturbations à l'entrée et à la sortie du canal inter-aubes. Le domaine physique de l'écoulement est subdivisé en éléments de volumes dans lesquels les équations d'Euler sont intégrées. Chaque élément de volume est limité par les surfaces suivantes:

* **Surface aube à aube:** C'est une surface qui est limitée par l'intrados de l'aube et l'extrados de l'aube voisine. Elle est générée par la rotation d'une pseudo-ligne de courant et est étendue jusqu'aux frontières d'entrée et de sortie du domaine.

* **Surface méridienne:** C'est la projection de l'aube sur le plan méridien (r, z) comme le montre la Figure 1.

<figure><img src="image_1.png" /><figcaption>Fig. 1 : Discretisation de la surface méridienne</figcaption></figure>

* **Surface orthogonale:** C'est une surface qui est située sur un plan perpendiculaire au plan méridien et qui est générée par la rotation d'une ligne quasi-orthogonale autour de l'axe de rotation.

Le volume de contrôle choisi pour la discrétisation des équations (2) est un élément bitrapézoïdal, constitué de deux éléments de volume consécutifs. Les variables de conservation sont stockées dans des noeuds situés aux centres des faces se trouvant sur des surfaces orthogonales (Figure 2).

<figure><img src="image_2.png" /><figcaption>Fig. 2 : Volume de contrôle</figcaption></figure>
Le système d'équations (2) est résolu par voie instationnaire associée à une discrétisation en volumes finis. Le schéma numérique utilisé est le schéma à viscosité corrigée qui est une extension du schéma de Lax-Fiedrichs [2] et qui a été introduit par Mc Donald et Couston [3,7]. Le schéma à viscosité corrigée se caractérise par une discrétisation des termes dépendants des coordonnées spatiales à l'aide d'une différence centrée et des termes dépendants du temps par une différence avant. Pour un noeud intérieur du domaine de l'écoulement, ce schéma s'écrit comme suit :

$$
\begin{align}
\sigma_{i,j,k}^{n+1} = \sigma_{i,j,k}^{n} & - \frac{\Delta t}{V} \left[ \iiint_V (\mathrm{div} \, \vec{F}) dV \right]^n \nonumber \\
& + \frac{1}{6} \left( \sigma_{i+1,j,k}^{n} + \sigma_{i-1,j,k}^{n} + \sigma_{i,j+1,k}^{n} + \sigma_{i,j-1,k}^{n} + \sigma_{i,j,k+1}^{n} + \sigma_{i,j,k-1}^{n} - 6\sigma_{i+1,j,k}^{n} \right) \nonumber \\
& - \frac{\alpha}{6} \left( \sigma_{i+1,j,k}^{0} + \sigma_{i-1,j,k}^{0} + \sigma_{i,j+1,k}^{0} + \sigma_{i,j-1,k}^{0} + \sigma_{i,j,k+1}^{0} + \sigma_{i,j,k-1}^{0} - 6\sigma_{i,j,k}^{0} \right) \tag{3} \\
\text{avec } \vec{F} = f(\vec{\sigma})I_r + g(\sigma)I_{\theta} + h(\sigma)I_z & \nonumber
\end{align}
$$

L'intégrale de volume est remplacée par une intégrale de surface de la fonction F grâce au théorème de Gauss. Cette dernière est égale à la somme des intégrales de surface à travers les six faces du volume de contrôle.

$$
\iiint_V \left( \operatorname{div} \vec{F} \right) dV = \iint_S \vec{F} \cdot d\vec{S} = \sum_{m=1}^{6} \left[ f_m \, dS_r + g_m \, dS_{\theta} + h_m \, dS_z \right] \quad (4)
$$

où :

$$
\vec{d} S = dS_r \quad \vec{I}_r + dS_\theta \quad \vec{I}_\theta + dS_z \quad \vec{I}_Z \text{ est l'élément de surface,}
$$

$f_m, g_m$, et $h_m$ sont les valeurs moyennes des termes convectifs $f, g$ et $h$ sur la face $m$
($m=1$ à 6) du volume de contrôle dans lequel le système d'équations (2) est intégré.

Les termes surmontés d'un zéro sont mis à jour toutes les NV itérations (NV=15....25), tandis
que le coefficient de viscosité artificielle α est calculé de la manière suivante [3,4]:

$$
\alpha = \beta \left( 1 - \frac{\rho_{i+1,j,k}^{n} + \rho_{i-1,j,k}^{n} + \rho_{i,j+1,k}^{n} + \rho_{i,j-1,k}^{n} + \rho_{i,j,k+1}^{n} + \rho_{i,j,k-1}^{n} - 6\rho_{i,j,k}^{n}}{6\rho_{01}} \right) \quad (5)
$$

où $\beta$ est un coefficient numérique tel que $0<\beta<1$. Le pas temporel $\Delta t$ est déterminé localement en chaque noeud à partir de la condition de stabilité de C.F.L. dont la formulation mathématique est comme suit:

$$
\Delta t \le \frac{\Delta l}{|W + \sqrt{\gamma RT}|} \qquad (6)
$$

où Δl est la distance minimale entre le noeud considéré et les six faces du volume de contrôle.
Les différents étapes de calcul sont résumés dans l'organigramme de la figure3.
<figure><img src="image_3.png" /><figcaption>Fig. 3 : Organigramme de Calcul</figcaption></figure>

## 3.1. CONDITIONS AUX LIMITES

Des conditions aux limites sont appliquées aux noeuds adjacents à une frontière solide, à une frontière périodique, à l'entrée ou à la sortie du domaine de l'écoulement. Au voisinage d'une paroi solide, la condition d'imperméabilité est appliquée tandis qu'au voisinage d'une frontière périodique, une condition de périodicité
de l'écoulement est appliquée. En ce qui concerne les noeuds adjacents à l'entrée ou à la sortie du domaine et étant donné que le système d'équations aux dérivées partielles (eq. 1) est du type hyperbolique par rapport au temps, le nombre de conditions à imposer est déterminé à l'aide de la méthode des caractéristiques. Ces conditions sont résumées comme suit :

<table><caption>**Tableau 1:** Nombre et type de conditions aux limites aux sections d'entrée et de sortie</caption><thead><tr><td>Type d'écoulement</td><td>Entrée</td><td>Sortie</td></tr></thead><tbody><tr><td>Subsonique (M &lt; 1)</td><td>4 conditions<br/>1 équation à résoudre</td><td>1 condition<br/>4 équations à résoudre</td></tr><tr><td>Supersonique (M &gt; 1)</td><td>5 conditions<br/>0 équation à résoudre</td><td>Aucune condition<br/>5 équations à résoudre</td></tr></tbody></table>

Un schéma excentré est appliqué aux volumes dont une face est située sur une des frontières solides du domaine de l'écoulement tandis que les noeuds adjacents à une frontière périodique sont traités de la même manière que les noeuds intérieurs.

Afin de commencer le calcul itératif pour la résolution du système d'équations (2), il est nécessaire d'adjoindre à celles-ci une solution initiale arbitraire du fait que le choix de cette dernière n'a aucune influence sur la solution finale étant donné qu'on ne s'intéresse qu'à la solution stationnaire de l'écoulement.

## 4. RESULTATS ET CONCLUSION

La méthode présentée dans cet article a été appliquée au calcul de l'écoulement tridimensionnel et non visqueux dans un rotor d'une turbine semi-axiale utilisée dans la suralimentation des moteurs à combustion interne et dont les caractéristiques géométriques sont données par Abidat et al [8].

<figure><img src="image_4.png" /><figcaption>Fig. 4: Evolution de l'erreur</figcaption></figure>
Les résultats présentés ont été obtenus pour les conditions de fonctionnement suivantes :

* Pression d'arrêt à l'entrée de la turbine : $P_0^* = 290$ KPa

* Température d'arrêt à l'entrée de la turbine : $T_0^* = 345$ K

* Pression statique à la sortie de la turbine : $P_s = 92$ KPa

* Vitesse de rotation : $N = 59740$ tr / mn

Les résultats sont présentés en termes d'évolution de l'erreur sur le débit et sur la vitesse relative (Figure 4), de distribution de la pression statique sur la surface de l'aube (Figures 5) et de champs des vitesses relatives sur différentes surfaces du canal inter - aubes (Figures 6, 7 et 8).

Les distributions de la pression statique sur la surface de l'aube présentées en figure 5, montrent que la charge de l'aube est plus élevée du coté du carter que du coté du moyeu. On remarque cependant une charge négative de l'aube à cause d'une très large négative incidence de l'écoulement près du carter au niveau du bord d'attaque. Les champs de vitesses représentés en figures 6 à 8 montrent une uniformité de l'écoulement principal (surface méridienne).

<figure><img src="image_5.png" /><figcaption>a - Prés du moyeu.</figcaption></figure>

<figure><img src="image_6.png" /><figcaption>a - Prés de l'entrée du rotor</figcaption></figure>

<figure><img src="image_7.png" /><figcaption>b - Au Rayon moyen.</figcaption></figure>

<figure><img src="image_8.png" /><figcaption>b - Milieu du canal</figcaption></figure>
<figure><img src="image_9.png" /><figcaption>Fig.5: Distribution de la pression statique
sur la surface de l'aube</figcaption></figure>

<figure><img src="image_10.png" /><figcaption>Fig.6: Champs de la vitesse relative
sur une surface orthogonale</figcaption></figure>

Les distributions de la pression statique sur la surface de l'aube présentées en figures 5, montrent que la charge de l'aube est plus élevée du coté du carter que du coté du moyeu. On remarque cependant une charge négative de l'aube à cause d'une très large négative incidence de l'écoulement près du carter au niveau du bord d'attaque. Les champs de vitesses représentés en figures 6 à 8 montrent une uniformité de l'écoulement principal (surface méridienne). La projection de la vitesse sur des surfaces orthogonales montre clairement l'effet des forces de Corriollis qui se traduisent par l'apparition d'un écoulement secondaire sous forme d'une zone tourbillonnaire dans la partie amont du canal inter-aubes.

La méthode de calcul présentée dans cet article peut donc être utilisée pour l'analyse des écoulements dans les turbomachines et peut constituer un outil appréciable lors de la conception de ce type de machines.

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

<figure><img src="image_12.png" /></figure>
<figure><img src="image_13.png" /><figcaption>Fig.7: Champs de la vitesse relative sur une surface méridienne</figcaption></figure>

<caption>Fig.8: Champs de la vitesse relative sur une surface aube à aube</caption>

## NOMENCLATURE

<table><tr><td>c<sub>v</sub></td><td>Chaleur spécifique à volume constant</td><td>j/kg K</td><td>V</td><td>Volume vitesse absolue</td><td>m<sup>3</sup> m/s</td></tr><tr><td>e</td><td>Energie totale</td><td>j/kg</td><td>W</td><td>Vitesse relative</td><td>m/s</td></tr><tr><td>f</td><td>Terme de convection suivant r</td><td>-</td><td>w</td><td>Composante de la vitesse relative</td><td>m/s</td></tr><tr><td>g</td><td>Terme de convection suivant θ</td><td>-</td><td>ρ</td><td>Masse volumique</td><td>kg/m<sup>3</sup></td></tr><tr><td>h</td><td>Terme de convection suivant z</td><td>-</td><td>γ</td><td>Rapport des chaleurs spécifiques</td><td>-</td></tr><tr><td>p</td><td>Pression statique</td><td>Pa</td><td>ω</td><td>Vitesse de rotation du rotor</td><td>tr/mn</td></tr><tr><td>r,θ,z</td><td>Coordonnées cylindriques</td><td>m, rad, m</td><td>α</td><td>Coefficient de viscosité numérique</td><td>-</td></tr><tr><td>t</td><td>Temps.</td><td>s</td><td>β</td><td>Coefficient</td><td>-</td></tr><tr><td>S:</td><td>Surface</td><td>m<sup>2</sup></td><td colspan="3"><strong>Indices</strong></td></tr></table>
<table><tr><td>Δt</td><td>Pas temporelle</td><td>s</td><td>i, j, k les indices de coordonnée</td><td>-</td></tr><tr><td>I</td><td>Vecteur unitaire</td><td>-</td><td></td><td></td></tr><tr><td>T</td><td>Température</td><td>K</td><td>n indice de temps</td><td>-</td></tr><tr><td>R</td><td>Constante des gaz</td><td>j/kg K</td><td></td><td></td></tr></table>

REFERENCES

[1] P. W. Mc Donald “The Computation of Transonic Flow Through Two Dimensional Gaz Turbine Cascades” ASME, Paper 71-GT-89, 1971.

[2] C. Hirsch “Numerical Computation of Internal and External Flows”, vol. 2, John Wiley & Sons, 1991.

[3] W. Van Hove. “Calculation of Three Dimensional Inviscid Rotational Flow in Axial Turbine Blade Cascade » ”ASME, Paper 83-GT-119, 1984.

[4] T. Arts. “Calculation of the Three Dimensional Steady, Inviscid Flow in a Transonic Axial Turbine Stage” I.Engrg Gaz Turbines and Power, vol 107, N°2 April 1985.

[5] M. Abidat, L. Merahi et M.K. Hamidou “Application d'une Méthode des Volumes Finis au Calcul des Ecoulements dans une Grille d'Aubes de Turbomachines” CIMASI'96, Casablanca du 14 au 16 Novembre 1996.

[6] M. Abidat, O. Hireche, L. Merahi et A. Azzi “Calcul de l'Ecoulement Visqueux dans une Grille d'Aubes Plane” 4ème Congrès de Mécanique, Mohammadia du 13 au 16 Avril 1999.

[7] M. Couston, P. W. Mc Donald and J. J. Smolderen “The Damping Surface Technique for Time Dependent Solutions to Fluid Dynamic Problems” VKI - TN 19, 1975.

[8] M. Abidat, H. Chen, N. C. Baines and M. R. Firth “Design of a Highly Loaded Mixed Flow Turbine.” Proc IMechE, Journal of Power and Energy, 1992.