Méthode de scission amont du terme d'advection

On considère l'équation hyperbolique scalaire suivante :

${\displaystyle {\frac {\partial u}{\partial t}}+{\frac {\partial f(u)}{\partial x}}=0}$

On peut l'écrire sous forme non conservative :

${\displaystyle {\frac {\partial u}{\partial t}}+a{\frac {\partial u}{\partial x}}=0\,,\quad a={\frac {\partial f}{\partial u}}}$

Il s'agit d'une équation de conservation pour ${\displaystyle u}$ où ${\displaystyle a}$ est la vitesse d'advection, constante le long d'une caractéristique. On peut ainsi calculer la solution exacte en remontant depuis le point de l'espace-temps ${\displaystyle (x_{j},t_{n+1})}$ pour chercher le point d'intersection de la caractéristique en ce point ${\displaystyle x_{a}=x_{j}\mp a\Delta t\,,\;\pm a>0}$ avec la ligne de temps ${\displaystyle t_{n}=const}$ :

${\displaystyle u(x_{j},t_{n+1})=u_{a}(x_{a},t_{n})}$

${\displaystyle x_{a}}$ ne coïncidant pas en général avec un point du maillage une interpolation est nécessaire, par exemple une interpolation linéaire respectant la condition CFL:

${\displaystyle u_{a}=\left(1-\alpha \right)u_{j}^{n}+\alpha u_{j\mp 1}^{n}\,,\quad \alpha =|a|{\frac {\Delta t}{\Delta x}}\,,\quad 0<\pm a{\frac {\Delta t}{\Delta x}}\leq 1}$

D'où le schéma de Courant-Isaacson-Rees (CIR) décentré amont en différences finies (1952):

${\displaystyle u_{j}^{n+1}=u_{j}^{n}-a{\frac {\Delta t}{\Delta x}}\left\{{\begin{array}{rcl}(u_{j}^{n}-u_{j-1}^{n})&si&a\geq 0\\(u_{j+1}^{n}-u_{j}^{n})&si&a<0\end{array}}\right.}$

Cette méthode n'est pas capable de capturer une discontinuité (choc, ligne de contact). Pour obtenir un tel résultat, Peter Lax a donné en 1954 une condition nécessaire d'entropie qui garantit l'existence d'une solution faible du problème de Rankine-Hugoniot^[7].

On prend pour exemple les équations d'Euler unidimensionnelles pour un gaz parfait :

${\displaystyle {\frac {\partial \mathbf {U} }{\partial t}}+{\frac {\partial \mathbf {F} (\mathbf {U} )}{\partial x}}=0}$

La version volumes finis d'un solveur de cette équation s'écrit :

${\displaystyle \mathbf {u} _{j}^{n+1}-\mathbf {u} _{j}^{n}=-{\frac {\Delta t}{\Delta x}}\left[\mathbf {f} (x_{j+1/2},t_{n})-\mathbf {f} (x_{j-1/2},t_{n})\right]}$

où ${\displaystyle \mathbf {u} }$ est l'approximation de ${\displaystyle \mathbf {U} }$ et ${\displaystyle \mathbf {f} (x_{k},t_{n})}$ le vecteur flux à l'interface ${\displaystyle x_{k}}$.

En suivant le principe de décentrage, ce flux sera écrit comme une fonction des quantités des cellules voisines à droite et à gauche ${\displaystyle (\mathbf {u} _{L},\mathbf {u} _{R})}$, par exemple :

${\displaystyle \mathbf {f} _{j+1/2}\equiv \mathbf {f} (x_{j+1/2},t_{n})=\mathbf {f} (\mathbf {u} _{j}^{n},\mathbf {u} _{j+1}^{n})\equiv \mathbf {f} (\mathbf {u} _{L},\mathbf {u} _{R})}$

Le calcul des états ${\displaystyle (\mathbf {u} _{L},\mathbf {u} _{R})}$ requiert la résolution d'un problème de Riemann. En l'absence de variables caractéristiques dans les équations on définit une variable ${\displaystyle \xi =(x-x_{j+1/2})/(t-t_{n})}$ jouant le rôle de caractéristique le long de laquelle la solution est constante. L'interface de cellule coïncide avec caractéristique ${\displaystyle \xi =0}$ et le flux numérique s'écrit :

${\displaystyle \mathbf {f} _{j+1/2}=\mathbf {f} (\mathbf {u} _{L},\mathbf {u} _{R})\mathbf {f} \left(\mathbf {w} (\xi =0;\mathbf {u} _{L},\mathbf {u} _{R})\right)}$.

Une telle solution n'est pas toujours possible, elle est coûteuse et sa précision remise en cause à chaque étape de l'intégration temporelle^[8].

Un exemple de résolution approchée est celle de Roe^[9]. Dans cette méthode les flux sont scindés dans deux directions suivant les valeurs propres de la matrice jacobienne ${\displaystyle A}$ évaluées dans un certain état moyen ${\displaystyle {\hat {\mathbf {u} }}(\mathbf {u} _{L},\mathbf {u} _{R})}$ :

${\displaystyle \mathbf {F} (\mathbf {u} _{R})-\mathbf {F} (\mathbf {u} _{L})=A({\hat {u}})(\mathbf {u} _{R}-\mathbf {u} _{L})=(A^{+}({\hat {\mathbf {u} }})+A^{-}({\hat {\mathbf {u} }}))(\mathbf {u} _{R}-\mathbf {u} _{L})}$

D'où le flux :

${\displaystyle {\begin{array}{rcl}\mathbf {f} _{j+1/2}&=&\mathbf {F} (\mathbf {u} _{L})+A^{-}({\hat {\mathbf {u} }})(\mathbf {u} _{R}-\mathbf {u} _{L})\\&=&\mathbf {F} (\mathbf {u} _{R})+A^{+}({\hat {\mathbf {u} }})(\mathbf {u} _{R}-\mathbf {u} _{L})\\&=&{\frac {1}{2}}\left[\mathbf {F} (\mathbf {u} _{L})+\mathbf {F} (\mathbf {u} _{R})+(A^{-}({\hat {\mathbf {u} }})-A^{+}({\hat {\mathbf {u} }}))(\mathbf {u} _{R}-\mathbf {u} _{L})\right]\end{array}}}$

Ce schéma représente des états amont et aval sont reliés par une simple discontinuité (choc ou discontinuité de contact) et ne respecte pas la condition d'entropie croissante : il peut conduire à une discontinuité de type « choc de détente » non physique.

Une autre approche, dans la logique de la méthode CIR, consiste à remonter les sources influençant le flux :

${\displaystyle \mathbf {f} (\mathbf {u} _{L},\mathbf {u} _{R})=\mathbf {F} ^{+}(\mathbf {u} _{L})+\mathbf {F} ^{-}(\mathbf {u} _{R})}$

Cela suppose la connaissance (l'approximation) des matrices jacobiennes associées à ${\displaystyle \mathbf {F} ^{+}}$ et ${\displaystyle \mathbf {F} ^{-}}$.

Parmi toutes les variantes proposées la méthode de van Leer^[3] est sans doute la plus efficace malgré les problèmes rencontrés pour le calcul des discontinuités de contact^[10].