Quelles méthodes peuvent garantir que les quantités physiques restent positives tout au long d'une simulation PDE?

Les quantités physiques comme la pression, la densité, l'énergie, la température et la concentration doivent toujours être positives, mais les méthodes numériques calculent parfois des valeurs négatives au cours du processus de solution. Ce n'est pas correct car les équations calculent des valeurs complexes ou infinies (écrasant généralement le code). Quelles méthodes numériques peuvent être utilisées pour garantir que ces quantités restent positives? Laquelle de ces méthodes est la plus efficace?

La méthode la plus courante consiste à réinitialiser les valeurs négatives à un petit nombre positif. Bien sûr, ce n'est pas une solution mathématiquement valable. Une meilleure approche générale qui peut fonctionner et qui est facile consiste à réduire la taille de votre pas de temps.

Les valeurs négatives surviennent souvent dans la solution des EDP hyperboliques, car l'apparition de chocs peut conduire à des oscillations, qui auront tendance à créer des valeurs négatives s'il y a des états de quasi-vide près du choc. L'utilisation d'une méthode de diminution de variation totale (TVD) ou d'une autre méthode non oscillatoire ( ENO, WENO ) peut réduire cette tendance. Ces méthodes sont basées sur l'utilisation de limiteurs non linéaires pour calculer les dérivées de la solution. Cependant, vous pouvez toujours obtenir des valeurs négatives pour plusieurs raisons:

• Si vous utilisez la méthode des lignes et appliquez un intégrateur de temps d'ordre élevé. La plupart des schémas TVD sont en fait TVD uniquement dans le sens semi-discret ou avec la méthode d'Euler. Pour une intégration temporelle d'ordre supérieur, vous devez utiliser une discrétisation temporelle à forte stabilité (SSP) ; ces schémas sont également appelés «contractuels» ou «préservant la monotonie». Il y a un livre récent sur le sujet par Sigal Gottlieb, Chi-Wang Shu et moi-même.
• Si vous n'utilisez pas la décomposition caractéristique locale pour les systèmes d'équations, votre solution ne sera pas TVD (les schémas TVD ne possèdent cette propriété que pour les problèmes scalaires). Il est donc préférable de reconstruire / interpoler en variables caractéristiques.
• Si vous avez un système non linéaire, des valeurs négatives peuvent survenir même si vous utilisez la décomposition caractéristique locale. Par exemple, tout solveur de Riemann linéarisé (tel qu'un solveur Roe) pour les équations en eau peu profonde ou les équations d'Euler peut générer des valeurs négatives dans des conditions suffisamment difficiles. Une solution consiste à utiliser un solveur HLL (ou une variante de HLL); certains d'entre eux sont prouvés positifs.
• Les schémas TVD ne sont que du second ordre; les schémas non oscillatoires d'ordre supérieur comme WENO ne satisfont pas strictement la TVD ou les principes maximaux. Mais une nouvelle modification de ces schémas d'ordre supérieur le fait; il est développé dans plusieurs articles récents par Xiangxiong Zhang (un étudiant de Chi-Wang Shu).

Il existe, bien sûr, de nombreuses autres approches spécialisées pour des équations particulières, comme dans le code GeoClaw de David George, qui utilise un solveur de Riemann avec des ondes non physiques supplémentaires pour imposer la positivité.

En supposant que nous résolvons des équations hyperboliques sans aucun terme source et en supposant que nous fournissons des conditions physiques initiales, nous assurer que le schéma numérique que nous utilisons est la réduction totale de la variation est un bon moyen d'assurer la «physicalité» de la solution calculée. Puisqu'un schéma TVD préserve la monotonie, aucun nouveau minimum ou maximum ne sera créé et la solution restera limitée par les valeurs initiales que nous espérons définir correctement. Bien sûr, le problème est que les schémas TVD ne sont pas les plus évidents. Parmi les schémas linéaires, seuls les schémas de premier ordre sont TVD (Godunov 1954). Ainsi, depuis les années 50, une variété de schémas TVD non linéaires ont été développés pour combiner haute précision et monotonie pour la résolution d'équations hyperboliques.

Pour mes applications, résolvant des équations de Navier-Stokes avec de grands gradients de pression / densité, nous utilisons un schéma hybride MUSCL- central pour capturer les grands gradients / discontinuités et en conserver une bonne précision. Le premier programme MUSCL (MUSCL signifie Monotone Upstream-center Schemes for Conservation Laws) a été conçu par Van Leer en 1979.

Si vous souhaitez en savoir plus sur ce sujet, veuillez consulter les œuvres de Harten, Van Leer, Lax, Sod et Toro.

Les réponses ci-dessus s'appliquent aux problèmes liés au temps, mais vous pouvez également exiger la positivité dans une simple équation elliptique. Dans ce cas, vous pouvez le formuler comme une inégalité variationnelle , en donnant des limites pour les variables.

Dans PETSc, il existe deux solveurs VI. On utilise une méthode à espace réduit, où les variables dans les contraintes actives sont supprimées du système à résoudre. L'autre utilise une méthode de Newton semi-lisse .

$A$

$$Au = b$$ $A$$A^{-1}$

$B\in\mathbb{R}^{n \times n}$$B\geq0$$B$

$$(B\geq 0) \qquad\Leftrightarrow\qquad (u\leq v ~\Rightarrow~ Bu\leq Bv, ~~\forall u,v \in \mathbb{R}^n )$$

$A$

$$0 \leq b ~\Rightarrow~ 0=A^{-1}0 \leq A^{-1}b = u$$ $b$$b\geq 0$

Généralement, les schémas de discrétisation qui conduisent à une matrice M sont appelés schémas monotones et ce sont ces schémas qui préservent la non-négativité.