PROPOSTA paula sesini blackoil

Prévia do material em texto

Simulação numérica de escoamentos tri-
fásicos em reservatórios de petróleo 
 
 
 
Proposta de Tese de Doutorado 
 
 
 
Paula Aida Sesini 
Orientador: Alvaro Luiz Gayoso de Azeredo Coutinho 
 
 
 
junho de 2004 
 1
 
1. Introdução 
 
A simulação de reservatórios de petróleo tem sido motivo de amplos estudos nos últimos anos 
devido a sua fundamental importância para a indústria petrolífera. O conhecimento do 
comportamento do reservatório durante sua exploração fornece dados determinantes como a 
localização de poços, injeção de água, tempo de produção, etc. Desta maneira é possível predizer 
a produtividade de um reservatório, desde um ponto de vista econômico. 
 
Os reservatórios são compostos por rochas porosas impregnadas com petróleo e em geral contêm 
camadas de gás natural, hidrocarbonetos e água sob grandes pressões. O petróleo pode ser 
extraído através de um processo chamado de recuperação primária, no qual, após a perfuração de 
um poço, a pressão interna impulsiona o gás, os hidrocarbonetos ou a água. No caso em que as 
pressões internas são baixas ou quando a pressão inicial diminui devido a um certo tempo de 
extração, podem ser usadas algumas técnicas de recuperação chamadas de recuperação avançada. 
Exemplos de estas técnicas são a injeção de água, a injeção de um componente químico miscível, 
etc. 
 
O comportamento de um reservatório de petróleo pode ser considerado um problema de 
escoamento multifásico em meios porosos apresentando três fases: água, óleo e gás. Diferentes 
modelos matemáticos tem sido desenvolvidos para descrever este tipo de escoamento. O modelo 
mais complexo é o chamado modelo composicional que considera o caso muito geral no qual 
existem N espécies químicas, ou componentes, as quais podem estar presentes em qualquer uma 
das três fases. Esta descrição do problema é necessária no caso de sistemas que contém óleo 
altamente volátil. 
 
A maioria dos reservatórios de petróleo contém óleo de baixa volatilidade e por este motivo, 
entre outros, o modelo mais freqüentemente usado é o chamado black-oil. Este modelo pode ser 
aplicado para descrever um sistema com óleo, composto essencialmente por metano e 
componentes pesados, usando dados extraídos de um teste de vaporização diferencial 
convencional aplicado a amostras de óleo de reservatório (Peaceman, 1977; Bear, 1972). O 
black-oil é um modelo composicional simplificado, de três componentes água, óleo e gás, que 
considera os efeitos de compressibilidade e transferência de massa entre as fases. Este modelo é 
aplicável tanto para a modelagem da recuperação primária (queda de pressão) como da 
recuperação avançada no caso da injeção de água (Trangenstein and Bell, 1989). 
 
O escoamento trifásico em reservatórios também pode ser representado por um modelo mais 
simplificado que o black-oil, o qual não considera os efeitos de compressibilidade e transferência 
de massa entre as fases. Este modelo tem sido motivo de numerosos estudos devido à 
necessidade de um melhor entendimento dos modelos matemáticos do escoamento no meio 
poroso, para otimizar os métodos de recuperação de óleo de reservatórios e de remoção de 
contaminantes de aqüíferos (Guzmán, 1995; Juanes, 2003; Abreu, 2003). 
 
O modelo black-oil, tanto como o modelo mais simplificado do escoamento trifásico, são 
descritos por um sistema de equações diferenciais parciais. Estas equações representam a 
conservação de massa de cada uma das componentes. As velocidades de cada uma das fases são 
calculadas usando a lei de Darcy estendida para o caso do escoamento multi-fásico. As 
saturações das fases estão relacionadas por uma equação de fechamento. As pressões de 
capilaridade tanto como as permeabilidades relativas são expressas por modelos dependentes das 
saturações das fases e as condições iniciais e de contorno são definidas dependendo de cada 
problema. 
 2
 
Os casos reais de escoamento multi-fásico em reservatórios se caracterizam por fatores como a 
heterogeneidade do meio, fraturas da rocha, geometrias complexas e condições de contorno 
dependentes do tempo. Devido à complexidade do problema se torna difícil o cálculo de uma 
solução analítica e, portanto surge a necessidade do uso de métodos numéricos que permitem 
obter uma solução aproximada. 
 
Os modelos matemáticos que descrevem o escoamento em meios porosos são compostos por 
equações diferenciais parciais não lineares, acopladas e dependentes do tempo. Por estes 
motivos, um problema importante em simulação de reservatórios é o desenvolvimento de 
métodos estáveis, eficientes, robustos e precisos para a solução deste tipo de sistemas de 
equações acopladas (Li et al., 2003). 
 
Existem diferentes métodos de solução em simulação de reservatórios. O método totalmente 
implícito, descrito em (Douglas et al., 1959), resolve todas as equações acopladas não-lineares 
simultaneamente. Este método é estável e pode usar grandes passos de tempo sem perder a 
estabilidade. Porém é um método que requer um alto custo computacional e grande espaço de 
memória. Uma outra alternativa consiste em resolver o sistema de equações de forma seqüencial 
(MacDonald and Coats, 1970). Isto é: resolver cada uma das equações do sistema separadamente 
e implicitamente. Este método é menos estável que o totalmente implícito mais por outro lado é 
mais eficiente devido a que requer menos espaço de memória e custo computacional. 
Recentemente Li et al., (2003) revisitaram o método seqüencial para resolver o modelo black-oil 
usando malhas não-estruturadas. Usaram o método CVFA (control volume function 
approximation) para discretizar as equações governantes. Concluíram que este método pode 
reduzir altamente o custo computacional e o tamanho da memória para simulação de 
reservatórios com malha não-estruturada. 
 
Por outro lado, Juanes (2003) apresentou uma alternativa diferente de integração no tempo. Ele 
propõe uma solução simultânea das duas equações de saturação (da água e do gás), usando o 
método multi-escala para a discretização das equações governantes. Obteve resultados 
satisfatórios para problemas unidimensionais. 
 
As equações que descrevem o escoamento em meios porosos, além de serem acopladas e não 
lineares, se caracterizam pelo transporte predominantemente advectivo, para o caso de efeito de 
capilaridade pequeno, e pela solução apresentando frentes de choques que se movem. É sabido 
que, devido a estas características do problema, os métodos numéricos clássicos aplicados ao 
problema do escoamento no meio poroso produzem soluções com oscilações espúrias no 
domínio todo e alta difusão numérica. Outro problema que apresenta o escoamento em meios 
porosos é o cálculo preciso da velocidade através da lei de Darcy. 
 
Nesta introdução não pretendemos fazer uma revisão completa da bibliografia, o qual é parte das 
atividades que serão realizadas até o final do trabalho. Portanto aqui citamos somente alguns 
trabalhos mais relevantes. 
 
Tradicionalmente a alternativa para resolver as dificuldades devidas ao transporte 
predominantemente advectivo e ao cálculo preciso da velocidade em meios porosos tem sido os 
métodos de elementos finitos mistos ou híbridos (Chavent et al., 1990; Brezzi and Fortin, 1991; 
Durlofsky, 1993), os quais permitem soluções satisfatórias e apresentam melhores soluções para 
a velocidade total que outras aproximações padrão (Durlofsky, 1994; Ewing, 1996). Alguns 
trabalhos mais recentes que estudam este método aplicado à simulação de reservatórios são os de 
Peszynska et al. (2000) e Wheeler and Peszynska (2002). 
 
 3
Outra alternativa para a simulação do problema em meios porosos tem sido os métodos que usam 
volumes finitos para a discretização das equações governantes. Exemplos destes métodos são o 
CVFE (controle volume finite element) (MacDonald and Coats, 1970; Fung et al., 1992) e o 
CVFA(control volume function approximation) (Li et al., 2003). A principal diferença entre 
estes dois métodos consiste em que, para a interpolação, o CVFA usa funções não polinomiais 
em vez de funções polinomiais usadas pelo CVFE. 
 
Nos últimos anos tem sido pesquisadas diversas alternativas para o problema da simulação de 
reservatórios, tornando mais conveniente a aplicação de métodos de elementos finitos 
estabilizados. A idéia básica destes métodos consiste na introdução de uma difusão artificial que 
age somente na direção das linhas de corrente. Estes métodos foram desenvolvidos para resolver 
o problema do transporte dominantemente advectivo em mecânica de fluidos (Hughes and 
Brooks, 1979). Exemplos de métodos de elementos finitos estabilizados são o método SUPG 
(streamline upwind Petrov-Galerkin) (Brooks and Hughes, 1982), o método GLS (Galerkin least-
squares) (Hughes et al., 1989) e o método ASGS (algebraic subgrid scale) (Hughes, 1995). Estes 
métodos são sempre acompanhados de uma técnica de captura de choque para controlar as 
instabilidades produzidas nas regiões de altos gradientes e de frentes de choque que se movem. 
 
Um exemplo de aplicação do método SUPG com operador de captura de choque e 
descontinuidades CAU ao escoamento em meios porosos é apresentado no trabalho de Coutinho 
e Alves (1999). Neste trabalho é simulado o efeito não-linear de viscous-fingering no 
escoamento miscível, obtendo resultados satisfatórios. Juanes and Patzek (2001) aplicaram o 
método de elementos finitos estabilizado multi-escala ao problema de injeção de traçadores 
(escoamento miscível) tanto como ao problema de águas subterrâneas (escoamento imiscível 
bifásico) para casos unidimensionais. Usaram uma difusão artificial no nível da sub-escala para a 
captura de choques, obtendo resultados razoáveis. Posteriormente Coutinho et al. (2003) 
realizaram uma comparação de esquemas de estabilização para a simulação por elementos finitos 
de escoamentos imiscíveis bifásicos em meios porosos. Estudaram exemplos bidimensionais 
obtendo bons resultados. Recentemente Juanes (2003) tem aplicado o método estabilizado multi-
escala ao problema do escoamento trifásico incompressível, sem considerar o efeito de 
transferência de massa entre as fases. Estudou casos unidimensionais, obtendo boas soluções. 
 
Neste trabalho usamos o método de elementos finitos para resolver o problema do escoamento 
multifásico em meios porosos no caso do modelo black-oil, tanto como no caso do modelo 
trifásico simplificado. Escolhemos este método devido principalmente a sua forte base 
matemática e à flexibilidade que apresenta para tratar geometrias complexas. 
 
O objetivo deste trabalho consiste em desenvolver um simulador de reservatórios de petróleo que 
trate problemas de escoamento trifásico. A primeira etapa consistirá em implementar um modelo 
simples que considere fluidos incompressíveis e que não ocorra transferência de massa entre as 
fases. Uma vez compreendidas e estudadas as dificuldades desta primeira etapa 
implementaremos o modelo black-oil. 
 
Em ambos os casos aplicaremos o método de Galerkin clássico para resolver a equação da 
pressão, enquanto que para resolver as equações de saturação, as quais apresentam transporte 
predominantemente advectivo, usaremos o método de elementos finitos estabilizados SUPG com 
operador de captura de choque e descontinuidades CAU. 
 
Quanto à integração no tempo propomos dois esquemas diferentes, com o intuito de comparar a 
performance computacional, estabilidade e precisão dos mesmos. Um destes esquemas é um 
método seqüencial muito similar ao apresentado em (Li et al., 2003). O segundo é um esquema 
que resolve inicialmente a equação da pressão, como no caso do método seqüencial, e em 
 4
seguida resolve as duas equações da saturação simultaneamente, de maneira similar ao esquema 
apresentado em (Juanes, 2003). 
Implementaremos estes esquemas de integração com base nos estúdios realizados sobre leis de 
conservação, apresentados nos capítulos 2 e 3 deste trabalho. Desta maneira, o esquema de 
integração das equações de saturação de forma simultânea pode ser implementado usando uma 
técnica muito similar àquela usada para a solução das equações de Euler (ver capítulo 2), devido 
à semelhança na estrutura das equações dos dois problemas. Da mesma forma, o esquema de 
integração sequencial pode ser implementado como uma extensão da implementação do 
escoamento bifásico (ver capítulo 3). 
 
Esta proposta de trabalho está organizada da seguinte forma. A seção 2 apresenta as equações 
governantes do modelo black-oil. Com objetivo de realizar comparações com outros trabalhos da 
literatura descrevemos dois sistemas diferentes de equações equivalentes. Um destes esquemas 
corresponde a uma formulação simultânea das equações de saturação enquanto que o outro 
sistema consiste em uma formulação seqüencial destas equações. Em seguida são apresentados 
também dois sistemas diferentes de equações equivalentes correspondentes a um modelo mais 
simples, o qual considera fluidos incompressíveis e que não ocorre transferência de massa entre 
as fases. 
 
Na seção 3 apresentamos a formulação semi-discreta de elementos finitos aplicada a este último 
modelo simplificado. Aplica-se esta formulação aos dois esquemas de equações equivalentes 
tratados. No caso do esquema com formulação simultânea das equações de saturação propomos 
uma formulação análoga àquela usada para resolver as equações de Euler, enquanto que no caso 
do esquema com formulação seqüencial das equações de saturação a aplicação da formulação de 
elementos finitos consiste numa extensão da técnica usada para resolver o caso do escoamento 
bifásico. 
 
Na seção 4 apresentamos os esquemas de integração propostos para resolver os sistemas de 
equações resultantes da aplicação da formulação de elementos finitos. Em ambos casos, 
propomos um algoritmo multi-corretor bloco-iterativo. A diferença consiste em que no caso do 
esquema de integração das equações de saturação de forma simultânea estas equações são 
resolvidas num único bloco usando uma técnica muito similar àquela usada para a solução das 
equações de Euler (ver capítulo 2). Por outro lado no caso do esquema de integração seqüencial 
cada equação de saturação é resolvida em um bloco separado como uma extensão da 
implementação do escoamento bifásico (ver capítulo 3). 
 
Finalmente a seção 5 apresenta o cronograma de atividades programadas para realizar os 
objetivos do trabalho proposto. 
 
 
 
 
 
 
2. Equações governantes 
 
Nesta seção apresentamos tanto as equações governantes do modelo black-oil como as do 
modelo trifásico simplificado. No modelo black-oil se assume que não existe 
transferência de massa entre a fase água e as outras duas fases, óleo e gás. O componente 
gás pode ser encontrado tanto na fase gás como dissolvido na fase óleo, enquanto que o 
componente óleo existe somente na fase óleo. O sistema de hidrocarbonetos é 
considerado de duas componentes: óleo e gás. A primeira é o líquido residual à pressão 
atmosférica depois de uma vaporização diferencial, entanto que a segunda é o fluido 
restante no meio poroso. 
 
O modelo black-oil é descrito pelas equações de conservação de massa de cada uma das 
componentes. Estas equações, em um domínio 2RΩ∈ com um contorno Γ em um 
intervalo de tempo [0, T] podem ser escritas da seguinte maneira (Peaceman,1977): 
 
)()( w
w
WS
ww
w
WS s
Bt
q
B
ρφρ ∂
∂=+⋅∇− v em Ω × [0, T] (1) 
)()( o
o
OS
oo
o
OS s
Bt
q
B
ρφρ ∂
∂=+⋅∇− v em Ω × [0, T] (2) 
⎥⎥⎦
⎤
⎢⎢⎣
⎡ +∂
∂=++⋅∇− )()( o
o
GSso
g
g
GS
go
o
GSso
g
g
GS s
B
Rs
Bt
q
B
R
B
ρρφρρ vv em Ω × [0, T] (3) 
 
onde os subíndicesw, o e g (em minúsculas) se referem às fases da água, óleo e gás 
enquanto que os subíndices W, O e G (em maiúsculas) se referem às componentes água, 
óleo e gás respectivamente. 
 
Nestas equações φ é a porosidade, ρIS é a densidade da componente I em condições 
padrão1 enquanto que Bi,vi, qi e si são o fator de volume de formação, a velocidade, a taxa 
de injeção de massa nos poços e a saturação correspondentes à fase i respectivamente. 
 
O fator volume de formação, Bi, expressa a relação entre o volume da fase i em condições 
do reservatório e o volume do componente I em condições padrão. Portanto: 
 
ISii VVTpB /),( = (
 
4) 
 solubilidade do gás, Rso, é definida como o volume de gás (medido em condições 
 
A
padrão) dissolvido para uma dada pressão e temperatura do reservatório em uma unidade 
de volume de óleo (medido em condições padrão). Isto é: 
 
 
1 Devido a que os reservatórios podem ter diferentes condições, todos os volumes são referidos a condições 
padrão as quais são: PSC=1 atm=14.5 psi, TSC=60oF = 15oC. 
 
 5
OSGSso VVTpR /),( = (5) 
 
A velocidade vi de cada fase é definida pela lei de Darcy como segue: 
 
( zKv ∇−∇−= gpk ii
i
ri
i ρµ ) i = w, o, g. (6) 
 
onde K é o tensor de permeabilidade absoluta dependente somente da posição e kri, µi, pi 
e ρi são a permeabilidade relativa, a viscosidade, a pressão e a densidade correspondentes 
à fase i respectivamente; g é a constante gravitacional e z∇ é o vetor unitário na direção 
da profundidade. 
 
As saturações da água, do óleo e do gás satisfazem a seguinte condição: 
 
1=++ gow sss (
 
7) 
s pressões de fase estão relacionadas com as pressões de capilaridade pcow e pcgo da 
 (8) 
 
anto as permeabilidades relativas como as pressões de capilaridade apresentam uma 
ubstituindo as equações de fechamento (5), (6) e (7) nas equações de conservação (1), 
este trabalho consideramos condições de contorno de não penetrabilidade do fluxo. 
i , i = w, o, g, 
A
seguinte maneira: 
 
wocow ppp −=
 (9) ogcgo ppp −=
T
dependência não linear devido ao fato que é aceito, como fato empírico, que estas 
propriedades são unicamente funções das saturações das fases. Para representar estas 
grandezas são usados modelos baseados em medições de laboratório. Curvas típicas no 
caso do escoamento bifásico são mostradas em (Peaceman,1977). Por outro lado, no caso 
do escoamento trifásico estas curvas são mais difíceis de obter e requerem interação entre 
experimentos de laboratório e de campo. 
 
S
(2) e (3) se obtém um sistema de três equações diferenciais parciais com três incógnitas, 
podendo se escolher diferentes conjuntos de variáveis independentes. Diversos estudos 
(Wu and Forsyth, 2001) mostram que a escolha das variáveis primitivas de um modelo 
multifásico tem um impacto significativo sobre a performance computacional. A seleção 
ótima destas variáveis depende, na maior parte das vezes, do problema particular que está 
sendo tratado. Portanto é conveniente realizar a escolha das variáveis levando em 
consideração: a eficiência computacional, a robustez e a simplicidade no cálculo de 
outras variáveis secundárias e da atualização das equações linearizadas. 
 
N
 
0. =n Γ∈x (10) 
 
v
 6
As condições iniciais dependem do estado do reservatório, o qual pode se encontrar em 
, (11) 
 
= , (12) 
 
0 xws= , (13) 
No estado satu onsi
orrespondentes são: 
 (14) 
O sistema de equações d
stema de equações equivalente, o qual é composto por uma equação para a pressão e 
as variáveis do 
roblema escolhidas são po, sw e sg: 
estado não-saturado ou saturado. O reservatório está em estado não-saturado quando todo 
o gás se dissolve na fase óleo e não existe fase gás. Ao contrário, o reservatório está em 
estado saturado quando as três fases água, óleo e gás existem simultaneamente. Para o 
estado não-saturado consideramos como variáveis do problema po, sw e pb, onde pb é a 
pressão do ponto de bolha2. As condições iniciais para este caso são: 
 
)()0,( 0 xx oo pp = Γ∈x
)()0,( 0 xx bb p Γ∈x 
(xws
p
)()0, Γ∈x
rado c deramos como incógnitas po, sw e so e as condições iniciais 
 
c
 
)()0,( 0 xx oo pp = , x Γ∈
 (15) 
 
)()0,( 0 xx ww ss = , Γ∈x
, (16) 
 
)()0,( 0 xx oo ss = Γ∈x
e conservação (1), (2) e (3) pode ser transformado em um 
 
si
duas equações para as saturações. Neste trabalho apresentamos diferentes esquemas para 
o sistema equivalente de equações governantes com o motivo de poder realizar 
comparações com métodos desenvolvidos em outros trabalhos. Estes são: 
 
• Esquema com formulação simultânea das equações de saturação no qual 
p
 
- Equação da pressão: 
( )soo
o
g
g
gg
o
oo
w
ww
o
tT t
cQ ∂−=⋅∇ vT φ
R
B
B
B
B
B
B
B
B
p
∇⋅−⎟⎟⎠
⎞
⎜⎜⎝
⎛∇⋅−⎟⎟⎠
⎞
⎜⎜⎝
⎛∇⋅−⎟⎟⎠
⎞
⎜⎜⎝
⎛∇⋅−
∂
vvvv 111
 (17) 
 
⎟⎟⎠
⎞
⎜⎜⎝
⎛ +++−=
o
so
o
go
g
g
g
g
o
o
o
o
w
w
w
w
t dp
dR
B
Bs
dp
dB
B
s
dp
dB
B
s
dp
dB
B
sc (18) 
 
2 Pressão de ponto de bolha: é a pressão mínima requerida para o gás dissolver no óleo. 
 
 7
GS
osog
GS
gg
OS
oo
WS
ww
T
qRBqBqBqBQ ρρρρ −++= 
 
 - Equação da saturação da água: 
 
 (19) 
 
[ −+⋅∇+∂ w ffs ((Kv ρρλφ ]
⎟⎟⎠
⎞
⎜⎜⎝
⎛∇⋅−⎟⎟⎠
⎞
⎜⎜⎝
⎛
∂
∂−=
⎥⎥⎦
⎤
⎢⎢⎣
⎡ ∇+∇+−⋅∇+
∇−+∂
w
ww
w
ww
WS
ww
g
g
cgo
gw
w
cwo
gow
gwgwooww
B
B
Bt
sBqB
s
ds
dp
s
ds
dpf
g
t
11
))((
))()T
v
K
z
φρ
λλλ
ρρλ
 (20) 
 
 
 - Equação da saturação do gás: 
 
[ +⋅∇+∂+∂ goo sogg ftBt vTφφ ]
⎟⎟⎠
⎞
⎜⎜⎝
⎛∇⋅−⎟⎟⎠
⎞
⎜⎜⎝
⎛∇⋅−⎟⎟⎠
⎞
⎜⎜⎝
⎛
∂
∂−⎟⎟⎠
⎞
⎜⎜⎝
⎛
∂
∂−=
⎥⎥⎦
⎤
⎢⎢⎣
⎡ ∇+−∇⋅∇+
⋅∇+∇−+−∂∂
o
so
og
g
gg
o
so
og
g
gg
GS
gg
g
g
cgo
oww
w
cwo
wg
o
o
sog
ogowgwg
B
RB
B
B
B
R
t
sB
Bt
sB
qB
s
ds
dp
s
ds
dpf
B
RB
gfs
RBs
vv
K
vzK
11
))((
))()((
φφρ
λλλ
ρρλρρλ
(21) 
 
As equações (18) e (19) podem também ser apresentadas usando uma formulação 
simultânea, obtida de maneira análoga ao caso do escoamentotrifásico incompressível 
uanes, 2003). (J
 
- Formulação simultânea das equações de saturação (18) e (19): 
 
QUDAUDFCFM =∇⋅∇+∇⋅∇+⋅∇+⋅∇+∂ )()( 2121t 
U∂ (22) 
 
⎥⎥
⎥
⎦
⎤
⎢⎢
⎢
⎣
⎡
−−=
o
gso
o
gso
B
BR
B
BR φφ
φ
)1(
0
M (23) 
 
 (24) ⎥⎦
⎤⎢⎣
⎡=
gB0
00
C
 
 8
⎥ (25) ⎦
⎤⎢⎣
⎡=
gg BB
00
A
 
 
⎥⎦
⎤⎢⎣
⎡=
g
w
s
s
U
⎦∇−+−+ zK gf ogowgwg ))()((T ρρλρρλ
 (26) 
⎤⎢⎣
⎡ ∇−+−+=
v
zKv
F
f
gff
g
gwgwooww ))()((T
1
ρρλρρλ
 (27) 
 
⎥
⎥⎦
∇z (28) ⎥
⎤
⎢⎢⎣
⎡
−+−+= KvF g
B
Rf
B
fR
gogwoo
o
so
o
o
oso ))()((
0
T2 ρρλρρλ
 
⎥⎥
⎥⎥
⎥
⎦
⎤
⎢⎢
⎢⎢
⎢
⎣
⎡
+−
+−
=
g
cgo
owg
w
cwo
wg
g
cgo
gw
w
cwo
gow
ds
dp
f
ds
dp
f
ds
dp
f
ds
dp
f
))((
))((
1
λλλ
λλλ
KK
KK
D (29) 
 
⎥⎥⎦
⎤
⎢⎢⎣
⎡
=
g
cgo
go
o
so
w
cwo
oo
o
so
ds
dp
f
B
R
ds
dpf
B
R λλ KKD
00
2
⎥⎥⎞⎛∂∂
⎠⎝⎠ w
RsB
B
1φ (31) 
 (30) 
⎥⎥
dos para o 
⎥
⎦
⎤
⎢⎢
⎢⎢
⎢
⎣
⎡
⎟⎟⎠⎜
⎜
⎝∂
−∂−⎟⎟⎠
⎞
⎜⎜⎝
⎛∇⋅−⎟⎟⎠
⎞
⎜⎜⎝
⎛
∂
∂−
⎟⎟
⎞
⎜⎜
⎛∇⋅−⎟⎟
⎞
⎜⎜⎝
⎛
∂
∂−
=
o
osog
so
o
og
g
gg
g
gg
GS
gg
ww
w
ww
WS
ww
Bt
sRB
tBB
B
Bt
sB
qB
B
Bt
sBqB
11
11
φφρ
φρ
v
v
Q
 
• Esquema com formulação seqüencial. Neste caso são considerados dois esta
reservatório de petróleo: 
 
 Estado saturado: existem as três fases água, óleo e gás. As variáveis escolhidas do 
q
o
problema são po, sw e so e a equação da pressão é a mesma equação (15) apresentada para 
o caso da formulação simultânea: 
 
- E uação da saturação da água: 
 
 9
[ ]
⎟⎟⎠
⎞
⎜⎜⎝
⎛∇⋅−⎟⎟⎠
⎞
⎜⎜⎝
⎛
∂
∂−=
⎥⎥⎦
⎤
⎢⎢⎣
⎡ ∇+∇+−⋅∇+
∇−+−+⋅∇+∂
∂
w
ww
w
ww
WS
ww
g
g
cgo
gw
w
cwo
gow
gwgwooww
w
B
B
Bt
sBqB
s
ds
dp
s
ds
dpf
gff
t
s
11
))((
))()((T
v
K
zKv
φρ
λλλ
ρρλρρλφ
 (
 
32) 
Equação da saturação do óleo: - 
 
[ ]
⎟⎟⎠
⎞
⎜⎜⎝
⎛∇⋅−⎟⎟⎠
⎞
⎜⎜⎝
⎛
∂
∂−=
⎥⎦
⎤⎢⎣
⎡ ∇+∇⋅∇+
∇−+−+⋅∇+∂
∂
o
oo
o
oo
OS
oo
o
o
cgo
gw
w
cwo
oo
gogwoooo
o
B
B
Bt
sBqB
s
ds
dp
s
ds
dpf
gff
t
s
11
)(
))()((T
v
K
zKv
φρ
λλ
ρρλρρλφ
 (33) 
 
stado não-saturado: existem somente duas fases água e óleo entanto que o 
omponente gás só existe dissolvido na fase óleo. As variáveis escolhidas do problema 
 
o E
c
são po, sw e pb aonde pb é a pressão do ponto de bolha. Neste caso a equação da pressão e 
a equação da saturação da água são iguais ao caso do estado saturado. A equação da 
pressão do ponto de bolha é obtida a partir da equação de conservação do componente 
gás nas condições do estado não-saturado: 
 
- Equação da pressão do ponto de bolha: 
 
[ ]
⎟⎟⎠
⎞
⎜⎜⎝
⎛∇⋅−∇⋅−⎥⎦
⎤⎢⎣
⎡ ∇⋅∇−
∇−+⋅∇−∂
∂−=
∂
∂⎟⎞⎜⎛∂ bbsoo pdpdRs 1φ ⎟⎠⎜⎝+∂
o
ooso
so
o
w
w
cwo
oo
woooo
o
ob
oo
bso
B
BR
R
s
ds
dpf
gff
t
s
tBdp
Bs
tdp
1)(
))((T
vvK
zKv
λ
ρρλφ
φ
 (34) 
 
Considerando fluidos incompressíveis e que não ocorre transferência de massa entre as 
ses o sistema de equações governantes pode ser simplificado da seguinte maneira: 
 (35) 
R
fa
 
• Esquema com formulação simultânea das equações de saturação: 
 
- Equação da pressão: 
 
TQ=⋅∇ Tv 
 10
 
ggoowwT qqqQ ρρρ ++= (36) 
 
Formulação simultânea das equações de saturação: - 
 
QUDFU =∇⋅∇+⋅∇+∂
∂ )(
t
 (37) 
 (38) 
 
⎥⎦
⎤⎢⎣
⎡=
g
w
s
s
U
⎥⎥
⎥⎥
⎦
⎤
⎢⎢
⎢⎢
⎣
⎡
∇−+−+
∇−+−+
=
zKv
zKv
F
g
ff
gff
ogowgw
gg
gwgwoo
ww
))()((
))()((
T
T
ρρλρρλφφ
ρρλρρλφφ (39) 
 
⎥⎥
⎥⎥
⎥
⎦
⎤
⎢⎢
⎢⎢
⎢
⎣
⎡
+
+
=
g
cgo
ow
g
w
cwo
w
g
g
cgo
g
w
w
cwo
go
w
ds
dpf
-
ds
dpf
ds
dpf
ds
dpf
-
)(
)(
λλφλφ
λφλλφ
KK
KK
D
⎡
 (40) 
 
wwqρ (41) 
 
• Esquema com formulação seqüencial. Neste caso a equação da pressão é a mesma 
quação (30) apresentada para o caso da formulação simultânea: 
⎥⎦⎢⎣
=
gg qρQ 
⎤
e
 
- Equação da saturação da água: 
 
 
[ oowww ffs ρλφ −+⋅∇+∂ ((T Kv ]
ww
g
g
cgo
gw
w
cwo
gow
gwgw
q
s
ds
dp
s
ds
dpf
g
t
ρ
λλλ
ρρλρ
=
⎥⎥⎦
⎤
⎢⎢⎣
⎡ ∇+∇+−⋅∇+
∇−+∂
))((
))()
K
z
 (42) 
 
- Equação da saturação do gás: 
 11
 
[ ]
gg
g
g
cgo
oww
w
cwo
wg
ogowgwgg
g
q
s
ds
dp
s
ds
dpf
gff
t
s
ρ
λλλ
ρρλρρλφ
=
⎥⎥⎦
⎤
⎢⎢⎣
⎡ ∇+−∇⋅∇+
∇−+−+⋅∇+∂
∂
))(
))()((T
K
zKv
 (43) 
 
Na seção seguinte veremos as formulações semidiscretas de elementos finitos aplicadas 
aos esquemas de equações governantes do modelo simplificado do escoamento trifásico. 
Estes casos serão os primeiros a ser implementados conforme é apresentado no 
cronograma de atividades. 
 12
3. Formulação semidiscreta de elementos finitos 
 
Nesta seção apresentamos a formulação variacional semidiscreta de elementos finitos 
aplicada às equações governantes que descrevem o escoamento trifásico. Porém neste 
documento mostramos somente os casos de escoamento trifásico incompressível: o caso 
com formulação simultânea e o caso com formulação seqüencial das equações de 
saturação da água e do gás. O desenvolvimento das restantes formulações variacionais 
correspondentes aos outros esquemas de equações governantes do modelo black-oil será 
realizado como parte das atividades programadas para o período seguinte. 
 
A formulação semidiscreta de elementos finitos consiste em uma discretização de 
elementos finitos no espaço seguida de uma discretização de diferenças finitas no tempo. 
 
3.1 Esquema com formulação simultânea das equações de saturação 
 
Consideramos um domínio computacional Ω dividido em nel elementos Ωe, e=1,2...,nel, 
onde Ω= Ωnele 1=∪ e e Ωi∩Ωj=ø. Os espaços das funções de interpolação para a pressão ph, 
para as saturações da água e do gás Sh e os espaços das funções de pesowh e Vh 
correspondentes a estas aproximações são respectivamente definidos como segue: 
 
ptpPpppp hee
hhhhh =Ω∈ΩΩ∈= )()],([/)],([/{ 1H in (44) }Γd
 
kk
e
e
hh gP =Ω∈ΩΩ∈= eUUHUUS .,)]([/,)]([/{ h21h21hh in (45) }Γgk
 
0)],([/)],([/{ h1hhh =Ω∈ΩΩ∈= wPwwww eehh H in }Γ (46) 
 
0.,)]([/,)]([/{ h21h21hh =Ω∈ΩΩ∈= keehh P eWWHWWV in }gkΓ (47) 
 
onde Hh(Ω) é um espaço de funções de dimensão finita sobre Ω, P1(Ωe) representa 
polinômios de primeira ordem em Ωe, H1h(Ω) é um espaço de funções de dimensão finita 
sobre Ω e P1(Ωe) representa polinômios de primeira ordem em Ωe. Γgk é o contorno de Ω 
com condições prescritas de Dirichlet. 
 
Considerando uma discretização padrão de Ω em elementos finitos a formulação de 
Galerkin para a equação da pressão é definida como, 
 
Ω=Ω⋅∇ ∫∫ ΩΩ dQwdw ThTh v (48) 
 
A equação da saturação (35) alternativamente pode ser escrita da seguinte maneira: 
 
0, =++ yyxxt ,, UAUAU em Ω × [0, T] (49) 
 
onde as matrizes jacobianas Ai, para i=x,y são definidas como se segue: 
 12
 
U
FA ∂
∂= ii (50) 
 
e o vetor Fi é definido da seguinte maneira: 
 
⎥⎥
⎥⎥
⎦
⎤
⎢⎢
⎢⎢
⎣
⎡
∇+∇−+−+
∇+∇−+−+
=
)())()((
)())()((
T
T
yiyxixogowgw
g
i
g
yiyxixgwgwoo
w
i
w
i
kkg
ff
kkgff
zzv
zzv
F
ρρλρρλφφ
ρρλρρλφφ (51) 
 
onde kix, kiy são os coeficientes do tensor de permeabilidade absoluta K, e são 
as componentes do vetor unitário na direção da profundidade e v
xz∇ yz∇
Tx e vTy são as 
componentes da velocidade total vT. 
 
Desta maneira, a formulação SUPG aplicada à equação da saturação (47) adicionando um 
operador de captura de choque (LeBeau & Tezduyar, 1991) é definida da seguinte forma, 
 
0)(]).[()(
).(
11
=Ω∂
∂
∂
∂∑ ∫+Ω∂
∂+∂
∂
∂
∂∑ ∫
+Ω∂
∂+∂
∂∫
= Ω= Ω
Ω
d
xx
d
xtx
d
xt
i
h
i
hnel
ei
h
h
i
h
k
hnel
e
th
k
i
h
h
i
h
h
ee
UWUAUWτA
UAUW
δ
 (52) 
 
Na equação (50) a primeira integral é a formulação de Galerkin, a primeira somatória de 
integrais no nível de elemento é o termo de estabilização SUPG e a segunda somatória 
representa o termo de captura de choque adicionado à formulação variacional para evitar 
a presença de oscilações espúrias em torno das regiões de choque. 
 
A matriz de estabilização τ pode ser definida por matrizes diagonais. Esta forma de 
estabilização foi inicialmente introduzida por Hughes & Tezduyar (1984) e mais tarde 
melhorada por Aliabadi & Tezduyar (1995). A matriz τ = τI depende do parâmetro τ que 
é definido como, 
 
)](,0max[ δττζττ −+= al (53) 
 
al CFL
τατ )21(3
2
+= (54) 
 
λτ 2
h
a = (55) 
 
2λ
δτδ = (56) 
 13
 
onde τl é o parâmetro de estabilização correspondente aos termos dependentes do tempo, 
τa o parâmetro de estabilização correspondente aos termos advectivos, τδ o parâmetro de 
estabilização para descontar os efeitos do operador de captura de choque e h um 
parâmetro dependente da malha definido como 1/2, onde  é a área do elemento. O 
parâmetro ζ é um coeficiente usado no algoritmo de integração no tempo e CFL é o 
número de Courant-Friedrichs-Lewy, definidos como, 
 
CFL
CFL
α
αζ
21
2
+= (57) 
 
h
tCFL ∆= λ (58) 
 
onde λ é o maior autovalor das matrizes Ai a ser determinado, α é um parâmetro que 
controla a estabilidade e precisão do algoritmo de integração no tempo e ∆t é o passo de 
tempo. Neste trabalho adotaremos α=0.5. O parâmetro de captura de choque δ pode ser 
definido de forma similar a como é apresentado em Le Beau et al., (1993) (δCD) ou como 
em Almeida & Galeão (1996) (δCAU) da seguinte maneira: 
 
2
1
2
2
⎥⎥
⎥⎥
⎥
⎦
⎤
⎢⎢
⎢⎢
⎢
⎣
⎡
∇
∂
∂
=
h
h
i
CD
i
U
UA
ξ
δ 
(59) 
 
hCAU U
UR
ξ
δ ∇=
)(
 (60) 
w
h
a
h
i
t
QUDvUUR +∇+⋅∇+∂
∂= )()( φ (61) 
y
y
x
x
y
y
x
x hhhhh
∂
∂
∂
∂+∂
∂
∂
∂+∂
∂
∂
∂+∂
∂
∂
∂=∇ UUUUU
2211 ξξξξξ (62) 
 
se ||∇ ξUh || ≠ 0, ao contrário δ=0. As componentes ∂xj/∂ξl são os termos da matriz de 
transformação entre as coordenadas físicas e as coordenadas locais dos elementos (ξl). 
 
3.2 Esquema com formulação seqüencial das equações de saturação 
 
Consideramos um domínio Ω dividido em nel elementos, Ωe, e=1, 2,…,nel, onde Ω= 
Ω
nel
e 1=∪
e e Ωi∩Ωj=ø. Os espaços das funções de interpolação para a pressão ph, a saturação da 
 14
água swh, a saturação do gás sgh e o espaço das funções de peso wh são respectivamente 
definidas por: 
 
ptpPpppp hee
hhhhh =Ω∈ΩΩ∈= )()],([/)],([/{ 11H in (63) }Γd
 
h
wi
h
w
e
e
h
w
hh
ww
h
w stsPssss =Ω∈ΩΩ∈= )()],([/)],([/{ 11h H in (64) }Γ i
 
h
gi
h
g
e
e
h
g
hh
g
h
g
h
g stsPssss =Ω∈ΩΩ∈= )()],([/)],([/{ 11H in (65) }Γi
 
0)],([/)],([/{ h1h1hh =Ω∈ΩΩ∈= wPwwww eehh H in }Γ (66) 
onde H1h(Ω) é um espaço de funções de dimensão finita sobre Ω, P1(Ωe) representa 
polinômios de primeira ordem em Ωe. Considerando uma discretização padrão de Ω em 
elementos finitos, a formulação de Galerkin para a equação da pressão, da mesma forma 
que no caso da formulação simultânea, é escrita como, 
Ω=Ω⋅∇ ∫∫ ΩΩ dQwdw ThTh v (67) 
A formulação SUPG para as equações de saturação de cada fase (água e gás) é: 
0)(
),()),((
nel
1e
nel
1e
*
=Ω∇∇+
Ω+Ω
∑∫
∑∫∫
= Ω
= ΩΩ
dsws
dsLwLdsLw
h
i
hh
i
h
ai
h
i
h
i
h
ii
h
ai
h
i
h
i
h
e
e
δ
τ vv
 i=w,g. (68) 
O operador diferencial para cada fase é definido por: ),( vhaihisL
 
i
h
gig
h
wiw
h
ai
h
ih
ai
h
ii Qsst
ssL +∇+∇+⋅∇+∂
∂= )(),( DDvv φ (69) 
 
onde as velocidades aparentes e são definidas respectivamente por: awv agv
 
D))()((T ∇−+−+= gff gwgwoowwaw ρρλρρλKvv (70) 
 
D))()((T ∇−+−+= gff ogowgwggag ρρλρρλKvv (71) 
 
as matrizes de difusão , , e são definidasrespectivamente por: wwD wgD gwD ggD
 
w
w
cwo
gowww sds
dpf ∇+−= )( λλKD (72) 
g
g
cgo
gwwg sds
dp
f ∇= λKD (73) 
w
w
cwo
wggw sds
dpf ∇= λKD (74) 
 15
g
g
cgo
owggg sds
dp
f ∇+−= )( λλKD (75) 
e o termo fonte para cada fase é definido por: iQ
 
iii qQ ρ−= (76) 
 
O operador diferencial para cada fase é definido como segue: )(* hi wL
 
hh
ai
h
i wwL ∇⋅= v* (77) 
 
De forma semelhante ao caso da formulação simultânea, na equação (65) a primeira 
integral é a formulação de Galerkin, a primeira somatória de integrais no nível de 
elemento é o termo de estabilização SUPG e a segunda somatória representa o termo de 
captura de choque adicionado à formulação variacional para evitar a presença de 
oscilações espúrias em torno das regiões de choque. 
 
O parâmetro de estabilização para cada fase iτ é definido por: 
 
1
2
)(
−
⎟⎟⎠
⎞
⎜⎜⎝
⎛=
h
s
c
h
i
h
ai
i
vτ (78) 
 
onde h é um parâmetro dependente da malha e c2 é um coeficiente constante. 
 
Neste trabalho adotaremos o parâmetro de captura de choque CAU introduzido por 
Galeão e Carmo4, o qual apresenta a seguinte forma: 
 
h
i
h
i
ii s
sR
h ∇=
)(αδ (79) 
 
onde i=w,g, αi é um parâmetro que depende do número de Peclet. é o resíduo para 
cada fase definido por: 
)( hisR
 
w
h
gig
h
wiw
h
ai
h
ih
i Qsst
ssR +∇+∇+⋅∇+∂
∂= )()( DDvφ (80) 
 
É importante notar que o parâmetro de captura de choque é nulo quando his∇ é zero. 
Na seção seguinte veremos diferentes esquemas propostos para a integração no tempo 
dos sistemas de equações obtidos da formulação de elementos finitos. Nesta proposta 
apresentamos somente os casos correspondentes ao modelo incompressível e sem 
transferência de massa entre as fases. 
 16
4. Esquemas de integração no tempo 
 
Como mostrado na seção 2 o objetivo é resolver diferentes esquemas de equações 
governantes. Nesta proposta descrevemos os esquemas de integração no tempo que 
aplicaremos aos diferentes sistemas de equações governantes do caso de fluidos 
incompressíveis. Estes esquemas serão os primeiros a serem implementados, enquanto os 
restantes serão tratados como parte das atividades seguintes à apresentação desta 
proposta. 
 
• Esquema de integração com solução simultânea das equações de saturação: 
 
As equações algébricas obtidas da formulação de elementos finitos são 
aproximadas pela regra trapezoidal generalizada (Hughes, 1987). Para resolver os 
sistemas de equações resultantes usaremos um algoritmo multi-corretor bloco-
iterativo composto por três blocos. 
Inicialmente se resolve a equação da pressão. No segundo bloco se calcula o 
campo de velocidades usando uma técnica de pós-processamento. Em seguida se 
resolve o sistema de equações da saturação usando um algoritmo Preditor-
Multicorretor muito parecido ao apresentado em Aliabadi and Tezduyar (1995). 
Uma vez calculadas a pressão e as saturações se atualizam as variáveis e este 
processo iterativo continua até se satisfazer algum critério de convergência. 
 
• Esquema de integração seqüencial: 
 
 Da mesma maneira que no caso do esquema anterior usaremos a regra trapezoidal 
generalizada (Hughes, 1987) para a discretização no tempo das equações. Como 
resultado desta aproximação se obtém o seguinte algoritmo multi-corretor bloco-
iterativo: 
 
• Bloco 1: Resolve a equação da pressão 
• Bloco 2: Calcula o campo de velocidades 
• Bloco 3: Resolve a equação da saturação da água 
• Bloco 4: Resolve a equação da saturação do gás 
• Atualiza as variáveis 
 
O processo iterativo continua até se satisfazer algum critério de convergência. 
 
Na seção seguinte apresentamos o cronograma de atividades que serão realizadas até o 
final do período para cumprir com os objetivos desta proposta de trabalho. 
 17
5. Cronograma de atividades 
 
Nesta proposta temos apresentado as equações governantes do modelo black-oil para 
simulação de reservatórios de petróleo. Temos descrito diferentes esquemas de 
equações governantes equivalentes com objetivo de realizar comparações da precisão, 
estabilidade e performance computacional dos mesmos. Para o caso do modelo black-
oil descrevemos uma formulação simultânea e uma formulação seqüencial. A sua vez, 
esta última foi decomposta em dois casos: estado saturado e estado não saturado do 
reservatório. De mesma forma para o caso do modelo trifásico incompressível sem 
transferência de massa entre as fases apresentamos uma formulação simultânea e uma 
formulação seqüencial. 
Temos aplicado somente a formulação de elementos finitos aos esquemas de 
equações governantes do modelo incompressível sem transferência de massa entre as 
fases. Da mesma forma temos apresentado os esquemas de integração no tempo para 
o caso deste modelo mais simples. 
Para a realização dos objetivos desta proposta de trabalho é necessário realizar uma 
série de atividades, que serão realizadas até o final do período. As mesmas são 
enunciadas na lista seguinte: 
 
• Aplicação da formulação de elementos finitos aos esquemas de equações 
equivalentes propostos para descrever o modelo black-oil – julho 2004. 
 
• Definir os esquemas de integração no tempo para todos os esquemas de 
equações equivalentes do modelo black-oil – julho 2004. 
 
• Implementação do modelo trifásico incompressível como extensão do código 
do miscível (descrito no capítulo 3-3.1.2) com resolução simultânea das 
equações de saturação – novembro 2004. 
 
• Implementação do modelo trifásico incompressível como extensão do código 
do miscível (descrito no capítulo 3-3.1.2) com resolução seqüencial das 
equações de saturação – novembro 2004. 
 
• Implementação do modelo black-oil como extensão do código do bifásico 
(descrito no capítulo 3-3.1.1) com resolução simultânea das equações de 
saturação – abril 2005. 
 
• Implementação do modelo black-oil como extensão do código do bifásico 
(descrito no capítulo 3-3.1.1) com resolução seqüencial das equações de 
saturação – abril 2005. 
 
• Exemplos numéricos para validação do código – setembro 2005. 
 
• Redação da tese – março 2006. 
 
• Data prevista para defesa – março 2006 
 19
 21
REFERÊNCIAS BIBLIOGRÁFICAS 
Bear J., Dynamics of Fluids in Porous Media, New York, Dover Publications Inc. (1972). 
Brezzi, F. and Fortin, M., Mixed and Hybrid Finite Element Methods, Vol. 15 of Springer Series 
in Computational Mathematics, Springer-Verlag, New York (1994). 
Brooks A.N. and Hughes T.J.R., Streamline Upwind/Petrov-Galerkin formulation for convection 
dominated flows with particular emphasis on the incompressible Navier-Stokes equations, 
Computer Methods in Applied Mechanics and Engineering, 32, 199-259 (1982). 
Chavent G., Cohen G., Jaffre J., Eymard R., Guerillot D. R. and Weill L., Discontinuous and 
Mixed Finite Elements for Two-Phase Incompressible Flow, SPERE, 567-575 (1990). 
Coutinho A.L.G.A. and Alves J.L.D., Finite element simulationof nonlinear viscous fingering in 
miscible displacements with anisotropic dispersion and nonmotonic viscosity profiles, 
Computational Mechanics, 23, 108-116 (1999). 
Coutinho A.L.G.A., Silva A.S. and Devloo P.R.B., Uma comparação de esquemas de 
estabilização para a simulação por elementos finitos de escoamentos imiscíveis bi-fásicos em 
meios porosos, Rev. Int. Mét. Num. Cálc. Dis. Ing., 19, 279-294 (2003). 
De Abreu E.C., Simulação Numérica de Escoamentos Trifásicos Água-óleo-gás em 
Reservatórios de Petróleo, Dissertação de Mestrado, IPRJ (2003). 
Douglas Jr. J., Peaceman D. W. and Rachford Jr. H. H., A Method for Calculating Multi-
Dimensional Immiscible Displacement, Trans.SPE AIME, 216, 297-306 (1959). 
Durlofsky L.J., A Triangle Based Mixed Finite Element-Finite Volume Technique for Modeling 
Two Phase Flow Through Porous Media, J. Comp. Physics, 105, 252-266 (1993). 
Durlofsky L.J., Accuracy of Mixed and Control Volume Finite Element Approximations to 
Darcy Velocity and Related Quantities, Water Resour. Res. 30(4): 965-973 (1994). 
Ewing R.E., Multidisciplinary Interactions in Energy and Environmental Modeling, J. Comput. 
Appl. Math. 74: 193-215 (1996). 
Franca L.P., Frey S.L. and Hughes T.J.R., Stabilized Finite Element Methods: I. Application to 
the Advective-Difusive Model, Comput. Methods Appl. Mech. Engrg. 95: 253-276 (1992). 
Fung L.S., Hiebert A.D. and Nghiem L.X., Reservoir Simulation with a Control-Volume Finite-
Element Method, SPERE, 349-357 (1992). 
Guzmán R.E., Mathematics of Three-Phase Flow, PhD Dissertation, Stanford University, 
Stanford (1995). 
Hughes T.J.R., Multiscale phenomena: Green’s functions, the Dirichlet-to-Neumann 
formulation, subgrid scale models, bubbles and the origins of stabilized methods, Comput. 
Methods Appl. Mech. Engrg. 127: 387-401 (1995). 
Hughes T.J.R. and Brooks A.N., A multidimensional upwind scheme with no crosswind 
difusion, in T.J.R. Hughes (ed.), Finite Element Methods for Convection Dominated Flows, 
ASME, New York, 19-35 (1979). 
 22
Hughes T.J.R., Franca L.P. and Hulbert G.M., A New Finite Element Formulation for 
Computational Fluid Dynamics: VIII. The Galerkin Least-Squares Method for Advective-
Difusive Equations, Comput. Methods Appl. Mech. Engrg. 73: 173-189 (1989). 
Juanes R. and Patzek T.W., Multiple-Scale stabilized finite elements for the simulation of tracer 
injections and waterflood, SPE Journal 75231 (2001). 
Juanes R., Displacement Theory and Multiscale Numerical Modeling of Three-phase Flow in 
Porous Media, Thesis, University of California Berkeley (2003). 
Li B., Chen Z. and Huan G., The Sequential Method for the Black-Oil Reservoir Simulation on 
Unstructured Grids, Journal of Computational Physics, 192 (2003) 36-72. 
MacDonald R. C. and Coats K.H., Methods for numerical simulation of water and gas coning, 
Trans. SPE AIME 249 (1970) 425-436. 
Peaceman D. W., Fundamentals of Numerical Reservoir Simulation, Amsterdam, Elsevier 
Scientific Publishing Company, 1977. 
Transgenstein J.A. and Bell J. B., Mathematical Structure of the Black-Oil Model for Petroleum 
Reservoir Simulation, SIAM J. Appl. Math., 49 (1989) 749-783. 
Wu Y. and Forsyth P. A., On the Selection of Primary Variables in Numerical Formulation for 
Modeling Multiphase Flow in Porous Media, Journal of Contaminant Hydrology 48 (2001) 277-
304.