Une bonne différence finie pour l'équation de continuité

Quelle serait une bonne discrétisation par différence finie pour l'équation suivante:

$\frac{\partial \rho}{\partial t} + \nabla \cdot \left(\rho u\right)=0$ ?

Nous pouvons prendre le cas 1D:

$\frac{\partial \rho}{\partial t} + \frac{d}{dx}\left(\rho u\right)=0$

Pour une raison quelconque, tous les schémas que je peux trouver concernent la formulation en coordonnées lagrangiennes. Je suis venu avec ce schéma pour le moment (ne tenez pas compte de l' indice j ):

$\frac{\rho^{n+1}_{i,j}-\rho^{n}_{i,j}}{\tau} + \frac{1}{h_x}\left(\frac{\rho^{n+1}_{i+1,j}+\rho^{n+1}_{i,j}}{2}u^{n}_{_xi+1/2,j}- \frac{\rho^{n+1}_{i,j}+\rho^{n+1}_{i-1,j}}{2}u^{n}_{_xi-1/2}\right)=0$

Mais semble être vraiment instable ou avoir une horrible condition de stabilité. Est-ce vrai?

La vitesse est en fait calculée par la loi de Darcy . De plus, nous avons l'équation d'état. Le système complet se compose également d'une équation énergétique et de l'équation d'état pour le gaz idéal. Les vitesses peuvent devenir négatives .$u=-\frac{k}{\mu}\nabla p$

Vous regardez l'équation de conservation de masse:

$\dfrac{dm}{dt}=0$

Lorsque l'on considère l'évolution de la masse par unité de volume, cela se résume à l'équation d'advection de la densité sous forme de flux:

$\dfrac{\partial \rho}{\partial t} = -\nabla \cdot (\rho u)$

La bonne chose à ce sujet est que ce n'est que l'équation d'advection d'un champ scalaire arbitraire (dans notre cas, il s'agit de la densité ) et qu'il est (relativement) facile à résoudre, à condition de disposer de schémas de différenciation temporelle et spatiale adéquats, et d'initiales et les conditions aux limites.$\rho$

Lors de la conception d'un schéma de différenciation finie, nous nous soucions de la convergence, de la stabilité et de la précision. Un schéma converge si lorsque . La stabilité des schémas garantit que la quantité reste finie lorsque . La précision formelle du schéma indique où se situe l'erreur de troncature dans la série d'expansion de Taylor de la dérivée partielle. Consultez un manuel CFD pour plus de détails sur ces propriétés fondamentales d'un schéma de différenciation. Δt→0At→$\dfrac{\Delta A}{\Delta t} \rightarrow \dfrac{\partial A}{\partial t}$$\Delta t \rightarrow 0$$A$$t \rightarrow \infty$

Maintenant, l'approche la plus simple consiste à passer directement à la différenciation amont du 1er ordre. Ce schéma est positif-défini, conservateur et efficace en termes de calcul. Les deux premières propriétés sont particulièrement importantes lorsque l'on modélise l'évolution d'une quantité toujours positive (c'est-à-dire masse ou densité).

Pour plus de simplicité, regardons le cas 1-D:

$\dfrac{\partial \rho}{\partial t} = -\dfrac{\partial(\rho u)}{\partial x}$

Il convient maintenant de définir le flux , de sorte que:$\Phi = \rho u$

$\dfrac{\partial(\rho u)}{\partial x} = \dfrac{\partial \Phi}{\partial x} \approx \dfrac{\Delta \Phi}{\Delta x} \approx \dfrac{\Phi_{i+1/2}-\Phi_{i-1/2}}{\Delta x}$

Voici un schéma de ce que nous simulons:

            u           u
|          -->         -->          |
|    rho    |    rho    |    rho    |
x-----o-----x-----o-----x-----o-----x
     i-1  i-1/2   i   i+1/2  i+1

Nous évaluons l'évolution de à la cellule . Le gain ou la perte nette provient de la différence entre ce qui entre, et ce qui sort, . C'est là que nous commençons à diverger de la réponse de Paul. Dans une véritable différenciation conservatrice en amont, la quantité au centre de la cellule est transportée par la vitesse à son bord de cellule, dans la direction de son mouvement. En d'autres termes, si vous imaginez que vous êtes la quantité advectée et que vous êtes assis au centre de la cellule, vous êtes transporté dans la cellule devant vous par la vitesse au bord de la cellule. L'évaluation du flux au bord de la cellule en tant que produit de la densité et de la vitesse, à la fois au bord de la cellule, n'est pas correcte et ne conserve pas la quantité advectée.i Φ $\rho$$i$ Φ i + 1 /$\Phi_{i-1/2}$$\Phi_{i+1/2}$

Les flux entrants et sortants sont évalués comme:

$\Phi_{i+1/2} = \dfrac{u_{i+1/2}+|u_{i+1/2}|}{2}\rho_{i} + \dfrac{u_{i+1/2}-|u_{i+1/2}|}{2}\rho_{i+1}$

$\Phi_{i-1/2} = \dfrac{u_{i-1/2}+|u_{i-1/2}|}{2}\rho_{i-1} + \dfrac{u_{i-1/2}-|u_{i-1/2}|}{2}\rho_{i}$

Le traitement ci-dessus de la différenciation des flux garantit une définition en amont. En d'autres termes, il ajuste la direction de différenciation en fonction du signe de la vitesse.

Le critère de stabilité de Courant-Friedrichs-Lewy (CFL), lors de la différenciation temporelle avec un simple premier ordre, la différenciation d'Euler vers l'avant est donnée comme suit:

$\mu = \dfrac{u\Delta t}{\Delta x} \leq 1$

A noter qu'en 2 dimensions, le critère de stabilité CFL est plus strict:

$\mu = \dfrac{c\Delta t}{\Delta x} \leq \dfrac{1}{\sqrt{2}}$

où est la magnitude de la vitesse, .$c$$\sqrt{u^{2}+v^{2}}$

Quelques choses à considérer. Ce schéma peut ou non être approprié pour votre application en fonction du type de processus que vous simulez. Ce schéma est hautement diffusif et convient aux écoulements très lisses sans gradients prononcés. Il est également plus diffusif pour des pas de temps plus courts. Dans le cas 1-D, vous obtiendrez une solution presque exacte si les gradients sont très petits et si . Dans le cas 2-D, cela n'est pas possible et la diffusion est anisotrope.$\mu = 1$

Si votre système physique prend en compte les ondes de choc ou les gradients élevés d'une autre sorte, vous devez rechercher une différenciation en amont d'ordre supérieur (par exemple, 3e ou 5e ordre). En outre, il peut être utile d'examiner la famille de systèmes de transport à flux corrigé (Zalesak, 1979, JCP); correction anti-diffusion pour le schéma ci-dessus par Smolarkiewicz (1984, JCP); Famille de schémas MPDATA par Smolarkiewicz (1998, JCP).

Pour la différenciation temporelle, la différenciation Euler de premier ordre peut être satisfaisante pour vos besoins. Sinon, examinez les méthodes d'ordre supérieur telles que Runge-Kutta (itérative), ou Adams-Bashforth et Adams-Moulton (multiniveau).

Il serait intéressant de consulter certains manuels de niveau supérieur de CFD pour un résumé des programmes mentionnés ci-dessus et bien d'autres.

Dans le cas 1D, vous ne voulez pas utiliser un schéma de différence directe ou centrale pour le terme dérivé spatial car ils sont numériquement instables. Au lieu de cela, il est préférable de discrétiser l'équation avec une différence finie explicite en arrière (au près) pour la dérivée spatiale:$(\frac{d}{dx})$

$\frac{\rho_i^{k+1}-\rho_i^k}{\Delta t} + \frac{\rho_i^k U_i^k-\rho_{i-1}^k U_{i-1}^k}{\Delta x} = 0$ .

Si les vitesses sont positives, ce schéma en arrière est stable. S'ils sont négatifs, une différence directe fonctionnera. Quoi qu'il en soit, il y a toujours une contrainte sur votre choix de et (nombre courant) pour rendre le schéma stable.$\Delta x$$\Delta t$