# Étude de schémas à haute précision pour la simulation d'écoulements transsoniques instationnaires ou visqueux.
## Application aux turbomachines

Stéphane Aubert

► To cite this version:

Stéphane Aubert. Étude de schémas à haute précision pour la simulation d'écoulements transsoniques instationnaires ou visqueux. Application aux turbomachines. Mécanique des fluides [physics.class-ph]. École Centrale de Lyon, 1993. Français. NNT: 1993ECDL0033. tel-04024181

HAL Id: tel-04024181

https://hal.science/tel-04024181

Submitted on 10 Mar 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 Attribution 4.0 International License
THESE

présentée devant

L'ECOLE CENTRALE DE LYON

pour obtenir

le titre de DOCTEUR

Spécialité MECANIQUE

par M. AUBERT Stéphane

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

ETUDE DE SCHEMAS A HAUTE PRECISION POUR LA
SIMULATION D'ECOULEMENTS TRANSSONIQUES
INSTATIONNAIRES OU VISQUEUX. APPLICATION AUX
TURBOMACHINES

Soutenue le 8 octobre 1993 devant la Commission d'Examen

Composition du jury

Professeur M.H. ATASSI (rapporteur)
Professeur A. BOLCS
Docteur W.N. DAWES
Docteur P. FERRAND
Professeur F. LEBOEUF
Professeur J. PIQUET (rapporteur)

Laboratoire de Mécanique des Fluides et d'Acoustique URA 263
Ecole Centrale de Lyon
Un homme qui ne commet jamais d'erreurs a depuis
longtemps cessé de faire du neuf. Un homme qui commet
toujours des erreurs est un homme condamné, un homme
aux ambitions démesurées. Mais celui qui saupoudre judi-
cieusement d'erreurs son succès apprend vite.

Reeho'na le o'Tghalie
dans *Les Mathématiques du Savoir*

A mes maîtres,
A mes asuthi
Remerciements

Cette étude a été réalisée au sein du Laboratoire de Mécaniques des Fluides et d'Acoustique, URA 263, de l'Ecole Centrale de Lyon, et de la société Metraflu. Elle s'inscrit dans le cadre des travaux de recherches théoriques du groupe Turbomachines. Elle a été cofinancée par le CNRS et la société Metraflu par une Bourse de Docteur Ingénieur. Aussi, je tiens à remercier Messieurs D.Jeandel et J. Bataille, du LMFA, et Messieurs G. Bois et E. Benarous, de Metraflu, pour leur accueil et leur soutien à ce travail.

Je suis reconnaissant envers Monsieur P. Ferrand, Chargé de Recherches au CNRS/LMFA, qui a assuré l'encadrement de cette thèse. Je tiens à le remercier pour sa grande disponibilité, son ouverture d'esprit, les intenses discussions scientifiques que nous avons partagées, ainsi que pour sa profonde humanité et sa gentillesse.

Ma gratidute va à Monsieur F. Leboeuf, Professeur à l'ECL, pour m'avoir incité à entreprendre ces travaux et pour m'avoir acceilli au sein du groupe Turbo-machines.

Mes plus vifs remerciements vont à Messieurs H.M. Atassi, Professeur à l'Université de Notre Dame, Indiana, Etats-Unis, et J. Piquet, Professeur à l'Ecole Centrale de Nantes, pour avoir accepté la charge de rapporteur de cette thèse.

La présence dans le jury de Messieurs W.N. Dawes, enseignant à l'Université de Cambridge, Grande-Bretagne, et A. Bölcs, Professeur à l'Ecole Polytechnique Fédérale de Lausanne, Suisse, m'honorent. Je les remercie d'autre part pour les contacts que nous avons eu dans le cadre de ce travail, et de ses extensions.

J'adresse enfin mes remerciements à toutes les personnes, condisciples et permanents, du LMFA et de Metraflu qui m'ont apporté leur soutien lors du déroulement de ce travail. Que soient assurés particulièrement ici de ma profonde gratitude Madame F. Maupas qui est une seconde mère pour nous, Monsieur G. Perrin
qui persiste sainement à programmer en C dans un monde inféodé au Fortran pour
ses programmes de visualisation et de maillage, et enfin Monsieur L. Hallo qui
restera, je le souhaite, un compagnon de route propre à stimuler mes travaux, et à
supporter mes extravagances et ma mauvaise humeur en période de stress.
<table><tbody><tr><td rowspan="5">MATHEMATIQUES<br/>INFORMATIQUE SYSTEMES</td><td>DAVID Bertrand</td><td>PROFESSEUR</td></tr><tr><td>MARION Martine</td><td>---</td></tr><tr><td>MAITRE Jean-François</td><td>---</td></tr><tr><td>THOMAS Gérard</td><td>---</td></tr><tr><td>MUSY François<br/>ROZINOER Jean</td><td>MAITRE DE CONFERENCES<br/>---</td></tr><tr><td rowspan="14">MECANIQUE DES FLUIDES</td><td>MATHIEU Jean</td><td>PROFESSEUR EMERITE<br/>PROFESSEUR</td></tr><tr><td>ALCARAZ Emilio</td><td>---</td></tr><tr><td>JEANDEL Denis</td><td>---</td></tr><tr><td>LEBOEUF Francis</td><td>---</td></tr><tr><td>SCOTT Jean</td><td>---</td></tr><tr><td>BATAILLE Jean</td><td>PROFESSEUR LYON I</td></tr><tr><td>GAY Bernard</td><td>---</td></tr><tr><td>GENCE Jean-Noël</td><td>---</td></tr><tr><td>LANCE Michel</td><td>---</td></tr><tr><td>BEN HADID Hamda</td><td>MAITRE DE CONFERENCES LYON I</td></tr><tr><td>BUFFAT Marc</td><td>---</td></tr><tr><td>HAMADICHE Mahmoud</td><td>---</td></tr><tr><td>MOREL Robert</td><td>PROFESSEUR INSA</td></tr><tr><td>CAMBON Claude</td><td>DIRECTEUR DE RECHERCHE CNRS</td></tr><tr><td>BERTOGLIO Jean-Pierre</td><td>CHARGE DE RECHERCHE CNRS</td></tr><tr><td>ESCUDIE Dany</td><td>---</td></tr><tr><td>FERRAND Pascal</td><td>---</td></tr><tr><td>HENRY Daniel</td><td>---</td></tr><tr><td rowspan="4">MECANIQUE DES SOLIDES</td><td>CAMBOU Bernard</td><td>PROFESSEUR</td></tr><tr><td>JEZEQUEL Louis</td><td>---</td></tr><tr><td>SIDOROFF François</td><td>---</td></tr><tr><td>SURRY Claude</td><td>PROFESSEUR ENISE</td></tr><tr><td rowspan="10">PHYSICOCHIMIE DES<br/>INTERFACES</td><td>CLECHET Paul</td><td>PROFESSEUR</td></tr><tr><td>JOSEPH Jacques</td><td>---</td></tr><tr><td>MARTELET Claude</td><td>---</td></tr><tr><td>MARTIN Jean-René</td><td>---</td></tr><tr><td>PORTE Louis</td><td>---</td></tr><tr><td>ROBACH Yves</td><td>MAITRE DE CONFERENCES</td></tr><tr><td>STREMSDOERFER Guy</td><td>---</td></tr><tr><td>JAFFREZIC Nicole</td><td>DIRECTEUR DE RECHERCHE CNRS</td></tr><tr><td>SOUTEYRAND Elyane</td><td>CHARGE DE RECHERCHE CNRS</td></tr></tbody></table>
PHOTACATALYSE

PICHAT Pierre
HERRMANN Jean-Marie

DIRECTEUR DE RECHERCHE CNRS

...

TECHNOLOGIE DES SURFACES

GEORGES Jean-Marie

PROFESSEUR

SABOT Jean

...

MARTIN Jean-Michel

...

MATHIA Thomas

DIRECTEUR DE RECHERCHE CNRS

KAPSA Philippe

...

LOUBET Jean-Luc

CHARGE DE RECHERCHE CNRS

LE BOSSE Jean-Claude

MAITRE DE CONFERENCES INSA

LOPEZ Jacques

MAITRE DE CONFERENCES UCB

ROUSSEAU Jacques

PROFESSEUR ENISE
# Résumé

L'étude de schémas numériques à haute précision appliqués à la simulation des phénomènes physiques régissant un écoulement transsonique instationnaire, voir turbulent, est l'objet principal de ce travail.

La méthode développée est bâtie sur les équations de Navier-Stokes moyennées, écrites sous forme conservative. Elles sont intégrées en espace dans le cadre d'une formulation volume fini en maillage curviligne structuré, mobile et déformable. L'intégration en temps est explicite. Les flux convectifs sont évalués par une approche MUSCL, couplée à un schéma mixte, résultant de l'hybridation en fonction du nombre de Mach d'un schéma centré et d'une technique de séparation de flux proposée par Van Leer. Les termes diffusifs sont quant à eux traiter par un schéma centré, la turbulence étant prise en compte à travers l'hypothèse de Boussinesq et le modèle algébrique de Michel, du type longueur de mélange. Enfin, les conditions aux limites géométriques sont gérées via l'extension du domaine de calcul par la création de noeuds supplémentaires, alors qu'un traitement assez sophistiqué des conditions physiques est mis en place à partir des équations de compatibilité.

Quelques configurations d'écoulements caractéristiques ont par ailleurs été simulées et analysées. Une description de grande qualité de discontinuités instationnaires et de leurs interactions a été obtenue. D'autre part, ont été simulés des écoulements multi-dimensionnels confinés stationnaires, recelant ou non des chocs, pour un fluide parfait ou en régime turbulent, et ce pour un grand nombre de nombres de Mach. Il en ressort une bonne capture des discontinuités, un traitement satisfaisant des conditions aux limites, une gestion correcte des effets visqueux, en particulier des interactions choc/couche limite, ce qui établit la qualité du schéma mixte. Une meilleure compréhension du comportement des méthodes numériques, entre autres de la technique de séparation des flux utilisée, en est un autre apport.

**Mots clefs:** Ecoulement compressible, Instationnaire, Visqueux, Turbomachines, Volumes finis, Séparation de flux, Hybridation, Relations de compatibilité
Abstract

The main purpose of this work is the study of highly accurate numerical schemes applied to simulate transonic unsteady, or turbulent, flows.

The method developed is based on the averaged conservative Navier-Stokes equations. These are integrated in space by means of a finite volume formulation over a moving curvilinear structured mesh. The time integration is explicit. Convective fluxes are evaluated by a MUSCL (Monotone Upstream-centred Scheme for Conservation Laws) approach, coupled with a mixed scheme. This is based on the hybridization according to Mach number of a centred scheme and of the flux vector splitting technique of Van Leer. Diffusive terms are computed using a centred scheme, and turbulence effects are taken into account by the Boussinesq hypothesis and by a classical mixing length model. Geometric boundary conditions are treated by extending the computational domain by two cells beyond the physical border. Physical boundary conditions are implemented by means of compatibility relations.

The method developed has been validated by simulating various steady and unsteady flows. Propagating shockwaves and contact discontinuities together with the interactions between them were accurately modelled. Confined multi-dimensional flows with and without discontinuities were also computed at various Mach numbers. These computations were performed with both inviscid and turbulent flow assumptions. Discontinuities were well captured, and the treatment of the boundary conditions proved satisfactory. Viscous effects were correctly predicted, including the complex shockwave/boundary layer interactions. One result of this work has been a better understanding of the numerical techniques adopted, especially the flux vector splitting scheme.

Key words: Compressible flow, Unsteady, Viscous, Turbomachinery, Finite volume, Flux splitting, Hybridization, Compatibility relations.
# Table des matières

<table><tr><td>1 Introduction</td><td>1</td></tr><tr><td>2 Modélisation d'un écoulement compressible turbulent</td><td>7</td></tr><tr><td>2.1 Introduction</td><td>7</td></tr><tr><td>2.2 Equations instantanées</td><td>7</td></tr><tr><td>2.3 Equations moyennées</td><td>9</td></tr><tr><td>2.4 Fermeture et modèles de turbulence</td><td>10</td></tr><tr><td>2.5 Système d'équations différentielles à résoudre</td><td>10</td></tr><tr><td>2.6 Conclusion</td><td>11</td></tr><tr><td>3 Résolution numérique des équations d'évolution</td><td>13</td></tr><tr><td>3.1 Introduction</td><td>13</td></tr><tr><td>3.2 Formulation volumes finis</td><td>14</td></tr><tr><td>3.2.1 Intégration en espace</td><td>14</td></tr><tr><td>3.2.2 Expression des flux</td><td>18</td></tr><tr><td>3.2.3 Equations discrétisées</td><td>21</td></tr><tr><td>3.3 Discrétisation temporelle</td><td>22</td></tr><tr><td>3.3.1 Remarques sur les schémas implicites en instationnaire</td><td>23</td></tr><tr><td>3.3.2 Schéma d'Euler explicite</td><td>23</td></tr><tr><td>3.4 Discrétisation spatiale des flux convectifs</td><td>24</td></tr><tr><td>3.4.1 Schémas décentrés amont</td><td>25</td></tr><tr><td>3.4.2 Discussion autour des limiteurs</td><td>32</td></tr><tr><td>3.4.3 Problèmes liés à la séparation de flux</td><td>37</td></tr><tr><td>3.4.4 Schémas centrés</td><td>53</td></tr><tr><td>3.4.5 Schémas mixtes</td><td>55</td></tr><tr><td>3.5 Discrétisation spatiale des flux diffusifs</td><td>56</td></tr><tr><td>3.6 Modélisation de la turbulence</td><td>57</td></tr><tr><td>3.6.1 Modèle de Prandtl</td><td>58</td></tr><tr><td>3.6.2 Modèle de Michel</td><td>59</td></tr><tr><td>3.6.3 Comparaisons et critiques</td><td>61</td></tr><tr><td>3.7 Résolution numérique</td><td>64</td></tr><tr><td>3.7.1 Progression en temps</td><td>64</td></tr><tr><td>3.7.2 Evaluation de la convergence</td><td>65</td></tr></table>
<table>
  <tbody>
    <tr>
      <td>3.8 Conclusion</td>
      <td>67</td>
    </tr>
    <tr>
      <td>4 Traitement des conditions aux limites</td>
      <td>69</td>
    </tr>
    <tr>
      <td>4.1 Introduction</td>
      <td>69</td>
    </tr>
    <tr>
      <td>4.2 Conditions géométriques</td>
      <td>71</td>
    </tr>
    <tr>
      <td>4.2.1 Extension du domaine de calcul</td>
      <td>72</td>
    </tr>
    <tr>
      <td>4.2.2 Différentes conditions géométriques</td>
      <td>73</td>
    </tr>
    <tr>
      <td>4.2.3 Gestion de maillages en O ou en C</td>
      <td>82</td>
    </tr>
    <tr>
      <td>4.2.4 Gestion d'un multi-domaines</td>
      <td>86</td>
    </tr>
    <tr>
      <td>4.3 Conditions physiques fortes</td>
      <td>87</td>
    </tr>
    <tr>
      <td>4.3.1 Equations d'évolution en variables caractéristiques</td>
      <td>88</td>
    </tr>
    <tr>
      <td>4.3.2 Imposition forte d'une condition aux limites</td>
      <td>93</td>
    </tr>
    <tr>
      <td>4.3.3 Différentes conditions physiques fortes</td>
      <td>96</td>
    </tr>
    <tr>
      <td>4.3.4 Avantages et inconvénients</td>
      <td>117</td>
    </tr>
    <tr>
      <td>4.4 Conditions physiques faibles</td>
      <td>118</td>
    </tr>
    <tr>
      <td>4.4.1 Formulation caractéristique multi-directionnelle</td>
      <td>118</td>
    </tr>
    <tr>
      <td>4.4.2 Imposition faible d'une condition aux limites</td>
      <td>119</td>
    </tr>
    <tr>
      <td>4.4.3 Différentes conditions physiques faibles</td>
      <td>123</td>
    </tr>
    <tr>
      <td>4.4.4 Conditions de non réflexion</td>
      <td>130</td>
    </tr>
    <tr>
      <td>4.4.5 Avantages et inconvénients</td>
      <td>132</td>
    </tr>
    <tr>
      <td>4.5 Résolution numérique</td>
      <td>133</td>
    </tr>
    <tr>
      <td>4.5.1 Gestion des noeuds supplémentaires</td>
      <td>133</td>
    </tr>
    <tr>
      <td>4.5.2 Gestion des noeuds frontières</td>
      <td>134</td>
    </tr>
    <tr>
      <td>4.6 Conclusion</td>
      <td>135</td>
    </tr>
    <tr>
      <td>5 Mouvements du maillage</td>
      <td>137</td>
    </tr>
    <tr>
      <td>5.1 Introduction</td>
      <td>137</td>
    </tr>
    <tr>
      <td>5.2 Adaptation statique</td>
      <td>138</td>
    </tr>
    <tr>
      <td>5.3 Adaptation dynamique</td>
      <td>145</td>
    </tr>
    <tr>
      <td>5.4 Conclusion</td>
      <td>146</td>
    </tr>
    <tr>
      <td>6 Simulation et analyse d'écoulements compressibles</td>
      <td>147</td>
    </tr>
    <tr>
      <td>6.1 Introduction</td>
      <td>147</td>
    </tr>
    <tr>
      <td>6.2 Propagation d'ondes mono-dimensionnelles</td>
      <td>148</td>
    </tr>
    <tr>
      <td>6.3 Interactions d'ondes mono-dimensionnelles</td>
      <td>154</td>
    </tr>
    <tr>
      <td>6.4 Interactions d'ondes bi-dimensionnelles</td>
      <td>163</td>
    </tr>
    <tr>
      <td>6.5 Tuyère supersonique</td>
      <td>171</td>
    </tr>
    <tr>
      <td>6.6 Tuyère transsonique</td>
      <td>181</td>
    </tr>
    <tr>
      <td>6.6.1 Ecoulements bi-dimensionnels non visqueux</td>
      <td>183</td>
    </tr>
    <tr>
      <td>6.6.2 Ecoulements bi-dimensionnels visqueux</td>
      <td>189</td>
    </tr>
    <tr>
      <td>6.6.3 Ecoulement bi-dimensionnel instationnaire</td>
      <td>202</td>
    </tr>
    <tr>
      <td>6.7 Grille d'aubes supersonique</td>
      <td>209</td>
    </tr>
    <tr>
      <td>6.8 Coude subsonique</td>
      <td>229</td>
    </tr>
  </tbody>
</table>
<table>
  <tbody>
    <tr>
      <td>6.9 Conclusion</td>
      <td>235</td>
    </tr>
    <tr>
      <td>7 Conclusions et perspectives</td>
      <td>241</td>
    </tr>
    <tr>
      <td>A Moyenne des équations de Navier-Stokes</td>
      <td>243</td>
    </tr>
    <tr>
      <td>A.1 Introduction</td>
      <td>243</td>
    </tr>
    <tr>
      <td>A.2 Equations instantanées</td>
      <td>243</td>
    </tr>
    <tr>
      <td>A.3 Moyennes employées</td>
      <td>244</td>
    </tr>
    <tr>
      <td>A.4 Moyenne des équations</td>
      <td>245</td>
    </tr>
    <tr>
      <td>A.5 Simplification du système d'équations</td>
      <td>246</td>
    </tr>
    <tr>
      <td>A.6 Conclusion</td>
      <td>248</td>
    </tr>
    <tr>
      <td>B Adimensionnalisation des équations de Navier-Stokes</td>
      <td>251</td>
    </tr>
    <tr>
      <td>C Formulation volumes finis</td>
      <td>255</td>
    </tr>
    <tr>
      <td>C.1 Introduction</td>
      <td>255</td>
    </tr>
    <tr>
      <td>C.2 Formulation intégrale en maillage mobile</td>
      <td>255</td>
    </tr>
    <tr>
      <td>C.2.1 Cas d'une équation scalaire</td>
      <td>255</td>
    </tr>
    <tr>
      <td>C.2.2 Cas d'une équation vectorielle</td>
      <td>257</td>
    </tr>
    <tr>
      <td>C.3 Changement de variables d'espace</td>
      <td>257</td>
    </tr>
    <tr>
      <td>C.4 Hypothèses des volumes finis</td>
      <td>259</td>
    </tr>
    <tr>
      <td>C.4.1 Première approche</td>
      <td>259</td>
    </tr>
    <tr>
      <td>C.4.2 Deuxième approche</td>
      <td>260</td>
    </tr>
    <tr>
      <td>C.4.3 Bilan</td>
      <td>261</td>
    </tr>
    <tr>
      <td>C.5 Equations discrétisées</td>
      <td>261</td>
    </tr>
    <tr>
      <td>C.6 Comparaison aux différences finies</td>
      <td>262</td>
    </tr>
    <tr>
      <td>D Géométrie différentielle en coordonnées curvilignes</td>
      <td>265</td>
    </tr>
    <tr>
      <td>D.1 Introduction</td>
      <td>265</td>
    </tr>
    <tr>
      <td>D.2 Base covariante</td>
      <td>265</td>
    </tr>
    <tr>
      <td>D.3 Base contravariante</td>
      <td>267</td>
    </tr>
    <tr>
      <td>D.4 Propriétés</td>
      <td>268</td>
    </tr>
    <tr>
      <td>D.5 Eléments différentiels</td>
      <td>268</td>
    </tr>
    <tr>
      <td>D.5.1 Elément différentiel linéaire</td>
      <td>268</td>
    </tr>
    <tr>
      <td>D.5.2 Elément différentiel surfacique</td>
      <td>269</td>
    </tr>
    <tr>
      <td>D.5.3 Elément différentiel volumique</td>
      <td>271</td>
    </tr>
    <tr>
      <td>D.6 Opérateurs différentiels</td>
      <td>271</td>
    </tr>
    <tr>
      <td>D.6.1 Divergence</td>
      <td>271</td>
    </tr>
    <tr>
      <td>D.6.2 Rotationnel</td>
      <td>271</td>
    </tr>
    <tr>
      <td>D.6.3 Gradient</td>
      <td>272</td>
    </tr>
    <tr>
      <td>D.6.4 Laplacien</td>
      <td>272</td>
    </tr>
    <tr>
      <td>D.7 Dérivées temporelles</td>
      <td>273</td>
    </tr>
    <tr>
      <td>D.8 Calcul numérique des métriques</td>
      <td>274</td>
    </tr>
    <tr>
      <td>D.8.1 Calcul des métriques en un noeud</td>
      <td>274</td>
    </tr>
  </tbody>
</table>
<table><tbody><tr><td>D.8.2</td><td>Calcul des métriques à une interface</td><td>275</td></tr><tr><td>E</td><td>Cas tests référencés</td><td>277</td></tr><tr><td>E.1</td><td>Ecoulement en tube à choc mono-dimensionnel</td><td>277</td></tr><tr><td>E.1.1</td><td>Description de la géométrie</td><td>277</td></tr><tr><td>E.1.2</td><td>Description des conditions initiales</td><td>277</td></tr><tr><td>E.1.3</td><td>Solution analytique</td><td>278</td></tr><tr><td>E.2</td><td>Interactions d'ondes mono-dimensionnelles</td><td>278</td></tr><tr><td>E.2.1</td><td>Description de la géométrie</td><td>278</td></tr><tr><td>E.2.2</td><td>Description des conditions initiales</td><td>278</td></tr><tr><td>E.2.3</td><td>Solution analytique</td><td>279</td></tr><tr><td>E.3</td><td>Interactions d'ondes bi-dimensionnelles</td><td>280</td></tr><tr><td>E.3.1</td><td>Description de la géométrie</td><td>280</td></tr><tr><td>E.3.2</td><td>Description des conditions initiales</td><td>280</td></tr><tr><td>E.4</td><td>Tuyère supersonique</td><td>281</td></tr><tr><td>E.4.1</td><td>Description de la géométrie</td><td>281</td></tr><tr><td>E.4.2</td><td>Description des conditions aux limites</td><td>281</td></tr><tr><td>E.5</td><td>Tuyère transsonique</td><td>283</td></tr><tr><td>E.5.1</td><td>Description de la géométrie</td><td>283</td></tr><tr><td>E.5.2</td><td>Description des conditions aux limites</td><td>284</td></tr><tr><td>E.5.3</td><td>Configuration instationnaire</td><td>285</td></tr><tr><td colspan="2">Bibliographie</td><td>292</td></tr></tbody></table>
Liste des figures

<table>
  <tbody>
    <tr>
      <td>3.1 Cellule d'intégration</td>
      <td>15</td>
    </tr>
    <tr>
      <td>3.2 Changement de variables d'espace</td>
      <td>17</td>
    </tr>
    <tr>
      <td>3.3 Stockage des variables conservatives aux noeuds</td>
      <td>18</td>
    </tr>
    <tr>
      <td>3.4 Stockage des variables conservatives au centre des mailles</td>
      <td>19</td>
    </tr>
    <tr>
      <td>3.5 Approche MUSCL</td>
      <td>27</td>
    </tr>
    <tr>
      <td>3.6 Tube à choc au premier ordre. Niveau de densité</td>
      <td>29</td>
    </tr>
    <tr>
      <td>3.7 Tube à choc au second ordre. Niveau de densité</td>
      <td>30</td>
    </tr>
    <tr>
      <td>3.8 Tube à choc au premier ordre. Niveau d'énergie totale</td>
      <td>30</td>
    </tr>
    <tr>
      <td>3.9 Tube à choc au deuxième ordre. Niveau d'énergie totale</td>
      <td>31</td>
    </tr>
    <tr>
      <td>3.10 Intéractions d'ondes en tube à choc. Niveau de densité</td>
      <td>32</td>
    </tr>
    <tr>
      <td>3.11 Tube à choc. Valeur initiale du limiteur de Mulder</td>
      <td>34</td>
    </tr>
    <tr>
      <td>3.12 Tube à choc. Limiteur de Mulder. ε = 10<sup>-6</sup></td>
      <td>35</td>
    </tr>
    <tr>
      <td>3.13 Tube à choc. Limiteur de Mulder. ε = 10<sup>-12</sup></td>
      <td>35</td>
    </tr>
    <tr>
      <td>3.14 Limiteur en cosinus</td>
      <td>36</td>
    </tr>
    <tr>
      <td>3.15 Limiteur en tangente hyperbolique</td>
      <td>37</td>
    </tr>
    <tr>
      <td>3.16 Choc cylindrique. Niveau de masse volumique en maillage radial</td>
      <td>38</td>
    </tr>
    <tr>
      <td>3.17 Choc cylindrique. Niveau de masse volumique en maillage cartésien</td>
      <td>38</td>
    </tr>
    <tr>
      <td>3.18 Choc cylindrique. Niveaux de masse volumique et de pression</td>
      <td>41</td>
    </tr>
    <tr>
      <td>3.19 Tuyère transsonique. Profils longitudinaux de pression totale</td>
      <td>43</td>
    </tr>
    <tr>
      <td>3.20 Tuyère transsonique. Profils longitudinaux de nombre de Mach</td>
      <td>43</td>
    </tr>
    <tr>
      <td>3.21 Tuyère transsonique. Profils longitudinaux de température totale</td>
      <td>44</td>
    </tr>
    <tr>
      <td>3.22 Tuyère faiblement subsonique. Nombre de Mach</td>
      <td>48</td>
    </tr>
    <tr>
      <td>3.23 Tuyère faiblement subsonique. Niveau de température totale. Van Leer</td>
      <td>48</td>
    </tr>
    <tr>
      <td>3.24 Tuyère faiblement subsonique. Niveau de température totale. Hänel</td>
      <td>49</td>
    </tr>
    <tr>
      <td>3.25 Ecoulement de Poiseuille. Convergence vers l'état stationnaire</td>
      <td>51</td>
    </tr>
    <tr>
      <td>3.26 Ecoulement de Poiseuille. Erreur sur la pression statique</td>
      <td>52</td>
    </tr>
    <tr>
      <td>3.27 Ecoulement de Poiseuille. Erreur sur la vitesse axiale</td>
      <td>52</td>
    </tr>
    <tr>
      <td>3.28 Molécule de calcul des flux visqueux</td>
      <td>57</td>
    </tr>
    <tr>
      <td>3.29 Evolution de la longueur de mélange</td>
      <td>58</td>
    </tr>
    <tr>
      <td>3.30 Modèle de Michel. Profils de vitesse</td>
      <td>61</td>
    </tr>
    <tr>
      <td>3.31 Modèle de Michel. Profil de vitesse amont</td>
      <td>62</td>
    </tr>
    <tr>
      <td>3.32 Modèle de Michel. Profil de vitesse aval</td>
      <td>62</td>
    </tr>
    <tr>
      <td>3.33 Modèle de Michel. Contrainte de cisaillement amont</td>
      <td>63</td>
    </tr>
  </tbody>
</table>
<table>
  <tbody>
    <tr>
      <td>3.34</td>
      <td>Modèle de Michel. Contrainte de cisaillement aval.</td>
      <td>63</td>
    </tr>
    <tr>
      <td>3.35</td>
      <td>Répartition statistique du résidu à convergence</td>
      <td>66</td>
    </tr>
    <tr>
      <td>4.1</td>
      <td>Molécule de calcul pour un noeud du coeur</td>
      <td>70</td>
    </tr>
    <tr>
      <td>4.2</td>
      <td>Molécule de calcul pour un noeud juste intérieur</td>
      <td>70</td>
    </tr>
    <tr>
      <td>4.3</td>
      <td>Molécule de calcul pour un noeud frontière</td>
      <td>70</td>
    </tr>
    <tr>
      <td>4.4</td>
      <td>Plans de maillage supplémentaires</td>
      <td>72</td>
    </tr>
    <tr>
      <td>4.5</td>
      <td>Symétrie plane</td>
      <td>74</td>
    </tr>
    <tr>
      <td>4.6</td>
      <td>Tuyère symétrique</td>
      <td>74</td>
    </tr>
    <tr>
      <td>4.7</td>
      <td>Périodicité.</td>
      <td>76</td>
    </tr>
    <tr>
      <td>4.8</td>
      <td>Grille d'aubes</td>
      <td>77</td>
    </tr>
    <tr>
      <td>4.9</td>
      <td>Périodicité angulaire</td>
      <td>78</td>
    </tr>
    <tr>
      <td>4.10</td>
      <td>Inducteur transsonique</td>
      <td>79</td>
    </tr>
    <tr>
      <td>4.11</td>
      <td>Extrapolation</td>
      <td>80</td>
    </tr>
    <tr>
      <td>4.12</td>
      <td>Construction d'un maillage en O</td>
      <td>83</td>
    </tr>
    <tr>
      <td>4.13</td>
      <td>Maillage en O</td>
      <td>84</td>
    </tr>
    <tr>
      <td>4.14</td>
      <td>Ecoulement autour d'un cylindre. Nombre de Mach</td>
      <td>84</td>
    </tr>
    <tr>
      <td>4.15</td>
      <td>Ecoulement autour d'un cylindre. Maillage</td>
      <td>85</td>
    </tr>
    <tr>
      <td>4.16</td>
      <td>Maillage en C</td>
      <td>85</td>
    </tr>
    <tr>
      <td>4.17</td>
      <td>Construction d'un multi-domaines</td>
      <td>86</td>
    </tr>
    <tr>
      <td>4.18</td>
      <td>Raccord entre deux domaines</td>
      <td>87</td>
    </tr>
    <tr>
      <td>4.19</td>
      <td>Choc sortant à l'atmosphère. Pression statique en bout de tube</td>
      <td>97</td>
    </tr>
    <tr>
      <td>4.20</td>
      <td>Choc sortant à l'atmosphère. Masse volumique en bout de tube</td>
      <td>97</td>
    </tr>
    <tr>
      <td>4.21</td>
      <td>Frontières libres</td>
      <td>98</td>
    </tr>
    <tr>
      <td>4.22</td>
      <td>Angle de calage dans le plan (x<sub>i</sub>, x<sub>j</sub>)</td>
      <td>101</td>
    </tr>
    <tr>
      <td>4.23</td>
      <td>Frontière glissante adiabatique</td>
      <td>102</td>
    </tr>
    <tr>
      <td>4.24</td>
      <td>Cylindre à M = 0.3. Erreur sur la pression totale</td>
      <td>105</td>
    </tr>
    <tr>
      <td>4.25</td>
      <td>Cylindre à M = 0.3. Erreur sur la température totale</td>
      <td>107</td>
    </tr>
    <tr>
      <td>4.26</td>
      <td>Frontière adhérente isotherme. Profil de vitesse</td>
      <td>109</td>
    </tr>
    <tr>
      <td>4.27</td>
      <td>Frontière adhérente isotherme. Température statique</td>
      <td>111</td>
    </tr>
    <tr>
      <td>4.28</td>
      <td>Frontière adhérente à flux de chaleur imposé.</td>
      <td>113</td>
    </tr>
    <tr>
      <td>4.29</td>
      <td>Extrapolation assurant l'adiabaticité</td>
      <td>116</td>
    </tr>
    <tr>
      <td>4.30</td>
      <td>Paroi glissante</td>
      <td>124</td>
    </tr>
    <tr>
      <td>4.31</td>
      <td>Arête glissante/glissante</td>
      <td>126</td>
    </tr>
    <tr>
      <td>4.32</td>
      <td>Coin glissant/glissant/sortie subsonique</td>
      <td>128</td>
    </tr>
    <tr>
      <td>4.33</td>
      <td>Interaction d'un choc et d'une frontière de non reflexion</td>
      <td>132</td>
    </tr>
    <tr>
      <td>4.34</td>
      <td>Interaction d'un choc et d'une frontière libre.</td>
      <td>133</td>
    </tr>
    <tr>
      <td>4.35</td>
      <td>Raccord de différentes conditions en un bord d'attaque</td>
      <td>135</td>
    </tr>
    <tr>
      <td>5.1</td>
      <td>Adaptation statique. Solution initiale.</td>
      <td>138</td>
    </tr>
    <tr>
      <td>5.2</td>
      <td>Adaptation statique. Limiteur après seuillage</td>
      <td>139</td>
    </tr>
    <tr>
      <td>5.3</td>
      <td>Adaptation statique. Limiteur après fenêtrage</td>
      <td>140</td>
    </tr>
    <tr>
      <td>5.4</td>
      <td>Adaptation statique. Limiteur après lissage</td>
      <td>140</td>
    </tr>
  </tbody>
</table>
<table>
  <tbody>
    <tr>
      <td>5.5</td>
      <td>Adaptation statique. Limiteur après intégration</td>
      <td>141</td>
    </tr>
    <tr>
      <td>5.6</td>
      <td>Adaptation statique. Solution finale</td>
      <td>143</td>
    </tr>
    <tr>
      <td>5.7</td>
      <td>Adaptation statique. Choc bi-dimensionnel</td>
      <td>143</td>
    </tr>
    <tr>
      <td>5.8</td>
      <td>Adaptation dynamique</td>
      <td>145</td>
    </tr>
    <tr>
      <td>6.1</td>
      <td>Tube à choc. Faisceau d'ondes</td>
      <td>148</td>
    </tr>
    <tr>
      <td>6.2</td>
      <td>Tube à choc. Niveau de masse volumique. Premier ordre</td>
      <td>149</td>
    </tr>
    <tr>
      <td>6.3</td>
      <td>Tube à choc. Niveau de pression statique. Premier ordre</td>
      <td>150</td>
    </tr>
    <tr>
      <td>6.4</td>
      <td>Tube à choc. Niveau de masse volumique. Troisième ordre</td>
      <td>150</td>
    </tr>
    <tr>
      <td>6.5</td>
      <td>Tube à choc. Niveau de pression statique. Troisième ordre</td>
      <td>151</td>
    </tr>
    <tr>
      <td>6.6</td>
      <td>Tube à choc. Profil de masse volumique. Premier ordre</td>
      <td>152</td>
    </tr>
    <tr>
      <td>6.7</td>
      <td>Tube à choc. Profil d'énergie totale. Premier ordre</td>
      <td>152</td>
    </tr>
    <tr>
      <td>6.8</td>
      <td>Tube à choc. Profil de masse volumique. Troisième ordre</td>
      <td>153</td>
    </tr>
    <tr>
      <td>6.9</td>
      <td>Tube à choc. Profil d'énergie totale. Troisième ordre</td>
      <td>153</td>
    </tr>
    <tr>
      <td>6.10</td>
      <td>Interactions d'ondes 1D. Faisceau d'ondes</td>
      <td>154</td>
    </tr>
    <tr>
      <td>6.11</td>
      <td>Interactions d'ondes 1D. Niveau de masse volumique</td>
      <td>155</td>
    </tr>
    <tr>
      <td>6.12</td>
      <td>Interactions d'ondes 1D. Niveau de pression statique</td>
      <td>156</td>
    </tr>
    <tr>
      <td>6.13</td>
      <td>Interactions d'ondes 1D. Niveau d'énergie totale</td>
      <td>157</td>
    </tr>
    <tr>
      <td>6.14</td>
      <td>Interactions d'ondes 1D. Limiteur associé à l'énergie totale</td>
      <td>159</td>
    </tr>
    <tr>
      <td>6.15</td>
      <td>Interactions d'ondes 1D. Profil de masse volumique à t=0.20s</td>
      <td>161</td>
    </tr>
    <tr>
      <td>6.16</td>
      <td>Interactions d'ondes 1D. Profil de masse volumique à t=0.40s</td>
      <td>161</td>
    </tr>
    <tr>
      <td>6.17</td>
      <td>Interactions d'ondes 1D. Profil de masse volumique à t=0.60s</td>
      <td>162</td>
    </tr>
    <tr>
      <td>6.18</td>
      <td>Interactions d'ondes 1D. Profil de masse volumique à t=0.70s</td>
      <td>162</td>
    </tr>
    <tr>
      <td>6.19</td>
      <td>Interactions d'ondes 2D. Niveau de masse volumique à t=0.00s</td>
      <td>164</td>
    </tr>
    <tr>
      <td>6.20</td>
      <td>Interactions d'ondes 2D. Niveau de masse volumique à t=0.20s</td>
      <td>165</td>
    </tr>
    <tr>
      <td>6.21</td>
      <td>Interactions d'ondes 2D. Niveau de masse volumique à t=0.30s</td>
      <td>165</td>
    </tr>
    <tr>
      <td>6.22</td>
      <td>Interactions d'ondes 2D. Niveau de masse volumique à t=0.40s</td>
      <td>166</td>
    </tr>
    <tr>
      <td>6.23</td>
      <td>Interactions d'ondes 2D. Niveau de masse volumique à t=0.50s</td>
      <td>166</td>
    </tr>
    <tr>
      <td>6.24</td>
      <td>Interactions d'ondes 2D. Niveau de masse volumique à t=0.70s</td>
      <td>167</td>
    </tr>
    <tr>
      <td>6.25</td>
      <td>Interactions d'ondes 2D. Profils radiaux de densité à 0.10s d'intervalle</td>
      <td>168</td>
    </tr>
    <tr>
      <td>6.26</td>
      <td>Interactions d'ondes 2D. Profils radiaux de densité à 0.10s d'intervalle</td>
      <td>169</td>
    </tr>
    <tr>
      <td>6.27</td>
      <td>Interactions d'ondes 2D. Profils radiaux de densité à 0.10s d'intervalle</td>
      <td>170</td>
    </tr>
    <tr>
      <td>6.28</td>
      <td>Interactions d'ondes 2D. Profils radiaux de densité à 0.10s d'intervalle</td>
      <td>170</td>
    </tr>
    <tr>
      <td>6.29</td>
      <td>Tuyère à rafales. Géométrie</td>
      <td>171</td>
    </tr>
    <tr>
      <td>6.30</td>
      <td>Tuyère à rafales. Nombre de Mach</td>
      <td>173</td>
    </tr>
    <tr>
      <td>6.31</td>
      <td>Tuyère à rafales. Nombre de Mach sur la paroi inférieure</td>
      <td>173</td>
    </tr>
    <tr>
      <td>6.32</td>
      <td>Tuyère à rafales. Nombre de Mach sur la paroi inférieure</td>
      <td>175</td>
    </tr>
    <tr>
      <td>6.33</td>
      <td>Tuyère à rafales. Vitesse longitudinale</td>
      <td>177</td>
    </tr>
    <tr>
      <td>6.34</td>
      <td>Tuyère à rafales. Vitesse transversale</td>
      <td>179</td>
    </tr>
    <tr>
      <td>6.35</td>
      <td>Tuyère à rafales. Erreur relative sur la pression totale</td>
      <td>179</td>
    </tr>
    <tr>
      <td>6.36</td>
      <td>Tuyère à rafales. Erreur relative sur la température totale</td>
      <td>181</td>
    </tr>
    <tr>
      <td>6.37</td>
      <td>Tuyère transsonique. Géométrie expérimentale</td>
      <td>183</td>
    </tr>
    <tr>
      <td>6.38</td>
      <td>Tuyère transsonique. Plan d'étude bi-dimensionnel</td>
      <td>184</td>
    </tr>
  </tbody>
</table>
<table>
  <tbody>
    <tr>
      <td>6.39</td>
      <td>Tuyère transsonique. Maillage bi-dimensionnel non visqueux</td>
      <td>184</td>
    </tr>
    <tr>
      <td>6.40</td>
      <td>Tuyère transsonique. Nombre de Mach</td>
      <td>185</td>
    </tr>
    <tr>
      <td>6.41</td>
      <td>Tuyère transsonique. Pression statique non visqueux</td>
      <td>187</td>
    </tr>
    <tr>
      <td>6.42</td>
      <td>Tuyère transsonique. Pression totale non visqueux</td>
      <td>188</td>
    </tr>
    <tr>
      <td>6.43</td>
      <td>Tuyère transsonique. Maillage bi-dimensionnel visqueux</td>
      <td>189</td>
    </tr>
    <tr>
      <td>6.44</td>
      <td>Tuyère transsonique. Nombre de Mach non visqueux</td>
      <td>190</td>
    </tr>
    <tr>
      <td>6.45</td>
      <td>Tuyère transsonique. Nombre de Mach laminaire. Schéma de Van Leer</td>
      <td>190</td>
    </tr>
    <tr>
      <td>6.46</td>
      <td>Tuyère transsonique. Nombre de Mach turbulent. Schéma de Van Leer</td>
      <td>191</td>
    </tr>
    <tr>
      <td>6.47</td>
      <td>Tuyère transsonique. Nombre de Mach turbulent. Schéma mixte ε' = 0.5191</td>
      <td></td>
    </tr>
    <tr>
      <td>6.48</td>
      <td>Tuyère transsonique. Nombre de Mach turbulent. Schéma mixte ε' = 0.1192</td>
      <td></td>
    </tr>
    <tr>
      <td>6.49</td>
      <td>Tuyère transsonique. Nombre de Mach turbulent. Schéma de Roe</td>
      <td>192</td>
    </tr>
    <tr>
      <td>6.50</td>
      <td>Tuyère transsonique. Pression totale laminaire. Schéma de Van Leer</td>
      <td>193</td>
    </tr>
    <tr>
      <td>6.51</td>
      <td>Tuyère transsonique. Pression totale</td>
      <td>194</td>
    </tr>
    <tr>
      <td>6.52</td>
      <td>Tuyère transsonique. Pression totale turbulent. Schéma de Van Leer</td>
      <td>194</td>
    </tr>
    <tr>
      <td>6.53</td>
      <td>Tuyère transsonique. Pression totale turbulent. Schéma de Roe</td>
      <td>194</td>
    </tr>
    <tr>
      <td>6.54</td>
      <td>Tuyère transsonique. Pression statique</td>
      <td>195</td>
    </tr>
    <tr>
      <td>6.55</td>
      <td>Tuyère transsonique. Pression statique</td>
      <td>197</td>
    </tr>
    <tr>
      <td>6.56</td>
      <td>Tuyère transsonique. Pression statique turbulent en paroi</td>
      <td>199</td>
    </tr>
    <tr>
      <td>6.57</td>
      <td>Tuyère transsonique. Pression statique turbulent en paroi</td>
      <td>199</td>
    </tr>
    <tr>
      <td>6.58</td>
      <td>Tuyère transsonique. Pression statique turbulent sur l'axe</td>
      <td>200</td>
    </tr>
    <tr>
      <td>6.59</td>
      <td>Tuyère transsonique. Vitesse au col visqueux</td>
      <td>201</td>
    </tr>
    <tr>
      <td>6.60</td>
      <td>Tuyère transsonique. Vitesse en sortie visqueux</td>
      <td>201</td>
    </tr>
    <tr>
      <td>6.61</td>
      <td>Ecoulement instationnaire. Pression statique pour P<sub>8</sub> = 0.6520</td>
      <td>203</td>
    </tr>
    <tr>
      <td>6.62</td>
      <td>Ecoulement instationnaire. Domaine d'étude</td>
      <td>205</td>
    </tr>
    <tr>
      <td>6.63</td>
      <td>Ecoulement instationnaire. Position instantanée du choc</td>
      <td>205</td>
    </tr>
    <tr>
      <td>6.64</td>
      <td>Ecoulement instationnaire. Pression statique instantanée de P<sub>1</sub> à P<sub>8</sub></td>
      <td>206</td>
    </tr>
    <tr>
      <td>6.65</td>
      <td>Ecoulement instationnaire. Pression statique instantanée en P<sub>3</sub></td>
      <td>207</td>
    </tr>
    <tr>
      <td>6.66</td>
      <td>Ecoulement instationnaire. Pression statique instantanée en P<sub>5</sub></td>
      <td>208</td>
    </tr>
    <tr>
      <td>6.67</td>
      <td>Grille d'aubes supersonique. Géométrie</td>
      <td>209</td>
    </tr>
    <tr>
      <td>6.68</td>
      <td>Grille d'aubes supersonique. Bord d'attaque</td>
      <td>210</td>
    </tr>
    <tr>
      <td>6.69</td>
      <td>Grille d'aubes supersonique. Bord de fuite</td>
      <td>210</td>
    </tr>
    <tr>
      <td>6.70</td>
      <td>Grille d'aubes supersonique. Nombre de Mach</td>
      <td>211</td>
    </tr>
    <tr>
      <td>6.71</td>
      <td>Grille d'aubes supersonique. Pression statique</td>
      <td>213</td>
    </tr>
    <tr>
      <td>6.72</td>
      <td>Grille d'aubes supersonique. Pression statique</td>
      <td>215</td>
    </tr>
    <tr>
      <td>6.73</td>
      <td>Grille d'aubes supersonique. Pression statique</td>
      <td>217</td>
    </tr>
    <tr>
      <td>6.74</td>
      <td>Grille d'aubes supersonique. Pression totale</td>
      <td>219</td>
    </tr>
    <tr>
      <td>6.75</td>
      <td>Grille d'aubes supersonique. Température totale</td>
      <td>221</td>
    </tr>
    <tr>
      <td>6.76</td>
      <td>Grille d'aubes supersonique. Géométrie grossière</td>
      <td>223</td>
    </tr>
    <tr>
      <td>6.77</td>
      <td>Grille d'aubes supersonique. Nombre de Mach en géométrie grossière</td>
      <td>225</td>
    </tr>
    <tr>
      <td>6.78</td>
      <td>Grille d'aubes supersonique. Pression totale en géométrie grossière</td>
      <td>227</td>
    </tr>
    <tr>
      <td>6.79</td>
      <td>Coude subsonique. Géométrie</td>
      <td>230</td>
    </tr>
    <tr>
      <td>6.80</td>
      <td>Coude subsonique. Vitesses</td>
      <td>232</td>
    </tr>
  </tbody>
</table>
<table>
  <tbody>
    <tr>
      <td>6.81</td>
      <td>Coude subsonique. Pression totale relative</td>
      <td>233</td>
    </tr>
    <tr>
      <td>6.82</td>
      <td>Coude subsonique. Pression statique</td>
      <td>237</td>
    </tr>
    <tr>
      <td>C.1</td>
      <td>Cellule d'intégration</td>
      <td>256</td>
    </tr>
    <tr>
      <td>C.2</td>
      <td>Changement de variables d'espace</td>
      <td>258</td>
    </tr>
    <tr>
      <td>C.3</td>
      <td>Stockage des variables conservatives aux noeuds</td>
      <td>261</td>
    </tr>
    <tr>
      <td>C.4</td>
      <td>Stockage des variables conservatives au centre des mailles</td>
      <td>262</td>
    </tr>
    <tr>
      <td>D.1</td>
      <td>Coordonnées cartésiennes et curvilignes</td>
      <td>266</td>
    </tr>
    <tr>
      <td>D.2</td>
      <td>Base covariante</td>
      <td>266</td>
    </tr>
    <tr>
      <td>D.3</td>
      <td>Base contravariante</td>
      <td>267</td>
    </tr>
    <tr>
      <td>D.4</td>
      <td>Elément différentiel linéaire</td>
      <td>269</td>
    </tr>
    <tr>
      <td>D.5</td>
      <td>Elément différentiel surfacique</td>
      <td>270</td>
    </tr>
    <tr>
      <td>D.6</td>
      <td>Elément différentiel volumique</td>
      <td>270</td>
    </tr>
    <tr>
      <td>D.7</td>
      <td>Calcul des métriques en un noeud</td>
      <td>274</td>
    </tr>
    <tr>
      <td>D.8</td>
      <td>Calcul des métriques à une interface</td>
      <td>275</td>
    </tr>
    <tr>
      <td>E.1</td>
      <td>Interactions d'ondes 1D. Faisceau d'ondes</td>
      <td>279</td>
    </tr>
    <tr>
      <td>E.2</td>
      <td>Tuyère supersonique. Description des parois</td>
      <td>281</td>
    </tr>
    <tr>
      <td>E.3</td>
      <td>Tuyère transsonique. Description des parois</td>
      <td>283</td>
    </tr>
  </tbody>
</table>
Liste des tableaux

<table><tr><td>3.1</td><td>Composantes du flux aux noeuds et aux interfaces</td><td>45</td></tr><tr><td>4.1</td><td>Nombre de conditions physiques pour une frontière libre</td><td>99</td></tr><tr><td>E.1</td><td>Interactions d'ondes 1D. Solution analytique</td><td>279</td></tr><tr><td>E.2</td><td>Interactions d'ondes 1D. Vitesse absolue des ondes</td><td>280</td></tr><tr><td>E.3</td><td>Interactions d'ondes 1D. Date et position des interactions</td><td>280</td></tr><tr><td>E.4</td><td>Tuyère supersonique. Description des parois</td><td>282</td></tr><tr><td>E.5</td><td>Tuyère transsonique. Position des sondes de pression statique</td><td>285</td></tr></table>
# Chapitre 1

# Introduction

Le développement des turbomachines actuelles nécessite une compréhension de plus en plus fine des écoulements qui y siègent. De nombreux outils d'analyse et de multiples campagnes d'essais ont permis, durant ces dernières décennies, d'en approfondir la connaissance.

Le mouvement d'un fluide au passage des différentes roues, fixes ou mobiles, est extrêmement complexe. L'écoulement est généralement visqueux, souvent *turbulent*. Les parois le guidant étant des surfaces présentant de fortes variations de courbure, il en résulte de nombreuses structures secondaires, tels les tourbillons de passsage, de coin, de bord d'attaque, ... S'y ajoutent le sillage de chaque aube. Ainsi, l'écoulement est *tri-dimensionnel, fortement inhomogène*. De plus, l'accroissement des vitesses de rotation, couplé aux effets de compressibilité, génère un régime transsonique, voir supersonique, s'accompagnant de chocs, d'ondes de détente, ou encore de surfaces de contact, qui se propagent, se réfléchissent, intéragissent entre elles ou avec les parois et les couches limites, ... Si bien que l'écoulement est le siège *d'importantes discontinuités*. D'autre part, la présence de parties fixes et d'autres mobiles, les variations de la vitesse de rotation, mais aussi les fluctuations des conditions en entrée ou en sortie de machine, sont à l'origine de phénomènes instationnaires. Enfin, des interactions fluide/structure, traduisant le couplage potentiel entre le mouvement du fluide et le déplacement des aubes, entrainent des déformations de géométrie.

La prise en compte de tous ces phénomènes lors de l'analyse ou du dessin d'une machine, est un projet très ambitieux. Cependant, une connaissance étendue de chacun, et de son influence sur les autres, permettrait de mieux comprendre les dégradations de performance constatées sous certaines conditions, voir la rupture et la destruction de la machine. Ainsi pourrait être définies, entre autres, des marges de sécurité moins empiriques. Deux méthodes d'investigation s'offre à nous, toutes deux nécessaires et complémentaires. D'une part, l'étude expérimentale permet d'obtenir une image de la "réalité physique". Cependant, elle souffre d'un problème d'acces-sibilité à toutes les informations, en particulier sans perturber l'écoulement, ainsi
que de coûts financiers importants, que la destruction potentielle de l'installation ne peut qu'alourdir. D'autre part, les simulations numériques sont moins chères et autorisent facilement un changement de paramètre, ou de géométrie. Mais elles sont contraintes par les hypothèses de modélisation, et par la capacité des moyens informatiques de calcul.

Cependant, la puissance des ordinateurs augmentant, et de nouvelles techniques voyant le jour, de nombreuses méthodes numériques de plus en plus performantes ont été développées pour simuler cette interaction complexe de phénomènes, en particulier au Laboratoire de Mécanique des Fluide et d'Acoustique de l'Ecole Centrale de Lyon. Ainsi y a débuté en 1987 au sein du groupe Turbomachines, un programme de recherche à long terme visant au calcul d'écoulements tri-dimensionnels turbulents et instationnaires. Dans ce cadre, des travaux fondés sur les théories linéarisées [20] et sur la théorie des disques d'action [18] ont conclu qu'actuellement, pour avancer dans la compréhension des mécanismes instationnaires, il était nécessaire de développer une méthode numérique non linéaire. Suite à quoi une étude [9] précédent le travail faisant l'objet du présent mémoire, a permis de construire un premier code de calcul. Il mettait en oeuvre des techniques éprouvées, mais aussi des développements récents, inemployés jusqu'alors au LMFA. Son but était de valider l'intégration du tout, d'évaluer ses performances, voir d'extrapoler celles des versions futures, et de mettre en évidence les points faibles.

La méthode développée par N. Benbouta reposait sur :

* la résolution des équations de Navier-Stokes moyennées.

* une formulation volumes finis en maillage structuré curviligne, couplée à une approche MUSCL.

* un schéma décentré amont pour l'évaluation des effets convectifs, bâti sur la séparation de flux de Van Leer.

* un schéma aux différences finies centrées pour le calcul des termes diffusifs.

* une approche implicite en temps, précise car ne recourant qu'à peu d'hypothèses simplificatrices.

Ce travail, orienté vers des applications stationnaires, a abouti à un code de calcul opérationnel, qui a permis d'établir la qualité de l'ensemble des méthodes numériques mises en oeuvre, et de leur intégration. Il a cependant révélé aussi un certain nombre de problèmes, dont les principaux sont :

* un mauvais comportement de la séparation de flux de Van Leer à très bas nombre de Mach, en particulier dans les couches limites.

* un traitement trop succinct et peu précis des conditions aux limites.
D'autre part, l'approche implicite en temps s'est révélée trop lourde à ce stade de recherche, principalement en occupation mémoire. Cette technique n'a pas été définitivement rejetée, mais juste jugée peu judicieuse aux vues des développements restant à accomplir.

L'objectif du travail présenté était de réaliser un nouveau code de calcul, bâti sur la structure de son prédécesseur, en tirant profit des enseignements qu'il avait apporté. Par rapport à ce dernier, il présente les caractéristiques supplémentaires suivantes:

* la prise en compte de la turbulence, d'une manière simple dans cette version à l'aide un modèle du type longueur de mélange (paragraphes 2.3, 2.4 et 3.6).

* une orientation résolument instationnaire, en particulier au niveau de la déformation du maillage (paragraphe 3.2 et chapitre 5).

* le retour à un schéma en temps explicite, dans un soucis de simplicité (paragraphe 3.3).

* l'introduction d'un schéma mixte, résultant de l'hybridation en fonction du nombre de Mach d'un schéma décentré amont et d'un schéma centré, qui permet de s'affranchir du mauvais comportement de la séparation de flux de Van Leer, qui est analysée et critiquée par ailleurs (paragraphe 3.4).

* un traitement sophistiqué des conditions aux limites

- géométriques, par introduction entre autres de la périodicité angulaire (paragraphe 4.2).

- physiques, par l'utilisation de deux formulations bâties sur les équations de compatibilité, alliant précision, cohérence avec le schéma interne, et gestion de l'instationnarité (paragraphes 4.3 et 4.4).

Il est à noter que ce n'est pas un code industriel, même si son industrialisation est possible, mais un code de recherche qui a été développé. Ainsi, les critères ayant présidé à sa réalisation ne sont pas des temps de calcul courts et une faible occupation mémoire, mais une grande précision par rapport aux problèmes instationnaires, et la possibilité de mener une étude physique de résultats non dénaturés par l'approche numérique. Un autre aspect important de son développement est sa modularité, qui permet une grande versatilité dans le choix des méthodes de résolution, mais aussi son utilisation en tant que laboratoire numérique où sont testés des techniques destinées à d'autres codes à finalité résolument industrielle quant à eux.

Ce travail a été validé par de nombreux cas tests, simulant des configurations élémentaires représentatives d'un problème spécifique, qui servent d'illustrations tout au long de la présentation de la méthode (chapitres 2 à 5). Des écoulements plus
complexes, proches de situations réelles en turbomachine, ont par ailleurs été calculés et analysés (chapitre 6).

La validation de l'approche présentée a été réalisée par des comparaisons à des solutions analytiques, numériques ou expérimentales, selon celles disponibles. Mais surtout, il a été fait grand usage des résultats fournis par un autre code de calcul développé au LMFA/METRAFLU [25]. Ce code est à la fois très proche de celui présenté ici, car reposant sur un autre schéma décentré amont proposé par Roe, en volumes finis, mais aussi très différent, car bâti sur un maillage non structuré, les effets de la turbulence étant introduits à l'aide d'un modèle $k-\epsilon$. De très nombreuses comparaisons ont pu être ainsi menées de manière poussée, grâce à une bonne connaissance des outils numériques utilisés.

Dans le chapitre 2, le système d'équations gouvernant le mouvement d'un fluide compressible est présenté. Puis les équations de Navier-Stokes y sont moyennées en utilisant la décomposition de Reynolds, avant que l'hypothèse de Boussinesq et un modèle de turbulence à zéro équation ne soient utilisées pour en assurer la fermeture. Elles sont alors adimensionnées.

Le chapitre 3 est consacré quant à lui à leur résolution numérique. Pour ce faire, elles sont tout d'abord intégrées en espace dans le cadre d'une formulation volumes finis en maillage curviligne structuré, mobile et déformable. Puis, l'intégration en temps est faite de manière explicite. Les flux convectifs sont évalués alors par une approche MUSCL, couplée à un schéma mixte, résultant de l'hybridation en fonction du nombre de Mach d'un schéma centré et d'une technique de séparation de flux proposée par Van Leer. Les termes diffusifs sont quant à eux discrétisés par des différences centrées au second ordre. Enfin, les effets de la turbulence sont introduits à travers un modèle algébrique, du type longueur de mélange, développé par Michel.

Le chapitre 4 expose les différents traitements des conditions aux limites qui ont été mis en oeuvre. Tout d'abord est abordée la gestion des conditions géométriques, telles que la périodicité, par l'extension du domaine de calcul et la création de noeuds fictifs supplémentaires. Puis sont présentées deux techniques bâties sur les équations de compatibilité visant à la prise en compte de contraintes extérieures physiques, potentiellement instationnaires.

Dans le chapitre 5, différents mouvements de maillage sont décrits. Ils s'appliquent à intégrer les effets de la rotation, mais aussi à décrire la capacité de la méthode développée à réaliser une adaptation statique, ou dynamique, du maillage.

Le chapitre 6 est consacré à la simulation et à l'analyse de configurations caractéristiques. Ainsi, une première série de résultats s'attache à la propagation et aux interactions d'ondes instationnaires dans un fluide parfait. Un second groupe vise à prédire des écoulements stationnaires confinés, recelant ou non des chocs, pour un fluide parfait ou en régime turbulent. Enfin, deux cas stationnaires proches de la réalité seront exposés. Le premier met en oeuvre une grille d'aubes supersonique, le second un coude tri-dimensionnel à bas nombre de Mach.
Enfin, une conclusion sur la présente étude ainsi que des perspectives d'utilisation et de développement sont données au chapitre 7.
# Chapitre 2

# Modélisation d'un écoulement
compressible turbulent

## 2.1 Introduction

Le système d'équations différentielles gouvernant le mouvement d'un fluide visqueux compressible est, au cours de ce chapitre, mis sous une forme susceptible d'être discrétisée, puis résolue numériquement par la suite. Dans un premier temps, les équations de Navier-Stokes instantanées, couplées à l'équation d'état des gaz parfaits, sont écrites sous forme conservative. Puis elles sont moyennées statistiquement afin que soient pris en compte les effets de la turbulence. La fermeture du nouveau système ainsi obtenu est assurée alors par un modèle de turbulence à zéro équation. Nous obtiendrons après une dernière étape d'adimensionnalisation, un ensemble d'équations différentielles qui modélisent les écoulements instationnaires compressibles turbulents que nous voulons étudier.

## 2.2 Equations instantanées

Le mouvement d'un fluide visqueux compressible est gouverné par les équations de Navier-Stokes, qui découlent de l'application des principes fondamentaux de mécanique et de thermodynamique. Elles expriment la conservation de la masse, de la quantité de mouvement et de l'énergie totale.

Un point du domaine spatio-temporel en lequel elles s'appliquent est repéré par ses coordonnées cartésiennes et une date, à savoir ($x_1, x_2, x_3, t$). Tout vecteur est projeté sur la base cartésienne ($\vec{e}_1, \vec{e}_2, \vec{e}_3$), qui est indépendante de ce point et du temps.

Le champ aérothermodynamique y est caractérisé par les variables dites conservatives, à savoir la masse volumique $\rho$, la quantité de mouvement $\rho\vec{V}$ et l'énergie totale $\rho E$.

Dans le cadre d'une description eulérienne et en l'absence de force extérieure, les équations de Navier-Stokes s'écrivent sous forme différentielle conservative comme
suit:

Conservation de la masse:

$$ \frac{\partial \rho}{\partial t} + \vec{\nabla} \cdot [\rho \vec{V}] = 0 \qquad (2.1) $$

Conservation de la quantité de mouvement:

$$ \frac{\partial \rho \vec{V}}{\partial t} + \vec{\nabla} \cdot [\rho \vec{V} \otimes \vec{V} + p \vec{I} - \vec{\tau}_l] = \vec{0} \qquad (2.2) $$

où:

- $\gamma$ étant le rapport des chaleurs spécifiques $C_p$ et $C_v$, la pression statique $p$ est explicitée à l'aide de l'équation d'état des gaz parfaits,

$$ p = \rho(\gamma - 1) \left[ E - \frac{1}{2} \vec{V}^2 \right] \qquad (2.3) $$

- $\lambda$ et $\mu$ étant les première et deuxième viscosités du fluide, $\vec{\tau}_l$ est le tenseur des contraintes visqueuses qui, pour un fluide newtonien isotrope en équilibre thermodynamique, est défini par:

$$ \vec{\tau}_l = \lambda \vec{\nabla} \cdot \vec{V} \vec{I} + \mu (\vec{\nabla} \vec{V} + \vec{V} \vec{V}^T) \qquad (2.4) $$

L'hypothèse de Stokes permet de relier $\lambda$ à $\mu$ par:

$$ 3\lambda + 2\mu = 0 \qquad (2.5) $$

Conservation de l'énergie:

$$ \frac{\partial \rho E}{\partial t} + \vec{\nabla} \cdot [(\rho E + p)\vec{V} + \vec{\Phi} - \vec{\tau}_l \cdot \vec{V}] = 0 \qquad (2.6) $$

où $\vec{\Phi}$ est le flux de chaleur qui est exprimé en fonction de la température statique $T$ et de la conductibilité thermique du fluide $\kappa$, par la loi de Fourier,

$$ \vec{\Phi} = -\kappa \vec{\nabla} T \qquad (2.7) $$

L'introduction du nombre de Prandtl $P_r$, et de l'énergie interne spécifique $e$, permet de reformuler ainsi cette équation:

$$ \frac{\partial \rho E}{\partial t} + \vec{\nabla} \cdot \left[ (\rho E + p) \vec{V} - \frac{\gamma \mu}{P_r} \vec{\nabla} e - \vec{\tau}_l \cdot \vec{V} \right] = 0 \qquad (2.8) $$

avec:

$$ P_r = \frac{\mu C_p}{\kappa} \qquad (2.9) $$
$$e \equiv C_v T \qquad (2.10)$$

$$= E - \frac{1}{2} \vec{V}^2 \qquad (2.11)$$

## 2.3 Equations moyennées

Le caractère turbulent de l'écoulement est représenté par le comportement fluctuant dans le temps et a priori désordonné des grandeurs qui le caractérisent. Les échelles temporelles et spatiales intervenant alors sont très disparates, et la résolution directe des équations instantanées de Navier-Stokes vues précédemment, est coûteuse en temps de calcul et en taille mémoire, compte tenu de l'importance des maillages à mettre en oeuvre [66]. C'est pourquoi une décomposition du champ aérothermodynamique en une composante moyenne et une composante fluctuante, est généralement utilisée [31, 33]. Introduite dans les équations instantanées qui sont moyennées statistiquement par la suite, cette décomposition conduit aux équations d'évolution de l'écoulement turbulent moyen. Cette transformation est décrite en détail annexe A.

Les variables définies par unité de masse ($\tilde{V}$, $E$, $e$) ayant été décomposées au sens de Favre, et les autres ($\rho$, $\bar{p}$) au sens de Reynolds, nous obtenons le système d'équations suivant, non fermé, admettant pour inconnues les valeurs moyennes, au sens turbulent, de la masse volumique $\bar{\rho}$, de la quantité de mouvement $\bar{\rho}\tilde{V}$ et de l'énergie totale $\bar{\rho}\tilde{E}$.

$$\frac{\partial \bar{\rho}}{\partial t} + \bar{\nabla} \cdot [\bar{\rho} \tilde{V}] = 0 \qquad (2.12)$$

$$\frac{\partial \bar{\rho} \tilde{V}}{\partial t} + \bar{\nabla} \cdot \left[ \bar{\rho} \tilde{V} \otimes \tilde{V} + \bar{\rho} \vec{I} - \overline{\vec{\tau}_{l+t}} \right] = \bar{0} \qquad (2.13)$$

$$\frac{\partial \bar{\rho} \tilde{E}}{\partial t} + \bar{\nabla} \cdot \left[ (\bar{\rho} \tilde{E} + \bar{p}) \tilde{V} - \left\{ \gamma \left( \frac{\bar{\mu}}{P_r} + \frac{\mu_t}{P_{rt}} \right) \bar{\nabla} \tilde{e} + \overline{\vec{\tau}_{l+t} \cdot \tilde{V}} + \frac{\mu_t}{\sigma_k} \bar{\nabla} \tilde{k} \right\} \right] = 0 \qquad (2.14)$$

avec:

$$\bar{p} = \bar{\rho}(\gamma - 1) \left[ \tilde{E} - \frac{1}{2} \tilde{V}^2 - \tilde{k} \right] \qquad (2.15)$$

$$\overline{\vec{\tau}_{l+t}} = (\bar{\lambda} + \lambda_t) \tilde{V} \cdot \tilde{V} \vec{I} + (\bar{\mu} + \mu_t) (\tilde{V} \tilde{V} + \tilde{V} \tilde{V}^t) - \frac{2}{3} \bar{\rho} \bar{k} \vec{I} \qquad (2.16)$$
## 2.4 Fermeture et modèles de turbulence

Trois variables supplémentaires, l'énergie cinétique turbulente moyenne $\tilde{k}$ et deux viscosités turbulentes $\lambda_t$ et $\mu_t$, sont apparues. De nouvelles équations, algébriques ou différentielles, devront donc être introduites pour évaluer $\tilde{k}$ et $\mu_t$, en supposant que la relation de Stokes relie $\lambda_t$ à $\mu_t$:

$$3\lambda_t + 2\mu_t = 0 \qquad (2.17)$$

Les équations assurant la fermeture du système précédent résultent d'une étape supplémentaire dans la modélisation de la turbulence faisant suite au processus de moyenne des équations instantanées et à l'hypothèse de Boussinesq [34].

Le choix d'un modèle doit être fait en particulier en fonction des écoulements étudiés, de sa facilité de mise en oeuvre, et du niveau de précision recherché dans la description de la turbulence. Dans le cadre de cette étude, nous sommes intéressés par les effets de la turbulence, entre autres au niveau des pertes qu'elle induit, mais non par son détail. D'autre part, une large majorité des modélisations développées l'ont été pour des écoulements stationnaires incompressibles. Et le prise en compte de la variation de densité, et a fortiori des fluctuations temporelles macroscopiques, en est encore à ses débuts. Or il s'agit précisément du type d'écoulement que nous voulons simuler, et il semblerait que les modèles les plus frustres soient les plus fiables dans ce cas. Enfin, nous ne voulons, nous ne pouvons pas investir d'emblée dans un modèle complexe, pour des raisons de coûts de développement, de temps de calcul et d'occupation mémoire. C'est pourquoi notre choix s'est porté vers un modèle du type longueur de mélange, simple, robuste, à défaut d'être très précis.

Nous ferons donc par la suite les hypothèses suivantes. Tout d'abord, l'énergie cinétique turbulente $\tilde{k}$ est négligeable et négligée, en particulier devant la pression statique et l'énergie cinétique moyenne. D'autre part, la viscosité turbulente $\mu_t$ est supposée être une fonction algébrique d'une échelle de longueur $l_m$, résultant d'une relation analytique faisant intervenir la géométrie et le champ aérodynamique, et d'une échelle de vitesse, bâtie sur le gradient de vitesse moyenne. Ainsi,

$$\tilde{k} = 0 \qquad (2.18)$$

$$\mu_t = \mu_t(l_m, \bar{\nabla} \tilde{V}) \qquad (2.19)$$

## 2.5 Système d'équations différentielles à résoudre

Après adimensionnalisation (annexe B), nous obtenons finalement les équations de Navier-Stokes moyennées suivantes, où $R_c$ est le nombre de Reynolds:
$$ \frac{\partial \rho}{\partial t} + \vec{\nabla} \cdot [\rho \vec{V}] = 0 \qquad (2.20) $$

$$ \frac{\partial \rho \vec{V}}{\partial t} + \vec{\nabla} \cdot [\rho \vec{V} \otimes \vec{V} + p \vec{I} - \frac{1}{R_e} \vec{\tau}_{l+t}] = \vec{0} \qquad (2.21) $$

$$ \frac{\partial \rho E}{\partial t} + \vec{\nabla} \cdot \left[ (\rho E + p) \vec{V} - \frac{1}{R_e} \left\{ \gamma \left( \frac{\mu}{P_r} + \frac{\mu_t}{P_{rt}} \right) \vec{V} e + \vec{\tau}_{l+t} \cdot \vec{V} \right\} \right] = 0 \qquad (2.22) $$

avec:

$$ p = \rho(\gamma - 1) \left[ E - \frac{1}{2}\vec{V}^2 \right] \qquad (2.23) $$

$$ = \rho(\gamma - 1)e \qquad (2.24) $$

$$ \vec{\tau}_{l+t} = (\lambda + \lambda_t) \vec{\nabla} \cdot \vec{V} \vec{I} + (\mu + \mu_t) (\vec{V}\vec{V} + \vec{V}\vec{V}^t) \qquad (2.25) $$

$$ \mu_t = \mu_t(l_m, \vec{\nabla}\vec{V}) \qquad (2.26) $$

Cette dernière relation sera explicitée à l'aide du modèle de turbulence retenu (paragraphe 3.6).

## 2.6 Conclusion

Ce chapitre a été consacré à l'édification d'un ensemble d'équations fermé, modélisant le mouvement instationnaire turbulent d'un fluide visqueux compressible. A partir des équations de Navier-Stokes instantanées, une moyenne statistique a été effectuée pour prendre en compte la turbulence. Après fermeture au moyen de l'hypothèse de Boussinesq et d'un modèle de turbulence à zéro équation, elles ont été adimensionnées. Les relations différentielles liant l'évolution temporelle et spatiale des variables conservatives ($\rho, \rho\vec{V}, \rho E$) ont ainsi été établies. Elles doivent main-tenant être discrétisées, avant d'être résolues numériquement. Tel est le sujet du prochain chapitre.
# Chapitre 3

# Résolution numérique des équations d'évolution

## 3.1 Introduction

Le mouvement du fluide simulé est modélisé par un ensemble d'équations aux dérivées partielles en temps et en espace, définies sur un domaine borné continu, dont la constitution à fait l'objet du chapitre précédent. Or le système ainsi obtenu n'admet pas de solution analytique dans de nombreux cas. C'est pourquoi nous nous posons de le résoudre numériquement, c'est à dire que nous ne recherchons le champ aérodynamique qu'en un nombre fini de points de l'espace, à un nombre fini d'instants. Cette solution discrète, approchée du problème initial, est alors d'autant plus précise que la segmentation de l'espace et du temps est fine et se rapproche d'un continuum.

Il existe différentes techniques de découpage spatio-temporel et d'approximation numériques des opérateurs différentiels. La sélection des uns ou des autres est à faire en fonction de la complexité de la géométrie, de l'écoulement, mais aussi du coût de mise en oeuvre et d'exploitation. Notre choix [9] s'est porté sur une formulation volumes finis en maillage structuré mobile et déformable en temps, dans laquelle les effets convectifs sont évalués à l'aide d'une technique de séparation de flux, proposée par Van Leer, couplée à une approche MUSCL; alors qu'un schéma aux différences finies centrées est utilisé pour calculer les termes visqueux et turbulents. Ceci pour la discrétisation spatiale. La discrétisation temporelle repose quant à elle sur un schéma explicite.

Les développements présentés dans ce chapitre concernent un noeud interne du domaine de calcul. Il est en particulier supposé être éloigné des frontières et posséder un voisinage suffisamment étendu pour mener à bien la discrétisation. Sont présentés par la suite les schémas numériques retenus et leur mise en oeuvre chronologique, accompagnés d'une critique de leur comportement.
## 3.2 Formulation volumes finis

L'idée de base de la formulation volumes finis est que l'espace est découpé en volumes élémentaires contigus, sur lesquels sont intégrées les équations différentielles. Les propriétés de l'écoulement sont supposées faiblement variables à l'intérieur de ces cellules d'intégration [54, 16]. Cette approche présente de nombreux avantages. Tout d'abord, elle induit une diminution de un de l'ordre tensoriel en espace du système, puisqu'un degré de dérivation disparaît. D'autre part, nous retrouvons le caractère conservatif des équations d'évolution, qui résultent initialement d'un bilan de masse, de quantité de mouvement ou d'énergie, sur un volume de contrôle. Nous mâtinons ainsi d'un peu de physique notre schéma numérique. Enfin, les discontinuités, telles que les ondes de chocs et les surfaces de contact, se retrouvent naturellement gérées en tant que solutions faibles du problème [55]. Nous tirons partie du fait que, bien que le champ aérodynamique puisse être discontinu, les flux à travers la surface limitant le volume d'intégration restent continus [39, p. 19]. Les développements mathématiques associés à cette transformation sont reproupés annexe C, et seules les grandes lignes directrices en seront exposées par la suite.

### 3.2.1 Intégration en espace

Rappelons que nous recherchons la solution approchée d'une équation différentielle ayant pour forme:

$$ \frac{\partial A}{\partial t} + \nabla \cdot \vec{B} = 0 \qquad (3.1) $$

où:

- A est une grandeur scalaire (resp. vectorielle)

- $\vec{B}$ est une grandeur vectorielle (resp. tensorielle)

Pour cela, cette dernière est intégrée, dans le repère cartésien, sur un volume élémentaire mobile et déformable en temps (figure 3.1). Ceci diffère d'une approche plus classique, où serait utilisé un maillage fixe, ou quasi-stationnaire, i.e. ne se déplaçant pas durant la phase d'évaluation du champ aérodynamique.

Plusieurs avantages découlent de cette déformation dans le temps des cellules d'intégration. Tout d'abord, il est possible d'adapter dynamiquement la taille et la position des mailles en fonction de la solution. Nous contrôlons ainsi la précision, en assurant une minimisation de l'erreur numérique. Cette technique correspond à un raffinage/déraffinage dans l'espace par déplacement des noeuds et non par création de nouveaux éléments. Une gestion informatique plus simple en découle puisque la topologie du maillage ne change pas. Ceci est donc particulièrement adapté en maillages structurés, à topologie intrinsèque figée. Dans la même veine d'idées, des perturbations, telles que des ondes de choc ou des sillages, peuvent être suivies lors de leur périgrination. Les noeuds ont en effet toute latitude de se déplacer
<figure><img src="image_3.png" /><figcaption>Figure 3.1: Cellule d’intégration</figcaption></figure>

avec elles. D'autre part, il est ainsi possible de gérer la rotation d'une roue sans introduire de changement de repère en coordonnées cylindriques, qui suppose une vitesse angulaire uniforme, interdisant toute simulation d'une montée progressive en régime. Il suffit de faire tourner tout le domaine de calcul. De plus, la prise en compte d'une déformation de la géométrie, résultant par exemple de la vibration des aubes, nécessite que soit introduite à la base dans les équations la notion de déformation et de déplacement de la surface physique du volume d'intégration. Ainsi se trouve grandement facilité le couplage avec un programme de calcul de structure si nous voulons réaliser une simulation de couplage fluide/structure. Un autre avantage découle du fait que la superposition de plusieurs mouvements, dus à une adaptation, une rotation, une déformation, ..., pourront être traités comme un seul mouvement résultant, peu important son origine et sa composition. Il est donc facile d'enrichir ces déplacements des noeuds. Enfin, contrairement à une approche classique, il n'est pas nécessaire d'interpoler la solution d'un instant à l'autre, le maillage et le champ aérodynamique évoluant en phase.

Cependant, il y a un certain nombre d'inconvénients à cette technique. Apparaissent entre autre de nouveaux termes dans les équations dus à la vitesse de la surface des cellules $\overline{W}$. Ils engendrent un faible surcoût en temps de calcul au niveau de l'évaluation du champ aérodynamique, mais la gestion du mouvement du maillage peut être chère, atteignant jusqu'à quelques 10% du CPU total. D'autre part, une contrainte supplémentaire existe lors de la discrétisation en temps, car il faut prendre en compte simultanément l'évolution de la solution et de la géométrie.

Nous obtenons alors l'équation suivante:
$$ \frac{d}{dt} \int_{\Omega(t)} Ad\Omega + \oint_{S(t)} \vec{C} \cdot \vec{n} dS = 0 \quad (3.2) $$

où:

- $A$ est une grandeur scalaire (resp. vectorielle)

- $\vec{C}$ est une grandeur vectorielle (resp. tensorielle)

- $S(t)$ est la surface de la cellule d'intégration $\Omega(t)$

- $\vec{n}$ est un vecteur unitaire perpendiculaire à $S(t)$ orienté vers l'extérieur

Il nous faut maintenant évaluer ces intégrales. Une solution serait de les calculer directement à partir des variables d'espace cartésiennes. Cependant, nous avons choisi de procéder d'abord à un passage en variables spatiales curvilignes, qui sont bien adaptées à la description d'un maillage structuré.

Ce type de maillage, couplé au changement de variables, permet une bonne discrétisation des frontières du domaine de calcul, et simplifie leur reconnaissance puisqu'une surface limite est alors un plan à coordonnée constante. D'autre part, nous réalisons ainsi le découplage de l'évaluation des grandeurs spatiales et aéro-dynamiques, ce qui nous offre la possibilité d'augmenter l'ordre de description de la géométrie indépendamment de celui retenu pour la solution. Si la forme du volume simulé est connue analytiquement, par des données de CAO ou des splines surfaciques, les paramètres géométriques peuvent même être calculés de manière exacte. Enfin, une gestion informatique simplifiée résulte de l'absence de table de connexion puisque la topologie d'un maillage structuré est intrinsèque et fixe.

Cependant, des géométries très complexes peuvent se révéler difficiles à mailler. De même, il est mal aisé de raffiner/déraffiner très localement, en agissant sur la densité des noeuds. D'autre part, de nouvelles grandeurs liées à la transformation géométrique, appelées les métriques, apparaissent. Elles induisent une légère complication des équations, mais entraînent surtout un surcoût de temps de calcul de part leur évaluation et une occupation mémoire accrue pour leur stockage.

Soient $(\xi^1, \xi^2, \xi^3, t)$ les variables curvilignes, $\sqrt{g} = \frac{D(x_1, x_2, x_3)}{D(\xi^1, \xi^2, \xi^3)}$ le jacobien de la transformation géométrique, et $\vec{a}^i = \vec{\nabla}\xi^i$ le vecteur de la base contravariante associé à la direction $\xi^i$. Les cellules d'intégration choisies sont des hexaèdres (figure 3.2), limités par des surfaces telles que $\xi^i$ soit constant sur $S_i$. Il en résulte que l'équation 3.2 se récrit:

$$ \int_{\omega} \frac{\partial}{\partial t} (\sqrt{g} A) \prod_{j=1}^{3} d\xi^j + \sum_{i=1}^{3} \left\{ \int_{S_i} \sqrt{g} \vec{C} \cdot \vec{a}^i \prod_{\substack{j=1 \\ j \neq i}}^{3} d\xi^j \right\} $$
<figure><img src="image_4.png" /><figcaption>Figure 3.2: Changement de variables d'espace</figcaption></figure>

$$ -\sum_{i=1}^{3} \left\{ \int_{S^i} \sqrt{g} \bar{C} \cdot \bar{a}^i \prod_{\substack{j=1 \\ j \neq i}}^{3} d\xi^j \right\} = 0 \qquad (3.3) $$

L'hypothèse principale de la formulation volumes finis est alors que les intégrands apparaissant dans l'équation 3.3 varient peu sur l'élément d'intégration. Elles peuvent donc être remplacées par une valeur moyenne par élément. Les approximations suivantes en découlent:

$$ \int_{\omega} \frac{\partial}{\partial t} (\sqrt{g} A) \prod_{j=1}^{3} d\xi^j \approx \frac{\partial}{\partial t} (\sqrt{g} A) \cdot \omega \qquad (3.4) $$

$$ \int_{S^i} \sqrt{g} \bar{C} \cdot \bar{a}^i \prod_{\substack{j=1 \\ j \neq i}}^{3} d\xi^j \approx \sqrt{g} \bar{C} \cdot \bar{a}^i \cdot S^i \qquad (3.5) $$

Si les cellules d'intégration sont choisies telles que $\xi^i = \xi^i \pm \frac{1}{2}, \omega = 1$ et $S^i = 1$. Si bien qu'une expression approchée de l'équation 3.3 est:

$$ \frac{\partial}{\partial t} (\sqrt{g} A) \Big|_{\xi^1, \xi^2, \xi^3} + \sum_{i=1}^{3} \left\{ \pm \sqrt{g} \bar{C} \cdot \bar{a}^i \Big|_{\xi^i \pm \frac{1}{2}} \right\} = 0 \qquad (3.6) $$

Par la suite, nous définirons le flux $F^i$ par:

$$ F^i = \sqrt{g} \bar{C} \cdot \bar{a}^i \qquad (3.7) $$
<figure><img src="image_5.png" /><figcaption>Figure 3.3: Stockage des variables conservatives aux noeuds</figcaption></figure>

Un certain nombre de contraintes sont attachées à cette approche. D'une part, le remplacement des intégrandes par une valeur moyenne n'autorise qu'un faible changement du champ aérodynamique à l'intérieur d'une maille. La dimension des structures décrites ne serait en particulier être inférieure à la taille des cellules d'intégration. Ainsi, pour un écoulement complexe, sujet à de fortes évolutions, des petites mailles seront nécessaires. Il est à noter que paradoxalement, cette technique, couplée à un calcul judicieux des flux, assure une bonne capture des chocs et des discontinuités ponctuelles. D'autre part, une certaine qualité de maillage est requise. Une évolution trop brutale de la taille des mailles, entrainant une discontinuité du jacobien, et des changements brusques de directions du maillage, produisant une discontinuité de la base contravariante, sont à proscire. Cependant, il n'existe pas de forte contrainte sur l'aspect de forme des cellules d'intégration.

La formulation volumes finis obtenue suppose que les variables conservatives sont connues aux noeuds du maillage, autour desquels sont bâtis les volumes élémentaires (figure 3.3). Une autre stratégie aurait consisté à stocker le champ aérodynamique au centre des cellules construites alors sur les noeuds (figure 3.4). L'approche retenue à l'avantage de positionner des noeuds sur la frontière du domaine de calcul, même si l'élaboration d'une maille dans ce cas est plus délicate.

### 3.2.2 Expression des flux

Rappelons que nous avons établi jusqu'à présent qu'une équation différentielle ayant pour forme:

$$ \frac{\partial A}{\partial t} + \nabla \cdot \vec{B} = 0 \qquad (3.8) $$
<figure><img src="image_6.png" /><figcaption>Figure 3.4: Stockage des variables conservatives au centre des mailles</figcaption></figure>

où:

- A est une grandeur scalaire (resp. vectorielle)

- $\vec{B}$ est une grandeur vectorielle (resp. tensorielle)

pouvait être approchée par:

$$ \left. \frac{\partial}{\partial t} (\sqrt{g} A) \right|_{\xi^1, \xi^2, \xi^3} + \sum_{i=1}^{3} \left\{ \pm F^i \left|_{\xi^i \pm \frac{1}{2}} \right\} = 0 \qquad (3.9) $$

avec:

$$ - F^i = \sqrt{g} \vec{C} \cdot \vec{a}^i $$

$$ - \vec{C} = \vec{B} - A \vec{W} \quad (\text{resp. } \vec{\vec{C}} = \vec{\vec{B}} - \vec{A} \otimes \vec{W}) $$

Il nous faut explicit er les flux $F^i$ dans le cas des équations de conservation de la masse 2.20, de la quantité de mouvement 2.21 et de l'énergie 2.22.

• Conservation de la masse

L'équation différentielle à considérer est:

$$ \frac{\partial \rho}{\partial t} + \vec{\nabla} \cdot [\rho \vec{V}] = 0 \qquad (3.10) $$

Suivant les notations utilisées, nous avons:
$$A = \rho \qquad (3.11)$$

$$\vec{B} = \rho \vec{V} \qquad (3.12)$$

$$\vec{C} = \rho(\vec{V} - \vec{W}) \qquad (3.13)$$

$$F^i = \sqrt{g \rho V^i} \qquad (3.14)$$

où, en remarquant que $-\vec{W} \cdot \vec{a}^i = \xi_i^i$ (annexe D.7),

$$V^i = \xi_{x_1}^i u_1 + \xi_{x_2}^i u_2 + \xi_{x_3}^i u_3 + \xi_i^i \qquad (3.15)$$

$V^i$ est la vitesse contravariante du fluide relative au maillage mobile dans la direction $\xi^i$.

## Conservation de la quantité de mouvement

L'équation différentielle à considérer est:

$$\frac{\partial \rho \vec{V}}{\partial t} + \vec{\nabla} \cdot [\rho \vec{V} \otimes \vec{V} + p \vec{I} - \frac{1}{R_e} \vec{\tau}_{l+t}] = \vec{0} \qquad (3.16)$$

Suivant les notations utilisées, nous avons:

$$\vec{A} = \rho \vec{V} \qquad (3.17)$$

$$\vec{B} = \rho \vec{V} \otimes \vec{V} + p \vec{I} - \frac{1}{R_e} \vec{\tau}_{l+t} \qquad (3.18)$$

$$\vec{C} = \rho \vec{V} \otimes (\vec{V} - \vec{W}) + p \vec{I} - \frac{1}{R_e} \vec{\tau}_{l+t} \qquad (3.19)$$

$$\vec{F}^i = \sqrt{g} \left\{ \rho \vec{V} \cdot V^i + p \cdot \vec{a}^i - \frac{1}{R_e} \vec{\tau}_{l+t} \cdot \vec{a}^i \right\} \qquad (3.20)$$

## Conservation de l'énergie

L'équation différentielle à considérer est:

$$\frac{\partial \rho E}{\partial t} + \vec{\nabla} \cdot \left[ (\rho E + p) \vec{V} - \frac{1}{R_e} \left\{ \gamma \left( \frac{\mu}{P_r} + \frac{\mu_t}{P_{rt}} \right) \vec{\nabla} e + \vec{\tau}_{l+t} \cdot \vec{V} \right\} \right] = 0 \qquad (3.21)$$

Suivant les notations utilisées, nous avons:

$$A = \rho E \qquad (3.22)$$

$$\vec{B} = (\rho E + p) \vec{V} - \frac{1}{R_e} \left\{ \gamma \left( \frac{\mu}{P_r} + \frac{\mu_t}{P_{rt}} \right) \vec{\nabla} e + \vec{\tau}_{l+t} \cdot \vec{V} \right\} \qquad (3.23)$$
$$
\vec{C} = (\rho E + p) \cdot (\vec{V} - \vec{W}) + p\vec{W} - \frac{1}{R_e} \left\{ \gamma \left( \frac{\mu}{P_r} + \frac{\mu_t}{P_{rt}} \right) \vec{\nabla} e + \vec{\tau}_{l+t} \cdot \vec{V} \right\} \quad (3.24)
$$

$$
F^i = \sqrt{g} \left\{ (\rho E + p) V^i - p \xi_t^i - \frac{1}{R_e} \left\{ \gamma \left( \frac{\mu}{P_r} + \frac{\mu_t}{P_{rt}} \right) \vec{\nabla} e + \vec{\tau}_{l+t} \cdot \vec{V} \right\} \cdot \vec{a}^i \right\} \quad (3.25)
$$

3.2.3 Equations discrétisées

Ainsi, le système d'équations différentielles à résoudre initial se réécrit sous forme
discrète en espace de la manière suivante.

$$
\left. \frac{\partial}{\partial t} (\sqrt{g} q) \right|_{\xi^1, \xi^2, \xi^3} + \sum_{i=1}^{3} \left[ F^i(\xi^i + \frac{1}{2}) - F^i(\xi^i - \frac{1}{2}) \right] = 0 \quad (3.26)
$$

avec:

$$
q = \begin{cases} \rho \\ \rho \vec{V} \\ \rho E \end{cases} \tag{3.27}
$$

$$
F^i = F_c^i - \frac{1}{R_c} F_v^i \qquad (3.28)
$$

$$
= \sqrt{g} \left| \begin{array}{l} \rho V^i \\ \rho \vec{V} V^i + p \vec{a}^i \\ (\rho E + p) V^i - p \xi_t^i \end{array} - \frac{1}{R_e} \sqrt{g} \left| \begin{array}{l} 0 \\ \vec{\tau}_{l+t} \\ \gamma \left( \frac{\mu}{P_r} + \frac{\mu_t}{P_{rt}} \right) \vec{\nabla} e + \vec{\tau}_{l+t} \cdot \vec{V} \end{array} \right| \right| \cdot \vec{a}^i \quad (3.29)
$$

$$
V^i = \xi_{x_1}^i u_1 + \xi_{x_2}^i u_2 + \xi_{x_3}^i u_3 + \xi_t^i \quad (3.30)
$$

$F_c^i$ (resp. $F_v^i$) représente le flux convectif (resp. diffusif) dans la direction $\xi^i$.

Il est possible d'établir qu'utiliser une formulation volumes finis en espace car-
tésien avec changement de variables d'espace en curviligne, tel que nous l'avons fait
jusqu'à présent, peut être assimilé à appliquer une approche différences finies en
repère curviligne (annexe C.6). L'apport principal de cette constatation réside dans
le fait que l'équation 3.26 peut être interprétée comme découlant d'une équation
différentielle qui serait:

$$
\frac{\partial}{\partial t} (\sqrt{g} q) + \sum_{i=1}^{3} \frac{\partial F^i}{\partial \xi^i} = 0 \quad (3.31)
$$

Cette remarque sera mise à profit lors du traitement des conditions aux limites.
## 3.3 Discrétisation temporelle

Il résulte de la formulation volumes finis 3.26 que l'évolution numérique du champ aérodynamique est décrite sous une forme semi-discrète en espace par l'équation:

$$ \frac{\partial}{\partial t}(\sqrt{g}q)\Big|_{\xi^1, \xi^2, \xi^3} = -R|_{\xi^1, \xi^2, \xi^3, t} \quad (3.32) $$

où R, dit le résidu, est défini par:

$$ R = \sum_{i=1}^{3} \left[ F^i(\xi^i + \frac{1}{2}) - F^i(\xi^i - \frac{1}{2}) \right] \quad (3.33) $$

L'intégration en temps de l'équation 3.32 nous permet alors d'obtenir la solution $q$ à l'instant $(n+1)$, $q^{(n+1)}$, à partir de sa valeur à un ou plusieurs instants antérieurs, $q^{(n)}$, $q^{(n-1)}$, $q^{(n-2)}$, ...

Il existe de nombreux schémas numériques réalisant cette opération. Un certain nombre d'entre eux peuvent être regroupés sous une formulation générale, dite *Delta de Pade* [7], qui s'exprime comme suit:

$$ \Delta (\sqrt{g}q)^n + \frac{\theta \Delta \tau}{1 + \beta} \Delta R^n = - \frac{\Delta \tau}{1 + \beta} R^n + \frac{\beta}{1 + \beta} \Delta (\sqrt{g}q)^{(n-1)} + O\left[ \left( \theta - \frac{1}{2} - \beta \right) \Delta \tau^2 + \Delta \tau^3 \right] \quad (3.34) $$

où:

* $\Delta \tau$ est le pas en temps séparant les instants $(n)$ et $(n+1)$:

$$ t^{(n+1)} = t^{(n)} + \Delta \tau \quad (3.35) $$

* $\Delta A^n$ est l'incrément de $A$ durant l'intervalle de temps $\Delta \tau$:

$$ A^{(n+1)} = A^{(n)} + \Delta A^n \quad (3.36) $$

* $\theta$ et $\beta$ sont deux paramètres fixant le caractère explicite ou implicite du schéma, et sa précision. Ainsi:

- $\theta = 0, \beta = 0$: Schéma d'Euler explicite du premier ordre.

- $\theta = \frac{1}{2}, \beta = 0$: Schéma trapezoidal implicite du second ordre.

- $\theta = 1, \beta = 0$: Schéma d'Euler implicite du premier ordre.

- $\theta = 1, \beta = \frac{1}{2}$: Schéma décentré en trois points implicite du second ordre.
### 3.3.1 Remarques sur les schémas implicites en instationnaire

Lorsque $\theta \neq 0$, l'équation 3.34 implique que la solution en un point à l'instant $(n+1)$ dépend du champ aérodynamique au même instant en d'autres noeuds, à cause du terme $\Delta R^n$, qui crée une dépendance spatiale. Il en résulte un système couplé, d'autant plus complexe que la molécule de calcul utilisée est étendue.

Cette classe de schémas autorise des pas en temps importants, particulièrement lorsque peu d'hypothèses simplificatrices sur l'évaluation du terme $\Delta R^n$ sont faites [9]. Elle est donc bien adaptée à la recherche d'un état stationnaire, pour lequel $t$ tend vers l'infini, ou la description de l'évolution d'un processus lent, tel que ce soit la vitesse de propagation de l'information physique qui impose le plus grand pas en temps admissible, et non pas des considérations de stabilité numérique du schéma.

Cependant, un certain nombre de défauts s'attachent à cette méthode. Tout d'abord, la qualité de la simulation d'un écoulement instationnaire est fortement tributaire des déphasages parasites qui sont introduits artificiellement dans la solution par le schéma en temps. Aussi, afin de les minimiser, aucun découplage en espace, tel que la technique des directions alternées [8], ne doit être utilisé. Il en résulte un système linéaire imposant. Par expérience [9], nous savons que l'occupation mémoire qui en découle peut devenir irréaliste par rapport à la capacité des ordinateurs actuels, atteignant 1 GigaOctets pour un maillage de 100 000 noeuds. D'autre part, l'évaluation du terme $\Delta R^n = R^{(n+1)} - R^{(n)}$ est primordiale. Ce dernier est le plus souvent modélisé à partir d'une linéarisation au premier ordre de $R^{(n+1)}$ autour de l'état $q^{(n)}$ [39]. Mais l'expression mathématiquement exacte du jacobien de $R$ est tellement complexe que seule une forme plus ou moins simplifiée en est traditionnellement utilisée. Une autre technique recourt à des sous-pas en temps numériques, entre les deux instants physiques ($n$) et ($n+1$), permettant d'approcher itérativement la valeur numérique de $\Delta R^n$, avant d'évaluer l'incrément de la solution physique $\Delta q^n$ [29]. Cependant, il semble aujourd'hui que ces méthodes soient limitées en précision à des domaines fréquentiels n'excédant pas cinq fois la fréquence réduite du phénomène étudié. Enfin, une approche implicite se révèle particulièrement lourde à gérer sur le plan informatique lors du développement et de l'extension d'un code de calcul, d'une part pour l'implémentation, d'autre part pour la validation de nouvelles techniques de discrétisation spatiale, ou de traitement des conditions aux limites par exemple.

### 3.3.2 Schéma d'Euler explicite

Suite à l'analyse des schémas implicites précédemment présentée, un schéma en temps explicite a été retenu pour cette étude. Le plus simple de tous, inclus dans la formulation de Pade, est sans doute le schéma d'Euler explicite. Il exprime que:

$$ \Delta (\sqrt{gq})^n = -\Delta\tau R^n + O(\Delta\tau^2) \quad (3.37) $$
Il est donc précis au premier ordre, et ne requiert que la connaissance de la solution à l'instant (n) pour évaluer Rⁿ, puis q(n+1).

Si sa stabilité peut être théoriquement établie sous une condition de CFL classique pour un schéma bâti sur une séparation de flux au premier ordre [39, p. 40], il n'en est rien pour une précision en espace plus élevée. Cependant, les schémas TVD, qui assurent partiellement la monotonie de la solution, se révèlent stables en pratique sous la même condition, qui déclare que:

$$CFL = (U + a) \frac{\Delta t}{\Delta x} \le 1 \qquad (3.38)$$

où:

* U est une échelle de vitesse caractéristique de l'écoulement au noeud considéré, par exemple le module de la vitesse du fluide $|| \vec{V} ||$.

* a est la vitesse du son.

* Δx est une échelle de longueur caractéristique de la géométrie au noeud considéré, par exemple la distance au point voisin le plus proche.

Le couple $U + a$ représente la vitesse de l'information physique se propageant le plus rapidement. Si bien que la condition 3.38 établit que cette onde ne peut pas parcourir en une itération une distance supérieure à celle séparant deux noeuds voisins. Ceci est le propre d'un schéma explicite: en un pas de temps élémentaire, un point n'influence que les noeuds qui lui sont reliés, contrairement à un schéma implicite où tous les points sont inter-dépendants.

Cette technique présente l'avantage d'une extrême simplicité, et ne requière que peu d'espace mémoire, de l'ordre de 140 MegaOctets pour 100 000 points dans la version actuelle du code. Cependant, la recherche d'un état stationnaire, ou la description d'une période d'un phénomène lent, peuvent nécessiter un temps de calcul important, puisque plus la vitesse est grande, et la taille des mailles petite, plus faible est Δt. C'est pourquoi des méthodes plus évoluées, telles qu'un schéma de Runge-Kutta [39] en instationnaire ou un lissage du residu en stationnaire[41], sont actuellement à l'étude. Elles présentent pour principales caractéristiques une précision en temps accrue, ainsi qu'une mise en oeuvre aisée à partir du schéma d'Euler explicite, qui offre donc une bonne base de travail.

## 3.4 Discrétisation spatiale des flux convectifs

Il nous faut maintenant évaluer le flux convectif $F_c^i$ apparu lors de la formulation volumes finis (équation 3.26). Il est défini par:
$$ F_c^i = \sqrt{g} \cdot \begin{vmatrix} \rho V^i & & & \\ \rho u_1 V^i & + p \xi_{x1}^i & & \\ \rho u_2 V^i & + p \xi_{x2}^i & & \\ \rho u_3 V^i & + p \xi_{x3}^i & & \\ (\rho E + p) V^i & - p \xi_i^i & & \end{vmatrix} \quad (3.39) $$

Et il doit être calculer à l'interface entre deux cellules d'intégration en lesquelles est connu le champ aérodynamique. La technique retenue pour ce faire doit assurer une bonne capture des discontinuités, potentiellement présentes dans un écoulement transsonique, et la stabilité du schéma numérique, en ayant pas ou peu recours à des paramètres ajustables en fonction de la solution.

### 3.4.1 Schémas décentrés amont

Les schémas décentrés amont répondent à ces deux critères. L'idée présidant à leur élaboration est la suivante. Il apparaît judicieux de s'inspirer de la réalité physique pour évaluer le flux convectif, en remarquant qu'il résulte de la combinaison d'informations transportées par des ondes se propageant dans différents sens. La discrétisation de ces dernières est alors réalisée en fonction de leur sens de déplacement. Cette technique nécessite d'importants développements mathématiques et physiques, mais elle offre en contrepartie une nouvelle approche guidée par la réalité très satisfaisante [27, pp. 408 – 492].

Plusieurs avantages sont attachés à cette classe de schémas. Tout d'abord, ils sont naturellement dissipatifs, et ne requièrent donc pas l'adjonction d'une dissipation artificielle supplémentaire. La stabilité numérique est ainsi assurée, sans avoir à régler de coefficients arbitraires. D'autre part, de par la philosophie suivie pour leur construction, ils ont la faculté intrinsèque de gérer des ondes telles que des chocs ou des surfaces de constact, stationnaires mais aussi et surtout instationnaires. Ainsi sont minimisés tous déphasages parasites dans la solution.

Cependant, ces méthodes plus sophistiquées entraînent un surcoût de temps de calcul et une implémentation plus complexe car, étant plus "intelligentes", elles nécessitent d'avantages de tests, ne facilitant pas en particulier la vectorisation des codes de calcul.

Notre choix [9] s'est porté, parmi les solutions existantes dans la littérature [27, pp. 438 – 472], sur la séparation de flux de Van Leer. Cette technique a initiallement été développée sur un maillage cartésien tri-dimensionnel [36]. Puis sa construction a été généralisée par Anderson, Thomas et Van Leer [3], à des applications utilisant un système de coordonnées curvilignes. Enfin, Parpia [49] en a proposé une formulation incluant le déplacement des noeuds du maillage. Un résumé succinct de ce dernier développement est présenté par la suite.

Le flux convectif $F_c^i$ est supposé résulter de la somme de deux contributions $F^{i+}$ et $F^{i-}$, telles que les matrices jacobiennes associées, $\partial F^{i+} / \partial(\sqrt{\bar{q}q})$ et $\partial F^{i-} / \partial(\sqrt{\bar{q}q})$,
aient respectivement, aucune valeur propre négative, ou positive. Chaque composante du flux est alors calculée en utilisant des variables décentrées amont ou aval, en fonction du signe des valeurs propres de la matrice jacobienne associée. Une infinité d'expressions de $F^{i+}$ et $F^{i-}$ répondent à ce critère. Celles proposées par Van Leer ont été déterminées sous les contraintes supplémentaires suivantes.

* $F^{i\pm}$ doit être continu, tel que $F^{i+} = F_c^i$ pour un nombre de Mach supérieur à 1, et $F^{i-} = F_c^i$ pour un nombre de Mach inférieur à -1.

* Les composantes de $F^{i\pm}$ doivent reproduire les symétries de $F_c^i$ par rapport au nombre de Mach, telles que $F^{i+}(M) = \pm F^{i-}(-M)$ si $F_c^i(M) = \pm F_c^i(-M)$.

* $\partial F^{i\pm}/\partial(\sqrt{\overline{gq}})$ doit être continue.

* $\partial F^{i\pm}/\partial(\sqrt{\overline{gq}})$ doit avoir une valeur propre qui s'annule pour des nombres de Mach subsoniques.

* $F^{i\pm}$, tout comme $F_c^i$, doit être un polynôme en nombre de Mach, du plus petit degré possible.

Il en résulte que le flux convectif dans la direction $\xi^i$ est décomposé en fonction du nombre de Mach contravariant $M^i$ de la manière suivante:

$$F_c^i = F^{i+} + F^{i-} \qquad (3.40)$$

$$M^i = \frac{V^i}{a \|\vec{a}^i\|} \qquad (3.41)$$

avec:

* Pour un écoulement supersonique ascendant: $M^i \le -1$

$$F^{i+} = 0 \qquad (3.42)$$

$$F^{i-} = F_c^i \qquad (3.43)$$

* Pour un écoulement subsonique: $-1 < M^i < 1$

$$F^{i\pm} = \sqrt{g} \|\vec{a}^i\| f_{\text{masse}}^{\pm} \cdot \begin{cases} 1 & \gamma \|\vec{a}^i\| \le M^i \pm 2 \quad \xi_{x_1}^i + u_1 \\ \frac{a}{\gamma \|\vec{a}^i\|} (-M^i \pm 2) \xi_{x_2}^i + u_2 & \gamma \|\vec{a}^i\| \le M^i \pm 2 \quad \xi_{x_3}^i + u_3 \\ \frac{a}{\gamma \|\vec{a}^i\|} (-M^i \pm 2) \xi_{x_3}^i + u_3 & \gamma \|\vec{a}^i\| \le M^i \pm 2 \quad \xi_{\text{energic}}^i + u_3 \end{cases} \qquad (3.44)$$
<figure><img src="image_7.png" /><figcaption>Figure 3.5: Approche MUSCL</figcaption></figure>

où:

$$f_{\text{masse}}^{\pm} = \pm \frac{1}{4} \rho a (M^{i} \pm 1)^{2} \qquad (3.45)$$

$$f_{\text{energie}}^{\pm} = \frac{a^2}{\gamma^2 - 1} \left[ -(\gamma - 1)M^{i2} \pm 2(\gamma - 1)M^{i} + 2 \right] + \frac{1}{2} (u_1^2 + u_2^2 + u_3^2) - \frac{a}{\gamma} \frac{\xi_t^i}{\|\tilde{a}^i\|} (-M^i \pm 2) \quad (3.46)$$

* Pour un écoulement supersonique descendant: $M^i \ge 1$

$$F^{i+} = F_c^i \qquad (3.47)$$

$$F^{i-} = 0 \qquad (3.48)$$

D'autre part, une technique MUSCL (Monotone Upstream-centered Scheme for Conservation Laws) est utilisée pour évaluer ces flux. Elle a été introduite pour la première fois par Van Leer [35] pour obtenir une précision spatiale d'ordre élevé, associée à un schéma TVD (Total Variation Diminishing). Dans cette approche, les variables conservatives sont interpolées à l'interface des mailles, puis les flux sont calculés; alors que par une méthode plus classique, les flux sont évalués au centre des cellules, puis interpolés sur leur frontière. Ainsi, pour calculer $F^{i+}$ (resp. $F^{i-}$), une variable décentrée amont $Q^-$ (resp. décentrée aval $Q^+$) est utilisée (figure 3.5). De plus, afin d'éviter toute oscillation à proximité des discontinuités de la solution, telles que les chocs, et pour assurer la monotonie du schéma, l'amplitude des gradients
locaux est restreinte entre deux bornes. Ceci est réalisé grâce à une fonction non linéaire $\Phi$, appelée un *limiteur*. En outre, les métriques, telles que $\sqrt{g}$ et $\|\vec{a}^i\|$, sont évaluées de manière exacte sur les surfaces d'intégration, et non pas interpolées à partir des valeurs connues au centre des mailles.

Finalement, nous obtenons pour $F_c^i(\xi^i - \frac{1}{2})$ par exemple:

$$F_c^i(\xi^i - \frac{1}{2}) = F^{i+}(Q_{\xi^i-\frac{1}{2}}^{-}) + F^{i-}(Q_{\xi^i+\frac{1}{2}}^{+}) \quad (3.49)$$

avec:

$$ 
\begin{aligned}
Q_{\xi^i-\frac{1}{2}}^{-} &= Q(q_{\xi^i}, q_{\xi^i-1}, q_{\xi^i-2}) \\
Q_{\xi^i-\frac{1}{2}}^{+} &= Q(q_{\xi^i-1}, q_{\xi^i}, q_{\xi^i+1})
\end{aligned}
\quad (3.50) 
$$

$$ 
\begin{aligned}
Q(q_0, q_1, q_2) &= q_1 + \frac{1}{4}\Phi[(1 + \kappa\Phi)(q_0 - q_1) + (1 - \kappa\Phi)(q_1 - q_2)] \quad (3.51) \\
&= \frac{1}{4}\Phi[1 + \kappa\Phi]q_0 + \left[1 - \frac{1}{2}\kappa\Phi^2\right]q_1 - \frac{1}{4}\Phi[1 - \kappa\Phi]q_2 \quad (3.52)
\end{aligned}
$$

où:

$$ \Phi = \Phi(q_0, q_1, q_2) \quad (3.53) $$

La précision spatiale est contrôlée par la fonction $\Phi$, comprise entre 0 pour un schéma au premier ordre et 1 pour un schéma d'ordre supérieur, et le paramètre $\kappa$ tel que nous ayons pour $\Phi = 1$:

* $\kappa = -1$: Schéma totalement décentré du second ordre

* $\kappa = 0$: Schéma de Fromm

* $\kappa = \frac{1}{3}$: Schéma partiellement décentré du troisième ordre

* $\kappa = 1$: Schéma centré du second ordre

Une valeur différente du limiteur $\Phi$ est calculée pour chaque composante de $Q^\pm$ en fonction de la variable conservative correspondante. Il en résulte que l'ordre d'interpolation en un point peut changer selon la sensibilité d'une grandeur particulière à une discontinuité spécifique, telle un choc, une surface de constact, une onde de détente, ...

Il est à noter que les développements présentés supposent une répartition de points équidistants, en particulier lors de l'évaluation de la précision spatiale (équation 3.51). Des modifications seraient à introduire pour prendre en compte une répartition irrégulière des noeuds, ainsi que des déformations géométriques, telles que des courbures [38]. Cependant, des résultats corrects sont obtenus si le maillage présente une évolution régulière. Nous remarquerons que, dans l'approche développée, le
<figure><img src="image_8.png" /><figcaption>Figure 3.6: Tube à choc au premier ordre. Niveau de densité</figcaption></figure>

découplage de l'évaluation des termes géométriques par rapport au champ aéro-dynamique rend très facile et peu coûteux en temps de calculer les métriques à un ordre élevé, et permet ainsi d'améliorer la description de géométries complexes, telles qu'un bord d'attaque rond par exemple.

Il existe de nombreux autres schémas décentrés amont, tels que la séparation de flux de Steger et Warming [62], ou encore les solveurs de Riemann approchés proposés par Roe [56] ou par Osher [46], ... Notre choix s'est porté sur celui de Van Leer car d'après la littérature, il produit de bons résultats, excepté il est vrai à très bas nombre de Mach, il est facile à mettre en oeuvre et peu coûteux d'utilisation. Ceci est dû à l'expression algébrique simple des contributions du flux, qui permet une analyse théorique du schéma et le rend accessible à quelques modifications. Il est d'autre part un des plus simples à discrétiser de manière implicite en temps. En d'autres termes, nous l'avons choisi car il conduit à des solutions honorables tout en offrant une bonne base de travail, raisonnablement ouverte.

A titre d'illustration des capacités de la séparation de flux de Van Leer, les résultats obtenus dans le cas d'un écoulement élémentaire en tube à choc, incluant des discontinuités instationnaires, sont présentés par la suite [4]. Ils sont comparés à la solution analytique (annexe E.1), ainsi qu'à celle obtenue par une méthode de Roe implémentée en maillage non structuré [25]. Les profils de masse volumique et d'énergie totale, au premier et au second ordre en espace après un temps significatif, sont visibles figures 3.6, 3.7, 3.8 et 3.9. Nous notons tout d'abord la faible différence
<figure><img src="image_9.png" /><figcaption>Figure 3.7: Tube à choc au second ordre. Niveau de densité</figcaption></figure>

<figure><img src="image_10.png" /><figcaption>Figure 3.8: Tube à choc au premier ordre. Niveau d'énergie totale</figcaption></figure>
<figure><img src="image_11.png" /><figcaption>Figure 3.9: Tube à choc au deuxième ordre. Niveau d'énergie totale</figcaption></figure>

entre les deux méthodes numériques, quelle que soit la précision spatiale retenue. D'autre part, une nette amélioration de la solution, en particulier au niveau de la surface de contact, est obtenue lorsque l'ordre en espace augmente. Ainsi les fronts de l'onde de détente ne sont alors que faiblement arrondis, pendant que l'onde de pression est étalée sur 2 ou 3 cellules, et la surface de contact sur 7 ou 8 mailles.

Le deuxième cas test présenté, toujours bâti sur un tube à choc mais mettant en oeuvre cette fois des interactions d'ondes instationnaires, révèle le même comportement d'une part vis-à-vis de la solution analytique (annexe E.2), d'autre part vis-à-vis de la méthode de Roe. Figure 3.10 est tracé le niveau de densité après que l'onde de pression centrale se soit réfléchie deux fois à l'extrémité du tube et une fois sur la surface de contact, en générant en particulier un choc transmis. Lors de ce calcul au second ordre, la position et l'amplitude des différentes ondes sont correctement évaluées. Tous les chocs restent étalés sur 3-4 cellules, alors que les surfaces de contact le sont sur 7-8 mailles. La dissipation des discontinuités n'augmente pas au cours du temps, ni lors des interactions puisque même le dernier choc central, de faible amplitude et de troisième génération, est correctement capturé.

Ces deux exemples établissent partiellement le bien fondé des choix faits en matière de discrétisation spatiale, puisqu'ils montrent la capacité du schéma retenu à gérer des ondes instationnaires, telles des chocs, même lorsque de nombreuses interactions se produident. Confronté à la méthode de Roe, seule une gestion des surfaces de contact légèrement meilleure peut être créditée au mérite de cette dernière.
<figure><img src="image_12.png" /><figcaption>Figure 3.10: Intéractions d'ondes en tube à choc. Niveau de densité</figcaption></figure>

### 3.4.2 Discussion autour des limiteurs

#### Définition et rôle du limiteur

Le limiteur $\Phi$, intervenant dans l'interpolation du champ aérodynamique à l'interface entre deux cellules (équation 3.51), est une fonction non linéaire détectant automatiquement les discontinuités de la solution. Il assure le caractère TVD du schéma en dégradant la précision spatiale au premier ordre dans ces zones, évitant ainsi toute oscillation parasite en espace [27, 39, pp. 552-556].

D'un point de vue mathématique, $\Phi$ est une fonction bornée entre 0 et 1, bâtie sur les gradients locaux de chaque côté de l'interface considérée, dans le sens du décentrement. Ainsi, pour un maillage régulier:

$$ \Phi(q_0, q_1, q_2) = \Phi(q_0 - q_1, q_1 - q_2) \quad (3.54) $$

Il est possible de proposer une infinité d'expressions pour $\Phi$. Cependant, un certain nombre de critères, en sus de ceux inclus dans sa définition, se révèlent importants lors du choix d'un limiteur. En particulier, il semble souhaitable qu'il soit continu et dérivable. Et ce afin d'éviter tout changement de précision brutal en un noeud lors du déplacement d'une discontinuité, surtout si elle oscille dans le temps. Ceci a une influence en instationnaire, mais aussi lors de la recherche d'un état stationnaire, où un "bon" limiteur facilitera le positionnement stable des discontinuités.
## Critiques de limiteurs classiques

Il existe de nombreux limiteurs dans la littérature. En voici deux parmi les plus populaires.

### • Minmod

Hänel [26] a montré que pour la résolution des équations de Navier-Stokes, de meilleurs résultats sont obtenus avec le limiteur minmod ainsi défini:

$$ \Phi(q_0, q_1, q_2) = \minmod(1, \omega \frac{q_1 - q_2}{q_0 - q_1}) \quad (3.55) $$

avec

$$ \minmod(a,wb) = sgn(a) \cdot \max \{0, \min \{|a|,wb \,sgn(a)\}\} \quad (3.56) $$

Le paramètre $\omega$, tel que $1 \le \omega \le \frac{3-\kappa}{1-\kappa}$, influence la sensibilité du limiteur. Plus sa valeur est grande, plus important est l'écart entre les états droit et gauche de chaque coté de l'interface, et plus franches sont les discontinuités.

Cependant, il n'est pas continûment dérivable, et Anderson [2] a établi qu'il en résultait un blocage de la convergence vers un état stationnaire. D'autre part, il conduit a une médiocre résolution des surfaces de contact [27, pp. 542].

### • Mulder

Ce limiteur est quant à lui, le fruit d'une modification du limiteur de Van Albada, qui est continûment dérivable, proposée par Mulder [45]:

$$ \Phi(q_0, q_1, q_2) = \frac{2(q_0 - q_1)(q_1 - q_2) + \varepsilon}{(q_0 - q_1)^2 + (q_1 - q_2)^2 + \varepsilon} \quad (3.57) $$

où $\varepsilon$ est une constante très petite ($\varepsilon \approx 10^{-6}$) introduite pour éviter une division par zéro dans une zone à gradient nul. Son choix n'est cependant pas sans conséquence, puisqu'une trop forte valeur de $\varepsilon$ conduit à une limitation peu efficace. Il peut donc être nécessaire d'ajuster ce paramètre en fonction de l'écoulement.

D'après les travaux d'Anderson [2], le limiteur de Mulder se comporte mieux que le limiteur minmod.

Cependant, il n'est pas sans défaut. Tout d'abord, il présente une forte sensibilité à la valeur de $\varepsilon$ pour un écoulement uniforme, ou faiblement variable. En effet, si un des deux gradients $(q_0 - q_1)$, ou $(q_1 - q_2)$, est faible devant cette constante, c'est le rapport du second à $\varepsilon$ qui fixe la valeur du limiteur. Ainsi, dans le cas extrême où $(q_1 - q_2)$ est nul,
<figure><img src="image_13.png" /><figcaption>Figure 3.11: Tube à choc. Valeur initiale du limiteur de Mulder</figcaption></figure>

$$ \Phi(q_0, q_1, q_2) = \frac{\epsilon}{(q_0 - q_1)^2 + \epsilon} \quad (3.58) $$

$$ \rightarrow 0 \quad \text{si} \quad (q_0 - q_1)^2 \gg \epsilon \quad (3.59) $$

$$ \rightarrow 1 \quad \text{si} \quad (q_0 - q_1)^2 \ll \epsilon \quad (3.60) $$

Ceci est illustré par les valeurs du limiteur présentées figures 3.11 à 3.13. Elles sont obtenues dans le cas d'un écoulement en tube à choc, initialisé à partir de la solution analytique à un instant où le faisceau d'ondes est pleinement développé. Y sont dessinés le niveau de masse volumique, ainsi que la valeur du limiteur associé, en fonction de l'abscisse du tube. Les discontinuités de la solution initiale sont parfaitement détectées figure 3.11, le limiteur y étant à 0. Figures 3.12 et 3.13 sont tracées les mêmes grandeurs après que la solution est évoluée en temps, pour deux valeurs de $\epsilon$, respectivement $10^{-6}$ et $10^{-12}$. Comme prédit par la relation 3.59, plus $\epsilon$ est petit, plus le limiteur est faible. Ainsi la précision globale du schéma en espace est-elle moins élevée dans le second cas. Nous noterons d'autre part une chute de $\Phi$ entre $X=2.5$ et $3$ pour $\epsilon = 10^{-12}$, alors que le niveau de masse volumique y est semble-t-il constant. Elle résulte du fait que pour certains points $(q_1 - q_2)$ est exactement nul, et $(q_0 - q_1)^2$ est supérieur à $10^{-12}$, $10^{-10}$ par exemple. Le passage à l'ordre un dans une telle région est sans importance sur la solution globale car la fonction est constante, mais il traduit quant même un comportement jugé peu satisfaisant.
<figure><img src="image_14.png" /><figcaption>Figure 3.12: Tube à choc. Limiteur de Mulder. $\epsilon = 10^{-6}$</figcaption></figure>

<figure><img src="image_15.png" /><figcaption>Figure 3.13: Tube à choc. Limiteur de Mulder. $\epsilon = 10^{-12}$</figcaption></figure>
<figure><img src="image_16.png" /><figcaption>Figure 3.14: Limiteur en cosinus</figcaption></figure>

D'autre part, nous noterons un effet de "mémoire" du limiteur qui maintient
dans les deux cas une précision faible autour du point X=3.7, correspondant
à la position initiale du choc, alors que ce dernier s'est déplacé au delà de
X=4.0. Ce phénomène, lié peut-être à la remarque précédente, n'a pas encore
été expliqué.

Bilan

Il ressort de ces deux exemples que le choix d'un limiteur n'est pas chose aisée, et semble-t-il fortement influencé par le type d'écoulements simulés, voir même par le cas spécifique étudié. Or c'est cette fonction qui conditionne le traitement des discontinuités, et la qualité de la solution obtenue. Si elle est trop sensible, de nombreuses régions seront sans raison calculées au premier ordre et une dissipation numérique excessive en résultera; si au contraire elle est trop peu sensible, des oscillations parasites apparaîtront. Le limiteur universel idéal ne semble pas encore avoir été découvert. Une attitude pragmatique, consistant pour chaque configuration étudiée à vérifier le comportement in situ de l'expression de Φ retenue, est conseillée.

Pour notre part, nous nous sommes essayé à la construction de nouveaux li-
miteurs, qui ne présentent pas les défauts de celui de Mulder pour un écoulement
uniforme. Le premier est bâti sur le cosinus de l'angle fait par les gradients (q₁ - q₂)
et (q₀ - q₁) (figure 3.14), le second sur l'écart entre l'état q₁ et celui prédit par une
interpolation linéaire entre q₀ et q₂ (figure 3.15). Leur principal défaut est qu'ils sont
sensibles à l'ordre de grandeur de la solution, ce qui nécessite l'emploi de constan-
tes d'adimensionnement à ajuster au cas par cas. Cependant, bien qu'ils ne soient
pas exploitables en l'état, ils semblent prometteurs, mais demandent d'avantage de
<figure><img src="image_17.png" /><figcaption>Figure 3.15: Limiteur en tangente hyperbolique</figcaption></figure>

développements.

### 3.4.3 Problèmes liés à la séparation de flux

Si les schémas décentrés amont apparaissent comme très élégants sur le plan physique et théorique, ils ne se révèlent pas irréprochables en pratique. Ainsi, différents défauts ont été constatés; ils sont exposés par la suite. Certains sont généraux, alors que d'autres sont plus spécifiques à la séparation de flux de Van Leer.

#### Sensibilité à l'orientation du maillage

Les schémas décentrés amont ont initialement été développés pour décrire des écoulements mono-dimensionnels de fluides parfaits. Dans ces conditions, la nature, la direction et le sens de propagation des ondes qui servent de support à la discrétisation sont non ambigus. Puis ces techniques ont été étendues, jusqu'à permettre le traitement des effets convectifs multi-dimensionnels pour un fluide visqueux. Or une partie sensible de la rigueur théorique initiale a fait les frais de cette généralisation. En particulier, une influence non négligeable de l'orientation du maillage par rapport aux discontinuités peut être mise en évidence dans quelques cas simples.

L'un d'entre eux met en oeuvre un écoulement bi-dimensionnel, résultant de la rupture d'un diaphragme cylindrique entre une chambre basse pression et une chambre haute pression (annexe E.3). Le même faisceau d'ondes de choc, de contact et de détente que celui rencontré en tube à choc se développe, si ce n'est que les différentes discontinuités sont maintenant circulaires. Les niveaux de masse volumique obtenus après un temps suffisant sur un maillage radial (resp. cartésien) sont présentés figure 3.16 (resp. figure 3.17). Alors que des plateaux réguliers séparent les ondes dans
<figure><img src="image_18.png" /><figcaption>Figure 3.16: Choc cylindrique. Niveau de masse volumique en maillage radial</figcaption></figure>

<figure><img src="image_19.png" /><figcaption>Figure 3.17: Choc cylindrique. Niveau de masse volumique en maillage cartésien</figcaption></figure>
le premier cas où les mailles sont alignées avec l'écoulement, de petites vagulettes indésirables sont visibles dans le second. Elles résultent de transferts angulaires de masse, qui semblent générés lors du déplacement des discontinuités à travers des cellules d'intégration mal orientées. Ces défauts ne sont pas propres à la séparation de flux de Van Leer, puisque le solveur approché de Riemann de Roe en maillage non structuré présente le même comportement, comme le montre les petits champignons visibles sur la figure 3.18.

Deux solutions existent à ce problème. La première consiste à adapter le maillage pour assurer l'alignement des discontinuités et des mailles. Elle n'est pas toujours applicable, en particulier lorsque d'autres contraintes géométriques, telle que la forme du canal inter-aubes, conditionnent fortement la génération du maillage. Le second remède quant à lui met en oeuvre une modification de la séparation de flux, qui n'est plus réalisée dans une direction normale à la surface d'une cellule d'intégration, mais perpendiculairement aux discontinuités [37, 57, 14]. Cette approche relativement récente n'a pas été testée dans le travail présent.

Mauvaise transition entre les zones subsoniques et supersoniques

Lors du calcul de configurations transsoniques stationnaires en tuyère symétrique (paragraphe 6.6), nous avons constaté la présence de sur-estimations ou de sous-estimations de certaines grandeurs telle que la pression totale, au niveau du choc, et ce même sur des résultats au premier ordre en espace réputés monotones (figure 3.19). Ce défaut semble être lié à la génération d'erreurs locales lorsque le régime d'écoulement passe de supersonique à subsonique. En effet dans cette région, le décentrement change et le raccord entre les deux zones n'apparaît pas sans défaut. D'une part, lorsque le nombre de Mach est supérieur à un, il n'y pas d'information voyageant à contre-courant, et toutes les ondes se déplacent vers l'aval. D'autre part, dans la région subsonique, des informations sont transportés dans les deux sens, vers l'amont et vers l'aval. Il existe donc une interface entre une maille supersonique et une maille subsonique où des informations supersoniques descendantes et des informations subsoniques remontantes sont mélangées d'une manière non physique. Or l'état stationnaire implique l'équilibre des flux à travers la frontière des cellules d'intégration. Pour satisfaire cette condition, les formulations volumes finis modifient le champ aérodynamique au centre des mailles. La solution à l'interface est alors correcte, au détriment de celle aux noeuds qui n'est plus physique. De ce fait, des oscillations très localisées peuvent apparaître.

Pour justifier cette analyse, il va être démontré par la suite qu'il n'est pas possible de conserver en même temps le flux à l'interface des cellules et aux noeuds, lorsque le régime d'écoulement change. Pour ce faire, nous considérons trois lignes de maillage extraites à différentes hauteur d'une tuyère transsonique, référencées par $j=1, j=8$ et $j=17$, pour un calcul en fluide parfait au premier ordre en espace. Les profils de nombres de Mach autour de la position du choc, sont présentés figure 3.20, alors que l'écart entre la température totale locale et la température totale imposée
<figure><img src="image_20.png" /><figcaption>Figure 3.18: Choc cylindrique. Niveaux de masse volumique et de pression</figcaption></figure>
<figure><img src="image_21.png" /><figcaption>Figure 3.19: Tuyère transsonique. Profils longitudinaux de pression totale</figcaption></figure>

<figure><img src="image_22.png" /><figcaption>Figure 3.20: Tuyère transsonique. Profils longitudinaux de nombre de Mach</figcaption></figure>
<figure><img src="image_23.png" /><figcaption>Figure 3.21: Tuyère transsonique. Profils longitudinaux de température totale</figcaption></figure>

en entrée, adimensionné par cette dernière, est visible figure 3.21. Cette grandeur est théoriquement nulle pour un écoulement adiabatique en fluide parfait. La discussion qui suit repose sur l'argumentation de Van Leer traitant de la gestion des discontinuités stationnaires que réalise sa méthode [36]. Les mêmes notations seront utilisées. Ainsi, l'état supersonique à l'amont du choc sera noté L, l'état subsonique à l'aval R, et les états intermédiaires P et Q, tels que la zone P soit supersonique et la zone Q subsonique. D'après la figure 3.21, nous constatons que seuls les noeuds P et Q présentent une erreur sensible. De plus, ce défaut est d'autant plus important que le noeud est proche de la ligne à Mach un. L'écart visible entre les états L et R est dû à la mauvaise précision spatiale, et n'est pas relié au sujet présentement discuté.

Les composantes du flux aux noeuds dans la direction de l'écoulement, $F^+$ et $F^-$, sont regroupées dans la table 3.1, ainsi que $F_l^+ + F_r^-$, le flux construit à chaque interface entre deux mailles. Pour un écoulement supersonique, $F^+$ est égale au flux complet $F$, et $F^-$ est nul.

Soit $\bar{F}$ le flux imposé par les conditions aux limites. Les états L et R sont supposés connus sans erreur tels que:

$$F(q_L) = F(q_R) = \bar{F} \qquad (3.61)$$

Cette hypothèse est valide si la différence de température totale entre les noeuds L et R est jugée négligeable (figure 3.21).

A l'interface PQ, nous avons:
<table><caption>Table 3.1: Composantes du flux aux noeuds et aux interfaces</caption><thead><tr><th>Noeud</th><th>F<sup>+</sup></th><th>F<sup>-</sup></th><th>F<sub>l</sub><sup>+</sup> + F<sub>r</sub><sup>-</sup></th></tr></thead><tbody><tr><td>L (M > 1)</td><td>F(q<sub>L</sub>)</td><td>0</td><td>F(q<sub>L</sub>)</td></tr><tr><td>P (M > 1)</td><td>F(q<sub>P</sub>)</td><td>0</td><td>F(q<sub>P</sub>) + F<sup>-</sup>(q<sub>Q</sub>)</td></tr><tr><td>Q (M < 1)</td><td>F<sup>+</sup>(q<sub>Q</sub>)</td><td>F<sup>-</sup>(q<sub>Q</sub>)</td><td>F<sup>+</sup>(q<sub>Q</sub>) + F<sup>-</sup>(q<sub>R</sub>)</td></tr><tr><td>R (M < 1)</td><td>F<sup>+</sup>(q<sub>R</sub>)</td><td>F<sup>-</sup>(q<sub>R</sub>)</td><td></td></tr></tbody></table>

$$F(q_P) + F^-(q_Q) = \bar{F} \qquad (3.62)$$

D'autre part au premier ordre,

$$F(q_R) = F^+(q_R) + F^-(q_R) = \bar{F} \qquad (3.63)$$

Si bien que la conservation du flux à l'interface QR,

$$F^+(q_Q) + F^-(q_R) = \bar{F} \qquad (3.64)$$

implique que:

$$F^+(q_Q) = F^+(q_R) \qquad (3.65)$$

Dans le cas d'un écoulement mono-dimensionnel, l'équation 3.65 se réduit de trois à deux équations scalaires indépendantes, le flux d'énergie de Van Leer étant fonction des flux de masse et de quantité de mouvement. Si bien que nous avons cinq équations, 3.62 et 3.65, pour six inconnues, $q_P$ et $q_Q$, admettant un paramètre qui est la localisation du choc entre les noeuds P et Q. Pour un écoulement multi-dimensionnel, la géométrie de la tuyère et les conditions aux limites d'entrée et de sortie, déterminent la position du choc. Les équations de quantité de mouvement dans les directions transversales remplacent alors le paramètre précédent. En trois dimensions, nous avons donc dix équations pour dix inconnues.

Si l'état P est tel que le flux soit conservé en ce noeud, $F(q_P) = \bar{F}$. D'où, d'après l'équation 3.62, $F^-(q_Q) = 0$, ce qui implique que $M(q_Q) > 1$. Or, ceci est en contradiction avec l'hypothèse selon laquelle la zone Q est subsonique. Aussi le flux ne peut être conservé en P, $F(q_P) \neq \bar{F}$, et par contrecoût en Q. Ainsi est-il établi qu'il n'est pas possible de conserver le flux numérique à l'interface des cellules et en même temps aux noeuds, à travers un choc.

Cette argumentation a été développée à partir de la séparation de flux de Van Leer. Cependant, sa conclusion est plus générale, et reste valable pour toute méthode construisant des états droit et gauche indépendants de chaque coté d'une interface, en fonction du nombre de Mach au centre des cellules. Il y aura toujours
génération d'erreurs au franchissement d'un choc stationnaire. D'autre part, ce défaut numérique semble occasionner une vibration longitudinale indésirable du choc. En effet, une instabilité transversale peut exister lorsqu'elle la position de la ligne à Mach un entre les noeuds P et Q change entre deux mailles successives dans la hauteur de la tuyère. L'erreur introduite alors n'étant pas régulièrement répartie, des oscillations parasites sont susceptibles d'apparaître.

## Non conservation de l'enthalpie totale

Hänel [26] a montré que l'enthalpie totale n'est pas conservée numériquement si l'expression des composantes du flux d'énergie développée par Van Leer est utilisée. Il propose à la place une autre formule, qui permet d'améliorer l'évolution de température totale, telle que:

$$f_{\text{energie}}^{\pm} = h_T \qquad (3.66)$$

où $h_T = E + \frac{p}{\rho}$ est l'enthalpie totale calculée à l'interface entre deux cellules à partir des champs aérodynamiques interpolés $Q^+$ et $Q^-$.

Cependant, avec cette formulation, le flux décomposé d'énergie ne change plus de signe pour un nombre de Mach égal à l'unité. Par conséquence, il y a risque d'étalement du choc dans le cas d'un calcul transsonique. En effet, le critère de la décomposition de Van Leer selon lequel une valeur propre nulle doit exister en régime subsonique n'est pas satisfait. Mais cette modification peut être intéressante pour des écoulements à faible nombre de Mach, près d'une paroi adhérente pour un fluide visqueux par exemple, comme l'a montré Hänel [26].

Pour notre part, nous avons obtenu une nette amélioration de la répartition de température totale au premier ordre dans une tuyère symétrique, traversée par un fluide parfait à des nombres de Mach compris entre 8.4 $10^{-2}$ et 9.6 $10^{-2}$ (figure 3.22). Nous constatons une augmentation (resp. une diminution) de $T_T$ lors de la compression (resp. la détente) si la technique de Van Leer est utilisée (figure 3.23). Par contre, il n'en est rien avec la correction de Hänel (figure 3.24), où l'écart entre la température totale imposée en entrée et la température totale locale calculée, est de l'ordre de la précision machine.

## Comportement non physique aux faibles nombres de Mach

La dissipation artificielle introduite par l'approche décentrée amont dépend du gradient local des variables conservatives, à travers l'évaluation des états gauche et droit $Q^-$ et $Q^+$, mais aussi de l'expression des composantes du flux $F^+$ et $F^-$ [21]. Pour preuve, considérons le bilan de flux de masse dans la direction $\xi^i$ par exemple.

Pour un schéma centré au second ordre, nous avons:

$$F_c^i(\xi^i + \frac{1}{2}) - F_c^i(\xi^i - \frac{1}{2}) = \frac{1}{2} [\rho \xi^i + 1 V_{\xi^i+1}^i - \rho \xi^i - 1 V_{\xi^i-1}^i] \qquad (3.67)$$

Alors qu'un schéma décentré amont au premier ordre conduit à:
<figure><img src="image_24.png" /><figcaption>Figure 3.22: Tuyère faiblement subsonique. Nombre de Mach</figcaption></figure>

<figure><img src="image_25.png" /><figcaption>Figure 3.23: Tuyère faiblement subsonique. Niveau de température totale. Van Leer</figcaption></figure>
<figure><img src="image_26.png" /><figcaption>Figure 3.24: Tuyère faiblement subsonique. Niveau de température totale. Hänel</figcaption></figure>

$$F_c^i(\xi^i + \frac{1}{2}) - F_c^i(\xi^i - \frac{1}{2}) = \rho_{\xi^i} V_{\xi^i}^i - \rho_{\xi^i-1} V_{\xi^i-1}^i \quad (3.68)$$

Par contre, la formulation de Van Leer donne:

$$F_c^i(\xi^i + \frac{1}{2}) - F_c^i(\xi^i - \frac{1}{2}) = \frac{1}{4} \left[ \rho_{\xi^i} a_{\xi^i} (M_{\xi^i+1}^i)^2 - \rho_{\xi^i+1} a_{\xi^i+1} (M_{\xi^i+1}^i - 1)^2 \right] \\ - \frac{1}{4} \left[ \rho_{\xi^i-1} a_{\xi^i-1} (M_{\xi^i-1}^i)^2 - \rho_{\xi^i} a_{\xi^i} (M_{\xi^i-1}^i - 1)^2 \right] \quad (3.69)$$

Prenons deux cas extrêmes où le nombre de Mach tend vers zéro ou vers un; si bien que nous trouvons les équivalents au premier ordre suivants:

$$(M^i + 1)^2 \approx 4M^i \quad \text{si} \quad M^i \to 1 \quad (3.70)$$

$$(M^i - 1)^2 \approx 0 \quad \text{si} \quad M^i \to 1 \quad (3.71)$$

$$\pm (M^i \pm 1)^2 \approx \pm 1 + 2M^i \quad \text{si} \quad M^i \to 0 \quad (3.72)$$

Reportés dans l'équation 3.69, ils nous permettent d'écrire:

$$F_c^i(\xi^i + \frac{1}{2}) - F_c^i(\xi^i - \frac{1}{2}) \approx \rho_{\xi^i} V_{\xi^i}^i - \rho_{\xi^i-1} V_{\xi^i-1}^i \quad \text{si} \quad M^i \to 1 \quad (3.73)$$

$$\approx \frac{1}{2} \underbrace{[\rho_{\xi^i+1} V_{\xi^i+1}^i - \rho_{\xi^i-1} V_{\xi^i-1}^i]}_{(a)}$$

$$+ \frac{1}{4} \underbrace{[\rho_{\xi^i} a_{\xi^i} - \rho_{\xi^i-1} a_{\xi^i-1}]}_{(b)}$$
<figure><img src="image_27.png" /><figcaption>Figure 3.25: Ecoulement de Poiseuille. Convergence vers l'état stationnaire</figcaption></figure>

$$ -\frac{1}{4} \underbrace{[\rho\xi^\prime + 1a\xi^\prime + 1 - \rho\xi^\prime a\xi^\prime]}_{(c)} \quad \text{si } M^\prime \to 0 \quad (3.74) $$

Ainsi, lorsque le nombre de Mach tend vers un, le schéma de Van Leer se rap-
proche d'un décentrement amont du flux complet (équation 3.68), ce qui est en
accord avec le sens physique. Par contre, lorsque le nombre de Mach tend vers zéro,
trois termes apparaissent. Le premier (a) correspond au bilan centré du flux com-
plet (equation 3.67), les deux autres (b) et (c) à la variation de pente du produit
ρa au point ξ'. Cette dérivée seconde constitue la viscosité numérique introduite
par la séparation de flux choisie. Si le produit ρa varie linéairement, (b) et (c) se
compensent et seul le terme centré demeure; à l'inverse, si le sens de variation de
ρa s'inverse en ξ', ces deux contributions se cumulent. Nous établissons ainsi que
la méthode de Van Leer génère, comparativement aux schémas non séparés, une
viscosité numérique qui varie avec le nombre de Mach.

Cette remarque permet d'expliquer les résultats peu physiques qui sont obtenus
lorsque le régime d'écoulement est subsonique faible, en particulier près d'une paroi
adhérente où la vitesse d'un fluide visqueux est proche de zéro [26, 40]. Dans ce cas,
les termes numériques (b) et (c) peuvent se révéler prépondérant devant l'information
physique (a), ρa étant très grand devant ρV' dans une telle région. Nous consta-
tons alors une bonne convergence du schéma, mais vers une solution médiocre, ces
deux phénomènes étant dus à une trop forte dissipation. A titre d'illustration, un
écoulement de Poiseuille à un nombre de Reynolds de 100 a été simulé en utilisant
la décomposition de Van Leer au troisième ordre d'une part, et un schéma centré
<figure><img src="image_28.png" /><figcaption>Figure 3.26: Ecoulement de Poiseuille. Erreur sur la pression statique</figcaption></figure>

<figure><img src="image_29.png" /><figcaption>Figure 3.27: Ecoulement de Poiseuille. Erreur sur la vitesse axiale</figcaption></figure>
précis au second ordre avec viscosité artificielle du quatrième ordre d'autre part (cf 3.4.4). L'évolution du résidu montre que la convergence est effectivement plus rapide dans le premier cas (figure 3.25). Cependant, le profil de vitesse théoriquement parabolique pour un écoulement incompressible établi [59, p. 84], et la perte de charge, théoriquement linéaire sous les mêmes conditions, sont déformés. Figure 3.26 est représenté l'écart entre les pressions statiques axiales théoriques et calculées, adimensionné par le niveau imposé en sortie ($X = 10$). Le résultat obtenu avec l'approche de Roe en maillage non structuré est fourni à titre de comparaison. Nous constatons tout d'abord que les trois courbes sont différentes du niveau analytique. Ceci peut s'expliquer par la nature faiblement compressible, et non totalement incompressible, du fluide simulé, le nombe de Mach maximum étant de 0.28. Au-delà de $X \approx 4$, le régime est considéré comme établi, et il est à noter qu'alors le schéma centré et celui de Roe conduisent à des solutions très proches, si l'erreur sur la pression statique en sortie du calcul non structuré est corrigée. Par contre, la séparation de flux de Van Leer produit une évolution de pression plus rapide, traduisant une perte de charge plus importante, semblable à celle obtenue avec un fluide d'avantage visqueux. Les écarts au profil transversal de vitesse dans la zone établie, adimensionnés par la vitesse axiale, laissent apparaître le même comportement (figure 3.27). La déformation du profil de vitesse obtenue par la méthode de Van Leer atteint 4% et ne peut être totalement imputée aux effets de compressibilité, qui n'affectent que de 1% les résultats centrés ou ceux de Roe. Nous mettons ainsi en évidence l'influence du choix des expressions des flux $F^+$ et $F^-$, qui déterminent l'évolution de la dissipation artificielle introduite par l'approche décentrée amont.

### Remarques complémentaires

La plupart de ces problèmes seraient évités si une séparation de flux exacte était utilisée. Elle permettrait alors l'évaluation de deux informations, $F^+$ et $F^-$, à partir de deux états discrets, $Q^-$ et $Q^+$, telles que le flux résultant de leur combinaison correspond exactement au flux physique. Or elle n'existe pas, et nous pouvons seulement faire en sorte que le bilan de flux s'approche de la réalité. La clef de ce problème est dans le choix du chemin numérique qui est suivi pour relier une interface entre deux cellules, aux centres des cellules qui la bordent. Ce chemin conditionne l'expression de la dissipation artificielle introduite. Dans cette perspective, la décomposition de Van Leer apparaît comme lourdement handicapée par le fait qu'elle repose trop peu sur des considérations physiques.

### 3.4.4 Schémas centrés

Les schémas centrés constituent une autre classe de méthodes numériques, qui contrairement aux schémas décentrés amont, sont plus bâties à partir de réflexions mathématiques. Communément utilisés pour des écoulements faiblements compresibles, ils présentent l'avantage d'être naturellement non dissipatifs, contrairement
aux schémas décentrés amont qui peuvent l'être excessivement. Cependant, ceci fait leur faiblesse, car ils se révèlent à l'usage fortement instables, en particulier en pré-sence de gradients importants. C'est pourquoi une dissipation artificielle doit leur être adjointe pour les stabiliser. Constitué généralement d'un terme du quatrième ordre, auquel se substitue un terme du second ordre à proximité des discontinuités, cet artifice est mal adapté aux régimes transsoniques, et non physique en écoulements supersoniques.

Cependant, les schémas décentrés amont sont très satisfaisants à ces nombres de Mach. Aussi nous ne leur recherchons une alternative qu'en régime subsonique faible, où les schémas centrés avec dissipation artificielle du quatrième ordre font merveilles. Voici comment une telle méthode numérique s'exprime en formulation MUSCL.

Le flux convectif à l'interface de deux cellules est évalué à l'aide de l'expression complète du flux, les variables conservatives étant quant à elles interpolées de manière centrée. Ainsi:

$$F_c^i(\xi^i - \frac{1}{2}) = F_c^i(Q_{\xi^i - \frac{1}{2}}) \qquad (3.75)$$

avec:

$$Q_{\xi^i - \frac{1}{2}} = \frac{1}{2} [q(\xi^i - 1) + q(\xi^i)] \qquad (3.76)$$

Il est à noter que cette expression est déjà incluse dans l'approche décentrée amont, où un choix judicieux de $\Phi = 1$ et $\kappa = 1$, permet de conserver le même formalisme.

D'autre part, est ajouté un terme dissipatif correspondant à une dérivée du troisième ordre des variables conservatives [30], tel que le flux convectif modifié $F_c^i*$ soit:

$$F_c^i*(\xi^i - \frac{1}{2}) = F_c^i(\xi^i - \frac{1}{2}) + \frac{\epsilon^i}{8} \sqrt{q} (V^i + a \|\vec{a}^i\|)_{\xi^i - \frac{1}{2}} [q_{\xi^i+1} - 3q_{\xi^i} + 3q_{\xi^i-1} - q_{\xi^i-2}] \qquad (3.77)$$

où $\epsilon^i \in [0; 2]$ est un coefficient à ajuster de manière à minimiser l'influence non physique de la dissipation artificielle tout en assurant la stabilité du schéma.

Si bien que le bilan de flux convectif dans la direction $\xi^i$ devient:

$$ \begin{aligned} & F_c^i*(\xi^i + \frac{1}{2}) - F_c^i*(\xi^i - \frac{1}{2}) \\ &= F_c^i(\xi^i + \frac{1}{2}) - F_c^i(\xi^i - \frac{1}{2}) + \frac{\epsilon^i}{8} \sqrt{q} (V^i + a \|\vec{a}^i\|)_{\xi^i} [q_{\xi^i+2} - 4q_{\xi^i+1} + 6q_{\xi^i} - 4q_{\xi^i-1} + q_{\xi^i-2}] \end{aligned} \qquad (3.78) $$

Afin de bien séparer les effets convectifs et dissipatifs, qu'ils soient physiques ou numériques, le second terme de cette équation, qui est une dérivée du quatrième
ordre des variables conservatives, sera par la suite considéré comme faisant partie du bilan des flux diffusifs et non plus du bilan des flux convectifs.

### 3.4.5 Schémas mixtes

Nous avons vu qu'un schéma décentré amont permettait une simulation satisfaisante des écoulements transsoniques et supersoniques, alors qu'un schéma centré était bien adaptés aux régimes faiblements compressibles. Une idée intéressante est de pouvoir utiliser chacun dans son domaine d'application optimum, une transition progressive de l'un à l'autre minimisant tout problème potentiel de raccord.

Sa mise est oeuvre se révèle simple dans une approche MUSCL. En effet, il suffit de faire évoluer le décentrement de l'interpolation des variables aux demi-mailles en fonction du nombre de Mach, l'expression des contributions amont et aval des flux restant par ailleurs inchangée. Une technique semblable a déjà été utilisée lors du traitement des discontinuités où la précision spatiale est ramené au premier ordre localement, si ce n'est que dans ce cas, nous n'agissons plus sur le limiteur $\Phi$ mais sur le paramètre $\kappa$. Ainsi, nous posons:

$$F_c^i(\xi^i - \frac{1}{2}) = F_c^{i+}(Q_{\xi^i - \frac{1}{2}}^-) + F_c^{i-}(Q_{\xi^i + \frac{1}{2}}^+) \quad (3.79)$$

où:

$$Q_{\xi^i - \frac{1}{2}}^- = Q(q_{\xi^i}, q_{\xi^i - 1}, q_{\xi^i - 2}) \quad (3.80)$$

$$Q_{\xi^i - \frac{1}{2}}^+ = Q(q_{\xi^i - 1}, q_{\xi^i}, q_{\xi^i + 1}) \quad (3.81)$$

avec:

$$Q(q_0, q_1, q_2) = \alpha^{\pm} Q^{\pm}(q_0, q_1, q_2) + (1-\alpha^{\pm}) Q^c(q_0, q_1, q_2) \quad (3.82)$$

et:

$$Q^{\pm}(q_0, q_1, q_2) = \frac{1}{4}\Phi[1 + \kappa\Phi]q_0 + \left[1 - \frac{1}{2}\kappa\Phi^2\right]q_1 - \frac{1}{4}\Phi[1 - \kappa\Phi]q_2 \quad (3.83)$$

$$Q^c(q_0, q_1, q_2) = \frac{1}{2}[q_0 + q_1] \quad (3.84)$$

$\alpha^{\pm}$ est une fonction du nombre de Mach calculé à partir du champ aérodynamique centré $Q^c$, qui assure la pondération entre les deux schémas. L'expression que nous utilisons fait intervenir un nombre de Mach de raccord $M_0^i$, au delà duquel seul le schéma décentré amont est utilisé:

$$ \begin{aligned} \alpha^{\pm} &= \frac{1}{2} \left[ 1 - \cos \left( \pi \frac{M^i}{M_0^i} \right) \right] & \text{si } |M^i| \le M_0^i \\ &= 1 & \text{si } |M^i| \ge M_0^i \end{aligned} \quad (3.85) $$
Cette fonction paire assure la continuité de $\alpha^{\pm}$ et de sa dérivée première au point de raccord et à l'origine. D'autre part, elle est telle que le schéma soit totalement centré lorsque $M^i = 0$.

La dissipation artificielle requise par le schéma centré subit la même pondération que $Q^c$, le schéma décentré amont incorporant ses propres effets dissipatifs. Aussi le coefficient $\epsilon^i$ apparaissant dans l'équation 3.78 est-il remplacé par $(1 - \alpha^{\pm}) \epsilon^i$.

## 3.5 Discrétisation spatiale des flux diffusifs

Rappelons que $F_v^i$, le flux diffusif à travers une face de la cellule d'intégration dans la direction $\xi^i$, est défini par (équation 3.28):

$$F_v^i = \sqrt{g} \vec{F}_v \cdot \vec{a}^i \qquad (3.86)$$

où:

$$\vec{F}_v = \begin{cases} 0 & \vec{r}_{l+t} \\ \vec{r}_{l+t} & \gamma \left( \frac{\mu}{F_r} + \frac{\mu_t}{F_{r,t}} \right) \vec{\nabla} e + \vec{r}_{l+t} \cdot \vec{V} \end{cases} \qquad (3.87)$$

$$\vec{r}_{l+t} = (\lambda + \lambda_t) \vec{\nabla} \cdot \vec{V} \vec{I} + (\mu + \mu_t) (\vec{\nabla} \vec{V} + \vec{\nabla} \vec{V}^t) \qquad (3.88)$$

$F_v^i$ est donc fonction des dérivées premières de l'énergie interne et de la vitesse, qui s'expriment comme suit en coordonnées curvilignes (annexe D.6):

$$\vec{\nabla} u_k = \sum_{j=1}^{3} \vec{a}^j \frac{\partial u_k}{\partial \xi^j} \qquad (3.89)$$

$$\vec{\nabla} \cdot \vec{V} = \sum_{j=1}^{3} \vec{a}^j \cdot \frac{\partial \vec{V}}{\partial \xi^j} \qquad (3.90)$$

$\vec{\nabla} \vec{V}$ (resp. $\vec{\nabla} \vec{V}^t$) est formé des vecteurs $\vec{\nabla} u_k$ écrits en ligne (resp. en colonne).

Ainsi l'expression du flux diffusif fait intervenir des termes géométriques, le champ aérodynamique et ses dérivées, qui tous doivent être connus à l'interface considérée. Dans un soucis de précision, le jacobien $\sqrt{g}$ et la base contravariante $\vec{a}^i$ sont calculés de manière exacte en cette demi-maille, et non pas interpolés. Par contre, les variables conservatives $q$ qui sont stockées aux noeuds, i.e. au centre des cellules, y sont évaluées par une moyenne centrée entre les états amont et aval, contrairement au flux convectif où une moyenne décentrée est utilisée. Nous posons donc:

$$Q_{\xi^i - \frac{1}{2}} = \frac{1}{2} [q(\xi^i - 1) + q(\xi^i)] \qquad (3.91)$$
<figure><img src="image_30.png" /><figcaption>Figure 3.28: Molécule de calcul des flux visqueux</figcaption></figure>

Les autres grandeurs, telles que les composantes de la vitesse ou l'énergie interne, en sont déduites. Enfin, les dérivées spatiales d'un scalaire A (équation 3.89) sont discrétisées en suivant une formulation différence finie en espace curviligne, d'où (figure 3.28):

$$ \bar{\nabla} A|_{\xi^i = -\frac{1}{2}} = \bar{a}^i|_{\xi^i = -\frac{1}{2}} [A(\xi^i) - A(\xi^i - 1)] + \sum_{\substack{j=1 \\ j \neq i}}^{3} \bar{a}^j|_{\xi^i = -\frac{1}{2}} \frac{A(\xi^j + 1) - A(\xi^j - 1)}{2} \quad (3.92) $$

Nous avons ainsi accès à toutes les données nécessaires à l'évaluation de $\bar{\nabla} \bar{V}$, $\bar{\nabla} \bar{V}^t$ et $\bar{\nabla} \cdot \bar{V}$, puis de $\vec{\tau}_{l+1}$ et $\bar{\nabla} e$, et enfin de $\vec{F}_v$ et de $F_v^i$.

## 3.6 Modélisation de la turbulence

Afin d'assurer la fermeture du système différentiel décrivant le mouvement moyen du fluide, un certain nombre d'hypothèses simplificatrices ont été faites (paragraphe 2.4). Rappelons qu'elles établissent que l'énergie cinétique turbulente moyenne $\tilde{k}$ est négligeable, et que la viscosité turbulente $\mu_t$ est une fonction algébrique d'une échelle de longueur $l_m$ et d'une échelle de vitesse, bâtie sur le gradient de vitesse moyenne. Deux modèles du type longueur de mélange explicitant cette fonction sont décrits par la suite. Reposant sur des hypothèses simples, du type couche limite bidimensionnelle turbulente établie en plaque plane, ils n'incluent aucune correction ou artifice visant à gérer la nature instationnaire de l'écoulement macroscopique. Cependant, aucun autre modèle réaliste, plus élaboré, proposé jusqu'à ce jour, n'en est capable. C'est pourquoi notre choix s'est porté vers une solution simple, offrant
<figure><img src="image_31.png" /><figcaption>Figure 3.29: Evolution de la longueur de mélange</figcaption></figure>

une bonne base de travail et d'extension, et d'un faible coût de mise en œuvre et d'exploitation. Cette partie traite d'un des aspects les plus importants et cependant les moins bien maîtrisés de la mécanique des fluides numérique. Un important travail de recherche est encore nécessaire.

### 3.6.1 Modèle de Prandtl

Prandtl [51] a proposé une expression de $\mu_t$, proportionnelle à une vitesse fluctuante moyenne $U$ et à une longueur, dite de mélange, $l_m$. Considérant un écoulement cisaillé ne comprenant qu'une seule composante du tenseur viqueux non nulle ($\overline{u_1u_2}$) et qu'une seule composante du gradient de vitesse non nulle ($\frac{\partial u_1}{\partial x_2}$), il a alors postulé que $U$ était égale au produit du gradient de vitesse moyenne et de la longueur de mélange [59, pp. 579-583]:

$$U = l_m \left| \frac{\partial u_1}{\partial x_2} \right| \qquad (3.93)$$

Si bien que la viscosité turbulente s'exprime alors:

$$\mu_t = \rho l_m^2 \left| \frac{\partial u_1}{\partial x_2} \right| \qquad (3.94)$$

L'étape suivante consiste en la détermination de la longueur de mélange. Dans le cas d'une couche limite bidimensionnelle, en plaque plane, sans gradient de pression, une fonction rampe telle que celle tracée figure 3.29, fournit une évolution satisfaisante de $l_m$. Patankar et Spalding [50] ont proposé pour valeurs des constantes arbitraires $K$ et $\lambda$:
$$K = 0.435 \qquad (3.95)$$

$$\lambda = 0.09 \qquad (3.96)$$

Ce modèle est sans doute un des plus simples qui soit. A ce titre, une grande précision ne saurait en être attendue, car il ne prend pas en compte, entre autres, la présence d'une zone très près de la paroi où l'écoulement à très faible vitesse est laminaire, connue sous le nom de sous-couche visqueuse laminaire. Cependant, il a été rappelé ici car il est à la base de nombreux autres modèles, tels que ceux de Cebeci-Smith [12] et de Baldwin-Lomax [6].

### 3.6.2 Modèle de Michel

L'expression de $\mu_t$ proposée par Michel, Quemard et Durant [42] s'inspire directement de ces travaux. Elle vise à en généraliser la formulation et à en corriger certains défauts. Ainsi, elle établie que la viscosité turbulente est définie par:

$$\mu_t = \rho l_m^2 F^2 \left| \frac{\partial V_t}{\partial n} \right| \qquad (3.97)$$

où:

* $l_m$ est la longueur de mélange, représentant l'échelle de longueur.

* $l_m \left| \frac{\partial V_t}{\partial n} \right|$, produit de la longueur de mélange et de la norme de la dérivée suivant la direction normale à la paroi de la composante tangentielle de la vitesse, est l'échelle de vitesse.

* F est une fonction d'amortissement, qui assure une décroissance rapide de $\mu_t$ très près de la paroi afin de retranscrire la sous-couche visqueuse laminaire.

Pour généraliser ce modèle à un écoulement tridimensionnel, Cambier et Escande [11] ont proposé de remplacer $\left| \frac{\partial V_t}{\partial n} \right|$ par la norme du rotationnel de vitesse, $\| \bar{\omega} \|$. Ainsi, l'équation 3.97 devient:

$$\mu_t = \rho l_m^2 F^2 \| \bar{\omega} \| \qquad (3.98)$$

La longueur de mélange $l_m$ est définie par:

$$l_m = 0.085 \delta \operatorname{th} \left( \frac{K}{0.085} \frac{d}{\delta} \right) \qquad (3.99)$$

où:

* $K = 0.41$ est la constante de Van Karman.
* $\delta$ est l'épaisseur de la couche limite.

* $d$ est la distance à la paroi du point considéré.

La tangente hyperbolique assure une transition douce entre les deux parties de la rampe 3.29, où:

$$ 
\begin{aligned}
l_m &\approx Kd \quad \text{si} \quad d \ll \delta \\
l_m &\approx 0.085\delta \quad \text{si} \quad d \approx \delta
\end{aligned} 
\quad (3.100) $$

La fonction F est telle que:

$$ F(\zeta) = 1 - \exp \left[ - \frac{\sqrt{\zeta}}{26 K} \right] \quad (3.101) $$

avec:

$$ \zeta = \rho l_m^2 \frac{\mu + \mu_t}{\mu^2} \| \bar{\omega} \| \quad (3.102) $$

Il est à noter que cette expression rend implicite la définition de $\mu_t$. Un processus itératif entre les équations 3.102, 3.101 et 3.98, $\rho$, $l_m$ et $\| \bar{\omega} \|$ étant par ailleurs fixés, est donc nécessaire. Par expérience, il faut alors entre une et quatre itérations au maximum pour obtenir une valeur stable de $\mu_t$.

Comme pour beaucoup de modèle de ce type, la pierre d'achoppement en est l'évaluation de l'épaisseur de la couche limite $\delta$. Ce problème est résolu ici de la manière suivante. Le long d'une normale à la paroi, le point en lequel la norme du rotationnel devient inférieure à un pourcentage fixé du rotationnel maximum, est considéré en frontière de couche limite. $\delta$ est donc telle que:

$$ \| \bar{\omega} \|_{d=\delta} \le \beta \| \bar{\omega} \|_{\max} \quad (3.103) $$

La valeur optimale de $\beta$ varie selon l'écoulement. Elle se situe cependant entre 2% et 5%.

Enfin, si la définition de $d$ est non ambigüe dans le cas d'une plaque plane isolée, il en est tout autre dans un canal. En présence de deux parois en vis-à-vis, configuration adaptée à l'étude d'une tuyère bidimensionnelle ou d'une grille d'aubes, nous adoptons une expression inspirée des travaux de Buleev [43]:

$$ d = \frac{d_1 d_2}{d_1 + d_2} \quad (3.104) $$

où $d_1$ et $d_2$ représentent les distances du point de calcul à chacune des parois.
<figure><img src="image_32.png" /><figcaption>Figure 3.30: Modèle de Michel. Profils de vitesse</figcaption></figure>

Ce modèle présente l'avantage de la simplicité, et peut même être étendu assez facilement au cas des écoulements turbulents tridimensionnels avec plusieurs parois, comme ceux se développant dans un canal inter-aubes par exemple [11]. Cependant, il repose sur un certain nombre d'hypothèses sur la géométrie (plaque plane) et le champ aérodynamique (écoulement établi) qui ne sont pas vérifiées en turbomachine. De plus, il recourt à des constantes tabulées dans des configurations simples, assez éloignées de la réalité. Une description précise des phénomènes turbulents ne serait donc en être attendue. Mais il permet en pratique d'introduire à faible coût une grande majorité de leur influence, ce qui est le but que nous poursuivons.

### 3.6.3 Comparaisons et critiques

Dans le but d'estimer la fiabilité du modèle de Michel, nous l'avons comparé d'une part à celui de Cebeci-Smith, d'autre part à des résultats expérimentaux [61]. L'expérience servant de référence [58] met en oeuvre une couche limite incompressible stationnaire se développant sur une plaque plane en présence d'un gradient adverse de pression, qui la ralentit sans qu'elle ne décolle cependant. Nous considérons par la suite deux stations de mesure, une amont et une autre aval, en laquelle le gradient de pression déforme sensiblement le profil de vitesse (figure 3.30). La répartition de vitesse expérimentale peut être approchée par une loi analytique proposée par Coles [32] (figures 3.31 et 3.32). Cette dernière nous fournit le champ de vitesse car nous n'effectuons pas ici un calcul couplé de la viscosité turbulente et du champ aérodynamique, qui est supposé constant et fixé en conséquence.
<figure><img src="image_33.png" /><figcaption>Figure 3.31: Modèle de Michel. Profil de vitesse amont</figcaption></figure>

<figure><img src="image_34.png" /><figcaption>Figure 3.32: Modèle de Michel. Profil de vitesse aval</figcaption></figure>
<figure><img src="image_35.png" /><figcaption>Figure 3.33: Modèle de Michel. Contrainte de cisaillement amont</figcaption></figure>

<figure><img src="image_36.png" /><figcaption>Figure 3.34: Modèle de Michel. Contrainte de cisaillement aval</figcaption></figure>
Figures 3.33 et 3.34 sont tracées les contraintes de cisaillement. Un bon accord est constaté à l'amont, où la couche limite n'est pas décélérée, entre les trois approches (figure 3.33). Le pic visible sur la solution de Cebeci-Smith est dû à un mauvais raccord entre les deux branches de la rampe décrivant la longueur de mélange $l_m$. L'utilisation d'une expression en tangente hyperbolique dans le modèle de Michel évite un tel comportement. Par contre à l'aval, où le gradient adverse de pression à un effet sensible, les deux méthodes numériques surestiment d'un facteur trois la contrainte de cisaillement (figure 3.34).

Il apparaît donc que lorsque l'écoulement s'éloigne de trop des hypothèses faites par Prandtl, les modèles bâtis sur ses travaux se révèlent peu fiables. Par contre, dans le cadre de ces approximations, ils sont sans reproche. Cette information est importante et doit être gardée à l'esprit lorsque des configurations complexes sont simulées. Cependant, aucune comportement franchement abherrant n'a été mis en évidence lors de ces tests, et nous jugeons qu'une utilisation pragmatique et circonspecte du modèle de Michel est possible dans un premier temps, même si des aménagements visant à intégrer l'effet de la courbure et de la pression sont souhaitables.

## 3.7 Résolution numérique

Finalement, nous avons transformé le système différentiel initial en un ensemble d'équations discrètes en temps et en espace dont la solution décrit l'évolution du champ aérodynamique. Il s'écrit:

$$ \Delta (\sqrt{g}q)^n = -\Delta\tau R^n \quad (3.105) $$

avec:

$$ R = \sum_{i=1}^{3} \left[ F^i(\xi^i + \frac{1}{2}) - F^i(\xi^i - \frac{1}{2}) \right] \quad (3.106) $$

$$ F^i = F_c^i - \frac{1}{R_c} F_v^i \quad (3.107) $$

L'expression des flux convectifs $F_c^i$ et diffusifs $F_v^i$ ayant fait l'objet des paragraphes précédents, il ne reste qu'à résoudre numériquement ce système.

### 3.7.1 Progression en temps

L'équation 3.105 décrit l'évolution temporelle de $\sqrt{g}q$, qui inclut deux informations que nous considérons découplées par ailleurs, une première géométrique, la seconde aérodynamique. Si bien que cette grandeur est transformée en:

$$ \begin{aligned} \Delta (\sqrt{g}q)^n &= \sqrt{g}^{(n+1)} q^{(n+1)} - \sqrt{g}^{(n)} q^{(n)} \\ &= \sqrt{g}^{(n+1)} [q^{(n+1)} - q^{(n)}] + [\sqrt{g}^{(n+1)} - \sqrt{g}^{(n)}] q^{(n)} \end{aligned} $$
$$= \sqrt{g}^{(n+1)} \cdot \Delta q^n + \Delta(\sqrt{g})^n \cdot q^n \quad (3.108)$$

Ainsi, l'équation 3.105 devient:

$$\sqrt{g}^{(n+1)} \cdot \Delta q^n = - \left[ \underbrace{\Delta \tau R^n - \sqrt{g}^{(n)} q^{(n)}}_{(a)} + \underbrace{\sqrt{g}^{(n+1)} q^{(n)}}_{(b)} \right] \quad (3.109)$$

La progression en temps, de l'instant $(n)$ à l'instant $(n+1)$, durant l'intervalle de temps $\Delta\tau$, se déroule alors comme suit:

* Le champ aérodynamique, $q$, et les métriques, $\tilde{a}$ et $\sqrt{g}$, étant connus à l'instant $(n)$, sont calculés les flux $F_c^i$, $F_v^i$, puis $F^i$ à l'interface des cellules, et enfin le résidu $R^n$ aux noeuds. Le terme (a) est alors évalué.

* Puis le maillage est modifié par déplacement des noeuds, ce qui permet le calcul de la vitesse de l'enveloppe des cellules $\tilde{W}$, et des métriques à l'instant suivant, en particulier $\sqrt{g}^{(n+1)}$. Le terme (b) peut alors être ajouté à (a).

* Enfin, il ne reste qu'à écrire que:

$$q^{(n+1)} = q^{(n)} - \frac{(a) + (b)}{\sqrt{g}^{(n+1)}} \quad (3.110)$$

* Nous disposons alors du champ aérodynamique et des métriques à l'instant $(n+1)$.

### 3.7.2 Evaluation de la convergence

Lors de la recherche d'un état stationnaire, la progression en temps se poursuit tant que la solution évolue. Dans la mesure où un calcul jusqu'à une date infinie est irréalisable, un critère doit être déterminé afin d'évaluer la qualité de l'état convergé. Il en existe de nombreux, bâts principalement sur le résidu. En effet, pour une composante du vecteur inconnu $q$, nous avons à convergence théorique:

$$\Delta (\sqrt{g}q)^n = 0 \quad (3.111)$$

et donc:

$$R^n = 0 \quad (3.112)$$

Le zéro parfait étant une réalité mathématique, mais un mythe numérique, cette dernière égalité est remplacée par:

$$\| R^n \| < \| R^n \|_{\min} \quad (3.113)$$

où:
<table><caption>Figure 3.35: Répartition statistique du résidu à convergence</caption><thead><tr><th colspan="3">Residu sur l'energie: L0=</th><th colspan="2">-.368 L1= -1.718 L2= -3.282</th></tr><tr><th>Residu</th><th>%</th><th>% cumule</th><th>Histogramme</th><th></th></tr></thead><tbody><tr><td>-3.282</td><td>2.392</td><td>2.392</td><td>**</td><td></td></tr><tr><td>-3.091</td><td>4.941</td><td>7.333</td><td>*****</td><td></td></tr><tr><td>-2.767</td><td>39.212</td><td>46.545</td><td>*****************************</td><td></td></tr><tr><td>-2.443</td><td>3.596</td><td>50.141</td><td>****</td><td></td></tr><tr><td>-2.119</td><td>10.163</td><td>60.303</td><td>*********</td><td></td></tr><tr><td>-1.795</td><td>16.807</td><td>77.111</td><td>*****************</td><td></td></tr><tr><td>-1.472</td><td>12.914</td><td>90.025</td><td>*********</td><td></td></tr><tr><td>-1.148</td><td>5.191</td><td>95.216</td><td>*****</td><td></td></tr><tr><td>-.824</td><td>4.190</td><td>99.406</td><td>****</td><td></td></tr><tr><td>-.692</td><td>.594</td><td>100.000</td><td>*</td><td></td></tr></tbody></table>

* $|| \cdot ||$ désigne une norme du vecteur $R^n$, qui comporte rappelons-le $N$ composantes, si $N$ est le nombre de noeuds.

* $|| R^n ||_{min}$ est une valeur seuil de la norme du résidu en dessous de laquelle nous considérons que l'état convergé est atteint.

La norme d'un vecteur peut être définie de nombreuses manières. Les suivantes en sont trois classiques:

$$L_0^n = \max \{|R_0^n|, i \in \{1; N\}\} \qquad (3.114)$$

$$L_1^n = \frac{1}{N} \sum_{i=1}^{N} |R_1^n| \qquad (3.115)$$

$$L_2^n = \frac{1}{N} \sqrt{\sum_{i=1}^{N} |R_2^n|^2} \qquad (3.116)$$

Par expérience, nous avons constaté qu'aucune des trois n'était suffisante. En effet, si tout le champ aérodynamique est convergé excepté un seul noeud, $L_0^n$ sera élevée alors que $L_1^n$ et $L_2^n$ seront faibles. Au contraire, il est possible que $L_0^n$ soit modérée, sans pour autant que $L_1^n$ et $L_2^n$ le soient, pour peu que de nombreux noeuds soient encore corrigés, même faiblement.

D'autre part, la convergence de chacune des cinq composantes de $q$ doit être vérifiée, car il apparaît que selon le régime d'écoulement, le comportement des unes par rapport aux autres change. En effet, si en transsonique toutes les variables évoluent plus ou moins en phase, à bas nombres de Mach, la masse volumique et l'énergie se singularisent par rapport à la quantité de mouvement.

Cependant, par expérience, c'est l'énergie totale qui semble être la plus sensible à tout écart. C'est pourquoi nous contrôlons la qualité de l'état convergé sur cette variable en construisant un histogramme de répartition du résidu entre $L_0^n$ et $L_2^n$ en
échelle logarithmique. Cette courbe, dont un exemple est visible figure 3.35, nous renseigne sur le pourcentage de points ayant atteints un certain niveau de convergence, et permet d'estimer l'importance des noeuds encore corrigés.

A ce problème de choix d'une norme appropriée, se superpose celui du seuil de référence $||R^n||_{min}$. Une technique classique consiste à adimensionner la norme du résidu à l'instant $(n)$ par sa valeur initiale [25]. Une décroissance relative d'un certain ordre de grandeur fixe est alors requise. Le piège de cette approche réside dans la composition du champ initial. S'il est très éloigné de la solution convergée, la norme initiale sera élevée et il est donc aisé de la faire rapidement décroître. Par contre, si l'état initial est peu différent de celui convergé, un gain même faible nécessitera de très nombreuses itérations.

C'est pourquoi cette approche n'a pas été retenue dans cette étude. Nous préférons poursuivre le calcul jusqu'à ce qu'un niveau sensiblement constant de la norme du résidu soit atteint, traduisant un plateau de convergence stable. Il faut alors évaluer la position de ce dernier. Sa valeur théorique est zéro. Cependant, la précision des calculs sur ordinateur étant finie, un incrément $\Delta (\sqrt{gq})^n$ significatif de la solution ne serait être inférieur à cette dernière. Remarquant que:

$$ ||R^n|| = \frac{||\Delta (\sqrt{gq})^n||}{\Delta \tau} \qquad (3.117) $$

et en notant $\epsilon$ la précision machine, telle que:

$$ ||\Delta (\sqrt{gq})^n|| \ge \epsilon \qquad (3.118) $$

nous en déduisons que:

$$ ||R^n|| \ge \frac{\epsilon}{\Delta \tau} \qquad (3.119) $$

La valeur minimale du résidu est donc $\frac{\epsilon}{\Delta \tau}$. Si elle est atteinte, nous déclarons la solution convergée au zéro machine. Un fait intéressant est que lors de l'adimensionalisation des équations, la longueur de référence $L_0$ n'intervient qu'au niveau des termes temporels, et donc de $\Delta \tau$. Il en résulte qu'une variation de $L_0$ changera artificiellement la position du plateau théorique, sans modifier pour autant la qualité de la solution convergée. Un changement du CFL aura le même effet. La position absolue du plateau est donc illusoire et non significative. Seules sa nature et sa position par rapport à la valeur théorique sont probantes.

## 3.8 Conclusion

Ce chapitre a été consacré à la description des méthodes numériques mises en oeuvre pour résoudre les équations différentielles initiales. Elles ont été tout d'abord intégrées en espace dans le cadre d'une formulation volumes finis en maillage curviligne structuré, mobile et déformable. Puis leur intégration en temps a été faite
explicitement. Les flux convectifs ont été alors évalués par une approche MUSCL, couplée à un schéma mixte, résultant de l'hybridation en fonction du nombre de Mach d'un schéma centré et d'une technique de séparation de flux proposée par Van Leer. Une étude critique de cette dernière a par ailleurs été présentée. Les termes diffusifs ont été quant à eux discrétisés par des différences centrées au second ordre. Enfin, les effets de la turbulence ont été introduits au moyen d'un modèle algébrique, du type longueur de mélange.

La dernière étape de la résolution numérique est la prise en compte des conditions aux limites. Celles-ci conditionnent l'existence et l'unicité de la solution calculée, et leur traitement influence fortement la qualité du champ aérodynamique obtenu. Il fait à ce titre l'objet du chapitre suivant.
# Chapitre 4

## Traitement des conditions aux limites

### 4.1 Introduction

L'objectif des développements présentés par la suite est la gestion des frontières du domaine de calcul. Nous considérons alors les noeuds qui le limitent, par opposition à ceux qui en forment le coeur. Une attention toute particulière doit être portée à ce traitement. En effet, rappelons que nous ne nous intéressons qu'à une portion de l'espace. Or dans la réalité, le reste de l'univers peut interagir avec cette boite que nous étudions. Il faut donc prendre en compte son environnement, jusqu'à l'infini, qui n'est malheureusement pas connu en détail car étant en dehors du domaine de calcul, il n'est pas discrétisé. Une mauvaise gestion de cette influence peut avoir des conséquences catastrophiques sur la solution globale obtenue.

Nous cherchons à atteindre plusieurs buts à travers ce traitement. Les uns sont physiques. Il s'agit en particulier d'imposer des contraintes extérieures sur l'écoulement, telles que la forme des parois guidant le fluide ou les fluctuations de certaines grandeurs, ainsi qu'un régime d'écoulement, qui peut être subsonique, transsonique, supersonique, stationnaire, instationnaire, ... Ceci est réalisé à travers la géométrie et des variables physiques, comme la pression statique, la température totale, ... D'autres objectifs poursuivis sont numériques. En effet, le schéma de discrétisation utilisé dans le coeur du domaine met en oeuvre, dans une direction, deux voisins de chaque côté d'un noeud, constituant ainsi une molécule de calcul tridimensionnelle à treize atomes (figure 4.1). Or un noeud frontière et un noeud juste intérieur n'ont pas un voisinage assez étendu, ayant respectivement zéro (figure 4.3) et un seul voisin (figure 4.2) dans la direction perpendiculaire à la surface limite. Deux solutions existent à ce problème: soit les noeuds manquants sont créés, soit un autre schéma, utilisant une molécule différente, doit être utilisé.

Afin de satisfaire à ces contraintes, nous définissons un traitement local, en fonction des conditions aux limites imposées au noeud frontière considéré. Trois con-
<figure><img src="image_37.png" /><figcaption>Figure 4.1: Molécule de calcul pour un noeud du coeur</figcaption></figure>

<figure><img src="image_38.png" /><figcaption>Figure 4.2: Molécule de calcul pour un noeud juste intérieur</figcaption></figure>

<figure><img src="image_39.png" /><figcaption>Figure 4.3: Molécule de calcul pour un noeud frontière</figcaption></figure>
figurations sont distinguées. La première, désignée par la suite comme mettant en oeuvre une frontière géométrique, concerne les noeuds appartenant à une surface limite du domaine de calcul présentant des propriétés géométriques remarquables, un plan de symétrie par exemple. Les deux autres, référencées comme étant des conditions aux limites fortes et des conditions aux limites faibles, sont appliquées à des noeuds en lesquels une contrainte extérieure est imposée. La distinction entre ces deux approches repose sur la manière dont ceci est réalisé. Dans la première, une condition physique, la tangence de la vitesse du fluide à la paroi par exemple, est introduite dans le système à résoudre en un noeud comme une équation à satisfaire. Elle remplace ainsi une équation du système initial. La condition est imposée à travers une valeur, ou à travers sa dérivée temporelle, ce qui permet de gérer des conditions aux limites instationnaires. D'un point de vue global, nous dirons que c'est la cause qui est imposée. Par contre, dans le cas du traitement faible, c'est l'effet qui est pris en compte, par exemple le fait que le gradient de pression statique perpendiculairement à une paroi de glissement plane soit nul. Toutes les équations initiales sont conservées, mais certains flux sont a priori fixés. La condition est ainsi imposée à travers une dérivée spatiale. A convergence d'un état stationnaire, les deux approches sont équivalentes, ce qui est plus douteux en instationnaire si aucun soin particulier n'est apporté au développement de la méthode faible.CACIFICATURALEMENT, nous pouvons dire que nous avons d'une part une technique dure, inflexible, le niveau imposé étant théoriquement satisfait, d'autre part une approche molle, tolérante, faisant preuve d'une plus grande souplesse. Le choix de l'une ou de l'autre se fera en fonction des besoins de stabilité ou de précision des techniques implémentées.

Dans tous les cas, une contrainte numérique majeure doit être prise en compte. Il s'agit de la compatibilité des schémas développés pour traiter les frontières avec celui mis en oeuvre dans le coeur, et qui repose sur un système de cinq équations différentielles du premier ordre en temps. Ceci afin d'utiliser en tout point la même discrétisation temporelle et d'éviter ainsi d'éventuels déphasages entre les différentes régions du domaine de calcul.

## 4.2 Conditions géométriques

Ce traitement s'applique aux noeuds appartenant à une frontière du domaine de calcul présentant une propriété géométrique remarquable. Il s'agit par exemple d'une symétrie plane, qui sépare une tuyère symétrique au niveau de son axe en deux demi-tuyères équivalentes, ou encore d'une périodicité angulaire, qui nous permet, en turbomachine, de ne considérer qu'un seul canal inter-aubes, en supposant que tous les autres formant la roue complète sont identiques à une rotation près à celui que nous simulons.

Dans les développements suivants, nous ferons l'hypothèse que l'écoulement, ainsi que les conditions aux limites imposées, satisfont aux mêmes propriétés géométriques
<figure><img src="image_40.png" /><figcaption>Figure 4.4: Plans de maillage supplémentaires</figcaption></figure>

que le domaine de calcul.

### 4.2.1 Extension du domaine de calcul

La clef de ce traitement est dans la création de nouveaux noeuds à partir de ceux existants, à l'aide des relations géométriques, par exemple par symétrie par rapport à un plan ou par translation angulaire. Nous construisons ainsi deux nouveaux plans de maillage, dits *supplémentaires*, sur chacune des six faces du domaine de calcul (figure 4.4). Ce nombre n'est pas innocent puisqu'il permet d'accroître le voisinage d'un noeud frontière jusqu'à reconstituer une molécule de calcul complète, identique à celle du coeur.

La stratégie présidant à cette extension du domaine de calcul est simple. Toutes les informations nécessaires en ces nouveaux noeuds sont calculées à partir de celles connues dans le coeur et à la frontière. Ceci vaut pour leur position *OM* et la base covariante $\bar{a}_i$, à partir de laquelle seront déduits le jacobien de la transformation géométrique $\sqrt{g}$ et la base contravariante $\bar{a}^i$. Le champ aérodynamique $\mathbf{q}$, quant à lui, ne peut être obtenu par résolution des équations de Navier-Stokes, car alors nous aurions simplement décalé le problème du voisinage. Il est lui aussi calculé à partir des relations géométriques et des variables conservatives internes ou frontières.

Cette technique n'est pas sans inconvénients, qui sont cependant exclusivement liés à son implémentation informatique. Tout d'abord, un coût mémoire supplémentaire en est le prix. En effet, soit $N_i$ le nombre de noeuds dans la direction $\xi^i$. Il y a $N_j \times N_k$ noeuds dans un plan à $\xi^i$ constant. Nous ajoutons deux tels plans sur
chacune des deux faces dans la direction $\xi^i$, soit $4N_jN_k$ noeuds au total. Ceci est à répéter dans les directions $\xi^j$ et $\xi^k$. $N_1N_2N_3$ étant le nombre total de noeuds du domaine de calcul non étendu, nous avons un accroissement mémoire relatif de:

$$ \frac{4N_2N_3 + 4N_3N_1 + 4N_1N_2}{N_1N_2N_3} = 4 \left[ \frac{1}{N_1} + \frac{1}{N_2} + \frac{1}{N_3} \right] \quad (4.1) $$

Il varie de quelques 1% pour un gros maillage (≈ 100 000 points) à quelques 10% pour un petit maillage (≈ 1 000 points). D'autre part, nous avons constaté un surcoût en temps de calcul de quelques 5% du coût total pour un maillage de 10 000 points. Ce pourcentage diminue lorsque la taille du maillage croît car la gestion des noeuds supplémentaires est peu coûteuse par rapport aux noeuds du coeur, en lesquels sont chèrement évalués les flux convectifs et diffusifs.

Par contre, de nombreux avantages résultent de cette approche. Premièrement, tous les noeuds du domaine de calcul ont alors suffisamment de voisins, y compris ceux des surfaces limites. Le schéma intérieur peut ainsi être utilisé partout, de manière totallement transparente. Aussi les noeuds d'une frontière géométrique seront-ils gérés dorénavant comme les noeuds internes. Deuxièmement, des traitements assez complexes peuvent être envisagés du fait de ce voisinage étendu, ce qui est un point important lors de la conception d'un code ouvert et versatile. Enfin, cette technique permet une gestion simple et naturelle de conditions traditionnellement difficiles à implémenter, en particulier dans une approche en temps implicite, telle que la périodicité.

### 4.2.2 Différentes conditions géométriques

Par la suite sont présentées quelques configurations intéressantes. Pour chacune seront données l'expression de la position des noeuds supplémentaires, de la base covariante et du champ aérodynamique.

#### Symétrie plane

* **Définition**

Soient $\vec{n}$ un vecteur unitaire perpendiculaire au plan de symétrie, et $P$, un point quelconque de la frontière (figure 4.5). Les noeuds supplémentaires $M_a$ et $M_b$ sont les points symétriques de, respectivement, $M_c$ et $M_d$ par rapport au plan défini par $\vec{n}$ et $P$. Il est à noter qu'un noeud supplémentaire ne dépend que d'un noeud du coeur du domaine, à savoir son symétrique.

* **Exemples d'utilisation**

Cette configuration se rencontre en particulier lors de l'étude de tuyères (figure 4.6) ou de grilles d'aubes présentant un ou deux plans de symétrie, ou dans le cas d'un profil symétrique à zéro degré d'incidence. Cette condition permet d'autre part de simuler une frontière glissante plane pour un fluide parfait (condition miroir).
<figure><img src="image_41.png" /><figcaption>Figure 4.5: Symétrie plane</figcaption></figure>

<figure><img src="image_42.png" /><figcaption>Figure 4.6: Tuyère symétrique</figcaption></figure>
• Relations géométriques

Nous avons:

- pour un scalaire (densité, énergie totale)

$$A^{\text{ext}} = A^{\text{int}} \qquad (4.2)$$

- pour un vecteur (position, quantité de mouvement)

$\vec{A}_t$ et $\vec{A}_n$ étant les composantes tangentielle et normale de $\vec{A}$,

$$\vec{A}_{t}^{\text{ext}} = +\vec{A}_{t}^{\text{int}} \qquad (4.3)$$

$$\vec{A}_{n}^{\text{ext}} = -\vec{A}_{n}^{\text{int}} \qquad (4.4)$$

ou encore,

$$\vec{A}^{\text{ext}} = \vec{A}^{\text{int}} - 2 (\vec{A}^{\text{int}} \cdot \vec{n}) \vec{n} \qquad (4.5)$$

• Position OM

En appliquant la relation 4.5,

$$PM^{\text{ext}} = PM^{\text{int}} - 2 (PM^{\text{int}} \cdot \vec{n}) \vec{n} \qquad (4.6)$$

d'où:

$$OM^{\text{ext}} = OM^{\text{int}} - 2 (OM^{\text{int}} \cdot \vec{n}) \vec{n} \qquad (4.7)$$

• Base covariante $\vec{a}_i$

D'après l'équation 4.7,

$$OM^{\text{ext}} = OM^{\text{int}} - 2 (OM^{\text{int}} \cdot \vec{n}) \vec{n} + 2 (OP \cdot \vec{n}) \vec{n} \qquad (4.8)$$

Il en résulte en différentiant par rapport à $\xi^i$, que:

$$\vec{a}_i^{\text{ext}} = \pm [\vec{a}_i^{\text{int}} - 2 (\vec{a}_i^{\text{int}} \cdot \vec{n}) \vec{n}] \qquad (4.9)$$

Le choix du signe dépend selon que la direction $\xi^i$ est respectivement tangente, ou sécante, à la frontière considérée.
<figure><img src="image_43.png" /><figcaption>Figure 4.7: Périodicité</figcaption></figure>

• Champ aérodynamique $q$

Il est possible de déduire des relations géométriques que:

$$ q^{\text{ext}} = \begin{bmatrix} 1 & 0 & 0 & 0 & 0 \\ 0 & 1 - 2n_1^2 & -2n_1n_2 & -2n_1n_3 & 0 \\ 0 & -2n_2n_1 & 1 - 2n_2^2 & -2n_2n_3 & 0 \\ 0 & -2n_3n_1 & -2n_3n_2 & 1 - 2n_3^2 & 0 \\ 0 & 0 & 0 & 0 & 1 \end{bmatrix} q^{\text{int}} \quad (4.10) $$

## Périodicité

• Définition

Soit $\bar{u}$ un vecteur (figure 4.7). Les noeuds supplémentaires $M_a$ et $M_b$ sont les points images de, respectivement, $M_g$ et $M_h$ par la translation dans le repère cartésien de vecteur $\bar{u}$, appelé le pas. Il est à noter qu'un nocud supplémentaire ne dépend que d'un noeud du coeur du domaine, à savoir son homologue. Ils sont tous deux équivalents, si ce n'est qu'ils sont distincts géométriquement par translation.

• Exemples d'utilisation

Cette condition permet de simuler une grille d'aubes plane (figure 4.8), ou toute autre rangée infinie d'obstacles périodiques. Elle suppose cependant, dans sa définition première présentée ici, que tous les canaux inter-obstacles sont identiques, sans déphasage temporel ou autre distorsion entre eux.
<figure><img src="image_44.png" /><figcaption>Figure 4.8: Grille d'aubes</figcaption></figure>

## • Relations géométriques

Nous avons:

- pour un scalaire (densité, énergie totale)

$$A^{\text{ext}} = A^{\text{int}} \qquad (4.11)$$

- pour un vecteur (quantité de mouvement)

$$\vec{A}^{\text{ext}} = \vec{A}^{\text{int}} \qquad (4.12)$$

## • Position OM

D'après la définition,

$$OM^{\text{ext}} = OM^{\text{int}} + \vec{u} \qquad (4.13)$$

## • Base covariante $\vec{a}_i$

En différentiant l'équation 4.13 par rapport à $\xi^i$, nous obtenons:

$$\vec{a}_i^{\text{ext}} = \vec{a}_i^{\text{int}} \qquad (4.14)$$
<figure><img src="image_45.png" /><figcaption>Figure 4.9: Périodicité angulaire</figcaption></figure>

* Champ aérodynamique $q$

Il est possible de déduire des relations géométriques que:

$$ q^{\text{ext}} = \begin{bmatrix} 1 & 0 & 0 & 0 & 0 \\ 0 & 1 & 0 & 0 & 0 \\ 0 & 0 & 1 & 0 & 0 \\ 0 & 0 & 0 & 1 & 0 \\ 0 & 0 & 0 & 0 & 1 \end{bmatrix} q^{\text{int}} \quad (4.15) $$

Cette condition géométrique est la plus simple de toutes. Son implémentation est remarquablement facile, grâce à l'adjonction des noeuds supplémentaires. Une contrainte y est cependant attachée. Elle suppose en effet qu'il y a colloca- tion des noeuds frontières, au pas près, sur deux faces homologues du domaine de calcul.

## Périodicité angulaire

* Définition

Soit $\Delta\theta$ un réel (figure 4.9). Les noeuds supplémentaires $M_a$ et $M_b$ sont les points images de, respectivement, $M_g$ et $M_h$ par la rotation dans le repère cartésien d'axe $Ox_1$ et d'angle $\Delta\theta$, appelé le *pas*. Il est à noter qu'un noeud supplémentaire ne dépend que d'un noeud du coeur du domaine, à savoir son homologue. Ils sont tous deux équivalents, si ce n'est qu'ils sont distincts
<figure><img src="image_46.png" /><figcaption>Figure 4.10: Inducteur transsonique</figcaption></figure>

géométriquement par rotation. Cette transformation correspond, dans le repère cylindrique, à une translation dans la direction angulaire.

* Exemples d'utilisation

Cette configuration est primordiale pour simuler une roue de turbomachine (figure 4.10). Elle permet en effet de ne calculer qu'un seul canal inter-aubes, dans la mesure où elle suppose comme précédement, que tous les canaux sont identiques, sans déphasage temporel ou autre distorsion entre eux.

* Relations géométriques

Nous avons:

- pour un scalaire (densité, énergie totale)

$$A^{ext} = A^{int} \qquad (4.16)$$

- pour un vecteur (position, quantité de mouvement)

$$ \left\{ \begin{aligned} A_1^{ext} &= A_1^{int} \\ A_2^{ext} &= A_2^{int} \cos\Delta\theta - A_3^{int} \sin\Delta\theta \\ A_3^{ext} &= A_3^{int} \sin\Delta\theta + A_3^{int} \cos\Delta\theta \end{aligned} \right. \qquad (4.17) $$

* Position OM

Par application directe de la relation 4.17, nous obtenons:
<figure><img src="image_47.png" /><figcaption>Figure 4.11: Extrapolation</figcaption></figure>

$$ egin{cases} x_1^{ext} = x_1^{int} \\ x_2^{ext} = x_2^{int} \cos\Delta\theta - x_3^{int} \sin\Delta\theta \\ x_3^{ext} = x_2^{int} \sin\Delta\theta + x_3^{int} \cos\Delta\theta \end{cases} \quad (4.18) $$

* Base covariante $\tilde{a}_i$

La différentiation de l'équation 4.18 par rapport à $\xi^i$ conduit à:

$$ egin{cases} a_{i1}^{ext} = a_{i1}^{int} \\ a_{i2}^{ext} = a_{i2}^{int} \cos\Delta\theta - a_{i3}^{int} \sin\Delta\theta \\ a_{i3}^{ext} = a_{i2}^{int} \sin\Delta\theta + a_{i3}^{int} \cos\Delta\theta \end{cases} \quad (4.19) $$

* Champ aérodynamique $q$

Il est possible de déduire des relations géométriques que:

$$ q^{ext} = \begin{bmatrix} 1 & 0 & 0 & 0 & 0 \\ 0 & 1 & 0 & 0 & 0 \\ 0 & 0 & \cos\Delta\theta & -\sin\Delta\theta & 0 \\ 0 & 0 & \sin\Delta\theta & \cos\Delta\theta & 0 \\ 0 & 0 & 0 & 0 & 1 \end{bmatrix} q^{int} \quad (4.20) $$

## Extrapolation

* Définition

Les noeuds supplémentaires $M_a$ et $M_b$ sont tels que le noeud $M_c$, qui est sur la frontière, soit le milieu du segment de droite les reliant à leur point source, à savoir $M_e$ et $M_d$ (figure 4.11). Il est à noter qu'un noeud supplémentaire dépend de deux noeuds, un du coeur et un de la bordure du domaine.
* Exemples d'utilisation

Cette condition est singulière dans la mesure où elle ne correspond à aucune configuration géométrique particulière. Cependant, elle permet de gérer les noeuds supplémentaires en lesquels il n'y a pas d'autre propriété géométrique remarquable, et qui seraient donc impossible à obtenir autrement.

* Relations géométriques

Nous avons:

- pour un scalaire (densité, énergie totale)

$$A^{\text{ext}} = 2A^c - A^{\text{int}} \quad (4.21)$$

- pour un vecteur (position, quantité de mouvement)

$$\vec{A}^{\text{ext}} = 2\vec{A}^c - \vec{A}^{\text{int}} \quad (4.22)$$

* Position OM

D'après cette dernière relation,

$$OM^{\text{ext}} = 2OM^c - OM^{\text{int}} \quad (4.23)$$

* Base covariante $\vec{a}_i$

En différentiant cette équation par rapport à $\xi^i$, nous obtenons:

$$\vec{a}_i^{\text{ext}} = 2\vec{a}_i^c - \vec{a}_i^{\text{int}} \quad (4.24)$$

* Champ aérodynamique $q$

Il est possible de déduire des relations géométriques que:

$$q^{\text{ext}} = \begin{bmatrix} 2 & 0 & 0 & 0 & 0 \\ 0 & 2 & 0 & 0 & 0 \\ 0 & 0 & 2 & 0 & 0 \\ 0 & 0 & 0 & 2 & 0 \\ 0 & 0 & 0 & 0 & 2 \end{bmatrix} q^c - \begin{bmatrix} 1 & 0 & 0 & 0 & 0 \\ 0 & 1 & 0 & 0 & 0 \\ 0 & 0 & 1 & 0 & 0 \\ 0 & 0 & 0 & 1 & 0 \\ 0 & 0 & 0 & 0 & 1 \end{bmatrix} q^{\text{int}} \quad (4.25)$$

Il est à noter que nous envisageons ici une extrapolation linéaire, qui ramène donc le schéma spatial au premier ordre au voisinage de la frontière. D'autre part, cette opération peut être apparentée à un décentrement forcé depuis l'intérieur du domaine, sans que le sens et le régime de l'écoulement soient pris en compte, comme ceci est le cas dans la séparation de flux. Aussi, une perte de précision est à escompter si les informations extrapolées sur les noeuds supplémentaires sont utilisées.
Remarques complémentaires

* Calculs d'écoulements mono ou bidimensionnels:

Trois noeuds uniquement, soient deux cellules, sont nécessaires dans les directions non significatives de tels écoulements. Une condition de symétrie plane, ou de périodicité, suffit alors à simuler une géométrie infiniment étendue dans une ou deux directions. En pratique, la première solution se révèle plus apte à absorber d'éventuelles fluctuations parasites d'origine numérique.

* Périodicité déphasée. Périodicité instationnaire

La définition des conditions de périodicité donnée précédemment suppose que tous les canaux sont équivalents aussi bien au niveau géométrique que sur le plan aérodynamique. Or cette hypothèse n'est pas toujours vérifiée, et un déphasage peut exister entre deux, ou plusieurs canaux. C'est par exemple le cas pour un étage où le stator et le rotor n'ont pas le même nombre de pales. La longueur du motif de périodicité est alors égale au plus petit multiple commun des deux pas. Le flottement et le décrochage tournant, où une ou plusieurs distorsions se déplacent de manière asynchrone devant une roue, sont d'autres configurations où la périodicité telle qu'elle a été présentée jusqu'ici est inutilisable.

La résolution exacte de ces problèmes passe par le calcul simultané d'autant de canaux que nécessaires pour retrouver une périodicité. Dans des cas extrèmes, il faudrait même calculer la roue complète. Or cette solution est très coûteuse en temps de calcul et en taille mémoire. Elle est encore raisonnablement inutilisable aujourd'hui, car il est impossible d'obtenir une bonne résolution de tous les canaux avec les moyens de calcul actuels.

Il existe cependant une solution approchée pour l'étude du flottement. Elle suppose un déphasage linéairement réparti entre tous les canaux du motif de périodicité. Cette distorsion est imposée sur les noeuds supplémentaires, qui sont alors équivalents à leur homologue à une translation spatiale et temporelle près. Le calcul peut ainsi être restreint à un seul canal. Il est à noter que l'extension du domaine de calcul facilite une fois de plus grandement l'implémentation de cette technique. En effet, la structure retenue permet la gestion de l'instationnarité de façon naturelle: elle résulte bien du voisinage de la frontière considérée.

### 4.2.3 Gestion de maillages en O ou en C

Il est à remarquer que si le pas de la périodicité angulaire est tel que $\Delta\theta = 2\pi$, le domaine de calcul se referme sur lui-même (figure 4.12). Nous obtenons alors un maillage dit en O, autour d'une surface de révolution qui correspond au moyeu en turbomachine. A partir de cette constatation, il apparaît possible de réaliser un solveur en O ou en C, à partir de celui en H que nous avons développé jusqu'ici, sans
<figure><img src="image_48.png" /><figcaption>Figure 4.12: Construction d'un maillage en O</figcaption></figure>

modifier en aucune façon le traitement du coeur et des frontières. En effet, il suffit de gérer de manière adaptée les noeuds supplémentaires.

Cette extension de la méthode offre une plus grande flexibilité dans la discrétisation du domaine de calcul qui n'est pas toujours aisée dans une approche structurée. Ainsi un type de maillage mieux adapté à la géométrie à traiter pourra être choisi, suivant qu'un canal, qu'un profil isolé, qu'un bord d'attaque ou de fuite rond ou effilé, seront étudiés.

### Gestion de maillages en O

Dans ce cas, il suffit, comme vu précédemment, que le domaine en H se referme parfaitement sur deux de ses faces opposées (figure 4.12). Une périodicité de vecteur nul, ou une périodicité angulaire à pas nul (figure 4.13), permettent alors de simuler cette configuration sans qu'aucune modification du solveur en H ne soit nécessaire. A titre d'illustration, est présenté figure 4.14 les iso-lignes de nombre de Mach autour d'un cylindre infini (figure 4.15), plongé dans un fluide parfait à basse vitesse, le nombre de Mach amont étant de 0.3. Il est à noter que le champ aérodynamique n'est en rien perturbé au niveau du raccord, qui apparaît dans le premier quadrant sous la forme d'une droite figure 4.14.

### Gestion de maillages en C

Il est possible de refermer le maillage en H en "soudant" non plus deux faces opposées, mais la même face sur elle-même (figure 4.16). Nous obtenons alors un domaine dit en C. La collocation des noeuds $M_a$ et $M_g$ (resp. $M_b$ et $M_h$) assure la continuité du domaine de calcul à travers le raccord. Cette configuration n'a pas été testée, mais elle apparaît comme parfaitement réalisable avec
<figure><img src="image_49.png" /><figcaption>Figure 4.13: Maillage en O</figcaption></figure>

<figure><img src="image_50.png" /><figcaption>Figure 4.14: Ecoulement autour d'un cylindre. Nombre de Mach</figcaption></figure>
<figure><img src="image_51.png" /><figcaption>Figure 4.15: Ecoulement autour d'un cylindre. Maillage</figcaption></figure>

<figure><img src="image_52.png" /><figcaption>Figure 4.16: Maillage en C</figcaption></figure>
<figure><img src="image_53.png" /><figcaption>Figure 4.17: Construction d'un multi-domaines</figcaption></figure>

de fortes chances de succès.

La possibilité d'étendre à très faible coût la méthode initiallement développée sur un maillage en H à des maillages en O ou en C est un autre argument favorable à l'extension du domaine de calcul. La gestion du type de géométrie est entièrement transparente pour le coeur et les frontières, car elle affecte exclusivement les noeuds supplémentaires.

### 4.2.4 Gestion d'un multi-domaines

Nous avons vu que la gestion des maillages en O ou en C se faisaient par recollement d'une frontière sur une autre frontière, pouvant être la même, du même domaine de calcul, les noeuds supplémentaires assurant la jonction entre les deux étant construits par collocation à partir des noeuds du coeur. Cette idée peut être généralisée, si ce n'est que les données nécessaires à l'extension du maillage considéré proviennent non plus de lui-même, mais d'un autre domaine (figure 4.17).

Il est alors possible de découper des géométries très complexes, difficiles à mailler en un seul bloc structuré, en plusieurs sous-domaines élémentaires. Ainsi, un canal inter-aubes en turbomachine pourra-t-il être décomposé en deux domaines en H à l'amont et à l'aval de la roue, en serrant un domaine en O ou en C autour de l'aube, un dernier domaine en H permettant de prendre en compte un éventuel jeu en bout de pale. Un autre exemple met en oeuvre deux domaines en H collés en emprisonnant un obstacle. Cette configuration permet de simuler un profil, avec ou sans incidence, dans une tuyère ou encore une pale intercalaire entre deux aubes.

La technique utilisée repose sur le fait que le domaine adjacent à celui considéré est vu à travers les noeuds supplémentaires. Il y a recollement bord à bord
<figure><img src="image_54.png" /><figcaption>Figure 4.18: Raccord entre deux domaines</figcaption></figure>

des domaines, avec recouvrement sur cinq noeuds collocatifs (figure 4.18). Une bonne continuité d'un domaine à l'autre est ainsi assurée, sans perte de précision ou d'introduction d'erreur par interpolation dans la zone de chevauchement. Les noeuds frontières voient effectivement les deux domaines simultanément.

Là encore, la simplicité d'un tel traitement est due à l'extension du domaine de calcul initial, qui se justifie donc une fois de plus.

## 4.3 Conditions physiques fortes

Ce traitement s'applique aux noeuds d'une frontière du domaine de calcul en lesquels une contrainte extérieure est imposée directement dans le système différentiel y décrivant l'évolution du champ aérodynamique. Cette modification des équations de conservation repose sur le sens de propagation des informations traversant la surface limite considérée, en exploitant la même idée que celle mise en oeuvre lors de la séparation des flux convectifs. Pour cela, une nouvelle approche du système initial est suivie afin de faire apparaître des dérivées tangentes et sécantes à la frontière. Les premières seront traitées comme dans le coeur du domaine, alors que les secondes feront l'objet d'un traitement particulier, en fonction de la nature entrante ou sorrante des ondes ainsi misces en évidence [27, pp. 362-364]. Le nombre de ces dernières nous renseigne en particulier sur le nombre de conditions aux limites physiques à imposer.

Il est à noter que la méthode présentée ici est, comme beaucoup d'autres, bâtie à partir des équations d'Euler, de nature hyperbolique. Dans le cas des équations de Navier-Stokes moyennées, qui sont incomplètement paraboliques, les termes visqueux
et turbulents sont traités comme des termes sources. D'autre part, contrairement aux conditions géométriques, cette technique concerne directement un noeud de la frontière, et non plus uniquement son voisinage étendu.

### 4.3.1 Equations d'évolution en variables caractéristiques

La première étape consiste en une réécriture des équations d'évolution sous une forme propice à des manipulations inspirées par le sens physique. Pour cela, nous exploitons le fait qu'une surface caractéristique de normale $\vec{n}$, a pour particularité qu'une combinaison linéaire appropriée des équations d'évolution conduit à un système ne contenant que des différentielles dans des directions tangentes à cette surface [27, pp. 152]. Ceci est réalisé lorsque les vecteurs propres gauches, associés à la matrice jacobienne du flux convectif $\vec{F}_c \cdot \vec{n}$, sont utilisés comme coefficients de cette combinaison linéaire. Il est possible alors d'interpréter le nouveau système d'équations comme décrivant des ondes se déplaçant le long de la normale. Leur sens de propagation est la clef du traitement numérique proposé.

Dans notre cas, une frontière correspond à une surface sur laquelle une des coordonnées curviligne, notée $\xi^i$, est constante. Si bien que $\pm \vec{a}^i$ est une normale à cette surface, pointant vers l'extérieur du domaine de calcul. Le gradient de $\xi^i$ étant orienté dans le sens des valeurs croissantes, le signe -- est choisi pour les frontières d'indice minimal, alors qu'aux indices maximaux est associé le signe +. Nous nous limiterons à ce dernier cas, l'autre étant traiter par remplacement de $\vec{a}^i$ par $-\vec{a}^i$. Dans ces conditions, $\vec{F}_c \cdot \vec{n}$, le flux convectif traversant la surface, est égal à $F_c^i$.

Par la suite sont présentés la matrice jacobienne qui est associée à ce dernier, les valeurs propres et vecteurs propres de celle-ci, ainsi que le nouveau système d'équations obtenu.

#### Matrice jacobienne des flux convectifs

La matrice jacobienne $K^i$ associée au flux convectif $F_c^i$ est ainsi définie:

$$K^i = [K_{ab}^i]_{(a,b)\in[1;5]^2} = [K_a^i]_{a\in[1;5]} \quad (4.26)$$

avec:

$$K_{ab}^i = \frac{\partial F_{ca}^i}{\partial (\sqrt{g_{bb}})} \quad (4.27)$$

$$K_a^i = \frac{\partial F_{ca}^i}{\partial (\sqrt{g_{qq}})} \quad (4.28)$$

Soient:
$$ K_1^i = \begin{vmatrix} \xi_{x_1}^i & & & \\ & \xi_{x_2}^i & & \\ & & \xi_{x_3}^i & \\ & & & 0 \end{vmatrix} \tag{4.29} $$

$$ K_2^i = \begin{vmatrix} -u_1(V^i - \xi_i^i) + \xi_{x_1}^i \Phi & & & \\ V^i - (\gamma - 2)\xi_{x_1}^i u_1 & u_1\xi_{x_2}^i - (\gamma - 1)\xi_{x_1}^i u_2 & & \\ u_1\xi_{x_3}^i - (\gamma - 1)\xi_{x_1}^i u_3 & u_2\xi_{x_2}^i - (\gamma - 1)\xi_{x_2}^i u_1 & & \\ (\gamma - 1)\xi_{x_1}^i & (\gamma - 1)\xi_{x_2}^i & & \end{vmatrix} \tag{4.30} $$

$$ K_3^i = \begin{vmatrix} -u_2(V^i - \xi_i^i) + \xi_{x_2}^i \Phi & & & \\ u_2\xi_{x_1}^i - (\gamma - 1)\xi_{x_2}^i u_1 & u_2\xi_{x_3}^i - (\gamma - 1)\xi_{x_2}^i u_3 & & \\ V^i - (\gamma - 2)\xi_{x_2}^i u_2 & (\gamma - 1)\xi_{x_2}^i & & \\ (\gamma - 1)\xi_{x_2}^i & (\gamma - 1)\xi_{x_3}^i & & \end{vmatrix} \tag{4.31} $$

$$ K_4^i = \begin{vmatrix} -u_3(V^i - \xi_i^i) + \xi_{x_3}^i \Phi & & & \\ u_3\xi_{x_1}^i - (\gamma - 1)\xi_{x_3}^i u_1 & u_3\xi_{x_2}^i - (\gamma - 1)\xi_{x_3}^i u_2 & & \\ V^i - (\gamma - 2)\xi_{x_3}^i u_3 & (\gamma - 1)\xi_{x_3}^i & & \\ (\gamma - 1)\xi_{x_3}^i & (\gamma - 1)\xi_{x_3}^i & & \end{vmatrix} \tag{4.32} $$

$$ K_5^i = \begin{vmatrix} (2\Phi - \gamma E)(V^i - \xi_i^i) & & & \\ (\gamma E - \Phi)\xi_{x_1}^i - (\gamma - 1)(V^i - \xi_i^i)u_1 & u_1\xi_{x_2}^i - (\gamma - 1)(V^i - \xi_i^i)u_2 & & \\ (\gamma E - \Phi)\xi_{x_2}^i - (\gamma - 1)(V^i - \xi_i^i)u_3 & (\gamma E - \Phi)\xi_{x_3}^i - (\gamma - 1)(V^i - \xi_i^i)u_3 & & \\ \xi_i^i + \gamma(V^i - \xi_i^i) & & & \end{vmatrix} \tag{4.33} $$

où:

$$ \Phi = \frac{1}{2}(\gamma - 1)\bar{V}^2 \tag{4.34} $$

Il est à noter que $V^i - \xi_i^i = \xi_{x_1}^i u_1 + \xi_{x_2}^i u_2 + \xi_{x_3}^i u_3$, et que donc, seuls deux termes supplémentaires sont introduits par la déformation en temps du maillage $\xi_i^i$, dans la première composante de $K_i^i$ et la dernière de $K_i^g$.

## Valeurs propres

Une valeur propre $\lambda$ de $K^i$ vérifie par définition:

$$ \det |K^i - \lambda I_5| = 0 \tag{4.35} $$

La résolution de cette équation fournit cinq valeurs réelles de $\lambda$, ce qui caractérise un système hyperbolique. Ce sont $\lambda = V^i$, d'ordre 3, et $\lambda = V^i \pm a \parallel \bar{a}^i \parallel$, toutes deux d'ordre 1.
## Vecteurs propres gauches et droits

Un vecteur propre gauche $l^i$, écrit en ligne, et un vecteur propre droit $r^i$, écrit en colonne, associés à la matrice jacobienne $K^i$ et à la valeur propre $\lambda$ sont définis par:

$$l^i K^i = \lambda l^i \quad (4.36)$$

$$K^i r^i = \lambda r^i \quad (4.37)$$

Soient $L^i$ la matrice formée des vecteurs propres gauches $l^i$ écrits en ligne, $R^i$ la matrice des vecteurs propres droits $r^i$ écrits en colonnes, et $\Lambda$ une matrice dont les valeurs propres constituent la diagonale, et qui est remplie de zéros par ailleurs. Elles vérifient les propriétés suivantes:

$$L^i R^i = I_5 \quad (4.38)$$

$$R^i L^i = I_5 \quad (4.39)$$

$$L^i K^i = \Lambda L^i \quad (4.40)$$

$$K^i R^i = R^i \Lambda \quad (4.41)$$

$$L^i K^i R^i = \Lambda \quad (4.42)$$

Pour chaque valeur propre, nous avons:

$$\lambda = V^i$$

$$L_1^i = \begin{bmatrix} [1 - \alpha] \hat{\xi}_{x_1}^i - \frac{1}{a} [\hat{\xi}_{x_3}^i u_2 - \hat{\xi}_{x_2}^i u_3] \\ \beta \hat{\xi}_{x_1}^i u_1 \\ \beta \hat{\xi}_{x_2}^i u_2 + \frac{1}{a} \hat{\xi}_{x_3}^i \\ \beta \hat{\xi}_{x_1}^i u_3 - \frac{1}{a} \hat{\xi}_{x_2}^i \\ -\beta \hat{\xi}_{x_1}^i \end{bmatrix} \quad R_1^i = \begin{bmatrix} \hat{\xi}_{x_1}^i \\ \hat{\xi}_{x_1}^i u_1 \\ \hat{\xi}_{x_1}^i u_2 + a \hat{\xi}_{x_3}^i \\ \hat{\xi}_{x_1}^i u_3 - a \hat{\xi}_{x_2}^i \\ \frac{1}{2} \hat{V}^2 \hat{\xi}_{x_1}^i + a [\hat{\xi}_{x_3}^i u_2 - \hat{\xi}_{x_2}^i u_3] \end{bmatrix} \quad (4.43) \quad (4.44)$$

$$L_2^i = \begin{bmatrix} [1 - \alpha] \hat{\xi}_{x_2}^i - \frac{1}{a} [\hat{\xi}_{x_1}^i u_3 - \hat{\xi}_{x_3}^i u_1] \\ \beta \hat{\xi}_{x_2}^i u_1 - \frac{1}{a} \hat{\xi}_{x_3}^i \\ \beta \hat{\xi}_{x_2}^i u_2 \\ \beta \hat{\xi}_{x_2}^i u_3 + \frac{1}{a} \hat{\xi}_{x_1}^i \\ -\beta \hat{\xi}_{x_2}^i \end{bmatrix} \quad R_2^i = \begin{bmatrix} \hat{\xi}_{x_2}^i \\ \hat{\xi}_{x_2}^i u_1 - a \hat{\xi}_{x_3}^i \\ \hat{\xi}_{x_2}^i u_2 \\ \hat{\xi}_{x_2}^i u_3 + a \hat{\xi}_{x_1}^i \\ \frac{1}{2} \hat{V}^2 \hat{\xi}_{x_2}^i + a [\hat{\xi}_{x_1}^i u_3 - \hat{\xi}_{x_3}^i u_1] \end{bmatrix} \quad (4.45) \quad (4.46)$$
$$
L_3^i = \begin{bmatrix} [1-\alpha]\hat{\xi}_{x_3}^i - \frac{1}{a}[\hat{\xi}_{x_2}^i u_1 - \hat{\xi}_{x_1}^i u_2] \\ \beta\hat{\xi}_{x_3}^i u_1 + \frac{1}{a}\hat{\xi}_{x_2}^i \\ \beta\hat{\xi}_{x_3}^i u_2 - \frac{1}{a}\hat{\xi}_{x_1}^i \\ \beta\hat{\xi}_{x_3}^i u_3 \\ -\beta\hat{\xi}_{x_3}^i \end{bmatrix} \qquad R_3^i = \begin{bmatrix} \hat{\xi}_{x_3}^i \\ \hat{\xi}_{x_3}^i u_1 + a\hat{\xi}_{x_2}^i \\ \hat{\xi}_{x_3}^i u_2 - a\hat{\xi}_{x_1}^i \\ \hat{\xi}_{x_3}^i u_3 \\ \frac{1}{2}\vec{V}^2\hat{\xi}_{x_3}^i + a[\hat{\xi}_{x_2}^i u_1 - \hat{\xi}_{x_1}^i u_2] \end{bmatrix} \qquad (4.47) \qquad (4.48)
$$

$$
\lambda = V^i + a \| \vec{a}^i \|
$$

$$
L_4^i = \begin{vmatrix} \alpha - \frac{1}{a \|\vec{a}^i\|} (V^i - \xi_i^i) \\ -\beta u_1 + \frac{1}{a} \hat{\xi}_{x_1}^i \\ -\beta u_2 + \frac{1}{a} \hat{\xi}_{x_2}^i \\ -\beta u_3 + \frac{1}{a} \hat{\xi}_{x_3}^i \\ \beta \end{vmatrix} \quad (4.49) \quad R_4^i = \frac{1}{2} \begin{vmatrix} 1 & & & \\ u_1 + a \hat{\xi}_{x_1}^i & u_2 + a \hat{\xi}_{x_2}^i & u_3 + a \hat{\xi}_{x_3}^i & \\ \frac{1}{2} \vec{V}^2 + \frac{1}{\beta} + \frac{a}{\|\vec{a}^i\|} (V^i - \xi_i^i) & & & \\ & & & \end{vmatrix} (4.50)
$$

$$
\lambda = V^i - a \| \vec{a}^i \|
$$

$$
L_5^i = \begin{vmatrix} \alpha + \frac{1}{a \|\vec{a}^i\|} (V^i - \xi_i^i) \\ -\beta u_1 - \frac{1}{a} \hat{\xi}_{x_1}^i \\ -\beta u_2 - \frac{1}{a} \hat{\xi}_{x_2}^i \\ -\beta u_3 - \frac{1}{a} \hat{\xi}_{x_3}^i \\ \beta \end{vmatrix} \quad (4.51) \quad R_5^i = \frac{1}{2} \begin{vmatrix} 1 & & & \\ u_1 - a \hat{\xi}_{x_1}^i & u_2 - a \hat{\xi}_{x_2}^i & u_3 - a \hat{\xi}_{x_3}^i & \\ \frac{1}{2} \vec{V}^2 + \frac{1}{\beta} - \frac{a}{\|\vec{a}^i\|} (V^i - \xi_i^i) & & & \\ & & & \end{vmatrix} \quad (4.52)
$$

où:

$$
\alpha = \frac{1(\gamma - 1)\vec{V}^2}{2 a^2} \qquad (4.53)
$$

$$
\beta = \frac{(\gamma - 1)}{a^2} \qquad (4.54)
$$

$$
\hat{\xi}_{x_j}^i = \frac{1}{\|\vec{a}^i\|} \xi_{x_j}^i \qquad (4.55)
$$

Nous remarquons que les valeurs propres sont modifiées pour prendre en compte
la vitesse de déplacement des noeuds du maillage, car $V^i = (\vec{V} - \bar{\vec{V}}) \cdot \vec{a}^i$, mais non les
vecteurs propres. D'autre part, le choix mathématique de ces derniers est arbitraire, à
la seule condition qu'ils forment une base du sous-espace propre associé. Ceci impose
qu'ils ne doivent jamais dégénérer, quelle que soit la configuration géométrique et
aérodynamique rencontrée. C'est en particulier le cas pour le sous-espace propre
de dimension trois, qui doit en plus satisfaire à une certaine symétrie afin de ne
privilégier aucune direction. Aussi, les bases ($L_1^i, L_2^i, L_3^i$) et ($R_1^i, R_2^i, R_3^i$) choisies sont-
elles identiques par permutation des indices. Enfin, d'un point de vue physique, les
vecteurs propres doivent être de dimensions homogènes, c'est-à-dire que les termes
d'une colonne de $L^i$ (resp. d'une ligne de $R^i$) ont la même unité physique, et si
possible, un ordre de grandeur semblable.

Les matrices $L^i$ et $R^i$ présentées précédemment satisfont à tous ces critères.

Formulation caractéristique mono-directionnelle

Rappelons que les vecteurs propres gauches de $K^i$ sont les coefficients de la combi-
naison linéaire qui permet de projeter les équations de conservation sur la surface
caractéristique. Le nouveau système ainsi obtenu forme les équations de compatibilité
associées à la normale $\bar{\epsilon}^i$. Voici comment est réalisée cette transformation.

Soit le système différentiel 3.31 associé aux équations de conservation en un point
de la frontière du domaine:

$$
\frac{\partial}{\partial t}(\sqrt{g}q) + \sum_{j=1}^{3} \frac{\partial F_c^j}{\partial \xi^j} + S = 0 \quad (4.56)
$$

$S$ contenant les termes sources, y compris les flux visqueux et turbulents.

Multiplié par la matrice des vecteurs propres gauches de $K^i$, il s'écrit alors:

$$
L^i \frac{\partial}{\partial t} (\sqrt{g}q) + L^i \left\{ \frac{\partial F_c^i}{\partial \xi^i} + \left[ \sum_{\substack{j=1 \\ j \neq i}}^3 \frac{\partial F_c^j}{\partial \xi^j} + S \right] \right\} = 0 \quad (4.57)
$$

Puis, l'introduction de la matrice jacobienne conduit à:

$$
L^i \frac{\partial}{\partial t} (\sqrt{g}q) + L^i K^i \frac{\partial}{\partial \xi^i} (\sqrt{g}q) + L^i \left[ \sum_{\substack{j=1 \\ j \neq i}}^3 \frac{\partial F_c^j}{\partial \xi^j} + S \right] = 0 \quad (4.58)
$$

Soient $w$, les variables dites caractéristiques, ainsi définies sous forme incrémentale:

$$
\delta (\sqrt{g}w) = L^i \delta (\sqrt{g}q) \tag{4.59}
$$

Nous introduisons de plus, un nouveau vecteur colonne $\ell^i$ tel que:

$$
\ell^i = L^i \frac{\partial F_c^i}{\partial \xi^i} \qquad (4.60)
$$

$$
= L^i K^i \frac{\partial}{\partial \xi^i} (\sqrt{g} q) \qquad (4.61)
$$

$$
= \Lambda L^i \frac{\partial}{\partial \xi^i} (\sqrt{g} q) \qquad (4.62)
$$

$$
= \Lambda \frac{\partial}{\partial \xi^i} (\sqrt{g} w) \qquad (4.63)
$$
Si bien que l'équation 4.58 se réécrit :

$$ \frac{\partial}{\partial t}(\sqrt{g}w) + \Lambda \frac{\partial}{\partial \xi^i}(\sqrt{g}w) + L^i \left[ \sum_{\substack{j=1 \\ j \neq i}}^{3} \frac{\partial F_c^j}{\partial \xi^j} + S \right] = 0 \quad (4.64) $$

Par contre, si nous la multiplions par la matrice des vecteurs propres droits $R^i$, afin de retrouver le système initial, nous obtenons :

$$ \frac{\partial}{\partial t}(\sqrt{g}q) + R^i \ell^i + \sum_{\substack{j=1 \\ j \neq i}}^{3} \frac{\partial F_c^j}{\partial \xi^j} + S = 0 \quad (4.65) $$

Le nouveau système 4.64 répond à notre besoin. Il fait clairement apparaître des termes tangents à la surface $\frac{\partial F_c^j}{\partial \xi^j}$, et d'autres sécants $\frac{\partial}{\partial \xi^i}(\sqrt{g}w)$. De plus, nous pouvons lui donner un sens physique en terme de propagation d'ondes. En effet, il décrit le transport des grandeurs caractéristiques $w$ dans la direction $\xi^i$, à la vitesse $\Lambda$, transport perturbé par un terme source incluant les variations perpendiculaires au cheminement et les effets dissipatifs. Ainsi, les valeurs propres de la matrice jacobienne $K^i$ peuvent être interprétées comme les vitesses des ondes se déplaçant perpendiculairement à la surface caractéristique, et donc colinéairement à la normale $\bar{a}^i$. Leur signe nous renseigne en particulier sur leur sens de propagation, et permet ainsi d'identifier les informations qui sortent du domaine de calcul, associées à une valeur propre positive, par opposition à celles qui y entrent à une vitesse négative. Il reste à exploiter ces renseignements pour discrétiser les équations à la frontière, et imposer les conditions aux limites.

### 4.3.2 Imposition forte d'une condition aux limites

Nous poursuivons deux buts, dont le premier est d'imposer une contrainte extérieure sur le domaine de calcul à travers ses noeuds frontières, et ce par valeur, par exemple un niveau de pression statique en sortie subsonique d'une roue. Le second est de discrétiser le système d'équations différentielles de façon cohérente sur le plan physique, en prenant en compte le sens de propagation des informations à travers la surface limite, et sur le plan numérique, si nous jugeons le voisinage étendu insufficient, ou en tous les cas peu précis à cause de l'extrapolation grossière qui y est réalisée.

L'idée présidant à la réalisation de ces objectifs repose sur la mise en évidence de deux familles d'ondes traversant la frontière. La première sort du domaine de calcul. Elle transporte une information depuis l'intérieur, qui est parfaitement connu puisque les équations de Navier-Stokes moyennées y sont résolues sur un voisinage satisfaisant. Un décentrement dans le sens de propagation est donc souhaitable. La seconde entre dans le domaine. Elle transporte une information qui traduit l'influence du reste de l'univers au point frontière considéré. Provenant d'une région non discrétisée, elle ne peut être connue avec précision, et est donc inutilisable. Par
contre, la physique, et la pratique, nous permettent de savoir quelle contrainte est imposée. C'est par exemple un niveau d'énergie, ou une direction d'écoulement. Les équations associées aux ondes entrantes sont donc retirées du système, pour être remplacées par des conditions physiques convenablement formulées.

## Mise en oeuvre

Le signe des valeurs propres $\Lambda$ apparaissant dans l'équation 4.64 permet de connaître le sens de propagation des informations à travers la frontière. Si $\lambda$ est positive, l'onde associée sort du domaine, et l'équation qui la décrit est à conserver. Par contre, si $\lambda$ est négative, l'équation concernée est supprimée. Le système 4.57, qui est une autre formulation de l'équation 4.64, devient alors:

$$ [P^{2i}] \frac{\partial}{\partial t} (\sqrt{g}q) + [P^{2i}] \left\{ \frac{\partial F_c^i}{\partial \xi^i} + \left[ \sum_{\substack{j=1 \\ j \neq i}}^{3} \frac{\partial F_c^j}{\partial \xi^j} + S \right] \right\} = 0 \quad (4.66) $$

où $[P^{2i}]$ est une matrice formée des lignes de $L^i$ correspondant aux valeurs positives de $\Lambda$, les autres étant remplacées par le vecteur ligne nul. Nous réalisons ainsi le filtrage du système ondulatoire.

Cet ensemble d'équations n'est plus fermé, puisqu'il contient autant de lignes remplies de zéros, qu'il y a d'ondes entrantes au domaine. Une condition physique doit alors être imposée pour chacune d'entre elles. Afin de conserver le même schéma en temps en tout point, cette dernière doit être exprimée sous forme d'une équation aux dérivées partielles du premier ordre en temps. Soit la relation fixant une certaine grandeur $B(q)$, telle que la pression statique ou la vitesse tangentielle, dont le niveau imposé, potentiellement variable dans le temps, serait $B^*(t)$:

$$ B(q) - B^*(t) = 0 \quad (4.67) $$

Elle est linéarisée au premier ordre en temps de la manière suivante:

$$ B^{(n+1)} - B^{* (n+1)} = \left\{ B^{(n)} + \left[ \frac{\partial B}{\partial(\sqrt{g}q)} \right]^{(n)} \left( \frac{\partial(\sqrt{g}q)}{\partial t} \right)^{(n)} \Delta t + O(\Delta t)^2 \right\} - B^{* (n+1)} \quad (4.68) $$

Si l'équation 4.67 est vérifiée à l'instant $n+1$,

$$ B^{(n+1)} - B^{* (n+1)} \approx 0 \quad (4.69) $$

il en résulte que:

$$ \left[ \frac{\partial B}{\partial(\sqrt{g}q)} \right] \frac{\partial}{\partial t} (\sqrt{g}q) = - \frac{1}{\Delta t} \left\{ B(q^n) - B^*(t^{n+1}) \right\} \quad (4.70) $$