Cette dernière équation est formellement semblable à l'équation 4.66, avec laquelle elle est fusionnée. Si bien qu'il est maintenant possible de remplacer le système 4.57 par:

$$[P^{1i}] \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\} = [P^{3i}] \quad (4.71)$$

où:

* $[P^{1i}]$ est formée des lignes de $L^i$ correspondant aux valeurs positives de $\Lambda$ (équation 4.66), les autres étant remplacées par la forme linéarisée de la condition physique imposée $B(q)$ (équation 4.70),

$$\frac{\partial B}{\partial(\sqrt{g}q)} \qquad (4.72)$$

* $[P^{2i}]$ est formée des lignes de $L^i$ correspondant aux valeurs positives de $\Lambda$ (équation 4.66), les autres étant remplacées par le vecteur ligne nul, car il n'y a pas de différentielles spatiales dans l'équation 4.70.

* $[P^{3i}]$ est un vecteur colonne rempli de zéros, puisque le second membre de l'équation 4.66 est nul, exceptées les lignes associées à une valeur négative de $\Lambda$, qui reçoivent l'écart au niveau imposé (équation 4.70),

$$-\frac{1}{\Delta t} \left\{ B(q^n) - B^*(t^{n+1}) \right\} \qquad (4.73)$$

Après multiplication par l'inverse de la matrice $[P^{1i}]$, l'équation 4.71 se discrétise comme suit en suivant une formulation différence finie en espace curviligne:

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

C'est cette équation, formellement proche de celle obtenue pour un noeud du coeur du domaine (équation 3.26), qui sera à résoudre en un noeud de la frontière.

Les termes dans la direction $\xi^i$ sont discrétisés de manière décentrée depuis l'intérieur du domaine, puisqu'ils sont associés à une information sortante. Dans une approche MUSCL, ceci consiste à utiliser l'expression complète du flux convectif $F_c^i$, mais en prenant un champ aérodynamique interpolé depuis l'intérieur. Les termes tangents à la surface, et les flux visqueux et turbulents, sont quant à eux traités comme ils le sont dans le coeur du domaine.

Ainsi les termes convectifs $F_c^j$ apparaissant dans l'équation 4.74 sont évalués par:
$$
\begin{align}
F_c^j(\xi^j \pm \frac{1}{2}) &= F_c^j(Q_{\xi^j \pm \frac{1}{2}}^{-}) && \text{si } j=i \nonumber \\
&= F_c^{j+}(Q_{\xi^j \pm \frac{1}{2}}^{-}) + F_c^{j-}(Q_{\xi^j \pm \frac{1}{2}}^{+}) && \text{si } j \neq i \tag{4.75}
\end{align}
$$

Remarques complémentaires

* L'interpolation décentrée perpendiculairement à la frontière peut être du premier ou du second ordre. D'après Gustafsson [23], ceci n'influence pas la précision globale du schéma en espace. En pratique, le premier ordre apparaît comme étant un tout petit peu plus stable.

* L'écart au niveau imposé 4.73 peut sembler superflu lors de la recherche d'un état stationnaire pour lequel $B^*(t)$ est constant. Cependant, ce terme se révèle tout aussi indispensable dans ce cas, qu'il l'est pour une simulation instationnaire. En effet, sans lui, la condition imposée $B(q)$ est susceptible de dériver au cours du temps par rapport à son niveau de référence, soit du fait d'accumulation d'erreurs numériques, soit sous l'action d'une perturbation générée au sein du domaine de calcul et qui interagit avec la frontière. Ceci peut, par exemple, être illustré dans le cas d'un tube à choc dont une extrémité est maintenue à pression statique constante et sur laquelle le choc vient impacter et se transformer en une onde de détente. La figure 4.19 présente l'évolution du profil de pression statique en ce bout de tube au cours du temps, avec et sans le terme $[P^{3i}]$; alors que le niveau de masse volumique y est visible figure 4.20. Si nous excluons l'interval très court après la date d'interaction, $t \simeq 0.31$, où le niveau réel s'écarte de sa valeur théorique, une amélioration indéniable est constatée. D'autre part, ce terme est la clef du contrôle temporel des conditions aux limites instationnaires, que se soient une fluctuation de pression totale ou la déformation d'une paroi solide.

* Un problème d'ordre numérique, inhérent à l'approche suivie, a été soulevé par ailleurs mais non résolu au cours de cette étude. Il résulte de la fusion des équations 4.66 et 4.70 qui nous conduit à mélanger des flux et des conditions physiques dans le second terme de l'équation 4.74, sans que l'ordre de grandeur des uns et des autres ne soit pris en compte. Sa résolution passe par l'étude précise de la matrice $[P^{1i}]$, en s'inspirant des remarques faites lors du choix des vecteurs propres, paragraphe 4.3.1.

4.3.3 Différentes conditions physiques fortes

Après avoir vu la méthode générale, nous allons passer en revue quelques confi-
gurations intéressantes. Pour chacune seront présentés le nombre et la nature des
conditions physiques à imposer ainsi que les formes linéarisées associées.
<figure><img src="image_55.png" /><figcaption>Figure 4.19: Choc sortant à l'atmosphère. Pression statique en bout de tube</figcaption></figure>

<figure><img src="image_56.png" /><figcaption>Figure 4.20: Choc sortant à l'atmosphère. Masse volumique en bout de tube</figcaption></figure>
<figure><img src="image_57.png" /><figcaption>Figure 4.21: Frontières libres</figcaption></figure>

## Frontière libre

* Définition

Un noeud d'une frontière libre est tel que du fluide entre ou sorte du domaine de calcul en ce point (figure 4.21).

* Exemples d'utilisation

Les faces d'entrée et de sortie d'une tuyère, ou encore l'amont et l'aval d'une roue de turbomachine correspondent à cette configuration. Elle permet d'autre part de simuler une frontière à l'infini qui entoure un profil isolé. Un jet de refroidissement à travers une aube de turbine est aussi pris en compte par cette condition.

* Nombre de conditions physiques à imposer

Les valeurs propres de $K^i$ s'écrivent sous forme générique, $\epsilon$ valant -1, 0 ou +1:

$$ \lambda = a \| \vec{a}^i \| (M^i + \epsilon) \quad (4.76) $$

Ainsi, selon le sens de l'écoulement à travers la frontière, et son régime, quatre situations existent. Elles sont regroupées dans le tableau 4.1, qui indique l'intervalle du nombre de Mach et le nombre de conditions physiques à imposer, i.e. le nombre de valeurs propres négatives.

* Conditions physiques imposées

Deux cas sont à distinguer selon qu'une sortie ou une entrée est considérée.
<table><caption>Table 4.1: Nombre de conditions physiques pour une frontière libre</caption><thead><tr><th></th><th>Régime subsonique</th><th>Régime supersonique</th></tr></thead><tbody><tr><th>Entrée</th><td>-1 ≤ M' ≤ 0<br>4 conditions</td><td>M' ≤ -1<br>5 conditions</td></tr><tr><th>Sortie</th><td>0 ≤ M' ≤ 1<br>1 condition</td><td>1 ≤ M' ≤ 1<br>0 condition</td></tr></tbody></table>

## - Traitement d'une sortie

Les informations franchissant une sortie supersonique ne peuvent remonter depuis l'extérieur du domaine de calcul. Ceci est cohérent avec le fait que les cinq valeurs propres soient positives, et que donc aucune condition physique doit être imposée. Ainsi, le champ aérodynamique en une telle frontière dépend uniquement du coeur du domaine.

Par contre en régime subsonique, une onde unique provient du reste de l'univers. Le choix de la condition physique destinée à la remplacer est vaste, mais non sans conséquence sur l'existence et l'unicité de la solution, ou encore la stabilité du schéma numérique. Ce problème est complexe, et il ne sera pas détaillé ici [27, pp. 349-353]. Traditionnellement en turbomachine, la pression statique de sortie est imposée. Nous avons ainsi:

$$B(q) = p \qquad (4.77)$$

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

la forme linéarisée associée étant:

$$\frac{\partial B}{\partial(\sqrt{gq})} = \frac{\gamma - 1}{\sqrt{g}} \begin{vmatrix} \frac{1}{2} \vec{V}^2 & & & \\ & -u_1 & & \\ & & -u_2 & \\ & & & -u_3 \\ & & & 1 \end{vmatrix} \qquad (4.79)$$

## - Traitement d'une entrée

Quatre ou cinq conditions physiques doivent être imposées. Encore une fois, leur choix ne peut être aléatoire car le fait que le problème soit bien posé en dépend. Un critère mathématique simple permet cependant d'évaluer la viabilité de l'ensemble des conditions physiques retenu. En effet, les formes linéarisées des conditions physiques introduites dans $[P^{1i}]$ doivent être une famille libre de vecteurs, afin de ne pas empêcher l'inversion de cette matrice. D'un point de vue pratique en turbomachine, il est courant que soient imposés en compressible les niveaux d'énergie et la direction de l'écoulement à travers une entrée.
L'énergie dynamique y est caractérisée par la pression totale $p_T$:

$$B(q) = p_T \quad (4.80)$$

$$= \left( \frac{\gamma - 1}{\gamma} \right)^{\frac{\gamma}{\gamma - 1}} p^{\frac{1}{\gamma - 1}} [p + \rho E]^{\frac{\gamma}{\gamma - 1}} \quad (4.81)$$

avec:

$$\frac{\partial B}{\partial(\sqrt{g}q_k)} = \left(\frac{\gamma - 1}{\gamma}\right)^{\frac{1}{\gamma - 1}} \left[1 + \frac{\rho E}{p}\right]^{\frac{1}{\gamma - 1}} \left\{ \left(1 - \frac{1}{\gamma}\left[1 + \frac{\rho E}{p}\right]\right) \frac{\partial p}{\partial(\sqrt{g}q_k)} + \frac{1}{\sqrt{g}} \delta_{k5} \right\} \quad (4.82)$$

où $\frac{\partial p}{\partial(\sqrt{g}q_k)}$ est donnée par l'équation 4.79, et $\delta_{ij}$ vaut 0 pour toutes les composantes $q_i$, exceptée la $j$-ième pour laquelle il vaut 1.

La température totale $T_T$ représente quant à elle l'énergie thermique:

$$B(q) = T_T \quad (4.83)$$

$$= \frac{1}{\gamma C_v \rho} [p + \rho E] \quad (4.84)$$

avec:

$$\frac{\partial B}{\partial(\sqrt{g}q_k)} = \frac{1}{\gamma C_v \rho} \left\{ \frac{\partial p}{\partial(\sqrt{g}q_k)} + \frac{1}{\sqrt{g}} \left( -\frac{1}{\rho} [p + \rho E] \delta_{k1} + \delta_{k5} \right) \right\} \quad (4.85)$$

A ces deux conditions, deux autres (resp. trois) fixant la direction de l'écoulement doivent être adjointes. Ceci peut se faire soit à l'aide du module des composantes de la vitesse, soit en conservant des angles de calage.

Dans le premier cas, soit $u_i$, composante de la vitesse $\vec{v}$ dans la direction cartésienne $x_i$, la grandeur à imposer. Nous avons:

$$B(q) = u_i \quad (4.86)$$

avec:

$$\frac{\partial B}{\partial(\sqrt{g}q)} = \frac{1}{\sqrt{g} \rho} \begin{vmatrix} -u_i & & \\ & \delta_{i1} & \\ & & \delta_{i2} \\ & & \delta_{i3} \\ & & 0 \end{vmatrix} \quad (4.87)$$
<figure><img src="image_58.png" /><figcaption>Figure 4.22: Angle de calage dans le plan $(x_i, x_j)$</figcaption></figure>

Dans le second cas, soit $\alpha_{ij}$ l'angle de calage formée par les composantes $u_i$ et $u_j$ de la vitesse $\vec{V}$ dans le plan cartésien $(x_i, x_j)$ (figure 4.22). La tangente de cet angle peut être utilisée comme condition physique, bien qu'il n'y ait pas d'équivalence stricte entre conserver un angle, défini sur $[0; 2\pi]$, et une tangente, correspondant à un intervalle angulaire de $[-\frac{1}{2}\pi; \frac{1}{2}\pi]$. La grandeur imposée est alors:

$$B(q) = \tan(\alpha_{ij}) \quad (4.88)$$

$$= \frac{u_j}{u_i} \qquad (4.89)$$

Il est à noter que cette expression n'est pas définie pour $u_i = 0$. C'est pourquoi il souhaitable dans ce cas de lui préférer l'expression précédente exprimant la conservation de $u_i$. La forme linéarisée associée est:

$$\frac{\partial B}{\partial(\sqrt{g}q_k)} = \frac{1}{\sqrt{g\rho u_i^2}} \{u_i\delta_{k(j+1)} - u_j\delta_{k(i+1)}\} \quad (4.90)$$

### • Remarques

La nature d'une frontière libre, entrée ou sortie, subsonique ou supersonique, n'est pas imposée, mais ce sont les conditions locales d'écoulement qui, à travers les valeurs propres, déterminent le traitement approprié. Ceci est nécessaire en instationnaire, où le champ aérodynamique à la frontière peut beaucoup changer, soit naturellement, lors du passage d'une perturbation, soit sous l'action d'une contrainte extérieure, par exemple une montée en régime du subsonique au supersonique.
<figure><img src="image_59.png" /><figcaption>Figure 4.23: Frontière glissante adiabatique</figcaption></figure>

## Frontière glissante adiabatique

* Définition

Une frontière glissante est telle que l'écoulement y est tangent en tout noeud (figure 4.23). De plus, le flux de chaleur perpendiculaire à cette surface est nul, assurant ainsi l'adiabaticité. Cette configuration concerne uniquement un fluide parfait, non visqueux et ne conduisant pas la chaleur, dont le mouvement est décrit par les équations d'Euler. Seuls demeurent les flux convectifs dans les équations de conservation.

* Exemples d'utilisation

Elle permet de simuler toute paroi solide guidant, déviant, ou faisant obstacle au fluide. C'est le parfait opposé d'une frontière libre. Cette condition est ainsi utilisée pour prendre en compte les parois d'une tuyère, le contour d'un profil, ou encore la présence d'un moyeu et d'un carter dans un canal inter-aubes.

* Nombre de conditions physiques à imposer

Par définition, la vitesse relative du fluide est tangente à la frontière, et donc perpendiculaire à sa normale. Si bien que:

$$V^i = (\vec{V} - \vec{W}) \cdot \vec{\sigma}^i = 0 \qquad (4.91)$$
Ainsi, quatres valeurs propres sont positives ou nulles, la dernière étant strictement négative. Nous sommes ici confronté à un cas limite de la méthode développée. En effet, d'après l'interprétation faite de l'équation 4.64, trois des ondes traversant la frontière sont de vitesse nulle, et stagnent donc là. Faut-il les considérer comme entrantes, et donc supprimer les équations associées, ou sortantes? La pratique montre qu'une seule contrainte doit être imposée. Nous allons donc agir comme si quatres ondes quittaient le domaine de calcul.

### • Conditions physiques imposées

L'unique condition physique à imposer est dictée par la définition: il s'agit de la tangence à la paroi de la vitesse. Ainsi,

$$B(q) = V^i \qquad (4.92)$$

$$B^*(t) = 0 \qquad (4.93)$$

dont la forme linéarisée est:

$$\frac{\partial B}{\partial(\sqrt{g}q)} = \frac{1}{\sqrt{g} \rho} \begin{vmatrix} -(V^i - \xi_i^i) \\ \xi_{x1}^i \\ \xi_{x2}^i \\ \xi_{x3}^i \\ 0 \end{vmatrix} \qquad (4.94)$$

### • Remarques

Le champ aérodynamique extrapolé sur les noeuds supplémentaires n'assure pas, en particulier, un flux de chaleur nul à travers la paroi. Ceci est une source potentielle d'erreurs si un décentrement d'ordre supérieur à un est utilisé perpendiculairement à la frontière. L'extrapolation pourrait être corrigée pour résoudre ce problème. Une meilleure solution serait d'imposer des conditions miroir locales à la place de l'extrapolation, de la même manière qu'une condition de symétrie permet de gérer une frontière glissante plane. Il faudrait considérer alors la normale au noeud, et non plus une normale unique pour toute la surface limite.

D'autre part, des variations erronées d'énergie très localisées près des parois glissantes ont été constatées. Ceci est clairement visible sur les figures 4.24 et 4.25 où sont représentés les écarts relatifs de pression totale et de température totale autour d'un cylindre. Un écoulement de fluide parfait à un nombre de Mach de 0.3 étant simulé dans ce cas, ces deux variables devraient être nulles. Plusieurs causes à ce comportement peuvent être avancées. Il peut résulter du traitement au premier ordre des frontières, de l'extrapolation simpliste faite sur les plans supplémentaires, d'un mauvais raccord entre le schéma
<figure><img src="image_60.png" /><figcaption>Figure 4.24: Cylindre à M = 0.3. Erreur sur la pression totale</figcaption></figure>
<figure><img src="image_61.png" /><figcaption>Figure 4.25: Cylindre à M = 0.3. Erreur sur la température totale</figcaption></figure>
<figure><img src="image_62.png" /><figcaption>Figure 4.26: Frontière adhérente isotherme. Profil de vitesse</figcaption></figure>

du coeur du domaine et celui utilisé sur les surfaces limites, ou encore de la séparation de flux de Van Leer [26]. Dadone [15] a simulé une configuration semblable, mettant en évidence le même comportement. Il a conclu quant à lui que la qualité du traitement d'une frontière glissante était fortement tributaire de l'extrapolation réalisée sur les noeuds fictifs. Il montre en particulier que plus cette opération est entachée d'erreurs, plus la solution est dépendante du maillage. L'extrapolation linéaire sur le maillage relativement grossier que nous utilisons là, a donc de fortes probabilités d'être à la source des défauts constatés.

## Frontière adhérente isotherme

* Définition

Une frontière adhérente isotherme est telle que d'une part la vitesse du fluide relativement à cette surface est nulle (figure 4.26), d'autre part la température de paroi y est imposée (figure 4.27).

* Exemples d'utilisation

Cette configuration permet de décrire le mouvement d'un fluide visqueux sur une paroi dont la température serait imposée. En particulier, une aube de turbine refroidie est simulable par cette condition.

* Nombre de conditions physiques à imposer

L'adhérence du fluide à la paroi fait que sa vitesse est égale à celle de la frontière. Ainsi:

$$ \bar{V} = \bar{W} \qquad (4.95) $$
<figure><img src="image_63.png" /><figcaption>Figure 4.27: Frontière adhérente isotherme. Température statique</figcaption></figure>

et donc:

$$V^i = 0 \qquad (4.96)$$

Aussi, comme pour une frontière glissante, quatres valeurs propres sont positives ou nulles, la dernière étant strictement négative. Encore une fois se pose la question du traitement des trois ondes stationnaires sur la surface. Mais dans ce cas, la pratique montre que ce sont quatre conditions physiques et non plus une seule, qui doivent être fixées.

### • Conditions physiques imposées

Le choix des conditions physiques à imposer est sans équivoque d'après la définition. Il s'agit des trois composantes de la vitesse du fluide et de la température statique.

Lors du traitement d'une frontière libre, nous avons vu comment fixer une composante de $\vec{V}$. La même technique est à appliquer ici, si ce n'est que le niveau de référence $B^*(t)$ est alors égal à la vitesse de la paroi $\vec{W}$. Ainsi,

$$B(q) = u_i \qquad (4.97)$$

$$B^*(t) = w_i \qquad (4.98)$$

La température statique se traite quant à elle comme la pression statique. Nous avons en effet:

$$B(q) = T \qquad (4.99)$$
<figure><img src="image_64.png" /><figcaption>Figure 4.28: Frontière adhérente à flux de chaleur imposé</figcaption></figure>

$$ = \frac{1}{r} \frac{p}{\rho} \qquad (4.100) $$

avec:

$$ \frac{\partial B}{\partial(\sqrt{g}q_k)} = \frac{1}{r} \frac{\partial p}{\partial \rho} \left\{ \frac{\partial p}{\partial(\sqrt{g}q_k)} - \frac{p}{\sqrt{g}\rho} \delta_{k1} \right\} \qquad (4.101) $$

## Frontière adhérente à flux de chaleur imposé

* Définition

Une frontière adhérente à flux de chaleur imposé est telle que d'une part la vitesse du fluide relativement à cette surface est nulle (figure 4.28), d'autre part le flux de chaleur qui la traverse perpendiculairement est fixé. Si ce dernier est nul, la surface est alors adiabatique.

* Exemples d'utilisation

Cette configuration est le pendant d'une frontière adhérente isotherme, si ce n'est que dans ce cas, ce n'est plus le niveau de température statique qui est imposé, mais son gradient spatial. Une paroi mouillée par un fluide visqueux et dont les échanges thermiques sont contrôlés peut être décrite par cette condition. Le principal usage qui en est fait concerne cependant des frontières adiabatiques, qui n'échangent pas, ou peu, de chaleur avec le milieu fluide. Cette approximation est très courante en turbomachine pour modéliser les surfaces solides, telles que les aubes, le moyeu et le carter.

* Nombre de conditions physiques à imposer
La thermique n'ayant aucune incidence sur les valeurs propres, la discussion faite pour une frontière adhérente isotherme est valable dans ce cas. Ainsi, nous en déduisons que quatre conditions physiques doivent être fixées.

## • Conditions physiques imposées

D'après la définition, les trois composantes de la vitesse du fluide et le flux de chaleur doivent être imposés. Par rapport à une frontière adhérente isotherme, seule la contrainte thermique change. Aussi, la vitesse sera traitée comme dans le cas précédent.

L'imposition du flux de chaleur $\Phi_w$ soulève quant à elle quelques problèmes. En effet, la méthode développée et utilisée jusqu'à présent permet de fixer une grandeur, mais non ses variations spatiales. Or il nous faut conserver ici le gradient de température perpendiculaire à la frontière tel que:

$$ -\kappa \frac{\partial T}{\partial \vec{n}} = \Phi_w \qquad (4.102) $$

En appliquant certaines relations de géométrie différentielle en coordonnées curvilignes (annexe D.5.2) et en remarquant que $\vec{n} = \vec{a}^i$, il est possible d'établir que:

$$ \frac{\partial T}{\partial \vec{n}} = \frac{1}{\| \vec{a}^i \|} \sum_{j=1}^{3} \vec{a}^i \cdot \vec{a}^j \frac{\partial T}{\partial \xi^j} \qquad (4.103) $$

Si bien que la grandeur à conserver $B(q)$ et son niveau de référence $B^*(t)$ sont définis par:

$$ B(q) = - \frac{\kappa}{\| \vec{a}^i \|} \sum_{j=1}^{3} \vec{a}^i \cdot \vec{a}^j \frac{\partial T}{\partial \xi^j} \qquad (4.104) $$

$$ B^*(t) = \Phi_w \qquad (4.105) $$

Or ces relations ne peuvent être fusionnées avec l'équation 4.66 puisqu'elles ne contiennent pas de différentielle temporelle, et d'autre part, leur linéarisation en temps fait apparaître des dérivées croisées temps-espace, impossibles à introduire elles-aussi dans le système d'équations. La difficulté augmente encore lorsque nous envisageons un maillage mobile et déformable. Aucune solution satisfaisante sur le plan théorique n'a été trouvée à ce problème. Cependant, quelques approximations permettent quant même de traiter une telle frontière à condition qu'elle soit adiabatique.

Tout d'abord, l'extrapolation du champ aérodynamique doit être modifiée afin d'assurer la nullité du flux de chaleur. Nous voulons donc que:
$$ \sum_{j=1}^{3} \vec{a}^i \cdot \vec{a}^j \frac{\partial T}{\partial \xi^j} = 0 \qquad (4.106) $$

Or l'expression de la température en fonction des variables conservatives étant:

$$ T = \frac{1}{C_v} \left\{ \frac{\rho E}{\rho} - \frac{1}{2} \frac{(\rho \vec{V})^2}{\rho^2} \right\} \qquad (4.107) $$

il est possible d'établir que:

$$ \frac{\partial T}{\partial \xi^j} = \frac{1}{C_v} \frac{1}{\rho} \left\{ [\vec{V}^2 - E] \frac{\partial \rho}{\partial \xi^j} - \vec{V} \frac{\partial \rho \vec{V}}{\partial \xi^j} + \frac{\partial \rho E}{\partial \xi^j} \right\} \qquad (4.108) $$

Cette équation se simplifie en remarquant que la vitesse est nulle à la paroi, d'où:

$$ \frac{\partial T}{\partial \xi^j} = \frac{1}{C_v} \frac{1}{\rho} \left\{ -E \frac{\partial \rho}{\partial \xi^j} + \frac{\partial \rho E}{\partial \xi^j} \right\} \qquad (4.109) $$

Ainsi, l'équation 4.106 est équivalente à:

$$ \sum_{j=1}^{3} \vec{a}^i \cdot \vec{a}^j \frac{\partial \rho E}{\partial \xi^j} - E \sum_{j=1}^{3} \vec{a}^i \cdot \vec{a}^j \frac{\partial \rho}{\partial \xi^j} = 0 \qquad (4.110) $$

Cette relation est impliquée par, A représentant ρ ou ρE :

$$ \sum_{j=1}^{3} \vec{a}^i \cdot \vec{a}^j \frac{\partial A}{\partial \xi^j} = 0 \qquad (4.111) $$

C'est cette équation qui est utilisée pour extrapoler le champ aérodynamique sur les noeuds supplémentaires. Il est possible de la dicrétiser comme suit (figure 4.29):

$$ \sum_{j=1}^{3} \vec{a}^i \cdot \vec{a}^j \frac{A (\xi^j + 1) - A (\xi^j - 1)}{2 \Delta \xi^j} = 0 \qquad (4.112) $$

Ce qui détermine pour nouvelle relation d'extrapolation, sachant que $\Delta \xi^j = 1$:

$$ A (\xi^i + 1) = A (\xi^i - 1) - \frac{1}{\vec{a}^i} \sum_{\substack{j=1 \\ j \neq i}}^{3} \vec{a}^j \cdot \vec{a}^i [A (\xi^j + 1) - A (\xi^j - 1)] \qquad (4.113) $$
<figure><img src="image_65.png" /><figcaption>Figure 4.29: Extrapolation assurant l'adiabaticité</figcaption></figure>

Appliquée à la densité et à l'énergie totale, la quantité de mouvement restant inchangée, elle assure ainsi que le flux de chaleur calculé perpendiculairement à la frontière est nul. Il est à noter que cinq noeuds du domaine de calcul sont utilisés, alors que deux uniquement interviennent dans l'extrapolation classique.

Par ailleurs, nous n'imposons que trois conditions physiques, les composantes de la vitesse, à la place de quatre. Le vecteur propre gauche de $K^i$ associé à la valeur propre $V^i - a \parallel \vec{a}^i \parallel$ est supprimé, ainsi que deux des trois vecteurs propres gauches associés à la valeur propre $V^i$. Les deux derniers vecteurs propres gauche à conserver sont choisis de telle façon que $[P^{1i}]$ soit inversible.

Ces deux modifications, d'une part de l'extrapolation, d'autre part de la construction du système d'équations, nous permettent de simuler une frontière adhérente adiabatique avec une précision raisonnable.

## • Remarques

Le fait que nous traitions une condition d'adhérence, propre à un fluide visqueux, avec une approche bâtie sur les équations d'Euler, n'est pas étrangère au fait que, premièrement, la thermique n'ait aucune incidence puisqu'un fluide parfait ne conduit pas la chaleur, deuxièmement, que les ondes stationnaires sur la surface limite ne soient pas traitées comme dans le cas d'une frontière glissante adiabatique bien que les valeurs propres soient les mêmes, enfin que nous ne sachions que faire du flux de chaleur, qui fait intervenir une relation de voisinage spatiale, alors que tous les autres cas mettent en jeu des conditions ponctuelles.
4.3.4 Avantages et inconvénients

La méthode développée présente de nombreux avantages par rapport à une approche plus classique et à première vue moins complexe.

Tout d'abord, le système d'équations discrétisées en espace obtenu à la frontière (équation 4.74) est formellement semblable à celui du coeur (équation 3.26). Ainsi, nous pouvons utiliser le même schéma en temps en tout noeud, et résoudre un problème uniformisé, la distinction en espace entre ces deux formulations se faisant à l'assemblage du système et non au moment de sa résolution. Ceci a peu d'importance lors d'un traitement temporel explicite, car il n'y a pas alors couplage entre les différents noeuds à un même instant. Mais c'est primordial pour un schéma implicite afin de maintenir la cohérence en temps du champ aérodynamique dans son ensemble.

D'autre part, un critère théorique simple a été dégagé pour choisir les équations à supprimer et celles à conserver, quelles que soient la géométrie et les conditions d'écoulement. Par exemple, une frontière glissante sera traitée de manière identique, qu'elle soit horizontale, verticale, ou en biais... Il en est de même pour une surface libre en laquelle les valeurs propres nous renseignent naturellement sur le sens et le régime de l'écoulement. Nous n'implémentons pas séparément ces différents cas.

Dans le même esprit d'uniformité, un formalisme unique est utilisé pour toutes les configurations. Ceci en facilite grandement l'implémentation puisque seules les matrices $[P^{1i}]$, $[P^{2i}]$ et $[P^{3i}]$ changent de l'une à l'autre.

Enfin le niveau d'une certaine grandeur $B(q)$ imposée est contrôlable de façon très précise en temps. Il peut être surveillé et modifié de manière totalement externe au traitement des conditions aux limites car il suffit d'agir sur sa valeur de référence $B^*(t)$. Nous disposons ainsi d'un moyen d'action d'une extrême facilité et d'une grande souplesse indispensable lors d'un calcul instationnaire.

Cependant, un certain nombre de reproches peuvent être formulés à l'encontre de ces conditions aux limites fortes.

En effet, elles se révèlent parfois trop rigides, gênant ou bloquant la convergence vers un état stationnaire. Ainsi, nous avons vu que les réflexions d'ondes sont très bien gérées. Ceci est un bon point en instationnaire lors de l'étude de l'impact d'un choc par exemple, mais cela devient un inconvénient lorsque des ondes parasites, dues souvent au champ aérodynamique initial, mettent très longtemps à se dissiper, en particulier pour un fluide parfait. Dans ce cas, une solution partielle consiste à débuter le calcul au premier ordre en espace afin de profiter d'un amortissement accru par une viscosité numérique excessive.

D'autre part, il existe une contrainte au niveau des conditions physiques qui sont parfois mal connues. Par exemple, la répartition bidimensionnelle de pression statique en sortie subsonique d'une roue est souvent remplacée par un niveau moyen uniforme faute de mieux. Le même problème se pose lorsqu'il faudrait laisser sortir une perturbation physique instationnaire, telle qu'un choc ou un sillage se positionnant.
Enfin, cette approche est inadaptée pour traiter les arêtes et les coins du domaine de calcul, en lesquels plusieurs conditions aux limites dans plusieurs directions devraient être prises en compte simultanément. A cause de la modification des équations initiales, en particulier due à la suppression de termes dans la direction sécante à la frontière, mais aussi dans les directions tangentielles, le système qui en résulte à des propriétés très différentes des équations d'Euler. Il n'est plus hyperbolique en autres. Il n'est donc pas possible d'appliquer la même méthode de suppression/remplacement une nouvelle fois dans une autre direction.

Il est à noter que nous nous intéressons à des ondes se propageant perpendiculairement à la surface caractéristique, colinéairement à la normale. Ce n'est donc pas forcément le long de la ligne de coordonnée sécante à la surface limite, sauf si le maillage est orthogonal. Une légère erreur lors du décentrement peut en résulter, conduisant à un problème d'alignement des ondes et des cellules d'intégration semblable à celui rencontré lors de la séparation de flux (paragraphe 3.4.3).

## 4.4 Conditions physiques faibles

Ce traitement s'applique aux noeuds en lesquels une contrainte extérieure est imposée de manière dite faible. Contrairement à ce qui est fait dans le cadre des conditions physiques fortes, nous ne cherchons plus à introduire de nouvelles équations fixant certaines grandeurs dans le système différentiel initial, mais nous déterminons de manière indirecte les informations en provenance de l'extérieur du domaine de calcul.

Il est nécessaire de développer une nouvelle formulation car une arête est l'intersection de deux frontières (resp. un coin est l'intersection de trois frontières). Des ondes se propageant dans deux (resp. trois) directions, chacune perpendiculairement à une surface limite, doivent alors être mises en évidence. Or l'approche forte ne peut être appliquée qu'une seule fois, après quoi le système dégénère et il ne décrit plus des ondes... Il faut donc définir une nouvelle stratégie.

Pour cela, nous suivons une approche bâtie sur les équations de compatibilité proche de celle présentée précédemment. Mais nous traitons différemment les ondes entrantes. Dans le premier cas, nous avons supprimé ce que nous ne connaissions pas et l'avons remplacé par une autre donnée. Ici, nous exprimons l'inconnu à l'aide des informations disponibles [65]. Si bien que nous tendons à imposer les effets du reste de l'univers plutôt que les causes originelles. C'est en ce sens que nous parlons de conditions faibles.

### 4.4.1 Formulation caractéristique multi-directionnelle

Notre première tache est donc de trouver une nouvelle technique qui ne dénature pas le système initial, pour pouvoir l'appliquer plusieurs fois afin de traiter les arêtes et les coins. L'idée directrice est que l'équation 4.65 renferme une part des informations qui ont été utilisées pour imposer les conditions physiques fortes. En effet, le terme $\ell^i$
est lié aux ondes qui traversent la frontière de normale $\vec{a}^i$. Cependant, cette équation n'étant qu'une autre écriture des équations de conservation, elle décrit aussi les ondes se propageant dans d'autres directions. D'autre part, l'approche ayant conduit à son obtention étant non destructive, au sens où n'en changeant pas les propriétés, il est possible de l'appliquer plusieurs fois.

Ainsi, le système:

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

peut-il être réécrit [65]:

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

où:

* $n$ est le nombre de frontières considérées.

$n=1$ pour une face

$n=2$ pour une arête

$n=3$ pour un coin

* $\ell^i$ est défini à partir des valeurs propres $\Lambda^i$ et des vecteurs propres gauches $L^i$ de la matrice jacobienne $K^i$ associée au flux convectif $F_c^i$ dans la direction $\xi^i$ par:

$$ \ell^i = L^i \frac{\partial F_c^i}{\partial \xi^i} \quad (4.116) $$

$$ = \Lambda^i \frac{\partial}{\partial \xi^i} (\sqrt{g} w^i) \quad (4.117) $$

Le système 4.115 n'est qu'une réécriture des équations des Navier-Stokes, mais il laisse apparaître la notion d'ondes se déplaçant dans $n$ directions tangentes aux $\vec{a}^i$, transportant les variables caractéristiques $w^i$.

### 4.4.2 Imposition faible d'une condition aux limites

Il faut maintenant reformuler le système 4.115 pour y inclure les contraintes extérieures. Pour cela, une stratégie proche dans l'esprit de celle mise en oeuvre pour les conditions physiques fortes est suivie. Elle repose elle-aussi sur la mise en évidence de deux familles d'ondes traversant une frontière dans une direction. 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 révèlent quelle condition physique sur la surface limite doit être imposée. Cette dernière va nous permettre de relier les deux familles d'ondes, et ainsi d'exprimer l'information entrante en fonction des informations sortantes.

Ainsi, pour chaque direction $\xi^i$, le signe des valeurs propres $\Lambda^i$ nous renseigne sur le sens de propagation des informations le long de la normale $\tilde{a}^i$. Soit $\lambda_k^i$, la $k$-ième valeur propre de $K^i$. Notre but est de déterminer l'expression de $\ell_k^i$ qui y est associé, pour la réintroduire dans l'équation 4.115.

* Si $\lambda_k^i$ est positive:

L'onde de vitesse $\lambda_k^i$ sort du domaine de calcul, et elle est de ce fait connue.
Nous avons:

$$ \ell_k^i = L_k^i \frac{\partial F_c^i}{\partial \xi^i} \qquad (4.118) $$

où $\frac{\partial F_c^i}{\partial \xi^i}$ est à calculer de manière décentrée depuis l'intérieur du domaine.

* Si $\lambda_k^i$ est négative:

L'onde associée apporte l'influence du reste de l'univers et elle n'est donc pas directement évaluable. Cependant, elle peut être reliée aux ondes sortantes par une condition physique judicieusement choisie. En effet, soient $B(q)$ une grandeur et $B^*(t)$ son niveau à imposer. D'après les développements faits autour des conditions physiques fortes, il est possible de linéariser au premier ordre en temps:

$$ B(q) - B^*(t) = 0 \qquad (4.119) $$

en:

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

Or 4.115 exprime que:

$$ - \frac{\partial}{\partial t} (\sqrt{gq}) = \sum_{j=1}^{n} R^j \ell^j + \sum_{j=n+1}^{3} \frac{\partial F_c^j}{\partial \xi^j} + S \qquad (4.121) $$

et donc:
$$ \sum_{j=1}^{n} \left[ \frac{\partial B}{\partial(\sqrt{g})} R^j \right] \ell^j = \frac{1}{\Delta t} \left\{ B(q^n) - B^*(t^{n+1}) \right\} - \left[ \frac{\partial B}{\partial(\sqrt{g})} \right] \left\{ \sum_{j=n+1}^{3} \frac{\partial F_c^j}{\partial \xi^j} + S \right\} \quad (4.122) $$

où les $\frac{\partial B}{\partial(\sqrt{g})} R^j$ sont $n$ vecteurs ligne multipliant $n$ vecteurs colonne $\ell^j$. $\ell_k^i$ est une composante de l'un d'eux. Sa valeur en fonction des autres $\ell^j$ est obtenue par résolution de cette équation linéaire.

### Conclusion:

Les équations 4.118 et 4.122 peuvent être regroupées sous une forme générique. Soient $C_j^{ik}$, $P_j^{1ik}$, $P^{2ik}$, des vecteurs ligne,

$$ C_j^{ik} = [C_{jl}^{ik}]_{l \in [1; 5]} \quad (4.123) $$

$$ P_j^{1ik} = [P_{jl}^{1ik}]_{l \in [1; 5]} \quad (4.124) $$

$$ P^{2ik} = [P_l^{2ik}]_{l \in [1; 5]} \quad (4.125) $$

et $P^{3ik}$ un scalaire tels que:

- Si $\lambda_k^i$ est positive:

$$ C_{jl}^{ik} = \delta_{ij} \delta_{kl} \quad (4.127) $$

$$ P_j^{1ik} = \delta_{ij} L_k^i \quad (4.128) $$

$$ P_l^{2ik} = 0 \quad (4.129) $$

$$ P^{3ik} = 0 \quad (4.130) $$

- Si $\lambda_k^i$ est négative:

$$ C_j^{ik} = \frac{\partial B}{\partial(\sqrt{g})} R^j \quad (4.131) $$

$$ P_{jl}^{1ik} = 0 \quad (4.132) $$

$$ P^{2ik} = \frac{\partial B}{\partial(\sqrt{g})} \quad (4.133) $$

$$ P^{3ik} = \frac{1}{\Delta t} \{ B(q^n) - B^*(t^{n+1}) \} \quad (4.134) $$

Nous avons alors:
$$ \sum_{j=1}^{n} C_j^{ik} \ell^j = \sum_{j=1}^{n} P_j^{1ik} \frac{\partial F_c^j}{\partial \xi^j} + P^{3ik} - P^{2ik} \left\{ \sum_{j=n+1}^{3} \frac{\partial F_c^j}{\partial \xi^j} + S \right\} \quad (4.135) $$

Le $\ell_k^i$ que nous recherchons, associé à la valeur propre $\lambda_k^i$, est une des inconnues de cette équation linéaire.

L'accumulation des équations 4.135 associées aux cinq valeurs propres (k) dans chacune des n directions (i) forme un système de 5n équations à 5n inconnues, les $\ell_k^i$. Son inversibilité est assurée par un choix de conditions aux limites physiques compatibles entre elles. Ce problème d'existence et d'unicité d'une solution a les mêmes racines que celui rencontré lors de la construction de la matrice $[P^{1i}]$ pour les conditions aux limites fortes. Le système ainsi obtenu s'écrit:

$$ [C] \ell = [P^1] \left[ \frac{\partial F_c^i}{\partial \xi^i} \right] + [P^3] - [P^2] \left\{ \sum_{j=n+1}^{3} \frac{\partial F_c^j}{\partial \xi^j} + S \right\} \quad (4.136) $$

où:

* $\ell$ est un vecteur colonne, formé de n vecteurs colonne de longueur 5, les $\ell^i_{i \in [1:n]}$.

* $[\frac{\partial F_c^i}{\partial \xi^i}]$ est un vecteur colonne, formé de n vecteurs colonne de longueur 5, les différentes des flux convectifs sécants aux frontières $\frac{\partial F_c^i}{\partial \xi^i}_{i \in [1:n]}$.

et:

* $[C]$ est une matrice carrée, formée de $n \times n$ blocs $5 \times 5$, les $[C_j^i]_{(i,j) \in [1:n]}$, qui résultent du regroupement des cinq lignes $[C_j^{1i}]_{k \in [1:5]}$.

* $[P^1]$ est une matrice carrée, formée de $n \times n$ blocs $5 \times 5$, les $[P_j^{1i}]_{(i,j) \in [1:n]}$, qui résultent du regroupement des cinq lignes $[P_j^{1ik}]_{k \in [1:5]}$.

* $[P^2]$ est une matrice colonne, formée de n blocs $5 \times 5$, les $[P^{2i}]_{i \in [1:n]}$, qui résultent du regroupement des cinq lignes $[P^{2ik}]_{k \in [1:5]}$.

* $[P^3]$ est un vecteur colonne, formé de n vecteurs colonne de longueur 5, les $[P^{3i}]_{i \in [1:n]}$, qui résultent du regroupement des cinq scalaires $[P^{3ik}]_{k \in [1:5]}$.

Les expressions des $\ell^i$ qui résultent de l'inversion du système 4.136 sont introduites dans l'équation 4.115, qui devient:
$$
\begin{align*}
& \frac{\partial}{\partial t}(\sqrt{gq}) + \\
& [R][C]^{-1}[P^1]\left[\frac{\partial F_c^i}{\partial \xi^i}\right] + \\
& (I_5 - [R][C]^{-1}[P^2])\left\{\sum_{j=n+1}^{3} \frac{\partial F_c^j}{\partial \xi^j} + S\right\} + \\
& [R][C]^{-1}[P^3] = 0 \tag{4.137}
\end{align*}
$$

où $[R]$ est une matrice ligne, formée de $n$ blocs $5 \times 5$, les $[R^i]_{i \in [1:n]}$.

Après développement des différents produits de matrices bloc, ce système peut
s'écrire sous forme synthétique:

$$
\frac{\partial}{\partial t}(\sqrt{gq}) + \sum_{j=1}^{3} \tilde{P}_j \frac{\partial F_c^j}{\partial \xi^j} + \tilde{S} = 0 \quad (4.138)
$$

Finalement, après discrétisation par différence finie en espace curviligne, l'équa-
tion que nous avons à résoudre en un noeud de la frontière est:

$$
\frac{\partial}{\partial t}(\sqrt{gq}) + \sum_{j=1}^{3} \tilde{P}_j \left[ F_c^j(\xi^j + \frac{1}{2}) - F_c^j(\xi^j - \frac{1}{2}) \right] + \tilde{S} = 0 \quad (4.139)
$$

Elle est formellement proche de celle obtenue pour un noeud du coeur du domaine
(équation 3.26).

Les termes dans les directions $\xi^i$ sécantes à une frontière sont discrétisés de manière décentrée depuis l'intérieur du domaine, puisqu'ils sont associés à une information sortante. Dans une approche MUSCL, ceci consiste à utiliser l'expression complète du flux convectif $F_c^i$, mais en prenant un champ aérodynamique interpolé depuis l'intérieur. Les termes tangents à la surface, et les flux visqueux et turbulents, sont quant à eux traités comme ils le sont dans le coeur du domaine.

Ainsi les termes convectifs $F_c^j$ sont évalués par :

$$
\begin{align}
F_c^j(\xi^j \pm \frac{1}{2}) &= F_c^j(Q_{\xi^j \pm \frac{1}{2}}^{-}) && \text{si } j \in [1; n] \nonumber \\
&= F_c^{j+}(Q_{\xi^j \pm \frac{1}{2}}^{-}) + F_c^{j-}(Q_{\xi^j \pm \frac{1}{2}}^{+}) && \text{si } j \in [n+1, 3] \tag{4.140}
\end{align}
$$

4.4.3 Différentes conditions physiques faibles

Il est difficile d'être exhaustif sur le plan analytique car le traitement possible de une
à trois directions simultanément génère de nombreux cas. Par contre, l'implémenta-
tion est relativement simple, en particulier parce que l'ordre dans lequel sont traitées
les différentes frontières est sans importance. Par ailleurs, la discussion du nombre et
de la nature des conditions physiques à imposer pour chaque cas a déjà été faite dans
le cadre des conditions aux limites fortes. De plus à cette occasion ont été exprimées
<figure><img src="image_66.png" /><figcaption>Figure 4.30: Paroi glissante</figcaption></figure>

les formes linéarisées qui y sont asssociées. Aussi, nous allons nous contenter par la
suite d'illustrer la méthode proposée sur trois configurations représentatives.

## Paroi glissante

* **Définition**

Une frontière glissante est telle que l'écoulement y est tangent en tout noeud.
De plus, le flux de chaleur perpendiculaire à cette surface est nul, assurant
ainsi l'adiabaticité. Cette configuration concerne uniquement un fluide parfait,
non visqueux et ne conduisant pas la chaleur, dont le mouvement est décrit par
les équations d'Euler. Seuls demeurent les flux convectifs dans les équations
de conservation.

* **Exemples d'utilisation**

Elle permet de simuler toute paroi solide guidant, déviant, ou faisant obstacle
au fluide (figure 4.30). Cette condition est ainsi utilisée pour prendre en compte
les parois d'une tuyère, le contour d'un profil, ou encore la présence d'un moyeu
et d'un carter dans un canal inter-aubes.

* **Nombre et nature des conditions physiques à imposer**

Nous sommes en présence d'une seule frontière (n = 1) à laquelle la vitesse relative du fluide est tangente. Si bien que comme lors du traitement de cette configuration par l'approche forte, une seule condition aux limites physiques est imposée. Soit $\bar{a}$ la normale à la surface considérée. Nous avons:
$$V^i = (\vec{V} - \vec{W}) \cdot \vec{a}^i = 0 \quad (4.141)$$

* Expression des $\ell^i$

Les quatres premières valeurs propres de $K^i$ sont positives ou nulles, seule la dernière étant strictement négative. Les différentes matrices intervenant dans l'évaluation des $\ell^i$ ont alors pour expression:

$$[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 & \frac{a}{2\sqrt{g\rho}} & -\frac{a}{2\sqrt{g\rho}} \end{bmatrix} \quad (4.142)$$

La dernière ligne de cette matrice provient du produit de $\frac{\partial V^i}{\partial(\sqrt{g\rho})}$ par $R^i$.

$$[P^1] = \begin{bmatrix} L_1^i \\ L_2^i \\ L_3^i \\ L_4^i \\ 0 \end{bmatrix} \quad (4.143)$$

$$[P^2] = \frac{1}{\sqrt{g} \rho} \begin{bmatrix} 0 & 0 & 0 & 0 & 0 \\ 0 & 0 & 0 & 0 & 0 \\ 0 & 0 & 0 & 0 & 0 \\ 0 & 0 & 0 & 0 & 0 \\ -(V^i - \xi_i^i) & \xi_{x1}^i & \xi_{x2}^i & \xi_{x3}^i & 0 \end{bmatrix} \quad (4.144)$$

$$[P^3] = \begin{bmatrix} 0 \\ 0 \\ 0 \\ 0 \\ \frac{1}{\Delta t} V^i \end{bmatrix} \quad (4.145)$$

## Arête glissante/glissante

* Définition

Cette configuration résulte de l'intersection de deux frontières glissantes adiabatiques.

* Exemples d'utilisation

Elle permet de simuler le raccord entre deux parois solides guidant le fluide (figure 4.31). Elle se rencontre en particulier à l'emplanture d'une aile ou d'une
<figure><img src="image_67.png" /><figcaption>Figure 4.31: Arête glissante/glissante</figcaption></figure>

aube, ou encore dans les arêtes d'un canal de section quadrangulaire, telle une tuyère.

* Nombre et nature des conditions physiques à imposer

Nous sommes en présence de deux frontières ($n=2$) telles que la vitesse relative du fluide est tangente à chacune d'elles. Une condition physique doit donc être imposée dans chaque direction. Soient $\vec{a}^i$ et $\vec{a}^j$ les normales aux frontières considérées. Nous avons:

$$V^i = (\vec{V} - \vec{W}) \cdot \vec{a}^i = 0 \qquad (4.146)$$

$$V^j = (\vec{V} - \vec{W}) \cdot \vec{a}^j = 0 \qquad (4.147)$$

* Expression des $\ell^i$

Les quatres premières valeurs propres de $K^i$ et de $K^j$ sont positives ou nulles, seules les dernières étant strictement négatives. Les différentes matrices inter- venant dans l'évaluation des $\ell^i$ ont alors pour expression:

$$[C] = \begin{bmatrix} C_1^1 & C_2^1 \\ C_1^2 & C_2^2 \end{bmatrix} \qquad (4.148)$$

avec:
$$ [C_1^1] = [C_2^2] = \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 & \frac{a}{2\sqrt{g\rho}} & -\frac{a}{2\sqrt{g\rho}} \end{bmatrix} \quad (4.149) $$

$[C_1^1]$ (resp. $[C_2^2]$) est remplie de zéros, exceptée la dernière ligne qui résulte du produit de $\frac{\partial V^i}{\partial(\sqrt{gq})}$ par $R^j$ (resp. de $\frac{\partial V^j}{\partial(\sqrt{gq})}$ par $R^i$).

$$ [P^1] = \begin{bmatrix} P_1^{11} & P_1^{12} \\ P_1^{12} & P_2^{12} \end{bmatrix} \quad (4.150) $$

avec:

$$ [P_1^{11}] = \begin{bmatrix} L_1^i \\ L_2^i \\ L_3^i \\ L_4^i \\ 0 \end{bmatrix} \quad (4.151) \qquad [P_2^{12}] = \begin{bmatrix} L_1^j \\ L_2^j \\ L_3^j \\ L_4^j \\ 0 \end{bmatrix} \quad (4.152) $$

$$ [P_2^{11}] = [P_1^{12}] = 0 \quad (4.153) $$

$$ [P^2] = \begin{bmatrix} P^{21} \\ P^{22} \end{bmatrix} \quad (4.154) $$

avec:

$$ [P^{21}] = \begin{bmatrix} 0 \\ 0 \\ 0 \\ 0 \\ \frac{\partial V^i}{\partial(\sqrt{gq})} \end{bmatrix} \quad (4.155) \qquad [P^{22}] = \begin{bmatrix} 0 \\ 0 \\ 0 \\ 0 \\ \frac{\partial V^j}{\partial(\sqrt{gq})} \end{bmatrix} \quad (4.156) $$

$$ [P^3] = \begin{bmatrix} P^{31} \\ P^{32} \end{bmatrix} \quad (4.157) $$

avec:

$$ [P^{31}] = \begin{bmatrix} 0 \\ 0 \\ 0 \\ 0 \\ \frac{1}{\Delta_t} V^i \end{bmatrix} \quad (4.158) \qquad [P^{32}] = \begin{bmatrix} 0 \\ 0 \\ 0 \\ 0 \\ \frac{1}{\Delta_t} V^j \end{bmatrix} \quad (4.159) $$
<figure><img src="image_68.png" /><figcaption>Figure 4.32: Coin glissant/glissant/sortie subsonique</figcaption></figure>

### Coin glissant/glissant/sortie subsonique

* **Définition**

Cette configuration résulte de l'intersection de deux frontières glissantes adiabatiques, et d'une frontière libre à travers laquelle le fluide quitte le domaine en régime subsonique.

* **Exemples d'utilisation**

Elle permet de simuler le raccord entre deux parois solides guidant le fluide et une surface de sortie (figure 4.32). Elle se rencontre en particulier à l'extrémité d'un canal de section quadrangulaire.

* **Nombre et nature des conditions physiques à imposer**

Nous sommes en présence de trois frontières (n = 3) telles que la vitesse relative du fluide est tangente à deux d'entre elles, et que le régime d'écoulement à travers la troisième est subsonique sortant. Une condition physique doit donc être imposée dans chaque direction. Soient $\vec{a}^i$, $\vec{a}^j$ et $\vec{a}^k$ les normales aux frontières considérées, et $p_k^*$ le niveau de pression statique de référence en sortie. Nous avons:

$$V^i = (\vec{V} - \vec{W}) \cdot \vec{a}^i = 0 \qquad (4.160)$$

$$V^j = (\vec{V} - \vec{W}) \cdot \vec{a}^j = 0 \qquad (4.161)$$

$$p - p_k^* = 0 \qquad (4.162)$$
• Expression des $\ell^i$

Les quatre premières valeurs propres de $K^i$, de $K^j$ et de $K^k$ sont positives ou nulles, seules les dernières étant strictement négatives. Les différentes matrices intervenant dans l'évaluation des $\ell^i$ ont alors pour expression:

$$ [C] = \begin{bmatrix} C_1^1 & C_2^1 & C_3^1 \\ C_1^2 & C_2^2 & C_3^2 \\ C_1^3 & C_2^3 & C_3^3 \end{bmatrix} \qquad (4.163) $$

avec:

$$ [C_1^1] = [C_2^2] = \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 & \frac{a}{2\sqrt{g\rho}} & -\frac{a}{2\sqrt{g\rho}} \end{bmatrix} \qquad (4.164) $$

$$ [C_3^3] = \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 & \frac{a^2}{2\sqrt{g}} & \frac{a^2}{2\sqrt{g}} \end{bmatrix} \qquad (4.165) $$

$[C_2^2]$ (resp. $[C_1^2]$) est remplie de zéros, exceptée la dernière ligne qui résulte du produit de $\frac{\partial V^i}{\partial(\sqrt{gq})}$ par $R^j$ (resp. de $\frac{\partial V^j}{\partial(\sqrt{gq})}$ par $R^i$).

$[C_3^1]$ (resp. $[C_2^1]$) est remplie de zéros, exceptée la dernière ligne qui résulte du produit de $\frac{\partial V^i}{\partial(\sqrt{gq})}$ par $R^k$ (resp. de $\frac{\partial V^j}{\partial(\sqrt{gq})}$ par $R^k$).

$[C_1^3]$ (resp. $[C_2^3]$) est remplie de zéros, exceptée la dernière ligne qui résulte du produit de $\frac{\partial p^i}{\partial(\sqrt{gq})}$ par $R^j$ (resp. de $\frac{\partial p^j}{\partial(\sqrt{gq})}$ par $R^i$).

$$ [P^1] = \begin{bmatrix} P_1^{11} & P_2^{11} & P_3^{11} \\ P_1^{12} & P_2^{12} & P_3^{12} \\ P_1^{13} & P_2^{13} & P_3^{13} \end{bmatrix} \qquad (4.166) $$

avec:

$$ [P_1^{11}] = \begin{bmatrix} L_1^i \\ L_2^i \\ L_3^i \\ L_4^i \\ 0 \end{bmatrix} \qquad [P_2^{12}] = \begin{bmatrix} L_1^j \\ L_2^j \\ L_3^j \\ L_4^j \\ 0 \end{bmatrix} \qquad [P_3^{13}] = \begin{bmatrix} L_1^k \\ L_2^k \\ L_3^k \\ L_4^k \\ 0 \end{bmatrix} \qquad (4.167) $$
$$ [P_2^{11}] = [P_3^{11}] = [P_1^{12}] = [P_3^{12}] = [P_1^{13}] = [P_2^{13}] = 0 \quad (4.170) $$

$$ [P^2] = \begin{bmatrix} P^{21} \\ P^{22} \\ P^{23} \end{bmatrix} \quad (4.171) $$

avec:

$$ [P^{21}] = \begin{bmatrix} 0 \\ 0 \\ 0 \\ 0 \\ \frac{\partial V^i}{\partial(\sqrt{g})} \end{bmatrix} \quad (4.172) $$ $$ [P^{22}] = \begin{bmatrix} 0 \\ 0 \\ 0 \\ 0 \\ \frac{\partial V^j}{\partial(\sqrt{g})} \end{bmatrix} \quad (4.173) $$ $$ [P^{23}] = \begin{bmatrix} 0 \\ 0 \\ 0 \\ 0 \\ \frac{\partial p}{\partial(\sqrt{g})} \end{bmatrix} \quad (4.174) $$

$$ [P^3] = \begin{bmatrix} P^{31} \\ P^{32} \\ P^{33} \end{bmatrix} \quad (4.175) $$

avec:

$$ [P^{31}] = \begin{bmatrix} 0 \\ 0 \\ 0 \\ 0 \\ \frac{1}{\Delta t} V^i \end{bmatrix} \quad (4.176) $$ $$ [P^{32}] = \begin{bmatrix} 0 \\ 0 \\ 0 \\ 0 \\ \frac{1}{\Delta t} V^j \end{bmatrix} \quad (4.177) $$ $$ [P^{33}] = \begin{bmatrix} 0 \\ 0 \\ 0 \\ 0 \\ \frac{1}{\Delta t} (p - p_k^*) \end{bmatrix} \quad (4.178) $$

### 4.4.4 Conditions de non réflexion

Une alternative au traitement fort ou faible des frontières libres réside dans les conditions dites de non réflexion. Elles s'appliquent principalement aux surfaces limites en lesquelles les contraintes physiques sont mal connues, telle que la distribution de pression statique en sortie subsonique d'une roue, ou variables dans le temps de manière incontrôlée, comme par exemple un sillage, un choc, ou toute autre perturbation convectée vers l'extérieur du domaine. Cette technique permet d'évacuer ces fluctuations à moindre mal.

Elle repose sur l'idée qu'il faut empêcher de rentrer les ondes inconnues en provenance du reste de l'univers qui sont susceptibles de perturber le champ aérodynamique calculé. Pour cela, toute information voyageant à vitesse négative est brutalement transformée en une onde stationnaire, en annulant sa vitesse. Nous parlons de condition de non réflexion car une onde interne se réfléchissant en sortie génère
un retour d'informations. N'autorisant aucune onde entrante, nous supprimons du même coup ce phénomène.

Ce traitement est un peu particulier dans la mesure où il ne fait pas intervenir de contraintes extérieures physiques. Il existe plusieurs formalismes pour le mettre en place [64, 27, pp. 384]. Celui présenté pour les conditions aux limites faibles s'y prête bien. Nous suivons la même approche, seul le traitement des ondes entrantes étant modifié.

Pour chaque direction $\xi^i$, le signe des valeurs propres $\Lambda^i$ nous renseigne sur le sens de propagation des informations le long de la normale $\bar{a}^i$. Soit $\lambda_k^i$, la k-ième valeur propre de $K^i$. Notre but est de déterminer l'expression de $\ell_k^i$ qui y est associé, pour la réintroduire dans l'équation 4.115.

* Si $\lambda_k^i$ est positive:

L'onde de vitesse $\lambda_k^i$ sort du domaine de calcul, et elle est de ce fait connue. Nous avons:

$$ \ell_k^i = L_k^i \frac{\partial F_c^i}{\partial \xi^i} \qquad (4.179) $$

où $\frac{\partial F_c^i}{\partial \xi^i}$ est à calculer de manière décentrée depuis l'intérieur du domaine. Comme précédemment:

$$ C_{jl}^{ik} = \delta_{ij} \delta_{kl} \qquad (4.180) $$

$$ P_j^{1ik} = \delta_{ij} L_k^i \qquad (4.181) $$

$$ P_l^{2ik} = 0 \qquad (4.182) $$

$$ P^{3ik} = 0 \qquad (4.183) $$

* Si $\lambda_k^i$ est négative:

L'onde associée apporte l'influence du reste de l'univers. C'est pourquoi elle est remplacée par une onde stationnaire, de vitesse $\lambda_k^i$ nulle. Si bien que:

$$ \ell_k^i = \lambda_k^i \frac{\partial}{\partial \xi^i} (\sqrt{g} w_k^i) \qquad (4.184) $$

$$ = 0 \qquad (4.185) $$

Et donc:

$$ C_{jl}^{ik} = \delta_{ij} \delta_{kl} \qquad (4.186) $$

$$ P_j^{1ik} = 0 \qquad (4.187) $$

$$ P_l^{2ik} = 0 \qquad (4.188) $$

$$ P^{3ik} = 0 \qquad (4.189) $$
<figure><img src="image_69.png" /><figcaption>Figure 4.33: Interaction d'un choc et d'une frontière de non reflexion</figcaption></figure>

C'est le cas le plus simple du formalisme développé car il ne fait pas intervenir de forme linéarisée associée à une condition physique. Il est à noter qu'il s'intègre parfaitement dans le traitement faible.

Cette méthode permet de répondre à plusieurs problèmes. Tout d'abord, elle est souvent indispensable lors d'une simulation instationnaire pour évacuer des perturbations générées au sein de l'écoulement, par exemple par le mouvement d'une pale, ou imposées sur une autre frontière, telles un sillage défilant en entrée d'un canal inter-aubes. D'autre part, utilisée en début de calcul, elle laisse sortir des perturbations parasites engendrées par le champ d'initialisation; et imposée en fin de convergence lorsque l'écoulement est établi, elle permet de stabiliser la solution si des ondes numériques se déplacent encore dans le domaine de calcul.

A titre d'illustration est présenté l'évolution au cours du temps du niveau de masse volumique dans un tube à choc dont les extrémités sont traitées soit par une condition de non reflexion (figure 4.33); soit par une condition de frontière libre (figure 4.34). Nous noterons dans le premier cas que le choc sort du domaine de calcul sans interaction notable avec la frontière, alors qu'il se transforme en une onde de détente dans le second.

### 4.4.5 Avantages et inconvénients

Cette approche présente les mêmes avantages que le traitement fort, en en corrigeant certains défauts. Tout d'abord, elle est générale, permettant de gérer aussi bien les faces, que les arêtes et les coins du domaine de calcul. De plus, elle est "commutative", au sens où toutes les directions sont équivalentes, si bien que l'ordre dans lequel elles sont traitées est sans importance. Enfin, elle assure elle aussi un contrôle
<figure><img src="image_70.png" /><figcaption>Figure 4.34: Interaction d'un choc et d'une frontière libre</figcaption></figure>

simple et précis des conditions aux limites au cours du temps.

Bien qu'elle est fait l'objet de moins de tests que l'imposition forte, cette technique apparaît comme viable et très prometteuse. Elle devrait même à terme supplanter sa rivale, bien qu'elle soit d'un coût CPU légèrement supérieur, l'implémentation de l'une et de l'autre étant peu différentes.

Le principal obstacle à son utilisation est le choix critique des conditions aux limites physiques devant être imposées à l'intersection de plusieurs frontières. En effet, elle est inutilisable si la matrice [C] se révèle non inversible. Mais ce problème n'est pas propre à cette approche, et sa résolution n'est sans doute pas liée au cadre dans lequel nous nous plaçons.

## 4.5 Résolution numérique

De la même manière que nous l'avons fait pour les noeuds du coeur du domaine de calcul, il nous faut maintenant évaluer numériquement le champ aérodynamique sur les plans supplémentaires et les frontières. Cette opération doit en particulier viser à ne pas introduire de déphasages éventuels entre ces différentes régions.

### 4.5.1 Gestion des noeuds supplémentaires

Il ressort du traitement des conditions aux limites géométriques qu'un noeud supplémentaire est lié à ceux du domaine de calcul par une relation linéaire ayant pour forme générique:

$$M_{ext} \cdot q^{ext} + M_{int} \cdot q^{int} + M_c \cdot q^c = 0 \quad (4.190)$$
où $M_{ext}$, $M_{int}$ et $M_c$ sont des matrices 5×5 déduites des équations 4.10, 4.15, 4.20 et 4.25.

Dans une approche Delta, nous écrivons que:

$$q^{(n+1)} = q^{(n)} + \Delta q^n \quad (4.191)$$

Si l'équation 4.190 est vérifiée à l'instant $(n+1)$, il en découle que:

$$M_{ext}^{(n+1)} \cdot \Delta q^{ext \ n} + M_{int}^{(n+1)} \cdot \Delta q^{int \ n} + M_c^{(n+1)} \cdot \Delta q^c \ n = \\ - [M_{ext}^{(n+1)} q^{ext \ (n)} + M_{int}^{(n+1)} q^{int \ (n)} + M_c^{(n+1)} q^c \ (n)] \quad (4.192)$$

Ainsi, la gestion des noeuds supplémentaires peut-elle être totalement implicite ou explicite. Dans ce dernier cas, la solution est évaluée dans un premier temps dans le domaine utile, produisant les valeurs de $\Delta q^{int \ n}$ et de $\Delta q^c \ n$, puis mise à jour sur les noeuds externes. Il est à noter qu'un multi-domaine ne pourra être gérer en instationnaire que d'une manière explicite, car l'approche implicite nécessiterait la résolution simultanée de tous les domaines à cause des noeuds communs, ce qui va à l'encontre d'une fragmentation du problème initial.

### 4.5.2 Gestion des noeuds frontières

L'équation 4.74, ou 4.139, semi-discrète obtenue en un noeud de la frontière est formellement proche en temps de celle résolue en un noeud du coeur (équation 3.26). Elle est donc discrétisée selon le même schéma, qui pour cette étude, est celui d'Euler explicite (paragraphe 3.3). Ainsi, la résolution numérique en un noeud frontière sera en tout point semblable à celle en un noeud interne (paragraphe 3.7).

Une réserve à l'encontre des conditions de périodicité doit être cependant émise. En effet, dans la méthode précédemment exposée, ce cas est traité de manière transparente au niveau des noeuds frontières puisque leur voisinage reconstruisant exactement cette situation, aucune autre disposition n'est théoriquement nécessaire. Cependant, au raccord d'une condition géométrique et d'une condition physique forte ou faible, peut apparaître un problème. Par exemple en un bord d'attaque (figure 4.35), la ligne de maillage frontière, périodique jusqu'alors, se divise en deux branches, une sur l'extrados, l'autre sur l'intrados. Le voisinage étendu assure que les flux perpendiculaires à cette ligne sont identiques. Par contre, les flux tangents diffèrent selon la branche suivie. Il en résulte un écart entre les solutions calculées au point de séparation. Cette erreur tend au cours du temps à se propager le long de la ligne, brisant la condition de périodicité. Un remède à ce problème est d'utiliser une condition dite de périodicité forte. Le champ aérodynamique n'est alors calculé de manière classique que sur une branche, puis recopier sur l'autre frontière, en laquelle les équations de Navier-Stokes ne sont plus résolues. Cependant, sur un maillage suffisamment fin, l'erreur commise reste faible et localisée.
<figure><img src="image_71.png" /><figcaption>Figure 4.35: Raccord de différentes conditions en un bord d'attaque</figcaption></figure>

## 4.6 Conclusion

Ce chapitre a été consacré au traitement des conditions aux limites. Tout d'abord, nous nous sommes attachés à la gestion des conditions géométriques, qui permettent entre autres de simuler un plan de symétrie, par extension du somme de calcul et création de noeuds supplémentaires. Les avantages d'une telle approche ont été mis en lumière, en particulier en terme de simplicité, mais aussi d'extension de la méthode initiallement développée sur un maillage en H, à des domaines en O ou en C, voir à des multi-domaines. Puis, un premier traitement instationnaire des conditions physiques, fondé sur les relations de compatibilité mono-dimensionnelles, a été présenté. Il s'est révélé de bonne qualité, mais peu apte à gérer les arêtes et les coins. C'est pourquoi une deuxième formulation, bâtie sur les relations de compatibilité multi-dimensionnelles cette fois, a été développée. Moins testée que la première, elle apparaît cependant prometteuse.

Ainsi s'achève la résolution numérique des équations d'évolution, dans le coeur et aux frontières du domaine de calcul. Le chapitre suivant s'attachera à décrire les mouvements éventuels du maillage.
# Chapitre 5

# Mouvements du maillage

## 5.1 Introduction

Ce chapitre s'intéresse aux possibilités offertes par la mise en oeuvre d'un maillage structuré mobile et déformable en temps. A ce titre, nous ne prétendons pas couvrir par la suite tous les aspects liés à cette spécificité, mais plutôt présenter quelques voies de développement jugées originales et prometteuses.

Une des premières utilisations, citée lors de l'intégration en espace, en est la gestion d'un mouvement de rotation. Si la vitesse angulaire est constante, le changement de variables spatiales de cartésiennes à curvillignes correspond alors exactement à la réécriture classique des équations différentielles initiales dans un repère cylindrique tournant [28, pp.16-21]. Mais un avantage indéniable de l'approche suivie est qu'elle ne suppose pas une vitesse de rotation uniforme, ni un déplacement en bloc des noeuds. En effet, tout autre mouvement, lié à la déformation du canal interaubes par exemple, pourra se superposer à la rotation sans que le formalisme ou les équations doivent être modifiés.

Un autre usage de cette faculté de déplacer les noeuds est fait lors de l'adaptation du maillage. Cette opération consiste à modifier la répartition de la taille des celules d'intégration afin de l'adapter au champ aérodynamique. Si ceci est intéressant dans les régions à fortes évolutions pour y augmenter la densité des noeuds, c'est par contre une quasi-nécessité au niveau des discontinuités. En effet, l'action des limiteurs dans le voisinage de ces dernières, dégrade au premier ordre le schéma de discrétisation spatiale. Pour assurer une précision assez uniforme de la solution, il faut donc que dans ces zones, les mailles soient relativement de très petites tailles. Dans la mesure où le nombre de noeuds reste constant, il n'y a pas création de nouveaux éléments, contrairement à une approche non structurée, et l'adaptation du maillage est faite par déplacement des points. Nous sommes donc bien dans le cadre d'un maillage mobile et déformable. Cette opération peut se réaliser en post-traitement, c'est à dire qu'à partir d'un champ aérodynamique un nouveau maillage est construit. Nous parlerons alors d'adaptation *statique*. Une autre approche consiste à coupler en temps le calcul de la solution et la déformation du maillage. Il s'agit
<figure><img src="image_72.png" /><figcaption>Figure 5.1: Adaptation statique. Solution initiale</figcaption></figure>

alors d'une adaptation dynamique. Ces deux techniques font l'objet des paragraphes suivants.

## 5.2 Adaptation statique

Les développements présentés par la suite sont inspirés des travaux de Fraga et Morris [22], dans lesquels une technique d'adaptation de maillage est mise en oeuvre pour capturer des solitons. Ce sont des ondes mobiles de forte amplitude et de petite longueur d'onde, qui ont la faculté de se croiser sans interagir. Les écoulements que nous simulons sont loin de cette description, comme le montre par exemple le profil de masse volumique extrait d'un tube à choc alors que le faisceau d'ondes est pleinement développé (figure 5.1). Cependant, si nous nous intéressons au limiteur associé à cette variable, tracé sur la même figure, nous voyons apparaître quatre ondes, d'amplitude unitaire et de faible longueur d'onde, correspondant à chacune des discontinuités de la solution. Ce seront nos "solitons". Le but de cette partie est d'exposer le cheminement qui va nous amener de la répartition uniforme de points visible figure 5.1, à celle mieux adaptée présentée figure 5.6. Lors de cette transformation, le champ aérodynamique est figé.

Tout d'abord, est déterminé le limiteur (équation 3.53) associé à chaque segment (i), compris entre les noeuds (i-1) et (i) (figure 5.1):
<figure><img src="image_73.png" /><figcaption>Figure 5.2: Adaptation statique. Limiteur après seuillage</figcaption></figure>

$$ \Phi_i^b = \min[\Phi(q_{i-1}, q_i, q_{i+1}), \Phi(q_i, q_{i-1}, q_{i-2})] \quad (5.1) $$

Puis, la fonction $\Phi^b$ étant continue entre 0.0 et 1.0, un seuillage la transforme en une fonction discontinue $\Phi^s$ telle que (figure 5.2):

$$ \Phi_i^s = 0.0 \quad \text{si} \quad \Phi_i^b \le 0.5 \quad (5.2) $$

$$ \Phi_i^s = 1.0 \quad \text{si} \quad \Phi_i^b > 0.5 \quad (5.3) $$

Ensuite, une opération de fenêtrage permet de créer une transition linéaire entre les deux extrema de $\Phi^s$, conduisant à la fonction $\Phi^f$ (figure 5.3). Un lissage est alors effectué, tel que (figure 5.4):

$$ \Phi_i^f = \frac{1}{4} [\Phi_{i-1}^f + 2\Phi_i^f + \Phi_{i+1}^f] \quad (5.4) $$

Toutes ces transformations ont eu pour but d'obtenir une fonction continue entre 0.0 et 1.0, qui soit la plus régulière possible, et qui soit minimale dans les régions où les plus petites mailles sont nécessaires. Elle est alors convertie en une information qui traduit le rapport $R$ entre le plus grand et le plus petit segment:

$$ \Phi_i^c = R + (1 - R)\Phi_i^f \quad (5.5) $$

Cette fonction est alors intégrée:
<figure><img src="image_74.png" /><figcaption>Figure 5.3: Adaptation statique. Limiteur après fenêtrage</figcaption></figure>

<figure><img src="image_75.png" /><figcaption>Figure 5.4: Adaptation statique. Limiteur après lissage</figcaption></figure>
<figure><img src="image_76.png" /><figcaption>Figure 5.5: Adaptation statique. Limiteur après intégration</figcaption></figure>

$$ \Phi_i^t = \Phi_{i-1}^t + \Phi_i^c[x_i - x_{i-1}] \quad (5.6) $$

puis normalisée, N étant le nombre de noeuds (figure 5.5):

$$ \Phi_i^n = \frac{N}{\Phi_i^c} \Phi_i^t \quad (5.7) $$

Cette dernière fonction représente la nouvelle répartition des points. Plus sa pente est importante, plus deux noeuds successifs sont proches. Le rapport entre sa pente la plus forte et sa pente la plus faible, fixe le rapport entre le segment le plus long et le segment le plus court. Le fenêtrage et le lissage assurent une évolution régulière entre les différentes régions.

Il ne reste maintenant plus qu'à lire l'indice des noeuds sur l'axe des ordonnées de la courbe 5.5, pour en déduire sur l'axe des abscisses leur position géométrique. Dans ce but, une simple interpolation linéaire conduit alors à la figure 5.6. Le resserrement des points de façon symétrique autour des discontinuités est à noter. Par ailleurs, une pondération par le gradient de masse volumique de la fonction $\Phi^s$ permettrait d'accroître la densité des noeuds dans l'onde de détente.

Cette technique, illustrée jusqu'à présent dans une configuration mono-dimensionnelle, a été appliquée avec succès à un écoulement transsonique faiblement bidimensionnel. Figure 5.7 est présenté le niveau de masse volumique, dans la région du choc, calculé sur un maillage bâti à partir d'une première solution obtenue avec une répartition régulière des noeuds. Un rapport de dix entre la plus grosse et la plus
<figure><img src="image_77.png" /><figcaption>Figure 5.6: Adaptation statique. Solution finale</figcaption></figure>

<figure><img src="image_78.png" /><figcaption>Figure 5.7: Adaptation statique. Choc bi-dimensionnel</figcaption></figure>
<figure><img src="image_79.png" /><figcaption>Figure 5.8: Adaptation dynamique</figcaption></figure>

petite maille a été choisi. Il en résulte une capture très fine du choc, avec une densité de noeuds tout à fait satisfaisante à la vue de l'évolution de la masse volumique.

## 5.3 Adaptation dynamique

Nous avons par ailleurs couplé la technique exposée jusqu'à présent avec le calcul de la solution analytique spatio-temporelle de l'écoulement dans un tube à choc. Il ne s'agit donc pas d'un couplage réel avec le reste du code de calcul développé, mais cette approche permet d'illustrer notre capacité à adapter dynamiquement le maillage à des discontinuités mobiles. La répartition des noeuds dans l'espace au cours du temps est visible figure 5.8. Les quatre discontinuités associées au choc, à la surface de contact et à l'onde de détente y sont clairement visibles. Un choix judicieux de la largeur de fenêtrage, ainsi que le lissage, assurent qu'à chaque pas de temps, les ondes se déplacent dans une zone finement maillée. En effet dans ce cas, il ne suffit pas de raffiner très locallement, il faut aussi prendre en compte le fait que les structures à calculer se déplacent. Le résultat obtenu est très encourageant. Cependant, il présente quelques défauts. Tout d'abord, au début de la simulation, les ondes sont très proches les unes des autres, et l'opération de seuillage/fenêtrage/lissage semble avoir quelques difficultés à les différencier. Il en résulte des déplacements importants des noeuds d'un instant à l'autre, donnant cet aspect chahuté. Puis une évolution monotone semble s'établir. Elle est cependant parfois perturbée par un brusque cisaillement qui est dû à une maille qui voit sa taille passée d'un extrème à l'autre.
Cette dilatation, ou cette contraction, influence tous les noeuds. Une évolution trop brutale de la fonction Φ" en est sans doute la cause.

## 5.4 Conclusion

Ce chapitre a été consacré à l'exploration de quelques voies de développement liées à la capacité du code de calcul présenté à gérer un maillage mobile et déformable en temps. Cette particularité peut être mise à profit pour prendre en compte les effets de la rotation du domaine de calcul, sans qu'une hypothèse restrictive de vitesse angulaire uniforme ou de mouvement en bloc ne soit nécessaire. D'autre part, cette faculté de déplacer les noeuds a été utilisée pour adapter le maillage à une solution stationnaire ou instationnaire.

Jusqu'à présent, aucun couplage entre le calcul du champ aérodynamique et le mouvement des noeuds n'a été réalisé. Cependant, chaque phase de développement de la méthode exposée, de la discrétisation en espace au traitement des conditions aux limites, intègre cette possibilité. Aussi aucune difficulté majeure ne devrait entraver cette opération. Dans un autre domaine, les critères de raffinement de maillage méritent sans doute une étude supplémentaire, mais ceux déjà mis en oeuvre ouvrent une voie prometteuse.
# Chapitre 6

# Simulation et analyse
d'écoulements compressibles

## 6.1 Introduction

Ce chapitre regroupe les résultats obtenus sur des configurations plus complexes que celles illustrant le développement et la validation de la méthode présentée jusqu'ici. Des écoulements en géométries réelles, ou en tous cas proches de la réalité des turbomachines, sont ainsi explorés. Nous nous attacherons en particulier à analyser d'un point de vue physique et numérique les solutions obtenues comme un ensemble, et non plus en tant que cas test de validation, ciblés sur un problème numérique spécifique.

Les sept configurations présentées par la suite peuvent être regroupées en trois familles, dont la complexité croît avec le réalisme.

La première s'attache à la propagation et aux interactions d'ondes instationnaires dans un fluide parfait. L'écoulement simulé le plus simple est un tube à choc isolé, mettant en évidence la capacité de la méthode développée à gérer des discontinuités mobiles mono-dimensionnelles (paragraphe 6.2). Puis nous nous intéresserons au paragraphe 6.3 aux interactions d'ondes qui peuvent se produire dans une telle configuration. Enfin, la simulation de discontinuités bi-dimensionnelles et de leurs réflexions fera l'objet du paragraphe 6.4.

La seconde famille de résultats vise à prédire des écoulements confinés entre deux ou plusieurs parois solides. Tout d'abord, une tuyère supersonique sera étudiée. Un fluide parfait y passe d'un régime subsonique à un régime supersonique, sans discontinuité du champ aérodynamique stationnaires (paragraphe 6.5). Puis un canal transsonique, recelant un choc, sera présenté. Des écoulements bi-dimensionnels, en fluide parfait ou visqueux, stationnaires ou instationnaires, y ont été calculés (paragraphe 6.6).

Enfin, deux cas stationnaires proches de la réalité seront exposés. Le premier met en oeuvre une grille d'aubes supersonique, dont la structure de chocs au bord
<figure><img src="image_80.png" /><figcaption>Figure 6.1: Tube à choc. Faisceau d'ondes</figcaption></figure>

de fuite est caractéristique (paragraphe 6.7), et le second un coude tri-dimensionnel à bas nombre de Mach, semblable à ceux présents entre deux roues d'une machine radiale (paragraphe 6.8).

## 6.2 Propagation d'ondes mono-dimensionnelles

Rares sont les configurations qui ont été autant étudiées expérimentalement et numériquement que l'écoulement dans un tube à choc [52]. Il présente la particularité d'être mono-dimensionnel, si l'influence des couches limites latérales est négligée, ce qui est souvent un gage de simplicité. Or dans ce cas, sa nature instationnaire et la présence de plusieurs discontinuités mobiles à des vitesses différentes, le rendent particulièrement complexe. Un autre paradoxe de ce problème est qu'il admet une solution analytique simple, alors qu'il est ardu à simuler numériquement. Il constitue donc une configuration de choix, à la fois bien documentée, simple d'apparence et très sensible à tout défaut de modélisation.

Un tube à choc peut être décrit comme un long tube horizontal, initiallement divisé en deux parties par un diaphragme (figure 6.1). Du côté gauche se trouve un gaz au repos à haute pression, et du côté droit le même gaz à une pression inférieure. La température est uniforme et identique dans les deux chambres, et le fluide est supposé parfait. Lorsque le diaphragme est rompu, une onde de pression se propage vers la droite dans la chambre basse pression, suivie par une surface de contact. Une ondes de détente remonte pendant ce temps dans la chambre haute pression. La discontinuité de contact, qui influence la valeur de la masse volumique, mais non les niveaux de pression et de vitesse, sépare les deux régions fluides ini-
<figure><img src="image_81.png" /><figcaption>Figure 6.2: Tube à choc. Niveau de masse volumique. Premier ordre</figcaption></figure>

tiales, qui ne se mélangent pas au cours du temps car il n'y a pas d'effets de diffusion.

Notre but à travers un tel exemple est d'estimer la capacité de la méthode présentée à suivre des discontinuités instationnaires en l'absence de toute interaction entre des ondes, ou avec des surfaces solides, et ce pour différents schémas en espace. La solution exacte de ce problème est connue, et utilisée comme référence [60, pp. 1007-1011], ainsi que celle obtenue en maillage non structuré avec la technique de Roe.

Les conditions initiales choisies sont telles qu'un rapport de pression statique et de densité de 20 existe entre les deux chambres (annexe E.1). D'autre part, la vitesse est nulle, et la température est telle que le rapport de la pression sur la masse volumique soit de 1. La longueur du tube est de 4.5, et il est divisé en 100 cellules régulières. La position du diaphragme est à $X=1.95$. Le pas en temps est constant en tout noeud et à chaque instant, et est fixé à $2.10^{-3}$ s. Dans ces conditions, 10 itérations sont nécessaires à l'onde de pression, qui est l'information physique la plus rapide, pour aller d'un point au suivant. Une grande précision en temps est ainsi obtenue. L'écoulement est transsonique, le nombre de Mach variant entre 0.00 et 1.36.

L'évolution au cours du temps des niveaux de masse volumique et de pression statique est présentée figures 6.2 et 6.3 pour un calcul au premier ordre en espace, et figures 6.4 et 6.5 pour un schéma au troisième ordre. Sur chacune, le faisceau
<figure><img src="image_82.png" /><figcaption>Figure 6.3: Tube à choc. Niveau de pression statique. Premier ordre</figcaption></figure>

<figure><img src="image_83.png" /><figcaption>Figure 6.4: Tube à choc. Niveau de masse volumique. Troisième ordre</figcaption></figure>
<figure><img src="image_84.png" /><figcaption>Figure 6.5: Tube à choc. Niveau de pression statique. Troisième ordre</figcaption></figure>

d'ondes prédit par la théorie est clairement visible. Le choc apparaît comme une marche, la plus à droite, ainsi que la surface de contact qui n'est pas visible sur les profils de pression statique. L'onde de détente forme un éventail s'ouvrant sur la gauche comme le temps croît, entre deux fronts rectilignes. La nature fortement dissipative du schéma au premier ordre est responsable du lissage des ondes figures 6.2 et 6.3, dissipation qui tend même à les faire disparaître. Par contre, leur amplitude et leur netteté se conservent très bien lorsque l'ordre en espace est augmenté (figures 6.4 et 6.5), ce qui traduit une faible dissipation numérique dans ce cas.

Les profils de densité et d'énergie totale à $t = 1s$ sont tracés figures 6.6 à 6.9. Tout d'abord, la similitude entre les deux solutions numériques est frappante.

Comme noté précédemment, le calcul au premier ordre lisse les discontinuités, à un point tel que la surface de contact disparaît presque sur les profils d'énergie (figure 6.7). Le saut visible au milieu de l'onde de détente ($X \approx 2.0$) est la trace de la position initiale du diaphragme (figures 6.6 et 6.7). Le schéma au premier ordre maintient cette discontinuité parasite située en un point immobile où le nombre de Mach est de 1. C'est là que le raccord entre les régions subsonique et supersonique s'opère mal. Ce défaut disparaît lorsque l'ordre en espace augmente, car l'onde de détente étant un phénomène continu, l'interpolation en ce point n'est plus d'ordre un.

Pour les schémas d'ordre supérieur, les résultats sont en accord avec la théorie. Les fronts de l'onde de détente sont juste un peu arrondis, pendant que l'onde de pression est étalée sur 2-3 mailles, et la surface de contact sur 7-8 cellules. A cause
<figure><img src="image_85.png" /><figcaption>Figure 6.6: Tube à choc. Profil de masse volumique. Premier ordre</figcaption></figure>

<figure><img src="image_86.png" /><figcaption>Figure 6.7: Tube à choc. Profil d'énergie totale. Premier ordre</figcaption></figure>
<figure><img src="image_87.png" /><figcaption>Figure 6.8: Tube à choc. Profil de masse volumique. Troisième ordre</figcaption></figure>

<figure><img src="image_88.png" /><figcaption>Figure 6.9: Tube à choc. Profil d'énergie totale. Troisième ordre</figcaption></figure>
<figure><img src="image_89.png" /><figcaption>Figure 6.10: Interactions d'ondes 1D. Faisceau d'ondes</figcaption></figure>

des limiteurs qui ramènent au premier ordre le schéma au niveau des discontinuités, il n'y a pas d'oscillations spatiales de la solution, mais une douce évolution. Dans ce cas, la surface de contact est préservée sur le profil d'énergie (figure 6.9).

A la vue de ces résultats, nous pouvons prétendre à une bonne résolution d'ondes instationnaires d'amplitudes disparates par la méthode développée, lorsqu'un schéma d'ordre élevé en espace est mis en oeuvre. Une remarque intéressante est que seules des différences peu significatives distinguent les solutions obtenues avec la séparation de flux de Van Leer, et la technique de Roe.

## 6.3 Interactions d'ondes mono-dimensionnelles

La configuration suivante met en oeuvre deux tubes à choc, partageant la même chambre basse pression (figure 6.10). A l'instant initial, les rapports de pression statique et de densité sont les mêmes à travers les deux diaphragmes. La température est uniforme, et le gaz au repos. Juste après la rupture des diaphragmes, le même faisceau d'ondes que celui décrit dans le cadre d'un seul tube à choc est visible. Les ondes de pression et les surfaces de contact se propagent dans la chambre basse pression, vers le milieu du tube, pendant que les ondes de détente se déplacent vers ses extrémités, dans les chambres haute pression. A cause des conditions initiales qui sont symétriques, le faisceau d'ondes l'est aussi par rapport au centre du tube. Entre les instants 0.20s et 0.40s, l'interaction des deux ondes de pression génère
<figure><img src="image_90.png" /><figcaption>Figure 6.11: Interactions d'ondes 1D. Niveau de masse volumique</figcaption></figure>

deux chocs réfléchis. Puis, autour de l'instant $t=0.50s$, ces derniers rencontrent les surfaces de contact, qui suivaient les ondes de pression initiales. Un choc réfléchi et un choc transmis sont alors créés, ainsi qu'une nouvelle discontinuité de contact. Les deux ondes de pression réfléchies se déplacent vers le centre du tube, où elles interagissent à nouveau. Etc, etc, ...

Nous avons vu dans le cas du tube à choc simple que la méthode développée est susceptible de suivre des ondes mobiles sans les dissiper, lorsqu'un schéma d'ordre élevé est utilisé. Avec un tel écoulement, nous cherchons à évaluer de quelle manière elle gère les interactions d'ondes de différente nature et d'amplitude variable. La solution exacte de ce problème est connue et nous sert de référence [44, 60, pp. 1018-1028].

Les conditions initiales choisies sont telles qu'un rapport de pression statique et de densité de 4 existe à travers les diaphragmes (annexe E.2). D'autre part, la vitesse est nulle et la température telle que le rapport de la pression sur la masse volumique soit de 1. Le tube s'étend de $X = -1.0$ à $X = 1.0$, et il est divisé en 200 cellules uniformes. La position initiale des diaphragmes est $X = -0.495$ et $X^* = 0.495$. Le pas en temps est constant, et fixé à $1.10^{-3}$ s. Dans ces conditions, il faut 6 itérations aux ondes de pression initiales pour parcourir une maille. Nous utilisons le schéma du troisième ordre en espace, ainsi que des conditions de non réflexion aux extrémités du tube afin d'évacuer les ondes de détente. Nous ne nous intéressons par la suite qu'aux trois premières interactions ayant lieu avant $t=0.70s$.
<figure><img src="image_91.png" /><figcaption>Figure 6.12: Interactions d'ondes 1D. Niveau de pression statique</figcaption></figure>

Figure 6.11 (resp. figure 6.12) est dessinée l'évolution au cours du temps de la masse volumique (resp. de la pression statique). Pour $t < 0.30s$, nous reconnaissons la même structure de marche que celle mise en évidence dans le cas d'un tube à choc isolé. Puis, à chaque réflexion, nous voyons les niveaux augmenter brutalement. Il est à noter que l'interaction de l'onde de pression initiale et de la surface de contact, qui est invisible sur le niveau de pression statique, engendre bien deux chocs, un réfléchi et un transmis qui eux, apparaissent clairement figure 6.12. Nous remarquons que les fronts d'onde conservent la même structure au cours de leur périgrination, sans se dissiper.

L'énergie totale, en fonction de l'abscisse du tube et du temps, est tracée figure 6.13. Son niveau étant constant entre deux discontinuités, sauf pour l'onde de détente, les plateaux qui séparent chaque front d'onde y forment une plage de couleur uniforme. Le motif qui en résulte est à rapprocher de la figure théorique 6.10. Si nous nous intéressons au limiteur associé à cette variable, les discontinuités apparaissent encore mieux sous la forme de deux traits parallèles (figure 6.14). Tracées à la même échelle, les figures 6.10 et 6.14 se superposent parfaitement, ce qui établit que la vitesse de propagation des différentes ondes est correctement prédite.

Afin de juger de l'exactitude de leur amplitude, nous avons extrait des profils de masse volumique à différents instants, $t=0.20s, 0.40s, 0.60s$ et $0.70s$, avant et après chaque interaction. Ils sont comparés figures 6.15 à 6.18 à la solution analytique, ainsi qu'à celle produite en maillage non structuré avec la méthode de Roe. Dans la mesure où le devenir de l'onde de détente n'est que de peu d'importance dans cette
<figure><img src="image_92.png" /><figcaption>Figure 6.13: Interactions d'ondes 1D. Niveau d'énergie totale</figcaption></figure>
<figure><img src="image_93.png" /><figcaption>Figure 6.14: Interactions d'ondes 1D. Limiteur associé à l'énergie totale</figcaption></figure>
<figure><img src="image_94.png" /><figcaption>Figure 6.15: Interactions d'ondes 1D. Profil de masse volumique à $t=0.20s$</figcaption></figure>

<figure><img src="image_95.png" /><figcaption>Figure 6.16: Interactions d'ondes 1D. Profil de masse volumique à $t=0.40s$</figcaption></figure>
<figure><img src="image_96.png" /><figcaption>Figure 6.17: Interactions d'ondes 1D. Profil de masse volumique à $t=0.60s$</figcaption></figure>

<figure><img src="image_97.png" /><figcaption>Figure 6.18: Interactions d'ondes 1D. Profil de masse volumique à $t=0.70s$</figcaption></figure>
simulation, la zone d'étude a été restreinte à l'intervalle $X \in [-0.5; 0.5]$. A l'instant $t=0.20s$ (figure 6.15), nous observons des résultats semblables à ceux exposés dans le cas d'un seul tube à choc. Nous reconnaissions les ondes de pression se déplaçant vers $X = 0.0$, suivies des surfaces de contact. Pour $t$ autour de 0.30s, se produit l'interaction entre les chocs, et deux ondes de pression réfléchies se propageant vers l'extérieur du tube en découlent. Elles sont visibles figure 6.16 à $t=0.40s$ de chaque côté de $X = 0.0$, alors que les discontinuités de contact continuent à avancer en direction du centre du tube. A une date proche de $t=0.50s$, il y a interaction entre les chocs réfléchis et les surfaces de contact. Une onde de pression réfléchie et une autre transmise apparaissent, ainsi qu'une nouvelle surface de contact. A $t=0.60s$, les chocs réfléchis sont proches de $X=0.0$, pendant que les discontinuité de contact sont à $X = \pm 0.2$ (figure 6.17). Les chocs transmis se propagent vers l'extérieur. La dernière interaction a lieu autour de l'instant $t = 0.65s$ entre les deux ondes de pression réfléchies. Deux petits chocs s'éloignant du centre du tube en résultent. Ils sont visibles figure 6.18 à $t=0.70s$ de chaque côté de l'origine.

Nous noterons que tout du long du calcul, les deux méthodes numériques prédissent la même amplitude et une position semblable pour chaque onde, même après plusieurs interactions. Les solutions ainsi obtenues sont très proches des résultats analytiques. Tous les chocs restent confinés sur 3-4 mailles, et les surfaces de contact sur 7-8 cellules. La dissipation des discontinuités est faible, et n'augmente ni au cours du temps, ni lors des interactions. Même les dernières ondes de pression, qui sont de faible amplitude, sont correctement capturées comparées à la théorie.

Nous établissons avec ces résultats que la méthode développée est capable de simuler des ondes instationnaires, y compris lorsque de nombreuses interactions de nature variable de produisent.

## 6.4 Interactions d'ondes bi-dimensionnelles

La simulation d'un choc cylindrique convergent est présentée par la suite. Dans ce problème bi-dimensionnel, un diaphragme cylindrique sépare initialement deux gaz parfaits au repos (figure 6.19). La pression et la densité sont plus faibles dans la chambre centrale ainsi délimitée. Aussi, le long d'un diamètre, les profils initiaux de masse volumique, de quantité de mouvement et d'énergie sont identiques à ceux rencontrés dans la configuration précédente, mettant en oeuvre deux tubes à choc en face l'un de l'autre. Mais, après que le diaphragme ait été rompu, ce sont des ondes cylindriques, et non plus planes, qui se propagent. Cependant, le schéma d'interactions d'ondes au cours du temps est le même. Seuls les dates et les niveaux changent sous l'influence du confinement angulaire qui n'existait pas dans le cas mono-dimensionnel.

Ce problème n'admet pas de solution analytique exacte, mais il a fait l'objet
<figure><img src="image_98.png" /><figcaption>Figure 6.19: Interactions d'ondes 2D. Niveau de masse volumique à t=0.00s</figcaption></figure>

de nombreuses simulations numériques. A titre de référence, nous utilisons la plus récente publiée, à notre connaissance, par Zhang et So [67]. Une telle configuration permet d'évaluer la résolution numérique du déplacement et de l'interaction d'ondes bi-dimensionnelles de natures diverses.

Les conditions initiales choisies sont telles qu'un rapport de pression statique et de densité de 4 existe à travers le diaphragme (annexe E.3). D'autre part, la vitesse est nulle et la température telle que le rapport de la pression sur la masse volumique soit de 1. Le domaine physique est un carré d'arête 2.0, centré sur l'origine. Le rayon initial du diaphragme est de 0.495. Le maillage cartésien de calcul couvre un quart de l'espace, $(X, Y) \in [0.0; 1.0]^2$, une condition de symétrie plane étant appliquée sur chaque axe. Des conditions de non réflexion sur les faces externes permettent d'évacuer les ondes de détente, et de simuler ainsi un quart de plan infini d'un point de vue aérodynamique. Il y a 101 noeuds dans chaque direction, pour un total de 10000 cellules cubiques et uniformes. Le pas en temps est constant et fixé à $1.10^{-3}$s, comme dans le cas précédent.

Les niveaux de masse volumiques prédits par la présente méthode numérique, à différents instants, sont présentés figures 6.19 à 6.24. Le diaphragme est rompu à $t = 0.00$s (figure 6.19). Le champ aérodynamique, initialement cylindrique, est correctement décrit, et ce bien que des volumes cubiques soient utilisés. Comme le temps passe, un choc et une surface de contact apparaissent, se déplaçant vers le centre, alors qu'une onde de détente remonte dans la chambre haute pression (figure 6.20).
<figure><img src="image_99.png" /><figcaption>Figure 6.20: Interactions d'ondes 2D. Niveau de masse volumique à $t=0.20s$</figcaption></figure>

<figure><img src="image_100.png" /><figcaption>Figure 6.21: Interactions d'ondes 2D. Niveau de masse volumique à $t=0.30s$</figcaption></figure>
<figure><img src="image_101.png" /><figcaption>Figure 6.22: Interactions d'ondes 2D. Niveau de masse volumique à t=0.40s</figcaption></figure>

<figure><img src="image_102.png" /><figcaption>Figure 6.23: Interactions d'ondes 2D. Niveau de masse volumique à t=0.50s</figcaption></figure>
<figure><img src="image_103.png" /><figcaption>Figure 6.24: Interactions d'ondes 2D. Niveau de masse volumique à $t=0.70s$</figcaption></figure>

Les fronts d'onde sont bien circulaires, mais de petites vaguelettes indésirables voient
le jour à la suite des ondes. Ce défaut, exposé au paragraphe 3.4.3, est aussi visible sur
la solution obtenue en maillage non structuré avec la méthode de Roe (figure 3.18),
ou encore sur les résultats de Zhang et So [67]. D'après notre expérience, il est lié
à l'orientation relative d'une discontinuité mobile, par rapport aux frontières de la
cellule d'intégration qu'elle traverse. Comme l'explique Hirsch. [27, pp. 475-483],
la position du maillage influence la vitesse de déplacement des ondes. Dans le cas
particulier d'une onde circulaire discrétisée sur un maillage orthogonal, à rayon const-
tant et pour un angle variant entre 0 et 90 degrés, toutes les orientations possibles
existent. Il en résulte que la vitesse du front d'onde n'est pas uniforme, ce qui con-
duit à des mouvements angulaires de masse, qui engendrent les vaguelettes mises en
évidence. Ces irrégularités de vitesse sont alors d'autant plus marquées que le rayon
est faible, l'orientation de la discontinuité dans la cellule d'intégration changeant
rapidement d'une maille sur l'autre. Pour t autour de 0.30s, le choc se réfléchit sur
lui-même à l'origine, qui est un point singulier où les niveaux de pression et de masse
volumique sont très élevés (figure 6.21). Là, naît un nouveau choc se déplaçant vers
l'extérieur. Dans la mesure où le choc convergent n'est plus circulaire au voisinage de
l'origine, à cause du phénomène exposé précédemment, le choc réfléchi ne l'est pas
lui aussi. Si bien que le niveau de densité le long du front d'onde n'est pas uniforme
(figure 6.22). Comme le temps augmente, la discontinuité de contact initiale et le
choc réfléchi interagissent (figure 6.23), générant un choc réfléchi convergent et un
autre transmis divergent, ainsi qu'une nouvelle surface de contact. Pour le dernier
instant présenté figure 6.24, il est à noter que les vaguelettes entre le choc diver-
<figure><img src="image_104.png" /><figcaption>Figure 6.25: Interactions d'ondes 2D. Profils radiaux de densité à 0.10s d'intervalle</figcaption></figure>

gent et la discontinuité de contact sont moins marquées que celles observées juste après que le diaphragme ait été rompu, et que le niveau de densité le long du front d'onde, devenu sensiblement circulaire, tend à s'uniformiser comme l'onde s'éloigne de l'origine.

Les profils de masse volumique extraits le long de la diagonale du domaine d'étude, sont comparés figures 6.25 et 6.26, à ceux prédits par la méthode de Roe en maillage non structuré. Ils sont à mettre en parallèle avec ceux présentés dans le cas mono-dimensionnel (figures 6.15 à 6.18). Ils traduisent le même processus de propagations et d'interactions d'ondes, mais dans ce cas, à cause du confinement angulaire, les niveaux entre les ondes ne sont pas constants. Les discontinuités voient en effet leur amplitude augmenter comme elles s'approchent de l'origine, puis décroître comme elles s'en éloignent. Même dans cette configuration assez sévère, les deux approches numériques produisent des solutions quasi-similaires, tant au niveau de l'amplitude que de la vitesse des ondes. Comme rapporté précédemment, des vaguelettes indésirables apparaissent dans les deux résultats. Elles sont particulièrement visibles aux instants $t = 0.30s$ et $t = 0.40s$ entre la surface de contact et le pied de l'onde de détente (figure 6.25). Il est à noter que la méthode non structurée maintient une surpression à l'origine douteuse, après que le choc convergent s'y soit réfléchi. Ce défaut est lié à l'implémentation de la condition aux limites de glissement, qui est utilisée dans ce cas pour simuler les axes de symétrie. L'approche mise en oeuvre repose sur une formulation faible, où le flux de masse, et non pas la vitesse normale, est imposé nul à la frontière. Ce traitement se révèle incapable de prédire
<figure><img src="image_105.png" /><figcaption>Figure 6.26: Interactions d'ondes 2D. Profils radiaux de densité à 0.10s d'intervalle</figcaption></figure>

le niveau de masse volumique instantané à l'origine, à cause des évolutions rapides du champ aérodynamique qui s'y produisent. La solution obtenue présente ainsi un déphasage manifeste. Par contre, la méthode structurée tire profit de l'extension du domaine de calcul, qui permet une implémentation exacte et sans défaut de la condition de symétrie.

Pour valider ces résultats, la même comparaison est réalisée avec ceux publiés par Zhang et So [67] figures 6.27 et 6.28. Le comportement général prédit est le même. Cependant, notre choc est plus franc et les fronts de l'onde de détente sont moins arrondis. Ceci traduit une dissipation numérique moindre. Dans l'approche de Zhang et So, la stabilité du schéma est assurée par l'emploi aux itérations paires d'un schéma décentré amont alterné explicite, et aux itérations impaires du schéma de MacCormack. Une précision du second ordre en espace est ainsi obtenue par une technique classique du type prédicteur/correcteur. Aussi, des opérateurs différents du second ordre, agissant comme termes dissipatifs, sont utilisés aux instants pairs. Notre propre approche inclus elle aussi de tels effets. Mais, d'autre part, une expression du flux complet, et non pas du flux décentré, est utilisée dans le schéma de MacCormack, alors que cette méthode conduit généralement à un lissage des discontinuités. Cette seconde étape est sans doute responsable de l'étalement du choc dans les résultats de Zhang et So. Cependant, grâce à un traitement spécial de la surface de contact, tirant profit de la continuité de la vitesse et de la pression à travers cette onde, ils maintiennent parfaitement son caractère discontinu, alors que la nôtre est étalée sur 3 ou 4 mailles.
<figure><img src="image_106.png" /><figcaption>Figure 6.27: Interactions d'ondes 2D. Profils radiaux de densité à 0.10s d'intervalle</figcaption></figure>

<figure><img src="image_107.png" /><figcaption>Figure 6.28: Interactions d'ondes 2D. Profils radiaux de densité à 0.10s d'intervalle</figcaption></figure>
<figure><img src="image_108.png" /><figcaption>Figure 6.29: Tuyère à rafales. Géométrie</figcaption></figure>

Ces résultats confirment, et étendent au cas bi-dimensionnel, les conclusions portées précédemment. La méthode développée se révèle apte à gérer des discontinuités instationnaires, qui se propagent et interagissent éventuellement entre elles. Une grande précision est ainsi atteinte dans la prédiction de la vitesse et de l'amplitude de chaque onde, qui se comparent très favorablement aux solutions numériques utilisées comme référence.

## 6.5 Tuyère supersonique

La première macro-configuration présentée vise à simuler une tuyère expérimentale, dite à rafales, installée au Laboratoire de Mécanique des Fluides et d'Acoustique de l'Ecole Centrale de Lyon (figure 6.29). De part sa géométrie, l'écoulement est supposé bi-dimensionnel et stationnaire. Il est subsonique en entrée et supersonique en sortie, l'étude expérimentale se déroulant à tuyère bloquée durant une dizaine de secondes. Une description précise de la géométrie et des conditions aux limites se trouve en annexe E.4. Nous ferons l'hypothèse que le fluide simulé est parfait.

Le maillage utilisé (figure 6.29) met en oeuvre 83 × 27 noeuds dans le plan d'écoulement, et deux mailles dans la direction non significative perpendiculaire à ce dernier, où une condition de symétrie plane est appliquée. Une légère raison géométrique de 1.1 augmente la densité des noeuds près des parois supérieure et inférieure, pour permettre une meilleure prise en compte des effets de courbure. Les conditions génératrices en entrée sont fixées par les grandeurs d'arrêt de l'atmosphère dans laquelle aspire la tuyère, et la sortie étant supersonique, aucune condition physique n'y est imposée. La géométrie de calcul a été rallongée à l'amont par une portion rectiligne qui n'existe pas dans la réalité. Elle nous permet d'imposer un écoulement horizontal mono-dimensionnel. Cette modification située avant les premiers points de mesure ne semble pas affecter les résultats. Enfin, le schéma au troisième ordre en espace est utilisé dans les trois directions.

La figure 6.30 présente le niveau de nombre de Mach. Ce dernier varie de faiblement subsonique, $M = 0.17$, à supersonique, $M = 1.95$. La courbure des iso-lignes
<figure><img src="image_109.png" /><figcaption>Figure 6.30: Tuyère à rafales. Nombre de Mach</figcaption></figure>

<figure><img src="image_110.png" /><figcaption>Figure 6.31: Tuyère à rafales. Nombre de Mach sur la paroi inférieure</figcaption></figure>
<figure><img src="image_111.png" /><figcaption>Figure 6.32: Tuyère à rafales. Nombre de Mach sur la paroi inférieure</figcaption></figure>

est à noter, ainsi que leur changement d'orientation après le col, où l'accélération la plus importante se situe sur la paroi inférieure rectiligne. Des prises de pression statique, reliées à un multi-manomètre à mercure, sont disposées sur cette frontière. Elles nous permettent d'accéder au nombre de Mach expérimental, qui est comparée à la solution calculée figure 6.31. Jusqu'à la section $X = -0.2$, l'accord entre les deux courbes est très bon, l'incertitude expérimentale étant estimée à ±1mmHg. La légère différence visible à l'amont a été imputée aux modifications apportées à la géométrie en entrée. Par contre, dans le reste de la tuyère, le nombre de Mach calculé se stabilise autour de 1.9 alors que la valeur expérimentale ne dépasse pas 1.8. Dans l'expérience, un profil est inséré dans l'écoulement en milieu de veine à $X = 0.0$. Le choc se développant à partir de son bord d'attaque vient impacter sur la paroi à l'aval de cette section et est responsable des fluctuations de pression visibles sur les mesures en fin de tuyère. La simulation numérique ne le prenant pas en compte prédit quant à elle une évolution monotone. Nous noterons par ailleurs que la position de cet objet est judicieuse, puisqu'elle le place dans un écoulement quasiment uniforme (figures 6.30 et 6.31).

Sont présentées figure 6.32 les solutions calculées d'une part par une méthode analytique supposant un écoulement mono-dimensionnel fonction de la contraction de veine, d'autre part par la méthode de Roe en maillage non structuré. L'accord entre les deux solutions numériques est excellent; elles reproduisent en particulier de la même manière la courbure de la paroi inférieure pour $X \in [-0.5; -0.4]$, ainsi que l'évolution du nombre de Mach dans la seconde moitié de la tuyère. Dans cette zone est retrouvé le niveau analytique de 1.9. Aussi l'écart de 0.1 constaté sur la
<figure><img src="image_112.png" /><figcaption>Figure 6.33: Tuyère à rafales. Vitesse longitudinale</figcaption></figure>

répartition du nombre de Mach figure 6.31 existe-t-il entre le résultat expérimental et chacune des trois solutions numériques. Ce phénomène s'explique par la couche limite se développant sur les parois latérales qui n'est pas prise en compte lors des calculs en fluide parfait. Son épaisseur croît tout au long de la tuyère, en particulier sur la paroi supérieure dans le divergent. Une diminution de la section de veine en résulte, entraînant ainsi une réduction du nombre de Mach supersonique.

Ce comportement a longtemps été tenu pour responsable du léger déplacement du col sonique visible figure 6.32 autour de $X = -0.3$ entre les solutions analytique mono-dimensionnelle et expérimentale. En fait, il ressort de cette simulation numérique que ce phénomène est d'avantage dû à la nature bi-dimensionnelle de l'écoulement comme le montre la figure 6.30 et les figures 6.33 et 6.34, qui établissent que bien que la vitesse transversale soit faible, elle n'est pas particulièrement négligeable devant la composante longitudinale, et ce même au voisinage de la paroi inférieure. Cette incapacité de la méthode analytique à gérer les courbures de géométrie est aussi clairement visible figure 6.32 pour $X \in [-0.5; -0.4]$.

Figure 6.35 sont présentés les écarts relatifs entre pression totale locale calculée et pression totale imposée en entrée, et figure 6.36 ceux existants entre la température totale locale calculée et la température totale imposée en entrée. Ces valeurs sont nulles en théorie pour un écoulement adiabatique de fluide parfait semblable à celui simulé. Nous constatons quelques divergences par rapport à cette affirmation, de l'ordre de un pour cent pour la pression totale et de un pour mille pour la température totale. Ces défauts sont concentrés près des parois, comme ceci a été mis en évidence dans le cas de l'écoulement autour d'un cylindre. Il est difficile de dire si ces erreurs sont uniformément générées, ou si elles sont produites par exemple dans le convergent à l'entrée, puis convectées. Cependant, elles apparaissent reliées dans une certaine mesure aux irrégularités des frontières. Nous noterons en
<figure><img src="image_113.png" /><figcaption>Figure 6.34: Tuyère à rafales. Vitesse transversale</figcaption></figure>

<figure><img src="image_114.png" /><figcaption>Figure 6.35: Tuyère à rafales. Erreur relative sur la pression totale</figcaption></figure>
<figure><img src="image_115.png" /><figcaption>Figure 6.36: Tuyère à rafales. Erreur relative sur la température totale</figcaption></figure>

particulier une sur-estimation de la température totale juste après le col, couplée à une sous-estimation de la pression totale au même point. Or, si nous nous référons à la figure 6.32 et à la solution théorique, nous constatons dans cette zone une rupture de pente du profil de nombre de Mach, qui ne peut s'expliquer que par une discontinuité dans l'évolution de la section de la tuyère. La normale à la paroi supérieure reflète alors ce phénomène, et peut être la cause, via la condition de glissement, d'un comportement local peu physique.

En conclusion, cette configuration illustre la capacité de la méthode à simuler un écoulement présentant une forte évolution du nombre de Mach, ainsi que des effets bi-dimensionnels modérés. L'intégration du traitement des conditions aux limites avec le schéma mis en oeuvre dans le coeur du domaine est satisfaisante, même si de légères erreurs sur les grandeurs d'arrêt sont encore mal expliquées. Mais un aspect important de cette étude est l'utilisation qui en est faite pour analyser des résultats expérimentaux. Ces derniers étant confinés à la paroi inférieure rectiligne, ils ne réflètent pas toute la réalité de l'écoulement, et la compréhension de leur évolution gagne à être faite à la lumière de la solution numérique obtenue. Ainsi, nous pouvons expliquer le déplacement du col sonique, à cause d'effets bi-dimensionnels, valider la position retenue pour insérer un profil dans une région où l'écoulement est établi et faiblement variable, et une simulation en fluide visqueux, peut-être turbulent, permettrait d'estimer l'influence réelle des couches limites.

## 6.6 Tuyère transsonique

L'écoulement dans un canal convergent-divergent présente l'avantage d'être souvent bi-dimensionnel, susceptible de receler un choc si les conditions aux limites
<figure><img src="image_116.png" /><figcaption>Figure 6.37: Tuyère transsonique. Géométrie expérimentale</figcaption></figure>

adéquates sont choisies, et d'être raisonnablement proche de celui rencontré dans un canal inter-aubes de turbomachine, tout en étant plus simple. Pour toutes ces raisons, une tuyère symétrique transsonique, étudiée expérimentalement à l'Ecole Polytechnique Fédérale de Lausanne en Suisse, dans le département du Professeur Bölcs par P. Ott [47], a fait l'objet des simulations numériques présentées par la suite.

Ce canal comporte une contraction de veine dans une direction de l'ordre de 10% due à une bosse symétrique (figure 6.37). Sa forme est une adaptation de celle proposée par Adamson et Liou [1], et est décrite en détail par Ott, Bölcs et Franson [48]. Sa largeur est constante et égale à la moitié de sa hauteur en entrée. Le nombre de Mach amont est de l'ordre de 0.65, et les grandeurs d'arrêt en entrée et la pression statique en sortie sont telles qu'un choc se positionne dans le divergent. Il résulte du retour au régime subsonique de l'écoulement qui a été accéléré après le col jusqu'à des nombres de Mach proches de 1.4. La géométrie exacte ainsi que les valeurs des conditions génératrices sont regroupées annexe E.5.

Nous avons simuler pour cette configuration géométrique des écoulements stationnaires ou instationnaires, bi-dimensionnels en fluide parfait ou visqueux, en régime laminaire ou turbulent [24, 5]. Ces différents cas sont présentés par la suite.

### 6.6.1 Ecoulements bi-dimensionnels non visqueux

La première configuration étudiée vise à simuler l'écoulement bi-dimensionnel dans le plan de symétrie qui sépare les deux parois latérales, dans la partie supérieure du
<figure><img src="image_117.png" /><figcaption>Figure 6.38: Tuyère transsonique. Plan d'étude bi-dimensionnel</figcaption></figure>

<figure><img src="image_118.png" /><figcaption>Figure 6.39: Tuyère transsonique. Maillage bi-dimensionnel non visqueux</figcaption></figure>
Calcul fluide parfait. Schéma de Van Leer

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

Calcul turbulent. Schéma mixte $\epsilon' = 0.1$

<caption>Figure 6.40: Tuyère transsonique. Nombre de Mach</caption>

canal (figure 6.38). Comme lors de l'étude de la tuyère à rafales, une légère progres-
sion géométrique de raison 1.1 permet de resserrer les mailles près de la paroi solide
(figure 6.39). Cependant, la principale adaptation que subit le maillage est autour de
la position escomptée du choc, la taille des cellules variant dans un rapport de 1 à 10
entre l'amont et le voisinage de cette discontinuité. Il comporte 82 × 26 × 3 noeuds.
Le fluide simulé est supposé parfait et l'écoulement uniforme, mono-dimensionnel
en entrée. Nous utilisons la séparation de flux de Van Leer au troisième ordre. Une
condition de symétrie plane dans la direction non significative assure le caractère
bi-dimensionnel de l'écoulement, et le même traitement sur la frontière inférieure
simule l'autre moitié de la tuyère.

Dans la partie supérieure de la figure 6.40 est visible la répartition du nombre de Mach. Nous noterons le choc situé à la fin du divergent, et la déformation des iso-lignes qui s'accentue comme cette grandeur augmente. Ainsi, alors que le choc est quasiment vertical, le nombre de Mach sur l'axe dépasse légèrement 1.3, mais il atteint 1.4 à la paroi. Cette accélération est due à la courbure générée par la bosse.

Les profils de pression statique extraits le long de la tuyère sur l'axe et la paroi sont comparés figure 6.41 à la solution numérique obtenue par la méthode de Roe en maillage non structuré. Les échelles en bas et à gauche se rapportent aux profils sur l'axe, alors que celles en haut et à droite font référence aux données pariétales. Un très bon accord existe entre les deux approches, le choc étant capturé à la même position sur 2 ou 3 mailles. Cette dernière se compare très favorablement sur l'axe à celle prédite par une méthode analytique mono-dimensionnelle. Cependant, si les
<figure><img src="image_120.png" /><figcaption>Figure 6.41: Tuyère transsonique. Pression statique non visqueux</figcaption></figure>

résultats sur l'axe sont indifférenciables, un léger écart est visible à la paroi dans les zones de forte courbure.

Ce comportement se retrouve au niveau de la répartition de pression totale (figu-
re 6.42). Il apparaît comme un défaut de l'approche non structurée, qui ne conserve
pas cette grandeur à l'amont du choc. Nous avons pu le mettre en relation avec un
transfert de fluide à travers la paroi, la condition de glissement n'étant dans ce cas
pas parfaitement vérifiée. Le traitement de manière faible de cette condition aux
limites dans l'approche non structurée n'y est pas étranger. Ce problème excepté,
les deux méthodes prédisent le même saut de pression totale à travers le choc, ainsi
que des niveaux constants à l'amont et à l'aval de ce dernier. Il est à noter en par-
ticulier que le nombre de Mach variant dans la hauteur de la tuyère de 1.3 à 1.4
à cause de la courbure imposée par la paroi, l'intensité du choc doit augmenter, ce
qui est fidèlement reproduit par une perte de pression totale qui passe de 2.5% à
3.5%. Ce résultat est excellent sur l'axe où l'écoulement est mono-dimensionnel et
où la théorie prévoit dans ce cas qu'à un nombre de Mach amont de 1.32 est associé
un saut de pression totale de 2.4%. Par contre, un écart de 3.9% devrait exister à
la paroi à la vue du nombre de Mach de 1.39. La légère sous-estimation constatée
est imputable à la nature bi-dimensionnelle de l'écoulement dans cette région. La
présence de sur- et de sous-estimations de la perte de pression totale très localisées
près du choc, a été expliquée par ailleurs comme résultant d'un mauvais raccord du
schéma entre les régions supersonique amont et subsonique aval (paragraphe 3.4.3).
Les deux approches présentent ce défaut. Enfin, de légères oscillations sont visi-
bles sur les profils obtenus avec le schéma de Van Leer. Elles semblent reliées au
<figure><img src="image_121.png" /><figcaption>Figure 6.42: Tuyère transsonique. Pression totale non visqueux</figcaption></figure>

problème cité à l'instant. En effet, si la pression statique aval est maintenue constante, nous constatons une mise en résonance d'une oscillation de forte fréquence et de faible amplitude du choc, qui se couple avec des fluctuations de pression statique se propageant jusqu'à la sortie, puis se réfléchissant. Ce phénomène disparaît quasiment si une condition de non réflexion est appliquée à la frontière aval. Une application possible à ce comportement est la suivante. La position du choc changeant à l'intérieure des mailles dans la hauteur de la tuyère, les erreurs qui y sont associées évoluent aussi. Elles tendent à déséquilibrer la solution dans la direction verticale, qui se corrige tant bien que mal en générant une bouffée de pression statique. Cette dernière est convectée jusqu'à la sortie, où elle interagit avec le niveau imposé constant et se réfléchit, en changeant potentiellement d'amplitude. Lorsqu'elle atteint le choc, ce dernier "croit" que la pression statique de sortie a changé et se déplace en conséquence. Ainsi s'amorce la résonance. Ce phénomène n'apparaît pas en non structuré où une intégration temporelle implicite est utilisée. Il en résulte l'utilisation d'un CFL plus élevé, conduisant à un pas en temps plus important. Ainsi sont filtrés des ondes à haute fréquence, qui sont au contraire parfaitement décrites avec un petit pas en temps explicite. Une simple technique de lissage du résidu [41], filtrant les fluctuations rapides, contribue d'ailleurs à améliorer la convergence dans ce cas.

Ces résultats nous permettent d'affirmer que dans le cas d'un écoulement bidimensionnel en fluide parfait, stationnaire et transsonique, la méthode développée en maillage structuré avec la séparation de flux de Van Leer produit une solution semblable à l'approche de Roe en maillage non structuré. D'autre part, cette simu-
<figure><img src="image_122.png" /><figcaption>Figure 6.43: Tuyère transsonique. Maillage bi-dimensionnel visqueux</figcaption></figure>

lation montre que le choc est correctement décrit vis-à-vis de la physique comparé à la solution analytique mono-dimensionnelle tant au niveau de sa position sur l'axe que du saut de pression totale qu'il engendre. L'influence de la courbure de la paroi est aussi fidèlement reproduite.

## 6.6.2 Ecoulements bi-dimensionnels visqueux

L'étape suivante a consisté à simuler le même écoulement bi-dimensionnel, mais en fluide visqueux cette fois. Le maillage utilisé comporte le même nombre de noeuds, seule la raison géométrique change et passe à 1.3 afin d'augmenter la résolution de la couche limite se développant sur la paroi solide (figure 6.43). Les conditions en entrée sont conservées. En particulier, un profil de pression totale uniforme est imposé. Nous retranscrivons ainsi le fait qu'il n'y a pas de couche limite à l'amont de la géométrie simulée. La pression statique en sortie quant à elle a été diminuée de 3% pour assurer la même position au choc sur l'axe. Comme dans le cas précédent, la solution obtenue par la méthode de Roe en maillage non structuré, couplée à un modèle de turbulence $k-ε$ et à une loi de paroi, est utilisée comme base de comparaison numérique. Pour notre part, nous avons mis en oeuvre quatre schémas différents. Le premier décrit un écoulement supposé laminaire, aucun effet turbulent n'étant pris en compte. Cette configuration n'a rien de physique, mais elle offre une base intéressante pour l'explication des mécanismes dissipatifs présents. Les trois autres cas, turbulents quant à eux, utilisent le modèle de Michel (paragraphe 3.6.2). Deux d'entre eux mettent en oeuvre le schéma mixte, dans les trois directions, avec une valeur du coefficient de viscosité artificielle $ε^i$ de, respectivement, 0.1 et 0.5, et un nombre de Mach de raccord de $M_0^i = 0.6$ (paragraphe 3.4.5). Seul le schéma de Van Leer avec la correction de Hänel sur le flux d'énergie est utilisé dans le dernier cas.

La présence de la couche limite turbulente influence sensiblement l'écoulement sur la paroi supérieure, particulièrement à l'aval du choc, mais peu son coeur (figure 6.40, bas). Le nombre de Mach à l'amont de l'interaction choc/couche limite étant plus faible que dans le cas fluide parfait (figures 6.44 et 6.46), de l'ordre de 1.35, nous sommes dans la situation critique où un décollement de la couche limite et un choc en lambda peuvent apparaître. En fait, nous sommes en deçà de la limite,
<figure><img src="image_123.png" /><figcaption>Figure 6.44: Tuyère transsonique. Nombre de Mach non visqueux</figcaption></figure>

<figure><img src="image_124.png" /><figcaption>Figure 6.45: Tuyère transsonique. Nombre de Mach laminaire. Schéma de Van Leer</figcaption></figure>
<figure><img src="image_125.png" /><figcaption>Figure 6.46: Tuyère transsonique. Nombre de Mach turbulent. Schéma de Van Leer</figcaption></figure>

<figure><img src="image_126.png" /><figcaption>Figure 6.47: Tuyère transsonique. Nombre de Mach turbulent. Schéma mixte $\epsilon' = 0.5$</figcaption></figure>
<figure><img src="image_127.png" /><figcaption>Figure 6.48: Tuyère transsonique. Nombre de Mach turbulent. Schéma mixte $\epsilon^i = 0.1$</figcaption></figure>

<figure><img src="image_128.png" /><figcaption>Figure 6.49: Tuyère transsonique. Nombre de Mach turbulent. Schéma de Roe</figcaption></figure>
<figure><img src="image_129.png" /><figcaption>Figure 6.50: Tuyère transsonique. Pression totale laminaire. Schéma de Van Leer</figcaption></figure>

et seul un épaississement important de la zone à basse vitesse est visible à l'aval du choc (figure 6.40, bas).

Figures 6.45 à 6.49 sont présentées les iso-lignes de nombre de Mach. Le résultat laminaire (figure 6.45) est totalement différent des autres, et il ne sera pas discuté ici. Par contre, ceux en maillage structurés sont proches les uns des autres. Nous noterons cependant que, passant du schéma de Van Leer (figure 6.46), au schéma mixte à $\epsilon^i = 0.5$ (figure 6.47) puis au schéma mixte à $\epsilon^i = 0.1$ (figure 6.48), les lignes à l'amont du choc sont de plus en plus anguleuses près de la couche limite, et correspondent ainsi sur une hauteur de la tuyère plus importante au résultat en fluide parfait (figure 6.44). Ceci traduit un comportement de moins en moins diffusif. La solution non structurée (figure 6.49) est proche du schéma mixte à $\epsilon^i = 0.5$ à l'amont du choc, sauf pour l'iso-ligne 1.3, et assez différente à l'aval, où elle prédit une zone à basse vitesse plus importante. Ces divergences peuvent s'expliquer par l'utilisation des lois de paroi faite dans l'approche $k-\epsilon$. La répartition de la vitesse à l'amont du choc près de la courbure s'en trouve modifiée. D'autre part, il a été nécessaire pour cette méthode d'ajouter un niveau de turbulence de 4% dans l'écoulement libre afin d'assurer la stabilité du schéma. Couplée au modèle $k-\epsilon$, cette modification par rapport au cas structuré tend à générer une dissipation accrue à travers le choc.

L'examen du niveau de pression totale dans le cas laminaire (figure 6.50) par rapport au schéma mixte à $\epsilon^i = 0.1$ (figure 6.51, bas) révèle pourquoi l'introduction de la turbulence modifie grandement l'écoulement. En effet, le choc interagissant avec la couche limite tend à la faire décoller, générant ainsi un bulbe de recirculation. En l'absence d'une dissipation accrue par la viscosité turbulente, ce décollement s'étend et gagne en épaisseur. Il entraîne ainsi une contraction de veine apparente, et un changement de la courbure locale des lignes de courant en frontière de couche limite, provoquant un déplacement du choc et un changement de la nature de l'interaction choc/couche limite, qui s'approche alors d'une structure en lambda. Nous noterons que les pertes dans le convergent sont plus faibles dans le cas laminaire que dans celui turbulent, ce qui est logique puisque dans ce dernier, la viscosité turbulente s'ajoute à la contribution laminaire et accroît ainsi la dissipation visqueuse.

Figure 6.51 est comparé le saut de pression totale à travers le choc pour un fluide parfait (haut) et un écoulement turbulent (bas). Il est proche dans la première moitié de la hauteur de la tuyère à partir de l'axe, puis change radicalement de sens
## Calcul fluide parfait. Schéma de Van Leer

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

## Calcul turbulent. Schéma mixte $\epsilon^i = 0.1$

<figure><img src="image_131.png" /><figcaption>Figure 6.52: Tuyère transsonique. Pression totale turbulent. Schéma de Van Leer</figcaption></figure>

d'évolution. Dans le premier cas, il continue à augmenter, alors qu'il diminue puis augmente brutalement dans la couche limite dans le second cas. Ceci tend à établir qu'au voisinage de l'axe, les effets visqueux sont négligeables, ce qui est normal. Par contre, dans le troisième quart de la hauteur, avant les pertes importantes dues à la couche limite qui atteignent 40%, l'intensité du choc semble plus faible. Si nous nous référons aux iso-lignes de nombre de Mach (figure 6.48), nous remarquons que cette zone correspond sensiblement à la poche à Mach 0.8. A cause de l'interaction choc/couche limite, le nombre de Mach à l'amont de cette région est plus faible que dans le cas fluide parfait. Si bien que le saut de pression totale diminue.

Les pertes à l'aval du choc hors couche limite ont la même structure avec le

<figure><img src="image_132.png" /><figcaption>Figure 6.53: Tuyère transsonique. Pression totale turbulent. Schéma de Roe</figcaption></figure>
<figure><img src="image_133.png" /><figcaption>Figure 6.54: Tuyère transsonique. Pression statique</figcaption></figure>

schéma de Van Leer (figure 6.52) et avec le schéma mixte à $\epsilon^i = 0.1$ (figure 6.51, bas), la solution du schéma mixte à $\epsilon^i = 0.5$, non présentée ici, se situant entre les deux. La séparation de flux de Van Leer se révèle la plus dissipative des trois, puisque les pertes continues à augmenter même après le choc. Le résultat obtenu avec le modèle $k-\epsilon$ est assez différent (figure 6.53). Ce comportement est à mettre en relation avec les écarts constatés sur le nombre de Mach.

Le niveau de pression statique pour un écoulement laminaire (haut) et un écoulement turbulent (bas), est présenté figure 6.54. Nous noterons le déplacement du choc, ainsi que la structure de l'interaction en lambda dans le premier cas, clairement visible car occupant près de la moitié de la hauteur. Un fait intéressant est que jusqu'aux chocs, les solutions en pression sont quasiment identiques, ce qui implique une faible remontée d'informations le long de la couche limite, au demeurant très peu épaisse jusque là. Ce phénomène s'explique par le taux de turbulence près de la paroi qui reste très faible, pour ne devenir significatif qu'à l'aval du choc. La même comparaison est faite figure 6.55, mais cette fois avec le résultat en maillage non structuré. Si nous excluons le pied du choc, aucune différence sensible n'est notable.

Les profils de pression statique sur la paroi sont tracés figure 6.56. Jusqu'à $X=0.6$, toutes les solutions se confondent, s'il est fait abstraction du premier point. En effet, c'est là qu'est visible la seule différence sensible entre les résultats structurés et non structurés. Elle est due au traitement du coin supérieur gauche de la géométrie en lequel deux conditions physiques devraient être prises en compte, d'une part une condition de frontière libre, correspondant à l'entrée, d'autre part une condition d'adhérence, associée à la paroi solide. L'approche non structurée combine ces
Calcul turbulent. Schéma de Roe

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

Calcul turbulent. Schéma mixte $\epsilon^i = 0.1$

<caption>Figure 6.55: Tuyère transsonique. Pression statique</caption>

deux informations, et les impose de manière faible. Par contre, une seule condition physique, l'adhérence, est prise en compte dans le cas structuré. Le niveau d'énergie n'étant plus directement fixé, il s'en suit une légère chute de pression statique. Dans le cas laminaire, la structure en lambda est associée à l'évolution lente, induite par la courbure, visible entre $X = 0.8$ et $1.2$. Le choc non structuré est légèrement plus amont, et suivi d'une évolution moins rapide de pression statique, ce qui est compatible avec la zone à faible vitesse visible figure 6.49. Tout ceci se retrouve sur l'agrandissement présenté figure 6.57. Une information supplémentaire intéressante est que lorsque la viscosité numérique $\epsilon^i$ est diminuée, nous nous rapprochons de la solution non structurée: le choc se déclenche plus tôt et la pression évolue plus lentement, ce qui traduit un épaississement aval plus important. Cependant, les trois solutions structurées restent quand même très proches.

Figure 6.58 sont comparés les profils de pression statique, sur l'axe cette fois, résultant de la simulation en fluide parfait (figure 6.41), de celle en écoulement turbulent avec le schéma mixte à $\epsilon^i = 0.1$, et enfin des résultats expérimentaux publiés par Ott [47]. Un très bon accord existe jusqu'au col entre les trois solutions, dans cette partie de l'écoulement quasi mono-dimensionnelle où les effets visqueux sont insignifiants. Par contre, la raideur et l'amplitude du choc diffèrent sensiblement. Les écarts constatés en sortie sont directement reliés aux pertes, puisque rappelons que la pression statique aval a été modifiée dans la simulation numérique pour placer le choc à sa position expérimentale. Ainsi le niveau requis pour le calcul turbulent bi-dimensionnel est légèrement plus faible que celui imposé dans le cas non visqueux, afin de contrebalancer l'influence dissipative de la couche limite se développant sur
<figure><img src="image_135.png" /><figcaption>Figure 6.56: Tuyère transsonique. Pression statique turbulent en paroi</figcaption></figure>

<figure><img src="image_136.png" /><figcaption>Figure 6.57: Tuyère transsonique. Pression statique turbulent en paroi</figcaption></figure>
<figure><img src="image_137.png" /><figcaption>Figure 6.58: Tuyère transsonique. Pression statique turbulent sur l'axe</figcaption></figure>

les parois supérieure et inférieure. La différence beaucoup plus importante entre les résultats numériques et expérimentaux a la même origine, mais ce sont les frontières latérales qui sont impliquées cette fois. En effet, la configuration expérimentale est tri-dimensionnelle et présente pour particularité que les couches limites sur ces parois latérales sont très épaisses, occupant jusqu'à la moitié du canal [48]. Les simulations bi-dimensionnelles menées ici, ne prenant pas en compte les pertes qu'elles engendrent, requièrent une pression statique en sortie supérieure. D'autre part, les mesures présentées ont été réalisées à l'aide de prises de pression statique pariétale montées sur une des faces latérales au niveau de l'axe. Le champ aérodynamique a ainsi été mesuré à travers une couche limite épaisse, avec laquelle le choc interagit, tout comme il le fait avec celle beaucoup plus fine se trouvant sur la paroi bombée (figure 6.55). Il en résulte un étalement de la remontée en pression juste après le choc, qui est à mettre en parallèle avec celui visible figure 6.56. Ces couches limites, ignorées dans la simulation numérique, ont donc une influence non négligeable et rendent toute comparaison à l'expérience peu significative à l'aval du choc. Un calcul tri-dimensionnel, ou une installation expérimentale plus proche d'une configuration bi-dimensionnelle, seraient nécessaires pour ce faire.

Le profil de vitesse axiale au col (figure 6.59) montre comme attendu des évolutions turbulentes plus rapides comparées à la solution laminaire. D'autre part, le résultat non structuré est assez différent des autres, sans doute à cause de la loi de paroi qui donne ici une vitesse pariétale moitié de la vitesse débitante. Nous noterons encore qu'un profil de plus en plus "carré" est obtenu lorsque la viscosité artificielle est diminuée. En sortie, figure 6.60, le profil de vitesse laminaire présente encore une
<figure><img src="image_138.png" /><figcaption>Figure 6.59: Tuyère transsonique. Vitesse au col visqueux</figcaption></figure>

<figure><img src="image_139.png" /><figcaption>Figure 6.60: Tuyère transsonique. Vitesse en sortie visqueux</figcaption></figure>
poche négative, ce qui montre l'étendue de la zone de recirculation dans ce cas. Les quatre solutions turbulentes sont assez proches, avec deux extrêmes. D'une part le résultat avec le schéma de Van Leer présente une assez bonne récupération après l'interaction choc/couche limite, d'autre part le profil non structuré montre encore un déficit de vitesse pour $Y \in [0.35; 0.39]$. Les deux schémas mixtes produisent des solutions intermédiaires, d'autant plus proches de celle obtenue avec la séparation de flux de Van Leer que le coefficient de viscosité numérique $\epsilon^i$ est élevé.

En conclusion, nous noterons tout d'abord qu'un écoulement bi-dimensionnel transsonique, visqueux, turbulent, a été simulé avec succès. Les solutions en maillage structuré et non structuré sont jugées peu différentes, et les écarts mis en évidence sont sans doute liés aux modèles de turbulence, ce qui est malheureusement une conclusion classique de ce genre de comparaison. Le comportement de la méthode présentée dans cette étude est très satisfaisant, puisqu'elle ne génère aucune aberration et qu'elle reproduit qualitativement l'évolution physique. Seules des mesures expérimentales dans la configuration étudiée, à savoir un écoulement uniquement bi-dimensionnel sans couche limite amont, permettraient cependant une comparaison quantitative fiable.

Sur le plan numérique, contrairement à ce que l'étude du comportement de la séparation de flux de Van Leer à très faible nombre de Mach laissait présager (paragraphe 3.4.3), aucune surestimation catastrophique des pertes par effet visqueux n'a été constatée, dans ce cas particulier. Cependant, la comparaison au schéma mixte, qui montre dans une telle configuration son applicabilité, nous apporte l'enseignement suivant. Des trois solutions structurées turbulentes, c'est celle obtenue avec le schéma de Van Leer qui est la plus diffusive, ce qui explique une récupération plus rapide après l'interaction choc/couche limite. Ainsi, d'une extrapolation des résultats produits avec $\epsilon^i = 0.1$ et $0.5$, nous pouvons déduire que la séparation de flux de Van Leer peut être assimilée à un schéma mixte dont le coefficient de viscosité artificielle serait supérieur à 0.5. Par ailleurs, il serait instructif de connaître la valeur minimale de ce paramètre assurant la stabilité numérique, et de mesurer ainsi le plus grand bénéfice apporté par cette technique d'hybridation en fonction du nombre de Mach d'un schéma décentré amont et d'un schéma centré. Cette étude, grande consommatrice de puissance de calcul car la vitesse de convergence diminue avec la dissipation numérique, n'a pas encore été réalisée.

### 6.6.3 Ecoulement bi-dimensionnel instationnaire

Pour finir, une configuration instationnaire a été étudiée. Elle correspond à une extension de celle présentée au paragraphe 6.6.1, dans la mesure où le fluide simulé est supposé parfait et l'écoulement bi-dimensionnel. Cependant dans le cas présent, la pression statique de sortie n'est plus constante dans le temps, mais fluctue sinusoïdalement de manière forcée. La pression totale amont étant fixée, et les parois indéformables, il en résulte une oscillation axiale du choc. Notre but ici est d'étudier
<figure><img src="image_140.png" /><figcaption>Figure 6.61: Ecoulement instationnaire. Pression statique pour $P_8 = 0.6520$</figcaption></figure>

ce déplacement périodique.

L'écoulement étant supersonique après le col, il n'y a pas de transfert d'informations depuis la région du choc vers l'entrée. Aussi des perturbations à l'aval modifieront la position du choc, mais non la solution amont. C'est pourquoi le domaine d'étude a été adapté pour ne recouvrir que la fin du divergent, ainsi qu'une portion aval de la géométrie, de section constante, plus longue que dans le cas précédent (figure 6.61). Le champ aérodynamique dans le nouveau plan d'entrée est supersonique, et déduit des résultats obtenus précédemment (paragraphe 6.6.1). Le niveau de pression statique est quant à lui imposé, via le terme $B^*(t)$ (équation 4.70), dans le plan de sortie, qui a été déplacé vers l'aval pour se conformer à la disposition expérimentale [48]. Une description précise de cette configuration est faite en an-nexe E.5.3. Le maillage utilisé est visible figure 6.62. Les repères $P_1$ à $P_8$ y désignent la position des sondes expérimentales de pression pariétale. Pour une fréquence d'excitation de 80 Hz, le choc se déplace entre les points $P_3$ et $P_5$. L'amplitude de cette oscillation est portée figure 6.61. Dans cette région, des mailles plus petites ont été mises en oeuvre, et elles grossissent progressivement dans la portion rectiligne. Leur taille cependant est bien supérieure à celle utilisée dans la configuration statique (figure 6.39), mais elle permet ici une bonne capture du choc et de ses déplacements, tout en autorisant un pas en temps plus conséquent.

Comme précédemment, le niveau de pression statique a été adapté pour assurer
au choc la même position axiale que celle relevée expérimentalement. L'emplacement
<figure><img src="image_141.png" /><figcaption>Figure 6.62: Ecoulement instationnaire. Domaine d'étude</figcaption></figure>

<figure><img src="image_142.png" /><figcaption>Figure 6.63: Ecoulement instationnaire. Position instantanée du choc</figcaption></figure>
<figure><img src="image_143.png" /><figcaption>Figure 6.64: Ecoulement instationnaire. Pression statique instantanée de $P_1$ à $P_8$</figcaption></figure>

instantané du point où le nombre de Mach est 1.0, est tracé figure 6.63 sur une période. Y sont comparés les résultats expérimentaux, ceux prédits par la présente étude, référencés Proust, ainsi que ceux obtenus numériquement par Ott [47]. Un bon accord existe entre les trois approches, avec toute fois, de légères différences de phase et d'amplitude. Les premières sont simplement dues à l'imprécision de l'origine des temps utilisée par Ott, cette dernière semblant variée entre ses mesures expérimentales et numériques, ainsi qu'entre chaque configuration étudiée. Si bien que la synchronisation des trois solutions, grâce au signal de pression statique aval en $P_8$, faite "à la main", est d'une précision douteuse, et un léger glissement en temps des trois courbes en résulte. D'autre part, la position des sondes est portée à droite du graphe 6.63. En nous référant à la figure 6.62, nous constatons que trois mailles uniquement sont présentes entre les sondes $P_2$ et $P_3$, ce qui permet de relativiser l'écart constaté dans cette zone entre les deux solutions numériques. En effet, de part la formulation volumes finis, le choc est capturé à l'intérieur d'une maille, et la précision avec laquelle est évaluée sa position axiale est de l'ordre de la taille de cette cellule. La différence d'amplitude constatée pourrait donc être simplement due au fait que les mailles utilisées par Ott sont trois à quatre fois plus longues que celles présentement mises en oeuvre.

Figure 6.64 sont tracés les niveaux de pression statique calculés en chaque sonde durant une période par la présente méthode. Comme établi par le graphe 6.63, l'écoulement en $P_1$ reste supersonique et le niveau de pression statique y est constant. Les sondes $P_3$ à $P_5$ sont alternativement à l'amont du choc, dans une zone supersonique à pression minimale constante, et à l'aval du choc, dans une zone sub-
<figure><img src="image_144.png" /><figcaption>Figure 6.65: Ecoulement instationnaire. Pression statique instantanée en $P_3$</figcaption></figure>

sonique à pression maximale. Une légère augmentation de pression est constatée en $P_2$, par rapport à un niveau uniforme par ailleurs. Comme le révèle la figure 6.63, cette sonde est très proche de la position extrême amont de la ligne à Mach 1.0, qui ne l'atteint pas cependant. Or, le choc est décrit numériquement comme une augmentation brutale de pression, non pas ponctuelle, mais étalée sur 2-3 mailles, au milieu desquelles est située cette ligne à Mach 1.0. La sonde $P_2$ se trouve manifestement dans cette région à fort gradient lorsque le choc atteint sa position extrême amont. Les courbes $P_3, P_4$ et $P_5$ présentent une certaine "symétrie", qui laisse envisagée une évolution du champ aérodynamique elle aussi "symétrique" lorsque le choc se déplace de l'aval vers l'amont, puis de l'amont vers l'aval. Cette opinion est confortée par la nature manifestement sinusoïdale qu'affiche l'emplacement du choc figure 6.63. Il n'existe donc pas de différence très sensible en une position fixe de la géométrie lors du passage du choc, qu'il remonte ou descende dans la tuyère. Enfin, l'écoulement en $P_6, P_7$ et $P_8$ reste subsonique. La pression statique y évolue sinusoïdalement, avec un déphasage entre les trois courbes qui augmente comme la distance séparant les sondes croît. Il correspond au temps nécessaire à la propagation de l'excitation.

Les résultats obtenus expérimentalement et numériquement par Ott en $P_3$ et $P_5$ sont comparés à la solution présente figures 6.65 et 6.66. Les deux solutions numériques sont considérées comme proches. La différence de phase a déjà été expliquée par ailleurs, et les changements de pression moins rapides prédits par Ott sont liés à la taille accrue de ses mailles. En effet, bien que dans les deux approches, le choc puisse être décrit sur le même nombre de cellules, le fait que les mailles
<figure><img src="image_145.png" /><figcaption>Figure 6.66: Ecoulement instationnaire. Pression statique instantanée en $P_s$</figcaption></figure>

mises en oeuvre dans la présente étude soient plus petites, confine l'évolution de la pression statique sur une plus courte distance, et conduit ainsi a une évolution plus rapide en temps, pour la même vitesse de déplacement du choc. Par contre, un écart sensible existe entre les solutions numériques et expérimentales. Il est normal que les niveau à l'aval du choc different, puisque rappelons que la pression statique en sortie a été modifiée pour assurer la même position axial au choc, et la même amplitude à son oscillation (figure 6.63). D'autre part, l'étalement important en temps des résulats expérimentaux est à mettre en relation avec les remarques faites au paragraphe 6.6.2 (figure 6.58). L'interaction du choc et de la couche limite épaisse sur les parois latérales, à travers laquelle sont effectuées les mesures, provoque une évolution de pression statique moins rapide d'une part, et plus étendue d'autre part. L'approche du choc est anticipée, et l'influence de son passsage persiste plus longtemps dans la couche limite.

Nous établissons avec ces résultats la capacité de la méthode développée à simuler le déplacement d'une discontinuité sous l'influence d'une contrainte extérieure imposée, dans ce cas le niveau de pression statique en sortie. La comparaison qualitative aux solutions expérimentales et numériques utilisées comme référence, montre la qualité des prédictions obtenues. Une investigation quantitative nécessiterait cependant un dispositif expérimental s'affranchissant des couches limites latérales perturbatrices, ou la prise en compte dans la simulation numérique de ces dernières. Un calcul tri-dimensionnel instationnaire et turbulent devrait alors être mené...
<figure><img src="image_146.png" /><figcaption>Figure 6.67: Grille d'aubes supersonique. Géométrie</figcaption></figure>

## 6.7 Grille d'aubes supersonique

Cette configuration stationnaire est bâtie sur un profil réel de pales du distributeur du deuxième étage de la turbine de la turbopompe à hydrogène du moteur VULCAIN développé par la SEP, à partir duquel est construite une grille d'aubes bi-dimensionnelle. Figure 6.67 est dessiné le canal inter-aubes que nous simulons. La déflexion de l'écoulement y atteint 110 degrés. Deux portions rectilignes, dont la longueur correspond à une demi-corde, ont été ajoutées à l'amont et à l'aval du canal proprement dit. Une condition de périodicité, non-angulaire, est appliquée sur les faces horizontales de ces deux domaines. Nous reconstituons ainsi un empilement infini d'aubes dans la direction verticale du dessin, pendant qu'une condition de symétrie plane permet de modéliser une envergure infinie perpendiculairement au plan d'écoulement. Bien qu'aucune donnée expérimentale sur cette géométrie particulière ne soit disponible, cette configuration permet de valider les potentialités de la méthode sur une situation réelle en turbomachine.

Le maillage utilisé (figure 6.67) comporte 133 x 42 noeuds, et deux mailles en envergure. Le bord d'attaque et le bord de fuite étant ronds, ces régions ont fait l'objet d'une attention toute particulière. Ainsi, un grand nombre de mailles ont été positionnées autour du bord d'attaque afin de décrire correctement la courbure, comme l'atteste la figure 6.68. Par contre, nous nous sommes contentés de placer un noeud unique au culot du bord de fuite (figure 6.69). L'écoulement uniforme en entrée est calé approximativement à 50 degrés, pour un nombre de Mach de l'ordre de 0.3. L'écoulement en sortie est supersonique, de l'ordre de Mach 1.5. Cependant, le nombre de Mach axial en sortie étant subsonique, de l'ordre de 0.6, un niveau de
<figure><img src="image_147.png" /><figcaption>Figure 6.68: Grille d'aubes supersonique. Bord d'attaque</figcaption></figure>

<figure><img src="image_148.png" /><figcaption>Figure 6.69: Grille d'aubes supersonique. Bord de fuite</figcaption></figure>
<figure><img src="image_149.png" /><figcaption>Figure 6.70: Grille d'aubes supersonique. Nombre de Mach</figcaption></figure>

pression statique doit être imposée sur cette frontière. Un premier calcul a été réalisé avec une telle condition aux limites, jusqu'à ce qu'une décroissance de trois ordres de grandeur du résidu soit atteinte. Puis, nous avons remplacé cette condition de frontière libre par une condition de non réflexion et mené la simulation à convergence. Ces deux étapes permettent d'établir le régime souhaité, tout en évitant la réflexion du sillage et des chocs de bord de fuite sur la frontière avale. Le fluide simulé est supposé parfait.

Figure 6.70 est représenté le nombre de Mach dans le canal. Nous noterons sa très faible valeur au bord d'attaque où se situerait le point d'arrêt en fluide visqueux, l'accélération subie par le fluide sur l'extrados du fait de la courbure, ainsi que l'éventail de trois ondes partant du bord de fuite. S'y trouvent deux chocs obliques de chaque côté du sillage géométrique du profil. Ces deux ondes de pression sont dues au changement brutal de direction que subit le fluide au culot, couplé à la nature supersonique de l'écoulement en sortie de canal. Si une condition de glissement est imposée en ce point, nous constatons une chute dramatique de la pression conduisant à une divergence du calcul. Afin de minimiser ce problème classique en fluide non visqueux de raccord de l'extrados et de l'intrados [27, pp. 395-401], nous nous sommes contentés de moyenner le champ aérodynamique au culot. La solution
<figure><img src="image_150.png" /><figcaption>Figure 6.71: Grille d'aubes supersonique. Pression statique</figcaption></figure>

obtenue est localement non physique, mais globalement satisfaisante.

Figure 6.71 est visible le niveau de pression statique, qui présente la même évolution que le nombre de Mach. En particulier, le choix d'une palette de couleurs adaptée à l'étude de la zone aval fait clairement apparaître les ondes naissant au bord de fuite (figure 6.72). Figure 6.73 a été reconstituée la périodicité sur quelques canaux. Il est à noter tout d'abord que le raccord entre le domaine de calcul et les autres canaux fantômes non simulés est indécelable, montrant ainsi la grande qualité de la condition de périodicité implémentée. Les deux chocs de bord de fuite sont parfaitement résolus, celui naissant sur l'intrados venant se réfléchir proprement sur l'extrados, alors que le deuxième de plus faible amplitude est rectiligne, faiblement incliné par rapport à l'axe. Grâce à la condition de non réflexion, ces ondes ne se réfléchissent pas en sortie, bien que le nombre de Mach perpendiculairement à cette frontière soit subsonique. Seule une légère distorsion du champ aérodynamique est visible sur la dernière rangée de mailles, liée vraisemblablement à l'extrapolation linéaire dans les plans fictifs (paragraphe 4.2.2).

Les niveaux d'erreurs sur la pression totale (figure 6.74) et la température totale (figure 6.75) montrent encore des défauts importants au voisinage des frontières de glissement, sur l'intrados et l'extrados (paragraphe 4.3.3). Au bord d'attaque en particulier, des sur- et des sous-estimations se concentrent. Une question non résolue
<figure><img src="image_151.png" /><figcaption>Figure 6.72: Grille d'aubes supersonique. Pression statique</figcaption></figure>
<figure><img src="image_152.png" /><figcaption>Figure 6.73: Grille d'aubes supersonique. Pression statique</figcaption></figure>
<figure><img src="image_153.png" /><figcaption>Figure 6.74: Grille d'aubes supersonique. Pression totale</figcaption></figure>
<figure><img src="image_154.png" /><figcaption>Figure 6.75: Grille d'aubes supersonique. Température totale</figcaption></figure>
<figure><img src="image_155.png" /><figcaption>Figure 6.76: Grille d'aubes supersonique. Géométrie grossière</figcaption></figure>

est de savoir si ces erreurs sont convectées le long des parois depuis cette région, ou si elles sont produites localement. Le niveau cependant est inférieur à un pour cent dans la majorité du canal. Nous noterons aussi les spots traversant du bord de fuite à l'extrados, qui correspondent à un des deux chocs faibles. Ils pourraient être dus au mauvais raccord du schéma à travers un choc (paragraphe 3.4.3). Enfin, le traitement succinct du culot n'est sans doute pas étranger aux erreurs présentent en sortie.

Nous avons simulé la même configuration avec un maillage comportant deux fois moins de mailles dans la direction aube à aube, soit $133 \times 21 \times 3$ noeuds (figure 6.76).

Si la répartition du nombre de Mach (figure 6.77) est globalement proche de celle obtenue sur la géométrie finie, en particulier en sortie, des différences sensibles voient le jour au voisinage du bord d'attaque, remarquables sur l'extrados. Elles sont dues à un changement du calage de l'écoulement dans cette région. Des mailles plus grosses, décrivant moins bien la courbure du bord d'attaque, étant utilisées, le fluide se répartit différemment entre les deux faces du profil, et conduit à ces écarts de comportement.

Ils sont clairement visibles sur le niveau de pression totale (figure 6.78), qui est sous-estimé maintenant à l'extrados alors qu'il est sur-estimé avec le maillage fin (figure 6.74). La région où une erreur notable existe reste cependant confinée aux parois. Elle apparaît plus épaisse dans ce cas, mais ceci est exclusivement dû à la taille accrue des mailles. L'erreur générée au bord d'attaque se propage donc le long du profil, et elle n'est pas produite localement, puisque rien n'a été changé au niveau de l'extrados, si ce n'est la densité des noeuds, et que dans une configuration
<figure><img src="image_156.png" /><figcaption>Figure 6.77: Grille d'aubes supersonique. Nombre de Mach en géométrie grossière</figcaption></figure>
<figure><img src="image_157.png" /><figcaption>Figure 6.78: Grille d'aubes supersonique. Pression totale en géométrie grossière</figcaption></figure>
il y a sur-estimation et dans l'autre sous-estimation. Cette conclusion est en accord avec celle proposée par Dadone [15], qui a établi l'influence du décentrement et de l'extrapolation sur les noeuds supplémentaires lors du traitement d'une frontière de glissement. Il a montré en particulier que plus cette dernière opération est peu précise, plus la taille des mailles joue un rôle important. En d'autres mots, une extrapolation de bonne qualité sera peu sensible au maillage, alors q'une médiocre le sera excessivement. Nous retrouvons à travers cette configuration ce comportement, puisque l'extrapolation linéaire utilisée, qui est sans doute la plus simple de toutes, conduit à des solutions, au niveau des pertes tout au moins, qui sont fortement tri-butaires de la taille des cellules. D'autre part, se superposent à ce phénomène les erreurs spécifiques à la séparation de flux de Van Leer [26] lorsque le nombre de Mach est faible, ce qui est précisément le cas au bord d'attaque. Ce défaut, combiné à celui introduit par l'extrapolation, contribue à un calcul erroné du champ aérodynamique local au début du profil, puis cette erreur est convectée le long des parois.

Les résultats satisfaisant obtenus en particulier sur la géométrie finement maillée tendent à montrer l'aptitude de la méthode développée à traiter des configurations de turbomachine, statiques il est vrai dans ce cas. Le comportement de l'écoulement est qualitativement correct et cohérent sur le plan physique. Seule une mauvaise répartition des grandeurs d'arrêt est à déplorer. Mais il apparaît qu'elle résulte d'avantage d'un phénomène d'erreurs locales qui sont convectées, que d'une production répartie le long des frontières. Une analyse plus soignée du bord d'attaque, et un traitement moins grossier du bord de fuite, contribueraient sans doute à améliorer grandement la solution.

## 6.8 Coude subsonique

Cette dernière macro-configuration met en oeuvre un coude tri-dimensionnel, semblable à ceux présents entre deux roues d'une machine radiale. Il a été dessiné au moyen de DegrN, un outils développé au Laboratoire de Mécanique des Fluides et d'Acoustique de l'Ecole Centrale de Lyon par Duchemin et al. [17]: Sa ligne moyenne est circulaire plane et présente une courbure de 90 degrés angulaires, entre les plans d'entrée et de sortie, alors que dans le même temps, sa section, carrée à l'entrée, se déforme et subit elle aussi une rotation de 90 degrés autour de cet arc de cercle (figure 6.79). L'écoulement est stationnaire, tri-dimensionnel, et très faiblement compressible, le nombre de Mach variant entre 0.07 et 0.27. Le fluide simulé est supposé parfait.

Le maillage utilisé (figure 6.79) met en oeuvre 31 × 15 × 15 nocuds, régulièrement répartis. Une condition physique forte de glissement est imposée sur les quatre faces latérales et un traitement explicite a posteriori des arêtes en assure "l'étanchéité". La formulation faible n'a pu être utilisée dans ce cas... faute de temps! Un gradient
<figure><img src="image_158.png" /><figcaption>Figure 6.79: Coude subsonique. Géométrie</figcaption></figure>
uniforme de pression totale, perpendiculaire au plan de la ligne moyenne, est im-
posé en entrée (figure 6.81), ainsi qu'un niveau uniforme de température totale. La
direction de l'écoulement y est supposée perpendiculaire au premier plan de calcul.
Par ailleurs, la sortie du canal étant subsonique, la répartition de pression statique y
est fixée (figure 6.82). La solution, purement incompressible, prédite par DegrN est
utilisée comme référence. Cet outil met en oeuvre un schéma au premier ordre en
espace le long de la ligne moyenne, et un développement de Taylor d'ordre 3 dans les
sections perpendiculaires à cette dernière. C'est pourquoi, dans ce cas, l'approche
 faisant l'objet de cette étude, référencée comme Proust par la suite, utilise le schéma
au premier ordre dans la direction principale, et celui au troisième ordre dans les
directions transversales.

Figure 6.80 est comparée la trace des vecteurs vitesses, dans les sections 1,11,21 et 28, prédits par les deux approches. La technique de visualisation est du type "brin de laine", et le symbole associé à chaque segment repère son origine, et non pas son extrémité. L'écoulement étant perpendiculaire au plan d'entrée, la projection de la vitesse y est nulle (figure 6.80, haut). Puis, l'influence combinée du gradient de pression totale et de la courbure, génère un tourbillon, apparenté au classique tourbillon de passage d'un canal inter-aubes, dont l'intensité va croissante (figure 6.80, de haut en bas). Les deux méthodes numériques prédissent qualitativement le même phénomène, et des niveaux de vitesse semblables. Cependant, le centre tourbillonnaire, en particulier dans le dernier plan, n'est pas situé au même endroit. Cette différence a été imputée aux traitements des conditions aux limites de glissement. En effet, de par son implémentation, DegrN ne satisfait pas la condition de tangence de la vitesse aux parois de la géométrie, mais minimise une fonctionnelle, bâtie sur le produit scalaire de la vitesse et de la normale, le long d'une courbe fermée autour de la ligne moyenne, incluse dans la section réelle. Il en résulte une "transpiration" du fluide au travers des frontières solides. Par contre, la condition de glissement est parfaitement respectée par Proust, et elle entraîne un déplacement du tourbillon, qui est repoussé vers le centre du canal par des parois totalement étanches cette fois.

Les écarts relatifs de pression totale, par rapport à la valeur sur la ligne moyenne dans le plan d'entrée, sont tracés figure 6.81. Cette grandeur étant imposée en entrée, ils y sont rigoureusement identiques. L'évolution prédite est qualitativement la même, les particules fluides à haute énergie restant proches de la paroi supérieure, alors que celles à basse énergie migrent vers le coin inférieur gauche. D'autre part, un bon accord quantitatif existe au voisinage de la ligne moyenne. Il se dégrade progressivement comme nous nous en éloignons, ainsi que près de la sortie. La forme en "S" des iso-lignes calculées par DegrN dans les sections intermédiaires, fournit une explication partielle à ces différences. En effet, elle est caractéristique d'une fonction polynômiale impaire de degrés 3, semblable à celle utilisée dans les directions transversales par DegrN... Si ce polynôme offre une précision satisfaisante autour de
<figure><img src="image_159.png" /><figcaption>Figure 6.80: Coude subsonique. Vitesses</figcaption></figure>
<figure><img src="image_160.png" /></figure>

<table><caption>Figure 6.81: Coude subsonique. Pression totale relative</caption><thead><tr><th>DegrN</th><th>Proust</th></tr></thead><tbody><tr><td>Minimum:</td><td>-0.015</td></tr><tr><td>Maximum:</td><td>0.025</td></tr><tr><td colspan="2">11 iso-lignes</td></tr></tbody></table>
l'origine, il présente par contre un comportement asymptotique peu réaliste, qui déforme sensiblement les résultats loin de la ligne moyenne, et conduit en particulier à une augmentation non physique de la pression totale dans le plan 28, où cette grandeur atteint 0.036 pour un niveau maximum en entrée de 0.022. D'autre part, le schéma au premier ordre utilisé par Proust est dissipatif, et responsable d'une diminution de la pression totale, en particulier dans le coin inférieur gauche, qui passe de -0.0125 à -0.014.

Les évolutions de la pression statique, présentées figure 6.82, appellent les mêmes remarques. Un accord certain est constaté au centre du canal, et il se dégrade comme nous nous éloignons de la ligne moyenne. Il est à noter cependant que, malgré un nombre de Mach très faible, les iso-lignes calculées par la présente méthode sont régulières, et n'ont pas ce comportement chahuté souvent rencontré lors du calcul d'un écoulement quasi incompressible par une méthode compressible. Il est classiquement dû à la faible variation relative du niveau de pression statique, qui est difficile à reproduire numériquement. Tel n'est pas le cas ici.

En conclusion, cette configuration illustre la capacité de la méthode développée à simuler un écoulement interne fortement tri-dimensionnel, à bas nombre de Mach. A la lumière des spécificités de chaque approche, il apparaît que le champ aérodynamique calculé se compare très favorablement à celui prédit par DegrN, utilisé comme référence. Un bon accord quantitatif a été constaté au voisinage de la ligne moyenne du canal, zone dans laquelle la précision des deux méthodes est proche. Ceci établit l'usage qui peut être fait de Proust en tant qu'outil d'analyse, combiné à des moyens de dessin.

## 6.9 Conclusion

La méthode numérique développée a ainsi été utilisée pour simuler et analyser plusieurs configurations d'écoulements.

Il ressort d'une première série de résultats en fluide parfait, qu'elle est apte à gérer des discontinuités instationnaires qui se propagent et interagissent éventuellement entre elles. Une grande précision est atteinte dans la prédiction de la vitesse et de l'amplitude de chaque onde, comparées aux solutions mono-dimensionnelles analytiques et numériques utilisées comme référence. Cette propriété est maintenue lorsque l'écoulement devient multi-dimensionnel, bien que dans ce cas, une sensibilité de la méthode à l'orientation relative des discontinuités et des mailles ait été mise en évidence.

L'étude d'écoulements confinés entre deux ou plusieurs parois solides, a révélé quant à elle la grande qualité du traitement des conditions aux limites mis en oeuvre. Ainsi, les frontières libres d'entrée/sortie sont parfaitement intégrées avec le reste du domaine de calcul, sans anomalie au voisinage de ces régions. D'autre part, les parois solides, glissantes ou adhérentes, sont correctement prises en compte. Un
<figure><img src="image_161.png" /></figure>

<table><caption>Figure 6.82: Coude subsonique. Pression statique</caption><thead><tr><th>DegrN</th><th>Proust</th></tr></thead><tbody><tr><td>Minimum:</td><td>23.65</td></tr><tr><td>Maximum:</td><td>24.33</td></tr><tr><td colspan="2">11 iso-lignes</td></tr></tbody></table>
choc stationnaire, bi-dimensionnel, a été simulé avec succès pour un fluide parfait ou visqueux. Il est résolu sur 2 ou 3 mailles uniquement. Enfin, des calculs turbulents multi-dimensionnels ont produit des résultats présentant un accord satisfaisant avec la solution numérique utilisée comme référence. Ils ont en particulier établi le bon comportement du schéma mixte, qui hybride en fonction du nombre de Mach un schéma centré et la séparation de flux de Van Leer, ce qui permet de remédier très facilement au mauvais comportement de cette dernière dans les couches limites. Par ailleurs, la simulation des oscillations forcées d'un choc bi-dimensionnel, en réponse à une fluctuation sinusoïdale de pression statique à l'aval, a été réalisée et les résultats obtenus, dans ce cas en fluide parfait, révèlent un comportement en temps et en espace très satisfaisant.

Enfin, deux configurations proches de la réalité en turbomachine ont montré la capacité de la méthode développée à traiter des géométries complexes. L'étude d'une grille d'aubes supersonique a confirmée la qualité de la gestion des frontières libres, puisqu'un domaine minimal de deux fois la longueur de la corde a pu être utilisé, ce qui est très peu. Et ce malgré la présence de deux chocs en bord de fuite, dont l'extension et la réflexion de l'un d'eux sur l'extrados ont par ailleurs été correctement simulées. Cette configuration a mis en évidence la nécessité d'un traitement approprié des noeuds supplémentaires, pour lesquels l'extrapolation linéaire mise en oeuvre est sans doute un peu trop frustre. Pour finir, le calcul de l'écoulement faiblement compressible dans un canal tri-dimensionnel a établi le bon comportement de la méthode à bas nombre de Mach.
# Chapitre 7

# Conclusions et perspectives

La compréhension des phénomènes physiques régissant un écoulement transsonique instationnaire, voir turbulent, ne peut être faite uniquement sur la base de résultats expérimentaux. C'est pourquoi l'objet principal de ce travail a été l'étude de schémas numériques à haute précision appliqués à une telle configuration.

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 traités 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é.

Le code de calcul ainsi construit a été utilisé pour simuler et analyser quelques configurations d'écoulements caractéristiques. Une description de grande qualité de discontinuités instationnaires et de leurs interactions a été obtenue. D'autre part, il a permis la simulation d'é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.

Cependant, un certain nombre de limitations ont été mises en évidence. D'une part, le traitement des flux convectifs s'est révélé sensible à l'orientation relative des discontinuités et des cellules d'intégration, ainsi que sujet à un comportement peu
physique au raccord de régions supersonique et subsonique. D'autre part, la qualité de la gestion des parois glissantes ou adhérentes apparaît comme fortement tributaire des nocuds supplémentaires introduits, et de la façon dont le champ aérodynamique y est déterminé.

A la vue de ces résultats, nous pouvons qualifier la méthode d'opérationnelle pour simuler des écoulements en fluide parfait ou visqueux, multi-dimensionnels, potentiellement instationnaires, et ce du régime faiblement compressible à superso-nique.

Les conditions de fiabilité numérique ayant été remplies, nous pouvons main-tenant envisager une utilisation de plus en plus appliquée à la physique du code développé. Ainsi, dans le prolongement de ce travail, la simulation du mouvement oscillant d'un choc dans une tuyère tri-dimensionnelle est en préparation en colla-boration avec le Professeur Bölcs et le Docteur Ott de l'EPFL. D'autre part, une première étude des phénomènes de couplage fluide/structure, en particulier du flotte-ment, est envisagée avec le Docteur Dawes de l'Université de Cambridge. Rappelons ici que notre but à long terme est une meilleure compréhension des mécanismes de transfert instationnaire de l'information au sein du milieu fluide, et de l'influence des effets non linéaires introduits par un choc par exemple.

Il faut cependant garder à l'esprit qu'une activité numérique devra être main-tenue. Elle pourra concerner entre autres des aspects géométriques, tels que le mouvement du maillage pour prendre en compte la rotation ou l'adaptation dy-namique, ainsi que l'extension à des domaines multiples, et des aspects temporels, avec l'accroissement de la précision en temps par utilisation d'un schéma de Runge-Kutta par exemple, ce dernier couplé à des techniques de lissage de résidu per-mettant par ailleurs une accélération de la convergence vers l'état stationnaire. Un travail complémentaire devrait être d'autre part accompli sur les flux convectifs, pour intégrer des limiteurs adaptés et réduire la dépendance géométrique, ainsi que sur les conditions aux limites, où un traitement plus élaboré des noeuds supplémentaires est à étudier, et où les potentialités de l'approche faible restent encore à explorer.
# Annexe A

## Moyenne des équations de Navier-Stokes

### A.1 Introduction

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, est coûteuse en temps de calcul et en taille mémoire, compte tenu de l'importance des maillages à mettre en oeuvre. C'est pourquoi une décomposition du champ aérothermodynamique en une composante moyenne et une composante fluctuante, est généralement utilisée. 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.

### A.2 Equations instantanées

Rappelons que les équations à moyenner s'écrivent:

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

$$ \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 (A.2) $$

$$ \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 (A.3) $$

avec:

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

$$ \vec{\tau}_l = \lambda \vec{\nabla} \cdot \vec{V} \vec{I} + \mu (\vec{\nabla} \vec{V} + \vec{V} \vec{V}^t) \qquad (A.5) $$
## A.3 Moyennes employées

Dans une approche très générale, la valeur moyenne $\bar{A}$ d'une variable A, est définie à partir de plusieurs réalisations successives de la même expérience. P étant la densité de probabilité de A, nous avons alors:

$$ \bar{A}(\vec{x}, t) = \int_{-\infty}^{+\infty} A(\vec{x}, t, n)P(\vec{x}, t, n)dn \quad (A.6) $$

et, $A'(\vec{x}, t)$ étant la partie fluctuante de A,

$$ A(\vec{x}, t) = \bar{A}(\vec{x}, t) + A'(\vec{x}, t) \quad (A.7) $$

Sur le plan numérique, mais aussi expérimental, cette moyenne statistique est difficile d'emploi. C'est pourquoi une moyenne dans le temps lui est souvent substituée. Trois échelles temporelles sont alors à considérer lors de l'étude d'un écoulement instationnaire. Tout d'abord, le comportement aléatoire fluctuant de la turbulence introduit un temps turbulent caractéristique $\tau'$. D'autre part, soit $\tau$ le temps caractéristique du phénomène instationnaire. Enfin, soit $\bar{\tau}$ la période durant laquelle est moyennée le champs aérothermodynamique. Dans la présente étude, nous supposerons que l'instationnarité de l'écoulement résulte de causes macroscopiques, souvent externes, qui sont sans lien direct avec la turbulence du milieu simulé. C'est le cas par exemple lors de la vibration d'une aube, ou au cours du défilement d'un sillage. Aussi, $\tau$ est supposé grand devant $\tau'$. Pour que la moyenne temporelle approche la moyenne statistique, sans pour autant masquer le phénomène instationnaire étudié, il faut alors que $\bar{\tau}$ soit comprise entre ces deux échelles de temps. Aussi,

$$ \tau' \ll \bar{\tau} \ll \tau \quad (A.8) $$

La moyenne de Reynolds, temporelle, remplaçant la moyenne statistique A.6, est ainsi définie:

$$ \bar{A}(\vec{x}, t) = \frac{1}{\bar{\tau}} \int_{-\frac{1}{2}\bar{\tau}}^{+\frac{1}{2}\bar{\tau}} A(\vec{x}, t+n)dn \quad (A.9) $$

Les variables décomposées en partie fluctuante et partie moyenne, sont introduites dans les équations de Navier-Stokes qui sont ensuite moyennées [10, 53]. Des termes inconnus, dus à la non linéarité des équations de départ, apparaissent alors. Ils expriment la corrélation entre les fluctuations de deux grandeurs, et ils doivent être modélisés pour fermer le système d'équations. Dans le cas d'un écoulement compressible, et afin de faire disparaître les corrélations impliquant des fluctuations de densité, une nouvelle moyenne dite de Favre, peut être utilisée [19]. Elle correspond à une moyenne de Reynolds pondérée par la masse volumique. Soient $\tilde{A}$ la partie moyenne au sens de Favre, et $A$ la partie fluctuante associée. Nous avons alors:

$$ \tilde{A}(\vec{x}, t) = \frac{\overline{\rho A(\vec{x}, t)}}{\overline{\rho}(\vec{x}, t)} \quad (A.10) $$
Si la moyenne de Reynolds A.9 est intuitivement proche de la moyenne statistique A.6, qui devrait être utilisée, il n'en est rien de la moyenne de Favre A.10. Son emploi peut donc être sujet à caution si nous nous réferons à l'esprit de moyenne statistique initial. Cependant, les décompositions de Reynolds et de Favre sont équivalentes lorsque la masse volumique est constante. D'autre part, Favre a démontré que pour un taux de turbulence modéré et un nombre de Mach inférieur à 5, la différence entre les deux approches est du second ordre. De plus, les résultats obtenus par certaines techniques expérimentales, comme l'anémométrie par fil chaud, sont proches des moyennes de Favre, tandis que d'autres méthodes, comme l'anémométrie Doppler laser, fournissent des moyennes de Reynolds. Aussi, nous utiliserons par la suite les deux et nous décomposerons au sens de Favre les variables définies par unité de masse ($\tilde{V}$, $E$, $e$), et au sens de Reynolds les autres ($\rho$, $p$).

## A.4 Moyenne des équations

Les variables ainsi décomposées en partie fluctuante et partie moyenne, sont donc introduites dans les équations de Navier-Stokes qui sont moyennées au sens de Reynolds. Il est à noter que la moyenne de Reynolds, tout comme la moyenne statistique, commute avec les opérateurs différentiels spatiaux et temporels. En utilisant les propriétés suivantes,

$$ \tilde{A} = A \qquad (A.11) $$

$$ \tilde{A}' = 0 \qquad (A.12) $$

$$ \tilde{A}'' = 0 \qquad (A.13) $$

$$ \overline{\rho A''} = 0 \qquad (A.14) $$

$$ \overline{\rho A} = \overline{\rho \tilde{A}} \qquad (A.15) $$

$$ \overline{\rho A \tilde{B}} = \overline{\rho \tilde{A} \tilde{B}} + \overline{\rho A'' \tilde{B}''} \qquad (A.16) $$

il est possible de transformer les équations instantanées A.1, A.2 et A.3, en:

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

$$ \frac{\partial \bar{\rho} \tilde{V}}{\partial t} + \nabla \cdot [\bar{\rho} \tilde{V} \otimes \tilde{V} + \bar{p} \vec{I} - \{\vec{\tau}_l - \overline{\rho \tilde{V}'' \otimes \tilde{V}''}\}] = 0 \qquad (A.18) $$

$$ \frac{\partial \bar{p} \tilde{E}}{\partial t} + \nabla \cdot \left[ (\bar{p} \tilde{E} + \bar{p}) \tilde{V} - \left\{ \frac{\gamma \rho}{P_r} \tilde{\nabla} e + \vec{\tau}_l \cdot \tilde{V} - \overline{(\rho E + p) \tilde{V}''} \right\} \right] = 0 \qquad (A.19) $$

avec:
$$ \bar{p} = \bar{p}(\gamma - 1) \left[ \tilde{E} - \frac{1}{2} \left\{ \tilde{\vec{V}}^2 + \tilde{V}^{\text{n}} \right\} \right] \quad (A.20) $$

$$ \tilde{\tau}_t = \lambda \tilde{\nabla} \cdot \tilde{V} \vec{I} + \bar{\mu} (\tilde{\nabla} \tilde{V} + \tilde{\nabla} \tilde{V}^t) \quad (A.21) $$

Pour obtenir entre autres cette dernière équation, nous avons fait l'hypothèse que $\lambda$ et $\mu$ sont proportionnelles à la densité $\rho$ et à la viscosité cinématique $\nu$, qui est supposée non fluctuante, de même que $\gamma$ et $P_r$. Si bien que:

$$ \bar{\mu} = \bar{p}\bar{\nu} = \bar{p}\nu \quad (A.22) $$

## A.5 Simplification du système d'équations

Le système ainsi obtenu, formé des équations A.17, A.18 et A.19, a pour particularité de contenir plus d'inconnues que d'équations. Il faut donc le simplifier, ou introduire de nouvelles équations pour faciliter sa fermeture. Voici comment nous résolvons chaque problème rencontré.

### Energie cinétique turbulente

L'expression de la pression moyenne A.20 contient un terme supplémentaire trauisant la corrélation des fluctuations de vitesse avec elles-mêmes. Soit $k$, l'énergie cinétique turbulente instantanée, définie par:

$$ k = \frac{1}{2} \tilde{V}^{\text{n}} \quad (A.23) $$

Nous avons alors:

$$ \bar{p} = \bar{p}(\gamma - 1) \left[ \tilde{E} - \frac{1}{2} \tilde{\vec{V}}^2 - \tilde{k} \right] \quad (A.24) $$

### Tenseur de Reynolds

Les inconnues de la forme $-\rho\tilde{V}'' \otimes \tilde{V}''$ constituent le tenseur de Reynolds. Elles traduisent la présence de flux turbulents de quantité de mouvement (équation A.18). Il est possible d'écrire les équations d'évolution de telles grandeurs. Cependant, apparaissent alors des corrélations d'ordre supérieur. Le système ainsi obtenu n'est toujours pas fermé et il est d'une complexité accrue. C'est pourquoi nous utiliserons une approche plus simple, fondée sur l'hypothèse de Boussinesq [10]. Dite de type premier gradient, elle suppose une dépendance linéaire entre le tenseur moyen du gradient de vitesse, et le tenseur de Reynolds. Ce dernier est alors formellement proche du tenseur visqueux $\tilde{\tau}_t$. Par analogie, deux viscosités turbulentes, $\lambda_t$ et $\mu_t$, sont introduites. Ainsi, nous posons:
$$ - \overline{\rho \vec{V}'' \otimes \vec{V}''} = \lambda_t \overline{\vec{\nabla} \cdot \vec{V}} \vec{I} + \mu_t (\overline{\vec{\nabla} \vec{V}} + \overline{\vec{\nabla} \vec{V}''}) - \frac{2}{3} \overline{\rho k} \vec{I} \quad (A.25) $$

Le dernier membre de cette expression assure que la trace du tenseur de Reynolds est proportionnelle à l'énergie cinétique turbulente moyenne, si $\lambda_t$ et $\mu_t$ sont reliées par l'hypothèse de Stokes.

## Equation de l'énergie

Le terme non linéaire $-(\rho E + p)\vec{V}''$ apparaissant dans l'équation A.19 introduit également des corrélations, qui s'écrivent de la façon suivante:

$$ - (\overline{\rho E + p}) \overline{u_i''} = -\gamma C_v \overline{\rho} \overline{T u_i''} + \sum_{j=1}^{3} ( -\overline{\rho u_i'' u_j''} \overline{u_j''} - \overline{\rho k} \overline{u_i''} ) \quad (A.26) $$

Les corrélations température/vitesse $\overline{Tu_i''}$ sont classiquement exprimées en fonction du gradient moyen de température par l'intermédiaire d'un nombre de Prandtl turbulent $P_{rt}$ comme suit [13]:

$$ \overline{\rho T \vec{V}} = - \frac{\mu_t}{P_{rt}} \overline{\nabla T} \quad (A.27) $$

Si bien que:

$$ \gamma C_v \overline{\rho} \overline{T \vec{V}} = - \frac{\gamma \mu_t}{P_{rt}} \overline{\nabla e} \quad (A.28) $$

La même approche est suivie pour les corrélations triples $\overline{ku_i''}$. D'où, $\sigma_k$ étant une constante,

$$ \overline{\rho k \vec{V}} = - \frac{\mu_t}{\sigma_k} \overline{\nabla k} \quad (A.29) $$

Le terme $-\overline{\rho u_i'' u_j''}$ est quant à lui un élément du tenseur de Reynolds, précédemment modélisé avec l'hypothèse de Boussinesq A.25.

Enfin, le développement de $\overline{\overline{\vec{r}_l} \cdot \vec{V}}$ laisse apparaître des corrélations $\overline{\overline{\vec{r}_l} \cdot \vec{V}''}$, qui expriment le travail moyen des contraintes visqueuses lors du déplacement engendré par les fluctuations de vitesse. Cet effet de dissipation est faible, et il sera par la suite négligé. Il ne reste alors que:

$$ \overline{\overline{\overline{r}_l} \cdot \vec{V}} = \overline{\overline{\overline{r}_l} \cdot \vec{V}} \quad (A.30) $$

## Moyenne de Favre et opérateurs différentiels

Contrairement à la moyenne statistique et à la moyenne de Reynolds, la décomposition de Favre ne commute par avec les opérateurs différentiels spatiaux. Aussi, $\overline{\overline{\nabla A}}$ sera différent de $\overline{\overline{\nabla A}}$. Il est en effet possible d'établir que:
$$ \overline{\overline{\nabla} A} = \frac{1}{\bar{\rho}} \left[ \overline{\nabla} (\bar{\rho} \tilde{A}) - \overline{A} \overline{\nabla} \rho \right] \quad (A.31) $$

alors que:

$$ \overline{\nabla} \tilde{A} = \frac{1}{\bar{\rho}} \left[ \overline{\nabla} (\bar{\rho} \tilde{A}) - \tilde{A} \overline{\nabla} \bar{\rho} \right] \quad (A.32) $$

Ces deux relations sont mathématiquement équivalentes à condition que soit vérifié le théorème du cancre énonçant que la moyenne du produit $\overline{A \tilde{\nabla} \rho}$ est égale au produit des moyennes $\tilde{A} \overline{\nabla} \bar{\rho}$. Cependant, des arguments physiques peuvent justifier l'identification de $\overline{\tilde{A} A}$ à $\overline{\tilde{A} \tilde{A}}$. Tout d'abord, les effets turbulents s'ajoutent aux termes diffusifs contenant de telles expressions. Or, le niveau d'approximation dont est entachée la modélisation de la turbulence est plus important que l'erreur commise alors. D'autre part, ces différentielles spatiales ne sont d'un ordre de grandeur sensible que dans des zones à faibles nombres de Mach, telles que les couches limites ou les sillages. Dans ces régions faiblement compressibles, les moyennes de Reynolds et de Favre se confondent, faisant ainsi disparaître toute ambiguïté mathématique. C'est pourquoi $\overline{\tilde{A} \tilde{A}}$ remplacera $\overline{\tilde{A} A}$ par la suite.

## A.6 Conclusion

Finalement, 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 \quad (A.33) $$

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

$$ \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 \quad (A.35) $$

avec:

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

$$ \overline{\vec{\tau}_{l+t}} = (\bar{\lambda} + \lambda_t) \bar{\nabla} \cdot \tilde{V} \vec{I} + (\bar{\mu} + \mu_t) (\tilde{\nabla} \tilde{V} + \tilde{\nabla} \tilde{V}^t) - \frac{2}{3} \bar{\rho} \tilde{k} \vec{I} \quad (A.37) $$

Trois variables supplémentaires, l'énergie cinétique turbulente moyenne $\tilde{k}$ et les 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 (A.38)$$
# Annexe B

## Adimensionnalisation des équations de Navier-Stokes

Le système d'équations différentielles fermé, formant les équations de Navier-Stokes moyennées, doit être adimensionné. Cette opération consiste à introduire des grandeurs de référence, afin de faire apparaître des variables réduites sans dimension. L'intérêt de cette manipulation est double. D'une part, sur le plan physique, elle permet de définir le principe de similitude, qui nous assure qu'une maquette à échelle réduite se comportera en soufflerie sensiblement comme comme un avion de taille réel en vol. Il est ainsi possible de classer le comportement des écoulements en fonction d'un certain nombre de paramètres, dits de similitude, tels que le nombre de Reynolds qui caractérise l'importance des effets convectifs par rapport aux effets diffusifs, ou encore le nombre de Prandtl, comparant la diffusion à la conduction thermique. D'autre part, sur le plan numérique, nous sommes confrontés, dans le système d'unités international, à une grande disparité d'ordres de grandeur des variables physiques, puisque traditionnellement la masse volumique est de quelques kilogrammes par mètre cube, alors que la pression statique est de quelques cent milles pascals. Or un adimensionnement judicieux permet de ramener les variables principales à un ordre de grandeur proche de l'unité. Une meilleure gestion de la précision finie des calculs sur ordinateur est ainsi rendue possible.

Dans ce but, nous définissons cinq grandeurs caratéristiques, en sus du rapport des chaleurs spécifiques $\gamma$ et du nombre de Prandtl $P_r$ déjà introduits :

* **Echelle géométrique**

Soit $L_0$ la longueur de référence.

* **Echelles cinématiques**

Soient $\rho_0$ et $V_0$, la masse volumique et la vitesse de référence. Il est possible de définir la pression de référence $P_0$ à la place de l'une de ces deux grandeurs, sachant que :
$$P_0 = \rho_0 V_0^2 \qquad (B.1)$$

## • Echelles thermodynamiques

Soient $r_0$ et $\mu_0$, la constante du gaz supposé parfait et la viscosité dynamique de référence. Il est possible de définir le nombre de Reynolds à la place de cette dernière par:

$$R_e = \frac{\rho_0 L_0 V_0}{\mu_0} \qquad (B.2)$$

Il est à noter qu'une conductibilité thermique de référence $\kappa_0$ peut être construite à partir du nombre de Prandtl car:

$$P_r = \frac{\gamma \frac{\mu_0 r_0}{\gamma - 1}}{\kappa_0} \qquad (B.3)$$

Ceci nous conduit aux variables sans dimension suivantes:

$$t^* = \frac{t}{\left(\frac{L_0}{V_0}\right)} \qquad (B.4)$$

$$x_i^* = \frac{x_i}{L_0} \qquad (B.5)$$

$$\rho^* = \frac{\bar{\rho}}{\rho_0} \qquad (B.6)$$

$$u_i^* = \frac{\tilde{u}_i}{V_0} \qquad (B.7)$$

$$E^* = \frac{\tilde{E}}{V_0^2} \qquad (B.8)$$

$$p^* = \frac{\bar{p}}{\rho_0 V_0^2} \qquad (B.9) \qquad T^* = \frac{\tilde{T}}{\left(\frac{V_0^2}{r_0}\right)} \qquad (B.10)$$

$$e^* = \frac{\tilde{e}}{V_0^2} \qquad (B.11)$$

$$\mu^* = \frac{\bar{\mu}}{\mu_0} \qquad (B.12) \qquad \tilde{\tau}^* = \frac{\tilde{\tau}}{\left(\frac{\mu_0 V_0}{L_0}\right)} \qquad (B.13)$$

Introduites dans les équations de Navier-Stokes moyennées et fermées, elles nous permettent d'écrire le système suivant dans lequel l'astérix désignant les grandeurs adimensionnées n'a pas été porté dans un soucis de lisibilité.

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

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

$$\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 (B.16)$$
avec:

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

$$\vec{\tau}_{l+t} = (\lambda + \lambda_t) \vec{\nabla} \cdot \vec{V} \vec{I} + (\mu + \mu_t) (\vec{\nabla} \vec{V} + \vec{V} \vec{\nabla}) \qquad (B.18)$$
# Annexe C

## Formulation volumes finis

### C.1 Introduction

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 est décomposée par la suite en trois étapes. Les équations différentielles sont tout d'abord intégrées en espace sur un maillage mobile. Puis, un changement de variables est réalisé, permettant de passer des variables spatiales cartésiennes $(x_1, x_2, x_3, t)$ aux variables curvilignes $(\xi^1, \xi^2, \xi^3, t)$. Enfin, un certain nombre d'hypothèses conduisent à la forme discrète recherchée des équations. Une comparaison au résultat obtenu dans le cadre d'une formulation différences finies clôture ces développements.

### C.2 Formulation intégrale en maillage mobile

Tout d'abord, 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 (C.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.

#### C.2.1 Cas d'une équation scalaire

$$ \frac{\partial A}{\partial t} + \nabla \cdot \vec{B} = 0 \qquad (C.2) $$
<figure><img src="image_162.png" /><figcaption>Figure C.1: Cellule d'intégration</figcaption></figure>

devient, par intégration sur la cellule $\Omega(t)$ (figure C.1),

$$ \int_{\Omega(t)} \frac{\partial A}{\partial t} d\Omega + \int_{\Omega(t)} \bar{\nabla} \cdot \vec{B} d\Omega = 0 \quad (C.3) $$

Le théorème de Leibnitz permet d'établir que:

$$ \int_{\Omega(t)} \frac{\partial A}{\partial t} d\Omega = \frac{d}{dt} \int_{\Omega(t)} A d\Omega - \oint_{S(t)} A(\vec{W} \cdot \vec{n}) dS \quad (C.4) $$

où:

* la dérivée temporelle du membre de gauche est calculée en un point fixe du repère cartésien.

* la dérivée temporelle du membre de droite est calculée en suivant le mouvement éventuel de la cellule d'intégration.

* $S(t)$ est la surface fermée limitant $\Omega(t)$.

* $\vec{n}$ est un vecteur unitaire perpendiculaire à $S(t)$ orienté vers l'extérieur de la cellule d'intégration.

* $\vec{W}$ est la vitesse de $S(t)$ exprimée dans le repère cartésien. $\vec{W} = \partial OM/\partial t$ où M est un point de $S(t)$.

D'autre part, le théorème de Gauss donne:

$$ \int_{\Omega(t)} \bar{\nabla} \cdot \vec{B} d\Omega = \oint_{S(t)} \vec{B} \cdot \vec{n} dS \quad (C.5) $$
Finalement, nous obtenons:

$$ \frac{d}{dt} \int_{\Omega(t)} Ad\Omega + \oint_{S(t)} [\vec{B} - A\vec{W}] \cdot \vec{n} dS = 0 \quad (C.6) $$

### C.2.2 Cas d'une équation vectorielle

De la même façon, l'équation :

$$ \frac{\partial \vec{A}}{\partial t} + \vec{\nabla} \cdot \vec{B} = \vec{0} \quad (C.7) $$

se transforme en:

$$ \frac{d}{dt} \int_{\Omega(t)} \vec{A} d\Omega + \oint_{S(t)} [\vec{\nabla} \cdot \vec{n} - \vec{A}(\vec{W} \cdot \vec{n})] dS = \vec{0} \quad (C.8) $$

En remarquant que :

$$ \vec{A}(\vec{W} \cdot \vec{n}) \equiv A \times (W^T \times n) = (A \times W^T) \times n \equiv (\vec{A} \otimes \vec{W}) \cdot \vec{n} \quad (C.9) $$

nous obtenons,

$$ \frac{d}{dt} \int_{\Omega(t)} \vec{A} d\Omega + \oint_{S(t)} [\vec{\nabla} - \vec{A} \otimes \vec{W}] \cdot \vec{n} dS = \vec{0} \quad (C.10) $$

## C.3 Changement de variables d'espace

Il nous faut maintenant évaluer les intégrales apparaissant dans les équations C.6 et C.10, qui ont pour forme générique :

$$ \frac{d}{dt} \int_{\Omega(t)} Ad\Omega + \oint_{S(t)} \vec{C} \cdot \vec{n} dS = 0 \quad (C.11) $$

où :

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

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

Cette équation, exprimée en fonction des variables cartésiennes $(x_1, x_2, x_3, t)$, doit l'être en fonction des variables curvilignes $(\xi^1, \xi^2, \xi^3, t)$. Il est à noter que la base sur laquelle sont projetées les grandeurs vectorielles reste dans les deux cas celle de départ, à savoir $(\vec{e}_1, \vec{e}_2, \vec{e}_3)$. Ceci a pour conséquence que, formellement, une équation vectorielle sera traitée comme une équation scalaire.

Les cellules d'intégration choisies sont des hexaèdres (figure C.2), limités par des surfaces telles que $\xi^i$ soit constant sur $S_i$. Ainsi :
<figure><img src="image_163.png" /><figcaption>Figure C.2: Changement de variables d'espace</figcaption></figure>

$$ \omega = \prod_{j=1}^{3} [\xi^{j-}, \xi^{j+}] \qquad (C.12) $$

$$ S^i = \prod_{\substack{j=1 \\ j \neq i}}^{3} [\xi^{j-}, \xi^{j+}] \qquad (C.13) $$

Par changement de variable, l'intégrale volumique devient:

$$ \frac{d}{dt} \int_{\Omega(t)} Ad\Omega = \frac{d}{dt} \int_{\omega} A \sqrt{g} \prod_{j=1}^{3} d\xi^{j} \qquad (C.14) $$

où $\sqrt{g} = \frac{D(x_1, x_2, x_3)}{D(\xi^2, \xi^2, \xi^3)}$ est le jacobien de la transformation géométrique.

Le maillage curviligne étant choisi fixe en temps, contrairement au maillage cartésien, la vitesse de la surface de $\omega$ est nulle. Nous obtenons alors:

$$ \frac{d}{dt} \int_{\Omega(t)} Ad\Omega = \int_{\omega} \frac{\partial}{\partial t}(\sqrt{g} A) \prod_{j=1}^{3} d\xi^{j} \qquad (C.15) $$

L'intégrale surfacique quant à elle se réécrit:

$$ \oint_{S(t)} \vec{C} \cdot \tilde{n} dS = \sum_{i=1}^{3} \left\{ \int_{S^{i+}} (\vec{C} \cdot \tilde{a}^i) \sqrt{g} \prod_{\substack{j=1 \\ j \neq i}}^{3} d\xi^j \right\} $$
$$ - \sum_{i=1}^{3} \left\{ \int_{S^{i-}} (\vec{C} \cdot \vec{a}^i) \sqrt{g} \prod_{\substack{j=1 \\ j \neq i}}^{3} d\xi^j \right\} \qquad (C.16) $$

où $\vec{a}^i = \vec{\nabla}\xi^i$ est le vecteur de la base contravariante associé à la direction $\xi^i$. Le choix du signe dépend de la position du volume $\omega$ par rapport à la surface $S^{i\pm}$.

Finalement, l'équation intégrale suivante, équivalente à l'équation C.11, est obtenue:

$$ \begin{aligned} & \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\} \\ & - \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\} = 0 \end{aligned} \qquad (C.17) $$

## C.4 Hypothèses des volumes finis

La formulation volumes finis permet d'évaluer numériquement les intégrales volumiques et surfaciques précédentes. Elle peut être exprimée de deux manières différentes.

### C.4.1 Première approche

L'hypothèse principale de la formulation volumes finis est que les intégrandes apparaissant dans l'équation C.17 varient peu sur l'élément d'intégration. Elles peuvent alors être remplacées par une valeur moyenne par élément.

Il en résulte les approximations suivantes:

$$ \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 (C.18) $$

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

Or, dans l'espace curviligne, les cellules d'intégration sont parallélépipédiques, d'où:
$$ \omega = \prod_{j=1}^{3} (\xi^{j+} - \xi^{j-}) \qquad (C.20) $$

$$ S^i = \prod_{\substack{j=1 \\ j \neq i}}^{3} (\xi^{j+} - \xi^{j-}) \Big|_{\xi^i = \text{const}} \qquad (C.21) $$

Si de plus, elles sont choisies cubiques telles que:

$$ \xi^{j\pm} = \xi^j \pm \frac{1}{2} \qquad (C.22) $$

nous avons:

$$ \omega = 1 \qquad (C.23) $$

$$ S^i = 1 \qquad (C.24) $$

Finalement, l'expression approchée de l'équation C.17 est:

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

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

$$ F^i = \sqrt{g} \vec{C} \cdot \vec{a}^i \qquad (C.26) $$

### C.4.2 Deuxième approche

L'approximation volumes finis précédente porte sur des intégrandes incluant des termes géométriques et des grandeurs aérodynamiques. Dans le cas où le maillage cartésien est fixe, il est possible d'obtenir une formulation découplée.

Ainsi:

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

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

Seules les grandeurs aérodynamiques sont moyennées par élément, l'intégrale des termes géométriques pouvant être évaluée indépendamment, par une méthode d'intégration de Gauss par exemple.
<figure><img src="image_164.png" /><figcaption>Figure C.3: Stockage des variables conservatives aux noeuds</figcaption></figure>

### C.4.3 Bilan

Dans le cas d'un maillage mobile, le coût de calcul supplémentaire introduit par la deuxième méthode est important. Il faut en effet évaluer les métriques à chaque pas en temps, puis les intégrer. C'est pourquoi nous suivrons la première approche.

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 C.3). Une autre stratégie aurait consisté à stocker le champ aérodynamique au centre des cellules construites alors sur les noeuds (figure C.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.

## C.5 Equations discrétisées

Nous avons donc établi qu'une équation différentielle ayant pour forme:

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

où:

- A est une grandeur scalaire (resp. vectorielle)

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

pouvait être approchée dans le cadre d'une formulation volumes finis par:
<figure><img src="image_165.png" /><figcaption>Figure C.4: Stockage des variables conservatives au centre des mailles</figcaption></figure>

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

avec:

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

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

## C.6 Comparaison aux différences finies

Les différences finies constituent une autre technique de discrétisation spatiale, au même titre que les volumes finis. Il est connu que ces deux méthodes sont équivalentes en maillage cartésien orthogonal régulier. Nous allons établir par la suite qu'utiliser une formulation volumes finis en espace carté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.

Partant de l'équation intégrale en repère curviligne C.17, nous faisons tendre vers zéro le volume d'intégration. Si bien que:

$$ \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 \prod_{j=1}^{3} \delta\xi^j \quad (C.31) $$

et:
$$ \begin{align} & \int_{S^i+} \sqrt{g} \vec{C} \cdot \vec{a}^i \prod_{\substack{j=1 \\ j \neq i}}^3 d\xi^j - \int_{S^i-} \sqrt{g} \vec{C} \cdot \vec{a}^i \prod_{\substack{j=1 \\ j \neq i}}^3 d\xi^j \nonumber \\ & \approx \int_{S^i} \left\{ [\sqrt{g} \vec{C} \cdot \vec{a}^i]^+ - [\sqrt{g} \vec{C} \cdot \vec{a}^i]^- \right\} \prod_{\substack{j=1 \\ j \neq i}}^3 d\xi^j \tag{C.32} \end{align} $$

$$ \approx \int_{S^i} \left\{ \frac{\partial}{\partial \xi^i} \left[ \sqrt{g} \vec{C} \cdot \vec{a}^i \right] \delta \xi^i \right\} \prod_{\substack{j=1 \\ j \neq i}}^3 d\xi^j \quad (C.33) $$

$$ \approx \frac{\partial}{\partial \xi^i} \left[ \sqrt{g} \vec{C} \cdot \vec{a}^i \right] \prod_{j=1}^{3} \delta \xi^j \quad (C.34) $$

D'où l'équation différentielle suivante, associée à l'équation C.17:

$$ \frac{\partial}{\partial t} (\sqrt{g} A) + \sum_{i=1}^{3} \frac{\partial}{\partial \xi^i} \left[ \sqrt{g} \vec{C} \cdot \vec{a}^i \right] = 0 \quad (C.35) $$

Dans une approche différences finies, les dérivées spatiales peuvent être discrétisées de manière centrée autour de $\xi^i$, si bien que:

$$ \frac{\partial}{\partial \xi^i} \left[ \sqrt{g} \vec{C} \cdot \vec{a}^i \right] \approx \frac{1}{\Delta \xi^i} \left[ \left. \sqrt{g} \vec{C} \cdot \vec{a}^i \right|_{\xi^i + \frac{1}{2}} - \left. \sqrt{g} \vec{C} \cdot \vec{a}^i \right|_{\xi^i - \frac{1}{2}} \right] \quad (C.36) $$

$$ \approx \sqrt{g} \vec{C} \cdot \vec{a}^i \left|_{\xi^i + \frac{1}{2}} - \sqrt{g} \vec{C} \cdot \vec{a}^i \right|_{\xi^i - \frac{1}{2}} \quad (C.37) $$

car $\Delta \xi^i = 1$, et donc:

$$ \sum_{i=1}^{3} \frac{\partial}{\partial \xi^i} \left[ \sqrt{g} \vec{C} \cdot \vec{a}^i \right] \approx \sum_{i=1}^{3} \left\{ \pm \sqrt{g} \vec{C} \cdot \vec{a}^i \left|_{\xi^i \pm \frac{1}{2}} \right. \right\} \quad (C.38) $$

Nous obtenons ainsi la même équation discrète que lors de l'approche volumes finis (équation C.25) à partir de l'équation C.35:

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

Ce qui établit la similitude entre les deux méthodes. Cette dernière s'explique par la forme hexaèdrique des cellules d'intégration et par le choix particulier des valeurs de $\xi^j\pm$, qui fixe $\omega$ et $S^i$ à 1.
# Annexe D

## Géométrie différentielle en coordonnées curvilignes

### D.1 Introduction

Les développements présentés par la suite visent à l'expression de certaines grandes géométriques, telles que le gradient d'un scalaire normal à une surface ou la divergence d'un vecteur, dans un système de coordonnées curvilignes [63]. En effet, un point du domaine de calcul spatio-temporel peut être repéré, soit par ses coordonnées catésiennes et une date, ($x_1, x_2, x_3, t$), soit par ses coordonnées curvilignes et une date, ($\xi^1, \xi^2, \xi^3, t$) (figure D.1). Il est à noter que dans les deux cas, tout vecteur est projeté sur la base catésienne ($\tilde{e}_1, \tilde{e}_2, \tilde{e}_3$), qui est indépendante de ce point et du temps. Or, dans le cadre d'une discrétisation de l'espace au moyen d'un maillage structuré, il est intéressant de formuler les éléments différentiels et les opérateurs différentiels, classiquement connus en fonction des coordonnées cartésiennes, à partir d'informations géométriques dérivées des coordonnées curvilignes. Pour ce faire, deux nouvelles bases vectorielles locales, dites covariante et contravariante, doivent cependant être introduites au préalable.

### D.2 Base covariante

Le vecteur $\tilde{a}_i$, de la base covariante, associé à la direction $\xi^i$, est défini par:

$$ \tilde{a}_i = \frac{\partial \vec{\tau}}{\partial \xi^i} \qquad (D.1) $$

Il a pour caractéristique d'être tangent en $\vec{\tau}$ à la ligne "$\xi^i$ variable/$\xi^j$ et $\xi^k$ constants" (figure D.2).

Les vecteurs de la base covariante ont pour propriété remarquable que:
<figure><img src="image_166.png" /><figcaption>Figure D.1: Coordonnées cartésiennes et curvilignes</figcaption></figure>

<figure><img src="image_167.png" /><figcaption>Figure D.2: Base covariante</figcaption></figure>