Chapter 16  turbke

• Fonction
• Discrétisation
• Mise en œuvre

0.1  Fonction

Le but de ce sous-programme est de résoudre le système des équations de k et ε de manière semi-couplée.
Le système d’équations résolu est le suivant :

P est le terme de production par cisaillement moyen :

G est le terme de production par gravité : G=−1/ρµ_t/σ_t ∂ρ/∂ x_ig_i

La viscosité turbulente est µ_t=ρ C_µk²/ε.

Les constantes sont :
C_µ=0,09 ; C_ε2=1,92 ; σ_k=1 ; σ_ε=1,3
C_ε3=0 si G0 (stratification instable) et C_ε3=1 si G0 (stratification stable).

Γ est un éventuel terme source de masse (tel que l’équation de conservation de masse devienne ∂ ρ/∂ t+div(ρu)=Γ). ϕ_i (ϕ=k ou ε) est la valeur de ϕ associée à la masse injectée ou retirée. Dans le cas où on retire de la masse (Γ<0), on a forcément ϕ_i=ϕ. De même, quand on injecte de la masse, on spécifie souvent aussi ϕ_i=ϕ. Dans ces deux cas, le terme disparaît de l’équation. Dans la suite du document, on qualifiera d’injection forcée les cas où on a Γ>0 et ϕ_i≠ϕ.

α_k, β_k, α_ε, β_ε sont des termes sources utilisateur éventuels, conduisant à une implicitation partielle, imposés le cas échéant par le sous-programme ustske .

0.2  Discrétisation

La résolution se fait en trois étapes, afin de coupler partiellement les deux variables k et ε. Pour simplifier, réécrivons le système de la façon suivante :

D est l’opérateur de convection/diffusion. S_k (resp. S_ε) est le terme source de k (resp. ε).

Première phase : bilan explicite

On résout le bilan explicite :

(le terme en Γ n’est pris en compte que dans le cas de l’injection forcée)

Deuxième phase : couplage des termes sources

On implicite les termes sources de manière couplée : soit

Le terme en div(ρu) n’est pas implicité car il est lié au terme D pour assurer que la matrice d’implicitation sera à diagonale dominante. Le terme en Γ et les termes sources utilisateur ne sont pas implicités non plus, mais ils le seront dans la troisième phase.

Et on écrit (pour ϕ=k ou ε)

On résout donc finalement le système 2× 2 :

Troisième phase : implicitation de la convection/diffusion

On résout le système : soit

Le terme en Γ n’est là encore pris en compte que dans le cas de l’injection forcée. Le terme en α n’est pris en compte que si α est négatif, pour éviter d’affaiblir la diagonale de la matrice qu’on va inverser.

Précisions sur le couplage

Lors de la phase de couplage, afin de privilégier la stabilité et la réalisabilité du résultat, tous les termes ne sont pas pris en compte. Plus précisément, on peut écrire :

en notant P = (∂ u_i/∂ x_j + ∂ u_j/∂ x_i)∂ u_i/∂ x_j −2/3(divu)²
et G = −1/ρσ_t ∂ρ/∂ x_ig_i

On a donc en théorie

En pratique, on va chercher à assurer k_ts>0 et ε_ts>0. En se basant sur un calcul simplifié, ainsi que sur le retour d’expérience d’ESTET, on montre qu’il est préférable de ne pas prendre en compte certains termes. Au final, on réalise le couplage suivant :

avec

(par définition de C_ε3, P+(1−C_ε3)G est toujours positif)

0.3  Mise en œuvre

• Calcul du terme de production
On appelle trois fois grdcel pour calculer les gradients de u, v et w. Au final, on a
TINSTK= 2(∂ u/∂ x)²+ 2(∂ v/∂ y)²+ 2(∂ w/∂ z)²+ (∂ u/∂ y+∂ v/∂ x)²+ (∂ u/∂ z+∂ w/∂ x)²+ (∂ v/∂ z+∂ w/∂ y)²
et
DIVU= ∂ u/∂ x+∂ v/∂ y +∂ w/∂ z

(le terme div(u) n’est pas calculé par divmas , pour correspondre exactement à la trace du tenseur des déformations qui est calculé pour la production)

• Lecture des termes sources utilisateur
Appel de ustske pour charger les termes sources utilisateurs. Ils sont stockés dans les tableaux suivants :
W7=Ωβ_k
W8=Ωβ_ε
DAM=Ωα_k
W9=Ωα_ε

Puis on ajoute le terme en (divu)² à TINSTK. On a donc
TINSTK=P

• Calcul du terme de gravité
Calcul uniquement si IGRAKE=1.
On appelle grdcel pour ROM, avec comme conditions aux limites COEFA=ROMB et COEFB=VISCB=0.
PRDTUR=σ_t est mis à 1 si on n’a pas de scalaire température.

G est calculé et les termes sources sont mis à jour :
TINSTK=P+G
TINSTE=P+Max[G,0] =P+(1−C_ε3)G

Si IGRAKE=0, on a simplement
TINSTK=TINSTE=P

• Calcul du terme d’accumulation de masse
On stocke W1=Ωdiv(ρu) calculé par divmas (pour correspondre aux termes de convection de la matrice).

• Calcul des termes sources explicites
On affecte les termes sources explicites de k et ε pour la première étape.
SMBRK=Ω(µ_t(P+G) −2/3ρ⁽ⁿ⁾ k⁽ⁿ⁾divu −ρ⁽ⁿ⁾ ε⁽ⁿ⁾)+Ω k⁽ⁿ⁾div(ρu)
SMBRE=Ωε⁽ⁿ⁾/k⁽ⁿ⁾ (C_ε1( µ_t(P+(1−C_ε3)G) −2/3ρ⁽ⁿ⁾ k⁽ⁿ⁾divu) −C_ε2ρ⁽ⁿ⁾ε⁽ⁿ⁾) +Ωε⁽ⁿ⁾div(ρu)

soit SMBRK=Ω S_k⁽ⁿ⁾+Ω k⁽ⁿ⁾div(ρu) et SMBRE=Ω S_ε⁽ⁿ⁾+Ωε⁽ⁿ⁾div(ρu).

• Calcul des termes sources utilisateur
On ajoute les termes sources utilisateur explicites à SMBRK et SMBRE, soit :
SMBRK=Ω S_k⁽ⁿ⁾+Ω k⁽ⁿ⁾div(ρu)+Ωα_k k⁽ⁿ⁾ +Ωβ_k
SMBRE=Ω S_ε⁽ⁿ⁾+Ωε⁽ⁿ⁾div(ρu) +Ωα_εε⁽ⁿ⁾ +Ωβ_ε

Les tableaux W7 et W8 sont libérés, DAM et W9 sont conservés pour être utilisés dans la troisième phase de résolution.

• Calcul des termes de convection/diffusion explicites
bilsc2 est appelé deux fois, pour k et pour ε, afin d’ajouter à SMBRK et SMBRE les termes de convection/diffusion explicites D(k⁽ⁿ⁾) et D(ε⁽ⁿ⁾). Ces termes sont d’abord stockés dans W7 et W8, pour être conservés et réutilisés dans la troisième phase de résolution.

• Termes source de masse
Dans le cas d’une injection forcée de matière, on passe deux fois dans catsma pour ajouter les termes en Ω Γ (k_i−k⁽ⁿ⁾) et Ω Γ (ε_i−ε⁽ⁿ⁾) à SMBRK et SMBRE. La partie implicite (ΩΓ) est stockée dans les variables W2 et W3, qui seront utilisées lors de la troisième phase (les deux variables sont bien nécessaires, au cas où on aurait une injection forcée sur k et pas sur ε, par exemple).

• Fin de la première phase
Ceci termine la première phase. On a
SMBRK=Ω ρ⁽ⁿ⁾k_e−k⁽ⁿ⁾/Δ t
SMBRE=Ω ρ⁽ⁿ⁾ε_e−ε⁽ⁿ⁾/Δ t

• Phase de couplage
(uniquement si IKECOU=1)

On renormalise SMBRK et SMBRE qui deviennent les seconds membres du système de couplage.
SMBRK=1/Ωρ⁽ⁿ⁾SMBRK =k_e−k⁽ⁿ⁾/Δ t
SMBRE=1/Ωρ⁽ⁿ⁾SMBRE =ε_e−ε⁽ⁿ⁾/Δ t
et DIVP23=2/3Max[div(u),0].

On remplit la matrice de couplage
A11=1/Δ t −2 C_µk⁽ⁿ⁾/ε⁽ⁿ⁾ Min[(P+G),0] +2/3Max[div(u),0]
A12= 1
A21= − C_ε1 C_µ(P +(1−C_ε3)G) − C_ε2(ε⁽ⁿ⁾/k⁽ⁿ⁾)²
A22=1/Δ t +2/3C_ε1Max[div(u),0] +2 C_ε2ε⁽ⁿ⁾/k⁽ⁿ⁾

On inverse le système 2× 2, et on récupère :
DELTK=k_ts−k⁽ⁿ⁾
DELTE=ε_ts−ε⁽ⁿ⁾

• Fin de la deuxième phase
On met à jour les variables SMBRK et SMBRE.
SMBRK=Ω ρ⁽ⁿ⁾k_ts−k⁽ⁿ⁾/Δ t
SMBRE= Ω ρ⁽ⁿ⁾ε_ts−ε⁽ⁿ⁾/Δ t

Dans le cas où on ne couple pas (IKECOU=0), ces deux variables gardent la même valeur qu’à la fin de la première étape.

• Calcul des termes implicites
On retire à SMBRK et SMBRE la partie en convection diffusion au temps n, qui était stockée dans W7 et W8.
SMBRK=Ω ρ⁽ⁿ⁾k_ts−k⁽ⁿ⁾/Δ t −Ω D(k⁽ⁿ⁾)
SMBRE= Ω ρ⁽ⁿ⁾ε_ts−ε⁽ⁿ⁾/Δ t −Ω D(ε⁽ⁿ⁾)

On calcule les termes implicites, hors convection/diffusion, qui correspondent à la diagonale de la matrice.
TINSTK=Ω ρ⁽ⁿ⁾/Δ t −Ωdiv(ρu)+ΩΓ+ΩMax[−α_k,0]
TINSTE=Ω ρ⁽ⁿ⁾/Δ t −Ωdiv(ρu)+ΩΓ+ΩMax[−α_ε,0]
(Γ n’est pris en compte qu’en injection forcée, c’est-à-dire qu’il est forcément positif et ne risque pas d’affaiblir la diagonale de la matrice).

• Résolution finale
On passe alors deux fois dans codits , pour k et ε, pour résoudre les équations du type :

TINST×(ϕ⁽ⁿ⁺¹⁾−ϕ⁽ⁿ⁾) = D(ϕ⁽ⁿ⁺¹⁾)+SMBR.

• clipping final
On passe enfin dans la routine clipke pour faire un clipping éventuel de k⁽ⁿ⁺¹⁾ et ε⁽ⁿ⁺¹⁾.