IBP2015_12 UM MÈTODO DOS VOLUMES FINITOS CENTRADO NA CÉLULA PARA SIMULAÇÃO DE ESCOAMENTO BIFÁSICO EM RESERVATÓRIOS DE PETRÓLEO HETEROGÊNEOS E ANISOTRÓPICOS Contreras, F.R.L.1, Cavalcante T. M.2, Carvalho, D. K. E.3, Lyra, P. R. M. 4
Copyright 2012, Instituto Brasileiro de Petróleo, Gás e Biocombustíveis - IBP Este Trabalho Técnico foi preparado para apresentação na Rio Oil & Gas Expo and Conference 2012, realizado no período de 17 a 20 de setembro de 2012, no Rio de Janeiro. Este Trabalho Técnico foi selecionado para apresentação pelo Comitê Técnico do evento, seguindo as informações contidas no trabalho completo submetido pelo(s) autor(es). Os organizadores não irão traduzir ou corrigir os textos recebidos. O material conforme, apresentado, não necessariamente reflete as opiniões do Instituto Brasileiro de Petróleo, Gás e Biocombustíveis, Sócios e Representantes. É de conhecimento e aprovação do(s) autor(es) que este Trabalho Técnico seja publicado nos Anais da Rio Oil & Gas Expo and Conference 2012.
Resumo No presente artigo, apresentamos uma metodologia numérica para discretização de equações que modelam escoamentos bifásicos (óleo-água), em reservatórios de petróleo, heterogêneos e anisotrópicos, usando um Cell-Centered Finite Volume Method (CCFVM). Os CCFVM têm, como principais vantagens, a sua fácil implementação e a relação de equivalência que existe entre célula e o volume de controle. As equações governantes são modeladas com o uso da metodologia IMPES (Implicit Pressure Explicit Saturation), que é um tipo de método segregado, fácil de implementar, em que o fenômeno é descrito a partir de uma equação elíptica de pressão na forma implícita, e de uma equação hiperbólica de saturação na forma explicita. Para discretizar a equação hiperbólica, foi usado o First Order Upwind Method e, para a equação de pressão, implementamos um método de CCFVM proposto por Zhiming Gao e Jiming Wu (ZHIMING, et al., 2010). Essa formulação, usada na equação de pressão, é um tipo de MPFA (Multipoint Flux Aproximation) que se encaixa num grupo de métodos conhecidos como Diamond Methods, cujo stencil para a construção da expressão do fluxo numa face consiste nos pontos de colocação das células que compartilham aquela face e nos nós que fazem parte dela. Nessa formulação, os nós que fazem parte de certa face são interpolados usando-se as soluções nos pontos de colocação das células que concorrem nesses nós, reduzindo-a assim a uma formulação completamente cellpr essão tem as seguintes características: taxas de convergência de centered . O método usado para discretizar a equação de pressão segunda ordem tanto para as soluções e de primeira pr imeira ordem para o fluxo, pode ser usado sobre malhas poligonais em e m geral, essencialmente o método tem uma boa acurácia e é computacionalemente eficiente e de baixo custo. Para validar a formulação, resolvemos alguns problemas de referência. Os resultados obtidos são promissores, indicando o potencial da formulação adotada.
Abstract In this paper, we present a numerical methodology for discretization of equations that model the two-phase fluid flow (oil-water) in oil reservoirs (considering anisotropic and heterogeneous domains), using a Cell-Centered Finite Volume Method (CCFVM). The CCFVM have, as its main advantages, an easy implementation and the equivalence relation that exists between the cell and the control volume. The governing equations are modeled using the IMPES methodology (Implicit Pressure and Explicit Saturation), which is a type of segregated segregated method that is easy to implement, in which the phenomenon is described by from an elliptical pressure equation, discretized implicitly, and a hyperbolic saturation, discretized explicitly. To discretize the hyperbolic equation the First Order Upwind Method was used. For the pressure equation was implemented a CCFVM method proposed by Zhiming Gao and Jiming Wu (ZHIMING, et al., 2010). The formulation, used in pressure equation, is a type of MPFA method (Multipoint Flux approximation) that fits into a group of methods known such as Diamond Methods, Methods, which stencil for the construction of the expression of the flow on one interface are the colocation points of the cells that share the interface and the nodes that are part of it, but, in this case, the nodes that are part of the interface are interpolated using the solutions at the colocation points of the cells that share each one of these nodes, reducing the scheme to completely cell-centered formulation. The method used to discretize the pressure equation has the following characteristics: it has has second order rates of convergence convergence for the solutions of the equation and for the velocities on the faces, can be used on general polygonal meshes, has essentially a good accuracy and low computational cost. To validate the formulation, some model problems are solved. The results are
______________________________ ¹ Mestrando, Dept. Engenharia Engenharia Mecânica Mecânica ² IC, Dept. Engenharia Mecânica ³ Doutor, Dept. Engenharia Mecânica 4 PhD, Dept. Engenharia Mecânica
Rio Oil & Gas Expo and Conference Conference 2012
promising, indicating the potential of the formulation described.
1. Introdução A simulação de escoamento de fluidos em reservatórios de petróleo tornou-se, nos últimos anos uma atividade de pesquisa importante, em virtude dos resultados promissores alcançados e do alto grau de versatilidade para lidar com configurações complexas. A simulação numérica de reservatórios consiste basicamente no uso de ferramentas computacionais para solucionar as equações diferenciais que governam os fenômenos de transporte de fluidos em meio porosos altamente heterogêneos e anisotrópicos. A metodologia que empregamos para discretizar as equações é um método de CCFV (Cell-Centered Finite Volume) que é localmente e globalmente conservativo, por ser uma formulação de volumes finitos. O método CCFV é bastante popular na área de engenharia e física computacional devido à facilidade de implementação e à robustez no tratamento de malhas não estruturadas e de geometrias complexas. As equações diferenciais que governam o escoamento bifásico (óleo-água) num reservatório de petróleo são equações elípticas, para a pressão, e equações hiperbólicas, para a saturação. A equação elíptica de pressão caracteriza-se pelo tensor de permeabilidade, que geralmente é heterogêneo e anisotrópico, devido às descontinuidades de camadas e falhas geológicas existentes nesse meio. A equação elíptica é discretizada através de um método de CCFV proposto por Zhiming Gao e Jiming Wu (ZHIMING, et al., 2010). São características de grande destaque desse método o fato dele se encaixar no grupo de métodos MPFA conhecido como Diamond Methods e o fato de haver o desenvolvimento de expressões explícitas para a interpolação das soluções nodais, por isso, ao longo desse artigo, nos referiremos a ele como MPFA-DEW (Multipoint MPFA-DEW (Multipoint Flux Approximation - Diamand with Explicit Weigths). A interpolação nodal é feita para que a formulação se torne completamente centrada na célula ( cell-centered ) a partir da combinação linear das soluções nos pontos de colocação das células que concorrem no nó em questão. Os coeficientes (pesos explícitos) dessa combinação linear são derivados de forma a tornar o método robusto e capaz de reproduzir uma solução linear por partes p artes com exatidão (linearity-preserving). Outros autores já propuseram diversas maneiras de calcular os pesos de interpolação, como (WU, et al., 2003), (HUANG, et al., 1998), (GX, et al., 2007), mas a maioria delas somente atende às expectativas (em termos de acurácia e convergência) quando o meio é isotrópico, além disso, dependem das descontinuidades do domínio e do tipo de malha, ao contrário dos pesos de interpolação derivados no presente artigo. Recentemente, alguns métodos CCFV linearity-preserving foram propostos em vários trabalhos (WU, 2005), (JIMING, et al., 2009), (ZHIMING, et al., 2010) e (WU, et al., 2011). A equação hiperbólica de saturação é caracterizada pelas curvas características em 2-D. Para discretizar a equação de saturação, o FOUM (First Order Upwind Method) foi usado, que é um método altamente versátil no tratamento das descontinuidades e que produz soluções não oscilatórias perto de regiões de choque (HIRSCH, 1994), porém com baixa taxa de convergência. A discretização temporal é derivada pelo EEM (Explicit Euler Method). As equações que governam o escoamento bifásico de óleo-água num meio poroso heterogêneo e anisotrópico serão resolvidos usando a metodologia IMPES ( Implicit Pressure and Exlicit Saturation), esta técnica resolve as equações de transporte de fluido de maneira segregada e sequencial; onde o campo de pressões é resolvido de maneira implícita a partir de uma saturação saturação inicial, daí obtemos as vazões vazões nas interfaces; então o campo de saturação é calculado explicitamente. Este procedimento é repetido até que se alcance um tempo ou PVI (Pore Volume Injected) final. Os exemplos foram adotados de literaturas ligadas à indústria de petróleo e mostram a boa acurácia e a eficiência do método MPFA-DEW. MPFA-DEW.
2. Formulação Matemática Nesta seção descrevemos as equações que governam o escoamento bifásico (água - óleo) num reservatório de petróleo, considerando que os fluidos são imiscíveis e incompressíveis e o meio poroso é um corpo rígido (porosidade constante). As equações governantes são obtidas a partir da Lei de Darcy e da Lei de Conservação da Massa. A velocidade de Darcy é obtida desprezando os efeitos de gravidade e capilaridade para cada fase i: vi i K P, i ki i
(1)
onde o tensor de permeabilidade é denotado por K . i e k i são a mobilidade e a permeabilidade relativa da fase i (=o,w), respectivamente. A equação da Lei de Conservação de Massa para um escoamento bifásico é dada por : i S i . vi qi t
(2)
2
Rio Oil & Gas Expo and Conference Conference 2012
onde , , i são a porosidade e a densidade, que são propriedades intrínsecas da rocha e do fluido, respectivamente; S i é a saturação da fase i e, por último, temos o termo de fonte específico, em massa, denotada por qi . Utilizando as equações anteriores, após uma série de manipulações, podemos chegar a:
v Q
onde
S w t
F Sw Q
F ( Sw ) f w (Sw )v ,
T o w f w w T
2.1.
(3)
representa o fluxo de escoamento bifásico
é a mobilidade total. A vazão volumétrica é dada por
(4) v vo vw K p T Q qo o qw w ,
é a velocidade total e a função fluxo fracional
vai depender da saturação da fase água.
Condições Iniciais e de Contorno
Seja o domínio computacional dividido em células finitas. As condições de contorno, para a equação de pressão, podem ser escritas da seguinte forma: p ( x, y ) g D , ( x, y) D T K p g N , ( x, y ) N
(5)
As condições iniciais e de contorno, para a equação de saturação, podem ser escritas da seguinte forma: Sw S w , sobre I S w ( x, 0) S w0 , sobre
(6)
onde I é o conjunto de elementos que correspondem ao poço injetor.
3. Formulação Numérica 3.1. Equação de Pressão Implícita Antes de começar a formulação discreta da pressão implícita, introduzimos as mobilidades na face de cada volume de controle como a média aritmética das mobilidades dos nós da mesma face. O cálculo das mobilidades dos nós é feito da seguinte forma: N ( I )
I
T (k )k
k 1
N
(7)
k
k 1
onde , I é a mobilidade no nó I, N(I) é número de elementos na vizinhança do nó I , k é o volume do k -ésimo -ésimo elemento e T ( k ) é a mobilidade total do k -ésimo -ésimo elemento, as mesmas especificações para a equação (8). N ( J )
J
T (k )k
k 1 N ( J )
(8)
k
k 1
onde, J é a mobilidade no nó J, N(J) é o número de elementos na vizinhança do nó J Finalmente a mobilidade na interface é dada por:
3
Rio Oil & Gas Expo and Conference Conference 2012
IJ
I J 2
(9)
Por tanto a nova permeabilidade é expressa como: K IJ K
Integrando a equação (3) na célula k e utilizando o teorema de divergência, temos:
v d
v Q
k
k
Q d
k
vn
ek e
k
k ,e
ds Q k
(10)
onde k representa o volume do k -ésimo -ésimo elemento em questão, nk ,e é o vetor normal a face e. O objetivo nesta seção é expressar a integral da equação (10) de forma discretizada, para cada face da malha computacional. Antes de prosseguirmos enunciamos um Lema fundamental proposto por (ZHIMING, et al., 2010). LEMA. Seja p uma função de pressão definida sobre o triangulo ABC (com vértices A, B e C ) ordenados em sentido anti-horário. Então: p
p A pB AB
AB
AB AB
2
(11)
p B pA cot BCA pC pA cot ABC ABC
Pode-se provar este lema assumindo que p é uma função é da forma for ma p ( x, y ) a x b y c no triângulo ABC , a prova deste Lema pode ser encontrada em (ZHIMING, et al., 2010). Considerando um modelo bidimensional introduzimos as notações, de acordo com a Figura 1. Dada uma face anti -horário, temos: IJ , de células cujos nós são ordenados em sentido anti-horário,
Figura 1. Representação de dois volumes de controle e a célula diamante. 1) L , e R são as células esquerda e direita, respectivamente. 2) As faces IJ e JI pertencem às células L e R , respectivamente. 3) Os vértices da face serão denotados por I e J . ˆ
ˆ
ˆ
ˆ
4) N IJ IJ e N JI JI são os vetores normais às faces IJ e JI , respectivamente, e é um operador que gira o vetor em 90º no sentido horário. 5) As alturas dos baricentros em relação à face são denotadas como h IJ , L e h JI , R para as células ˆ
respectivamente. 6) K ( Lˆ ) K Lˆ é a permeabilidade da célula
ˆ
L ˆ
e
R ˆ
,
L. ˆ
7) p( L) p L e p ( R ) pR são as pressões nos baricentros das células. ˆ
ˆ
ˆ
8)
p( I ) p I e p( J ) pJ
ˆ
, são as pressões nos vértices da face. 4
Rio Oil & Gas Expo and Conference Conference 2012
Com estes notações o método pode ser descrito, conforme o que se segue, segundo o trabalho de (ZHIMING, et al., 2010).
3.2. Desenvolvimento das Expressões Discretas para as Vazões na Interface Em princípio, derivamos o fluxo na face ( n)
L
q IJ KIJ ˆ
onde
( p
ˆ
ˆ
ˆ
IJL : ˆ
(t ) lIJ KIJ ( p J pI )
pl ) cot IJl ( pJ - pl ) cot
I
q I LJ K Lp N IJ
IJ , usando o lema anterior no triângulo
ˆ
ˆ
(12)
é o fluxo na face IJ .
ˆ
( t ) ( n) As constantes K IJ , K IJ podem ser obtidos através da seguinte combinação linear:
(t ) ( n) K LT ( N IJ ) K IJ IJ KIJ ( NIJ )
(13)
ˆ
Multiplicando a equação (13) pela direita por
IJ
obtemos:
T (t ) ( N IJ ) K L ( IJ ) K IJ 2 IJ
(14)
ˆ
E multiplicando a equação (13) pela direita por IJ obtemos: T ( n) ( N IJ ) K L ( N IJ ) K IJ 2 IJ
(15)
ˆ
Logo, aplicando a integral de linha sobre a face IJ na equação (9) e substituindo os valores das cotangentes, que no fundo estão relacionadas com as alturas e ao comprimento da face IJ , obtemos: L ˆ
h IJ
( n) K IJ IJ
L q IJ dN ˆ
(t ) K IJ L ( JL JI ) ( IL IJ ) pJ pI (n) hIJ ( pI p L ) L ( pJ - pL ) L IJ K IJ JI JI h h IJ IJ
1
ˆ
ˆ
ˆ
ˆ
ˆ
De maneira análoga, calculamos a integral de linha na face JI do triângulo R ˆ
h JI
(n)
K JI
R
q JI dN ˆ
JI
(16)
ˆ
ˆ
JIl
ˆ
:
(t ) K ( JR JI ) ( IR IJ ) ( p I p R ) pJ pI (JI n) h RJI ( pJ - pR ) R R JI K JI IJ hIJ JI hJI
1
ˆ
(17)
ˆ
ˆ
ˆ
ˆ
ˆ
ˆ
Pela continuidade de fluxo temos:
q
ˆ L IJ
(18)
ˆ
R dN qJI dN
IJ
JI
L R Substituindo a integral q JI dN por q IJ dN em equação (17), se tem: ˆ
ˆ
JI
R h JI ˆ
(n) K JI JI
L dN q JI ˆ
IJ
( t ) K JI ( JR JI ) ( IR IJ ) ( p I p R ) pJ pI ( n) ( pJ - pR ) R R JI K JI h h I J J I JI JI
1
ˆ
ˆ
ˆ
ˆ
ˆ
ˆ
R ˆ
hJI
(19)
5
Rio Oil & Gas Expo and Conference Conference 2012
Finalmente, tomando a média ponderada das equações (16) e (19), temos a vazão na interface IJ : F IJ
q LIJ dN K (IJ1) p R pL D (IJ1) pJ pI
(20)
ˆ
ˆ
ˆ
IJ
onde K IJ(1)
(n) (n ) K IJ K JI (n ) R (n ) L K IJ hJI K JI hIJ ˆ
IJ e
ˆ
DIJ(1)
LR IJ ˆ
ˆ
2
IJ
(t ) 1 K IJ
IJ
h L ( n ) IJ K IJ ˆ
(t ) K JI
hR ( n ) JI K JI ˆ
.
(1)
Note que a parte tangencial de D IJ desaparece quando o tensor de permeabilidade é isotrópico, então a expressão (20) fica idêntica àquela representada em (JIMING, et al., 2009). Lembremos que todas as aproximações anteriores são obtidas de forma que o método seja Linearity-Preserving Linearity-Preserving. O tratamento das vazões numa face que esteja sobre o contorno de Dirichlet é dada, a partir da equação (16): ( n)
F IJ
K IJ L h IJ ˆ
JL IJg IJ ˆ
D
( I ) LI IJg D (J ) IJ
2
ˆ
(21)
(t )
p L g D (J ) g D (J ) K IJ ˆ
onde os nós ( g D ( J ) pJ , g D ( I ) pI ) tem pressão conhecida e não é necessário interpolá-los. Quando as faces coincidem com uma fronteira de Neumann, a vazão é expressa da seguinte forma: L q IJ ds g N ds
(22)
ˆ
IJ
IJ
L
onde K L p nIJ g N . Geralmente, no presente trabalho, estes tipos de fluxos serão nulos ( g N qIJ 0) . ˆ
ˆ
3.3. Tratamento dos Vértices Na equação (20), usamos as pressões nos vértices que compõem a face em que calculamos o fluxo, mas para tonar a formulação completamente centrada na célula, devemos escrever essas pressões como combinações lineares das pressões nos pontos de colocação das células na vizinhança dos nós, para um nó I : p I w1 p1ˆ ... wk 1 pk 1 w ˆ p ˆ k
k
N ( I )
w p i 1
iˆ
iˆ
(23)
Existem duas formas de obtenção desses pesos de interpolação. Inicialmente, definimos duas constantes geométricas, baseadas na Figura 2, considerando k como o ponto médio da face Ik : ˆ k ,1
k ,2 ˆ
Ik ˆ k h Ik
I k 1 h Ik ,k 1
(24)
(25)
ˆ
6
Rio Oil & Gas Expo and Conference Conference 2012
(a)
(b)
Figura 2. (a) Elementos na vizinhança do vértice I , (b) representação dos angulos e alturas na célula k , adpatado do artigo (ZHIMING, et al., 2010) . Calculo dos pesos explícitos de tipo 1(MPFA-DEW1)
Os pesos explícitos de tipo 1 podem ser definidos da seguinte forma: wk ˆ
k ˆ
1, 2, 2, ..., N ( I )
, k
N ( I )
(26)
ˆ k
k 1
onde: ˆ K (ˆn ) ˆ
( kˆ ) k ,1 k ,1
k
2
K (ˆn ) ˆ (k 1) K ˆ(n ) cot( ˆ ˆ k ,2 k ,2
i 1
k ,i
k ,i
k, i
)
K ˆ(t ) K ˆ(t ) k ,1
k ,2
(27)
e: K k( n-1),2 cot k -1,2 K ( n,1) cot ,1 K k ( t-)1,2 - K ( t,1)
k k ( k ) (n) ( n) K k -1,2 cot k -1,2 K k ,1 cot k ,1 ˆ
k (t ) (t ) K k -1,2 - K k ,1
(28)
e:
K ( t ),i
T
ˆ
ˆ
ˆ
ˆ
ˆ
( n)
K ˆ ,i
i K 2 i
N i
N ˆ i ˆ
ˆ
ˆ
(29)
ˆ
T
K ˆ N ˆ i ˆ 2 ˆ ˆ i
, com i
1, 2
(30)
Calculo dos pesos explícitos de tipo 2 (MPFA-DEW2)
Os pesos explícitos de tipo 2 podem ser definidos da seguinte forma:
7
Rio Oil & Gas Expo and Conference Conference 2012
w k ˆ
k
(31)
( n) k ,2 k ,2
(32)
, k 1, 2, ..., N ( I )
ˆ
N ( I )
k ˆ
k 1
onde:
k
K
( n) (k ) ( k k ,1 k ,1 ˆ
ˆ
1) K
ˆ
ˆ
ˆ
e: K
( k )
(n) k 1
( n) K k 1, 2
K
(t ) k
ˆ
cot
(n)
K k 1 cot k 1,1 K K
k 1, 2
( n) k ,1 ˆ
cot k ,1 ˆ
( n)
c ot
k ( t ) K k 1, 2 ˆ
k ,2
(33)
ˆ
K
(t ) k ,1
ˆ
e:
( n) K
T
T
N 1
N 1
(t ) K
K N 1 2 1
(34)
K 1
1
(35)
2
3.4. Equação de Saturação Explicita Integrando a equação (4) sobre o um volume de controle k e usando o teorema de Gauss, podemos escrever:
S w( k )
k
t
d k
( f w v ) nk d k
k
Q d w
k
(36)
k
Logo, usando o Explicit Euler Method na discretização temporal da equação (36) temos: t n 1 n S w ( k ) S w( k ) Q w f w F k k , k k
(37)
onde Fk , qk d e qk v nk , representam a vazão e fluxo na face , respectivamente. ˆ
ˆ
As saturações no instante n e n+1 são denotadas por S wn ( k1)
e S wn ( k ) , respectivamente. O cálculo do passo de tempo
t
é feito respeitando o critério de CFL (Courant-Friedrichs-Lewy Condition). Aproximamos os termos ( f w F k , ), da somatória da equação (37), pelo Firts Order Upwind Method da seguinte forma: f w F k ,
f ˆ Fk , , w,L f F , w,Rˆ k ,
0 F k ,
(38)
0 F k ,
onde f w,L e f w, R representam o fluxo fracional do elemento a esquerda e a direita, respectivamente. Note que se F k , é ˆ
ˆ
negativo, então a informação vem da célula da esquerda para a célula da direita, no caso contrário, a informação vem da célula da direita para a da esquerda.
8
Rio Oil & Gas Expo and Conference Conference 2012
4. Resultados Numéricos Nesta seção serão mostrados dois problemas que foram resolvidos empregando a formulação desenvolvida, os resultados obtidos serão comparados com outros já existentes na literatura. Para o caso do escoamento bifásico consideramos um problemas clássicos de “five-spot ” e comparamos os resultados com os de (DURLOFSKY, 1993).
4.1. Estudo do Campo de Pressões num Domínio Heterogêneo e Anisotrópico O domínio é um quadrado [-1,1]x[-1,1], com condições de contorno de Dirichlet, dadas pela aplicação da solução analítica do problema nos contornos. A solução analítica é dada por:
2 sin( y) cos( y) x sin( y) ; x 0 x e sin( y) ; x 0
u
(39)
As permeabilidades são dadas da seguinte forma:
1 0 0 1 2 1 1 2
para x 0
(40) para x 0
De forma que o problema pode ser escrito da maneira que se segue: u f ( x, y )
(41)
Com termo de fonte :
2 sin( y) cos( y) x sin( y) ; x 0 x 2e cos( y) ; x 0
f ( x, y )
(42)
Será necessária a definição das normas dos erros a partir das quais foram calculadas as taxas de convergência. Definamos primeiramente u como a solução analítica do problema e u como a solução numérica. A partir daí, definimos a norma do erro máximo da seguinte forma: E MAX MAX max u u
(43)
E a norma do erro RMS como:
N u u 2 E RMS i 1 N
0.5
(44)
onde N é o número de células da malha. A norma dos erros nas velocidades também deve ser definida, para que se defina a taxa de convergência das velocidades. Sendo v a velocidade calculada analiticamente e v a velocidade calculada numericamente, temos:
9
Rio Oil & Gas Expo and Conference Conference 2012
N 2 v v A i i 1 N A i i 1
0.5
A
E VEL
(45)
A
onde N A é o número de faces que existem na malha e Ai é a soma dos volumes das células que compartilham a i-ésima face do domínio. Conhecida a notação que define as malhas (N:MXP), podemos definir as taxas de convergência da seguinte forma:
E TN MXP E TN 2MX2P
R TN log 2
(46)
onde TN é o tipo de norma de erro (MAX, RMS ou VE L). Tabela 1. Resultados pelo método MPFA-DEW2 MPFA-DEW2, com α=1. N 8X8 16X16 32X32 64X64
EMÁX 0.0156 0.0046 0.0013 3.43x10^(-4)
RMÁX ------#--1.7618 1.8231 1.9222
ERMS 0.0062 0.0016 3.95x10^(-4) 9.89x10^(-5)
RRMS ----#-----1.9542 2.0181 1.9978
Nesse exemplo, foi usado o MPFA-DEW2, com pesos explícitos de tipo 2, sobre uma malha triangular estruturada. É possível observar as altas taxas de convergência, tanto para a solução como para a velocidade, o que é ainda mais positivo. Mas comparemos com alguns resultados de (CRUMPTON, et al., 1994) e (HYMAN, et al., 1997), em que também foi resolvido esse problema: Tabela 2. Taxa de convergência na velocidade. N 8X8 16X16 32X32 64X64
Taxa de Convergência na Velocidade -------------#-----------------------1.5538 1.6415 1.7681
Tabela 3. Resultado obtidos por Hyman e Crumpton para N 8X8 16X16 32X32 64X64
(HYMAN, et al., 1997), com α=1 EMÁX RMÁX ------#-------------#-------3.74x10^(-3) ------#-------9.66x10^(-4) 1.9529 2.45x10^(-4) 1.9792
α=1.
(CRUMPTON, et al., 1994), com α =1 ERMS RRMS 3.33x10^(-3) --------#-----------9.37x10^(-4) 1.8294 2.45x10^(-4) 1.9353 6.25x10^(-5) 1.9709
Além dessas altas taxas de convergência nos erros das velocidades, é possível observar que o MPFA-DEW apresenta resultados bastante semelhantes a de (CRUMPTON, et al., 1994) e (HYMAN, et al., 1997) , tanto com respeito aos erros nas soluções quanto com respeito às taxas de convergência desses erros, o que a utilidade do método diante desse tipo de domínio (heterogêneo e anisotrópico).
4.2. Escoamento Bifásico num ¼ de Cinco Poços com Razão de Viscosidade Moderadamente Adversa M = 4.0. Este exemplo, adaptado de (DURLOFSKY, 1993) é uma versão adimensionalizada do clássico problema de ¼, de cinco poços. Para este problema, as saturações residuais de água e óleo são S ro ro=S rw rw=0 e o meio poroso é considerado 10
Rio Oil & Gas Expo and Conference Conference 2012
homogêneo e isotrópico com K = I em todo o domínio. Assumimos ainda que a porosidade é constante, sendo que seu valor real não é relevante já que a utilizamos apenas para adimensionalizar o tempo. As viscosidades da água e do óleo são, respectivamente, w 0.1 e o 0.4 de modo que a razão entre as viscosidades é M =4. =4. As condições de contorno são S w,i 1.0 no poço injetor e as pressões no canto superior esquerdo e no canto inferior direito são Pse Pid 0 . Na Figura 5, observamos que, para a malha alinhada, o resultado obtido no presente trabalho apresenta um tempo de irrupção de água no poço produtor (Breakthrough) um pouco antecipado se comparado com os resultados obtidos com a malha não-estruturada. Por outro lado, os resultados obtidos com a malha transversal, tanto no presente texto quanto no trabalho de (DURLOFSKY, 1993), apresentam um pequeno atraso para a irrupção de água no poço pro dutor. 0.7 Durlofsky Alinhada
0.6
Transversal Não Estruturada Quadrilatero
0.5 o d a l u 0.4 m u c A o 0.3 e l Ó
0.2
0.1
0
0
0.1
0.2
0.3
0.4
0. 5 VPI
0.6
0.7
0. 8
0. 9
1
Figura 5. Óleo acumulado para o problema ¼ de cinco poços com razão de viscosidades M = 4.0, para uma malha com (16x16) subdivisões, este resultado foi obtido pelo MPFA-DEW2.
1.1 Durlofsky Alinhada Transversal Não Estruturada Quadrilatero
1
0.9
0.8
o 0.7 d a r e p u c 0.6 e R o e l Ó 0.5
0.4
0.3
0.2
0.1
0
0. 2
0.4
0.6
0. 8
1
VPI
Figura 6. Óleo recuperado, para o problema ¼ de cinco poços com razão de viscosidades M = 4.0, para uma malha com (16x16) subdivisões, este resultado foi obtido pelo MPFA-DEW2.
5. Conclusão Todas as aproximações feitas pelo MPFA-DEW são “linearity-preserving”, trata-se de um método robusto e que apresenta taxas de convergência de segunda ordem para a solução do campo de pressões e de primeira ordem para o fluxo 11
Rio Oil & Gas Expo and Conference Conference 2012
em domínios altamente heterogêneos e anisotrópicos. No caso do escoamento bifásico, resolvemos um problema de ¼ de cinco poços (five-spot) e logo comparamos os resultados obtidos pelo MPFA-DEW com os resultados de (DURLOFSKY, 1993) num meio isotrópico. O caso analisado foi feito usando malha triangular transversal, alinhada, não estruturada e malha quadrilateral, porém o MPFA-DEW apresentou, de maneira geral, soluções tão boas ou melhores do que as encontradas na literatura. No futuro próximo, pretendemos utilizar um método de discretização de ordem mais alta para a equação de saturação, além de avaliar o método sobre malhas distorcidas.
6. Agradecimento Agradecimento especial a Fundação de Amparo a Ciência e Tecnologia do Estado de Pernambuco (FACEPE),pelas bolsas de estudos concedidas durante o período de desenvolvimento da pesquisa .
7. Referências CRUMPTON P.J., SHAW G.J. and WARE A. F. Discretisation and Multigrid Solution of Elliptic Equations With Mixed Derivative Terms and Strongly Discontinuous Coefficients. Journal Computation Physics . - 1994. DURLOFSKY L. J. A. Triangle Based Mixed Finite Element-Finite Volume Technique for Modeling Two Phase Flow Through Porous Media. Journal of Computational Computational Physics. - 1993. - 252-266 : Vol. 105. GX LV, LG SHEN and ZJ. SHEN Numerical Methods for Energy Flux of Diffusion Equation on Unstructured Grids , Chinese Journal of Computational Physics, 2007. - 379 – 386 386 : Vol. 24(4). HIRSCH C. Numerical Computation of Internal and External Flows. Brussels- Belgium : A Wiley-lnterscience publication, 1994. HUANG and KAPPEN A Study of Cell-Centered Finite Volume Methods for Diffusion Equations. Mathematics Research Report. – Lawrence.1998. HYMAN JAMES, SHASHKOV MIKHAIL and STEINBERG STANLY The Numerical Solution of Diffusion Problems in Strongly Heterogeneous Non-isotropic Materials. Journal of computational physics. - New Mexico, 1997. Vol. 132. JIMING WU, ZHIMING GAO Linearity Preserving Nine-Point Schemes for Diffusion Equation on Distorted Quadrilateral Meshes. Journal of Computational Physics. - Bejing, PR China, 2009. - Vol. 229. JIMING WU ZHIMING GAO A Nine Point Scheme with Explicit Weights for Diffusion Equations on Distorted Meshes. Applied Numerical Mathematics. - 2010. JISHENG KOU SHUYU SUN. A New Treatment of Capillarity to Improve the Stability of IMPES Two-Phase Flow Formulation. Computational Transport Phenomena Laboratory. - 2010. KUZMIN DMITRI Linearity Prserving Flux Correction and Convergence Acceleration for Constrained Galerkin Schemes. Applied Mathematics III. - 2011. M M. BASKO J. MARUHN and TAUSCHWITZ. An Efficient Cell- Centered Diffusion Scheme for Quadrilateral Grids. Journal Computational Physics . - 2009. Q. Y. CHEN and J WAN Y.YANG and R. T. MIFFLIN Enriched Multi-point Flux Approximation for General Grids. Journal Computational Physics . - 2008. WU JIMING and GAO ZHIMING A. Nine-Point Scheme with Explicit Weights for Diffusion Equations on Distorted Meshes. Applied Numerical Mathematics. - Bejing, China : Applied Numerical Mathematics, 2011. - Vol. 61. WU JM Linearity Exact Method and the Difference Scheme for Diffusion Equation on Quadrilateral Meshes. Beijing. Institute of Applied Physics and Computational Mathematics, 2005. - Vols. 156 – 169. 169. WU JM, FU S. W. and SHEN L. J. A. Difference Scheme with High Resolution for the Numerical Solution of a Nonlinear Diffusion Equation. Journal on Numerical Methods and Computer Applications. - 2003. - Vol. 24. ZHIMING and WU JIMING A. Linearity-Preserving Cell-Centered Scheme for the Heterogeneous and Anisotropic Diffusion Equations on General Meshes. Laboratory of Computational Physics . – Bejing. 2010.
12
