Méthodes de relaxation

Dans ce chapitre, nous présentons des méthodes numériques dont le but est d’approcher les solutions du système d’équations diphasiques isotherme compositionnel ½ ∂t(ρck) + ∂x(ρ k LRk L vL + ρ k GRk GvG) = 0, k = 1, K, ∂t(ρv) + ∂x(ρLRLv 2 L + ρGRGv 2 G + p) = 0. (10.1) Nous présenterons deux stratégies de résolution : une méthode avec découplage et une méthode couplée. Nous analyserons ces méthodes du point de vue du respect de la positivité des fractions massiques des constituant ainsi que de leur coût. Nous proposerons deux variantes pour la méthode avec découplage. Nous constaterons que la première (fondée sur une extension naturelle de la méthode de Larrouturou) ne permet pas d’assurer le principe du maximum dans le cadre d’un modèle `a 2 vitesse. Nous constaterons également que la seconde (fondée sur le système de Born-Infeld) ne permet pas de vérifier le principe du maximum pour toutes les lois de glissement, mais seulement pour certaines d’entre elles. Nous considérerons alors un schéma de relaxation fondé sur le modèle couplé et qui a l’avantage de garantir la positivité des fractions massiques quelque soit la loi de glissement. L’idée clé est de construire un système de relaxation sur la base de l’étude du système d’équilibre en coordonnées Lagrangiennes. Nous proposons de relaxer la pression ainsi que les fonctions non linéaires associées au transport de chaque constituant. L’étude de cette méthode de relaxation est facilitée par l’introduction de seulement deux coefficients de relaxation (indépendamment du nombre de constituants). Cette méthode de relaxation permet de diminuer le coˆut d’évaluation du flux numérique. En effet, il n’est désormais plus nécessaire de calculer numériquement les valeurs propres de la matrice Jacobienne, comme c’était le cas dans la méthode VFRoe–Tacite. 

Relaxation avec découplage 

Le système compositionnel étant de grande taille, on cherche `a construire un système de petite taille par découplage. Dans ce but, on approche le système compositionnel par les solutions du système gaz-liquide complété par des équations de transport sur les fractions massiques des constituants. Rappelons en effet que ρck = ρ k LR k L + ρ k GR k G et ρv = ρ k LR k L vL + ρ k GR k GvG. (10.2) L’hypothèse d’équilibre (égalité des potentiels chimiques) permet de calculer ρ k LRk L et ρ k GRk G en fonction de ρck. Nous proposons d’abandonner cette hypothèse et de considérer des lois d’évolution 137 ethodes de relaxation CHAPITRE 10. METHODES ´ DE RELAXATION sur ρ k LRk L et ρ k GRk G. On procède à une “relaxation chimique” en découplant chaque équation de conservation en deux et on considère le système    ∂t ¡ ρ k LRk L ¢ + ∂x ¡ ρ k LRk L vL ¢ = λ Qk, k = 1, K, ∂t ¡ ρ k GRk G ¢ + ∂x ¡ ρ k GRk GvG ¢ = −λ Qk, k = 1, K, ∂t(ρv) + ∂x(ρLRLv 2 L + ρGRGv 2 G + p) = 0, (10.3) avec Qk des termes de relaxation chimiques et λ un paramètre de relaxation. Nous allons approcher les solutions du système à l’équilibre (10.1) par les solutions du système étendu (10.3). Notons u = (ρc1, . . . , ρcK, ρv) T et supposons connue à l’instant t n une solution u(x,t n ). On se propose de réactualiser cette solution en 2 pas de temps dictés par une décomposition en opérateurs. 1er pas On considère la donnée initiale v(x,t n ) = ³© ρ k LRk L ª k=1,K , © ρ k GRk G ª k=1,K , ρv´T construit à partir de u(x,t n ) sous l’hypothèse d’équilibre chimique. On résout le système homogène    ∂t ¡ ρ k LRk L ¢ + ∂x ¡ ρ k LRk L vL ¢ = 0, k = 1, K, ∂t ¡ ρ k GRk G ¢ + ∂x ¡ ρ k GRk GvG ¢ = 0, k = 1, K, ∂t(ρv) + ∂x(ρLRLv 2 L + ρGRGv 2 G + p) = 0. (10.4) 2è pas On résout    ∂t ¡ ρ k LRk L ¢ = λ Qk, k = 1, K, ∂t ¡ ρ k GRk G ¢ = −λ Qk, k = 1, K, ∂t(ρv) = 0, (10.5) dans la limite d’un paramètre de relaxation λ infini avec pour donnée initiale v(x,t n + ∆t) du problème précédent. Ceci revient à construire ρck par ρck = ρ k LRk L + ρ k GRk G. Il est attendu que la procédure d’approximation proposée soit stable dès que la condition de Whitham, imposant que la vitesse du son du système en déséquilibre est plus grande que la vitesse du son du système à l’équilibre, est respectée. On se concentre désormais sur le premier pas, c’est à dire sur la résolution du système (10.4). Nous allons voir pourquoi le premier pas peut être compris comme la résolution de    ∂t(ρLRL) + ∂x(ρLRLvL) = 0, ∂t(ρGRG) + ∂x(ρGRGvG) = 0, ∂t(ρv) + ∂x(ρLRLv 2 L + ρGRGv 2 G + p) = 0, (10.6) complété de la résolution de ½ ∂t(ρXXk) + ∂x(ρXXkv + Xkσ) = 0, k = 1, K − 1, ∂t(ρY Yk) + ∂x(ρY Ykv − Ykσ) = 0, k = 1, K − 1. (10.7) Soulignons que le découplage est intéressant car le système (10.6) est le système gaz-liquide tel qu’il a été étudié dans la première partie. Pour eclairer cette approche, considérons plus précisément le système (10.4). • Si on somme les K premières équations de (10.4) entre elles d’une part, et les K dernières équations de (10.4) entre elles d’autre part, on obtient les équations de conservation de la masse de liquide et de la masse de gaz : ∂t(ρLRL) + ∂x(ρLRLvL) = 0, ∂t(ρGRG) + ∂x(ρGRGvG) = 0, (10.8) par définition de ρLRL et ρGRG. Méthodes de relaxation 138 Michael¨ Baudin CHAPITRE 10. METHODES ´ DE RELAXATION • Par ailleurs, en utilisant les égalités ρ k LRk L = ρLRLXk et ρ k GRk G = ρGRGYk, il est facile de mettre (10.4) sous la forme : ∂t(ρLRLXk) + ∂x(ρLRLXkvL) = 0, k = 1, K, ∂t(ρGRGYk) + ∂x(ρGRGYkvG) = 0, k = 1, K. (10.9) • En reportant les définitions ρX = ρLRL et ρY = ρGRG dans (10.9), on obtient ∂t(ρXXk) + ∂x(ρXXkvL) = 0, k = 1, K, ∂t(ρY Yk) + ∂x(ρY YkvG) = 0, k = 1, K. (10.10) Or les flux peuvent peuvent se réécrire ρXXkvL = ρXXkv + Xkσ, ρY YkvG = ρY Ykv − Ykσ, (10.11) par application des égalités (9.4) et avec la définition σ = ρXY ϕ. Nous obtenons alors ½ ∂t(ρXXk) + ∂x(ρXXkv + Xkσ) = 0, k = 1, K, ∂t(ρY Yk) + ∂x(ρY Ykv − Ykσ) = 0, k = 1, K. (10.12)

 Remarque 10.1.1

En introduisant, dans les équations de conservation (10.9), les équations de conservation de chaque phase (10.8), on obtient des équations de transport : ∂tXk + vL ∂xXk = 0, k = 1, K ∂tYk + vG ∂xYk = 0, k = 1, K, (10.13) qui mettent en évidence 2i invariants de Riemann forts en coordonnées Eulériennes. L’algorithme proposé devra respecter la positivité des fractions massiques sous une condition CFL. La principale difficulté vient du fait que le transport s’effectue aux vitesses vL et vG et non à la vitesse v moyenne. Nous allons présenter maintenant deux méthodes numériques fondées sur le système découplé. Nous verrons qu’une extension de la méthode de Larrouturou ne fonctionne pas. Nous proposerons alors une méthode numérique fondée sur le système de Born-Infeld. Mais il s’avère que cette méthode est limitée par une condition de Whitham qui n’est pas vérifiée par toutes les lois de glissement. 

Une extension de la méthode de Larrouturou

 On s’intéresse ici au système    ∂t(ρLRL) + ∂x(ρLRLvL) = 0, ∂t(ρGRG) + ∂x(ρGRGvG) = 0, ∂t(ρv) + ∂x(ρLRLv 2 L + ρGRGv 2 G + p) = 0, ∂t(ρLRLXk) + ∂x(ρLRLvLXk) = 0, k = 1, K − 1, ∂t(ρGRGYk) + ∂x(ρGRGvGYk) = 0, k = 1, K − 1. (10.14) 10.1.1.1 Construction On procède à une méthode inspirée des travaux de B. Larrouturou [42]. Nous traitons séparément les 3 premières équations des suivantes. Nous supposons qu’`a l’instant n, toutes les quantités conservatives sont connues sur chaque maille : (ρc1, . . . , ρcK, ρv) n i . (10.15) L’algorithme se compose alors des étapes suivantes. Michael¨ Baudin 139 Méthodes de relaxation CHAPITRE 10. METHODES ´ DE RELAXATION 1. Décomposition des variables Nous appelons le modèle thermodynamique de telle sorte que l’on peut calculer les composantes “diphasiques” (ρLRL, ρGRG, ρv) n i (10.16) ainsi que les composantes “compositionnelles” ³ (ρLRLXk)k=1,K , (ρGRGYk)k=1,K´n i , (10.17) pour toutes les mailles i. 2. Equations ´ diphasiques. Nous connaissons un schéma explicite capable de procéder à la mise à jour des trois premières composantes de telle sorte que l’on peut calculer les composantes “diphasiques” (ρLRL, ρGRG, ρv) n+1 i (10.18) ainsi que les flux numériques “diphasiques” H1 i+1/2 = (ρLRLvL) n i+1/2 , H2 i+1/2 = (ρGRGvG) n i+1/2 , (10.19) pour toutes les arêtes i + 1/2. 3. Transport des constituants. On définit les flux numériques H4 i+1/2 = H1 i+1/2 × ( (Xk) n i , si H1 i+1/2 ≥ 0, (Xk) n i+1, si H1 i+1/2 ≤ 0, H5 i+1/2 = H2 i+1/2 × ( (Yk) n i , si H2 i+1/2 ≥ 0, (Yk) n i+1, si H2 i+1/2 ≤ 0. La mise à jour des fractions massiques des constituants dans chaque phase se fait alors grˆace au schéma d’évolution (ρLRLXk) n+1 i = (ρLRLXk) n i − λ ³ H4 i+1/2 − H4 i−1/2 ´ , (10.20) (ρGRGYk) n+1 i = (ρGRGYk) n i − λ ³ H5 i+1/2 − H5 i−1/2 ´ , (10.21) pour les constituant k = 1, K − 1 et pour chaque maille i, avec λ = ∆t n ∆x . 4. Recomposition des variables La réactualisation des fractions massiques des constituants se fait par (ρck) n+1 i = (ρLRLXk) n+1 i + (ρGRGYk) n+1 i , (10.22) pour les constituants k = 1, K −1 et pour toutes les mailles i. Le dernier constituant est évalué par ρ n+1 i = (ρLRL) n+1 i + (ρGRG) n+1 i , (10.23) (ρcK) n+1 i = ρ n+1 i − K X−1 k=1 (ρck) n+1 i . (10.24) Méthodes de relaxation 140 Michael¨ Baudin CHAPITRE 10. METHODES ´ DE RELAXATION Ainsi présenté, le schéma est explicite, c’est à dire précis, vis-`a-vis du transport des constituants. Nous constatons alors un avantage très important de cette technique : le gain en temps calcul. En effet, la taille du système pour lequel il est nécessaire de calculer des flux numériques est réduite : quelque soit le nombre de constituants, nous ne calculons des flux numériques que pour le système gaz-liquide de taille 3. A ce stade, il reste à démontrer que l’étape 4. de mise à jour des fractions massiques des constituants préserve leur positivité. C’est ce point, essentiel, que nous allons développer par la suite. Nous allons constater qu’il n’est pas possible d’assurer cette contrainte, qui est pour nous minimale, dans le cadre d’un modèle à deux vitesses. 

Positivité des fractions massiques des constituants 

Le schéma de B. Larrouturou est construit dans le but de préserver la positivité des fractions massiques dans le cadre d’un système diphasique à une vitesse. Il est important de noter que le système traité par B. Larrouturou est le même que notre système gaz-liquide, si l’on considère le cas particulier du glissement nul, c’est à dire lorsqu’on considère une seule vitesse v = vL = vG. Dans son article [42], il fait la démonstration du fait que les mises à jours successives de la fraction massique de gaz respectent le principe du maximum de telle sorte que, pour toute maille i et pour tout pas de temps n ≥ 0, on a l’égalité min j Y 0 j ≤ Y n i ≤ max j Y 0 j sous réserve que le pas de temps vérifie une condition de type CFL. Si on veut comprendre l’application que l’on en fait ici, il est indispensable d’analyser dans le détail la démonstration qu’il propose. Elle est fondée sur la structure de la mise à jour de la densité et de la fraction massique de gaz. On suppose que le schéma numérique est tel que ρ n+1 i > 0. Il est alors possible de démontrer que la fraction massique à l’itéré n + 1 est une combinaison linéaire positive des fractions massiques à l’itéré n, sous la condition CFL. Donc, si les fractions massiques sont positives à l’itéré n, elles le sont également à l’itéré n + 1. Pour nous, supposons qu’`a l’itéré n, nous avons : ρ n i ≥ 0, Y n i ≥ 0, Xn i ≥ 0, (Yk) n i ≥ 0, (Xk) n i ≥ 0, (ck) n i ≥ 0, (10.25) pour toutes les mailles i et tous les constituants k = 1, K. On suppose que le schéma numérique “diphasique” est tel que ρ n+1 i ≥ 0, Y n+1 i ≥ 0, X n+1 i ≥ 0. (10.26) Pour nous, il y a deux possibilités pour démontrer que c n+1 i = (XXk) n+1 i + (Y Yk) n+1 i ≥ 0. (10.27) Soit on démontre que chaque membre est positif, d’ou` il suit que la somme est également positive, soit on démontre directement cette inégalité. – Empruntons la première voie et tentons de démontrer que (Xk) n+1 i ≥ 0 et (Yk) n+1 i ≥ 0. Par introduction de la définition ρXXk = ρLRLXk dans (10.20), il vient (ρXXk) n+1 i = (ρXXk) n i − λ ³ H4 i+1/2 − H4 i−1/2 ´ , (10.28) Michael¨ Baudin 141 Méthodes de relaxation CHAPITRE 10. METHODES ´ DE RELAXATION et de même pour les fractions massiques de gaz. Supposons maintenant que H1 i±1/2 ≥ 0. Alors, on a (XXk) n+1 i = h (ρX) n i − λ H1 i+1/2 i (Xk) n i + λ H1 i−1/2 (Xk) n i−1 ρ n+1 i (10.29) et de même pour (Y Yk) n+1 i . Nous sommes donc amenés à imposer (ρX) n i − λ H1 i+1/2 ≥ 0 c’est à dire λ H1 i+1/2 ρ n i ≤ Xn i , (10.30) de telle sorte que (XXk) n+1 i ≥ 0. Il faut maintenant détailler le flux numérique sous la forme H1 i+1/2 = (ρXv + σ)i+1/2 ou` σ = ρXY ϕ avec ϕ = vL − vG est le glissement. Nous utilisons le schéma de relaxation défini dans la première partie de cette thèse pour le système gaz-liquide, pour lequel nous avons démontré qu’il préserve la positivité des fractions massiques de liquide et de gaz : Y n+1 i ≥ 0 et X n+1 i ≥ 0. Alors le flux numérique est celui d’un schéma de Godunov H1 i+1/2 = ρ ?X? v ? + Σ ? , (10.31) avec ρ ? ≥ 0 et X? ≥ 0. Supposons dans un premier temps que le glissement est nul. Nous allons vérifier que, dans le cas particulier ϕ = 0 (qui est celui de B. Larrouturou), la condition de positivité (10.30) est λ ρ ?v ? ρ n i ≤ 1, (10.32) avec ρ ? et v ? les états intermédiaires définis par (3.50). Comme le glissement est nul, σ = 0 donc, par définition, Σ ? = 0. Par ailleurs, par hypothèse, H1 i+1/2 ≥ 0 donc v ? ≥ 0. Par conséquent, la fraction massique de liquide X? est “décentrée amont” : X? = Xn i . Ainsi, le flux numérique se simplifie en H1 i+1/2 = ρ ?Xn i v ? . Il est alors évident que si Xn i = 0, alors la condition (10.30) est trivialement vérifiée, puisque le flux H1 i+1/2 est nul. Par opposition, si Xn i 6= 0, alors on peut simplifier (10.30) pour trouver la condition (10.32). Notons que, sur le plan de la stricte logique, il serait désormais nécessaire de compléter cette démonstration dans le cas ou` les flux H1 i±1/2 ne sont pas tous positifs. Mais notre but n’est pas de faire cette démonstration, mais bien de trouver un contre-exemple dans le cas ou` le glissement est non nul : c’est ce que nous allons faire. Supposons maintenant que le glissement est non nul : désormais, Σ ? 6= 0. Nous allons montrer que la condition (10.30) ne peut être vérifiée si Xn i = 0 et Xn i+1 6= 0. En effet, il n’est plus possible de démontrer que le flux numérique est nul. En particulier, si v ? − bτ ? L < 0 et v ? > 0 (avec τ le co-volume défini par τ = 1/ρ), la solution du problème de Riemann est X? = Xn i +Xn i+1 2 + Xn i −Xn i+1 2 b 6= 0. – Empruntons la seconde voie. Supposons maintenant que H 1,2 i±1/2 sont positifs ou nuls. Il vient (ck) n+1 i = h (ρX) n i − λ H1 i+1/2 i (Xk) n i + h (ρY ) n i − λ H2 i+1/2 i (Yk) n i + +λ H1 i−1/2 (Xk) n i−1 + λ H2 i−1/2 (Yk) n i−1 ρ n+1 i , (10.33) Méthodes de relaxation 142 Michael¨ Baudin CHAPITRE 10. METHODES ´ DE RELAXATION = (ρXXk) n i + (ρY Xk) n i − λ h H1 i+1/2 (Xk) n i + H2 i+1/2 (Yk) n i i + +λ h H1 i−1/2 (Xk) n i−1 + H2 i−1/2 (Yk) n i−1 i ρ n+1 i , (10.34) En introduisant la notation (ρck) n i = (ρXXk) n i + (ρY Xk) n i , nous sommes donc amenés à imposer (ρck) n i − λ h H1 i+1/2 (Xk) n i + H2 i+1/2 (Yk) n i i ≥ 0 c’est à dire λ H1 i+1/2 (Xk) n i + H2 i+1/2 (Yk) n i ρ n i ≤ (ck) n i . (10.35) Sous cette condition, on est assuré que (ck) n+1 i ≥ 0. De même que tout à l’heure, si le glissement est nul, la démonstration est possible et si le glissement est non nul, elle ne l’est pas. 

Conclusion 

Nous avons montré qu’une extension naturelle du schéma de B. Larrouturou est possible mais conduit à un schéma qui ne préserve pas la positivité des fractions massiques des constituants. La raison principale est que le système traité par B. Larrouturou concerne un mélange de différents gaz (de telle sorte que ϕ = vL−vG = 0) et non un mélange gaz-liquide (pour lequel ϕ = vL−vG 6= 0). Comme la positivité est pour nous une contrainte centrale, nous abandonnons cette stratégie et nous en proposons une autre qui, comme nous allons le voir, a elle aussi ses limites. 10.1.2 Relaxation découplée par Born-Infeld On considère le système    ∂t(ρLRL) + ∂x(ρLRLvL) = 0, ∂t(ρGRG) + ∂x(ρGRGvG) = 0, ∂t(ρv) + ∂x(ρLRLv 2 L + ρGRGv 2 G + p) = 0, ∂t(ρXXk) + ∂x(ρXXkv + Xkσ) = 0, k = 1, K − 1, ∂t(ρY Yk) + ∂x(ρY Ykv − Ykσ) = 0, k = 1, K − 1. (10.36) Il a déj`a été vu que le système précédent peut se mettre sous la forme    ∂t(ρ) + ∂x(ρv) = 0, ∂t(ρv) + ∂x(ρv2 + P) = 0, ∂t(ρY ) + ∂x(ρY v − σ) = 0, ∂t(ρXXk) + ∂x(ρXXkv + Xkσ) = 0, k = 1, K − 1, ∂t(ρY Yk) + ∂x(ρY Ykv − Ykσ) = 0, k = 1, K − 1, (10.37) avec σ = ρXY ϕ et P = p + ρXY ϕ 2 . En coordonnées Lagrangiennes, définies par le changement de variables, dy = ρdx − ρvdt, le système devient    ∂tτ − ∂yv = 0, ∂tY − ∂yσ = 0, ∂tv + ∂yP = 0, ∂t(XXk) + ∂y(Xk σ) = 0, k = 1, K − 1, ∂t(Y Yk) + ∂y(−Yk σ) = 0, k = 1, K − 1. (10.38) Michael¨ Baudin 143 Méthodes de relaxation CHAPITRE 10. METHODES ´ DE RELAXATION Il est facile de démontrer que : X ∂tXk + σ ∂yXk = 0, k = 1, K, Y ∂tYk − σ ∂yYk = 0, k = 1, K. (10.39) Ces deux dernières équations permettent d’exhiber les invariants forts du précédent système, ainsi que les valeurs propres associées. Nous pouvons faire une analogie entre les systèmes (10.13) et (10.39). Elle montrent que les vitesses de transport des fractions massiques Xk et Yk sont vL et vG en coordonnées Eulériennes et sont σ X et − σ Y en coordonnées Lagrangiennes. On constate que les vitesses caractéristiques − σ Y et σ X peuvent donner lieu à des valeurs propres indéterminées lors de l’évaluation de la solution du problème de Riemann. En effet, si nous relaxons le système de la manière de la première partie de cette thèse, il vient que la solution du problème de Riemann fait apparaˆıtre les états intermédiaires − Σ? Y ? et Σ? X? . Ainsi, si Y ? = 0 ou X? = 0, il devient impossible de définir ces fractions. C’est pourquoi nous proposons d’utiliser, pour relaxer l’équation de conservation de la fraction massique de gaz, d’utiliser la variante de Born-Infeld. Cette proposition se fonde sur l’article en annexe C qui définit le schéma dans le cas d’une équation de conservation scalaire ainsi que dans le cas du système gaz-liquide. Nous introduisons les variables supplémentaires w et z qui représentent − σ Y et σ X , c’est à dire les vitesses des phases liquide et gaz en coordonnées Lagrangiennes. Tous les champs de ce système sont linéairement dégénérés. On peut envisager d’appliquer un solveur de Godunov. 

Conclusion

 Il a été démontré (voir l’annexe C) que cette stratégie nous contraint à vérifier que la loi de glissement vérifie une condition de Whitham. Or, nous cherchons une méthode numérique qui permette d’assurer la robustesse du schéma quelque soit la loi de glissement. C’est pourquoi cette stratégie a été abandonnée.