Brevet FR2756045A1 — publié le 19980522

La présente invention concerne une méthode pour former un modèle de simulation d'écoulements diphasiques dans des conduites et notamment des écoulements d'hydrocarbures dans des conduites pétrolières. Par le brevet US A-5 550 761, on connaît une méthode pour modéliser des écoulements de mélanges diphasiques dans des conduites où l'on utilise un système d'équations aux dérivées partielles de modélisation comportant des équations de conservation de la masse pour chaque phase liquide et gaz du mélange ainsi que de la quantité de mouvement totale du mélange, qui s'écrivent: #t[#GRG]+#x[#GRGVG]=0 #/#t[#LRL]+#/#x[#LRLVL]=0 (1) #/#t[#GRGVG+#LRLVL]+#/#x[#GRGVG2+#LRLVL2+P]=Tw-#gsin# où RG (resp. RL) est la fraction volumique de gaz (resp. liquide), VG (resp. VL) est la vitesse de la phase gaz (resp. liquide), P est la pression. p = #GRG + PLRL est la densité moyenne du mélange, T est une donnée fonction des inconnues qui représente le terme de frottement pariétal, g est l'accélération de la pesanteur et 0 est l'angle que fait la conduite avec l'horizontale. On y ajoute des équations de fermeture de trois sortes 1) une loi de fermeture "conservative" qui s'exprime sous la forme: RL + RG = 1 2) une loi de fermeture hydrodynamique qui s'écrit formellement: #(VM,XG,#(P,T),dV,x)=0 où: - VM est la vitesse moyenne du mélange diphasique, #GRG - XG est la fraction massique de gaz (XG = ), #GRG+#LRL - r(P, T) représente l'ensemble des propriétés physiques des fluides, - T est la température, - dV=VG-VL est la vitesse de glissement entre les deux phases. 3) des lois de fermeture thermodynamiques Le système d'équations aux dérivées partielles (1) s'écrit sous forme vectorielle conservative : #/#t(W)+#/#x(F(W,x)=Q(W,x) où F et Q s'expriment en fonction de W par l'intermédiaire des lois hydrodynamiques et thermodynamiques mais nous ne disposons pas d'expression analytique de F(W,x) et de Q(W,x). La matrice jacobienne du système en W est notée : A(W,x)=#F/@@(W,x) aw Les valeurs propres #1, #2 et X3 de A(W,x) sont réelles et vérifient A1 (W) O. En outre, #1 et #3 sont d'au moins un ordre de grandeur supérieur à X2. La modélisation des écoulements diphasiques demande la résolution d'un système d'équations aux dérivées partielles hyperbolique non linéaire relativement complexe. Il est connu pour résoudre ce type de système hyperbolique non linéaire d'utiliser un schéma soit explicite, soit implicite. L'inconvénient d'un schéma implicite est son manque de précision. Un schéma explicite conduit à des résultats beaucoup plus précis mais en contrepartie, il requiert une limitation sur la grandeur du pas de temps qui est difficile à appliquer dans le cadre d'une application à des écoulements diphasiques dans des conduites en raison de contraintes particulières. Comme indiqué ci-dessus, le système d'équations possède la caractéristique d'avoir des valeurs propres dont les ordres de grandeur sont très différents. Les valeurs propres X1 et B3 sont liées aux vitesses des ondes de pression dans la conduite qui sont de l'ordre de quelques dizaines à quelques centaines de m/s, beaucoup plus grandes donc que celle des fluides qui est de l'ordre de quelques m/s. L'autre valeur propre B2 est liée à la vitesse de propagation de la fraction gazeuse dans la conduite. Avec un schéma explicite, le pas de temps est réglé par la plus grande valeur propre et il est donc relativement petit en l'occurrence dans l'application traitée. Cette contrainte n'est pas compatible avec une résolution explicite si l'on veut que la valeur du pas de temps ne soit pas trop petite pour rester dans des durées de calcul raisonnables. Une solution pour éviter la difficulté peut consister, comme il est connu, à adopter un schéma de résolution implicite qui permet de choisir un pas de temps plus grand. La contrepartie de cet avantage est que le schéma implicite introduit une diffusion numérique qui rend les résultats obtenus moins précis, notamment pour les ondes lentes. La méthode de modélisation d'écoulements de mélanges diphasiques dans des conduites selon l'invention permet d'éviter les difficultés mentionnées ci-dessus. Elle comporte la résolution d'équations de conservation de masse pour chacune des phases du mélange ainsi qu'une équation de conservation de quantité de mouvement du mélange, en respectant des relations de fermeture. Elle est caractérisée en ce qu'elle comporte l'utilisation d'un schéma de résolution explicite sur les ondes lentes correspondant à la propagation de la fraction gazeuse dans la conduite et d'un schéma de résolution implicite sur les ondes rapides correspondant à la propagation des ondes de pression, de façon à permettre un contrôle du pas de temps par les ondes lentes tout en gardant la précision sur ces mêmes ondes inhérente au schéma explicite. D'autres caractéristiques et avantages de la méthode selon l'invention, apparaîtront à la lecture de la description ci-après de modes de réalisation décrits à titre d'exemples non limitatifs, en se référant aux dessins annexés où - la Fig. l montre des courbes de variation en fonction du temps t du débit massique de liquide, en différents points d'abcisse x d'une conduite, quand on impose à l'entrée (amont) une certaine fluctuation du débit massique de gaz, en adoptant pour la modélisation un schéma implicite; - la Fig.2 montre dans des conditions identiques, les courbes obtenues en adoptant un schéma de calcul explicite; et - la Fig.3 montre toujours dans les mêmes conditions, les courbes que l'on obtient en adoptant le schéma mixte de la méthode selon l'invention. DISCRETISATION Pour la modélisation de l'écoulement, on discrétise la conduite en un ensemble de I mailles. La maille i est notée Mi et ses extrémités sont les points xi-1.2 et xi+1/2. On pose #xi = xi+1/2 - xi-1/2, xi = y2(xi+1/2 - xi-1/2). Les inconnues discrètes Wi, i=1,...,I, sont localisées aux centres des mailles. Un schéma de résolution par volumes finis est obtenu en intégrant le système (2) sur chaque maille Mi : d/dt(Wi)#xi + Hi+1/2 - Hi-1/2 = Qi#xi où Hi+1/2 est le flux numérique en xi+1/2 à l'instant t, approximation de F en xi+1/2 : Hi+1/2=F(Wi-p+1,..., Wi+p-1, Wi+p). On utilise, pour la modélisation, le schéma suivant: Hi+1/2 = 1/2(F(Wi) + F(Wi+1)) + KDi+1/2(Wi-Wi+1) (3) Dans cette relation, on considère le flux comme la somme d'un terme centré (F(Wi)+ F(Wi+1))/2 et d'un terme KDi+1/2(Wi - Wi+1) où Di+112 est une matrice de diffusion, construite à partir des matrices jacobiennes Aj=A(Wi) et Ai+l=A(Wi+l) et K est un paramètre réel positif. Après cette discrétisation en espace, le système d'équations s'écrit formellement: d dt(Wi)=Li(Wj,J=1..I) (4) Il reste à effectuer la discrétisation en espace. Etape explicite Le schéma d'Euler explicite que l'on peut utiliser pour réaliser cette discrétisation temporelle, et qui s'écrit Win+L - wn = Li(Wj,j = dt requiert une limitation sur le pas de temps de type CFL contrôlée par la plus grande des valeurs propres. Pour le schéma de résolution de la relation (3), cette limitation s 'écrit t CFL Min( Min(Axi (5) yKMax 0 où CFL représente un nombre positif inférieur ou égal à 1. Etape implicite La condition imposée par la relation (5) est trop restrictive en pratique et conduit à des temps de calcul importants, car le pas de temps y est réglé par la plus grande valeur propre k3, elle-même liée à la vitesse des ondes de pression dans la conduite qui est de l'ordre de quelques dizaines à quelques centaines de m/s alors que la vitesse des fluides dans la conduite est de l'ordre de quelques mètres par seconde. Cette contrainte n'est pas compatible avec une résolution explicite si l'on veut que la valeur du pas de temps ne soit pas trop petite pour rester dans des durées de calcul raisonnables. Pour éviter d'avoir à impliciter le schéma numérique qui certes autorise un pas de temps plus grand mais produit en contrepartie des résultats bien moins précis, la méthode selon l'invention propose donc d'adopter une version explicite sur les ondes lentes associées à la valeur propre #2 correspondant à la propagation de la fraction gazeuse dans la conduite et implicite sur les ondes rapides de pression associées à #1 et B3. Schéma d'ordre 1 en temps, En adoptant un schéma d'Euler implicite, la relation (4) s'écrit: Win+1 - Win = L:(W:n+1, i=-1..I) At Cette discrétisation n'est pas soumise à une condition de type CFL mais elle demande à chaque pas de temps la résolution d'un système non linéaire entre les inconnues(Wjn+l ,j = 1..I). Comme cette résolution est en général coûteuse, on utilise des schémas linéairement implicites consistant à approcher L; (Wjn+ ,j = 1.. I) par une combinaison linéaire de (Wjn+1,j=1..I). Le système d'équations (2) discrétisé à l'ordre 1 en espace s'écrit: Li(Wjn+1, j=1..I) = -1/#x:(Hi+1/2(Win+1, Wi+1n+1) - Hi-1/2(Wi-In+1, Win+1))+Qin+1 H(Win+1, Wi+1n+1) est linéarisé en faisant un développement de Taylor à l'ordre 1 par rapport à W1n et Ce schéma est linéairement implicite pour toutes les ondes. Pour garder une bonne précision sur les ondes lentes, on n'implicite dans H(Win+1, Wi+1n+1) que les quantités relatives aux ondes rapides. On ne garde alors que les termes liées aux valeurs propres #1 et X3 dans les dérivées partielles de H. Le schéma ainsi obtenu est soumis à une condition CFL beaucoup moins restrictive que celle de la relation (5) car seule l'onde lente intervient: (dx,,dx,,l . At =I, ,I | KMax(2i) IJ Le pas de temps est donc beaucoup plus grand qu'avec le schéma explicite et permet ainsi d'obtenir une meilleure précision sur les fronts se propageant avec la vitesse X2. Ordre 2 en temps La généralisation du schéma à l'ordre 2 en temps se fait de manière connue par un schéma de Runge-Kutta (schéma de Henn) que l'on trouve décrit notamment par Le Veque R.J. 1992. Numerical methods for conservaflon laws. Lectures in Mathematics. ETH Zurich. Birkhauser. Résultats de simulation La méthode selon l'invention a été testée en réalisant des simulations d'écoulements dans des conduites. Les résultats présentés ci-après à titre d'exemple illustratif, sont relatifs à l'effet produit sur le débit de gaz à la sortie d'une conduite horizontale quand on impose les conditions aux limites suivantes: Limite amont: QL(t)=20kg/s pour t 2 Os QG(t)=0,2 + 0,02t pour 0 QG(t)=0,4 pour t > 10s QL et QL représentant respectivement le débit massique de liquide et de gaz. Limite aval P(t) = 1 MPa pour t 2 Os Le pas d'espace utilisé est constant égal à 200m. La quantité qui intéresse les pétroliers est le débit liquide en sortie de conduite. Les courbes représentent l'évolution du débit liquide au cours du temps en différents points de la conduite et en particulier à la limite aval. On compare les résultats obtenus pour le schéma explicite, le schéma implicite et le schéma explicite sur les ondes de taux de vide (propagation de la fraction gazeuse) et implicite sur les ondes de pression. Les temps calculs sont les suivants Schéma explicite : 8704s Schéma implicite : 205s Schéma explicite/implicite : 2 1 8s. Comme on le voit, le schéma explicite/implicite requiert un temps calcul qui est du même ordre de grandeur que le temps calcul du schéma implicite mais plus de 40 fois inférieur au temps calcul du schéma explicite. On observe aussi sur les courbes des Fig. 1 à 3 que la précision obtenue est équivalente que l'on utilise le schéma explicite ou le schéma mixte explicite/implicite de la méthode selon l'invention, et qu'elle est bien supérieure à la précision que l'on peut obtenir avec un schéma implicite. REVENDICATIONS 1) Méthode pour constituer un modèle de simulation d'écoulements de mélanges diphasiques dans une conduite comportant la résolution d'équations de conservation de masse pour chacune des phases du mélange ainsi qu'une équation de conservation de quantité de mouvement du mélange, en respectant des relations de fermeture, caractérisée en ce qu'elle comporte l'utilisation d'un schéma de résolution explicite sur les ondes lentes correspondant à la propagation de la fraction gazeuse dans la conduite, et d'un schéma de résolution implicite sur les ondes rapides correspondant à la propagation des ondes de pression, de façon à permettre un contrôle du pas de temps par les ondes lentes tout en gardant la précision sur ces mêmes ondes inhérente au schéma explicite.