Chapter 6  covofi

• Fonction
• Discrétisation
• Mise en œuvre
• Points à traiter
• Annexe 1 : Inversibilité de la matrice EM _n
• Annexe 2 : Remarques à propos du respect du principe du maximum discret

0.1  Fonction

Dans ce sous-programme, on résout :
soit l’équation de convection-diffusion d’un scalaire en présence de termes sources : Ici a représente la valeur instantanée du scalaire en approche laminaire ou, en approche RANS, sa moyenne de Reynolds a. Les deux approches étant exclusives et les équations obtenues similaires, on utilisera le plus souvent aussi la notation a pour a.
soit, dans le cas d’une modélisation RANS, la variance de la fluctuation d’un scalaire en présence de termes sources¹ : a"² représente ici la moyenne du carré des fluctuations² de a. K, Γ, T_s^imp et T_s^exp représentent respectivement le coefficient de diffusion, la valeur du terme source de masse, les termes sources implicite et explicite du scalaire a ou a"². µ_t et σ_t sont respectivement la viscosité turbulente et le nombre de Schmidt ou de Prandtl turbulent, ε est la dissipation de l’énergie turbulente k et R_f définit le rapport constant entre les échelles dissipatives de k et de a"² (R_f est constant selon le modèle assez simple adopté ici).
On écrit les deux équations précédentes sous la forme commune suivante : avec, pour f=a ou f=a"² :

Le terme ∂ (ρ f)/∂ t est décomposé de la sorte : En utilisant l’équation de conservation de la masse (cf. preduv ), l’équation précédente s’écrit finalement :

0.2  Discrétisation

Pour intégrer l’équation (??), une discrétisation temporelle de type θ-schéma est appliquée à la variable résolue³ :

L’équation (??) est discrétisée au temps n+θ en supposant les termes sources explicites pris au temps n+θ_S, et ceux implicites en n+θ. Par souci de clarté, on suppose, en l’absence d’indication, que les propriétes physiques Φ (K, ρ,...) et le flux de masse (ρ u) sont pris respectivement aux instants n+θ_Φ et n+θ_F, où θ_Φ et θ_F dépendent des schémas en temps spécifiquement utilisés pour ces grandeurs⁴.

où : Le terme de production affecté d’un indice n+θ_S est un terme source explicite et il est donc traité comme tel :
L’équation (??) s’écrit : avec : On rappelle que, pour un scalaire f, le sous-programme codits résout une équation du type suivant f_s^exp représente les termes sources discrétisés de manière explicite en temps (hormis contributions de la convection diffusion explicite provenant du θ-schéma) et f_s^imp fⁿ⁺¹ représente les termes linéaires en fⁿ⁺¹ dans l’équation discrétisée en temps.
On réécrit l’équation (??) sous la forme (??) qui est ensuite résolue par codits .

0.3  Mise en œuvre

On distingue deux cas suivant le type de schéma en temps choisi pour les termes sources :
• Si les temres sources ne sont pas extrapolés, toutes les contributions du second membre vont directement dans le vecteur SMBRS.
• Sinon, un vecteur supplémentaire est nécessaire afin de stocker les contributions du pas de temps précedent (PROPCE). Dans ce cas :

• - le vecteur PROPCE sert à stocker les contributions explicites du second menbre au temps n−1 (pour l’extrapolation en n+θ_S).
• - le vecteur SMBRS est complété au fur et à mesure.
  

L’algorithme de ce sous-programme est le suivant :

• mise à zéro des vecteurs représentants le second membre (SMBRS) et de la diagonale de la matrice (ROVSDT).
• calcul des termes sources du scalaire définis par l’utilisateur en appelant le sous-programme ustssc .
  
  ⋆ Si les termes sources sont extrapolés, SMBRS reçoit −θ_S fois la contribution au temps n−1 des termes sources qui sont extrapolés (stockés dans PROPCE). La contribution des termes sources utisateurs (au pas temps courant) est répartie entre PROPCE (pour la partie T_s^exp qui est à stocker en vue de l’extrapolation) et SMBRS (pour la partie explicite provenant de l’utilisation du θ schéma pour T_s^imp). La contribution implicite est alors mise dans ROVSDT (après multiplication par θ) quel que soit son signe, afin de ne pas utiliser des discrétisations temporelles différentes entre deux pas de temps successifs, dans le cas par exemple où T_s^imp change de signe⁵.
  
  ⋆ Sinon la contibution de T_s^exp est directement mise dans SMBRS. Celle de T_s^imp est ajoutée à ROVSDT si elle est positive (de manière à conserver la dominance de la diagonale), ou explicitée et mise dans le second membre sinon.
  
• prise en compte des physiques particulières : arc électrique, rayonnement, combustion gaz et charbon pulvérisé. Seuls les vecteurs ROVSDT et SMBRS sont complétés (schéma d’ordre 1 sans extrapolation).
• ajout des termes sources de masse en Γ (f_i−f) par appel au sous-programme catsma .
  
  ⋆ Si les termes sources sont extrapolés, le terme explicite en Γ f_i est stocké dans PROPCE. Le θ-schéma est appliqué au terme implicite, puis les contributions implicite et explicite réparties entre ROVSDT et SMBRS.
  
  ⋆ Sinon, la partie implicte en −Γ f va dans ROVSDT, et tout le reste dans SMBRS.
  
• calcul du terme d’accumulation de masse en div(ρ u) par appel à divmas et ajout de sa contribution dans SMBRS, et dans ROVSDT après multiplication par θ⁶.
• ajout du terme instationnaire à ROVSDT.
• calcul des termes de production (2 µ_t/σ_t(grad a)²) et de dissipation ( − ρ ε/R_f ka"²) si on étudie la variance des fluctuations d’un scalaire avec un modèle de turbulence de type RANS. Ce calcul s’effectue en calculant préalablement le gradient du scalaire f par appel au sous-programme grdcel .
  
  ⋆ Si les termes sources sont extrapolés, la production est mise dans PROPCE puis l’énergie cinétique k et la dissipation turbulentes ε sont calculées (XK et XE) en fonction du modèle de turbulence utilisé. SMBRS reçoit − ρ ε/R_f ka"² au temps n et ROVSDT le coefficient d’implicitation ρ ε/R_f k après multiplication par THETAP = θ.
  
  ⋆ Sinon, la production va dans SMBRS, et la dissipation est répartie de la même manière que précédemment avec THETAP = 1.
  
• une fois la contribution de tous les termes sources calculés, le second membre est assemblé, et le vecteur PROPCE ajouté après multiplication par 1+θ_S à SMBRS, dans le cas où les termes sources sont extrapolés.
• calcul du coefficient de diffusion K au centre des cellules, et des valeurs aux faces par appel au sous-programme viscfa .
• résolution de l’équation complète (avec les termes de convection diffusion) par un appel au sous-programme codits avec f_s^exp=SMBRS et f_s^imp=ROVSDT.
• ajustement (clipping) du scalaire ou de la fluctuation du scalaire en appelant le sous-programme clpsca .
• impression du bilan explicite d’expression ||E_n(fⁿ) − ρⁿ/Δ t ( f ⁿ⁺¹ − fⁿ )|| , où || . || désigne la norme euclidienne.
  
  

On résume dans les tableaux ?? et ?? les différentes contributions (hors convection-diffusion) affectées à chacun des vecteurs PROPCE, SMBRS et ROVSDT suivant le schéma en temps choisi pour les termes sources. En l’absence d’indication, les propriétés physiques ρ,µ, ... sont supposées prises en au temps n+θ_Φ, et le flux de masse ( ρ u) pris au temps n+θ_F, les valeurs de θ_F et de θ_Φ dépendant du type de schéma sélectionné spécifiquement pour ces grandeurs⁷.

Avec extrapolation des termes sources :

Sans extrapolation des termes sources :

0.4  Points à traiter

• Intégration du terme de convection-diffusion
Dans ce sous-programme, les points litigieux sont dus à l’intégration du terme de convection-diffusion. On renvoie donc le lecteur au sous-programme bilsc2 qui les explicite.

0.5   Annexe 1 : Inversibilité de la matrice EM _n

Dans cette section, on va étudier plus particulièrement l’inversibilité de la matrice EM _n, matrice du système linéaire à résoudre associée à EM_n pour le cas d’un schéma en temps de type Euler implicite d’ordre un (θ=1). Pour toutes les notations, on se reportera à la documentation sur le sous-programme covofi . On va montrer que la démarche adoptée permet de s’assurer que la matrice des systèmes de convection-diffusion dans les cas courants est toujours inversible.

0.5.1  Introduction

Pour montrer l’inversibilité, on va utiliser le fait que la dominance stricte de la diagonale l’implique⁸. On cherche donc à déterminer sous quelles conditions les matrices de convection diffusion sont à diagonale strictement dominante.

On va montrer qu’en incluant dans la matrice le terme en div(ρ  u) issu de ∂ ρ/∂ t, on peut établir directement et exactement⁹ la propriété. Par contre, si ce terme n’est pas pris en compte dans la matrice, il est nécessaire de faire intervenir le pendant discret de la relation : Cette relation n’est cependant vérifiée au niveau discret qu’à la précision du solveur de pression près (et, en tous les cas, ne peut être approchée au mieux qu’à la précision machine près). Il paraît donc préférable de s’en affranchir.

Avant d’entrer dans les détails de l’analyse, on rappelle quelques propriétés et définitions.

Soit C une matrice carrée d’ordre N, d’élément générique C_ij. On a par définition :

Définition : La matrice C est à diagonale strictement dominante ssi

On convient de dire que C est à diagonale simplement dominante ssi l’inégalité n’est pas stricte, soit :

Remarque : Si, sur chaque ligne, la somme des éléments d’une matrice est nulle, que les éléments extradiagonaux sont négatifs et que les éléments diagonaux sont positifs, alors la matrice est à diagonale simplement dominante. Si la somme est strictement positive, la diagonale est strictement dominante.

On a l’implication suivante :

Propriété : Si la matrice C est à diagonale strictement dominante, elle est inversible.
Cette propriété¹⁰ se démontre simplement si l’on admet le théorème de Gerschgörin ci-dessous :

Théorème : Soit B une matrice carrée d’ordre N dans lC × lC, d’élément générique B _ij, les valeurs propres λ_l de B sont, dans le plan complexe, telles que ||λ_l − B_ii|| _lC ∑_{j=1, j≠ i}^j=N||B_ij|| _lC

Si B est à éléments réels, on écrira ||λ_l − B_ii|| _lC ∑_{j=1, j≠ i}^j=N|B_ij|

Démontration de la propriété précédente :

Soit C à diagonale strictement dominante à éléments réels. On montre qu’il est possible d’inverser le système CX=S d’inconnue X, quel que soit le second membre S. Pour cela, on décompose C en partie diagonale (D) et extradiagonale (−E) soit :

C étant à diagonale strictement dominante, tous ses éléments diagonaux sont non nuls. D est donc inversible (et les élements de l’inverse sont réels). On considère alors la suite¹¹ :

On peut écrire :

Cette somme converge si le rayon spectral de B=D⁻¹E est strictement inférieur à 1. Or, la matrice C est à diagonale strictement dominante. On a donc pour tout i∈IN (à partir de la relation (??) et en divisant par |C_ii|) :

ce qui s’écrit encore :

ou bien :

d’où, avec le théorème de Gerschgörin, une relation sur les valeurs propres λ_l de B :

et comme B_ii=0 :

en particulier, la valeur propre dont la norme est la plus grande vérifie également cette équation. Ceci implique que le rayon spectral de B est strictement inférieur à 1. La suite (Xⁿ)_n∈IN converge donc (et la méthode de Jacobi converge). Il existe donc une solution à l’équation CX=S. Cette solution est unique¹² et la matrice C est donc inversible.

0.5.2  Avec prise en compte des termes issus de ∂ ρ/∂ t dans EM _n

Introduction

Pour montrer que la matrice EM _n est inversible, on va montrer qu’elle est à diagonale strictement dominante. Pour cela, on va considérer successivement les contributions :
   - des termes différentiels d’ordre 0 linéaires en δ f ^n+1,k+1,
   - des termes issus de la prise en compte de ∂ ρ/∂ t,
   - des termes différentiels d’ordre 1 (convection),
   - des termes différentiels d’ordre 2 (diffusion).

Pour chacune de ces contributions, on va examiner la dominance de la diagonale de l’opérateur linéaire associé. Si, pour chaque contribution, la dominance de la diagonale est acquise, on pourra alors conclure à la dominance de la diagonale pour la matrice (somme) complète¹³ EM _n et donc à son inversibilité.

Contributions des termes différentiels d’ordre 0 linéaires en δ f ^n+1,k+1

L’unique contribution est sur la diagonale : il faut donc vérifier qu’elle est strictement positive.

Pour chaque ligne I, f_s ^imp_I (cf. (??)) contient au minimum la quantité strictement positive¹⁴ ρ_Iⁿ |Ω_i|/Δ t. Les autres expressions, (− |Ω_i| (T_s ^imp)_I , +|Ω_i|  Γ_I , − |Ω_i| (T_s ^pd,imp) _I), lorsqu’elles existent, contribuent toutes positivement¹⁵.

L’opérateur linéaire associé à ces contributions vérifie donc bien la dominance stricte de la diagonale (propriété 1). Ce n’est cependant pas vrai si on extrapole les termes source, à cause de T_s ^imp. Il en résulte une contrainte sur la valeur du pas de temps.

Contributions des termes différentiels d’ordre 1 et des termes issus de la prise en compte de ∂ ρ/∂ t

Les termes considérés sont au nombre de deux dans (??) :
   - la contribution issue de la prise en compte de ∂ ρ/∂ t se retrouve dans f_s ^imp_I (équation ??),
   - la contribution du terme de convection linéarisé.

Après intégration spatiale, la somme de ces deux termes discrets s’écrit :

Pour chaque ligne I, on va chercher les propriétés de dominance de la diagonale en traitant séparément les faces internes (équation (??)) et les faces de bord (équation (??)).

 • la contribution des faces internes ij (facteur de δ f_I ^n+1,k+1) à la diagonale est positive ; la contribution aux extradiagonales est négative (facteur de δ f_J ^n+1,k+1) et la somme de ces contributions est exactement nulle (équation  (??)). Si l’on note C_IJ les coefficients de la matrice issus de la contribution de ces termes, on a donc |C_II| ∑_{J=1, J≠ I}^J=N|C_IJ| qui traduit la dominance “simple” (l’inégalité n’est pas “stricte”) de la diagonale et règle la question des contributions des faces internes.

 • la contribution des faces de bord doit être réécrite en utilisant l’expression des conditions aux limites sur f pour préciser la valeur de δ f _bik (on omet l’exposant ( n+1,k+1) pour alléger les notations) :
     - pour une condition de Dirichlet : δ f _bik = 0,
     - pour une condition de Neumann : δ f _bik = δ f_I,
     - pour une condition mixte (f _bik = α + β f_i) : δ f _bik = β δ f_I.
   Pour la contribution des faces de bord, il faut alors considérer deux cas de figure possibles.

• Le flux de masse au bord est positif ou nul ( m _bikⁿ = (ρ u) _bikⁿ . S _bik 0). Cette situation correspond par exemple aux sorties standards (fluide sortant du domaine), aux symétries ou aux parois étanches (flux de masse nul). Les contributions aux faces de bord sont alors toutes nulles, quelles que soient les conditions aux limites portant sur la variable f. En conséquence, la diagonale issue de ces contribution est simplement dominante.  
• Le flux de masse au bord est strictement négatif. Cette situation correspond à une entrée de fluide dans le domaine. Les contributions considérées s’écrivent alors : Il convient alors de distinguer plusieurs situations, selon le type de condition à la limite portant sur f :
  si la condition à la limite de f est de type Dirichlet, seule subsiste une contribution positive ou nulle à la diagonale, qui assure donc la dominance simple : si la condition à la limite de f est de type Neumann, la somme des contributions dues aux faces de bord est alors nulle, ce qui assure donc la dominance simple.
  si la condition à la limite de f est de type mixte, la contribution des faces de bord est sur la diagonale et vaut : On ne peut pas se prononcer quand à la dominance de la diagonale, à cause de la présence de (1−β) (la valeur de β est fixée par l’utilisateur) et la démarche adoptée ici ne permet donc pas de conclure. Il faut néanmoins noter que cette situation est rare dans les calculs standards. Elle demande un complément d’analyse et sera pour le moment exclue des considérations exposées dans le présent document.
  

On peut conclure, quand il n’y a pas de condition à la limite de type mixte, que la matrice associée aux contributions des termes différentiels d’ordre 1 (convectifs) et à la prise en compte des termes issus de ∂ ρ/∂ t et est à diagonale simplement dominante.

Contributions des termes différentiels d’ordre 2

On va considérer enfin les contributions des termes différentiels d’ordre 2 (issus du terme
− div (Kⁿ grad δ f ^n+1,k+1)). Pour ces termes, la contribution à la diagonale est positive¹⁶, négative aux extradiagonales¹⁶, compte tenu de :

Considérons deux cas :
   - la cellule courante I n’a que des faces internes au domaine de calcul (pas de faces de bord). La somme des contributions est nulle. On a donc dominance simple de la diagonale.
   - la cellule courante I a des faces de bord. La somme des contributions diagonales et extradiagonales est positive quand on a une condition à la limite de type Dirichlet ou de type Neumann sur f. La diagonale est alors strictement dominante. Lorsqu’il y a des conditions à la limite de type mixte, il n’est plus possible de conclure (situation écartée précédemment).
On peut conclure, quand il n’y a pas de condition à la limite de type mixte, que la matrice associée aux contributions des termes différentiels d’ordre 2 est au moins à diagonale simplement dominante.

Conclusion

En travaillant sur des maillages non pathologiques (à transmittivité positive, voir la note de bas de page numéro 0.5.2) et en n’imposant pas de condition à la limite de type mixte sur les variables, on peut donc conclure que EM _n est la somme de matrices à diagonale simplement dominante et d’une matrice à diagonale strictement dominante (paragraphe 0.5.2). Elle est donc à diagonale strictement dominante, et donc inversible (de plus, la méthode itérative de Jacobi converge).

0.5.3  Sans prise en compte des termes issus de ∂ ρ/∂ t dans EM _n

Introduction

Pour identifier les cas dans lesquels la matrice EM _n est inversible, on va rechercher les conditions qui assurent la dominance de la diagonale. Par rapport à l’analyse présentée au paragraphe 0.5.2, seules diffèrent les considérations relatives aux contributions des termes différentiels d’ordre 1, puisqu’elles sont traitées au paragraphe 0.5.2 avec les termes issus de la prise en compte de ∂ ρ/∂ t.

Contributions des termes différentiels d’ordre 1

La contribution du terme de convection est la seule à prendre en compte. Elle s’écrit, d’après les équations (??) et la discrétisation explicitée pour le sous-programme covofi :

On constate que pour chaque ligne I, la contribution des faces internes (facteur de δ f_I ^n+1,k+1) à la diagonale est positive et qu’elle est négative aux extradiagonales (facteur de δ f_J ^n+1,k+1). Cependant, contrairement au cas présenté au paragraphe 0.5.2, la somme de ces contributions n’est pas nulle dans le cas général. Pour obtenir un résultat quant à la dominance de la diagonale, il faut faire intervenir la version discrète de la propriété (??) rappelée ci-dessous :

Soit, sous forme discrète :

Il n’est donc pas possible d’analyser séparément les contributions des faces internes et celles des faces de bord (contrairement à la situation rencontrée au paragraphe 0.5.2). On se place ci-après dans le cas général d’une cellule qui a des faces internes et des faces de bord (si elle n’a que des faces internes, la démonstration est la même, mais plus simple. On peut l’écrire en considérant formellement que la cellule “a zéro faces de bord”, c’est à dire que γ_b(i) est l’ensemble vide).
Il faut alors considérer deux cas de figure, selon la valeur du flux de masse aux faces de bord (éventuelles) de la cellule :

• Le flux de masse au bord est positif ou nul ( m _bikⁿ = (ρ u) _bikⁿ . S _bik 0). Cette situation correspond à des cellules qui ont des faces de bord de sortie standard (fluide sortant du domaine), de symétrie ou de paroi étanche (flux de masse nul). Les contributions s’écrivent alors : Dans ce cas, la somme des contributions à la diagonale est positive, les contributions aux extradiagonales sont négatives et, avec la relation (??), on vérifie que la somme des contributions diagonales et extradiagonales est nulle. On a donc dominance simple de la diagonale.
• Le flux de masse au bord est strictement négatif. Cette situation correspond à des cellules qui ont des faces de bord d’entrée standard (entrée de fluide dans le domaine). Les contributions considérées s’écrivent alors : Il convient alors de distinguer plusieurs situations, selon le type de condition à la limite portant sur f (on omet l’exposant ( n+1,k+1) pour alléger les notations) :
  - pour une condition de Dirichlet : δ f _bik = 0,
  - pour une condition de Neumann : δ f _bik = δ f_I,
  - pour une condition mixte (f _bik = α + β f_I) : δ f _bik = β δ f_I.
  Selon le cas on se trouve dans une des situations suivantes :
  si la condition à la limite de f est de type Dirichlet, la contribution des faces de bord est nulle dans la matrice. La contribution des faces internes à la diagonale est positive, la contribution aux extradiagonales négative et la somme de ces contributions vaut ∑_{j∈ Vois(i)} m _ijⁿ, soit avec la relation (??) :

  ∑ j∈ Vois(i)

  m ijn=−

  ∑ k∈ γb(i)

  m bikn

  Elle est strictement positive et la diagonale est donc strictement dominante.
  si la condition à la limite de f est de type Neumann, la contribution des faces de bord se réduit au terme : ∑_{k∈ γb(i)} m _bikⁿ  δ f_I ^n+1,k+1 . La somme des contributions à la diagonale est alors SC_i:

  SCi=

  1 2

  ∑ j∈ Vois(i)

  ( + m ijn + | m ijn| ) +

  ∑ k∈ γb(i)

  m bikn

  En utilisant deux fois la relation (??), on obtient donc pour la diagonale :

  SCi=

  1 2

  ⎡ ⎢ ⎢ ⎣

  ∑ j∈ Vois(i)

  | m ijn| +

  ∑ k∈ γb(i)

  m bikn

  ⎤ ⎥ ⎥ ⎦

  =

  1 2

  ⎡ ⎢ ⎢ ⎣

  ∑ j∈ Vois(i)

  ( | m ijn|− m ijn )

  ⎤ ⎥ ⎥ ⎦

  Cette grandeur est positive et égale à l’opposé de la somme des termes extradiagonaux qui sont tous négatifs. La diagonale est donc simplement dominante.
  si la condition à la limite de f est de type mixte, la somme des contributions dues aux faces de bord est :
  On ne peut donc pas conclure quant au signe de cette contribution, le facteur β étant choisi librement par l’utilisateur. Cette situation a été écartée dans le paragraphe 0.5.2.
  

On peut donc conclure, quand il n’y a pas de condition à la limite de type mixte, que la matrice associée aux contributions des termes différentiels d’ordre 1 (convectifs) est à diagonale simplement dominante, à condition que la relation (??) soit vérifiée exactement.

Conclusion

En travaillant sur des maillages non pathologiques (à transmittivité positive, voir la note de bas de page numéro 0.5.2) et en n’imposant pas de condition à la limite de type mixte sur les variables, on peut donc conclure que EM _n est à diagonale strictement dominante, donc inversible (et la méthode itérative de Jacobi converge) à condition que la relation (??) soit vérifiée exactement. Ce n’est généralement pas le cas (la précision du solveur de pression et la précision machine sont finies). Même si la contribution diagonale en ρ_Iⁿ |Ω_i|/Δ t peut suffire à assurer la dominance, on a cependant souhaité, dans Code Saturne, s’affranchir du problème potentiel en prenant en compte les termes issus de ∂ ρ/∂ t dans la matrice.

0.6   Annexe 2 : Remarques à propos du respect du principe du maximum discret

0.6.1  Introduction

Les considérations exposées ici sont relatives au fait que, en continu, une variable qui n’est que convectée par un champ de débit à divergence nulle doit rester dans les bornes minimales et maximales définies par les conditions initiales et par les conditions aux limites en espace. Ainsi, les valeurs d’un scalaire passif initialement nul dont les conditions aux limites sont des conditions de Neumann homogène et des conditions de Dirichlet de valeur 1 devront nécessairement rester comprises dans l’intervalle [0 ; 1]. C’est ce que l’on entend ici par “principe du maximum”.

Soient u un champ de vitesse figé et connu et t un réel positif. On considère le problème modèle P de convection des variables scalaires ρ et ρ f, défini par : avec une condition initiale f⁰ donnée ainsi que des conditions aux limites associées sur f de type Dirichlet ou Neumann.
Dans Code Saturne, la deuxième équation de P est réécrite en continu, en utilisant la première, sous la forme :
et discrétisée temporellement comme suit : Dans un premier temps, on va étudier la discrétisation spatiale associée à (??), qui correspond donc à la prise en compte de la contribution de ∂ ρ/∂ t dans l’équation en continu (et se traduit par la présence de −div((ρ u)ⁿ) dans l’expression de f_s ^imp_I ), puis dans un deuxième temps, la discrétisation spatiale de l’expression ;
qui correspond à un problème de convection pure classique.

On étudiera ensuite un exemple simplifié (monodimensionnel à masse volumique constante).

Les considérations présentes mériteraient d’être approfondies.

0.6.2  Cas général

Discrétisation spatiale de ρⁿ   fⁿ⁺¹ − fⁿ/Δ t − fⁿ⁺¹ div (ρ u)ⁿ +  div ((ρ u)ⁿ fⁿ⁺¹) = 0

En intégrant sur une cellule Ω_i à l’aide de la formulation volumes finis habituelle, on obtient : où f_f ij ⁿ⁺¹ et f_f b ik ⁿ⁺¹ sont les valeurs de f aux faces internes et de bord issues du choix du schéma convectif.

En reprenant les notations précédentes, en imposant un schéma décentré amont au premier membre (i.e. en exprimant δ f _ij ^n+1,k+1 et δ f _bik ^n+1,k+1) et en raisonnant en incréments (cf. sous-programme navsto ), on aboutit à : avec : et (f ^n+1,k)_{k∈ IN} suite convergeant vers f ⁿ⁺¹, n entier donné, solution de (??) .

Discrétisation spatiale de ρⁿ   fⁿ⁺¹ − fⁿ/Δ t + div ((ρ u)ⁿ fⁿ⁺¹) = 0

En procédant de façon analogue et en adoptant les mêmes hypothèses, on obtient : (où f_f ij ⁿ⁺¹ et f_f b ik ⁿ⁺¹ sont les valeurs de f aux faces internes et de bord issues du choix du schéma convectif)

avec : et (f ^n+1,k)_{k∈ IN} suite convergeant vers f ⁿ⁺¹, n entier donné, solution de (??) .

0.6.3  Exemple pour le principe du maximum

On va maintenant se placer en monodimensionnel, sur un maillage régulier formé de trois cellules de pas h constant (figure II.0.1) et étudier le comportement du premier membre pour les deux types d’expressions, entre le pas de temps n Δ t et le pas de temps (n+1) Δ t, avec, comme condition initiale f₁⁰ = f₂⁰ =f₃⁰ = 0 et comme conditions aux limites, une de type Dirichlet et l’autre de type Neumann homogène :

On supposera de plus que :
    le schéma convectif utilisé est le schéma upwind
    la masse volumique est constante
    (ρ u) _b1ⁿ > 0, (ρ u) ₁₂ⁿ > 0 , (ρ u) ₂₃ⁿ >  0 , (ρ u) _b2ⁿ > 0 et S _b1 < 0.

---

[width=3.5cm,angle=-90]./graphics/domaine1d

Figure II.0.1: Définition du domaine de calcul unidimensionnel considéré.

---

On s’intéresse à l’influence sur le respect du principe du maximum discret de la précision avec laquelle est vérifiée sous forme discrète la relation :

soit, ici :

Prise en compte de la contribution de ∂ ρ/∂ t dans la matrice

Le système à résoudre est alors, en omettant pour simplifier l’exposant ( n+1,k+1) : ce qui donne : d’où : On obtient donc bien une solution qui vérifie le principe du maximum discret, même pour des grands pas de temps Δ t, et ce, quelle que soit la précision avec laquelle est vérifiée, à l’étape de correction, la forme discrète (??) de la conservation de la masse ∫_Ωidiv (ρ u) dΩ = 0 dont on ne s’est pas servi ici.

Sans la contribution de ∂ ρ/∂ t dans la matrice

On obtient de même :

soit :

Ici, on constate que le respect du principe du maximum discret :
est équivalent à la condition : Contrairement à la situation du paragraphe 0.6.3, on ne peut obtenir ici un résultat qu’en faisant intervenir l’égalité (??), forme discrète de la conservation de la masse. On obtient bien alors, à partir du système précédent :

Si l’on s’intéresse à la cellule Ω₁ et que l’on suppose (ρ u) ₁₂ⁿ . S ₁₂=− (ρ u) _b1ⁿ . S _b1−ε (ρ u) ₁₂ⁿ . S ₁₂ (où ε est la précision locale relative pour l’équation de conservation de la masse discrète), on constate que l’on obtient δ f₁ > f_b1 = 1 (valeur non admissible) dès lors que :

c’est-à-dire dès que le nombre de CFL local (ρ u) ₁₂ⁿ . S ₁₂Δ t/ρ₁|Ω₁| excède l’inverse de la précision relative locale ε.

0.6.4  Conclusion

Prendre en compte la contribution de ∂ ρ/∂ t dans la matrice permet un meilleur respect du principe du maximum discret, lorsque la précision de ∫_Ωidiv (ρ u) dΩ = 0 n’est pas exactement vérifiée.

1

Davroux A. et Archambeau F. : Calcul de la variance des fluctuations d’un scalaire dans le solveur commun. Application à l’expérience du CEGB dite “Jet in Pool”, HE-41/99/043.

2

a et a"², sous forme discrète en espace, correspondent donc en fait à des vecteurs dimensionnés à NCELET de composantes a_I et a"²_I respectivement, I décrivant l’ensemble des cellules.

3

Si θ=1/2, ou qu’une extrapolation est utilisée, le pas de temps Δ t est supposé uniforme en temps et en espace.

4

cf. introd

5

cf. preduv

6

cette opération est faite quel que soit le schéma en temps de façon à rester cohérent avec ce qui est fait dans bilsc2

7

cf. introd

8

Ce faisant, on choisit cependant une condition forte et la démonstration n’est probablement pas optimale.

9

Hormis dans le cas de conditions aux limites mixtes, qu’il conviendrait d’examiner plus en détail.

10

Lascaux, P. et Théodor, R. : Analyse Numérique Matricielle Appliquée à l’art de l’Ingénieur, Tome 2, Ed. Masson.

11

On reconnaît la méthode de Jacobi

12

On peut le voir “par l’absurde”. En effet, supposons qu’il existe deux solutions distinctes X₁ et X₂ à l’équation CX=S. Alors Y=X₂−X₁ vérifie CY=0, soit DY=−EY, donc D⁻¹EY=−Y. Ceci signifie que Y (qui n’est pas nul, par hypothèse) est vecteur propre de D⁻¹E avec λ=−1 pour valeur propre associée. Or, le rayon spectral de D⁻¹E est strictement inférieur à 1 et λ=−1 ne peut donc pas être une valeur propre de D⁻¹E. En conséquence, il ne peut exister qu’une seule solution à l’équation CX=S.

13

Ce raisonnement n’est pas optimal (la somme de valeurs absolues étant supérieure à la valeur absolue de la somme), mais permet d’obtenir des conclusions dans le cas présent (condition suffisante).

14

Ceci permettra de conclure à la stricte dominance de la diagonale de la matrice somme complète EM _n.

15

Le terme de dissipation ρ1/R_f ε/k, spécifique à l’étude de la variance des fluctuations, est positif par définition et ne remet donc pas en cause la conclusion.

16

Ceci n’est en fait pas toujours vrai. En effet, pour chaque face ij, la transmittivité Kⁿ/I′J′S_ij fait intervenir la mesure algébrique du segment I′J′, où I′ et J′ sont les projetés orthogonaux sur la normale à la face du centre des cellules voisines. Cette grandeur est une valeur algébrique et peut théoriquement devenir négative sur certains maillages pathologiques, contenant par exemple des mailles non convexes. On pourra se reporter au dernier point à traiter du sous-programme matrix .