Chapter 4  codits

• Fonction
• Discrétisation
• Mise en œuvre
• Points à traiter

0.1  Fonction

Ce sous-programme, appelé entre autre par preduv , turbke , covofi , resrij , reseps , ..., résout les équations de convection-diffusion d’un scalaire a avec termes sources du type : où ρ u, f_s^exp et f_s^imp désignent respectivement le flux de masse, les termes sources explicites et les termes linéarisés en aⁿ⁺¹. a est un scalaire défini sur toutes les cellules¹. Par souci de clarté on suppose, en l’absence d’indication, les propriétes physiques Φ (viscosité totale µ_tot,...) et le flux de masse (ρ u) 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².
L’écriture des termes de convection et diffusion en maillage non orthogonal engendre des difficultés (termes de reconstruction et test de pente) qui sont contournées en utilisant une méthode itérative dont la limite, si elle existe, est la solution de l’équation précédente.

0.2  Discrétisation

Afin d’expliquer la procédure utilisée pour traiter les difficultés dues aux termes de reconstruction et de test de pente dans les termes de convection-diffusion, on note, de façon analogue à ce qui est défini dans navsto mais sans discrétisation spatiale associée, E_n l’opérateur : L’équation (??) s’écrit donc : La quantité E_n(aⁿ⁺¹) comprend donc :
f_s ^imp aⁿ⁺¹, contribution des termes différentiels d’ordre 0 linéaire en aⁿ⁺¹,
θ   div((ρu) aⁿ⁺¹) − θ   div(µ _tot grad aⁿ⁺¹), termes de convection-diffusion implicites complets (termes non reconstruits + termes de reconstruction) linéaires³ en aⁿ⁺¹,
f_s ^exp− f_s ^imp aⁿ et (1−θ)  div((ρ u) aⁿ) − (1−θ)   div(µ _tot grad aⁿ) l’ensemble des termes explicites (y compris la partie explicite provenant du schéma en temps appliqué à la convection diffusion).

De même, on introduit un opérateur EM_n approché de E_n, linéaire et simplement inversible, tel que son expression contient :
la prise en compte des termes linéaires en a,
la convection intégrée par un schéma décentré amont (upwind) du premier ordre en espace,
les flux diffusifs non reconstruits.
Cet opérateur permet donc de contourner la difficulté induite par la présence d’éventuelles non linéarités introduites par l’activation du test de pente lors du schéma convectif, et par le remplissage important de la structure de la matrice découlant de la présence des gradients propres à la reconstruction.
On a la relation⁴, pour toute cellule Ω_I de centre I : On cherche à résoudre : Soit : On va pour cela utiliser un algorithme de type point fixe en définissant la suite (a^n+1, k)_{k∈ IN}⁵: où δ a^n+1, k+1 est solution de : Soit encore, par linéarité de EM_n :

Cette suite, couplée avec le choix de l’opérateur E_n, permet donc de lever la difficulté induite par la présence de la convection (discrétisée à l’aide de schémas numériques qui peuvent introduire des non linéarités) et les termes de reconstruction. Le schéma réellement choisi par l’utilisateur pour la convection (donc éventuellement non linéaire si le test de pente est activé) ainsi que les termes de reconstruction vont être pris à l’itération k et traités au second membre via le sous-programme bilsc2 , alors que les termes non reconstruits sont pris à l’itération k+1 et représentent donc les inconnues du système linéaire résolu par codits ⁶.
On suppose de plus que cette suite (a^n+1, k)_k converge vers la solution aⁿ⁺¹ de l’équation (??), i.e. lim_k→∞ δ a^n+1, k = 0, ceci pour tout n donné.
(??) correspond à l’équation résolue par codits . La matrice EM _n, matrice associée à EM_n est à inverser, les termes non linéaires sont mis au second membre mais sous forme explicite (indice k de a^n+1, k) et ne posent donc plus de problème.

Remarque 1

La viscosité µ _tot prise dans EM_n et dans E_n dépend du modèle de turbulence utilisé. Ainsi on a µ _tot=µ _laminaire + µ _turbulent dans EM_n et dans E_n sauf lorsque l’on utilise un modèle R_ij−ε, auquel cas on a µ _tot=µ _laminaire.
Le choix de EM_n étant a priori arbitraire (EM_n doit être linéaire et la suite (a^n+1, k)_k∈IN doit converger pour tout n donné), une option des modèles R_ij−ε (IRIJNU=1) consiste à forçer µ _totⁿ dans l’expression de EM_n à la valeur µ _laminaireⁿ + µ _turbulentⁿ lors de l’appel à codits dans le sous-programme navsto , pour l’étape de prédiction de la vitesse. Ceci n’a pas de sens physique (seul µ _laminaireⁿ étant censé intervenir), mais cela peut dans certains cas avoir un effet stabilisateur, sans que cela modifie pour autant les valeurs de la limite de la suite (a^n+1, k)_k.

Remarque 2

Quand codits est utilisé pour le couplage instationnaire renforcé vitesse-pression (IPUCOU=1), on fait une seule itération k en initialisant la suite (a^n+1, k)_k∈IN à zéro. Les conditions de type Dirichlet sont annulées (on a INC = 0) et le second membre est égal à ρ |Ω_i|. Ce qui permet d’obtenir une approximation de type diagonal de EM_n nécessaire lors de l’étape de correction de la vitesse⁷.

0.3  Mise en œuvre

L’algorithme de ce sous-programme est le suivant :
- détermination des propriétés de la matrice EM_n (symétrique si pas de convection, non symétrique sinon)
- choix automatique de la méthode de résolution pour l’inverser si l’utilisateur ne l’a pas spécifié pour la variable traitée. La méthode de Jacobi est utilisée par défaut pour toute variable scalaire a convectée. Les méthodes disponibles sont la méthode du gradient conjugué, celle de Jacobi, et le bi-gradient conjugué stabilisé (BICGStab) pour les matrices non symétriques. Un préconditionnement diagonal est possible et utilisé par défaut pour tous ces solveurs excepté Jacobi.
- prise en compte de la périodicité (translation ou rotation d’un scalaire, vecteur ou tenseur),
- construction de la matrice EM_n correspondant à l’opérateur linéaire EM_n par appel au sous-programme matrix ⁸. Les termes implicites correspondant à la partie diagonale de la matrice et donc aux contributions différentielles d’ordre 0 linéaires en aⁿ⁺¹,(i.e f_s^imp), sont stockés dans le tableau ROVSDT (réalisé en amont du sous-programme appelant codits ).
- création de la hiérarchie de maillage si on utilise le multigrille ( IMGRP >0 ).
- appel à bilsc2 pour une éventuelle prise en compte de la convection-diffusion explicite lorsque θ ≠ 0.
- boucle sur le nombre d’itérations de 1 à NSWRSM (appelé NSWRSP dans codits ). Les itérations sont représentées par k appelé ISWEEP dans le code et définissent les indices de la suite (a^n+1, k)_k et de (δ a^n+1, k)_k.
Le second membre est scindé en deux parties :
    un terme, affine en a^n+1, k−1, facile à mettre à jour dans le cadre de la résolution par incrément, et qui s’écrit :
    les termes issus de la convection/diffusion (avec reconstruction) calculée par bilsc2 .

La boucle en k est alors la suivante :

• Calcul du second membre, hors contribution des termes de convection-diffusion explicite SMBINI; le second membre complet correspondant à E_n(a^n+1, k−1) est, quant à lui, stocké dans le tableau SMBRP, initialisé par SMBINI et complété par les termes reconstruits de convection-diffusion par appel au sous-programme bilsc2 .
  À l’itération k, SMBINI noté SMBINI ^k vaut :
  
  • Avant de boucler sur k,un premier appel au sous-programme bilsc2 avec THETAP=1−θ permet de prendre en compte la partie explicite des termes de convection-diffusion provenant du schéma en temps. Avant de boucler sur k, le second membre SMBRP ⁰ est stocké dans le tableau SMBINI ⁰ et sert pour l’initialisation du reste du calcul.
  • pour k = 1, On a donc : et SMBRP ¹ est complété par un second appel au sous-programme bilsc2 avec THETAP=θ, de manière à ajouter dans le second membre la partie de la convection-diffusion implicite. • pour k = 2,
  de façon analogue, on obtient : Soit : l’appel au sous-programme bilsc2 , étant systématiquement fait par la suite avec THETAP=θ, on obtient de même : où • pour l’itération k+1,
  Le tableau SMBINI ^k+1 initialise le second membre complet SMBRP ^k+1 auquel vont être rajoutées les contributions convectives et diffusives via le sous-programme bilsc2 .
  on a la formule : Puis suit le calcul et l’ajout des termes de convection-diffusion reconstruits de − E_n(a^n+1, k), par appel au sous-programme bilsc2 . On rappelle que la convection est prise en compte à cette étape par le schéma numérique choisi par l’utilisateur (schéma décentré amont du premier ordre en espace, schéma centré du second ordre en espace, schéma décentré amont du second ordre S.O.L.U. ou une pondération (blending) des schémas dits du second ordre (centré ou S.O.L.U.) avec le schéma amont du premier ordre, utilisation éventuelle d’un test de pente).
  Cette contribution (convection-diffusion) est alors ajoutée dans le second membre SMBRP ^k+1 (initialisé par SMBINI ^k+1).
• Résolution du système linéaire en δ a^n+1, k+1 correspondant à l’équation (??) par inversion de la matrice EM_n, en appelant le sous programme invers . On calcule a^n+1, k+1 grâce à la formule : Soit :
• Traitement de la périodicité et du parallélisme.
• Test de convergence :
  Il porte sur la quantité ||SMBRP ^k+1|| < ε ||EM_n(aⁿ) + SMBRP ¹|| , où || . || représente la norme euclidienne. Si le test est vérifié, la convergence est atteinte et on sort de la boucle sur les itérations. La solution recherchée est a ⁿ⁺¹ = a^n+1, k+1.
  Sinon, on continue d’itérer dans la limite des itérations imposées par NSWRSM dans usini1 .
  En “continu” ce test de convergence s’écrit aussi : Si bien que sur maillage orthogonal avec schéma de convection upwind et en l’absence de terme source, la suite converge en théorie en une unique itération puisque par construction :

Fin de la boucle.

0.4  Points à traiter

• Approximation EM_n de l’opérateur E_n
D’autres approches visant soit à modifier la définition de l’approximée, prise en compte du schéma centré sans reconstruction par exemple, soit à abandonner cette voie seraient à étudier.

• Test de convergence
La quantité définissant le test de convergence est également à revoir, éventuellement à simplifier.

• Prise en compte de T_s^imp
Lors de la résolution de l’équation par codits , le tableau ROVSDT a deux fonctions : il sert à calculer la diagonale de la matrice (par appel de matrix ) et il sert à mettre à jour le second membre à chaque sous-itération de la résolution en incréments. Or, dans le cas où T_s^imp est positif, on ne l’intègre pas dans ROVSDT, afin de ne pas affaiblir la diagonale de la matrice. De ce fait, on ne l’utilise pas pour mettre à jour le second membre, alors que ce serait tout à fait possible. Au final, on obtient donc un terme source traité totalement en explicite (T_s^exp+T_s^impaⁿ), alors que la résolution en incréments nous permettrait justement de l’impliciter quasiment totalement (T_s^exp+T_s^impa^{n+1,k_fin−1}, où k_fin est la dernière sous-itération effectuée).
Pour ce faire, il faudrait définir deux tableaux ROVSDT dans codits .

1

a, sous forme discrète en espace, correspond à un vecteur dimensionné à NCELET de composante a_I, I décrivant l’ensemble des cellules.

2

cf. introd

3

Lors de la discrétisation en espace, le caractère linéaire de ces termes pourra cependant être perdu, notamment à cause du test de pente.

4

On pourra se reporter au sous-programme matrix pour plus de détails relativement à EM_disc, opérateur discret agissant sur un scalaire a.

5

Dans le cas ou le point fixe en vitesse-pression est utilisé (NTERUP> 1) a^n+1,0 est initialisé par la dernière valeur obtenue de aⁿ⁺¹.

6

cf. le sous-programme navsto .

7

cf. le sous-programme resolp .

8

On rappelle que dans matrix , la convection est traitée, quelque soit le choix de l’utilisateur, avec un schéma décentré amont d’ordre 1 en espace et qu’il n’y a pas de reconstruction pour la diffusion. Le choix de l’utilisateur quant au schéma numérique pour la convection intervient uniquement lors de l’intégration des termes de convection de E_n, au second membre de (??) dans le sous-programme bilsc2 .