Método das Diferenças Finitas Aplicado na Solução Numérica da Equação de Difusão Hidraúlica Linear do Óleo unidimensional

Método das Diferenças Finitas Aplicado na Solução Numérica da Equação de Difusão Hidraúlica Linear do Óleo unidimensional Júlia Borges dos Santos Universidade Federal de Uberlândia - Faculdade de Engenharia Civil Graduando em Engenharia Civil- PROMAT/FAMAT/UFU

ju_ jbs@ hotmail. com

Santos Alberto Enriquez-Remigio Universidade Federal de Uberlândia - Faculdade de Matemática Professor Adjunto I

santos@ famat. ufu. br No presente trabalho resolve-se numericamente a equação de difusão hidráulica linear do óleo pelo seguinte esquema numérico baseado no método das diferenças nitas: Euler no tempo e aproximação para a derivada segunda que considera estêncil de três pontos. Devido à extensão dos reservatórios de óleo é costume o uso de malha mista: na, perto do poço e grossa longe deste. Neste caso, considerou-se uma malha com espaçamento homogêneo perto do poço e heterogêneo, baseado numa progressão geométrica crescente, após o m da homogênea. O código implementado é validado resolvendo-se a equação de difusão linear com solução manufaturada e logo aplicado na solução da equação de difusão hidráulica linear do óleo. Resultados numéricos mostram: 1) bom comportamento para a aproximação da equação com solução manufaturada, em que convergência de segunda ordem no espaço é obtida; 2) erros pequenos na aplicação do mesmo esquema na resolução numérica da equação do óleo, embora parte da malha usada seja heterogênea com espaçamentos crescentes. Resumo:

1 Introdução A exploração de petróleo deve ser realizada com planejamento e conhecimento da área, para obterse a viabilidade econômica. Para isto é necessário o estudo das propriedades das rochas, propriedades dos uidos, e das leis físicas que regem o movimento dos uidos no interior das rochas. Com isso, a equação da difusividade hidráulica se torna a síntese do problema em questão, pois sua solução nos permite observar como a pressão no reservatório se comporta nas três direções e à medida que o tempo passa, possibilitando um estudo a respeito da produtividade do reservatório em questão [7]. Porém, esta é uma equação diferencial de segunda ordem, e sua solução analítica (a depender do tipo de regime do uido que está sendo considerado), é complexa como a maioria das equações que regem fenômenos de importância prática. Neste contexto, os métodos numéricos se tornam essenciais, especialmente na área de reservatórios de petróleo na qual os problemas são muito complexos para serem analisados por técnicas clássicas [7, 6]. Devido a importância da simulação numérica no estudo de diferentes fenômenos, vários trabalhos vem sendo desenvolvidos em todo o Brasil [6, 8, 7]. Diante da relevância da solução da equação da difusividade hidráulica para a área de reservatórios de petróleo, este trabalho tem como objetivo solucionar numericamente esta equação para regimes especícos, os quais possuem solução analítica conhecida, a m de comparar os resultados e mostrar a precisão das soluções numéricas para que, em trabalhos futuros, possamos solucionar com segurança problemas que não possuem solução analítica, contribuindo para um maior conhecimento e entendimento destes fenômenos físicos.

2 Equação de Difusão Hidráulica do Óleo Linear 1D Considerando-se um reservatório de óleo como mostrado na Figura 2.1 e admitindo-se as seguintes hipóteses: meio poroso homogêneo e isotrópico, permeabilidade constante, compressibilidade e vis-

1

Figura 2.1: Fluxo linear 1D do óleo.

cosidade do uido constantes, e somente um uido no meio; a equação da difusão não linear do óleo é dada por [7]:

∂2p cρ 2 + c2 ρ ∂x



∂p ∂x

2

∂2p + cρ + c2 ρ ∂y



cp ∂y

2

onde c é a compressibilidade do uido, p é a pressão, do meio,

µ

∂2p + cρ 2 + c2 ρ ∂z



∂p ∂z

2 =

φµct ∂p cρ , k ∂t

(2.1)

ρ é a massa especíca do uido, φ é a porosidade

é a viscosidade do uido, k é a permeabilidade do meio, e t é o tempo.

Admitindo-se que a compressibilidade e os gradientes de pressão para o óleo são valores muito pequenos, os termos

c2 ρ



∂p ∂x

2

c2 ρ



∂p ∂x

2

c2 ρ



∂p ∂x

2

da equação (2.1) podem ser considerados desprezíveis

[7] e dessa forma a Eq. (2.1) ca linear e é dada por:

φµct ∂p ∂2p ∂2p ∂2p + 2+ 2 = 2 ∂x ∂y ∂z k ∂t

(2.2)

Considerando-se unicamente variação da pressão na direção x e denindo-se

η =

k φµct , então a

equação (2.2) pode ser escrita:

∂p ∂2p =η 2 ∂t ∂x

2.1

(2.3)

Condições auxiliares para a equação de difusão Para a equação de difusão (4.1) ser resolvida, precisa-se das seguintes condições auxiliares: condição

inicial e condições de contorno. A condição inicial é dada por:

CI :

p(x, t = 0) = p0 (x),

e as condições de contorno consideradas no presente trabalho, são:

(I)

( ∂p (x = a, t) = p1 (t) CCI : ∂x CCE : p(x = b, t) = p2 (t)

ou (II)

( CCI : CCE :

∂p ∂x (x ∂p ∂x (x

= a, t) = p1 (t) = b, t) = p3 (t)

(2.4)

As condições de contorno dadas em (I) são do tipo Neumann em x= a e Dirichlet em x=b, enquanto que as dadas em (II) são do tipo Neumann em x=a e x=b. No caso de reservatório de óleo 1D, podem ser consideradas as seguintes expressões para

p1 (t), p2 (t)

e

p3 (t)

[7]:

p0 (x) = pi onde

q , µ, ct , A

e

p0 (x),

pi

p1 (t) =

qµ , kA

p2 (t) = pi ,

p3 (t) = 0.

são, respectivamente, a vazão, a viscosidade, compressibilidade total, área e

pressão inicial do reservatório. As condições de contorno dadas em (I) e (II) denem, na área do reservatório de óleo, regimes de uxo conhecidos, respectivamente como permanente a pseudo-permanente.

2

2.2

Solução analítica A solução analítica para a equação de difusão linear (2.3), sujeita à condição inicial e as condições

de contorno dadas em (I), é (Ver [5]):

 p(x, t) = pi −

qw µxf  x 8 1− − 2 kA xf π

∞ X

−(nπ)2

n=impar

1 e n

ηt 4x2 f

 cos(

nπ  x) , 2xf

(2.5)

e sujeita as condições de contorno (II) é:

p(x, t) = pi + C(x −

∞ Cxf 2Cxf X exp(−λ2n ηt) Cη x2 )− t− + cos(λn x), 2xf xf 3 π2 n2

(2.6)

n=1

sendo:

λn =

nπ xf e

C=

qw µxf kA .

Na Figura 2.2 e 2.3, apresenta-se o esboço dos grácos associados ao regime permanente e pseudopermanente, respectivamente, e associados aos parâmetros indicados na Tabela 1:

2)

k(m

p

−16 20.9.86923.10

N m2




4 4905.10

φ

N m2 .s −3 0.18.10

0.18






µ

q

m3 s



−3 0.1.0.1667.10

cf 0

cr



m2 N



1 −8 9.81 .10

2

A(m )

xf(m)

60.60

1000

Tabela 1: Parâmetros considerados para o reservatório de óleo.

Figura 2.3: tempo nal = 1728000s (20dias).

Figura 2.2: tempo nal = 8640000s (100dias).

3 Método das Diferenças Finitas É um dos métodos usados para aproximar derivadas. Tais aproximações são obtidas usando-se a série de Taylor [3]. As aproximações usadas para a derivada temporal e espacial da Eq. (2.3) são, respectivamente:

pn+1 − pni ∂p ≈ i ∂t ∆t onde o índice i indica o ponto

xi ,

∂2p 4xi pi−1 + 4xi pi+1 − (4xi−1 + 4xi )pi ≈ , 2 ∂x 4xi−1 4xi ( 4xi−1 +4xi ) 2

e n, indica o tempo

Quando a malha é homogênea, isto é

4xi = 4xi−1

para derivada segunda (ordem de aproximação dois):

3

tn (Veja [3]). = 4xi+1 , chegamos

à seguinte aproximação

∂2p pi−1 − 2pi + pi+1 ≈ 2 ∂x 4x2 A expressão discreta usada para resolver a Eq. (2.3) é:

pn+1 − pni i ≈α 4t

4xi pi−1 + 4xi pi+1 − (4xi−1 + 4xi )pi

!

4xi−1 4xi ( 4xi−12+4xi )

e para aproximar a condição de contorno do tipo Neumann (Veja (2.4)) foi usada aproximação de segunda ordem não centradas (Veja [1, 2, 3]).

4 Resultados Numéricos Nas seguintes seções apresentam-se: 1) resultados numéricos da vericação numérica da ordem de convergência do método, garantindo dessa forma a correta implementação do código; 2) resultados da solução da equação de difusão linear do óleo para o regime permanente e pseudo-permanente no tempo t=100s. A malha considerada é mista: homogêna em todo o comprimento  d e heterogênea no restante do domínio. Na Figura 4.1 pode-se ver um exemplo do tipo de malha usada no presente trabalho.

Figura 4.1: Exemplo de malha mista usada no presente trabalho: homogênea perto do poço (x0 ) e heterogênea, com espaçamento dado por um progressão geométrica crescente, longe do poço.

φe e aproximado φa é calculado usando-se a seguinte v uN −1 uX =t (φe (xi , t) − φa (xi , t))2 ∆x2i ,

A norma dos erro entre o valor da função exata fórmula:

erro = kφe − φa kL2

i=0 onde

N

é o número de subdivisões do intervalo de cálculo

[0, xf ].

O método numérico usado é explícito, portanto o passo no tempo, condição:

4.1

∆t <

∆x2 2η .

∆t,

deve satisfazer a seguinte

Vericação numérica da ordem de convergência do método implementado Impondo-se a seguinte função

p(x, t) = 0.001 cos(t)arctg(x)

como solução da seguinte equação de

difusão linear:

obtem-se a seguinte expressão

∂p ∂2p = 0.084 2 + f (x, t) ∂t ∂x   xcos(t) f (x, t) = 0.001 −sen(t)arctan(x) + 0.168 (1+x . 2 )2

(4.1) Com essas infor-

mações resolveu-se numericamente a equação de difusão (4.1) para as condições de contorno dadas por (I) e (II). Nas seguintes seções N1 e N2 indicam o número de subintervalos com espaçamento homogêneo e heterogêneo, respectivamente, na qual foi dividida o domínio unidimensional,

4

xf .

4.1.1

Condições de contorno do tipo (I)

Observando os dados da Tabela 2 e 3, percebe-se que confere a ordem dois do método numérico implementado para as malhas indicadas, tanto para um tempo muito pequeno, como para outro um pouco maior.

d N1

h=

1 N1

=

N2(r=1)

N2(r=1.2)

4

156

19

0.125

0.0002380

8

312

22

0.0625

0.0000605

3.934245

16

624

26

0.03125

0.0000151

32

1248

30

0.015625

0.0000037

64

2496

34

0.0078125

0.00000095

d N1

1 N1

=

q=

erro(r=1)

erro(r=1.2)

q=

erro(h) erro(h/2)

0.0002381 0.0000606

3.927642

3.983368

0.0000152

3.969591

4.001695

0.00000386

3.952403

3.996047

0.00000102

3.767017

x ∈ [0, 20], d = 0.5

Tabela 2: CC tipo (I):

h=

erro(r=1)

erro(h) erro(h/2)

N1

e tempo nal =1s.

q=

erro(h) erro(h/2)

N2(r=1)

N2(r=1.2)

4

156

19

0.125

0.0002671

8

312

22

0.0625

0.0000663

4.022900

0.0000662

4.017795

16

624

26

0.03125

0.0000165

4.003421

0.0000165

3.989590

32

1248

30

0.015625

0.0000041

4.000868

0.0000042

3.949687

64

2496

34

0.0078125

0.0000010

4.000295

0.0000011

3.757087

q=

0.0002660

x ∈ [0, 20], d = 0.5

Tabela 3: CC tipo (I):

erro(r=1.2)

erro(h) erro(h/2)

N1

e tempo nal = 20s.

4.1.2 Condições de contorno tipo II A ordem de convergência do método também é conferida para este regime. Porém, estes tempos são muitos pequenos em comparação à prática. Resta saber se, para tempos maiores, está ordem será mantida, podendo assim ser aplicado o método com segurança para problemas práticos.

h=

d N1

=

1 N1

q=

N2(r=1.2)

4

156

19

0.125

0.0002400

8

312

22

0.0625

0.0000610

3.929841

0.0000606

3.927258

16

624

26

0.03125

0.0000153

3.981324

0.0000152

3.969039

32

1248

30

0.015625

0.0000038

4.000716

0.0000038

3.950598

64

2496

34

0.0078125

0.0000009

3.995570

0.0000010

3.761288

x ∈ [0, 20], d = 0.5

5

erro(r=1.2)

q=

erro(h) erro(h/2)

N2(r=1)

Tabela 4: CC (II):

erro(r=1)

erro(h) erro(h/2)

N1

0.0002383

e tempo nal=1s.

h=

d N1

=

1 N1

N2(r=1.2)

4

156

19

0.125

0.0002683

8

312

22

0.0625

0.0000667

4.021641

0.0000662

4.018851

16

624

26

0.03125

0.0000166

4.002189

0.0000166

3.989414

32

1248

30

0.015625

0.0000041

4.000097

0.0000042

3.947987

64

2496

34

0.0078125

0.0000001

3.999871

0.0000012

3.751465

4.2

q=

erro(r=1.2)

q=

erro(h) erro(h/2)

N2(r=1)

Tabela 5: CC (II):

erro(r=1)

erro(h) erro(h/2)

N1

0.0002663

x ∈ [0, 20], d = 0.5

e tempo nal=20s.

Solução Numérica da Equação Linear do Óleo Foi calculado a solução numérica da Eq. (2.3) nos tempos

tf = 5s, 100s, 200s

para os regimes

permanente e pseudo-permanente, sendo que os parâmetros para o reservatório de óleo são os indicados na Tabela 1. Nos primeiro d=100m do comprimento do reservatório usou-se uma malha homogênea e no resto uma malha heterogênea com espaçamento crescente denido por uma progressão geométrica de razão r=1.2. Na Tabela 6, pode ser observado que os erros obtidos nos tempo

tf = 5s,100s,

e

200s

para o regime permanente são aparentemente grandes, enquanto que para o regime pseudo-permanente são pequenos. Tal armação não é verdadeira, pois os erros são insignicantes quando comparados com a ordem de grandeza da pressão ( aproximadamente

N1

N2

RP

2000

44

RPP

2000

44

N 106 m 2

).

∆t

kpe − pa kL² (tf = 5s)

kpe − pa kL² (tf = 100s)

kpe − pa kL² (tf = 200s)

5e-02

2.091569e-03

3.077825

1.449694

1.218571

5e-02

2.091569e-03

1.195117e-07

1.199403e-07

1.681794e-07

h=

Tabela 6: Erro nos tempos

d 2000

t = 5s, 100s, 200s

para o regime permanente (RP) e regime pseudo-

permanente (RPP) usando-se o método de Euler no tempo com estêncil que usa três pontos para a derivada segunda.

5 Conclusões No presente trabalho, resolveu-se a equação de difusão linear (2.3) usando-se o método de Euler para a derivada temporal e diferenças nitas que usa estêncil de três pontos para a segunda derivada espacial. A malha usada foi homogênea no início do domínio, e heterogênea no restante, com espaçamento dado por uma progressão geométrica. Embora a malha toda não seja igualmente espaçada, e nem o método usado seja de alta ordem, os resultados se mostraram satisfatórios. Observando-se os resultados para a solução manufaturada, verica-se o comportamento de um método de segunda ordem, mostrando-se um bom meio para se solucionar o problema, para as condições especicadas.

Em relação aos resultados da solução numérica, percebe-se grandes erros

para o regime permanente, e pequenos para regime pseudo-permanente. Apesar desta divergência, os dois regimes apresentaram bons resultados. Isto devido ao fato de que os erros de aproximação, para a equação de difusão linear do óleo, se tornam desprezíveis quando comparados à ordem de grandeza dos valores da pressão do óleo para a solução aproximada. A despeito disso, é importante ressaltar que apesar deste bom comportamento, isto não implica que se obterão bons resultados para a solução da equação de difusão linear de uma forma genérica. As simulações foram realizadas para regimes especícos, com parâmetros físicos também particulares, sendo os resultados intimamente ligados a estes fatores.

6

Referências [1] Aziz, K. Settari, A.,  Petroleum Reservoir Simulation , Applied Science Publisher, 1979. [2] Ferziger, J.H. and Peric, M., Computational Methods for Fluid Dynamics, Springer, 2002. [3] Fortuna, A. O., Técnicas Computacionais para Dinâmicas dos Fluidos. Conceitos Básicos e Aplicações, EDUSP, 2000. [4] Neide, B.F., Cálculo Numérico, EDUSP, (2007), [5] LENEP-UENF, Fluxo Monofásico de Fluidos em Meios Porosos, Equipe de Análise de Testes de Pressão, 2008. [6] Oliveira, B.C., Simulação de reservatórios de petróleo em ambiente MPSoC, Dissertação de mestrado, Universidade Federal do Rio Grande do Norte, 2009. [7] Rosa, A.J., Carvalho, R.S., Xavier, J.A.D.,  Engenharia de Reservatórios de Petróleo , Interciência, Brasil, 2006. [8] Silva, S.A., Solução numérica da equação de difusão de calor em geometrias tridimensionais, Dissertação de mestrado, Universidade Federal de Pernanbuco, 2004.

7