# Au sujet du Gaussien Anisotrope

Adam Larat

► To cite this version:

Adam Larat. Au sujet du Gaussien Anisotrope. [Rapport de recherche] CNRS. 2016. hal-01272438

HAL Id: hal-01272438

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

Submitted on 15 Feb 2016

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

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

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

Distributed under a Creative Commons CC BY 4.0 - Attribution - International License
# Au sujet du Gaussien Anisotrope

Adam Larat

15 février 2016

## 1 Introduction Générale

Dans ce document, on considère un nuage de particules solides transporté par un support gazeux. Pour le moment, on ne s'intéresse pas au problème de la distribution en taille de ces particules. On cherche une description purement statistique du problème et l'on considère pour cela la fonction de répartition $f(t, \vec{x}, \vec{c})$ du nombre probable de particule se déplaçant à la vitesse $\vec{c}$, à la position $\vec{x}$ et au temps $t$. L'équation de Williams-Boltzmann régit le transport de cette fonction de répartition dans l'espace des phase $(t, \vec{x}, \vec{c})$ :

$$ \partial_t f + \partial_{\vec{x}} (\vec{c} f) + \partial_{\vec{c}} (\vec{F} f) = 0 \qquad (1) $$

Les deux premiers termes expriment le transport des particules à leur vitesse respective, tandis que le troisième terme exprime un couplage avec le champ de vitesse du gaz sous-jacent. En fonction de son inertie, la particule relaxe plus ou moins vite vers la vitesse du gaz. Ce phénomène est caractérisé par le temps de relaxation $\tau_p$ et la force de couplage peut être modélisée sous sa forme la plus simple comme :

$$ \vec{F} = (\vec{u}_g - \vec{c})/\tau_p. \qquad (2) $$

Afin d'éviter la simulation numérique de cette équation dans tout l'espace des phases (de taille $2d + 1$, $d$ étant la dimension en espace, et potentiellement de taille infinie dans les directions de vitesse), on souhaite dorénavant ne transporter que les premiers moments en vitesse de la fonction de distribution en nombre :

$$ \mathcal{M}_k = \int_{\mathbb{R}^d} (\otimes^k \vec{c}) f(t, \vec{x}, \vec{c}) d\vec{c}, \qquad (3) $$

où l'opérateur $\otimes^k$ représente la puissance $k$-ième du produit tensoriel,

$$ \otimes^0 \vec{c} = 1, \quad \otimes^1 \vec{c} = \vec{c}, \quad \otimes^2 \vec{c} = \vec{c} \otimes \vec{c}, \quad \text{etc...} $$

En multipliant (1) par $\otimes^k \vec{c}$ et en intégrant sur tout l'espace des vitesses, on obtient l'équation d'évolution des moments d'ordre exactement $k$, représentée sous forme tensorielle :

$$ \partial_t \mathcal{M}_k + \int_{\mathbb{R}^d} (\otimes^k \vec{c}) \partial_{\vec{x}} (\vec{c} f) d\vec{c} = - \int_{\mathbb{R}^d} (\otimes^k \vec{c}) \partial_{\vec{c}} (\vec{F} f) d\vec{c}. \qquad (4) $$
Par symétrie du tenseur obtenu par l'opérateur $\otimes^k$ (produit tensoriel multiple d'un même vecteur), on montre aisément que l'intégrande du second terme devient

$$ \partial_{\vec{x}} (\otimes^{k+1} \vec{c} f). $$

On peut alors permuter l'intégrale et la dérivée partielle en espace pour faire ressortir le tenseur des moments d'ordre supérieur. De même, on peut regarder le terme $(i_1, \dots, i_k)$ du membre de droite. En intégrant par partie et en faisant bien attention aux indices (ou certainement en invoquant un puissant théorème d'intégration par partie dans l'espace des tenseurs symétriques), on peut montrer, voir Annexe A, que l'intégrande devient au signe près :

$$ k f \left[ (\otimes^{k-1} \vec{c}) \odot \vec{F} \right]. $$

L'opérateur $\odot$ désigne ici le produit tensoriel symétrique : si $\mathbb{B}$ et $\mathbb{D}$ sont des tenseurs d'ordres respectifs $p$ et $q$, avec $p + q = n$, alors

$$ (\mathbb{B} \odot \mathbb{D})_{(i_1, \dots, i_n)} = \frac{1}{n!} \sum_{\sigma \in S(n)} (\mathbb{B} \otimes \mathbb{D})_{\sigma(i_1, \dots, i_n)}. \qquad (5) $$

Dans ce contexte, d'un point de vue très général, le transport du tenseur des $C_{d+k-1}^k$ (voir Annexe B) moments d'ordre $k$ s'écrit :

$$ \partial_t \mathcal{M}_k + \partial_{\vec{x}} \mathcal{M}_{k+1} = \frac{k (\mathcal{M}_{k-1} \odot \vec{u}_g - \mathcal{M}_k)}{\tau_p}. \qquad (6) $$

Puisque l'équation régissant l'évolution des moments d'ordre $k$ fait intervenir le tenseur des moments d'ordre $k+1$, il s'agit maintenant d'apporter une fermeture adéquate, afin de ne transporter qu'un nombre fini de moments. Une manière de faire consiste en une hypothèse sur la forme de la fonction de distribution en vitesse. Par exemple, si la distribution dans l'espace est supposée être un Dirac, c'est à dire qu'on associe une unique vitesse à chaque point de l'espace et du temps, $f(t, \vec{x}, \vec{c}) = \rho(t, \vec{x})\delta(\vec{c} - \vec{u}(t, \vec{x}))$, les équations sur les deux premiers moments deviennent :

$$ \left\{ \begin{array}{l} \partial_t \rho + \partial_{\vec{x}} (\rho \vec{u}) = 0, \\ \\ \partial_t (\rho \vec{u}) + \partial_{\vec{x}} (\rho \vec{u} \otimes \vec{u}) = \frac{\rho \vec{u}_g - \rho \vec{u}}{\tau_p}, \end{array} \right. \qquad (7) $$

système qui est parfaitement clos. On remarquera que pour être parfaitement définie, la fonction de distribution nécessite la connaissance de $d+1$ quantités ($\rho$ et les $d$ composantes de la vitesse) en espace et en temps. Ces quantités sont simplement les solutions des $d+1$ équations de transport des deux premiers ordres de moments.

On peut chercher à aller plus loin et ne fermer le système de transport qu'à l'ordre 2. On considère alors $\frac{(d+1)(d+2)}{2}$ équations sur les moments. Par analogie avec ce qui vient d'être dit, la fonction de distribution choisie devra présenter autant de degrés de liberté. Une possibilité consiste à choisir une distribution Gaussienne et anisotrope en vitesse :

$$ f(t, \vec{x}, \vec{c}) = \frac{\rho(t, \vec{x})}{\sqrt{(2\pi)^d \det \Sigma}} \exp \left( -\frac{1}{2} (\vec{c} - \vec{u}(t, \vec{x}))^T \Sigma^{-1} (\vec{c} - \vec{u}(t, \vec{x})) \right), \qquad (8) $$

où $\Sigma$ est un tenseur symétrique de covariance.
# 2 Système d'équation pour la distribution Gaussienne anisotrope

## 2.1 Moments de la distribution Gaussienne

On commence par donner un rapide rappel sur le calcul des différents moments de la fonction de répartition Gaussienne. En ce qui concerne la moyenne, autrement dit le moment d'ordre zéro, il suffit de remarquer que la fonction de répartition Gaussienne appartient à $\mathcal{S}(\mathbb{R}^{2d+1})$, tous les théorèmes sur les intégrales s'appliquent donc :

$$ 
\begin{aligned}
\left(\int_{-\infty}^{+\infty} e^{-X^2} dX\right)^2 &= \left(\int_{-\infty}^{+\infty} e^{-X^2} dX\right) \left(\int_{-\infty}^{+\infty} e^{-Y^2} dY\right) = \int_{\mathbb{R}^2} e^{-(X^2+Y^2)} dX dY \\
&= \int_{\rho>0} \int_{\theta=0}^{2\pi} \rho e^{-\rho^2} d\rho d\theta = 2\pi \int_{\rho>0} \left(\frac{1}{2} e^{-\rho^2}\right)' d\rho = \pi.
\end{aligned}
$$

D'où

$$ \int_{-\infty}^{+\infty} e^{-X^2} dX = \sqrt{\pi} \quad \text{et} \quad \int_{\mathbb{R}} e^{-\frac{x^2}{\alpha}} dx = \sqrt{\alpha\pi} \qquad (9) $$

On peut maintenant calculer aisément le moment d'ordre zéro de notre répartition anisotrope (8) :

$$ 
\begin{aligned}
\mathcal{M}_0 &= \int_{\mathbb{R}^d} f(t, \vec{x}, \vec{c}) d\vec{c} = \frac{\rho(t, \vec{x})}{\sqrt{(2\pi)^d \det \Sigma}} \int_{\mathbb{R}^d} \exp \left( -\frac{1}{2} \left[ \Sigma^{-1/2} (\vec{c} - \vec{u}) \right]^T \left[ \Sigma^{-1/2} (\vec{c} - \vec{u}) \right] \right) d\vec{c} \\
&= \frac{\rho(t, \vec{x})}{(2\pi)^{d/2}} \int_{\mathbb{R}^d} \exp \left( -\frac{\|\vec{X}\|^2}{2} \right) d\vec{X} \\
&= \rho(t, \vec{x}).
\end{aligned}
$$

Ensuite, il faut remarquer que par antisymétrie, tous les moments centrés d'ordre impair de la Gaussienne sont nuls :

$$ \int_{\mathbb{R}} x^{2k+1} e^{x^2} dx = 0. \qquad (10) $$

Ainsi, le moment d'ordre 1 de $f$ s'obtient facilement :

$$ 
\begin{aligned}
\int_{\mathbb{R}^d} \vec{c} f(t, \vec{x}, \vec{c}) d\vec{c} &= \mathcal{M}_0 \vec{u} + \frac{\rho(t, \vec{x})}{\sqrt{(2\pi)^d \det \Sigma}} \int_{\mathbb{R}^d} (\vec{c} - \vec{u}) \exp \left( -\frac{1}{2} \|\Sigma^{-1/2} (\vec{c} - \vec{u})\|^2 \right) d\vec{c} \\
&= \rho \vec{u}(t, \vec{x})
\end{aligned}
$$

En ce qui concerne les moments d'ordre deux, c'est un peu plus complexe. Je me demande si on peut faire plus simple. On souhaite évaluer

$$ \mathcal{M}_2 = \frac{\rho(t, \vec{x})}{\sqrt{(2\pi)^d \det \Sigma}} \int_{\mathbb{R}^d} \vec{c} \otimes \vec{c} \exp \left( -\frac{1}{2} \|\Sigma^{-1/2} (\vec{c} - \vec{u})\|^2 \right) d\vec{c}. \qquad (11) $$
Pour cela, comme précédemment, on injecte la vitesse moyenne dans le produit tensoriel des vitesses :

$$ \vec{c} \otimes \vec{c} = (\vec{c} - \vec{u}) \otimes (\vec{c} - \vec{u}) + 2\vec{u} \odot (\vec{c} - \vec{u}) + \vec{u} \otimes \vec{u}. $$

Par antisymétrie, le terme du milieu donnera une intégrale nulle. La seule difficulté réside dans l'évaluation de l'intégrale associée au premier des trois termes. Considérons le terme $(i, j) \in [1, d]^2$ de ce tenseur :

$$ \left( \int_{\mathbb{R}^d} \left[ \Sigma^{1/2} \vec{X} \otimes \Sigma^{1/2} \vec{X} \right] \exp \left( -\frac{\|\vec{X}\|^2}{2} \right) d\vec{X} \right)_{ij} = \int_{\mathbb{R}^d} \left( \Sigma^{1/2} \vec{X} \right)_i \left( \Sigma^{1/2} \vec{X} \right)_j \exp \left( -\frac{\|\vec{X}\|^2}{2} \right) d\vec{X} \\ = \int_{\mathbb{R}^d} \left[ \left( \Sigma^{1/2} \right)_{ik_1} X^{k_1} \right] \left[ \left( \Sigma^{1/2} \right)_{jk_2} X^{k_2} \right] \exp \left( -\frac{X_k X^k}{2} \right) d\vec{X} $$

Toujours par l'argument d'antisymétrie, tous les termes qui font intervenir des indices muets $k_1$ et $k_2$ différents sont nuls. Il ne reste donc plus que les termes de la forme

$$ \begin{aligned} & \sum_{k=0}^{d} \left\{ \int_{\mathbb{R}^d} (\Sigma^{1/2})_{ik} (\Sigma^{1/2})_{jk} (X^k)^2 e^{-X_k^2/2} e^{-(X_{j \neq k} X^{j \neq k})/2} d\vec{X} \right\} \\ &= \sum_{k=0}^{d} \left\{ (\Sigma^{1/2})_{ik} (\Sigma^{1/2})_{jk} \sqrt{2\pi} \int_{\mathbb{R}^{d-1}} e^{-(X_{j \neq k} X^{j \neq k})/2} \right\} \\ &= (2\pi)^{d/2} \sum_{ij}. \end{aligned} $$

Attention ! Il est extrêmement important de remarquer que le passage à la dernière ligne n'est possible que sous une hypothèse de symétrie sur $\Sigma^{1/2}$. En général, la matrice de covariance $\Sigma$ est symétrique définie positive et possède donc une unique racine symétrique définie positive. Une étude plus poussée dans le cas général est nécessaire, notamment afin de poursuivre la remarque suggérée plus loin dans la sous-section 2.2. On peut enfin conclure au sujet du moment d'ordre deux, la dernière égalité faisant office de définition :

$$ M_2 = \rho \vec{u} \otimes \vec{u} + \rho \Sigma = 2\rho \mathbb{E} \quad (12) $$

Enfin, nous aurons besoin du moment d'ordre trois de cette distribution anisotrope pour fermer le système. Toujours en effectuant le changement de variable $\Sigma^{1/2} \vec{X} = \vec{c} - \vec{u}$, on a :

$$ \int_{\mathbb{R}^d} (\vec{c} \otimes \vec{c} \otimes \vec{c}) \, f d\vec{c} = \frac{\rho}{(2\pi)^{d/2}} \int_{\mathbb{R}^d} \otimes^3 (\Sigma^{1/2} \vec{X} + \vec{u}) e^{-\frac{\|\vec{X}\|^2}{2}} d\vec{X}. \quad (13) $$

Le développement de la puissance troisième du produit tensoriel donne après un petit calcul

$$ \otimes^3 (\Sigma^{1/2} \vec{X} + \vec{u}) = \otimes^3 (\Sigma^{1/2} \vec{X}) + 3 (\Sigma^{1/2} \vec{X} \otimes \Sigma^{1/2} \vec{X}) \odot \vec{u} + 3 (\vec{u} \otimes \vec{u}) \odot (\Sigma^{1/2} \vec{X}) + \otimes^3 \vec{u} $$

Le produit tensoriel $\otimes$ et le produit tensoriel symétrique $\odot$ étant des applications bili- néaires sur l'espace de tenseurs, on a :

$$ M_3 = \frac{\rho(t, \vec{x})}{(2\pi)^{d/2}} \left\{ \begin{aligned} & \int_{\mathbb{R}^d} \otimes^3 (\Sigma^{1/2} \vec{X}) e^{-\frac{\|\vec{X}\|^2}{2}} d\vec{X} + 3\vec{u} \odot \int_{\mathbb{R}^d} (\Sigma^{1/2} \vec{X} \otimes \Sigma^{1/2} \vec{X}) e^{-\frac{\|\vec{X}\|^2}{2}} d\vec{X} \\ & + 3(\vec{u} \otimes \vec{u}) \odot \int_{\mathbb{R}^d} \Sigma^{1/2} \vec{X} e^{-\frac{\|\vec{X}\|^2}{2}} d\vec{X} + (2\pi)^{d/2} \otimes^3 \vec{u} \end{aligned} \right\} $$
La première intégrale est nulle car dans chacun des termes du tenseur d'ordre trois il y a une composante de $\vec{X}$ qui est à un ordre impair. La seconde intégrale est connue et vaut $\Sigma$. La troisième est nulle par asymétrie. Finalement, on a donc :

$$ M_3 = \rho(t, \vec{x}) (3\vec{u} \odot \Sigma + \otimes^3 \vec{u}) = (3\rho\Sigma + \rho\vec{u} \otimes \vec{u}) \odot \vec{u} = 2(\rho E + P) \odot \vec{u}, \quad (14) $$

où l'on a défini le tenseur de pression comme

$$ P = \rho\Sigma \quad (15) $$

## 2.2 Système à trois moments

D'après ce qui vient d'être démontré, on peut concaténer tous les résultats pour montrer que le transport des trois premiers moments de l'équation (1) donne le système suivant :

$$ \left\{ \begin{array}{l} \partial_t \rho + \partial_{\vec{x}} (\rho \vec{u}) = 0 \\[1em] \partial_t (\rho \vec{u}) + \partial_{\vec{x}} (\rho \vec{u} \otimes \vec{u} + P) = \rho \frac{\vec{u}_g - \vec{u}}{\tau_p} \\[1em] \partial_t (\rho E) + \partial_{\vec{x}} ((\rho E + P) \odot \vec{u}) = \frac{\rho \vec{u}_g \odot \vec{u} - 2\rho E}{\tau_p} \end{array} \right. \quad (16) $$

Grâce à la définition (12) et la "loi d'état" (15), ce système est fermé pour les inconues tensorielles $\rho$, $\vec{u}$ et $\Sigma$.

**Remarque Hypothétique :** Les calculs précédents ont été réalisés sous l'hypothèse que $\Sigma$ est un tenseur symétrique. Cette hypothèse est valide si les particules du brouillard sont supposées ponctuelles. En revanche, lorsque ces particules prennent une taille finie non nulle, elles acquièrent $\frac{d(d-1)}{2}$ moments inertiels en rotation supplémentaires. Un transfert d'énergie du moment cinétique vers le moment angulaire est alors envisageable lors des collisions. La trâce de ce transfert croisé de l'énergie cinétique vers l'énergie de rotation apparaîtrait dans la structure non nécessairement symétrique du tenseur de covariance $\Sigma$. Cette configuration doit être étudiée plus précisément.

## 2.3 Étude dans le cas bi-dimensionnel

Dorénavant, on fixe $d = 2$. Si on note

$$ E = \begin{pmatrix} e_{11} & e_{12} \\ e_{12} & e_{22} \end{pmatrix} \quad \text{et} \quad P = \begin{pmatrix} p_{11} & p_{12} \\ p_{12} & p_{22} \end{pmatrix}, \quad (17) $$

le système (16) devient

$$ \partial_t \mathbf{W} + \partial_x \mathbf{F}(\mathbf{W}) + \partial_y \mathbf{G}(\mathbf{W}) = \mathbf{S}(\mathbf{W}), \quad (18) $$
avec

$$
\mathbf{W} = \begin{pmatrix} \rho \\ \rho u \\ \rho v \\ \rho e_{11} \\ \rho e_{12} \\ \rho e_{22} \end{pmatrix}
$$

et

$$
\vec{\mathcal{F}}(\mathbf{W}) = (\mathbf{F}(\mathbf{W}), \mathbf{G}(\mathbf{W})) = \left\{ \begin{pmatrix} \rho u \\ \rho u^2 + p_{11} \\ \rho u v + p_{12} \\ \frac{1}{2}(\rho u^3 + 3p_{11}u) \\ \frac{1}{2}(\rho u^2 v + p_{11}v + 2p_{12}u) \\ \frac{1}{2}(\rho v^2 u + p_{22}u + 2p_{12}v) \end{pmatrix}, \begin{pmatrix} \rho v \\ \rho u v + p_{12} \\ \rho v^2 + p_{22} \\ \frac{1}{2}(\rho u^2 v + p_{11}v + 2p_{12}u) \\ \frac{1}{2}(\rho v^2 u + p_{22}u + 2p_{12}v) \\ \frac{1}{2}(\rho v^3 + 3p_{22}v) \end{pmatrix} \right\}
$$

2.3.1 Splitting Dimensionnel

On souhaite maintenant extraire la structure d'onde d'un tel système. Le but est
de montrer que ce système est hyperbolique et d'expliciter les valeurs propres, vecteurs
propres et les formes linéaires propres de la matrice Jacobienne généralisée du système.
Les Jacobiennes dans chacune des directions spatiales s'écrivent :

$$
A = \frac{\partial \mathbf{F}}{\partial \mathbf{W}} =
\begin{pmatrix}
0 & 1 & 0 & 0 & 0 & 0 \\
0 & 0 & 0 & 2 & 0 & 0 \\
0 & 0 & 0 & 0 & 2 & 0 \\
\frac{u}{2}(u^2 - 3c_{11}^2) & \frac{3}{2}(c_{11}^2 - u^2) & 0 & 3u & 0 & 0 \\
\frac{u^2 - c_{11}^2}{2}v - c_{12}^2u & c_{12}^2 - uv & \frac{c_{11}^2 - u^2}{2} & v & 2u & 0 \\
\frac{v^2 - c_{22}^2}{2}u - c_{12}^2v & \frac{c_{22}^2 - v^2}{2} & c_{12}^2 - uv & 0 & 2v & u
\end{pmatrix}
\quad (19)
$$

$$
B = \frac{\partial \mathbf{G}}{\partial \mathbf{W}} =
\begin{pmatrix}
0 & 0 & 1 & 0 & 0 & 0 \\
0 & 0 & 0 & 0 & 2 & 0 \\
0 & 0 & 0 & 0 & 0 & 2 \\
\frac{u^2 - c_{11}^2}{2}v - c_{12}^2u & c_{12}^2 - uv & \frac{c_{11}^2 - u^2}{2} & v & 2u & 0 \\
\frac{v^2 - c_{22}^2}{2}u - c_{12}^2v & \frac{c_{22}^2 - v^2}{2} & c_{12}^2 - uv & 0 & 2v & u \\
\frac{v}{2}(v^2 - 3c_{22}^2) & 0 & \frac{3}{2}(c_{22}^2 - v^2) & 0 & 0 & 3v
\end{pmatrix}
\quad (20)
$$

Les quantités $c_{ij}^2$ s'apparentent ici à des vitesses de propagation de l'information et sont
en fait les composantes du tenseur $\Sigma$, que l'on notera également $\mathbb{C}$ :

$$
c_{ij}^2 = \frac{p_{ij}}{\rho} = \sigma_{ij}, \quad \forall i, j \in \{1, 2\}. \tag{21}
$$

On remarquera que par positivité des termes diagonaux de Σ les quantités c₁₁ et c₂₂ ont un sens mais pas c₁₂. D'ailleurs cette variable apparaîtra toujours au carré dans le reste du rapport : c²₁₂.
On cherche maintenant à expliciter les vecteurs propres $X_A = (x_1, \dots, x_6)^t$ associés aux valeurs propres $\lambda_A$ de la matrice $A$. Les équations sur les composantes de $X_A$ obtenues se découpulent assez bien. Voici ce que l'on obtient, ligne par ligne :

$$
\left\{
\begin{aligned}
x_4 &= \frac{\lambda}{2}x_2 = \frac{\lambda^2}{2}x_1 \\
x_5 &= \frac{\lambda}{2}x_3 \\
(\lambda - u) &((\lambda - u)^2 - 3c_{11}^2) x_1 = 0 \\
((\lambda - u)^2 - c_{11}^2) (vx_1 - x_3) + 2c_{12}^2(\lambda - u)x_1 = 0 \\
((\lambda - u)(c_{22}^2 - v^2) - 2c_{12}^2v) \frac{x_2}{2} + ((\lambda - u)v + c_{12}^2) x_3 = (\lambda - u)x_6
\end{aligned}
\right.
$$

On en déduit assez rapidement l'existence de 6 valeurs propres

$$
\lambda_A = u \pm \sqrt{3}c_{11}, u \pm c_{11}, u, \qquad (22)
$$

la dernière étant double, que l'on classe comme d'habitude par ordre croissant. En ré-
solvant complètement le système précédent, on trouve la matrice des vecteurs propres à
droite pour A :

$$
R_A =
\begin{pmatrix}
1 & 0 & 1 & 0 & 1 \\
u - \sqrt{3}c_{11} & 0 & u & 0 & u + \sqrt{3}c_{11} \\
v - \sqrt{3}\frac{c_{12}^2}{c_{11}} & 1 & v & 0 & v + \sqrt{3}\frac{c_{12}^2}{c_{11}} \\
\frac{1}{2}\left(u - \sqrt{3}c_{11}\right)^2 & 0 & \frac{u^2}{2} & 0 & \frac{1}{2}\left(u + \sqrt{3}c_{11}\right)^2 \\
\frac{1}{2}\left(u - \sqrt{3}c_{11}\right)\left(v - \sqrt{3}\frac{c_{12}^2}{c_{11}}\right) & \frac{u - c_{11}}{2} & \frac{uv}{2} & 0 & \frac{1}{2}\left(u + \sqrt{3}c_{11}\right)\left(v + \sqrt{3}\frac{c_{12}^2}{c_{11}}\right) \\
\frac{1}{2}\left(v - \sqrt{3}\frac{c_{12}^2}{c_{11}}\right)^2 + \frac{|\mathbb{C}|}{2c_{11}^2} & v - \frac{c_{12}^2}{c_{11}} & \frac{v^2}{2} & 1 & \frac{1}{2}\left(v + \sqrt{3}\frac{c_{12}^2}{c_{11}}\right)^2 + \frac{|\mathbb{C}|}{2c_{11}^2}
\end{pmatrix}
\tag{23}
$$

et la matrice de vecteurs propres à gauche, également appelés "formes linéaires propres",
qui est donnée en annexe pour des questions de mise-en-page. Annexe C, équations (101),
(102), (103). Dans une stratégie de splitting dimensionnel en deux dimensions, on aura
également besoin de la structure d'onde pour la matrice Jacobienne B, équation (20), que
l'on donne ici. On trouvera en annexe C la matrice des formes linéaires propres associées
(102), ainsi que le code Maxima qui a permis de vérifier tous ces calculs.

$$
R_B =
\begin{pmatrix}
1 & 0 & 1 & 0 & 1 \\
u - \sqrt{3} \frac{c_{12}^2}{c_{22}} & 1 & u & 0 & u + \sqrt{3} \frac{c_{12}^2}{c_{22}} \\
v - \sqrt{3} c_{22} & 0 & v & 0 & v + \sqrt{3} c_{22} \\
\frac{1}{2} \left( u - \sqrt{3} \frac{c_{12}^2}{c_{22}} \right)^2 + \frac{|\mathbb{C}|}{2c_{22}^2} & u - \frac{c_{12}^2}{c_{22}} & \frac{u^2}{2} & 1 & u + \frac{c_{12}^2}{c_{22}} & \frac{1}{2} \left( u + \sqrt{3} \frac{c_{12}^2}{c_{22}} \right)^2 + \frac{|\mathbb{C}|}{2c_{22}^2} \\
\frac{1}{2} \left( v - \sqrt{3} c_{22} \right) \left( u - \sqrt{3} \frac{c_{12}^2}{c_{22}} \right) & \frac{v - c_{22}^2}{2} & \frac{uv}{2} & 0 & \frac{v + c_{22}^2}{2} & \frac{1}{2} \left( v + \sqrt{3} c_{22} \right) \left( u + \sqrt{3} \frac{c_{12}^2}{c_{22}} \right) \\
\frac{1}{2} \left( v - \sqrt{3} c_{22} \right)^2 & 0 & \frac{v^2}{2} & 0 & 0 & \frac{1}{2} \left( v + \sqrt{3} c_{22} \right)^2
\end{pmatrix}
\tag{24}
$$
### 2.3.2 Hyperbolicité dans le cas général

Dans un cadre totalement non-structuré, il nous faut connaître la structure d'onde dans une direction unitaire $\vec{n} = (n_x, n_y)$ quelconque. Avec Maxima, voire annexe C, on montre que la matrice $J = A.n_x + B.n_y$ est diagonalisable dans toutes les directions, de valeur propres

$$ \lambda_J = \vec{u} \cdot \vec{n} \pm \sqrt{3}c_{nn}, \quad \vec{u} \cdot \vec{n} \pm c_{nn}, \quad \vec{u} \cdot \vec{n}, \quad (25) $$

où l'on a utilisé la notation

$$ c_{nn} = \sqrt{\mathbb{C}(\vec{n}) \cdot \vec{n}} = \sqrt{c_{11}^2 n_x^2 + 2c_{12}^2 n_x n_y + c_{22}^2 n_y^2}, \quad (26) $$

et la dernière valeur propre est toujours double. Dans ce contexte, on peut écrire la matrice des vecteurs propres à droite dans une forme peu orthodoxe :

$$ R = (\mathbf{r}_1, \mathbf{r}_2, \mathbf{r}_3, \mathbf{r}_4, \mathbf{r}_5, \mathbf{r}_6) \quad (27) $$

comme

$$
\begin{align*}
\mathbf{r}_{6^+/1^-} &= \begin{pmatrix} 1 & , & \vec{u} \pm \frac{\sqrt{3}\mathbb{C}(\vec{n})}{\sqrt{\mathbb{C}(\vec{n}) \cdot \vec{n}}} & , & \frac{1}{2} \begin{pmatrix} \vec{u} \pm \frac{\sqrt{3}\mathbb{C}(\vec{n})}{\sqrt{\mathbb{C}(\vec{n}) \cdot \vec{n}}} \end{pmatrix}^{\otimes 2} + \frac{c_{nn}^2 \mathbb{C} - \mathbb{C}(\vec{n}) \otimes \mathbb{C}(\vec{n})}{2c_{nn}^2} \end{pmatrix}^t \\
\mathbf{r}_{5^+/2^-} &= \begin{pmatrix} 0 & , & \vec{t} & , & \begin{pmatrix} \vec{u} \pm \frac{\mathbb{C}(\vec{n})}{\sqrt{\mathbb{C}(\vec{n}) \cdot \vec{n}}} \end{pmatrix} \odot \vec{t} \end{pmatrix}^t \\
\mathbf{r}_3 &= \begin{pmatrix} 1 & , & \vec{u} & , & \frac{\vec{u} \otimes \vec{u}}{2} \end{pmatrix}^t \\
\mathbf{r}_4 &= \begin{pmatrix} 0 & , & \vec{0} & , & \vec{t} \otimes \vec{t} \end{pmatrix}^t
\end{align*}
$$

où $\vec{t}$ est le vecteur tangent unitaire tel que $(\vec{n}, \vec{t})$ forme une base orthonormée :

$$ \vec{t} = (-n_y, n_x)^t. $$

De même, on écrit avec la même convention les formes linéaires propres, voir équation (103) dans l'annexe C. Pour vérifier, si l'on note les vecteurs à droite et à gauche comme :

$$ \mathbf{r}_k = (r_k, \vec{r}_k, \mathbb{R}_k)^t \quad \text{et} \quad \mathbf{l}_k = (l_k, \vec{l}_k, \mathbb{L}_k), \quad (28) $$

on a la relation d'orthonormalité du repère propre :

$$ < \mathbf{l}_i, \mathbf{r}_j > = l_i r_j + \vec{l}_i \cdot \vec{r}_j + \mathbb{L}_i :: \mathbb{R}_j = \delta_{ij}. $$

Pour cela, on aura besoin des formules simples suivantes pour des vecteurs 2D $\vec{w}, \vec{x}, \vec{y}$ et $\vec{z}$ quelconques :

$$ (\vec{x} \otimes \vec{x}) :: (\vec{y} \otimes \vec{y}) = (\vec{x} \cdot \vec{y})^2, \quad (\vec{x} \otimes \vec{x}) :: (\vec{y} \odot \vec{z}) = (\vec{x} \cdot \vec{y})(\vec{x} \cdot \vec{z}), \quad (29) $$

$$ (\vec{w} \odot \vec{x}) :: (\vec{y} \odot \vec{z}) = \frac{1}{2} ((\vec{w} \cdot \vec{y})(\vec{x} \cdot \vec{z}) + (\vec{w} \cdot \vec{z})(\vec{x} \cdot \vec{y})), \quad (30) $$

$$ \vec{x}^{\perp} \cdot \vec{y} = -\vec{x} \cdot \vec{y}^{\perp} \quad \text{et} \quad \vec{x}^{\perp} \cdot \vec{y}^{\perp} = \vec{x} \cdot \vec{y}, \quad (31) $$

et des relations pour le tenseur $\mathbb{C}$ :

$$ \mathbb{C} :: (\vec{n} \otimes \vec{n}) = c_{nn}^2, \quad (\mathbb{C}(\vec{n})^{\perp} \odot \vec{n}) :: \mathbb{C} = 0, \quad (\mathbb{C}(\vec{n})^{\perp} \otimes \mathbb{C}(\vec{n})^{\perp}) :: \mathbb{C} = |\mathbb{C}|c_{nn}^2. $$
### 2.3.3 Structure d'onde en variables primitives

On souhaite maintenant réaliser la même étude dans les variables

$$ \mathbf{V} = (\rho, \vec{\mathbf{u}}, \mathbb{P})^t. \qquad (32) $$

Si l'on note $\mathcal{D}_t = \partial_t + (\vec{\mathbf{u}} \cdot \nabla)$, la dérivée particulaire, les équations (18) peuvent se réécrire comme

$$ \left\{ \begin{array}{l} \mathcal{D}_t \rho + \rho \operatorname{div}(\vec{\mathbf{u}}) = 0 \\[1em] \mathcal{D}_t \vec{\mathbf{u}} + \frac{1}{\rho} \operatorname{div}(\mathbb{P}) = 0 \\[1em] \mathcal{D}_t \mathbb{P} + \mathbb{P} \operatorname{div}(\vec{\mathbf{u}}) + \sum_k (\mathbb{P}_{ik} \partial_{x_k} u_j + \mathbb{P}_{jk} \partial_{x_k} u_i) = 0 \end{array} \right. \qquad (33) $$

Les matrices Jacobiennes s'écrivent

$$ A = \begin{pmatrix} u & \rho & 0 & 0 & 0 & 0 \\ 0 & u & 0 & \frac{1}{\rho} & 0 & 0 \\ 0 & 0 & u & 0 & \frac{1}{\rho} & 0 \\ 0 & 3p_{11} & 0 & u & 0 & 0 \\ 0 & 2p_{12} & p_{11} & 0 & u & 0 \\ 0 & p_{22} & 2p_{12} & 0 & 0 & u \end{pmatrix} \qquad (34) $$

$$ B = \begin{pmatrix} v & 0 & \rho & 0 & 0 & 0 \\ 0 & v & 0 & 0 & \frac{1}{\rho} & 0 \\ 0 & 0 & v & 0 & 0 & \frac{1}{\rho} \\ 0 & 2p_{12} & p_{11} & v & 0 & 0 \\ 0 & p_{22} & 2p_{12} & 0 & v & 0 \\ 0 & 0 & 3p_{22} & 0 & 0 & v \end{pmatrix} \qquad (35) $$

$$ J = A \cdot nx + B \cdot ny = \begin{pmatrix} \vec{\mathbf{u}} \cdot \vec{\mathbf{n}} & \rho nx & \rho ny & 0 & 0 & 0 \\ 0 & \vec{\mathbf{u}} \cdot \vec{\mathbf{n}} & 0 & \frac{nx}{\rho} & \frac{ny}{\rho} & 0 \\ 0 & 0 & \vec{\mathbf{u}} \cdot \vec{\mathbf{n}} & 0 & \frac{nx}{\rho} & \frac{ny}{\rho} \\ 0 & 3p_{11}nx + 2p_{12}ny & p_{11}ny & \vec{\mathbf{u}} \cdot \vec{\mathbf{n}} & 0 & 0 \\ 0 & 2p_{12}nx + p_{22}ny & p_{11}nx + 2p_{12}ny & 0 & \vec{\mathbf{u}} \cdot \vec{\mathbf{n}} & 0 \\ 0 & p_{22}nx & 2p_{12}nx + 3p_{22}ny & 0 & 0 & \vec{\mathbf{u}} \cdot \vec{\mathbf{n}} \end{pmatrix} \qquad (36) $$

Grâce à Maxima et beaucoup de calcul, on retrouve la structure d'onde trouvée pré-
cédemment dans le cadre des variables conservatives, avec

$$ \lambda_{6/1}^{\pm} = \vec{\mathbf{u}} \cdot \vec{\mathbf{n}} \pm \sqrt{3} c_{nn}, \quad \lambda_{5/2}^{\pm} = \vec{\mathbf{u}} \cdot \vec{\mathbf{n}} \pm c_{nn}, \quad \lambda_3 = \lambda_4 = \vec{\mathbf{u}} \cdot \vec{\mathbf{n}}. \qquad (37) $$

$$ \begin{align} \vec{r}_{6/1}^{\pm} &= \left( \begin{array}{c} 1 \\ \pm \frac{\sqrt{3}\mathbb{C}(\vec{n})}{\rho \mathbb{C}(\vec{n}) \cdot \vec{n}} \end{array} , \begin{array}{c} \frac{3\mathbb{C}(\vec{n}) \otimes \mathbb{C}(\vec{n}) + |\mathbb{C}| \vec{t} \otimes \vec{t}}{\mathbb{C}(\vec{n}) \cdot \vec{n}} \\ \pm \frac{2\rho \mathbb{C}(\vec{n}) \odot \vec{t}}{c_{nn}} \end{array} \right)^t, \nonumber \\ \vec{r}_{5/2}^{\pm} &= \left( \begin{array}{c} 0 \\ \vec{t} \end{array} , \begin{array}{c} \rho \mathbb{C}(\vec{n}) \cdot \vec{n} \end{array} \right)^t, \nonumber \\ \vec{r}_3 &= \left( \begin{array}{c} 1 \\ \vec{0} \end{array} , \begin{array}{c} \mathbb{O} \end{array} \right)^t, \nonumber \\ \vec{r}_4 &= \left( \begin{array}{c} 0 \\ \vec{0} \end{array} , \begin{array}{c} \vec{t} \otimes \vec{t} \end{array} \right)^t \tag{38} \end{align} $$
et les formes linéaires propres

$$
\begin{equation}
\begin{aligned}
\vec{l}_{6/1}^{\pm} &= \begin{pmatrix} 0 & , & \pm \frac{\rho}{2\sqrt{3}c_{nn}} \mathbf{n} & , & \frac{\mathbf{n} \otimes \mathbf{n}}{6c_{nn}^2} \\ & & \mathbb{C}(\mathbf{n})^{\perp} & , & \end{pmatrix}, \\
\vec{l}_{5/2}^{\pm} &= \begin{pmatrix} 0 & , & \frac{\mathbb{C}(\mathbf{n})^{\perp}}{2\mathbb{C}(\mathbf{n}) \cdot \mathbf{n}} & , & \pm \frac{\mathbb{C}(\mathbf{n})^{\perp} \odot \mathbf{n}}{2 * \rho c_{nn}^3} \\ & & \mathbb{C}(\mathbf{n})^{\perp} & , & \end{pmatrix}, \\
\vec{l}_{3} &= \begin{pmatrix} 1 & , & \mathbf{0} & , & -\frac{\mathbf{n} \otimes \mathbf{n}}{3c_{nn}^2} \\ & & \mathbb{C}(\mathbf{n})^{\perp} & , & \end{pmatrix}, \\
\vec{l}_{4} &= \begin{pmatrix} 0 & , & \mathbf{0} & , & \frac{3\mathbb{C}(\mathbf{n})^{\perp} \otimes \mathbb{C}(\mathbf{n})^{\perp} - |\mathbb{C}|\mathbf{n} \otimes \mathbf{n}}{3c_{nn}^4} \\ & & \mathbb{C}(\mathbf{n})^{\perp} & , & \end{pmatrix}
\end{aligned}
\tag{39}
\end{equation}
$$

2.4 Résolution du problème de Riemann dans une direction quel-
conque

**2.4.1 Transformation invariante**

On rappelle d'abord succinctement le problème, toujours en 2 dimensions, mais les résultats suivants se généralisent. Il s'agit du système de loi de conservation

$$
\frac{\partial \mathbf{W}}{\partial t} + \operatorname{div} (\vec{\mathcal{F}}(\mathbf{W})) = 0, \quad (40)
$$

où, en notation tensorielle,

$$
\mathbf{W} = (\rho, \rho\vec{u}, \rho\mathbb{E})^t, \quad \vec{\mathcal{F}}(\mathbf{W}) = (F(\mathbf{W}), G(\mathbf{W})) = (\rho\vec{u}, \rho\vec{u} \otimes \vec{u} + \mathbb{P}, \rho\mathbb{H} \odot \vec{u})^t. \quad (41)
$$

On se réfèrera au paragraphe 2.3 pour plus de détails.

Pour tout vecteur unitaire **n** = (cos θ, sin θ), on définit la matrice de rotation

$$
\mathbf{R} = \begin{pmatrix} \cos \theta & \sin \theta \\ -\sin \theta & \cos \theta \end{pmatrix} \tag{42}
$$

et les opérateurs algébriques linéaires inversibles suivants, projetant respectivement les
variables et les flux dans la direction **n** :

$$
\mathcal{T} : \left\{
\begin{array}{l@{\quad}c@{\quad}l@{\quad}c@{\quad}l}
\mathbb{R}_{\mathbf{W}}^{6} & \to & \mathbb{R}_{\mathbf{W}}^{6} & & \\
(\rho, \vec{u}, \mathbf{T}) & \mapsto & (\rho, \mathbf{R}\vec{u}, \mathbf{R}\mathbf{T}\mathbf{R}^{-1}) & &
\end{array}
\right.
\tag{43}
$$

et

$$
\mathcal{R} : \left\{
\begin{array}{lcl}
(\mathbb{R}^n)^2 & \to & (\mathbb{R}^n)^2 \\[1em]
\vec{\mathcal{F}} = (\mathbf{F}, \mathbf{G}) & \mapsto & (\vec{\mathcal{F}} \cdot \mathbf{n}, \vec{\mathcal{F}} \cdot \vec{t}) = (\cos \theta \mathbf{F} + \sin \theta \mathbf{G}, -\sin \theta \mathbf{F} + \cos \theta \mathbf{G})
\end{array}
\right.
\tag{44}
$$

Si l'on note $(d\vec{n}, d\vec{t})$ les formes linéaires canoniques dans le repère lié à $(\vec{n}, \vec{t})$, on a
évidemment

$$
(d \vec{n}, d \vec{t})^t = \mathbf{R} \cdot (dx, dy)^t
$$

et il vient que pour tout flux bi-dimensionnel $\vec{\mathcal{F}} = (\mathbf{F}, \mathbf{G})$

$$
\nabla_{(x,y)} \cdot (\mathbf{R}^{-1}(\vec{\mathcal{F}})) = \partial_x (\cos \theta \mathbf{F} + \sin \theta \mathbf{G}) + \partial_y (-\sin \theta \mathbf{F} + \cos \theta \mathbf{G}) = \nabla_{(\vec{n}, \vec{t})} \cdot \vec{\mathcal{F}} \quad (45)
$$
Or, dans le cas du système avec pression anisotrope, on peut vérifier la relation suivante, à la main ou en utilisant Maxima (cf. feuille de calcul en annexe) :

$$ \mathcal{R}(\vec{\mathcal{F}}) = \mathcal{T}^{-1} [\vec{\mathcal{F}}(\mathcal{T}\mathbf{W})]. \qquad (46) $$

Ainsi, si (40) est vérifié dans un repère cartésien arbitraire $(x, y)$, pour toute direction unitaire $\vec{n}$, on a :

$$ \begin{align*} \partial_t \mathbf{W} + \nabla_{(x,y)} \cdot (\vec{\mathcal{F}}(\mathbf{W})) &= 0 \\ \Leftrightarrow \partial_t \mathbf{W} + \nabla_{(x,y)} \cdot (\mathcal{R}^{-1} [\mathcal{T}^{-1} \vec{\mathcal{F}}(\mathcal{T}\mathbf{W})]) &= 0 \\ \Leftrightarrow \partial_t \mathcal{T}\mathbf{W} + \nabla_{(\vec{n},\vec{t})} \cdot \vec{\mathcal{F}}(\mathcal{T}\mathbf{W}) &= 0. \end{align*} $$

Une fois ceci démontré, le problème de Riemann à travers une surface $\Sigma$ de normale $\vec{n}$, s'écrit

$$ \partial_t \mathcal{T}\mathbf{W} + \partial_{\vec{n}} \mathcal{F}(\mathcal{T}\mathbf{W}) = 0, $$

où les dérivées le long de la direction $\vec{t}$ on été omises par symétrie du problème. Par multiplication à gauche par $\mathcal{T}^{-1}$ et en utilisant (46), le problème considéré est bien équivalent à :

$$ \frac{\partial \mathbf{W}}{\partial t} + \frac{\partial \vec{\mathcal{F}}(\mathbf{W}) \cdot \vec{n}}{\partial \vec{n}} = 0. \qquad (47) $$

On est ainsi ramené à un problème de Riemann 1D dans le repère $(\vec{n}, t)$, où le flux est bien $\vec{\mathcal{F}} \cdot \vec{n}$, avec la structure d'onde préliminairement étudiée lors des paragraphes 2.3.2 et 2.3.3. Dans la suite, par abus de notation, la variable $x$ représentera la coordonnée spatiale dans la direction normale $\vec{n}$.

On se ramène donc à un problème de Riemann monodimensionnel pour le système (47) avec un état gauche $\mathbf{W}_L$ et un état droit $\mathbf{W}_R$. Vu la dimension du problème, dans un cadre très général, la solution est composée de 6 ondes séparant 7 états constants, ce qui signifie $5 \times 6 = 30$ inconnues. En réalité la solution est bien plus simple, voir figure 2, et on peut grandement réduire le nombre d'inconnues, ce que nous allons faire durant les prochains paragraphes.

### 2.4.2 Nature de chacune des ondes

**Ondes 3 et 4** : Les ondes 3 et 4 sont toutes les deux associées à la valeur propre

$$ \lambda^{\otimes} = \vec{u} \cdot \vec{n}. $$

On a alors :

$$ \nabla_{\mathbf{W}} \lambda^{\otimes} = \left( -\frac{\vec{u} \cdot \vec{n}}{\rho}, \frac{\vec{n}}{\rho}, \mathbf{0} \right)^t. $$

Il vient immédiatement :

$$ \overline{\nabla_{\mathbf{W}} \lambda^{\otimes} \cdot \vec{r}_3} = \overline{\nabla_{\mathbf{W}} \lambda^{\otimes} \cdot \vec{r}_4} = 0. \qquad (48) $$

Les champs 3 et 4 sont donc **Linéairement Dégénérés**.
**Ondes 2 et 5** : Les ondes 2 et 5 sont associées aux valeurs propres

$$ \lambda_{\pm}^{ld} = \vec{u} \cdot \vec{n} \pm c_{nn}. $$

On peut alors calculer

$$ \nabla_w \lambda_{\pm}^{ld} = \frac{1}{\rho c_{nn}} \left( -\vec{u} \cdot \vec{n} c_{nn} \pm \left( \vec{u} \cdot \vec{n}^2 - \mathbb{E}(\vec{n}) \cdot \vec{n} \right) , (c_{nn} \mp \vec{u} \cdot \vec{n}) \vec{n} , \pm \vec{n} \otimes \vec{n} \right)^t, $$

pour voir immédiatement que

$$ \overline{\nabla_w \lambda_{-}^{ld} \cdot \vec{r}_2} = \overline{\nabla_w \lambda_{+}^{ld} \cdot \vec{r}_5} = 0, \quad (49) $$

et les champs 2 et 5 sont également **Linéairement Dégénérés**.

**Ondes 1 et 6** : Les ondes 1 et 6 sont associées aux valeurs propres

$$ \lambda_{\pm} = \vec{u} \pm \sqrt{3} c_{nn}. $$

Par analogie avec le calcul précédent,

$$ \nabla_w \lambda_{\pm} = \frac{1}{\rho c_{nn}} \left( -\vec{u} \cdot \vec{n} c_{nn} \pm \sqrt{3} \left( \vec{u} \cdot \vec{n}^2 - \mathbb{E}(\vec{n}) \cdot \vec{n} \right) , (c_{nn} \mp \sqrt{3} \vec{u} \cdot \vec{n}) \vec{n} , \pm \sqrt{3} \vec{n} \otimes \vec{n} \right)^t. $$

Si $\vec{r}_{\pm}$ désigne respectivement $\vec{r}_6$ et $\vec{r}_1$, alors on peut écrire :

$$ \overline{\nabla_w \lambda_{\pm} \cdot \vec{r}_{\pm}} = \pm 2 \sqrt{3} \frac{c_{nn}}{\rho}, \quad (50) $$

ce qui justifie le caractère **Vraiment Non-Linéaire** des champs 1 et 6.

### 2.4.3 Relations de saut pour le Gaussien anisotrope

<figure><img src="image_2.png" /><figcaption>FIGURE 1 – Représentation du choc droit se déplaçant à la vitesse $\sigma \vec{n}$ en deux dimensions d'espace. Définition de la cellule de contrôle $C$.</figcaption></figure>
Soit $\Sigma$ une discontinuité séparant des états constants $W_1$ et $W_2$, de normale unitaire $\vec{n}$ et se déplaçant à la vitesse $\vec{\sigma} = \sigma\vec{n}$. Voir Figure 1. Dans la suite, il est pratique de considérer la vitesse du fluide dans le référentiel du choc

$$ \vec{v} = \vec{u} - \vec{\sigma}. \qquad (51) $$

Par la relation de Rankine-Hugoniot pour l'équation (47), on obtient :

— Masse :

$$ [[\rho \vec{u} \cdot \vec{n}] ] = \sigma [[\rho] ] \Leftrightarrow [[\rho(\vec{u} \cdot \vec{n}) - \sigma)]] = [[j]] = 0 \qquad (52) $$

— Quantité de Mouvement :

$$ [[(\rho \vec{u} \otimes \vec{u} + P) \{\vec{n}\}]] = \sigma [[\rho \vec{u}]] \Leftrightarrow j[[\vec{u}]] + [[P(\vec{n})]] = 0 \qquad (53) $$

— Énergie :

$$ [[\rho (\vec{u} \odot H(\vec{u}, C)) \{\vec{n}\}]] = \sigma [[\rho E]] \qquad (54) $$

où

$$ 2\rho H(\vec{u}, C) = \rho \vec{u} \otimes \vec{u} + 3P \Rightarrow 2\rho \vec{u} \odot H(\vec{u}, C) = \rho \vec{u}^{\otimes 3} + 3\vec{u} \odot P. $$

De plus

$$ \begin{aligned} (3\vec{u} \odot P) \{\vec{n}\} &= \sum_k (\vec{u} \odot P)_{ijk} n_k \\ &= \sum_k (P_{ij}u_k + P_{ik}u_j + P_{jk}u_i) n_k \\ &= P\vec{u} \cdot \vec{n} + 2\vec{u} \odot P(\vec{n}), \end{aligned} $$

donc en prenant deux fois (54), on obtient :

$$ j[[\vec{u} \otimes \vec{u} + P/\rho]] + 2[[\vec{u} \odot P(\vec{n})]] = 0 \qquad (55) $$

Il est important de noter que l'on peut utiliser $\vec{v}$ à la place de $\vec{u}$ dans les relations de saut (53) et (54). La justification pour la première relation est immédiate, puisque $[[\sigma]] = 0$. C'est en revanche beaucoup plus délicat pour la seconde. On peut d'une part s'en convaincre par le calcul, mais on peut également appliquer la conservation de l'énergie totale sur la cellule de contrôle $C$ définie sur la Figure 1. Comme cette cellule est liée au référentiel de la discontinuité, on a par intégration de l'équation tensorielle d'énergie :

$$ [[\rho (\vec{v} \odot H(\vec{v}, C)) \{\vec{n}\}]] = 0. $$

On raisonne comme précédemment en utilisant en particulier la relation

$$ (3\vec{v} \odot P) \{\vec{n}\} = \vec{v} \cdot \vec{n} P + 2\vec{v} \odot P(\vec{n}) $$

et le fait que $[[\vec{v} \cdot \vec{n}] = [[j]] = 0$, pour obtenir

$$ j[[\vec{v} \otimes \vec{v} + P/\rho]] + 2[[\vec{v} \odot P(\vec{n})]] = 0. $$

En projetant cette équation tensorielle sur sa composante $(\vec{n}, \vec{n})$, on remarque alors que l'enthalpie totale normale se conserve dans le référentiel du choc, c'est-à-dire

$$ [[H(\vec{v}, C) \{\vec{n}\} \cdot \vec{n}] ] = 0. \qquad (56) $$
Attention, cette relation est fausse quant à l'enthalpie totale du fluide dans le référentiel absolue. En effet $\mathbb{H}(\vec{u}, \mathbb{C})$ subit un travail de pression lors du passage du choc, ce que l'on peut voir en regardant la composante $(\vec{n}, \vec{n})$ de l'équation tensorielle (54) :

$$j[\mathbb{H}(\vec{u}, \mathbb{C})\{\vec{n}\}.\vec{n}] + 2\sigma[\mathbb{P}_{nn}] = 0$$

### 2.4.4 Relation de choc droit pour le Gaussien anisotrope :

On se place dans le cas où la discontinuité est un choc, c'est à dire que $j = \rho \vec{v} \cdot \vec{n} \neq 0$. On rappelle les trois relations scalaires de saut :

Masse :

$$\rho_1 \vec{v} \cdot \vec{n}_1 = \rho_2 \vec{v} \cdot \vec{n}_2, \qquad (57)$$

QdM :

$$\rho_1(\vec{v} \cdot \vec{n}_1)^2 + \mathbb{P}_{nn1} = \rho_2(\vec{v} \cdot \vec{n}_2)^2 + \mathbb{P}_{nn2}, \qquad (58)$$

Énergie :

$$(\vec{v} \cdot \vec{n}_1)^2 + 3c_{nn1}^2 = (\vec{v} \cdot \vec{n}_2)^2 + 3c_{nn2}^2, \qquad (59)$$

ainsi que les relations suivantes :

$$\mathbb{P} = \rho \mathbb{C} \quad \text{et} \quad c_{nn} = \sqrt{\mathbb{C}(\vec{n}) \cdot \vec{n}}.$$

On définit enfin le nombre de Mach normal comme

$$M_{nn} = \frac{\vec{v} \cdot \vec{n}}{\sqrt{3c_{nn}}}. \qquad (60)$$

En divisant les équations (58) et (59) par $3c_{nn1}^2$ et en utilisant extensivement la conservation de la masse (57) pour éliminer les occurences de $\vec{v} \cdot \vec{n}_2$, on obtient les relations

$$\left(\frac{c_{nn2}}{c_{nn1}}\right)^2 = \frac{\rho_1}{\rho_2} (3M_{nn1}^2 + 1) - 3 \left(\frac{\rho_1}{\rho_2}\right)^2 M_{nn1}^2 \qquad (61)$$

et

$$\left(\frac{c_{nn2}}{c_{nn1}}\right)^2 = M_{nn1}^2 + 1 - \left(\frac{\rho_1}{\rho_2}\right)^2 M_{nn1}^2. \qquad (62)$$

La combinaison de ces deux relations fournie cette équation de degré 2 en $\rho_1/\rho_2$

$$2M_{nn1}^2 \left(\frac{\rho_1}{\rho_2}\right)^2 - (3M_{nn1}^2 + 1)\frac{\rho_1}{\rho_2} + M_{nn1}^2 + 1 = 0,$$

dont le discriminant est $\Delta = (M_{nn1}^2 - 1)^2$ et les deux racines sont 1 (peu intéressant) et

$$\frac{\rho_2}{\rho_1} = \frac{2M_{nn1}^2}{M_{nn1}^2 + 1}. \qquad (63)$$

Par (61) ou (62), il vient

$$\frac{c_{nn2}^2}{c_{nn1}^2} = \frac{(M_{nn1}^2 + 1)(3M_{nn1}^2 - 1)}{4M_{nn1}^2} \qquad (64)$$
et par suite

$$ \frac{\mathbb{P}_{nn2}}{\mathbb{P}_{nn1}} = \frac{\rho_2 c_{nn2}^2}{\rho_1 c_{nn1}^2} = \frac{1}{2} (3M_{nn1}^2 - 1) \quad (65) $$

puis par l'équation de conservation de la masse

$$ \frac{\vec{v} \cdot \vec{n}_2}{\vec{v} \cdot \vec{n}_1} = \frac{\rho_1}{\rho_2} = \frac{M_{nn1}^2 + 1}{2M_{nn1}^2}. \quad (66) $$

On remarquera qu'on obtient exactement les mêmes relations de choc que pour Euler, lorsque $\gamma$ vaut 3 et que la vitesse du son vaut $c = \sqrt{3}c_{nn}$. On se servira de cette propriété par la suite.

Enfin, les relations de saut supplémentaires :

$$ \left\{ \begin{array}{l} j\left[\left\|\vec{u} \cdot \vec{t}\right\|\right] + \left\|\mathbb{P}_{nt}\right\| = 0 \\ j\left[\left\|\left(\vec{u} \cdot \vec{n}\right)(\vec{u} \cdot \vec{t})\right\| + \mathbb{P}_{nt}/\rho\right] + \left\|\left[(\vec{u} \cdot \vec{n})\mathbb{P}_{nt} + (\vec{u} \cdot \vec{t})\mathbb{P}_{nn}\right]\right\| = 0 \\ j\left[\left\|\left(\vec{u} \cdot \vec{t}\right)^2\right\| + \mathbb{P}_{tt}/\rho\right] + 2\left\|\left[(\vec{u} \cdot \vec{t})\mathbb{P}_{nt}\right]\right\| = 0 \end{array} \right. \quad (67) $$

permettent de trouver les inconnues restantes.

### 2.4.5 Cas de la détente

Comme dans [1], on montre que $s_3 = \frac{\mathbb{P}_{nn}}{\rho^3}$ est une entropie pour le système (47). Pour cela, on considère la conservation de la masse

$$ \partial_t \rho + (\vec{u} \cdot \vec{n}) \partial_x \rho + \rho \partial_x \vec{u} \cdot \vec{n} \quad (68) $$

et la combinaison suivante des lignes 4, 5 et 6 :

$$ n_x^2 L_4 + 2n_x n_y L_5 + n_y^2 L_6 \Rightarrow \partial_t \mathbb{P}_{nn} + \vec{u} \cdot \vec{n} \partial_x \mathbb{P}_{nn} + 3 \mathbb{P}_{nn} \partial_x \vec{u} \cdot \vec{n} = 0 \quad (69) $$

Si on divise cette dernière relation par $\rho^3$ et qu'on lui retranche $\frac{\mathbb{P}_{nn}}{\rho^4}$ fois la précédente, on obtient alors :

$$ \partial_t \left( \frac{\mathbb{P}_{nn}}{\rho^3} \right) + \vec{u} \cdot \vec{n} \partial_x \left( \frac{\mathbb{P}_{nn}}{\rho^3} \right) = 0. \quad (70) $$

On trouve ici la conservation d'une des entropies du système (47) : $s_3 = \frac{\mathbb{P}_{nn}}{\rho^3}$, qui se conserve à travers les détentes. TODO : réf vers Berthon, lien avec l'entropie micro et Th. H.

Ensuite, en projetant le système exprimé en variables primitives (33) sur chacune de ses formes linéaires propres (39), on en déduit les invariants de Riemann pour chacune des courbes caractéristiques. L'annexe D fournie plus de détails. Pour les deux champs vraiment non-linéaires, on a :

$$ \partial_{\epsilon_{\pm}} (\mathbb{P}_{nn}) \pm \sqrt{3} \rho c_{nn} \partial_{\epsilon_{\pm}} (\vec{u} \cdot \vec{n}) = 0. \quad (71) $$

Pour les champs 5 et 2, on a :

$$ \left( \partial_{\epsilon_{\pm}} (\mathbb{P}(\vec{n})) \pm \rho c_{nn} \partial_{\epsilon_{\pm}} (\vec{u}) \right) \cdot \mathbb{C}(\vec{n})^{\perp} = 0. \quad (72) $$
Enfin pour le 3-champs :

$$ \partial_{\epsilon_0} \rho - \frac{1}{3c_{nn}^2} \partial_{\epsilon_0} (\mathbb{P}_{nn}) = 0 \qquad (73) $$

et pour le 4-champs :

$$ (\partial_{\epsilon_0} \mathbb{P}) :: \{3\mathbb{C}(\vec{n})^{\perp} \otimes \mathbb{C}(\vec{n})^{\perp} - |\mathbb{C}|\vec{n} \otimes \vec{n}\} = 0. \qquad (74) $$

Il s'agit maintenant d'intégrer ces équations dans le cas d'une transformation régulière (onde simple). Commençons par la relation (73). Puisque $c_{nn}^2 = \mathbb{P}_{nn}/\rho$, on retrouve la conservation de l'entropie citée ci-dessus

$$ \partial_{\epsilon_0} s_3 = 0. $$

On remarquera d'ailleurs que la projection des équations en variables primitives sur la forme propre $\vec{l}_3$ correspond à la combinaison des équations proposée pour obtenir $s$. Passons à la relation différentielle (74). Dans son article de 2006 [1], C. Berthon propose une seconde entropie $s' = \frac{p_{11}p_{22} - p_{12}^2}{\rho^4} = \frac{|\mathbb{P}|}{\rho^4}$ pour le système dans la direction $\vec{n} = (1, 0)$. Ce n'est pas cette entropie que nous allons trouver ici, mais les trois entropies sont liées au sens d'une certaine algèbre multiplicative (TODO : Détailler !). Réécrivons (74) comme

$$ 3 (\mathbb{P}(\vec{n})^{\perp} \otimes \mathbb{P}(\vec{n})^{\perp} + |\mathbb{P}|) :: \partial_{\epsilon_0} \mathbb{P} - 4|\mathbb{P}|\partial_{\epsilon_0} (\mathbb{P}_{nn}) = 0. $$

Un rapide calcul montre que le premier terme se réduit simplement à $3 (\mathbb{P}_{nn}) \partial_{\epsilon_0} |\mathbb{P}|$, ce qui nous donne donc la forme de l'entropie conservée le longs des 4-caractéristiques :

$$ s_4 = \frac{|\mathbb{P}|^3}{(\mathbb{P}_{nn})^4}. \qquad (75) $$

Puisque les 3- et 4-caractéristiques sont confondues, $s_3$ et $s_4$ sont conservés au travers des mêmes régions. Dans ce cas, $s' = \sqrt[3]{\frac{s_4}{s_3^4}}$ est également une entropie du système.

Les ondes de détentes étant isentropiques, on en déduit en particulier que

$$ \frac{\mathbb{P}_{nn}}{\rho^3} = \left( \frac{c_{nn}}{\rho} \right)^2 = \bar{s}_3, $$

ce qui nous permet maintenant d'intégrer la relation (71). En effet,

$$ d\mathbb{P}_{nn} = 3\bar{s}_3\rho^2 d\rho $$

d'où

$$ d(\vec{u} \cdot \vec{n}) \pm \sqrt{3\bar{s}_3} d\rho = 0 $$

et les invariants de Riemann pour les $\mathscr{C}_{\pm}$

$$ J_{\pm} = \vec{u} \cdot \vec{n} \pm \sqrt{3} c_{nn} \qquad (76) $$

Il reste maintenant à intégrer les deux relations pour les ondes linéairement dégénérées associées aux valeurs propres 2 et 5. En écrivant $\mathbb{P}(\vec{n})$ et $\vec{u}$ comme

$$ \mathbb{P}(\vec{n}) = (\mathbb{P}_{nn}) \cdot \vec{n} + (\mathbb{P}_{nt}) \cdot \vec{t} \quad \text{et} \quad \vec{u} = (\vec{u} \cdot \vec{n}) \cdot \vec{n} + (\vec{u} \cdot \vec{t}) \cdot \vec{t}, $$
et en utilisant les relations (31), on peut réécrire les deux relations (72) comme

$$-\mathbb{P}_{nt}d(\mathbb{P}_{nn}) + \mathbb{P}_{nn}d(\mathbb{P}_{nt}) \pm \rho c_{nn}\mathbb{P}_{nn}d(\vec{\mathbf{u}}\cdot\vec{\mathbf{t}}) \mp \rho c_{nn}\mathbb{P}_{nt}d(\vec{\mathbf{u}}\cdot\vec{\mathbf{n}}) = 0.$$

Si on considère de plus qu'on se situe à travers une détente '+', alors l'invariant de Riemann $J_-$ est uniformément conservé à travers la détente et

$$d(\vec{\mathbf{u}}\cdot\vec{\mathbf{n}}) = \sqrt{3}d(c_{nn}), \qquad (77)$$

soit

$$\frac{\rho c_{nn}}{\mathbb{P}_{nt}} d(\vec{\mathbf{u}}\cdot\vec{\mathbf{t}}) = \sqrt{3} \frac{d(c_{nn})}{c_{nn}} \pm \frac{d(\mathbb{P}_{nn})}{\mathbb{P}_{nn}} \mp \frac{d(\mathbb{P}_{nt})}{\mathbb{P}_{nt}}. \qquad (78)$$

**Tour de Magie :** Dans le cas d'une onde simple, éventail de $\mathcal{C}_+$ (resp. $\mathcal{C}_-$), les deux invariants de Riemann exprimés dans (78), correspondant aux champs 2 et 5, se conservent le long des $\mathcal{C}_+^{ld}$. Or dans le cadre d'un problème de Riemann, les deux états externes $\mathbf{W}_R$ et $\mathbf{W}_L$ sont uniformes. Il en résulte que ces invariants de Riemann sont également uniformes dans l'éventail. Malgré le fait que les formes différentielles était initialement exprimées dans les directions particulières des $\mathcal{C}_+^{ld}$, elles sont maintenant valables dans toutes les directions, d'où la notation $d(\cdot)$. On peut donc sommer et retrancher ces deux relations et il vient :

$$\frac{d(\mathbb{P}_{nn})}{\mathbb{P}_{nn}} - \frac{d(\mathbb{P}_{nt})}{\mathbb{P}_{nt}} = 0 \quad \Rightarrow \quad \frac{\mathbb{P}_{nn}}{\mathbb{P}_{nt}} = c_{ste} \qquad (79)$$

et par suite, dans le cas d'un éventail de $\mathcal{C}_+$,

$$d(\vec{\mathbf{u}}\cdot\vec{\mathbf{t}}) - \sqrt{3} \frac{\mathbb{P}_{nt}}{\mathbb{P}_{nn}} d(c_{nn}) = 0 \quad \Rightarrow \quad \vec{\mathbf{u}}\cdot\vec{\mathbf{t}} - \sqrt{3} \frac{\mathbb{P}_{nt}}{\mathbb{P}_{nn}} c_{nn} = c_{ste}, \qquad (80)$$

et dans le cas d'un éventail de $\mathcal{C}_-$,

$$d(\vec{\mathbf{u}}\cdot\vec{\mathbf{t}}) + \sqrt{3} \frac{\mathbb{P}_{nt}}{\mathbb{P}_{nn}} d(c_{nn}) = 0 \quad \Rightarrow \quad \vec{\mathbf{u}}\cdot\vec{\mathbf{t}} + \sqrt{3} \frac{\mathbb{P}_{nt}}{\mathbb{P}_{nn}} c_{nn} = c_{ste}. \qquad (81)$$

On a maintenant toutes les relations nécessaires à la résolution complète du problème de Riemann pour le système du Gaussien Anisotrope dans une direction $\vec{n}$ quelconque.

### 2.4.6 Résolution complète du problème de Riemann pour le système Gaussien Anisotrope

On cherche à résoudre (47), pour $x \in \mathbb{R}$ et $t \ge 0$, où les variables conservatives et leurs flux sont donnés par (41). La condition initiale est

$$\mathbf{W}(t=0, x) = (x > 0)?\mathbf{W}_R : \mathbf{W}_L.$$

Par l'étude de la structure d'onde, on sait que le problème est autosimilaire, faisant apparaître 5 ondes (puisque les ondes 3 et 4 sont confondus), séparant 6 états constants, qui dans l'ordre sont :

$$\mathbf{W}_L, \mathbf{W}_L^*, \mathbf{W}_L^0, \mathbf{W}_R^0, \mathbf{W}_R^*, \mathbf{W}_R.$$

Dans un premier temps, nous avons donné les structures propres associées à ces ondes, puis nous avons montré que les trois ondes internes sont linéairement dégénérées, tandis que les deux ondes externes sont soit des chocs, soit des détentes. On rappelle également les trois relations de saut, $\vec{v} = \vec{u} - \vec{\sigma}$, $j = \rho\vec{v}\cdot\vec{n}$ :
— Masse : $[j] = 0;$

— QdM : $j[\vec{u}] + [P(\vec{n})] = 0;$

— Énergie : $j[\vec{u} \otimes \vec{u} + C] + 2[\vec{u} \odot P(\vec{n})] = 0.$

On peut maintenant grandement réduire le nombre d'inconnues du problème qui sont pour l'instant au nombre de $4 \times 6 = 24$.

1. Pour $\mathcal{C}_0$: par l'expression de $\vec{r}_3$ et $\vec{r}_4$, voir (38), on déduit que

$$[\vec{u}] = 0 \qquad (82)$$

puis que

$$[j] = 0 \qquad (83)$$

et

$$[P(\vec{n})] = 0 \qquad (84)$$

Il en résulte qu'à travers la discontinuité centrale, seuls $\rho$ et $P_{tt}$ subissent un saut.

2. Pour $\mathcal{C}_\pm^{ld}$ : par la forme de $\vec{r}_\pm^{ld}$, il vient

$$[\rho] = 0 \quad \text{et} \quad [\vec{u} \cdot \vec{n}] = 0. \qquad (85)$$

En projetant la relation de saut sur la quantité de mouvement sur $\vec{n}$, on obtient subsidiairement

$$[P_{nn}] = 0. \qquad (86)$$

Les projections de la relation de saut sur la quantité de mouvement sur $\vec{n}$ et de celle sur l'énergie sur $(\vec{n}, \vec{t})$ donnent un système en $[\vec{u} \cdot \vec{t}]$ et $[P_{nt}]$, qu'on peut montrer singulier. La projection de la relation de saut en énergie sur $(\vec{t}, \vec{t})$ fournit une relation explicite entre $[P_{tt}]$ et $[\vec{u} \cdot \vec{t}]$ et $[P_{nt}]$.

À travers les deux ondes linéairement dégénérées conjuguées, on a un saut de $\vec{u} \cdot \vec{t}$, $P_{nt}$ et $P_{tt}$. De plus, puisque les valeurs propres $\lambda_\pm^{ld}$ se conservent à travers ces ondes, les caractéristiques associées sont alignées avec la discontinuité et on peut montrer que $j = \mp \rho c_{nn}$.

On a alors réduit le nombre d'inconnues à 14, comme indiqué sur la Figure 2. De manière classique, on va maintenant pouvoir résoudre itérativement l'intersection des courbes de 1– et 6–choc pour déterminer $\vec{u} \cdot \vec{n}^*$ et $P_{nn}^*$ qui sont constants dans l'éventail intérieur du problème de Riemann. Puis on en déduira de proche en proche les autres quantités inconnues.

3. Pour $\mathcal{C}_\pm$ :

a) Cas du choc : Puisque les relations de saut sont les mêmes que pour les équations d'Euler avec $\gamma = 3$ et $c = \sqrt{3}c_{nn}$, la relation qui relie $\vec{u} \cdot \vec{n}^*$ et $P_{nn}^*$ dans le cas d'un $6-1$-choc est donnée par :

$$ \begin{aligned} \vec{u} \cdot \vec{n}^* &= \vec{u} \cdot \vec{n}_{R/L} \pm \frac{c_{nnR/L}}{\sqrt{P_{nnR/L}}} \frac{P_{nn}^* - P_{nnR/L}}{\sqrt{2P_{nn}^* + P_{nnR/L}}} \\ &= \vec{u} \cdot \vec{n}_{R/L} \pm \frac{1}{\sqrt{\rho_{R/L}}} \frac{P_{nn}^* - P_{nnR/L}}{\sqrt{2P_{nn}^* + P_{nnR/L}}}, \end{aligned} \qquad (87) $$
sous les conditions de validité $\mathbb{P}_{nn}^* > \mathbb{P}_{nnL}$ pour la 1-onde et $\mathbb{P}_{nn}^* < \mathbb{P}_{nnR}$ pour la 6-onde. Une fois le problème résolu en $\vec{u} \cdot \vec{n}^*$ et $\mathbb{P}_{nn}^*$, on en déduit $\rho_R^*$ et $\rho_L^*$ par la relation de saut

$$ \frac{\rho_{R/L}^*}{\rho_{R/L}} = \frac{2\mathbb{P}_{nn}^* + \mathbb{P}_{nnR/L}}{\mathbb{P}_{nn}^* + 2\mathbb{P}_{nnR/L}} \qquad (88) $$

b) **Cas de la détente** : Dans ce cadre, la solution restant continue, on sait que $J_\mp$ est un invariant de Riemann à travers l'onde $\mathscr{C}_\pm$ et que les deux entropies $s_3$ et $s_4$ se conservent. On en déduit donc que

$$ \frac{c_{nnR/L}^*}{c_{nnR/L}} = \frac{\rho_{R/L}^*}{\rho_{R/L}} = \left( \frac{\mathbb{P}_{nn}^*}{\mathbb{P}_{nnR/L}} \right)^{1/3} \qquad (89) $$

et

$$ \begin{aligned} \vec{u} \cdot \vec{n}^* &= \vec{u} \cdot \vec{n}_{R/L} \pm \sqrt{3}(c_{nnR/L}^* - c_{nnR/L}) \\ &= \vec{u} \cdot \vec{n}_{R/L} \pm \sqrt{3}c_{nnR/L} \left( \left( \frac{\mathbb{P}_{nn}^*}{\mathbb{P}_{nnR/L}} \right)^{1/3} - 1 \right). \end{aligned} \qquad (90) $$

Les conditions de validité de ces relations sont complémentaires de celles du choc, c'est-à-dire : $\mathbb{P}_{nn}^* < \mathbb{P}_{nnL}$ pour la 1-onde et $\mathbb{P}_{nn}^* > \mathbb{P}_{nnR}$ pour la 6-onde.

De la même manière que dans le cas des équations d'Euler, on peut montrer que les raccords des 1- et 6-transformations sont $\mathcal{C}^2$ en $\mathbb{P}_{nn}^* = \mathbb{P}_{nnR/L}$ et que les quantités $\vec{u} \cdot \vec{n}^*$ et $\mathbb{P}_{nn}^*$ peuvent être trouvées par une recherche itérative de type Newton-Raphson. Une fois le problème résolu, on connaît alors la nature de ces 2 ondes et on peut en déduire les quatres quantités.* manquantes :

c) **Cas de la détente** : Par (89), on a $\rho_{R/L}^*$. On utilise alors les invariants de Riemann des $\mathscr{C}_\pm^{ld}$ pour déterminer $\vec{u} \cdot \vec{t}_{R/L}^*$ et $\mathbb{P}_{ntR/L}^*$. (79) donne immédiatement

$$ \mathbb{P}_{ntR/L}^* = \frac{\mathbb{P}_{nn}^*}{\mathbb{P}_{nnR/L}} \mathbb{P}_{ntR/L} \qquad (91) $$

et (80) et (81) donnent :

$$ \vec{u} \cdot \vec{t}_{R/L}^* = \vec{u} \cdot \vec{t}_{R/L} \pm \sqrt{3} \frac{\mathbb{P}_{nt}}{\mathbb{P}_{nn}} (c_{nnR/L}^* - c_{nnR/L}). \qquad (92) $$

On remarquera que la quantité $\frac{\mathbb{P}_{nt}}{\mathbb{P}_{nn}}$ étant constante à travers la détente, on a omis les indices. Enfin, par conservation de l'entropie $s_4$, on obtient $\mathbb{P}_{tt}$.

c) **Cas du choc** : Par la relation de saut (88), on obtient $\rho_{R/L}^*$. On utilise ensuite les relations de saut supplémentaires (67). On voit immédiatement que les deux premières relations sont indépendantes de la troisième. On écrit alors

$$ \llbracket j \vec{u} \cdot \vec{t} + \mathbb{P}_{nt} \rrbracket = 0 \quad \text{et} \quad \llbracket (j \vec{u} \cdot \vec{n} + \mathbb{P}_{nn}) \vec{u} \cdot \vec{t} + \left( \frac{j}{\rho} + \vec{u} \cdot \vec{n} \right) \mathbb{P}_{nt} \rrbracket = 0. \qquad (93) $$
Ce système s'écrit

$$
\mathcal{M}_{R/L}^{*} \mathbf{W}_{R/L}^{*} = \mathcal{M}_{R/L} \mathbf{W}_{R/L},
$$

avec

$$
\mathcal{M} = \begin{pmatrix} j & 1 \\ j \vec{u} \cdot \vec{n} + \mathbb{P}_{nn} & \frac{j}{\rho} + \vec{u} \cdot \vec{n} \end{pmatrix}.
$$

Le déterminant de ce système vaut

$$
\begin{align}
\det \mathcal{M}^* &= \frac{j}{\rho^*} - \mathbb{P}_{nn}^* \quad \left( j = \rho(\vec{u} \cdot \vec{n} - \sigma) = \mp \sqrt{3\rho c_{nn}} M_{nn} \right), \nonumber \\
&= \rho^* (c_{nn}^*)^2 (3(M_{nn}^*)^2 - 1) \quad \left( 3(M_{nn}^*)^2 - 1 = \frac{4}{3M_1^2 - 1} = 2 \frac{\mathbb{P}_{nnR/L}}{\mathbb{P}_{nn}^*} \right), \nonumber \\
&= 2\mathbb{P}_{nnR/L} > 0. \tag{94}
\end{align}
$$

Le système est toujours inversible et de plus, par la relation de saut de la quantité
de mouvement sur la direction $\vec{n}$, le terme en bas à gauche est constant à travers
la transformation. On trouve alors facilement

$$
\vec{u} \cdot \vec{t}_{R/L}^* = \vec{u} \cdot \vec{t}_{R/L} + \frac{\vec{u} \cdot \vec{n}^* - \vec{u} \cdot \vec{n}_{R/L}}{\mathbb{P}_{nnR/L}} \mathbb{P}_{ntR/L}, \quad (95)
$$

et

$$
\mathbb{P}_{ntR/L}^* = \mathbb{P}_{ntR/L} \frac{\mathbb{P}_{nn}^*}{\mathbb{P}_{nnR/L}}. \quad (96)
$$

Enfin, la troisième relation de saut de (67) donne $\mathbb{P}_{ttR/L}^*$.

4. État $\mathbf{W}_{R/L}^\circ$ : On a dorénavant complètement caractérisé l'état $\mathbf{W}_{R/L}^*$ et par suite les quantités $\rho_{R/L}^*$, $\vec{u} \cdot \vec{n}^*$ et $\mathbb{P}_{nn}^*$ qui restent constantes à travers les $\mathscr{C}_\pm^{ld}$. Il reste donc à calculer $\vec{u} \cdot \vec{t}^\circ$, $\mathbb{P}_{nt}^\circ$ et $\mathbb{P}_{ttR/L}^\circ$. Pour $\vec{u} \cdot \vec{t}^\circ$ et $\mathbb{P}_{nt}^\circ$, comme dans le cas des $\mathscr{C}_\pm$, on a les deux relations de saut (93). Cependant, ces deux relations sont liées puisque $j = \mp \rho c_{nn}$ et par suite $\det \mathcal{M} = 0$. Il en résulte alors que $\vec{u} \cdot \vec{t}^\circ$ et $\mathbb{P}_{nt}^\circ$ sont solutions du sytème

$$
j_{\pm} \vec{u} \cdot \vec{t}^{\circ} + \mathbb{P}_{nt}^{\circ} = j_{\pm} \vec{u} \cdot \vec{t}_{R/L}^{*} + \mathbb{P}_{nt_{R/L}}^{*}, \quad (97)
$$

qui est évidemment inversible et donne :

$$
\vec{u} \cdot \vec{t}^{\circ} = \frac{j_{+} \vec{u} \cdot \vec{t}_{R}^{*} - j_{-} \vec{u} \cdot \vec{t}_{L}^{*} + \mathbb{P}_{ntR}^{*} - \mathbb{P}_{ntL}^{*}}{j_{+} - j_{-}}, \quad (98)
$$

$$
\mathbb{P}_{nt}^{\circ} = \frac{j_{+}j_{-}(\vec{u} \cdot \vec{t}_{L}^{*} - \vec{u} \cdot \vec{t}_{R}^{*}) + j_{+}\mathbb{P}_{ntL}^{*} - j_{-}\mathbb{P}_{ntR}^{*}}{j_{+} - j_{-}}. \quad (99)
$$

Enfin, la dernière relation de saut utilisée à travers les ondes $\mathscr{C}_\pm^{ld}$ :

$$
j \ll \frac{\mathbb{P}_{tt}}{\rho} + \ll \vec{u} \cdot \vec{t} (j \vec{u} \cdot \vec{t} + 2\mathbb{P}_{nt}) \gg = 0 \quad (100)
$$

donne les deux dernières inconnues, $\mathbb{P}_{ttR/L}^\circ$.

Références

[1] Christophe Berthon. Numerical approximations of the 10-moment gaussian closure. *Mathematics of Computation*, 75(256) :pp. 1809–1831, 2006.
<figure><img src="image_3.png" /><figcaption>FIGURE 2 – Représentation schématique de la solution du problème de Riemann 1D pour le système du Gaussien Anisotrope après réduction des inconnues du problème. L'onde centrale en pointillé noir est une discontinuité de contact. Les deux ondes rectilignes en bleu sont également des discontinuités de contact associées aux valeurs propres 2 et 5. Enfin les deux ondes extérieures sont associées aux champs 1 et 6, vraiments non-linéaires; elles peuvent être soit des chocs, soit des détentes.</figcaption></figure>
# A Annexe : Moments Succesifs du Terme de Couplage

On s'intéresse au terme de droite de l'équation (4). C'est un tenseur d'ordre $k$ dont on peut regarder le terme $(i_1, \dots, i_k)$ :

$$ 
\begin{align*}
\left(-\int_{\mathbb{R}^d} (\otimes^{k} \vec{c}) \, \partial_{\vec{c}} (\vec{F} f) \, d\vec{c}\right)_{i_1, \dots, i_k} &= -\int_{\mathbb{R}^d} \prod_{q=1}^{k} c_{i_q} \sum_{j=1}^{d} \frac{\partial F_j f}{\partial c_j} d\vec{c} \\
&= \sum_{j=1}^{d} \int_{\mathbb{R}^d} \frac{\partial (\prod_{q=1}^{k} c_{i_q})}{\partial c_j} \cdot F_j f \, d\vec{c} \\
&= \int_{\mathbb{R}^d} \sum_{p=1}^{k} \left(\prod_{q \neq p} c_{i_q}\right) \cdot (F_{i_p} f) \, d\vec{c} \\
&= \int_{\mathbb{R}^d} \frac{1}{(k-1)!} \sum_{\sigma \in S(k)} \left( (\otimes^{k-1} \vec{c}) \otimes \vec{F} f \right)_{\sigma(i_1, \dots, i_k)} \, d\vec{c} \\
&= k \int_{\mathbb{R}^d} \left( (\otimes^{k-1} \vec{c}) \odot \vec{F} f \right)_{i_1, \dots, i_k} \, d\vec{c}
\end{align*}
$$

# B Annexe : Tenseurs Totalement Symétriques

On cherche ici à calculer la dimension du sous-espace vectoriel des tenseurs totalement symétriques. Dans un tenseur totalement symétrique d'ordre $k$ et de dimension $d$, deux composantes dont les indices ne diffèrent que par une permutation ont même valeur :

$$ \forall (i_1, \dots, i_k) \in [1; d]^k, \forall \sigma \in S(k), \quad \mathbb{D}_{i_1, \dots, i_k} = \mathbb{D}_{\sigma(i_1, \dots, i_k)} $$

On peut alors ranger tous les indices d'une même composante les uns à la suite des autres, dans un ordre croissant. Par exemple :

$$ \mathbb{D}_{31311} = \mathbb{D}_{11133}. $$

Pour dénombrer les composantes distinctes d'un tel tenseur, il faut donc créer des séquences $n_1, \dots, n_k$, où $n_i$ désigne le nombre d'occurence de l'indice $i$ dans les indices de la composante. On peut alors représenter chaque composante indépendante du tenseur par une chaîne constituée uniquement de 0 et de 1 de la manière suivante : on écrit autant de 1 que l'indice 1 apparaît de fois, puis on écrit 0 qui joue le rôle de séparateur, puis on écrit autant de 1 que l'indice 2 apparaît de fois, puis on écrit 0 pour séparer, etc... Il faut donc placer $k$ chiffres 1 et $d-1$ séparateurs 0, dans un total de $d+k-1$ cases. Les 0 sont placés de manière aléatoire, si deux 0 apparaissent consécutivement, cela veut dire qu'un indice est absent, comme l'indice 2 dans l'exemple ci-dessus. Le nombre de combinaisons possibles est donc le nombre de possibilités pour ranger $d+k-1$ objets (des 0 et des 1) dans $d+k-1$ cases, soit $(d+k-1)!$, mais tous les 0 et tous les 1 étant équivalents, le nombre de combinaisons est :

$$ \dim = \frac{(d+k-1)!}{k!(d-1)!} = C_{d+k-1}^{k}. $$

C'est le nombre de composantes indépendantes d'un tenseur totalement symétrique de rang $k$ dans un espace à $d$ dimensions.
C Annexe : Formes linéaires propres des matrices Jacobiennes

$$
3c_{11}^2 L_A = \begin{pmatrix}
\frac{u}{2} (u + \sqrt{3}c_{11}) & -u - \frac{\sqrt{3}}{2} c_{11} & 0 & 1 & 0 & 0 \\
-\frac{3c_{11}}{2} (u + c_{11}) \left(v - u - \frac{c_{12}^2}{c_{11}^2}\right) & \frac{3c_{11}}{2} \left(v - 2u - \frac{c_{12}^2}{c_{11}^2} - \frac{c_{12}^2}{c_{11}}\right) & \frac{3c_{11}}{2} (u + c_{11}) & \frac{3c_{12}^2}{c_{11}} & -3c_{11} & 0 \\
3c_{11}^2 - u^2 & 2u & 0 & -2 & 0 & 0 \\
\frac{1}{c_{11}^2} (u|C| + 3c_{12}^2 (vc_{11}^2 - uc_{12}^2)) & 3 (uc_{12}^2 - vc_{11}^2) & \frac{3c_{11}}{2} (u - c_{11}) & \frac{3c_{12}^2 - |C|}{c_{11}^2} & -6c_{12}^2 & 3c_{11}^2 \\
\frac{3c_{11}}{2} (u - c_{11}) \left(v - u - \frac{c_{12}^2}{c_{11}^2}\right) & -v + 2u \frac{c_{12}^2}{c_{11}^2} - \frac{c_{12}^2}{c_{11}} & -\frac{3c_{11}}{2} (u - c_{11}) & -\frac{3c_{12}^2}{c_{11}} & 3c_{11} & 0 \\
\frac{u}{2} (u - \sqrt{3}c_{11}) & -u + \frac{\sqrt{3}}{2} c_{11} & 0 & 1 & 0 & 0
\end{pmatrix} \quad (101)
$$

$$
3c_{22}^2 L_B = \begin{pmatrix}
\frac{v}{2} (v + \sqrt{3}c_{22}) & 0 & -v - \frac{\sqrt{3}}{2} c_{22} & 0 & 0 & 1 \\
-\frac{3c_{22}}{2} (v + c_{22}) \left(u - v - \frac{c_{22}^2}{c_{22}^2}\right) & \frac{3c_{22}}{2} (v + c_{22}) & \frac{3c_{22}}{2} \left(u - 2v - \frac{c_{22}^2}{c_{22}^2} - \frac{c_{22}^2}{c_{22}}\right) & 0 & -3c_{22} & \frac{3c_{12}^2}{c_{22}} \\
\frac{1}{2c_{22}^2} (3(vc_{12}^2 - uc_{22}^2)^2 - v^2|C|) & 3(vc_{12}^2 - uc_{22}^2) & \frac{1}{c_{22}^2} (v|C| + 3c_{12}^2 (uc_{22}^2 - vc_{12}^2)) & 3c_{22}^2 & -6c_{12}^2 & \frac{3c_{12}^2 - |C|}{c_{12}^2} \\
\frac{3c_{22}}{2} (v - c_{22}) \left(u - v - \frac{c_{12}^2}{c_{22}^2}\right) & -\frac{3c_{22}}{2} (v - c_{22}) & \frac{3c_{22}}{2} \left(-u + 2v - \frac{c_{12}^2}{c_{22}^2} - \frac{c_{12}^2}{c_{22}^2}\right) & 0 & 3c_{22} & -\frac{3c_{12}^2}{c_{22}} \\
\frac{v}{2} (v - \sqrt{3}c_{22}) & 0 & -v + \frac{\sqrt{3}}{2} c_{22} & 0 & 0 & 1
\end{pmatrix} \quad (102)
$$
$$
\begin{align}
\mathbf{r}_{6+1/2}^{-} &= \left( \begin{array}{c} 1, \\ \vec{u} \pm \frac{\sqrt{3}\mathcal{C}(\vec{n})}{\sqrt{\mathcal{C}(\vec{n})\cdot\vec{n}}}, \end{array} \right. \left( \begin{array}{c} \frac{1}{2} \left( \vec{u} \pm \frac{\sqrt{3}\mathcal{C}(\vec{n})}{\sqrt{\mathcal{C}(\vec{n})\cdot\vec{n}}} \right)^{\otimes 2} \\ + \frac{c_{nn}^2 \mathcal{C} - \mathcal{C}(\vec{n}) \otimes \mathcal{C}(\vec{n})}{2c_{nn}^2} \end{array} \right)^t \nonumber \\
&= \left( \begin{array}{c} 0, \\ \vec{t} \end{array} \right. , \left. \begin{array}{c} \vec{u} \pm \frac{\mathcal{C}(\vec{n})}{\sqrt{\mathcal{C}(\vec{n})\cdot\vec{n}}} \odot \vec{t} \end{array} \right)^t \tag{27}
\end{align}
$$

$$
\begin{align*}
\mathbf{r}_3 &= \left( \begin{array}{c} 1, \\ \vec{u} \end{array} \right. , \left. \begin{array}{c} \frac{\vec{u} \otimes \vec{u}}{2} \\ \vec{t} \otimes \vec{t} \end{array} \right)^t \\
\mathbf{r}_4 &= \left( \begin{array}{c} 0, \\ \vec{0} \end{array} \right. , \left. \begin{array}{c} \vec{t} \otimes \vec{t} \end{array} \right)^t
\end{align*}
$$

$$
\begin{align}
l_{6+1}^{-} &= \frac{1}{3c_{nn}^2} \begin{pmatrix} \frac{\vec{u} \cdot \vec{n}}{2} (\vec{u} \cdot \vec{n} \mp \sqrt{3c_{nn}}) & -\left(\vec{u} \cdot \vec{n} \mp \frac{\sqrt{3}}{2} c_{nn}\right) \vec{n} & \vec{n} \otimes \vec{n} \\ \frac{1}{2c_{nn}^3} (\mathcal{C}(\vec{n}) \cdot \vec{u}^{\perp}) (c_{nn} \mp \vec{u} \cdot \vec{n}) & (c_{nn} \mp \vec{u} \cdot \vec{n}) \mathcal{C}(\vec{n})^{\perp} \pm (\mathcal{C}(\vec{n}) \cdot \vec{u}^{\perp}) \vec{n} & \pm 2 \mathcal{C}(\vec{n})^{\perp} \odot \vec{n} \\ \frac{1}{3c_{nn}^2} 3 \mathcal{C}(\vec{n}) \cdot \vec{n} - (\vec{u} \cdot \vec{n})^2 & 2(\vec{u} \cdot \vec{n}) \vec{n} & -2(\vec{n} \otimes \vec{n}) \end{pmatrix} \nonumber \\
l_{5+1/2}^{-} &= \frac{1}{2c_{nn}^3} \begin{pmatrix} \frac{1}{3c_{nn}^2} (3(\mathcal{C}(\vec{n}) \cdot \vec{u}^{\perp})^2 - (\vec{u} \cdot \vec{n})^2 |C|) & (\vec{u} \cdot \vec{n} |C|) \vec{n} + 3 (\mathcal{C}(\vec{n}) \cdot \vec{u}^{\perp}) \mathcal{C}(\vec{n})^{\perp} & 3 (\mathcal{C}(\vec{n})^{\perp} \otimes \mathcal{C}(\vec{n})^{\perp}) - |C| \vec{n} \otimes \vec{n} \end{pmatrix} \tag{103} \\
l_3 &= \frac{1}{3c_{nn}^2} \begin{pmatrix} \frac{1}{3c_{nn}^2} (3(\mathcal{C}(\vec{n}) \cdot \vec{u}^{\perp})^2 - (\vec{u} \cdot \vec{n})^2 |C|) & (\vec{u} \cdot \vec{n} |C|) \vec{n} + 3 (\mathcal{C}(\vec{n}) \cdot \vec{u}^{\perp}) \mathcal{C}(\vec{n})^{\perp} & 3 (\mathcal{C}(\vec{n})^{\perp} \otimes \mathcal{C}(\vec{n})^{\perp}) - |C| \vec{n} \otimes \vec{n} \end{pmatrix} \nonumber \\
l_4 &= \frac{1}{3c_{nn}^4} \begin{pmatrix} \frac{1}{2} (3(\mathcal{C}(\vec{n}) \cdot \vec{u}^{\perp})^2 - (\vec{u} \cdot \vec{n})^2 |C|) & (\vec{u} \cdot \vec{n} |C|) \vec{n} + 3 (\mathcal{C}(\vec{n}) \cdot \vec{u}^{\perp}) \mathcal{C}(\vec{n})^{\perp} & 3 (\mathcal{C}(\vec{n})^{\perp} \otimes \mathcal{C}(\vec{n})^{\perp}) - |C| \vec{n} \otimes \vec{n} \end{pmatrix} \nonumber
\end{align}
$$
# D Annexe : Invariants de Riemann

On considère un système hyperbolique de loi de conservation en 1D, mis sous une forme quasi-linéaire :

$$ \partial_t \mathbf{W} + A(\mathbf{W}) \partial_x \mathbf{W} = 0, \quad \mathbf{W} \in \Omega \subset \mathbb{R}^d. \qquad (104) $$

Puisque $A(\mathbf{W})$ est partout diagonalisable, on peut multiplier à gauche par les formes linéaires

$$ \mathcal{L}(\mathbf{W}) \partial_t \mathbf{W} + \Lambda(\mathbf{W}) \mathcal{L}(\mathbf{W}) \partial_t \mathbf{W} = 0, $$

soit, pour chacun des champs propre $i$, $i \in \int 1, d$,

$$ \vec{l}_i (\partial_t \mathbf{W}) + \lambda_i \vec{l}_i (\partial_x \mathbf{W}) = 0 \Rightarrow \sum_k l_{ik} \partial_{\mathcal{C}_i} \mathbf{W}_k = 0 \Rightarrow \vec{l}_i (\partial_{\mathcal{C}_i} \mathbf{W}) = 0 \qquad (105) $$

Si l'on peut intégrer cette dernière relation, on obtient alors l'invariant de Riemann pour le $i$-champs $J_i$, tel que

$$ \partial_{\mathcal{C}_i} (J_i) = \vec{l}_i (\partial_{\mathcal{C}_i} \mathbf{W}) = \partial_{\mathbf{W}} J_i \cdot \partial_{\mathcal{C}_i} \mathbf{W}, \qquad (106) $$

la dernière égalité étant obtenue par dérivation composée. Comme cette dernière relation est valable dans tous l'espace des phase $\Omega$, il en résulte par le théorème de Riesz que la forme linéaire $\vec{l}_i$ s'identifie au vecteur $\partial_{\mathbf{W}} J_i$ et par suite :

$$ \forall \mathbf{W} \in \Omega, \quad \vec{\nabla}_{\mathbf{W}} J_i \cdot \vec{r}_k = \vec{l}_i(\vec{r}_k) = \delta_{ik}, $$

ce qui signifie bien que $J_i$ se conserve le long de toutes les courbes intégrales autres que $\mathcal{I}_i$, et donc à travers toutes les ondes simples autres que celle associée à son propre champs.