APLICAÇÃO DO MÉTODO DOS ELEMENTOS DE CONTORNO …objdig.ufrj.br/60/teses/coppe_m/ConradoKeidel.pdf · iii Keidel, Conrado Aplicação do Método dos Elementos de Contorno na Modelagem

APLICAÇÃO DO MÉTODO DOS ELEMENTOS DE CONTORNO NA MODELAGEM DE

TESTES DE PRESSÃO EM POÇOS DE PETRÓLEO

Conrado Keidel

Dissertação de Mestrado apresentada ao Programa

de Pós-graduação em Engenharia Civil, COPPE,

da Universidade Federal do Rio de Janeiro, como

parte dos requisitos necessários à obtenção do

título de Mestre em Engenharia Civil.

Orientadores: José Antonio Fontes Santiago

José Claudio de Faria Telles

Rio de Janeiro

Março de 2011



Conrado Keidel

DISSERTAÇÃO SUBMETIDA AO CORPO DOCENTE DO INSTITUTO ALBERTO LUIZ

COIMBRA DE PÓS-GRADUAÇÃO E PESQUISA DE ENGENHARIA (COPPE) DA

UNIVERSIDADE FEDERAL DO RIO DE JANEIRO COMO PARTE DOS REQUISITOS

NECESSÁRIOS PARA A OBTENÇÃO DO GRAU DE MESTRE EM CIÊNCIAS EM

ENGENHARIA CIVIL.

Examinada por:

________________________________________________

Prof. José Antonio Fontes Santiago, D.Sc.

________________________________________________

Prof. José Claudio de Faria Telles, Ph.D.

________________________________________________

Prof. Alvaro Marcello Marco Peres, Ph.D.

________________________________________________

Prof. Ricardo Alexandre Passos Chaves, D.Sc.

RIO DE JANEIRO, RJ - BRASIL

MARÇO DE 2011

iii

Keidel, Conrado

Aplicação do Método dos Elementos de Contorno na

Modelagem de Testes de Pressão em Poços de Petróleo/

Conrado Keidel. – Rio de Janeiro: UFRJ/COPPE, 2011.

XV, 137 p.: il.; 29,7 cm.



Dissertação (mestrado) – UFRJ/ COPPE/ Programa de

Engenharia Civil, 2011.

Referencias Bibliográficas: p. 125-127.

1. Método dos Elementos de Contorno. 2. Testes de

Pressão em Poços de Petróleo. 3. Espaço de Laplace. 4.

Algoritmo de Stehfest. I. Santiago, José Antonio Fontes et al.

II. Universidade Federal do Rio de Janeiro, COPPE, Programa

de Engenharia Civil. III. Título.

iv

À memória de meu querido pai, Carlos,

Ao vigor de meu querido filho, Gustavo.

v

Agradecimentos

À minha mãe, Selene. Às minhas irmãs, Gisela e Lívia.

À minha esposa, Renata e sua família de origem.

Ao Antonio Carlos Decnop Coelho e, por extensão, à PETROBRAS, por terem me

concedido um período de reflexão jamais experimentado anteriormente.

À COPPE e todo seu corpo docente, pelos ensinamentos transmitidos. Em especial, aos

orientadores de dissertação, por terem me fornecido de bandeja um conhecimento por mim tão

almejado.

Ao orientador acadêmico e professor Álvaro Coutinho.

Aos colegas de jornada, especialmente aos excepcionais Edivaldo Fontes Jr. e Marlúcio

Barbosa.

Aos colegas de Laboratório de Mecânica Computacional (LAMEC). Principalmente, ao

Cleberson Dors.

vi

Resumo da Dissertação apresentada à COPPE/UFRJ como parte dos requisitos necessários

para a obtenção do grau de Mestre em Ciências (M.Sc.).



Conrado Keidel

Março/2011



Programa: Engenharia Civil

O fluxo de fluidos pouco compressíveis em reservatórios de petróleo pode ser modelado

analiticamente através de equações integrais compostas por funções de Green. O Método dos

Elementos de Contorno (MEC) corresponde à solução numérica de tais equações. Para

elaboração desta dissertação são implementados programas computacionais em VBA,

baseados no MEC, capazes de lidar com fluxo monofásico para poços verticais em domínios

simples e compostos de múltiplas regiões homogêneas. As soluções são obtidas via campo de

Laplace e são invertidas para o campo real por meio do algoritmo de Stehfest. Os resultados

gerados a partir de modelos simples são analisados em confronto com soluções analíticas

constantes da literatura especializada sobre análise de testes em poços de petróleo.

vii

Abstract of Dissertation presented to COPPE/UFRJ as a partial fulfillment of the requirements

for the degree of Master of Science (M.Sc.)

APPLYING THE BOUNDARY ELEMENTS METHOD FOR MODELING

PETROLEUM WELL TESTS

Conrado Keidel

March/ 2011

Advisors: José Antonio Fontes Santiago


Department: Civil Engineering

Flow of slightly compressible fluids through petroleum reservoirs can be analytically

modeled by means of integral equations composed by Green’s functions. The Boundary

Element Method (BEM) corresponds to the numerical solutions of such equations. For the

preparation of this dissertation they are implemented VBA computing codes, based on BEM,

capable of dealing with single phase flow to vertical wells located in simple and multi-region

composed domains. The solutions are obtained via Laplace space and are inverted back to the

real field with the aid of Stehfest’s algorithm. Simple models are run and their results are then

compared with analytical solutions taken from specialized literature about well test analysis.

viii

Sumário 1. Introdução........................................................................................................................... 1

1.1. Considerações Preliminares........................................................................................ 1

1.2. Revisão Bibliográfica ................................................................................................. 5

1.3. Objetivos e Conteúdos................................................................................................ 6

1.4. Acerca da Nomenclatura Adotada.............................................................................. 8

2. Relações Básicas............................................................................................................... 10

2.1. Conceitos Matemáticos Essenciais........................................................................... 10

2.2. Operações no Espaço de Laplace ............................................................................. 11

2.3. Equações de Dinâmica dos Fluidos .......................................................................... 13

2.3.1. Lei da Conservação de Massas (ou Equação da Continuidade) ....................... 14

2.3.2. Lei de Darcy ..................................................................................................... 16

2.3.3. Equações de Estado .......................................................................................... 19

2.4. Equação da Difusão Hidráulica ................................................................................ 21

2.5. Condição Inicial, Condições de Contorno e Condições de Interface Empregadas .. 24

2.5.1. Condição Inicial................................................................................................ 24

2.5.2. Condições de Contorno .................................................................................... 25

2.5.3. Condições de Interface ..................................................................................... 25

2.6. Superposição no Espaço ........................................................................................... 26

2.7. Superposição no Tempo (Convolução) .................................................................... 28

2.8. Soluções Analíticas para Escoamento de Fluidos em Meios Porosos Baseadas em

Funções de Green ................................................................................................................. 30

2.8.1. Conceitos Preliminares – Delta de Dirac.......................................................... 30

2.8.2. Funções de Green ............................................................................................. 32

2.8.3. Equação Integral Ponderada da Difusão – Ponto-fonte Pertencente ao Domínio

33

2.8.4. Equação Integral Ponderada de Difusão – Ponto-fonte Pertencente ao Contorno

39

2.9. Soluções Fundamentais ............................................................................................ 42

2.9.1. Considerações Preliminares.............................................................................. 42

ix

2.9.2. Funções de Green Básicas (ou Soluções Fundamentais) ................................. 43

3. MEC Aplicado ao Escoamento em Meios Porosos .......................................................... 51

3.1. Aspectos Gerais ........................................................................................................ 51

3.2. Variação Temporal Aproximada por Série de Taylor (Diferenças Finitas) ............. 51

3.3. Solução no Campo Real com a Solução Fundamental Dependente do Tempo....... 53

3.4. Solução Via Campo de Laplace................................................................................ 53

3.5. Procedimento Numérico ........................................................................................... 55

3.5.1. Domínios Simples ............................................................................................ 55

3.5.2. Domínios Compostos ....................................................................................... 58

3.6. Considerações a Respeito das Implementações ....................................................... 64

3.6.1. Generalidades ................................................................................................... 64

3.6.2. Domínios Simples ............................................................................................ 66

3.6.3. Domínios Compostos ....................................................................................... 74

4. Aplicações Numéricas ...................................................................................................... 79

4.1. Considerações Iniciais .............................................................................................. 79

4.1.1. Metodologia...................................................................................................... 79

4.1.2. Base de Comparação ........................................................................................ 81

4.2. Discussão Acerca do Consumo de Tempo ............................................................... 82

4.3. Avaliação da Exatidão dos Códigos Gerados........................................................... 91

4.4. Limitações do Algoritmo de Stehfest ..................................................................... 118

5. Conclusões e Sugestões.................................................................................................. 121

5.1. Conclusões.............................................................................................................. 121

5.2. Sugestões ................................................................................................................ 122

Referências Bibliográficas...................................................................................................... 125

Apêndice A ............................................................................................................................. 128

Apêndice B ............................................................................................................................. 129

Apêndice C ............................................................................................................................. 131

Apêndice D............................................................................................................................. 132

Apêndice E ............................................................................................................................. 135

Apêndice F.............................................................................................................................. 136

x

Lista de Figuras

Figura 1.1 - Contexto em que se insere o estudo em objeto....................................................... 2

Figura 1.2 – O problema inverso de interpretação de testes em poços de petróleo.................... 4

Figura 1.3 - "Fio condutor” desta revisão bibliográfica. ............................................................ 6

Figura 1.4 - Domínio típico de interesse. ................................................................................... 8

Figura 2.1 – Ilustração para o teorema da divergência............................................................. 11

Figura 2.2 - Aplicação do Método das Imagens....................................................................... 26

Figura 2.3 - Barreiras ao fluxo: em ângulo de 90º e em ângulo qualquer. ...............................27

Figura 2.5 – Convolução com o delta de Dirac. ....................................................................... 31

Figura 2.6 – Domínio Considerado. ......................................................................................... 33

Figura 2.7 - Fonte nas proximidades de um vértice de contorno. ............................................ 39

Figura 2.8 - Fonte pontual nas proximidades de contorno suave. Sistema equivalente. .......... 41

Figura 2.9 – Soluções para domínios ortotrópicos obtidas por produtos de Newman. ............ 50

Figura 3.1– Representação em elementos de contorno............................................................ 55

Figura 3.3 – Exemplo esquemático de um domínio composto................................................. 58

Figura 3.4 – Porções do domínio composto destacadas. .......................................................... 59

Figura 3.5 – Aspecto das matrizes............................................................................................ 60

Figura 3.6 – Exemplo esquemático de domínio composto com região interna........................ 63

Figura 3.7 – Porções destacadas do domínio composto com região interna. ........................... 63

Figura 3.8- Fluxograma de cálculo para MEC aplicado ao Espaço de Laplace com inversão

pelo algoritmo de Stehfest. ....................................................................................................... 65

Figura 3.9 – Domínio simples hipotético. ................................................................................ 66

Figura 3.10 – Parâmetros de entrada 1/2. ................................................................................. 67

Figura 3.11 – Parâmetros de entrada 2/2. ................................................................................. 67

Figura 3.12 - Integração Numérica – quando x e ξ não pertencem ao mesmo elemento. ....... 68

Figura 3.13 - Integral da Função de Bessel Modificada de 2ª Espécie e Ordem 0................... 71

Figura 3.14 – Exemplo para Domínio Composto: conectividade. ...........................................74

Figura 3.15 – Exemplo: tabelas de entrada 1/3. ....................................................................... 75

Figura 3.16 – Exemplo: dados de entrada 2/3. ......................................................................... 75

Figura 3.17 – Exemplo: dados de entrada 3/3. ......................................................................... 76

xi

Figura 3.18 – Sentido positivo de percurso para as regiões consideradas de forma isolada. ... 76

Figura 3.19 – Aspecto das matrizes globais e do vetor global para o exemplo considerado. .. 77

Figura 4.1 - Avaliação do consumo de tempo. ......................................................................... 83

Figura 4.2 - Comportamento de K0 e K1. ................................................................................ 84

Figura 4.3 - K0(x) com aproximação da fonte pontual na direção ortogonal ao elemento

(unidades SI)............................................................................................................................. 85

Figura 4.4 - K0(x) para variação da posição da fonte na direção do elemento (unidades: SI). 86

Figura 4.5 – Representação do reservatório circular com 6 elementos. ................................... 88

Figura 4.6 – Reservatório circular com poço centralizado (NE=6, Ninter=8, NSteh = 2). ...... 89

Figura 4.7 – Reservatório circular com poço centralizado (NE=6, Ninter=8, NSteh=6). ........ 89

Figura 4.8 - Reservatório circular com poço centralizado (NE=6, Ninter=8, NSteh=10)........ 90

Figura 4.9 – Exemplo 1– NE= 4, NPGL=8. ............................................................................. 92

Figura 4.10 – Exemplo 1– NE= 4, NPGL=32. ......................................................................... 93

Figura 4.11 - Exemplo 1– NE= 16, NPGL=8........................................................................... 94


Figura 4.13 - Exemplo 1– NE= 64 não uniformes, NPGL=8................................................... 96

Figura 4.14 – Exemplo 1: Quedas de pressão em outros pontos do reservatório..................... 97

Figura 4.15 - Exemplo 2– NE= 4, NPGL=8............................................................................. 98

Figura 4.16 - Exemplo 2– NE= 4, NPGL=8............................................................................. 99


Figura 4.18 - Exemplo 2– NE= 7 não uniformes, NPGL=8................................................... 100

Figura 4.19 - Exemplo 3 – NE=18 não uniformes, NPGL=8.................................................101

Figura 4.20 - Exemplo 3 – NE=80, NPGL=8......................................................................... 102

Figura 4.21 - Exemplo 4 – NE=20 não uniformes, NPGL=8.................................................103

Figura 4.22 – Exemplo 5 - NE=25, NPGL=8......................................................................... 104

Figura 4.23 – Exemplo 5 - NE=49, NPGL=8......................................................................... 105

Figura 4.24 - Exemplo 6 – NE=12, NPGL=8......................................................................... 107

Figura 4.25 - Exemplo 6 – NE=42, com comprimento uniforme na direção x, NPGL=8. .... 107

Figura 4.26 - Exemplo 6 – NE=84,com comprimento uniforme na direção x, NPGL=8. ..... 108

Figura 4.27 – Exemplo 7 – NF=21, NPGL=8. ....................................................................... 109

Figura 4.28 - Exemplo 7 – NF=61, NPGL=8........................................................................ 110

xii

Figura 4.29 – Exemplo 7 – NF=121, NPGL=8. ..................................................................... 110

Figura 4.30 – Exemplo 8 - NE de interface = 16, NPGL=8. .................................................. 112

Figura 4.31 – Exemplo 8 - NE de interface = 16, NPGL=8. .................................................. 113

Figura 4.32 - Exemplo 8 - NE de interface = 16, NPGL=8.................................................... 114

Figura 4.33 - Exemplo 8 - NE de interface = 16, NE de fronteira externa=6, NPGL=8........ 114

Figura 4.34 - Exemplo 9 - NE de interface = 8, NPGL=8...................................................... 115

Figura 4.35 - Exemplo 9 - NE de interface = 2x16 + 8, NPGL=8. ........................................ 116

Figura 4.36 - Exemplo 10 - NE de interface = 25, sendo o menor com 10m......................... 117

Figura 4.37 - Exemplo 11 – modelo composto de 3 regiões. ................................................. 118

Figura 4.38 - Exemplo sintético do uso do algorítmo de Stehfest. Caso de vazões variáveis.119

Figura D.1 – Fonte Planar....................................................................................................... 132

xiii

Lista de Tabelas Tabela 4-1 – Regimes de fluxo característicos, observados através da curva de derivadas de

queda de pressão, utilizados neste trabalho. ............................................................................. 80

Tabela 4-2- Parâmetros adotados para os estudos iniciais. ...................................................... 82

Tabela 4-3 – Tempo consumido em função do número de pontos de Gauss-Legendre

Considerado apenas 1 ponto interno, representando o poço. ................................................... 87

Tabela 4-4 – Exemplos 1, 2, 3 e 4: parâmetros considerados. ................................................. 92

Tabela 4-5 – Exemplo 5: parâmetros considerados................................................................ 104

Tabela 4-6 – Exemplo 6: parâmetros considerados................................................................ 106

Tabela 4-7 - Exemplo 7: parâmetros considerados. ............................................................... 109

Tabela 4-8 - Exemplo 8: parâmetros considerados. ............................................................... 111

Tabela 4-9 –Exemplo 10: parâmetros considerados............................................................... 116

Tabela 4-10 –Exemplo 11: parâmetros considerados............................................................. 117

Tabela A-1: Coeficientes de Stehfest. .................................................................................... 128

xiv

Nomenclatura Letras gregas

ξ - posição do ponto-fonte. η - constante de difusão hidráulica ( µφη tck= ).

φ - porosidade de reservatório.

µ - viscosidade dinâmica de fluido.

Ω - domínio.

Γ - contorno.

Letras latinas

c - compressibilidade.

h - altura da formação.

k - permeabilidade da formação.

L - operador diferencial linear; transformada de Laplace; comprimento.

p – pressão.

q– vazão.

q~ - vazão específica, por unidade de superfície ou de comprimento de fonte.

r - distância.

t - tempo.

u - queda de pressão, p∆ , exceto onde indicado.

u - queda de pressão especificada (conhecida).

u - queda de pressão no espaço de Laplace.

x - posição do ponto de integração de elemento (ponto-campo).

Subscritos

0 – inicial.

c – constante.

e – fronteira externa. Contorno.

f – fratura.

i – contador associado à posição da fonte.

xv

j – contador associado ao elemento considerado.

k – contador associado aos pontos de integração de Gauss-Legendre.

r – rocha.

t – total.

w – fonte ou poço (well).

x – direção do eixo coordenado x.

y - direção do eixo coordenado y.

z - direção do eixo coordenado z.

Siglas/ Abreviaturas

ECDS – Programa computacional desenvolvido com base no MEC para Domínio Simples.

ECDC - Programa computacional desenvolvido com base no MEC para Domínio Composto.

MDF – Método das Diferenças Finitas.

MEC – Método dos Elementos de Contorno.

MEF – Método dos Elementos Finitos.

NE – Número de Elementos.

NPGL – Número de pontos para integração numérica por Gauss-Legendre.

NPT – Número de passos de tempo.

NSteh – Número de passos de tempo de acordo com o algoritmo de Stehfest.

CAPÍTULO 1: INTRODUÇÃO

1

1. Introdução

1.1. Considerações Preliminares

A fonte de energia dominante em todo o globo desde o princípio do século XX tem sido

oriunda de materiais fósseis. Preponderantemente, de óleo e de gás.

Estes compostos são gerados no subsolo do planeta a partir do soterramento progressivo

de rochas ricas em matéria orgânica1. Tal soterramento acarreta um aumento gradual de

temperatura e de pressão no sistema portador e redunda em algumas transformações químicas

na matéria original, que na ausência de grandes quantidades de substâncias oxidantes, acaba

por se transformar em petróleo.

Devido a fatores diversos e ainda controvertidos (ver THOMAS, 2001) o petróleo gerado

posteriormente é expulso da rocha inicial e ascende, migrando até que alguma barreira

geológica o impossibilite de prosseguir seu curso, permanecendo, assim, alojado em uma

rocha sedimentar, porosa e permeável2, denominada reservatório.

Milhões de anos após este processo, o homem intervém: prospectando, perfurando poços

e desenvolvendo planos de extração destes recursos energéticos.

Da necessidade de ordenar e otimizar a explotação do petróleo vem o imperativo de se

entender e simular o movimento dos fluidos em meios porosos.

Este deslocamento, quando a taxas relativamente baixas, se dá preponderantemente por

difusão e é governado pela equação da difusão hidráulica, cuja forma compacta é:

t

uu

∂∂=∇2η . (1.1)

O operador Laplaciano 2∇ corresponde ao divergente do gradiente e é função do sistema

coordenado adotado. A grandeza escalar u representa a queda de pressão3 e é dependente de

posição e de tempo, bem como de características de rocha e fluido, agrupadas pela constante

de difusão η . 1 Destacadamente, folhelhos. 2 Mormente, arenitos e carbonatos. 3 Que está intimamente relacionada ao potencial de fluxo.


2

Ao longo do tempo, vários métodos analíticos têm sido utilizados em busca de soluções

específicas para condições iniciais e de contorno simplificadas. Dentre os quais têm se

destacado os que se utilizam de transformadas de Fourier e, principalmente, os que lançam

mão de transformadas de Laplace4 (ver Figura 1.1, caminho 1).

Figura 1.1 - Contexto em que se insere o estudo em objeto.

Outros métodos analíticos, bastante conhecidos, porém dominados por um número bem

menos expressivo de profissionais da indústria do petróleo, consistem na construção de

soluções a partir do impulso fundamental de Lorde Kelvin.

De posse de soluções fundamentais5, neste contexto, há basicamente duas alternativas para

solução de determinado problema:

4 Tendo como marco o extraordinário artigo de VAN EVERDINGEN e HURST (1949). 5 São muitas vezes referidas como funções de Green. A de Lorde Kelvin, por exemplo, é uma delas.

PROBLEMA FÍSICO-MATEMÁTICO DE

DIFUSÃO HIDRÁULICA +

CONDIÇÃO INICIAL + CONDIÇÕES DE

CONTORNO

SOLUÇÃO FUNDAMENTAL

(IMPULSO)

EQUAÇÃO INTEGRAL DA

DIFUSÃO HIDRÁULICA

EQUAÇÃO DIFERENCIAL DA DIFUSÃO HIDRÁULICA

SOLUÇÃO. MÉTODO DAS IMAGENS PARA CONSIDERAÇÃO DOS LIMITES DE DOMÍNIO.

SOLUÇÃO BASEADA NA INTEGRAÇÃO EM TEMPO E EM ESPAÇO.

MÉTODO DAS IMAGENS PARA CÁLCULO DAS

FRONTEIRAS DE RESERVATÁORIO.

SOLUÇÃO.

1

2

3

F L E X I B I L I D A D E


3

a) ou se integra diretamente em tempo e em espaço, de acordo com suas especificidades.

Neste caso, as fronteiras, caso existam, são computadas com o auxílio do método das

imagens (ver Figura 1.1, caminho 2);

b) ou se busca resolver a equação integral governante, obtida a partir de ponderação da

equação diferencial original pela solução fundamental específica para o problema de

interesse6 (ver Figura 1.1, caminho 3).

A opção disposta em b), acima, é, em teoria, a que confere maior flexibilidade à solução

analítica de problemas, desde que as integrais existentes no processo possam ser calculadas.

Ou seja, aparentemente é a que pode ser aplicada a uma maior variedade de condições de

contorno.

Mas em certos casos, a complexidade a ser avaliada pode ser tal que a obtenção de

soluções analíticas elegantes e explícitas se torne inviável.

Nestas circunstâncias, geralmente se recorre a métodos numéricos de resolução. Dentre os

quais, no âmbito da Engenharia de Petróleo têm se destacado o de diferenças finitas (MDF),

preferido pelas empresas que comercializam os mais difundidos programas específicos para

simulação de escoamento de fluidos em meios porosos, e o de elementos finitos (MEF),

freqüentemente adotado nos meios acadêmicos.

Porém, para ambos é necessária a discretização do domínio do problema em células. E em

decorrência disto podem haver: a) dificuldades no que concerne à orientação de malha – o

fluxo para o poço é radial e a malha geralmente é cartesiana – e b) perda de precisão devido às

aproximações que ocorrem no domínio do problema, entre outros.

O Método dos Elementos de Contorno (MEC), em determinadas conjunturas, constitui-se

excelente opção em relação aos outros dois, já que tende a apresentar como vantagens

(WROBEL et al, 1984; KIKANI et al, 1992; KRYUCHKOV et al, 2001):

- a possibilidade de redução da dimensão do problema7;

- fácil conformação às fronteiras do reservatório;

- eliminação de malha de domínio e, conseqüentemente, ausência de problemas de

orientação de malha e de dispersão numérica;

- excelente acurácia no trato com singularidades;

6 Tal equação é explorada mais adiante e também pode ser diretamente obtida das identidades de Green. 7 Caso a condição inicial seja de potencial constante e uniforme no domínio do problema e caso não haja fontes distribuídas no reservatório.


4

- ausência de aproximação da equação governante no domínio do problema;

- capacidade de lidar com fronteiras móveis e espaços infinitos.

Com efeito, o MEC corresponde à solução numérica da equação da difusão hidráulica

pelo caminho de nº 3, representado na Figura 1.1.

Modelagem de Testes de Pressão em Poços de Petróleo

Testes de pressão em poços de petróleo são operações administradas com o objetivo de se

inferir as propriedades dinâmicas de um sistema poço-reservatório.

Ao sistema, posicionado em subsuperficie e, portanto, inacessível por meios diretos, é

aplicada uma perturbação na forma de retirada de massa, q(t), e é medida sua resposta, em

termos de queda de pressão, u(t).

A partir dos registros q(t) e u(t) se busca uma função de transferência que possa relacioná-

los satisfatoriamente8 (ver Figura 1.2). Habitualmente, este processo, denominado

interpretação, se dá por regressão não-linear entre dados medidos e dados simulados, a partir

da adoção do modelo mais apropriado para o conjunto de informações que se tem a

disposição.

O elenco de modelos analíticos disponíveis é restrito. Um método numérico simples,

robusto e potencialmente veloz como o MEC, neste contexto, pode ampliar e muito as

possibilidades de modelagem.

Figura 1.2 – O problema inverso de interpretação de testes em poços de petróleo.

8 Por definição, um problema inverso.

Sistema (?)

q(t) u(t)

Função de transferência


5

1.2. Revisão Bibliográfica

O fio condutor que dá ensejo à construção da presente dissertação pode ser resumido no

texto abaixo e na Figura 1.3.

CARSLAW e JAEGER (1959) divulgaram e deram ênfase a idéia, na ocasião de cerca de

80 anos de idade, atribuída a Lorde Kelvin, do uso sistemático de soluções fundamentais em

problemas de condução de calor.

Posteriormente, amparados pela extrema similaridade com problemas daquela natureza,

GRINGARTEN e RAMEY (1973), em seu notável artigo, transpuseram aqueles conceitos

para o campo de escoamento em meios porosos9. Ademais, organizaram e sistematizaram o

uso de funções de Green para diversas possibilidades de configuração de poço e de

reservatório de petróleo.

Alguns anos depois, OZKAN e RAGHAVAN (1991a, 1991b) e RAGHAVAN (1993)

ampliaram o escopo do trabalho anterior e apresentaram as correspondentes soluções

analíticas no campo de Laplace.

No que tange ao uso do MEC para problemas de difusão de um modo geral, RIZZO e

SHIPPY (1970) utilizaram-no para aplicações em problemas de condução de calor via

transformadas de Laplace. LIGGET e LIU (1979) resolveram problemas referentes a aqüíferos

subterrâneos via equações integrais no campo real. Um compêndio destas e outras aplicações

pode ser fartamente encontrado em BREBBIA, TELLES e WROBEL (1984).

Especificamente na área de simulação de escoamento em reservatórios de petróleo, o

maior desenvolvimento se deu na Universidade de Stanford. Por MATSUKAWA e HORNE

(1988) e, particularmente, por KIKANI e HORNE (1992, 1993). Estes últimos chegaram a

implementar rotinas para o cômputo de escoamento em reservatórios compostos de duas

regiões homogêneas, sendo uma delas interna a outra.

9 Ambos são regidos por equação diferencial parcial de mesma espécie.


6

Figura 1.3 - "Fio condutor” desta revisão bibliográfica.

1.3. Objetivos e Conteúdos

Os principais objetivos desta dissertação correspondem ao estudo do uso sistemático de

funções de Green para a solução de problemas de escoamento de fluidos pouco compressíveis

em meios porosos e à solução numérica10, através do Método dos Elementos de Contorno

10 Através de programas computacionais, escritos em Visual Basic for Applications (VBA), tendo como ponto de partida os disponibilizados por BREBBIA(1980) em seu livro fundador, formulados para problemas de potencial em regime permanente.

tempo

1959 1973 1991 1992 1993

CARSLAW e JAEGGER: codificação e sistematização do uso de soluções fundamentais em solução de problemas de condução de calor.

GRINGARTEN e RAMEY: transposição dos conceitos para o campo de escoamento em meios porosos.

1984

BREBBIA, TELLES e WROBEL: investigação, ampliação e consolidação do conhecimento sobre o MEC à época. Publicação de livro com parte dedicada à solução da equação da difusão, particularmente.

HORNE e outros: aplicação do MEC para solução de problemas de escoamento em meios porosos.

OZKAN e RAGHAVAN: obtenção e organização de soluções da equação da difusão baseadas em funções de Green no espaço de Laplace.

1988


7

(MEC), da equação da difusão hidráulica, em sua forma integral, quando submetida às

seguintes premissas11:

• domínio bidimensional, dividido em múltiplas (sub-)regiões homogêneas (rocha e

fluido)12;

• escoamento monofásico de fluido pouco compressível;

• contorno de geometria qualquer, sujeito a condição de queda de pressão ou de fluxo

prescrito;

• uma ou mais fontes presentes no domínio, produzidas por esquema de vazões

variáveis, sem inclusão de efeitos locais13;

• solução via espaço de Laplace, com inversão pelo tradicional e consagrado14 algoritmo

de STEHFEST (1970), resumido no Apêndice A.

No capítulo 2 é apresentado todo o ferramental físico-matemático envolvido na completa

definição analítica do problema supracitado.

A solução numérica via Método dos Elementos de Contorno e aspectos de sua

implementação computacional são descritos e explorados no capítulo 3.

No capítulo 4 é demonstrada a validade dos códigos computacionais gerados, e, por

extensão, seus limites, através do confronto com soluções analíticas já consagradas na

literatura especializada15 (veja-se EARLOUGHER JR., 1977; RAGHAVAN, 1993 e

BOURDET, 2001).

Finalmente, no capítulo 5 são oferecidas conclusões acerca dos resultados obtidos e

recomendações para futuros trabalhos.

11 Tipicamente, problemas semelhantes ao apresentado na Figura 1.4. 12 Esta dissertação é, portanto, uma expansão do estudo desenvolvido por KIKANI e HORNE (1992, 1993). 13 Quedas de pressão registradas a partir de produção de poços de petróleo estão sujeitas a dois efeitos locais: de estocagem e de película. O primeiro dos quais está relacionado à descompressão do fluido no interior da coluna, no momento imediatamente posterior a abertura do poço. O segundo é associado a danos à formação causados em fases anteriores à produção propriamente dita. Estes efeitos não fazem parte do escopo deste trabalho. 14 Largamente utilizado no campo da Engenharia de Petróleo. 15 Daqui por diante, sempre em que há referência à “literatura especializada” significa a relativa à análise de testes em poços de petróleo.


8

Figura 1.4 - Domínio típico de interesse.

1.4. Acerca da Nomenclatura Adotada

As teorias pertinentes ao escoamento de fluidos em meios porosos e ao MEC compõem

dois campos de conhecimento maduros e bem estruturados, que, por conseguinte, possuem

nomenclaturas próprias, ajustadas às culturas e às necessidades particulares de cada qual.

A exploração destes dois campos em conjunto, num mesmo documento, conduz, assim, a

certas divergências concernentes aos símbolos.

Neste aspecto, nesta dissertação, diferente do que é geralmente encontrado em livros de

Engenharia de Petróleo, a queda de pressão é tratada como u ao invés de p∆ .

Por outro lado, diversamente do que aparece nos livros de MEC, u , ou qualquer símbolo

que contenha a barra superior, está relacionado à variável no espaço de Laplace. Os valores

conhecidos no contorno são representados por “^”. Assim, ( ) uu ˆ=x significa que a queda de

pressão na posição x é especificada.

Os termos ponto-fonte (source point) e ponto-campo (field point) se referem a posição em

que se encontram, respectivamente, a fonte - solução fundamental para o problema em

x

y

w1

w2

u

n

u

∂∂ ˆ

n

u

∂∂ ˆ

u)

u

n

u

∂∂ ˆ

( )3,,, φµ tck

( )2,,, φµ tck

( )1,,, φµ tck

u)

n

u

∂∂ ˆ

w3

w4


9

análise, representada por ξ - e a variável de integração que percorre os elementos de contorno,

especificada por x.

As soluções fundamentais são classificadas de acordo com sua espécie (geometria): fonte

planar (plane source), fonte linear (line source) e fonte pontual (point source).

Os poços de petróleo são aludidos como fontes internas ao domínio.

CAPÍTULO 2: RELAÇÕES BÁSICAS

10

2.Relações Básicas

2.1. Conceitos Matemáticos Essenciais

Dois conceitos de cálculo vetorial são de fundamental importância para a consumação

do presente estudo e, em decorrência, merecem menção. As simples expressões que se

seguem constituem alicerce de praticamente todo o desenvolvimento subseqüente.

O teorema da divergência16 enuncia que se F é um campo vetorial de classe C1 para o

domínio Ω , limitado por uma superfície fechada Γ , orientada por vetores normais

apontados para exterior, então:

Γ⋅=Ω⋅∇ ∫∫ΓΩ

dd nFF . (2.1)

O símbolo ∇ representa o operador de Hamilton, que para um domínio cartesiano

tridimensional é definido como:

kjizyx ∂

∂+∂∂+

∂∂=∇ , (2.2)

sendo i, j e k os vetores unitários nas direções x,y e z, respectivamente.

O produto escalar (.) entre ∇ e F mede a divergência vetorial em uma determinada

posição do domínio.

Assim, a integral ao longo de todo o domínio do operador divergente aplicado a F

corresponde à integral de superfície da projeção de F na direção de seu vetor normal n,

conforme Figura 2.1.

16 Também conhecido como teorema de Gauss, teorema de Ostrogradski ou teorema de Gauss-Ostrograski. Veja-se, por exemplo, PISKUNOV (S/A).


11

Figura 2.1 – Ilustração para o teorema da divergência.

Um segundo conceito igualmente necessário e importante é a regra de

multiplicação de derivadas aplicada a domínios multidimensionais. Por exemplo, se

sobre determinada região estão definidas as funções escalares ( )xu e ( )xw , uma das

possibilidades de uso de tal regra é:

( ) uwuwuw 2∇+∇⋅∇=∇⋅∇ . (2.3)

2.2. Operações no Espaço de Laplace

O problema de interesse do presente estudo é de caráter transitório e, como

conseqüencia, sua solução, analítica ou numérica, é naturalmente dependente do tempo.

Em certas circunstâncias, o banimento de tal dependência pode ser vantajoso. As

transformadas de Laplace, por meio da expressão abaixo, oferecem tal oportunidade.

( ) ( ) ( )dttfesftfL st ,,,0

xxx ∫∞

−== . (2.4)

Com efeito, a integral acima redunda no desaparecimento da variável temporal, t. A

função transformada passa a estar relacionada com a variável s. É caro notar que a parcela

Ω

n F

Γ

F.n


12

da função relacionada às coordenadas espaciais, x, não sofre nenhuma alteração quando da

passagem do espaço original para o espaço transformado.

Entretanto, para que haja transformada de Laplace, evidentemente, é preciso que ( )tf ,x

seja seccionalmente contínua no intervalo ∞<< t0 e que os limites para 0→t e para

∞→t existam e possam ser calculados.

No espaço de Laplace, algumas operações matemáticas tidas, em variáveis tempo-

espaço, como de algum grau de complexidade, estão correlacionadas a simples operações.

Por exemplo, a diferenciação com relação ao tempo corresponde a:

( ) ( ) ( )+−=

0,,, xxx fsfst

dt

dfL . (2.5)

Inversamente, para a integração no intervalo [ ]t,0 , se tem que:

( ) ( )sfs

dfLt

,1

,0

xx =

∫ ττ . (2.6)

A convolução, outra operação de grande valia em alguns campos da física-

matemática17, também sofre notável simplificação.

( ) ( ) ( ) ( ) ( ) ( )sgsfdtgtfLtgtfLt

,,,.,,*,0

xxxxxx =

−= ∫ ττ . (2.7)

Além disso, cumpre notar que as transformadas de Laplace respeitam propriedade de

linearidade. Isto é, se a e b são constantes então:

( ) ( ) ( ) ( )sgbsfatbgtafL ,,,, xxxx +=+ . (2.8)

Com a facilidade propiciada pelas propriedades aludidas, naturalmente os cálculos

algébricos no espaço transformado tendem a se traduzir em operações mais simples.

Contudo, ao término de tais operações é necessário que se converta a resposta

novamente para o campo real de origem. Esta transformação inversa, caso não seja

encontrada sob a forma tabelada, deve ser calculada pela expressão:

17 Em particular, para o desenvolvimento desta dissertação.


13

( ) ( ) ( )∫∞+

∞−

− ==ia

ia

st dssfei

tfsfL ,2

1,,1 xxx

π. (2.9)

Como é patente, tal avaliação envolve integração sujeita a limites complexos18.

Por outra, em casos que tal conversão seja problemática ou não possa ser obtida de

forma explícita, comumente se lança mão de algoritmos numéricos. Dentre tais, destaca-se

o de STEHFEST (1970), representado pela seguinte equação:

( )

=≈ ∑=

it

sfVt

tfNSteh

ii

2ln,

2ln,

1

xx . (2.10)

Em poucas palavras, uma aproximação da função ( )tf ,x no espaço real pode ser obtida

através da avaliação em NSteh pontos de ( )sf ,x no espaço transformado. É um

procedimento de certa forma similar ao de integração numérica por Gauss-Legendre19.

Os coeficientes iV dependem exclusivamente do valor NSteh adotado e, por extensão,

podem ser tabelados (ver Apêndice A).

Maiores detalhes analíticos a respeito de transformadas de Laplace podem ser

amplamente encontrados em ABRAMOVITZ e STEGUN (1972) e CHURCHILL (1972).

No Apêndice C deste documento são apresentadas as transformadas utilizadas neste estudo.

2.3. Equações de Dinâmica dos Fluidos

O escoamento de fluidos pouco compressíveis em meios porosos é regido pela equação

da difusão hidráulica, correspondente à associação dos conceitos e expressões a saber:

• Lei da Conservação de Massas;

• Lei de Darcy e

• Equações de estado, dependentes das espécies de fluido e rocha sob

consideração.

18 A rigor, de forma geral s é complexa. 19 No sentido em que a aproximação da integral correspondente a transformação inversa é realizada através de um somatório baseado em pesos e pontos amostrais pré-definidos. A semelhança é somente esta. O procedimento de Gauss-Legendre é utilizado para aproximação de integrais espaciais. De forma diversa, o de STEHFEST(1970), para transformadas temporais.


14

2.3.1. Lei da Conservação de Massas (ou Equação da Continuidade)

Seguindo-se o procedimento de PISKUNOV (S/A): seja um volume, Ω , cuja

porosidade possa ser representada por ( )t,xφ 20 e a permeabilidade por ( )xk , limitado por

uma superfície Γ , descrita através do vetor normal unitário a ela própria externamente

orientado ( )xn , saturado, por simplificação, de um único fluido móvel de massa específica

( )t,xρ que se desloca à velocidade ),( txv . Admita-se, também, a presença de fontes (ou

sumidouros) dispersas no domínio, representadas pelo conjunto de suas vazões específicas

volumétricas ( )tq w ,~ x . Seja conhecida a massa que este sistema porta num instante de

tempo t.

Considerando-se que a influência ocasionada por algum presumível campo

gravitacional possa ser desprezível para fins de cálculo, pode-se estabelecer, pela lei da

conservação de massas, que:

0=∆+∆+∆=∆∑ ρMMMM

we. (2.11)

Ou seja, se nenhuma massa é criada ou extinta durante um período de tempo t∆ ,

necessariamente, a soma das diversas variações consideradas deve ser nula. Portanto, ao

balanço de entrada e saída de fluidos pelas fronteiras do corpo, e

M∆ , e pelas fontes

possivelmente existentes no domínio, w

M∆ , corresponde a variação de massa no interior

do corpo, ρ

M∆ .

O primeiro termo da equação (2.11) pode ser calculado por:

dtdMtt

te

Γ⋅=∆ ∫ ∫∆+

Γ

nvρ , (2.12)

onde Γd se refere a um elemento infinitesimal de superfície.

Pode-se verificar, por simples inspeção de coerência entre as unidades, que o resultado

da integral acima reflete efetivamente uma quantidade de massa.

20 A porosidade representa um papel curioso. Em princípio é considerada como função dependente da queda de pressão e, por conseguinte, do tempo. Mais adiante, essa dependência é abandonada, para que a equação da difusão hidráulica seja tornada linear.


15

A aplicação do teorema da divergência - expressão (2.1) - à equação (2.12) redunda em:

( ) dtdMtt

te

Ω⋅∇=∆ ∫ ∫∆+

Ω

vρ . (2.13)

A variação de massa devida à concorrência de fluidos provenientes de fontes internas

ao domínio pode ser calculada diretamente por:

dtdqMtt

tw

Ω=∆ ∫ ∫∆+

Ω

~ρ . (2.14)

Com efeito, a vazão específica volumétrica, q~ , pode ser transformada em vazão

mássica através de sua simples multiplicação pela massa específica do fluido em objeto,ρ .

Naturalmente, a integral deve ser calculada em porções do domínio, Ω , onde haja fontes. É

importante ressaltar, ainda, que este resultado será positivo quando a taxa de injeção for

maior do que a de retirada de fluidos21.

Por seu turno, a variação da quantidade de massa no interior do domínio pode ser

expressa por:

( ) dtdt

Mtt

t

Ω∂∂=∆ ∫ ∫

∆+

Ω

ρφρ22,23. (2.15)

Finalmente, com o estabelecimento da igualdade (2.11) chega-se a:

( ) ( ) 0~ =Ω∂∂+Ω+Ω⋅∇ ∫ ∫∫ ∫∫ ∫

∆+

Ω

∆+

Ω

∆+

Ω

dtdt

dtdqdtdtt

t

tt

t

tt

t

ρφρρv . (2.16)

Ou, agrupando-se as integrais:

21 Freqüentemente e por simples convenção, na literatura especializada, se encontra a retirada como positiva. Mais adiante, o sinal desta variação está de acordo com a orientação do vetor normal em relação à superfície considerada. 22 A variação associada à porção sólida (rochosa) é habitualmente negligenciada. 23 A rigor, para o cálculo dessa quantidade só é necessário o conhecimento prévio das massas internas inicial e

da final do sistema, posto que: ( ) ( ) ( ) Ω−=Ω∂∂

∫∫ ∫Ω

∆+Ω

∆+

dddtt ttt

tt

t

ρφρφρφ .


16

( ) ( ) 0~ =Ω

∂∂++⋅∇∫ ∫

∆+

Ω

dtdt

qtt

t

ρφρρv . (2.17)

Mas como tanto o intervalo de tempo considerado nos cálculos quanto a forma a que se

restringe o domínio de avaliação são arbitrários resulta que o integrando de (2.17) deve ser

nulo. Portanto:

( ) ( ) 0~ =∂∂++⋅∇ ρφρρt

qv . (2.18)

A equação (2.18) é a lei de conservação de massas para escoamento de fluidos pouco

compressíveis em meios porosos.

Cumpre notar, ademais, que para a completa definição do problema existe a

necessidade de equações complementares que correlacionem a velocidade de fluxo, a massa

específica do fluido móvel e a porosidade do meio com a grandeza física passível de

medição em campo24 – a queda de pressão.

2.3.2. Lei de Darcy

A lei de Darcy, que estabelece a relação, para um meio poroso, entre a velocidade

macroscópica, v, e o potencial de fluxo entre duas posições, Φ , pode ser expressa, em

termos gerais, por:

Φ∇−= kvµρ

, (2.19)

onde os parâmetros µ e ρ são, respectivamente, a viscosidade dinâmica e a massa

específica do fluido em movimento. O sinal negativo indica que o escoamento se dá da

posição de maior potencial para a de menor.

O tensor de segunda ordem k corresponde às permeabilidades do meio poroso para o

caso geral anisotrópico, conforme a relação:

24 Ao lado da vazão das fontes (ou sumidouros).


17

=

zzzyzx

yzyyyx

xzxyxx

kkk

kkk

kkk

k . (2.20)

O potencial de fluxo, segundo RAGHAVAN (1993), pode ser expresso por:

( )2

' 2

00

mp

p

vzzg

dp +−+=Φ ∫ ρ, (2.21)

onde os valores 0p e 0z correspondem, respectivamente, à pressão e à cota num nível de

referência. A velocidade média microscópica é representada por mv e a aceleração da

gravidade por g. Os três termos do segundo membro da igualdade acima estão relacionados,

nesta ordem, a: energia decorrente de forças viscosas, energia potencial e energia cinética.

Na grande maioria dos problemas de escoamento de fluidos pouco compressíveis em

meios porosos e na totalidade dos casos analisados neste trabalho as forças viscosas

dominam completamente a resposta, sendo, portanto, natural se desprezar o 2º e o 3º termos

de (2.21).

Utilizando-se o teorema do valor médio e tendo-se em mente as considerações acima

expostas, o gradiente da expressão (2.21) fica reduzido a:

u∇≈Φ∇ρ~1

, (2.22)

com ρ~ respresentando uma massa específica média e com a queda de pressão sendo

definida como

0ppu −= 25, (2.23)

onde 0p é a pressão inicial do reservatório.

Substituindo-se a expressão (2.22) na equação (2.19) e admitindo-se, em adição,

existência de pequenos gradientes de pressão no meio poroso, chega-se à lei de Darcy

generalizada num dos formatos de interesse:

25 Note-se que tal definição é coerente com a de injeção de fluido positiva, citada anteriormente.


18

u∇−=µk

v . (2.24)

A expressão acima, escrita em termos matriciais, para um domínio bidimensional, por

exemplo, é:

∂∂∂∂

−=

yu

xu

kk

kk

v

v

yyyx

xyxx

y

x

µ1

. (2.25)

Evidentemente, se o meio poroso é homogêneo e isotrópico26, o tensor k pode ser

tratado como uma grandeza escalar. Assim sendo, a equação (2.24) também pode ser escrita

como:

uk ∇−=µ

v . (2.26)

Considerando-se, ainda, que vazão infinitesimal, dq, e velocidade de um fluido que

perpassa um elemento de área infinitesimal,wdΓ , se relacionam por

wddq Γ⋅= nv , (2.27)

a expressão (2.26) ainda pode assumir a forma

∫Γ

Γ⋅∇−=w

wduk

q nµ

. (2.28)

Esta expressão determina a vazão de fluido que atravessa a superfície delimitada por wΓ

- por exemplo, a de um poço de petróleo.

Se, adicionalmente, a transposição de fluido a esta superfície se dá de maneira

uniforme, a expressão (2.28) ainda pode ser simplificada para27:

26 Ou mesmo ortotrópico, com a devida transformação no sistema de coordenadas.

27 Observe-se que n

uu

∂∂=⋅∇ n e µφη tck= .


19

wn

u

c

q

t Γ∂∂−= η

φ

~. (2.29)

A vazão específica (ou vazão por unidade de área de fonte interna) deve ser, neste caso,

definida por:

∫Γ

Γ=w

wdqq~ . (2.30)

Por motivos elementares, na igualdade (2.29) a derivada nu ∂∂ deve ser avaliada na

superfície interna ao domínio considerada. O parâmetro tc representa a compressibilidade

total do sistema e é definido adiante.

2.3.3. Equações de Estado

fluido

A compressiblidade isotérmica28 de um fluido é definida como a variação relativa entre

seu volume, V, e pressão a ele aplicada, p, e pode ser expressa por:

Tp

V

Vc

∂∂−= 1

, (2.31)

ou, ainda, considerando a definição (2.23) e que 1=∂∂ pu :

Tu

V

Vc

∂∂−= 1

. (2.32)

O sinal negativo indica que uma diminuição da pressão aplicada induz um aumento do

volume de fluido e vice-versa.

Da definição de massa específica de um fluido se tem que:

ρm

V = . (2.33)

28 As variações de temperatura durante a operação de um poço de petróleo são habitualmente ignoradas. Todas as propriedades consideradas, por simplificação, como constantes, são avaliadas à pressão e à temperatura iniciais do reservatório.


20

Sendo a massa, m, constante, a derivada da igualdade (2.33) deve ser calculada como:

u

m

p

V

∂∂−=

∂∂ ρ

ρ 2. (2.34)

Finalmente, a substituição das expressões (2.33) e (2.34) em (2.32) fornece:

Tuc

∂∂= ρ

ρ1

. (2.35)

Portanto, a compressibilidade de um fluido também pode ser calculada pela variação

relativa entre sua massa específica e a queda de pressão a ele aplicada.

Rocha

Analogamente, a compressibilidade isotérmica efetiva de uma rocha reservatório, cr, é

definida por:

T

r p

V

Vc

∂∂= 1

, (2.36)

ou, pelos mesmos motivos usados na definição da equação (2.32),

Tr u

V

Vc

∂∂= 1

. (2.37)

Todavia, diferente do que ocorre em (2.31), a relação entre variação de pressão e

variação de volume é direta. Isto é, a diminuição da pressão aplicada à rocha conduz a

diminuição do espaço poroso.

Por definição, o volume poroso de uma rocha relaciona-se ao seu volume total, Vt, por:

tVV φ= . (2.38)

Sendo Vt aproximadamente constante, a derivada de (2.38) em relação à queda de

pressão vale:

uV

u

Vt ∂∂≈

∂∂ φ

. (2.39)


21

Pela substituição de (2.38) e (2.39) em (2.32) vem:

ucr ∂

∂≈ φφ1

. (2.40)

2.4. Equação da Difusão Hidráulica

Conforme mencionado anteriormente, a equação da difusão hidráulica resulta da

composição entre as equações abaixo resumidas.

• Equação da continuidade

( ) ( ) 0~ =∂∂++⋅∇ ρφρρt

qv ,

• Lei de Darcy

u∇−=µk

v ,

• Equações de estado

Tuc

∂∂= ρ

ρ1

e

ucr ∂

∂≈ φφ1

.

Como tal, carrega em seu seio todas as hipóteses simplificadoras anteriormente

relacionadas.

A obtenção de sua forma sintética desejada pode ser estabelecida através do

procedimento abaixo descrito.

Primeiramente, substitui-se (2.26) em (2.18) e transfere-se o terceiro termo para o

segundo membro, resultando em:

( )ρφρρµ t

qu∂∂−=+

∇⋅∇− ~k. (2.41)


22

Expandindo-se as derivadas dos produtos de ambos os lados da equação anterior e

multiplicando-a inteiramente por (-1), vem:

ttquu

∂∂+

∂∂=−∇⋅∇+

∇⋅∇ ρφφρρρµµ

ρ ~kk. (2.42)

Dividindo-se (2.42) por ρ , tendo em conta que uu∇∂∂=∇ ρρ , e colocando-se φ

em evidência no segundo membro, chega-se a:

( )

∂∂+

∂∂=−∇

∂∂+

∇⋅∇tt

qup

uρ

ρφ

φφρ

ρµµ11~1 2kk

. (2.43)

Como φ e ρ dependem da pressão, as seguintes igualdades se aplicam:

∂∂

∂∂=

∂∂

∂∂

∂∂=

∂∂

.t

u

ut

t

u

utρρ

φφ

(2.44)

Tendo em vista as transformações acima e introduzindo-se as equações de estado (2.35)

e (2.40) vem:

( )t

ucqucu t ∂

∂=−∇+

∇⋅∇ φµµ

~2kk (2.45)

Onde ct é a compressibilidade total do sistema e seu valor resulta da soma da

compressibilidade do fluido com a da rocha29.

rt ccc += (2.46)

Convém notar que o segundo termo da equação (2.45) pode ser desprezado, pois

corresponde à multiplicação entre duas quantidades pequenas em relação às demais: a

compressibilidade de um fluido pouco compressível e o quadrado do gradiente da queda de

pressão.

29Para o caso de um único fluido saturando os poros da rocha.


23

Ademais, admitindo-se que a viscosidade do fluido possa ser considerada

aproximadamente constante30 e com a divisão de toda a equação (2.45) por tcφ , chega-se a

uma versão de interesse da equação da difusão hidráulica, a seguir:

( )t

u

c

qu

c tt ∂∂=−∇⋅∇

φφµ

~1k . (2.47)

Em adição, se o domínio de interesse é regionalmente homogêneo e ortotrópico e se a

orientação cartesiana dos eixos coordenados coincide com as direções principais de

permeabilidade, a equação (2.47) pode ainda assumir a forma:

t

u

c

q

y

u

c

k

x

u

c

k

tt

y

t

x

∂∂=−

∂∂+

∂∂

φφµφµ

~2

2

2

2

. (2.48)

Por fim, pode-se mostrar com certa facilidade, que se se operar em (2.48) uma

transformação de coordenadas do tipo ( )

→ y

k

kx

k

kyx

yx

,, , com k sendo um valor de

referência (por conveniência: yxkkk = ) esta expressão será conduzida a forma

correspondente ao caso isotrópico (ver Apêndice B). Ou seja,

t

u

c

qu

t ∂∂=−∇

φη

~2 , (2.49)

com

tc

k

φµη = (2.50)

sendo a constante de difusão hidráulica.

Caso não haja atuação de fonte ou sumidouro no domínio do problema31, a expressão

(2.49) assume a compacta forma:

30 Essa simplificação é de crucial importância para a linearidade da equação. 31 O escoamento em reservatórios de petróleo é causado pela produção de (injeção em) poços. Deste modo, na prática, sempre há sumidouros (fontes) presentes no domínio. A incorporação matemática de tais elementos pode se dar pela consideração da equação (2.49), diretamente, ou por condição de contorno aplicada à expressão (2.51), de forma indireta.


24

02 =∂∂−∇

t

uuη . (2.51)

As expressões (2.47), (2.49) e (2.51) são, notadamente, equações diferenciais parciais

de 2ª ordem.

Com efeito, a variável de interesse quando da solução da equação de difusão hidráulica

– a queda de pressão - é dependente do tempo e do espaço.

2.5. Condição Inicial, Condições de Contorno e Condições de

Interface Empregadas

2.5.1. Condição Inicial

Na totalidade dos casos práticos investigados neste trabalho, a pressão inicial do

reservatório é considerada constante e uniforme para um nível especificado de referência32.

As causas que fundamentam tal escolha são basicamente três:

a) Ora, quando um poço é posto pela primeira vez em produção num domínio ainda

virgem, é justamente nesta condição que o parâmetro é encontrado. É suposto que o

reservatório encontre-se em equilíbrio, atingido após longo processo em escala geológica;

b) As posições dos poços e seus históricos de produção são – ou, ao menos deveriam ser –

registrados com toda a diligência presumível. Com essas informações, a partir da

superposição dos efeitos no tempo e no espaço33 , caso se conheça com exatidão o sistema,

é possível obter o potencial em posição e tempo especificados;

c) Ademais, não é plausível se proceder de outra forma, posto que a única maneira de se

acessar ( )tu ,x é através de registros em poços.

Deste modo, a condição inicial é expressa por:

( ) 00, ptp ==x (2.52)

Ou ainda,

32 Cota ou profundidade. 33 Como já se viu, a equação da difusão hidráulica, com as hipóteses simplificadoras adotadas, constitui-se numa expressão linear e como tal admite superposição dos efeitos.


25

( ) ( ) 00,0, 0 =−=== ptptu xx . (2.53)

2.5.2. Condições de Contorno

As condições de contorno formuladas para fins deste trabalho são:

a) Queda de pressão prescrita (condição de contorno de Dirichlet)

( ) utu e ˆ, =x ; (2.54)

b) Derivada da queda de pressão na direção normal ao contorno prescrita (condição de

contorno de Neumann)

( ) ( )n

tu

n

tu ee

∂∂=

∂∂ ,ˆ, xx

. (2.55)

Destas, a mais largamente utilizada na prática é ( )

0,ˆ

=∂

∂n

tu ex, para representação dos

limites físicos do reservatório.

Quando há presença de um aqüífero de potência considerável nas cercanias da formação

também é comum se lançar mão de ( ) 0, =tu ex na porção correspondente a tal ocorrência.

2.5.3. Condições de Interface

Caso o meio poroso seja dividido em regiões de propriedades distintas, duas condições

de continuidade devem ser asseguradas em suas interfaces. Quais sejam:

a) De compatibilidade entre as quedas de pressão – ou seja, pressões na interface,

calculadas em função dos parâmetros de cada região a que ela pertença, devem ser

idênticas;

b) De equilíbrio entre os fluxos – os fluxos na interface, calculados em função dos

parâmetros de cada região a que ela pertença, devem ter mesma magnitude e sentidos

opostos34.

34 Considerando-se que as regiões correspondentes a interface, avaliadas de forma isolada, sejam orientadas por vetores normais unitários apontados para o exterior.


26

2.6. Superposição no Espaço

O Princípio da Superposição no Espaço preconiza que a queda de pressão numa posição

x de um reservatório de petróleo corresponde à soma entre os efeitos causados, naquela

posição, pela produção dos poços operados no domínio e pelas barreiras ao fluxo ali

existentes, representados genericamente por suas respectivas localizações xj . Ou,

matematicamente:

( ) ( )∑=

−=N

jj tutu

1

,, xxx , (2.56)

onde N representa o total de ocorrências35.

Método das Imagens O Método das Imagens constitui-se num artifício matemático freqüentemente adotado

para o cálculo de limites de reservatório.

Por exemplo, seja um poço de petróleo produzido à vazão q, situado a uma distância L

de uma barreira selante de extensão infinita para o tempo considerado nos cálculos, tal qual

a representada na Figura 2.2a.

Figura 2.2 - Aplicação do Método das Imagens.

35 Poços e barreiras ao escoamento.

a) sistema real b) sistema equivalente


27

De acordo com o Método das Imagens, a queda de pressão numa posição qualquer do

sistema apresentado na Figura 2.2a pode ser calculada através do sistema equivalente

revelado na Figura 2.2b.

Tal assertiva é baseada na configuração das idênticas linhas de corrente desenvolvidas

nos dois sistemas. É importante notar que um poço imaginário, posicionado a 2L do poço

real, também deve ser considerado como produzido à taxa q.

Assim sendo, no caso em tela, lançando-se mão do Método das Imagens, a queda total

de pressão num ponto x qualquer pode ser calculada por (2.56), tomando-se N=2 (1 poço

real e 1 imaginário), neste caso.

Seguindo-se o mesmo raciocínio, a queda de pressão numa posição x, derivada da

produção de um poço situado nas cercanias de uma dupla de barreiras que perfazem um

ângulo reto entre si (ver Figura 2.3a e b), também pode ser calculada através de (2.56),

sendo, desta vez, N=4 (1 poço real e 3 imaginários).

Ora, por extensão deste exercício ou por simples inspeção é possível perceber que a

generalização do raciocínio correspondente aos exemplos anteriores implica em que a

quantidade de termos da somatória em (2.56), N, pode ser calculada por:

θπ2=N , (2.57)

sendo θ o ângulo interno da dupla de barreiras (ver Figura 2.3c).

Figura 2.3 - Barreiras ao fluxo: em ângulo de 90º e em ângulo qualquer.

a) a 90º: sistema real b) a 90º: sistema equivalente c) em ângulo qualquer


28

Para a consideração da expressão (2.56), o valor N deve ser obrigatoriamente inteiro36.

2.7. Superposição no Tempo (Convolução)

Seguindo-se o procedimento de RAGHAVAN (1993): seja um reservatório de

petróleo isotrópico, inicialmente em equilíbrio37, submetido à produção de fluido a

partir de um único poço, regida pela equação da difusão hidráulica, conforme

( ) ( )0

,,2 =

∂∂−∇

t

tutu

xxη , (2.58)

ou, no espaço transformado de Laplace,

( ) ( ) 0,,2 =−∇ sussu xxη . (2.59)

Admitindo-se, inicialmente, que a referida produção se dê por taxa constante, no

próprio poço deve ser observada a condição dada pela lei de Darcy (2.29), que

transformada para o espaço de Laplace pode ser escrita como:

( )w

n

suk

s

q cc

Γ∂∂

−=,~ x

µ. (2.60)

Admita-se que se conheça a solução do problema até aqui exposto, representada

por ( )suc ,x , com o índice c fazendo referência à vazão constante.

Se ( )suc ,x é solução do problema descrito, conseqüentemente ( ) ( )susu c ,, xx =

respeita a igualdade (2.59).

Agora, por outra, imaginando-se um esquema de vazões variáveis em lugar da

produção à taxa constante, a condição de contorno no poço passa a:

( ) ( )w

n

suksq

Γ∂∂−= ,~ x

µ. (2.61)

Para este segundo caso, pode-se provar que com

36 Neste caso, as possibilidades de θ se restringem aos divisores de π2 . 37 Portanto, com ( ) 00, ==tu x .


29

( ) ( ) ( )sussqq

su cc

,~~1

, xx = (2.62)

a equação diferencial (2.59) é satisfeita38. Isto é, a expressão (2.62) é solução da

igualdade (2.59), para a condição de contorno (2.61).

Assim, a resposta de um sistema submetido a vazões variáveis se relaciona com a

com a resposta deste mesmo sistema operado à taxa constante por (2.62).

Ademais, observando-se que

( ) ( )t

tususL c

c ∂∂=− ,

,1 xx , se ( ) 00, ==tuc x , (2.63)

chega-se ao teorema da convolução39,40 no espaço real original,

( ) ( ) ( ) ( ) ( ) τττdt

t

u

q

q

t

tu

q

tqtu

tc

c

c

c

−∂

∂=∂

∂= ∫ ,,

*,0

xx

x . (2.64)

Portanto, a queda de pressão num ponto qualquer de um reservatório de petróleo,

num dado instante, pode ser calculada pela convolução – representada pelo símbolo (*)

- entre a derivada temporal da resposta do sistema à vazão constante e o histórico de

vazões considerado.

Se o tal histórico de vazões varia sob a forma de N degraus41 então a expressão

(2.64) pode ser simplificada para:

( ) ( ) ττ dtt

uq

qtu

j

j

t

t

cN

jj

c

−∂

∂= ∫∑−

=

,1

,1

1

xx . (2.65)

38 Curiosamente, o quociente entre (2.61) e (2.60) seguido de rearranjo dos termos produz:

( ) ( ) ( )n

sussq

qn

su c

c ∂∂

=∂

∂ ,~~1, xx

.

39 Também conhecido como Princípio de Duhamel. 40 Note-se que ( ) ( ) cc qtqqtq =~~ . 41 Geralmente as mudanças substanciais de vazão em poços de petróleo são ocasionadas por troca de estranguladores de fluxo (chokes) alocados nas linhas de produção. Como estes elementos têm diâmetro de passagem discretos, sua troca acarreta mudança brusca de vazão. Portanto, a vazão costuma variar em degraus.


30

A integral anterior pode ser prontamente calculada, com o artifício de uma

transformação de variáveis do tipo:

ατ =−t . (2.66)

Neste caso, a expressão (2.65) passa a ser:

( ) ( ) ( )[ ]∑=

− −−−=N

jjcjcj

c

ttuttuqq

tu1

1 ,,1

, xxx , (2.67)

com ( ) 0=− ttuc , evidente.

Por outra, expandindo-se o somatório e reagrupando-se os termos passa-se a ter que:

( ) ( )∑=

−− −−=N

jjcjj

c

ttuqqq

tu1

11 ,1

),( xx , (2.68)

com 00 =q , N sendo o período em exame e j sendo o contador dos períodos (degraus).

2.8. Soluções Analíticas para Escoamento de Fluidos em

Meios Porosos Baseadas em Funções de Green

2.8.1. Conceitos Preliminares – Delta de Dirac

O delta de Dirac, δ , ilustrado na Figura 2.4, é definido como um impulso unitário e

instantâneo.

Figura 2.4 – delta de Dirac.

0 t,τ

δ ( ) 10

=−∫∞

ττδ dt

ε

ε1

0→ε


31

No domínio do tempo, sua principal propriedade é

( ) ( ) ( )tfdtft

∫ =−0

ττδτ . (2.69)

Adicionalmente, é interessante notar pela simples inspeção da Figura 2.5, que

( ) ( ) ( )tfdtft

∫+

=−ε

ττδτ0

(2.70)

e

( ) ( ) 00∫−

=−ε

ττδτt

dtf . (2.71)

Figura 2.5 – Convolução com o delta de Dirac.

Analogamente a (2.69), para o domínio espacial se pode escrever que:

( ) ( ) ( )ξxξx fdf =Ω∫Ω

,δ . (2.72)

Além disso, pela combinação de (2.69) e de (2.72), para um domínio espaço-temporal

se tem:

0

ε+t t τ

ε−t

( )τf

τ−t

( )τδ −t


32

( ) ( ) ( ) ( )tfddtft

,,,0

ξxξx =Ω−∫ ∫Ω

ττδδτ . (2.73)

A própria definição do delta de Dirac neste último caso, estabelece que:

( )

≠−≠==−=∞→

−+

+

0,0

0,,,

ττ

τtou

tetf

ξx

ξxxξ . (2.74)

Onde o índice “+” sobre o zero indica que o impulso já tenha sido efetivado42.

2.8.2. Funções de Green

Sem muito rigor matemático, seja uma equação diferencial não homogênea típica das

manipuladas no presente estudo:

( ) ( )tftLu ,, xx = , (2.75)

com L representando um operador linear e u e f funções escalares quaisquer, dependentes

de tempo e espaço.

A função de Green para (2.75), ( )τ−tg ,,xξ , de interesse para este trabalho, é aquela

que honra a seguinte igualdade:

( ) ( ) ( )τδξδτξ −−=− ttgL xx ,,,* , (2.76)

com L* sendo operador diferencial adjunto43 a L e δ o delta de Dirac.

Convém notar que ( )τ−tg ,,xξ também é um impulso e, como tal, respeita a

propriedade (2.73). Em conseqüência:

( ) ( ) ( )tfddtgft

,,,,0

ξxξx =Ω−∫ ∫Ω

τττ . (2.77)

Além do mais, se por definição o impulso é unitário, então para todo e qualquer t,

obrigatoriamente:

42 Até mesmo ( ) 00,, =−xξf . 43 Operador adjunto é aquele que permite estabelecer: gLugLu *,, = .


33

( ) 1,, =Ω∫Ω

dtg xξ .

As funções de Green, de fato, são distribuições e também são conhecidas como

soluções fundamentais.

2.8.3. Equação Integral Ponderada da Difusão – Ponto-fonte Pertencente

ao Domínio

Seja um domínio bidimensional isotrópico, Ω , limitado externa e internamente por

curvas representadas por eΓ e wΓ , respectivamente, com a fronteira externa denotando os

limites físicos do reservatório e a interna fazendo o papel de poços.

A Figura 2.6 é um desenho esquemático representativo da espécie de domínio em

questão. Note-se que as fronteiras externas podem estar sujeitas a dois tipos de condição de

contorno: queda de pressão prescrita ou derivada normal de queda de pressão prescrita. As

fronteiras internas (fontes) são admitidas exclusivamente como de fluxo prescrito.

Figura 2.6 – Domínio Considerado.

Seja o escoamento neste domínio governado pela equação diferencial (2.58), com o

tempo sendo aqui denotado por τ , ou seja

( ) ( )0

,,2 =

∂∂−∇

t

uu

ττη xx . (2.78)

ξξξξ

x


34

Ponderando-se (2.78) por uma função de natureza conhecida ( )τ,xw 44 e integrando-se

no domínio do problema passa-se a ter:

( ) ( ) ( ) ( ) 0,

,,0

2 =Ω

∂∂−∇∫ ∫

Ω

t

ddt

uuw τττητ x

xxx , (2.79)

ou separando-se os termos,

( ) ( ) ( ) ( ) ( ) ( ) 0,

,,,00

2 =Ω∂

∂−Ω∇ ∫ ∫∫ ∫ΩΩ

tt

ddt

uwdduw ττττττη x

xxxxx . (2.80)

Com o objetivo de transferir os operadores lineares da função incógnita, u, para a

função de ponderação, w, pode-se aplicar sucessivamente a regra de multiplicação de

derivadas, mencionada no capítulo anterior. Neste caso:

( ) uwuwuw ∇⋅∇−∇⋅∇=∇2 (2.81)

e

( ) wuwuuw 2∇−∇⋅∇=∇⋅∇ . (2.82)

Mormente, da combinação de (2.81) e (2.82) se tem que

( ) ( ) wuwuuwuw 22 ∇+∇⋅∇−∇⋅∇=∇ . (2.83)

Da mesma forma, para a derivada temporal,

( )t

wuwu

tt

uw

∂∂−

∂∂=

∂∂

. (2.84)

Aplicando-se, agora, a (2.80) as transformações (2.83) e (2.84) chega-se a:

( ) ( )[ ] ( ) ( )[ ] ( )

( ) ( ) ( ) ( )

( ) ( ) ( ) .0,,

,,,

,,,,

0

0

2

0

=Ω

∂∂−

−Ω

∂∂+∇+

+Ω∇⋅∇−∇⋅∇

∫ ∫

∫ ∫

∫ ∫

Ω

Ω

Ω

τττ

τττητ

τττττη

ddt

uw

ddt

wwu

ddwuuw

t

t

t

xxx

xx

xx

xxxxx

(2.85)

44 Função obrigatoriamente de classe igual ou superior a de u.


35

Observe-se que o 2º termo da expressão (2.85) contém a equação adjunta a (2.78). Para

que este termo possa ser avaliado de forma simples e direta convém admitir que ( )τ,xw

seja substituída pela função de Green para o problema em estudo, ( )τ−tg ,,xξ .

Assim sendo, como por definição a função ( )τ−tg ,,xξ satisfaz a equação adjunta ao

problema original (2.76) então aplica-se:

( ) ( ) ( ) ( )τδδττη −−=∂

−∂+−∇ tt

tgtg xξ

xξxξ ,

,,,,2 . (2.86)

Pela aplicação de ( )τ−tg ,,xξ em lugar de ( )τ,xw associada à consideração de (2.86),

a expressão (2.85) passa a ser expressa por:

( ) ( )[ ] ( ) ( )[ ] ( )

( ) ( ) ( )[ ] ( )

( ) ( ) ( ) .0,,,

,,

,,,,,,

0

0

0

=Ω

∂−∂−

−Ω−−+

+Ω−∇⋅∇−∇−⋅∇

∫ ∫

∫ ∫

∫ ∫

−

Ω

−

Ω

−

Ω

τττ

ττδδτ

τττττη

ε

ε

ε

ddt

utg

ddtu

ddtguutg

t

t

t

xxxξ

xxξx

xxξxxxξ

(2.87)

Note-se que para que seja evitada a avaliação das integrais temporais em condição de

singularidade45, como artifício matemático subtrai-se46 um curto período de tempo a seus

limites superiores, ε . Posteriormente, se extrai o limite para 0→ε .

O 1º termo de (2.87) pode ser calculado, de acordo com o teorema da divergência

(equação (2.1)), como:

( ) ( )[ ] ( ) ( )[ ] ( )

( ) ( ) ( ) ( ) ( ) .,,

,,

,,

,,,,,,lim

0

0

0

∫ ∫

∫ ∫

Γ+Γ

−

Ω→

Γ

∂−∂−

∂∂−=

=Ω−∇⋅∇−∇−⋅∇

t

t

we

ddn

tgu

n

utg

ddtguutg

τττττ

τττττε

ε

xxξ

xx

xξ

xxξxxxξ

(2.88)

Considerando-se a propriedade (2.71), o 2º termo de (2.87) pode ser calculado como47:

45 Seus limites superiores remetem a condição em que o impulso é infinito 46 Ou soma-se. 47 Valendo-se, adicionalmente, da condição expressa por (2.74).


36

( ) ( ) ( ) ( ) .0,,lim0

0 =Ω−− ∫ ∫−

Ω→

ε

ε ττδδτt

ddtu xxξx (2.89)

Finalmente, tendo em mente as condições (2.69) e (2.72), do 3º termo de (2.87) se tem

que:

( ) ( )[ ] ( ) ( )[ ] ( )

( ) ( ) ( )[ ] ( ).,,,,

,,,,,,lim

0

00

∫

∫

Ω=

=Ω

−=→

Ω−−=

=Ω−−−

xξxξξ

xxxξxxξ

dutgtu

dutgutg t

τ

τετε

ττ

ττττ (2.90)

Substituindo-se as transformações (2.88), (2.89) e (2.90) na expressão original, (2.87),

chega-se a:

( ) ( ) ( ) ( ) ( )

( ) ( ) ( )[ ] ( ) .0,,,,

,,,

,,,

0

0

=Ω+−

−Γ

∂−∂−

∂∂−

∫

∫ ∫

Ω=

Γ+Γ

xξxξξ

xxξ

xx

xξ

dutgtu

ddn

tgu

n

utg

t

we

ττ

τττττη (2.91)

Ou, reorganizando-se os membros:

( ) ( ) ( ) ( ) ( ) ( )

( ) ( )[ ] ( ).,,,

,,,

,,,,

0

0

∫

∫ ∫

Ω=

Γ+Γ

Ω+

+Γ

∂−∂−

∂∂−=

xξxξ

xxξ

xx

xξξ

dutg

ddn

tgu

n

utgtu

t

we

ττ

τττττη (2.92)

O primeiro termo do segundo membro de (2.92) ainda pode ser subdividido entre

fronteiras externas e internas, ou precisamente, limites do reservatório e fontes internas ao

domínio. Isto é:

( ) ( ) ( ) ( ) ( ) ( )

( ) ( ) ( ) ( ) ( )

( ) ( )[ ] ( ).,,,

,,,

,,,

,,,

,,,,

0

0

0

∫

∫ ∫

∫ ∫

Ω=

Γ

Γ

Ω+

+Γ

∂−∂−

∂∂−+

+Γ

∂−∂−

∂∂−=

xξxξ

xxξ

xx

xξ

xxξ

xx

xξξ

dutg

ddn

tgu

n

utg

ddn

tgu

n

utgtu

t

t

e

w

ττ

τττττη

τττττη

(2.93)


37

Este resultado está em consonância com o princípio da superposição, apresentado no

item 2.6. A expressão (2.93) pode ser traduzida matematicamente ou compactada para:

( ) ( ) ( ) ( )0,,,, ξutututuew

++= ξξξ , (2.94)

onde o ( )w

tu ,ξ refere-se à queda de pressão total causada pelas fontes existentes no

domínio; o termo ( )e

tu ,ξ está associado à influência dos limites externos e ( )0,ξu

determina a distribuição inicial dos potenciais.

Portanto,

( ) ( ) ( ) ( ) ( ) ( )∫ ∫Γ

Γ

∂−∂−

∂∂−=

t

w

w

ddn

tgu

n

utgtu

0

,,,

,,,, τττττη x

xξx

xxξξ . (2.95)

Adotando-se, daqui por diante, unicamente fontes de dimensões infinitesimais e de

vazão especificada48, a expressão (2.95) se torna:

( ) ( ) ( ) ( )∫ ∫Γ

Γ∂

∂−=t

w

w

ddn

utgtu

0

,,,, τττη x

xxξξ , (2.96)

Em adição, considerando-se que a admissão de fluido pelas fontes se dê de modo

uniformemente distribuído, a igualdade (2.96) se reduz a:

( ) ( ) ( ) ( )∫ ∫Γ

Γ−∂

∂=t

w

w

ddtgn

utu

0

,,, τττη xxξξ . (2.97)

Ademais, pela introdução da lei de Darcy no formato de (2.29) em (2.97) se tem que

( ) ( ) ( ) ( )∫ ∫Γ

Γ−=t

tw

w

ddtgqc

tu0

,,~1, τττ

φxxξξ , (2.98)

uma vez que

tck φµη 1= . (2.99)

48 Como é hábito na indústria do petróleo.


38

Interpondo-se o resultado (2.98) na expressão inicial (2.93) chega-se a:

( ) ( ) ( ) ( ) ( ) ( )

( ) ( ) ( ) ( ) ( )[ ] ( ).,,,,,~1

,,,

,,,,

00

0

∫∫ ∫

∫ ∫

Ω=

Γ

Γ

Ω+Γ−+

+Γ

∂−∂−

∂∂−=

xξxξxxξ

xxξ

xx

xξξ

dutgddtgqc

ddn

tgu

n

utgtu

t

t

t

w

e

τττττφ

τττττη

(2.100)

Esta equação indica que para se obter a queda de pressão total num ponto qualquer do

domínio, ( )tu ,ξ , eventualmente há necessidade de se conhecer a distribuição inicial das

quedas de pressão ao longo do reservatório. Contudo, na quase totalidade dos casos práticos

observados na indústria petrolífera e na totalidade dos casos investigados neste trabalho,

essa distribuição é assumida como constante e uniforme para 0=t , resultando que a

integral a ser calculada sobre o domínio desaparece. Restam apenas integrais temporais e de

contorno, conforme a expressão49:

( ) ( ) ( ) ( ) ( ) ( )

( ) ( ) ( ) .,,~1

,,,

,,,,

0

0

∫ ∫

∫ ∫

Γ

Γ

Γ−+

+Γ

∂−∂−

∂∂−=

t

t

t

w

e

ddtgqc

ddn

tgu

n

utgtu

τττφ

τττττη

xxξ

xxξ

xx

xξξ

(2.101)

Esta, finalmente, constitui expressão valiosíssima para desenvolvimento deste estudo.

Revela que para que se conheça o potencial em qualquer ponto do domínio, ξ , em um

determinado tempo, t , basta que se conheçam as propriedades do sistema, as vazões a ele

impostas e as quedas de pressão e suas derivadas em cada ponto do contorno –

representados por x50.

49 Por inspeção, pode-se notar que xe ξ são intercambiáveis. Observe-se também que o mesmo resultado

poderia ser obtido com a consideração de fontes já na equação da difusão, a partir de (2.49). 50 Ponto-campo. Para a determinação do potencial numa posição de cálculo, ξ , o ponto-campo deve percorrer

todo o contorno Γ .


39

No entanto, estas duas quantidades – quedas de pressão e suas derivadas no contorno –

não são a priori inteiramente conhecidas51. É necessário, pois, calculá-las. Para tanto, pode-

se posicionar o ponto-fonte, ξ , no próprio contorno.

Para se proceder desta maneira, é necessário, antes, examinar como (2.101) se comporta

no contorno de Ω .

2.8.4. Equação Integral Ponderada de Difusão – Ponto-fonte Pertencente

ao Contorno

Esteja o ponto-fonte posicionado a uma pequena distância, ε , de um vértice de

contorno qualquer definido por um ângulo θ (ver Figura 2.7a).

Figura 2.7 - Fonte nas proximidades de um vértice de contorno.

Considerando-se o Princípio da Superposição no Espaço e o Método das Imagens, já

aludidos anteriormente, o efeito causado por uma fonte aplicada em ξ , distante ε do

vértice, pode ser calculado através do sistema equivalente (ver Figura 2.7b), pela expressão:

( ) ( )∑=

Σ=

N

jj tgtg

1

,,,, xξxξ . (2.102)

51 Para cada porção do contorno do problema, se conhece, inicialmente, apenas uma delas– condições de contorno. Observe-se a Figura 2.6.

a) sistema real b) sistema equivalente


40

O símbolo ( )Σ

tg ,,xξ representa o efeito combinado das N fontes no sistema

equivalente e jξ a posição de cada uma delas.

É importante notar que a magnitude da distância entre cada fonte no sistema

equivalente e o vértice é sistematicamente a mesma (ver Figura 2.7, parâmetro ε )

Fazendo-se com que a fonte e, por conseguinte, suas imagens se aproximem do vértice,

(2.102) passa a ser:

( ) ( )

= ∑

=→Σ

N

jjε

tgtg1

0 ,,lim,, xξxξ . (2.103)

A consideração expressa pelo limite acima implica em que todas as fontes do sistema

equivalente passem a ocupar a mesma posição Vξ . Assim sendo, a expressão (2.103) resulta

em:

( ) ( )tNgtg V ,,,, xξxξ =Σ

, (2.104)

com Vξ sendo a coordenada do vértice e

θπ2=N ,

conforme já discutido em 2.652.

Tendo-se em vista a equação (2.57) e aplicando-se agora a solução fundamental,

correspondente ao desenvolvimento anterior, na equação integral (2.101) passa-se a ter:

( ) ( ) ( ) ( ) ( ) ( ) ( )

( ) ( ) ( )∫ ∫

∫ ∫

Γ

Γ

Γ−+

+Γ

∂−∂−

∂∂−=

t

t

t

w

e

ddtgqc

ddn

tgu

n

utgtuc

0

0

,,~1

,,,

,,,,

τττφ

τττττη

xxξ

xxξ

xx

xξξξ

(2.105)

com

52 Apesar deste resultado ter sido obtido para N inteiro, pode-se provar que esta expressão á válida para qualquer ângulo θ .


41

( ) πθ 2=ξc e Γ∈ξ . (2.106)

Em particular, em contornos ditos suaves, como o apresentado na Figura 2.8

( ) 2/1=ξc , se Γ∈ξ . (2.107)

Figura 2.8 - Fonte pontual nas proximidades de contorno suave. Sistema equivalente.

Observe-se que a equação (2.101) pode ser representada por (2.105) fazendo-se

( ) 1=ξc .

ε

Fonte pontual Fonte Imaginária

Contorno

ε

x ξ Iξ

x

θ


42

2.9. Soluções Fundamentais

Este item é uma espécie de interlúdio neste encadeamento de assuntos. Da

filosofia delineada em 2.8 – particularmente 2.8.3 e 2.8.4 - se poderia partir diretamente

para a obtenção e aplicação do Método dos Elementos de Contorno.

Entretanto, conforme evidenciado principalmente em (2.86) e (2.89), é necessário

discorrer brevemente sobre as já mencionadas soluções fundamentais.

2.9.1. Considerações Preliminares

Conforme já demonstrado anteriormente, para um reservatório de extensão infinita, pelo

Princípio da Superposição no Tempo, a queda de pressão pode ser calculada como:

( ) ( ) ( ) τττdt

t

u

q

qtu

tc

c

−∂

∂= ∫ ,,0

xx .

Por seu turno, da expressão (2.98), ligeiramente modificada, com x em lugar de ξ 53, tal

parâmetro também pode advir de:

( ) ( ) ( ) ( )∫ ∫Γ

Γ−=t

t

ddtgqc

tuw0

,,~1, τττ

φxxξx . (2.108)

Comparando-se ambas as expressões anteriores e tendo-se em conta que as igualdades

devem valer para todo e qualquer tempo, t, resulta que os integrandos de cada qual devem

produzir resultados idênticos.

Considerando-se, ainda, que para as espécies de fontes abordadas neste trabalho, de

dimensões infinitesimais54 vale a igualdade:

( ) ( ) ( )ττ −=Γ−∫Γ

tgdtg w

w

,,,, xξxxξ , (2.109)

então, resulta que:

53 Sem prejuízo algum de seu significado original. As coordenadas espaciais são, de fato, intercambiáveis. 54 Fonte pontual, fonte linear e fonte planar.


43

( ) ( )τφ

τ −=−∂

∂tg

c

qt

t

uw

t

cc ,,~

, xξx . (2.110)

Isto significa que a queda de pressão num sistema sujeito a um impulso unitário pode

ser prontamente obtida pela derivada temporal da solução para o mesmo sistema submetido

à vazão constante.

No espaço transformado de Laplace, a equação (2.110) pode ser escrita como

( ) ( )sgc

qsus

t

cc ,,

~, xξx

φ= , (2.111)

desde que o sistema esteja inicialmente em equilíbrio.

2.9.2. Funções de Green Básicas (ou Soluções Fundamentais)

Com as considerações acima se pode resolver a equação adjunta

( ) ( ) ( ) ( )τδδη −−=∂

∂+∇ tt

tgtg xξ

xx ,

,,2 , (2.112)

partindo-se de

( ) ( )0

,,2 =

∂∂+∇

t

tgtg

xxη , (2.113)

com posterior aplicação de (2.110) ou (2.111).

Fonte Pontual (de Lorde Kelvin)

A primeira solução à (2.112) a ser examinada, e provavelmente a mais importante

delas, é a fonte pontual de Lorde Kelvin55.

Seja um meio poroso ilimitado em todas as direções, inicialmente em equilíbrio –

o que implica em ( ) 00, ==tu x – submetido a uma retirada contínua e pontual de fluido

à taxa q .

Pela lei de Darcy para fluxo esférico isso equivale a escrever que

55 Para maiores detalhes, consulte-se CARSLAW e JAEGER (1959).


44

∂∂−=

Γ

∂∂−= →

=Γ→ ∫ r

ukd

n

ukq

r

w

wµ

πεµ ε

ε

ε2

00 4limlim . (2.114)

O parâmetro ε representa o raio de uma esfera, que deve tender a 0 de forma a que

esta represente um ponto. O sinal negativo advém dos sentidos opostos dos vetores n e r .

Sendo o sistema coordenado esférico e isotrópico, a equação (2.111) deve ser

escrita como:

012

2

2

=∂∂+

∂∂+

∂∂

t

u

r

u

rr

u

η. (2.115)

Ou ainda, com uma transformação de variáveis do tipo56

ruv = , (2.116)

a equação (2.115) se torna

( ) ( )0

,,2

2

=∂

∂+∂

∂t

trv

r

trvη . (2.117)

A lei de Darcy, representada em (2.114), por sua vez, passa a ser expressa por

( )tvk

r

v

r

v

rr

kq ,04

1lim4

22

0 µπ

µπ ε =

−∂∂−= → . (2.118)

Há algumas maneiras de se resolver o problema estabelecido por (2.117) e

(2.118). Uma das mais simples é a solução via transformadas de Laplace.

Tendo-se em conta que ( ) 00, ==trv , a igualdade (2.117), no espaço de Laplace,

pode ser expressa por

( ) ( ) 0,1,

2

2

=+∂

∂srvs

r

srv

η. (2.119)

Já a transformada da equação (2.118), rearranjada, deve ser

56 Veja-se CARSLAW e JAEGER (1959).


45

( )k

q

ssv

πµ

4

1,0 = . (2.120)

onde s é a variável de Laplace.

A solução de equações da espécie de (2.119) é da forma

( ) rsrs

eCeCsrv ηη −−= 21, . (2.121)

Como, pelo fato do domínio em referência ser infinito, aplica-se

( ) 0,lim =∞→ srvr , (2.122)

e então, obrigatoriamente,

01 =C . (2.123)

Pela substituição de (2.120) em (2.121), chega-se a

k

q

sC

πµ

4

12 −= (2.124)

e, em decorrência

( ) rs

ek

q

ssrv η

πµ −

=4

1, . (2.125)

Ou, levando-se em conta a transformação inicial de variáveis, dada por (2.116), a queda

de pressão no espaço transformado vale:

( ) rs

ekr

q

ssru η

πµ −

=4

1, . (2.126)

De acordo com a transformada inversa de (2.126) , com auxílio de (C.5), a expressão

equivalente no espaço real é:

( )

=

t

rerfc

kr

qtru

ηπµ

24, , (2.127)


46

onde erfc é a função erro complementar e é definida conforme

( ) .2 2

∫∞

−=z

u duezerfcπ

(2.128)

Agora, admitindo-se uma retirada instantânea de fluido, de volume unitário, pela

aplicação do enunciado na igualdade (2.110) e a consideração de equilíbrio inicial pode-se

estabelecer que:

( ) ( )t

r

t

etc

qtru η

πηφ4

23

2

4

1~,

−= (2.129)

onde tcq φ~ é a intensidade da fonte e

( ) ( )t

r

et

tg η

πη4

23

2

4

1,,

−=xξ . (2.130)

Esta é a consagrada expressão de Lorde Kelvin, que transformada para o espaço de

Laplace corresponde a (ver Apêndice C):

( ) η

πηsre

rsg −=

41

,,xξ.

(2.131)

Fonte Linear

Imaginando uma seqüência de fontes pontuais dispostas na direção z, para aquele

mesmo espaço infinito definido anteriormente, se obtém prontamente a solução para a fonte

linear57.

Matematicamente, esta consideração pode ser alcançada pela integração da equação

(2.130) em relação z.

( ) ( )∫∞

∞−

−= ζ

πηη de

ttg t

r

423

2

4

1,,xξ , (2.132)

57 Aqui, e como em todo este trabalho, portanto, a fonte linear é ortogonal ao plano cartesiano em que se define o problema bidimensional (x,y).


47

com

( ) ( ) ( )2222 ζψξ −+−+−= zyxr . (2.133)

Retirando-se os termos não relacionados a z da integral, se tem que

( ) ( )( ) ( ) ( )

∫∞

∞−

−−−+−−= dzee

ttg t

z

t

yx

ηζ

ηψξ

πη44

23

222

4

1,,xξ . (2.134)

Operando-se, agora, a mudança de variáveis definida por

( )t

zu

ηζ

4

22 −= , (2.135)

o diferencial dz passa a ser

dutdz η2= (2.136)

e, então, a integral (2.134) passa a produzir (ver ABRAMOVITZ e STEGUN, 1972):

( )

( ) terfctdze t

z

πηπηηζ

22

24

2

=∞−=∫∞

∞−

−−. (2.137)

Portanto, a função de Green para fonte linear é estabelecida por:

( )( ) ( )

t

yx

et

tg ηψξ

πη4

22

41

,,−+−−

=xξ . (2.138)

Ou, transformando-a para o espaço de Laplace, por (C.7),

( ) ( )ηπη

srKsg 021

,, =xξ . (2.139)

Este impulso é de suma importância para o desenvolvimento deste trabalho.

0K é a função de Bessel modificada de 2ª espécie e ordem zero e pode ser calculada

por:


48

( ) ( ) ( )( )∑ ∑

∞

= =

−+

+

−=1 1

12

2

00!

21

2

1ln

r

r

m

r

mr

xxIxxK γ , (2.140)

com

( ) ( )( )∑

∞

=

=0

2

2

0!

4

j

j

j

xxI . (2.141)

O parâmetro γ é a constante de Euler e vale

577215665,0lnlim1

1 ≈

−= ∑=

−∞→ nk

n

knγ . (2.142)

Neste caso, a distância em relação à fonte é definida por:

( ) ( )222 ψξ −+−= yxr (2.143)

Fonte Planar

Seguindo-se o mesmo procedimento anterior, pode-se posicionar uma seqüência de

fontes lineares alinhadas na direção de um dos eixos coordenados e se obter a solução para

a fonte planar58.

Integrando-se, por simples escolha, a expressão (2.138) em y passa-se a ter:

( )( ) ( )

∫∞

∞−

−+−−= dye

ttg t

yx

ηψξ

πη4

22

41

,,xξ . (2.144)

Extraindo-se as constantes da integral chega-se a:

( )( ) ( )

∫∞

∞−

−−−−= dyee

ttg t

y

t

x

ηψ

ηξ

πη44

22

41

,,xξ , (2.145)

cuja integral tem a mesma forma e fornece o mesmo resultado de (2.137). Portanto,

58 O plano assim definido, é de extensão infinita em y e z. É, portanto, ortogonal ao plano (x,y).


49

( )( )

t

x

et

tg ηξ

πη4

2

2

1,,

−−=xξ , (2.146)

com correspondência no espaço de Laplace a (ver transformação C.4):

( ) ηξ

ηsxetg −−=

2

1,,xξ (2.147)

Produtos de Newman

Dependendo da complexidade do problema em análise, em princípio pode ser difícil se

determinar a função de Green apropriada para solução analítica da integral governante

(2.100).

Contudo, NEWMAN (1936) mostrou que para determinadas geometrias e condições

iniciais, tal complexidade pode ser substancialmente diminuída, pois é possível se separar

pura e simplesmente problemas bi ou tridimensionais em problemas unidimensionais

correspondentes. Mais do que apenas isto, o mesmo autor revelou que, em termos

geométricos, a solução obtida por este procedimento corresponde à interseção dos sistemas

unidimensionais relacionados.

No que concerne a este trabalho, pode-se afirmar que aos impulsos instantâneos

referentes ao espaço real podem ser aplicados produtos de Newman.

Uma conseqüência imediata do que NEWMAN (1936) observou é que todo o

desenvolvimento apresentado neste item 2.9.2 pode ser obtido de forma inversa, para

domínios ortotrópicos, a partir da fonte planar (veja-se o Apêndice D).

Pelo exposto, a solução de fonte linear para o caso ortotrópico pode ser obtida pela

consideração de duas fontes planares em direções ortogonais de propriedades distintas (ver

Figura 2.9a). E a fonte pontual, de forma semelhante, pela interseção de uma fonte linear

com uma planar (ver Figura 2.9b).


50

Figura 2.9 – Soluções para domínios ortotrópicos obtidas por produtos de Newman.

Assim, a fonte linear para o caso ortotrópico é

( )( ) ( ) ( ) ( )

t

y

t

x

yx

t

y

y

t

x

x

yxyx et

et

et

tg ηψ

ηξ

ηψ

ηξ

ηηππηπη4444

2222

4

1

2

1.

2

1,,

−−−−

−−−−

==xξ , (2.148)

e a fonte pontual

( )( ) ( ) ( )

( )

( ) ( ) ( )

.2

1

2

1.

4

1,,

444

23

444

222

222

t

z

t

y

t

x

zyx

t

z

z

t

y

t

x

yx

zyx

xyx

et

et

et

tg

ηζ

ηψ

ηξ

ηζ

ηψ

ηξ

ηηηπ

πηηηπ−−−−−−

−−−−−−

=

==xξ

(2.149)

Pode-se notar, nas expressões acima, que caso as constantes de difusão tenham o

mesmo valor em todas as direções (2.148) se torna (2.138) e (2.149) se torna (2.130).

De forma marcante, estes resultados confirmam que a transformação de coordenadas

disposta no Anexo B redunda num sistema isotrópico equivalente.

x

y

a) obtenção da fonte linear. b) obtenção da fonte pontual.

x

y

z z

CAPÍTULO 3: MEC APLICADO AO ESCOAMENTO EM MEIOS POROSOS

51

3. MEC Aplicado ao Escoamento em Meios

Porosos

3.1. Aspectos Gerais

Neste capítulo primeiramente são apresentadas as possibilidades mais evidentes de

solução da equação integral (2.105) via Método dos Elementos de Contorno. Tais

possibilidades estão intimamente relacionadas com a forma de manusear o caráter temporal

do problema.

Na seqüência são apresentadas justificativas para a escolha e são concentrados esforços

na solução via espaço de Laplace.

São explorados elementos de geometria linear e variáveis constantes. Método da

colocação.

3.2. Variação Temporal Aproximada por Série de Taylor

(Diferenças Finitas)

A primeira e aparentemente mais simples possibilidade corresponde à aproximação do

caráter temporal da equação da difusão hidráulica por série de Taylor.

Assim sendo, dando um passo atrás e aplicando-se este conceito à equação da difusão,

(2.78), passa-se a ter:

( ) ( ) ( )0

,,,2 =

∆+

∆∆+−∆+∇

t

tu

t

ttuttu

xxxη . (3.1)

onde ( )tu ,x é uma quantidade conhecida. Trata-se da distribuição de potencial no passo de

tempo anterior ao ora em análise.

Daqui por diante se podem seguir os mesmos passos elencados no item 2.8.3.

Multiplicando-se (3.1) por uma função de ponderação ( )xξ,g , integrando-se no

domínio do problema e separando-se os termos passa-se a ter:


52

( ) ( ) ( ) ( ) ( ) ( )

( ) ( ) ( ) .0,,1

,,1

,, 2

=Ω∆

+

+Ω∆+∆

−Ω∆+∇

∫

∫∫

Ω

ΩΩ

xxxξ

xxxξxxxξ

dtugt

dttugt

dttug

η

η (3.2)

Em adição, com o uso associado da regra de multiplicação de derivadas e do teorema

da divergência, já delineados anteriormente, a 1ª integral de (3.2) pode ser representada por:

( ) ( ) ( ) ( ) ( ) ( )

( ) ( ) ( ) ( ) ( )xxξ

xx

xξ

xxξxξxxξ

∫

∫∫

Γ

ΩΩ

Γ

∂∂∆+−

∂∆+∂+

+Ω∇∆+=Ω∆+∇

dn

gttu

n

ttug

dgttudttug

,,

,,

,,,, 22

(3.3)

e se a função de ponderação ( )xξ,g é solução da equação adjunta ao problema

( ) ( ) ( )xξxξxξ ,,1

,2 δη −=∆

−∇ gt

g , (3.4)

então

( ) ( ) ( ) ( ) ( )ttudgt

gttu ∆+−=Ω

∆−∇∆+∫

Ω

,,1

,, 2 ξxxξxξx η . (3.5)

Com as considerações oriundas de (3.3) e (3.5) aplicadas a equação de partida deste

desenvolvimento (3.1), a equação integral ponderada passa a ser59:

( ) ( ) ( ) ( ) ( ) ( )

( ) ( ) ( ).,,1

,,

,,,

xxxξ

xxξ

xx

xξξ

Ω∆

+

+Γ

∂∂∆+−

∂∆+∂=∆+

∫

∫

Ω

Γ

dtugt

dn

gttu

n

ttugttu η

(3.6)

Pelo exame da equação (3.6) é possível perceber que há necessidade do cálculo de

integral de domínio referente ao passo de tempo anterior, para atualização da distribuição

dos potenciais60. Este tipo de abordagem elimina por completo uma das principais

59 Nesta dedução foram desprezadas, por simplificação, fontes no domínio. 60 Com exceção do primeiro, se a condição de pressão inicial uniforme em todo o domínio for observada.


53

vantagens do MEC em relação a outros métodos numéricos: a falta de necessidade de se

discretizar o domínio do problema em células. Além disto, ocorre que, por haver vinculação

entre os passos de tempo, as discrepâncias em relação à solução analítica são

sucessivamente propagadas. Ou seja, os erros são acumulados.

Com efeito, conforme explicitam BREBBIA, TELLES e WROBEL(1984), para a

obtenção de resultados com boa exatidão a partir desta abordagem, os passos de tempo

devem ser bastante reduzidos.

3.3. Solução no Campo Real com a Solução Fundamental

Dependente do Tempo

Uma segunda possibilidade é a solução numérica de (2.105), reproduzida abaixo.

( ) ( ) ( ) ( ) ( ) ( ) ( )

( ) ( ) ( )∫ ∫

∫ ∫

Γ

Γ

Γ−+

+Γ

∂−∂−

∂∂−=

t

t

t

w

e

ddtgqc

ddn

tgu

n

utgtuc

0

0

,,~1

,,,

,,,,

τττφ

τττττη

xxξ

xxξ

xx

xξξξ

Como puderam observar BREBBIA, TELLES e WROBEL (1984) a solução numérica

da expressão (2.105) só é possível com o estabelecimento de uma marcha de tempo,

evidentemente mais frouxa do que a da abordagem anterior.

Segundo KIKANI e HORNE (1992), esta prática origina a necessidade de se trabalhar

com convolução de matrizes, permanentemente armazenadas em memória. Isto acarreta seu

uso extensivo e, em conseqüência, razoável consumo de tempo. Ademais, uma má escolha

de marcha de tempo pode redundar em propagação/ acumulação de erros.

A despeito disto, ainda segundo os mesmos autores, tal técnica é capaz de produzir

excelentes resultados.

3.4. Solução Via Campo de Laplace

Considerando-se a propriedade de convolução, (2.7), a conversão para o campo de

Laplace e o rearranjo de (2.105) originam:


54

( ) ( ) ( ) ( ) ( ) ( ) ( ) ( )

( ) ( ),,,~1

,,,

,,,,

sgqc

dn

susgd

n

sgsusuc

wt

ee

xξ

xx

xξxxξ

xξξ

τφ

ηη

+

+Γ∂

∂=Γ∂

∂+ ∫∫ΓΓ (3.7)

para 1 única fonte interna, ou para NF fontes internas verticais:

( ) ( ) ( ) ( ) ( ) ( ) ( ) ( )

( ) ( ).,,~1

,,,

,,,,

1∑

∫∫

=

ΓΓ

+

+Γ∂

∂=Γ∂

∂+

NF

mmm

t

sgsqc

dn

susgd

n

sgsusuc

ee

xξ

xx

xξxxξ

xξξ

φ

ηη (3.8)

A solução numérica de (3.8), limitada às características relacionadas no item 1.3, é

efetivamente o escopo deste trabalho.

No espaço de Laplace, devido à eliminação da dependência temporal de u não ocorre

acumulação de erros61.

Ademais, com esquemas de vazão constantes é possível obter a resposta do sistema

exclusivamente no período desejado. Ou seja, não é preciso operar os cálculos desde um

curto tempo inicial.

Mormente, com a adoção de algoritmos simples e robustos, como o de STEHFEST

(1970), é possível se converter, com boa precisão, a solução do campo de Laplace para o

campo real.

KIKANI e HORNE (1993) o classificaram como vantajoso em relação ao de cômputo

por soluções dependentes do tempo.

Mais recentemente, surgiu quem advogasse o contrário (KRYUCHKOV, SANGER e

BARDEN, 2001). Tais autores afirmaram, com certa razão, que no caso de produção de

poços sujeitos a históricos de vazões variáveis, a solução via espaço de Laplace não é a

mais adequada62.

61 Em outras palavras, a determinação da resposta do sistema em um determinado instante independe dos cálculos nos passos de tempo anteriores. 62 Isto, de fato, só foi percebido na fase final de elaboração desta dissertação. O motivo que conduz a esta assertiva é mencionado mais à frente.


55

3.5. Procedimento Numérico

3.5.1. Domínios Simples

Como a grande maioria dos métodos numéricos pressupõe, a solução numérica de uma

equação integral implica na partição do meio continuo sob escrutínio em elementos

discretos e sua – da equação integral - posterior aproximação por um sistema linear.

O Método dos Elementos de Contorno, em particular, resulta da discretização do

contorno do problema, apenas.

Para sua consumação, considerando-se elementos de geometria linear e parâmetros

constantes63, deve-se:

a) Aproximar o contorno do meio físico original por NE elementos, conforme a

relação abaixo64 e a Figura 3.1.

∑=

Γ≅ΓNE

jj

1

(3.9)

Figura 3.1– Representação em elementos de contorno.

Procedendo-se desta maneira, a expressão (3.8) passa a ser escrita como:

63 Segundo a nomenclatura habitual do MEC (cf. com elementos lineares). O que se quer dizer é que a todas as posições elementares correspondem os mesmos valores u e du/dn. O termo constante confere idéia de ausência de variação com o andamento do tempo. Em problemas transitórios, a rigor, isto não ocorre. 64 Ao se dividir o contorno em elementos, o ponto-fonte ao invés de percorrer continuamente o contorno passa a ocupar posições discretas, representadas pelo centro destes elementos.

a) sistema físico original. b) aproximação por MEC.


56

( ) ( ) ( ) ( ) ( )

( ) ( ) ( ) ( ) ( ) ...1,,~1,,,

,,,

,

11

1

NEisgsqcn

sudsg

sudn

sgsuc

NF

miimim

t

jjNE

jjiij

NE

jjjj

jiijiiii

e

j

=+∂

∂

Γ=

=

Γ

∂∂

+

∑∑ ∫

∑ ∫

== Γ

= Γ

ξx

xxξ

xxxξ

ξξ

φη

η

(3.10)

Note-se que os contadores i e j devem variar de 1 a NE. O contador i corresponde à

posição do ponto-fonte (solução fundamental), j se refere ao elemento sob inegração e m às

fontes de domínio (poços).

Observe-se também que ( )su jj ,x e ( )n

su jj

∂∂ ,x

deixam de fazer parte das integrais, pois,

segundo premissa já estabelecida, por definição, não apresentam variação no interior de

cada elemento.

Definindo-se

( ) ( ) ( )

( ) ( )

≠Γ∂

∂

=Γ∂

∂+

=

∫

∫

Γ

Γ

jidn

sg

jidn

sgc

h

j

j

jjiij

jjiij

i

ij

,,,

,,,

xxξ

xxξ

ξ

η

η

; (3.11)

( )suu jjj ,x= ; (3.12)

( ) ( )∫Γ

Γ=j

dsgg jiijij xxξ ,,η ; (3.13)

( )n

sup jj

j ∂∂

=,x

e (3.14)

( ) ( )∑=

=NF

miimim

tj sgsq

cb

1

,~1ξ

φ, (3.15)

a equação (3.10) assume a forma sintética matricial:

bpGuH += . (3.16)


57

Note-se que:

i.os coeficientes ijh e ijg , correspondentes às matrizes H e G , respectivamente,

dependem exclusivamente das propriedades petrofísicas e da geometria do sistema. Podem,

portanto, ser prontamente calculados.

ii.o vetor b , de coeficientes jb , está relacionado à produção das fontes internas ao

domínio (poços). Desta forma, também é conhecido a priori.

iii.entretanto, para cada elemento, inicialmente, se conhece apenas ju ou jp , que são as

condições de contorno inerentes ao problema sob avaliação.

b) Posicionar o ponto-fonte em cada centro de elemento65 e montar o sistema

linear designado pela expressão (3.16), tendo-se como objetivo calcular a incógnita de

contorno restante - ju ou n

u j

∂∂

.

Nesta etapa, pelo fato do ponto-fonte estar alocado sobre o contorno

( ) 2/1=iic ξ . (3.17)

c) Rearranjar o sistema (3.16). Isto significa enviar todos os valores conhecidos

para um dos lados daquela igualdade66. Tal lado da igualdade, de posteriormente é operado

e transformado num vetor de termos independentes (y ) e o sistema original é modificado

para:

yxA = , (3.18)

com x sendo o vetor relacionado aos parâmetros incógnitos no contorno e A a matriz

de coeficientes correspondentes.

65 Denominados de forma corriqueira, nos livros sobre MEC, nós funcionais (cf. nós geométricos). O ponto-fonte é posicionado, então, sobre o nó funcional i de determinado elemento e o ponto-campo percorre elemento a elemento (índice j), integralmente o contorno, dando origem à construção de uma linha i do

sistema linear bpGuH += . Em seguida, o ponto-fonte é deslocado para o elemento seguinte, e o ponto-

campo realiza o percurso novamente. E assim sucessivamente, até que o ponto-fonte tenha sido posicionado em todos os nós funcionais existentes, e como conseqüência, o sistema linear tenha sido completamente definido. 66 Na prática, isto corresponde a mudanças de colunas entre as matrizes H e G .


58

d) Uma vez resolvido o sistema linear e, portanto, com u e n

u

∂∂

já determinados

ao longo de todo o contorno, operar (3.10) – desta feita com i=1 e ( ) 1=iic ξ - e avaliar

quedas de pressão em um (ou mais) ponto(s) interno(s) ao domínio67,68,69.

3.5.2. Domínios Compostos

É ocasião de se expandir os conceitos anteriores para o caso de domínios compostos de

múltiplas regiões homogêneas.

Seja um domínio hipoteticamente dividido em n regiões homogêneas. Por conveniência

de explanação, seja n=3. Por exemplo, seja tal domínio composto aquele apresentado na

Figura 3.5.

Figura 3.2 – Exemplo esquemático de um domínio composto.

Este domínio, para efeito de raciocínio, pode ser repartido em suas regiões

componentes, conforme a Figura 3.3.

Observe-se que cada uma delas possui contornos propriamente ditos e interfaces com as

demais regiões.

67 Importante notar que c) não adiciona erro analítico à solução. 68 Também é possível se determinar as derivadas direcionais das quedas de pressão em pontos internos ao domínio. Desta forma, podem ser percebidas as linhas de corrente do problema. Isto foi explorado por MATSUKAWA e HORNE (1988) e não faz parte do escopo deste trabalho. 69 Aqui há somente multiplicação de vetores. As integrais, por seu turno, precisam ser calculadas.

x

y

1Ω

2Ω

3Ω


59

Figura 3.3 – Porções do domínio composto destacadas.

Suprimindo-se o travessão indicativo do espaço de Laplace, por economia, e adotando-

se os índices sobrescritos das matrizes abaixo para fazerem referência, respectivamente, à

região sob investigação e à outra região a que pertence o elemento de interface se pode

escrever70, para 1Ω que:

[ ] [ ] 1,1

3,1

2,1

1,1

3,12,11,1

3,1

2,1

1,1

3,12,11,1 b

p

p

p

GGG

u

u

u

HHH +

=

. (3.19)

Analogamente, para 2Ω e 3Ω , avaliados de forma isolada, aplicam-se, respectivamente

[ ] [ ] 2,2

3,2

2,2

1,2

3,22,21,2

3,2

2,2

1,2

3,22,21,2 b

p

p

p

GGG

u

u

u

HHH +

=

(3.20)

e

70 Se o elemento faz parte do contorno do domínio composto como um todo, os índices sobrescritos devem se

equivaler. Assim, ii ,H , ii ,G , ii ,u e ii ,p estão relacionados ao contorno do reservatório propriamente

dito. ji ,H , ji ,G , ji ,u e ji ,p , com ji ≠ , estão relacionados às interfaces entre as regiões ( ji Ω∩Ω ). Os

termos relacionados às fontes internas a cada região, ii ,b , evidentemente, não têm relação com outras regiões.

x

y

1Ω

2Ω

3Ω


60

[ ] [ ] 3,3

3,3

2,3

1,3

3,32,31,3

3,3

2,3

1,3

3,32,31,3 b

p

p

p

GGG

u

u

u

HHH +

=

. (3.21)

As matrizes expressas por (3.19), (3.20) e (3.21) apresentam o aspecto revelado na

Figura 3.4. Para construção de cada qual, o ponto-fonte é posicionado em cada elemento

constituinte da região sob consideração e o ponto-campo navega apenas no contorno

correspondente a ji Ω∩Ω . A dimensão de cada uma delas é, portanto, dada pelo total de

elementos pertencentes à região sob análise (m) e pelo total de elementos referentes à

interface (ou, se i=j , contorno) sob avaliação (n).

Figura 3.4 – Aspecto das matrizes.

Mas de acordo com o disposto no item 2.5.3, as interfaces entre as regiões devem

obrigatoriamente obedecer:

a) a compatibilidade entre os potenciais;

ijji ,, uu = (3.22)

e

b) o equilíbrio entre os fluxos, baseado na lei de Darcy.

ji

i

iij

j

j kk ,, ppµµ

−= (3.23)

ou, definindo-se o multiplicador escalar,

ij

jiji

k

kF

µµ

=, , (3.24)

( )ji Ω∩Ω


61

jijiij F ,,, pp −= . (3.25)

Aplicando-se as condições (3.22) e (3.25) aos sistemas expressos por (3.19) a (3.21) e

agrupando-os em um arranjo global71 passa-se a ter:

.3,3

2,2

1,1

3,3

3,2

2,2

3,1

2,1

1,1

3,32,31,21,31,3

3,22,21,21,2

3,12,11,1

3,3

3,2

2,2

3,1

2,1

1,1

3,32,31,3

3,22,21,2

3,12,11,1

+

−−−=

=

b

b

b

p

p

p

p

p

p

GG0G00

0GG0G0

000GGG

u

u

u

u

u

u

HH0H00

0HH0H0

000HHH

FF

F

(3.26)

Considerando-se, ainda, que nas interfaces os potenciais e sua derivadas são a priori

desconhecidos, este sistema global ainda pode ser modificado para:

71 Ou, alternativamente, como fez DORS (2002), se pode partir para um método de solução iterativo, como por exemplo, gradientes bi-conjugados.


62

.3,3

2,2

1,1

3,3

2,2

1,1

3,3

2,2

1,1

3,3

3,2

3,2

2,2

3,1

3,1

2,1

2,1

1,1

3,32,32,32,31,31,31,3

3,23,22,21,21,21,2

3,13,12,12,11,1

+

=

=

−−−

b

b

b

p

p

p

G00

0G0

00G

u

p

u

u

p

u

p

u

u

HGH0GH000

0GHH00GH0

0000GHGHH

FF

F

(3.27)

Este procedimento é geral e, como tal, pode ser aplicado a um sistema de n regiões

homogêneas72.

As regiões, de forma geral, possuem elementos de contorno e de interface.

O caso em que há regiões inteiramente contidas em outras, portanto, é particular.

No sistema apresentado na Figura 3.5 (e repartido na Figura 3.6) , por exemplo, a região

3Ω está contida na 2Ω . Assim sendo, não há interface de 3Ω com qualquer região que não

seja 2Ω . Em conseqüência, o arranjo global dado por (3.28) é compactado para:

.3,3

2,2

1,1

3,3

2,2

1,1

3,3

2,2

1,1

3,3

3,2

3,2

2,2

2,1

2,1

1,1

3,32,32,32,3

3,23,22,21,21,21,2

2,12,11,1

+

=

=

−−

b

b

b

p

p

p

G00

0G0

00G

u

p

u

u

p

u

u

HGH0000

0GHHGH0

0000GHH

F

F

(3.28)

72 A rigor, a (3.27) ainda devem ser interpostas as condições de contorno.


63

Figura 3.5 – Exemplo esquemático de domínio composto com região interna.

Figura 3.6 – Porções destacadas do domínio composto com região interna.

x

y

1Ω 3Ω

2Ω

1Ω

x

y

2Ω

3Ω 2Ω

3Ω


64

3.6. Considerações a Respeito das Implementações

3.6.1. Generalidades

Para efetuação dos estudos propostos foram concebidos dois programas

computacionais. O primeiro, mais simples, para reservatórios homogêneos, denominado

ECDS73; o segundo, para domínios compostos, referido como ECDC74. Evidentemente, a

rigor, o segundo é capaz de resolver problemas da espécie tratada pelo primeiro.

A lógica do código para a simulação de domínios compostos é essencialmente a mesma

de domínios simples. No entanto, há substanciais diferenças no que tange a entrada de

dados75 e a montagem do sistema linear.

A construção de ambos segue o fluxograma apresentado na Figura 3.7.

Após leitura de todos os parâmetros e configurações que cercam o problema em estudo,

incluindo-se aí quantidade de passos de tempo e sua marcha e a quantidade de passos de

Stehfest (NSteh), são percorridos sucessivamente os dois laços destacados no fluxograma.

O laço menor corresponde ao algoritmo de Stehfest. Note-se que para cada passo de

tempo especificado, são montados e resolvidos NSteh sistemas lineares.

O laço maior refere-se aos passos de tempo especificados pelo usuário.

O cálculo de queda de pressão para os pontos internos só é efetuado após completa

definição, para o passo de tempo em curso, das variáveis de contorno do problema.

A seguir são tecidos alguns comentários no que tange os blocos destacados em cinza

(na Figura 3.7), primeiramente relativos a ECDS e posteriormente a ECDC.

73 Acrônimo de Elementos de Contorno, Domínio Simples. 74 Acrônimo de Elementos de Contorno, Domínio Composto. 75 A entrada de dados é extremamente simplificada em ESDS em relação a ECDC.


65

Figura 3.7- Fluxograma de cálculo para MEC aplicado ao Espaço de Laplace com inversão pelo

algoritmo de Stehfest.

Formação do Sistema Linear de Equações (SL)

Modificação do SL (condições de contorno)

Solução do SL

Leitura de Dados

NStehz < Acumulação Solução

1+z

ftt < tt ∆+

1=z

Armazenamento da Solução

1,0 =zt

Impressão de Resultados

Pontos Internos


66

3.6.2. Domínios Simples

Leitura de dados A leitura de dados para domínios simples envolve a definição de 6 blocos de

informação.

A título de ilustração, suponha-se que se queira avaliar o reservatório apresentado na

Figura 3.8.

Figura 3.8 – Domínio simples hipotético. Neste caso, os blocos de entrada para ECDS devem ser (veja-se a Figura 3.9 e a Figura

3.10):

a) coordenadas dos nós geométricos – a conectividade entre eles é

automaticamente assumida como seqüencial, isto é, de nó i para nó i+1. O último nó é

conectado internamente ao primeiro, para fechamento do contorno76;

b) condições de contorno – que vinculam o elemento, representado por seu

centro. São dispostas de forma igualmente seqüencial, de acordo com a ordem assumida

76 O domínio pretendido deve obedecer a convenção de estar a esquerda de quem caminha sobre seu contorno.

1

1

2

2

3

4 3

4

Legenda: Nós Elementos Poços Pts. internos


67

anteriormente. São discriminadas pelos códigos 0 e 1 queda de pressão ou derivada normal

de queda de pressão prescrita, respectivamente;

c) fontes internas – são informadas suas coordenadas e as vazões constantes

correspondentes;

d) pontos internos de cálculo– referenciados por suas coordenadas;

e) propriedades de rocha e fluido – relativas ao domínio sob avaliação;

f) configurações – onde são apontados o número de passos de Stehfest por

passo de tempo (NSteh), o número de pontos de integração numérica por elemento (NPGL),

o tempo inicial de cálculo (t0), o tempo final de cálculo (tf), o número de passos de tempo

(NPT) e sua progressão (cartesiana ou logarítmica).

Figura 3.9 – Parâmetros de entrada 1/2.

Figura 3.10 – Parâmetros de entrada 2/2.

Formação do sistema linear

Com o abastecimento destas informações se pode proceder às sucessivas montagens do

sistema linear, conforme discutido em 3.5.1. Observe-se que são geradas, por passo de

Stehfest, duas matrizes cheias e não simétricas.

Considerando-se poços representados por fontes lineares, cujas soluções fundamentais

no campo de Laplace sejam do tipo (2.139), os coeficientes de H apontados em (3.11)

podem ser calculados como:


68

( ) ( ) ( )dr

r

sgd

n

sgh

jj

jijji

jji

ij ∫∫ΓΓ

⋅∇∂

∂=Γ

∂∂

= nrxξ

xxξ ,,,,

ηη . (3.29)

Figura 3.11 - Integração Numérica – quando x e ξ não pertencem ao mesmo elemento.

Admitindo-se a avaliação da integral baseada no procedimento de Gauss-Legendre, a

equação anterior pode ser modificada para (veja-se a Figura 3.11):

( )

jjik

NPGL

k

kikij Jr

r

sgwh n

xξ⋅∇

∂∂

= ∑=1

,,η , (3.30)

onde jJ é o determinante do Jacobiano da transformação de coordenadas do sistema

global ( )yx, para o sistema local( )0,υ e, neste caso, vale 2/jL .

A derivada embutida em (3.30) deve ser expressa por:

( ) ( )ηπη

srKs

r

sgik

ki12/3

2/1

2

,, −=∂

∂ xξ. (3.31)

Por seu turno, a magnitude do vetor que mede a distância entre o ponto localizado no

elemento sob integração e a posição da fonte é definida por

( ) ( )22ikikik yxr ψξ −+−= (3.32)

e seu gradiente vale:

αυ

αυ

α

α

cos2

2

)(cos

)(

12

12

jmk

jmk

j

j

Lyy

senL

xx

L

yy

L

xxsen

+=

−=

−=

−−=

α

),( mm yx n Posição da fonte

Elemento sob integração

* ),( ii ψξ

),( 11 yx

),( 22 yx

ψ,y

ikr

ξ,x

υ

12 =υ

11 −=υ

jΓ


69

( ) ( )[ ]ji ikikik

ik yxr

r ψξ −+−=∇ 1. (3.33)

Já o vetor normal ao elemento considerado pode ser representado por seus cossenos

diretores, isto é

jin αα sen+= cos , (3.34)

ondeα é o ângulo que o vetor normal ao elemento faz com o eixo x (ver Figura 3.11).

Os valores numéricos de tais cossenos diretores podem ser obtidos pela semelhança entre

os triângulos destacados na figura já apontada.

Portanto, aplicando-se (3.31) a (3.34) em (3.30) se chega a:

( ) ( ) ( )[ ]αψαξηηπ

senyxr

srKwsL

h ikikik

erN

kikk

jij −+−−= ∑

=

cos1

4

int

11 . (3.35)

Mas se o elemento sob integração é o mesmo em que se encontra o ponto-fonte então77:

( ) ( )2/1

,,22/1

,,22/1

2

0

2

0

=⋅∇∂

∂+=∂

∂+= ∫∫ drr

sgdr

n

sgh

L

iiiii

L

iiii nr

xξxξ ηη . (3.36)

Analogamente, os coeficientes de (3.13) devem ser calculados por:

( ) ( )jjiij dsggj

xxξ Γ= ∫Γ

,,η , (3.37)

ou, substuindo-se a função de Green pertinente,

( ) ( )∫Γ

Γ=j

jijij dsrKg xηπ 02

1, (3.38)

ou ainda, por Gauss-Legendre

( )∑=

=NPGL

kikk

jij srKw

Lg

104

ηπ

. (3.39)

Se, entretanto, o ponto-fonte está situado no próprio elemento sob integração, devido à

simetria, a expressão (3.37) pode ser avaliada analiticamente por: 77 Note-se que, neste caso nr ⊥ .


70

( )drsggL

ii ∫=2

0

,,2 xξη . (3.40)

Ou, aplicando-se a solução fundamental:

( )∫=2

0

0

1 L

ii drsrKg ηπ

. (3.41)

Operando-se sobre (3.41) a transformação abaixo relacionada

drsdu

sru

η

η

=

=, (3.42)

chega-se finalmente a

( )∫=η

ηπ

sL

ii duuKs

g2

0

0

1. (3.43)

Integrais da espécie de (3.43) podem ser calculadas (ver ABRAMOVITZ e STEGUN,

1972) através da expressão:

( ) ( )( ) ( )

( )( ) ( )

( )( ) ( )∑∑

∑∑∫

=

∞

=

∞

=

∞

=

++

++

++

+=

k

nk

k

k

k

k

kx

nkk

x

kk

xx

kk

xx

xduuK

102

2

022

2

02

2

0

0

1

12!

2

12!

2

12!

2

2ln γ

. (3.44)

Conforme Figura 3.12, pode-se perceber que K0 é monótona e decrescente e, em

conseqüência, sua integral é limitada ao valor abaixo.

( )20

0

π=∫∞

duuK (3.45)


71

Figura 3.12 - Integral da Função de Bessel Modificada de 2ª Espécie e Ordem 0.

Com relação ao termo correspondente às fontes internas, (3.15), pode-se estabelecer

que:

( )sgcsh

qb w

ti ,,xξ

φ= (3.46)

para q constante, tendo-se em conta que para fontes lineares hqq /~ = .

Substituindo-se a expressão correspondente a ( )sg w,,xξ se tem que:

( )ηπµ

srKskh

qb iwi 0

12

= , (3.47)

com

( ) ( )22wiwiiwr ψψξξ −+−= . (3.48)


72

Considerando agora, hipoteticamente, um reservatório “infinito” para o tempo total de

cálculo78, com (3.47) em (3.10) passa-se a ter:

( )ηπµ

srKskh

qu iwi 0

12

=v. (3.49)

A expressão (3.49) corresponde à solução analítica para produção à taxa constante de

fonte linear em reservatório “infinito”79. Como conseqüência imediata e direta, se pode

concluir que até que haja influência de algum limite ou interface de reservatório o caráter

analítico da solução é preservado.

Modificação do sistema de equações

A modificação do sistema de equações corresponde a transformação do sistema

bpGuH += em yxA = pela consideração das condições de contorno do problema em

objeto.

Tal etapa consiste simplesmente em mudanças de colunas entre as matrizes H e G ,

com todos os valores conhecidos sendo posicionados do lado direito da primeira igualdade,

posterior operação deste lado da igualdade por produto matriz-vetor e subseqüente soma

com o vetor corresponde às fontes existentes no domínio, conforme procedimento descrito

em 3.5.1.

Solução do sistema de equações Uma vez obtido o sistema linear yxA = procede-se à sua solução pelo método direto de

eliminação de Gauss, com pivoteamento parcial80.

Calculo nos pontos internos O fluxograma disposto na Figura 3.7 mostra que o cálculo em pontos internos é

efetuado após a completa definição dos valores das incógnitas no campo real. Em outras

78 Em outras palavras, assumindo que neste período não há qualquer influência de limites ou barreiras ao fluxo. 79 A adimensionalização de variáveis utilizada em KIKANI e HORNE (1992,1993) não é compatível com esta observação. 80 Métodos diretos são mais eficientes no trato com problemas relativamente pequenos.


73

palavras, somente após o término e conseqüente acumulação dos valores numéricos obtidos

através dos passos de Stehfest.

Procedimento diferente, com o cálculo para os pontos internos efetuado ainda no laço

menor, foi tentado sem sucesso.


74

3.6.3. Domínios Compostos

A manipulação de domínios compostos segue, em linhas gerais, as mesmas regras

adotadas para domínios simples. A mais significativa modificação corresponde à

necessidade do trato de mais de 1 sistema bpGuH += , referente a cada região.

A leitura de dados sofre impacto significativo.

Leitura de dados

À guisa de ilustração, seja o reservatório representado na Figura 3.13. Como se pode

observar, diversamente do que ocorre para domínios simples, passa a haver necessidade de

se informar a conectividade entre os nós.

Assim, inicialmente se deve definir a numeração das regiões existentes. Isto posto, é

necessário que se estabeleça um critério para a incidência dos elementos de interface. Para

confecção do programa computacional ECDC foi adotado o critério de conectividade

informada como se elemento pertencesse à região de menor valor.

Evidentemente, é preciso que se esclareçam as duas regiões as quais pertence o

elemento (veja-se a Figura 3.14).

Figura 3.13 – Exemplo para Domínio Composto: conectividade.

1 2

3 4 5

6

7

8

9

10

11

1

2

3 4

5

6

7

8

9

10

11

12 13 1Ω

2Ω

3Ω

x

y


75

Assim, por exemplo, a conectividade do elemento 11, que pertence às regiões 1 e 2

deve se dar do nó 10 para o nó 1.

Em ECDC, a quantidade de blocos de entrada passa de 6, caso de ECDS, para 7.

Com efeito, o campo de coordenadas dos nós geométricos deve ser complementado pelo de

conectividade dos elementos.

Figura 3.14 – Exemplo: tabelas de entrada 1/3. Os demais dados de entrada em muito se assemelham aos referidos anteriormente.

Para as fontes internas e pontos internos de cálculo passa a ser necessário que se

informe a região a que pertencem (conforme Figura 3.15).

Figura 3.15 – Exemplo: dados de entrada 2/3. Evidentemente, também há necessidade de se especificar as propriedades por região, de

forma seqüencial, conforme Figura 3.16.


76

Figura 3.16 – Exemplo: dados de entrada 3/3.

Formação do sistema linear

Com toda a massa de dados lida, procede-se à montagem do sistema linear, de

acordo com diversas ocasiões demandadas.

Figura 3.17 – Sentido positivo de percurso para as regiões consideradas de forma isolada. As regiões são isoladas e as correspondentes matrizes são organizadas sob a forma

de matrizes globais, cuja topologia, para o exemplo em tela, corresponde a Figura 3.18.

1

1Ω 2

3

8

9

10

11

8

9

2Ω

3Ω

11

10 12 13

6

7

5

4

12 13


77

Esta é justamente a etapa de maior dificuldade na construção de códigos para domínios

compostos. Com base em (3.27) se deve alocar os coeficientes num arranjo global esparso.

Para esta função, lançou-se mão da construção de vetores apontadores81, conforme os

designados abaixo:

- vetor apontador do nó funcional em que se encontra o ponto-fonte;

- vetor apontador da região sob avaliação;

- vetor apontador do elemento sob integração para a matriz H modificada;

- vetor apontador do elemento sob integração para a matriz G modificada.

Figura 3.18 – Aspecto das matrizes globais e do vetor global para o exemplo considerado.

No terceiro vetor acima aludido as repetições seqüenciais indicam colunas

representativas de coeficientes g e h de interface.

O cálculo dos coeficientes, representados por “x” na Figura 3.18, é baseado na mesma

formulação disposta em 3.6.2.

81 Facilmente obtidos dos registros reg1 e reg2, informados pelo usuário. Para maiores detalhes consulte-se o Apêndice F.


78

Note-se que, neste caso sintético e ilustrativo, a região 3 não possui fontes.

O caso em que há regiões internas a outras é particular e, deste modo, também pode ser

manipulado pelo programa computacional implementado.

Modificação e Solução do sistema de equações

Estas duas etapas são idênticas as executadas pelo programa ECDS.

Entretanto, ao final do processo, para ECDC, há necessidade de se separar as variáveis

por região, incluindo-se aí a troca de sinal nas interfaces, quando necessária.

Calculo nos pontos internos

Isto posto, se pode proceder ao cálculo para os pontos internos solicitados.

A queda de pressão em um ponto interno à determinada região é função exclusiva do

contorno dela própria.

CAPÍTULO 4: APLICAÇÕES NUMÉRICAS

79

4.Aplicações Numéricas

4.1. Considerações Iniciais

Com o objetivo de avaliar os programas computacionais escritos com vistas à

consolidação do presente estudo, na seqüência são apresentadas algumas comparações com

soluções analíticas já consagradas na literatura especializada82.

4.1.1. Metodologia

Idealmente, tais comparações devem abarcar toda a faixa de valores de parâmetros

práticos e devem contemplar toda sorte de geometrias imagináveis. Na prática, um estudo

com toda essa abrangência e em curto período não é factível.

Nesta dissertação, a metodologia adotada é a de iniciar as verificações a partir de

modelos simples e incorporar gradativamente complexidades, acumulando aprendizado

exemplo após exemplo83. Desta forma, é possível estabelecer diretrizes práticas de uso dos

códigos, ainda que não comprovadas à exaustão.

O limitante deste estudo são, na prática, as soluções analíticas existentes. A rigor, tais

soluções em grande parte estão sujeitas a geometrias e a condições de contorno

simplificadas - basicamente 0ˆ =u ou 0ˆ

=∂∂n

u no contorno.

Fundamentalmente, o termo de comparação é solução no poço84. Exceto onde indicado,

as soluções analíticas são as extraídas de um programa comercial largamente utilizado na

indústria de petróleo, para fins de interpretação de testes em poços, denominado

PANSYSTEM85. Neste programa, de forma geral, as fronteiras são contabilizadas pelo

método das imagens.

82 Em uma palavra, validação dos códigos gerados. 83 Método científico de indução. Examinar casos particulares sucessivamente no intuito de se estabelecer generalizações. 84 Os potenciais, ou quedas de pressão, são medidos exclusivamente nos poços de petróleo, com auxílio de registradores eletrônicos. A rigor, não é possível mensurá-los diretamente em outras posições do sistema. 85 Versão 3.5. Para maiores informações, consultar o sítio: http://www.ep-solutions.com/solutions/eps/PanSystem.htm.


80

Recursos da Teoria de Interpretação de Testes em Poços de Petróleo

Quando um sistema de petróleo é posto em produção, é gerado um transitório de

pressões que viaja pelo reservatório, do poço para adiante, sendo alterado de acordo com os

obstáculos e modificações nas características físicas do meio poroso atravessado.

O consagrado gráfico de diagnóstico é um recurso dos mais largamente utilizados com

fins de interpretação de testes em poços de petróleo. Este aparato corresponde à

representação em eixos logarítmicos das curvas de queda de pressão e derivada de queda de

pressão.

A de queda de pressão é prontamente obtida por meios analíticos e numéricos86. Já a de

derivada da queda de pressão, definida por BOURDET et al (1983), para um único período

de fluxo é calculada através da expressão:

td

duu

ln'= . (4.1)

Esta curva possibilita a identificação dos regimes de fluxo característicos perpassados

pelo transitório de pressões, como, por exemplo os apresentados na Tabela 4-1. Além do

mais, em relação à curva de quedas de pressão, a derivada apresenta variações

amplificadas.

Tabela 4-1 – Regimes de fluxo característicos, observados através da curva de derivadas de queda de pressão, utilizados neste trabalho.

REGIME ASPECTO

APRESENTADO

OBSERVAÇÃO

Radial ( )tu ln∝ Patamar Reservatório infinito

Barreiras em ângulo

Modelo radial composto

Linear ( )tu ∝ Reta de Inclinação ½ Poço fraturado

Barreiras paralelas.

Cartesiano ( )tu ∝ Reta de inclinação 1 Fronteiras seladas

Constante( )cteu = Derivada cai abruptamente Fronteira realimentada

86 Como o MEC.


81

Neste trabalho a validação dos códigos computacionais gerados está quase que

integralmente assentada neste tipo de representação gráfica.

Uma primeira verificação, de cunho qualitativo, é exercida justamente sobre a curva de

derivada de queda de pressão. Consiste em observar sua forma. A segunda, quantitativa,

corresponde ao cálculo da discrepância relativa entre a resposta analítica e a advinda do

MEC.

Nos gráficos apresentados a seguir, as soluções analíticas são sistematicamente

representadas por linhas cheias e as soluções advindas da aplicação do MEC por círculos.

Os erros relativos87 são dispostos no eixo secundário das ordenadas e são determinados de

acordo com a seguinte expressão:

ANALÍTICAMECR uu−= 1ε (4.2)

São utilizadas unidades brasileiras de campo. As conversões para o Sistema

Internacional são achadas no Apêndice E.

4.1.2. Base de Comparação

O julgamento da eficácia das implementações baseadas no MEC, nas especulações a

seguir, é inteiramente baseado no binômio custo-benefício.

O termo custo está relacionado ao tempo total de conclusão dos cálculos. Já a palavra

benefício, à exatidão da resposta. Ambos são influenciados sobremaneira:

• pela quantidade de elementos de contorno necessários à representação adequada

do problema físico-matemático;

• pela quantidade de passos de Stehfest, adotada para transformação inversa da

resposta ao campo real;

• pela quantidade imposta de pontos de integração numérica de Gauss-Legendre.

O custo é influenciado, ainda, pela quantidade de passos de tempo desejada.

87 Ou discrepâncias relativas. Daqui por diante estes termos são usados de forma indistinta.


82

De antemão, se pode assegurar que as principais fontes de inexatidão do processo são as

seguintes:

• a discretização do contorno em elementos de parâmetros constantes;

• a integração numérica sobre os elementos de contorno de parâmetros constantes;

• os erros acarretados pelo processo numérico de eliminação de Gauss, utilizado

para determinação das incógnitas do sistema linear;

• a aproximação das funções de Bessel por polinômios racionais;

• os erros de truncamento associados aos cálculos. Com efeito, há muitas

operações com números grandes.

4.2. Discussão Acerca do Consumo de Tempo

Os programas computacionais implementados podem ser resumidos em 5 grandes

blocos: leitura de dados, montagem do sistema linear associado, resolução do sistema

linear, cálculo em pontos internos e impressão de saída.

Dentre estes, montagem e resolução do sistema linear são os computacionalmente mais

intensos.

Para ilustrar esta discussão acerca do consumo de tempo, considere-se, por

simplificação, um reservatório circular com limites externos selados88, com um único poço

ao centro, de acordo com os parâmetros apresentados na Tabela 4-2, abaixo.

Tabela 4-2- Parâmetros adotados para os estudos iniciais. Parâmetro Valor

k [mD] 1000

φ 15%

µ [cP] 1,5

ct [cm²/kgf] 1,5.10-4

h [m] 25

q [m³/d] 350

re [m] 10.000

88 O que equivale a dizer 0ˆ

=∂∂n

uem todo o contorno.


83

Avaliação dos Tempos de Montagem e Resolução do Sistema Linear

Mantendo-se fixos a quantidade de passos de Stehfest (6), a quantidade de pontos de

integração de Gauss-Legendre (8), o número de passos de tempo (100) e o número de

pontos internos de cálculo (1, na parede do próprio poço) procede-se à construção do

gráfico abaixo.

Figura 4.1 - Avaliação do consumo de tempo.

A Figura 4.1 mostra claramente que nas circunstâncias descritas, com até 100 elementos

de cálculo o tempo gasto na montagem do Sistema Linear é superior a 85% do tempo total

por rodada do programa.

Com algo em torno de 1000 elementos, os tempos de montagem e resolução do sistema

linear passam a se equivaler.

Ademais, na mesma figura é possível perceber, através das inclinações das retas

demarcadas, que a montagem do sistema linear possui dependência quadrática com a


84

quantidade de elementos e a resolução deste mesmo sistema, dependência

aproximadamente cúbica89.

Avaliação da Natureza e do Tempo Gasto pela Integração Numérica

Para a consecução dos cálculos do MEC é necessário se calcular integrais do tipo

( )∫Γ

Γj

dsrK ij η0 e ( )∫Γ

Γ⋅∇j

drsrK ijij nη1 . Para tanto, nos casos em que o ponto-fonte se

situe numa posição não pertencente ao elemento sob integração, é aplicada a Quadratura de

Gauss-Legendre90.

A função de Bessel de 2ª espécie modificada 0K e sua derivada, 1K , são as dispostas na

Figura 4.2. Como se pode observar, são monótonas e assumem valores significativos para

argumentos pequenos. A função K1 é mais aguda.

Figura 4.2 - Comportamento de K0 e K1.

89 Para domínios compostos de múltiplas regiões NE faz o papel de elementos não nulos da matriz H(ou G). 90 Também conhecida como Quadratura Gaussiana. O método fornece resposta exata para integração de

polinômios de ordem 2n-1. N é o número de pontos necessários para cômputo de ( ) ( )∑∫=−

≈n

iii fwdxxf

1

1

1

ξ .

Os valores de iw e iξ provêm dos polinômios de Legendre e são freqüentemente encontrados sob a forma

tabelada.


85

O argumento das funções de Bessel para o problema em análise é ηsrij .

A constante η depende exclusivamente da definição das características físicas do

problema e tipicamente assume valores entre 10-2 a 102, em unidades internacionais.

O valor ijr depende da geometria do problema – comprimento dos elementos e distância

das posições assumidas pela fonte em relação ao elemento sob integração - e a variável de

Laplace, s, diminui à medida que o tempo aumenta, posto que Lt=1/s².

Fixando-se o comprimento de um hipotético elemento de integração em 100=jL e

tomando-se, por exemplo, o conjunto 10,0=ηs se pode avaliar a influência da posição

relativa do ponto-fonte em relação ao centro do elemento. Na Figura 4.3 se pode observar

que à medida que a fonte se aproxima do elemento, a função ( )ηsrK ij0 assume maiores

valores de pico e revela formas mais acentuadas, o que acarreta maior dificuldade de

integração.

Figura 4.3 - K0(x) com aproximação da fonte pontual na direção ortogonal ao elemento (unidades SI).

Por outra, adotando-se uma posição ortogonal ao elemento fixa e transladando-se o

ponto-fonte na mesma direção do elemento, conforme Figura 4.4, se vê que a função


86

preserva seu valor de pico até que a fonte se afaste da região delimitada pelos vetores

ortonormais aos extremos do elemento sob integração. Com a diminuição do valor de pico

e a função passando a ser monótona, passa a haver menor dificuldade de integração.

Figura 4.4 - K0(x) para variação da posição da fonte na direção do elemento (unidades: SI).

Claro fica, portanto, que a distância ao centro do elemento deva ser um critério para a

escolha da quantidade de pontos de integração. E, intuitivamente, se espera que o

comprimento do elemento também influencie no processo.

Em resumo quanto menor j

ij

L

sr η, maior tende – apenas tende – a ser a exatidão da

integração numérica para um mesmo NPGL.

Maiores quantidade de números de pontos integração é necessária à medida que:

a) a fonte se aproxima do centro do elemento em objeto. Particularmente, se a

fonte estiver posicionada no espaço compreendido entre as projeções

ortonormais às suas extremidades;

b) o tempo avança.


87

Como o avanço do tempo é inexpugnável, a tarefa de obtenção de um critério universal

de seleção automática – portanto, interna aos programas; sem intervenção do usuário - do

número ideal de pontos de integração numérica (NPGL) é dificultada.

Quanto ao tempo gasto na operação de tais integrações numéricas se pode dizer que:

com base ainda no reservatório definido anteriormente, sujeito às mesmas condições já

especificadas, foi concebida a Tabela 4-3. Por esta tabela, se pode perceber que o tempo

gasto na montagem do sistema linear é função quadrática do número de elementos

utilizados. Em outras palavras, se a 2ª, 3ª e 4ª colunas são normalizadas em relação ao

quadrado do número de elementos, para cada quantidade de pontos de Gauss-Legendre

considerada na primeira coluna, as linhas terão valores muito próximos entre si.

Tabela 4-3 – Tempo consumido em função do número de pontos de Gauss-Legendre Considerado apenas 1 ponto interno, representando o poço.

Tempo de Montagem do Sistema Linear[s] Número de pontos

de Gauss-Legendre NE=6 NE=12 NE=24

4 5,02 18,77 75,41

8 5,53 19,72 85,31

16 6,25 24,8 101,83

32 8,22 32,91 136,18

64 12,37 49,37 206,84

100 16,39 67,61 281,05

Mais importante ainda é constatar que o acréscimo de número de pontos de Gauss-

Legendre corresponde a um aumento do consumo de tempo inferior a relação quadrática

entre os casos considerados.

Por exemplo, a um aumento do número de pontos de Gauss-Legendre de 4 para 100

(25 vezes) corresponde um acréscimo no tempo de montagem do sistema linear inferior a 5.

Avaliação da quantidade de vezes em que se monta o sistema linear

Conforme exposto acima, para problemas relativamente pequenos, com até poucas

centenas de elementos, o tempo de cálculo é majoritariamente consumido em montagens do

sistema linear. Estas, por sua vez, são função do número de elementos a que se lança mão


88

para a representação do problema físico, NE, da quantidade de passos de tempo

especificada (NPT) e da quantidade de passos de Stehfest por passo de tempo (NSteh).

Considerando que NE seja função da exatidão esperada da resposta e que a quantidade

de passos de tempo esteja de acordo com o gosto do usuário dos códigos (NPT), deve-se

buscar limitar a quantidade de passos de Stehfest (NSteh).

Avaliação da quantidade de passos de Stehfest

Para fins de ilustração, seja o reservatório circular com poço centralizado, apresentado

na Figura 4.5 (em vermelho), representado para efeitos de cálculo pelo hexágono em azul.

Portanto, NE=6.

Considere-se que todas as suas fronteiras sejam seladas91.

Fixando-se o número de pontos de integração numérica em 8, se pode avaliar

unicamente a influência da quantidade de passos de Stehfest (NSteh) na exatidão da

resposta.

A Figura 4.6, a Figura 4.7 e a Figura 4.8 mostram os resultados obtidos,

respectivamente, para NSteh=2, 6 e 10.

Figura 4.5 – Representação do reservatório circular com 6 elementos.

91 Portanto, com .0ˆ

=∂∂n

u


89

Figura 4.6 – Reservatório circular com poço centralizado (NE=6, Ninter=8, NSteh = 2).

Figura 4.7 – Reservatório circular com poço centralizado (NE=6, Ninter=8, NSteh=6).


90

RESERVATÓRIO CIRCULAR COM LIMITES SELADOSPOÇO CENTRALIZADO

0.1

1

10

100

0.01 0.1 1 10 100 1000 10000 100000

∆∆∆∆t [h]

u,u'

[kgf

/cm

²]

-7%

-6%

-5%

-4%

-3%

-2%

-1%

0%

εε εε R

u analítica u' analíticau MEC u' MECeR

Figura 4.8 - Reservatório circular com poço centralizado (NE=6, Ninter=8, NSteh=10).

A curva analítica em vermelho, a de derivada de pressões, revela que até cerca de 300h

as fronteiras não impõem contribuições significativas à solução. Até este instante, portanto,

o escoamento é desenvolvido de maneira radial e infinita. A partir de então, há uma curta

transição até que o transitório de pressões encontre concomitantemente os limites selados

do sistema, fazendo com que se desenvolva o regime pseudo-permanente92 – inclinação

unitária de derivada de pressões93.

Há, portanto, dois regimes de fluxo característicos a serem observados pelo MEC.

Com gritante evidência, é possível perceber que NSteh = 2 não é capaz de representar o

segundo deles (Figura 4.6). A rigor, o caráter da resposta naquele trecho é preservado mas

erros superiores a 100% são obtidos ao final do período.

Com NSteh=6, mantidas as demais condições, há substancial redução da discrepância

para a casa dos 6% na transição entre os dois regimes aludidos, com posterior redução do

erro à medida que o tempo avança (Figura 4.7).

92 Regime em que a queda de pressão passa a ter igual intensidade em todas as posições do domínio. 93 Consulte-se a Tabela 4.1.


91

Com NSteh=10 não há melhoria significativa da exatidão da resposta, em comparação

com os resultados advindos da aplicação de NSteh=6, que possam justificar preferência

por este valor. Com efeito, as discrepâncias relativas praticamente não se alteram.

Este simples exemplo é apenas ilustrativo e serve apenas para revelar, como número

ideal para os cálculos via MEC, a escolha por NSteh = 6. Na realidade, foram examinados

diversos casos e rigorosamente todos apontaram para a mesma preferência94.

Outros valores, em princípio possíveis, curiosamente redundam em resultados espúrios

(NSteh=4, NSteh=8, NSteh=12).

Valores maiores podem induzir erros relacionados à precisão de máquina. Este tema é

explorado em item específico, mais adiante.

4.3. Avaliação da Exatidão dos Códigos Gerados

A seguir são apresentados exemplos inspirados na metodologia disposta em 4.1.1.

É avaliado inicialmente um modelo bem simples e posteriormente são interpostas

complexidades, com o intuito de explorar as funcionalidades e os limites dos códigos

concebidos. O aprendizado é incorporado nas simulações subseqüentes.

Os poços são considerados como produzidos a taxas constantes.

O julgamento da exatidão corresponde ao confronto das soluções geradas pelo MEC

com as analíticas. É necessário, pois, que: 1) a curva de derivada de queda de pressão

represente o caráter esperado da resposta e, 2) se sim, se avalie a discrepância relativa

associada, calculada com base nas curvas de queda de pressão.

Exemplo 1: Reservatório Quadrado – Fronteiras Seladas – Poço Centralizado

Seja, portanto, inicialmente um reservatório quadrado com fronteiras seladas,

conforme o quadro menor da Figura 4.9.

Sejam as principais propriedades do sistema, as dispostas na Tabela 4-4.

Para este reservatório são, naturalmente, esperados dois regimes de fluxo: nos

tempos iniciais, o radial infinito, diagnosticado pela curva de derivadas como um patamar e

94 O trabalho desenvolvido por KINKANI e HORNE (1992) apresentou a mesma conclusão.


92

nos tempos finais, o pseudo-permanente, representado por uma reta de inclinação unitária95.

Entre eles há uma transição.

Tabela 4-4 – Exemplos 1, 2, 3 e 4: parâmetros considerados.

PARÂMETRO VALOR

k [mD] 800,0000

fi 0,2000

mi [cP] 1,0000

ct [cm²/kgf] 0,00010

h [m] 50,0000

q [m³/d] 600 rw [m] 0,10

Adotando-se 4 elementos de contorno e 8 pontos de integração numérica é obtido o

resultado apresentado na Figura 4.9.

Figura 4.9 – Exemplo 1– NE= 4, NPGL=8.

O caráter da resposta esperada é inteiramente reproduzido, mesmo com este

reduzido número de elementos.

95 Veja-se a Tabela 4.1.


93

As discrepâncias são insignificantes no primeiro regime atingido, devido à

característica analítica da solução até aquele momento. Quando a influência das fronteiras

passa a produzir efeito, ocorre um período de transição, ligeiramente diferente da reposta

analítica, com aumento do módulo da discrepância relativa até cerca de 8% e, uma vez

atingido este valor máximo, posterior redução.

Imaginando-se que a discrepância obtida anteriormente pudesse estar associada à

imprecisão na avaliação das integrais, componentes das matrizes H e G , procedeu-se ao

aumento dos pontos de integração numérica de Gauss-Legendre de 8 para 3296.

Figura 4.10 – Exemplo 1– NE= 4, NPGL=32.

Como se pode observar, através da Figura 4.10, tal modificação não produziu o efeito

desejado. O erro máximo se manteve-se na mesma faixa.

Restou, então, aumentar a quantidade de elementos.

96 Se esta hipótese viesse a se confirmar, seria de bom grado, posto que o aumento deste parâmetro acarreta menor consumo de tempo computacional que o aumento no número de elementos de contorno.


94

Representando-se, agora, o mesmo problema físico através de 16 elementos de

comprimento uniforme, o erro máximo de transição passa a 6% (ver Figura 4.11). Esta

quantidade de elementos corresponde a um comprimento, por elemento, de 1000m.

Seguindo-se na mesma direção e ampliando-se o número de elementos para 64, de

comprimento uniforme, não são obtidos ganhos substanciais (ver Figura 4.12).

Figura 4.11 - Exemplo 1– NE= 16, NPGL=8.


95


Por outra, tendo-se em mente o discutido anteriormente a respeito da integração das

funções de Bessel, se poderia imaginar que uma distribuição não uniforme dos elementos

como a apresentada na Figura 4.13 pudesse causar diminuição no erro relativo.

No entanto, como se pode observar, isto não se verificou.

Assim, imagina-se que este erro típico de transição entre períodos de fluxo

característicos deva ser intrínseco ao próprio método, pelo menos da maneira que tenha

sido implementado, com elementos de parâmetros constantes.


96

Figura 4.13 - Exemplo 1– NE= 64 não uniformes, NPGL=8.

Com efeito, em outras posições do reservatório, o método comporta-se da mesma

forma: representa fielmente os valores esperados em tempos cutos e longos, mas deixa a

desejar em tempos intermediários (veja-se a Figura 4.14). A posição de número 5, mais

próxima ao canto é aquela que apresenta maior discrepância, mesmo em tempos curtos.


97

Figura 4.14 – Exemplo 1: Quedas de pressão em outros pontos do reservatório.


98

Exemplo 2: Reservatório Quadrado – Fronteiras Seladas – Poço Descentralizado

Dando-se seqüência ao estudo, deseja-se neste exemplo, avaliar a influência da posição

do poço em relação aos limites do reservatório e sua relação com a exatidão da resposta.

Assumindo-se, inicialmente, a discretização mínima para a geometria em objeto e poço

posicionado a 1000m do limite mais próximo, se obtém o resultado apresentado na Figura

4.15.

A transição passa a ser mais longa do que a do exemplo anterior e, como se pode

admitir, mesmo com esta pobre representação, o método numérico é capaz de bem

representar os regimes de fluxo perpassados. As maiores discrepâncias permanecem na

região de transição entre eles.


Posicionando-se, posteriormente, o poço a 250m do limite mais próximo, conforme a

Figura 4.16, o código computacional ainda é capaz de produzir resultados satisfatórios. Ao

menos compatíveis com o tipo e faixa de erros identificados anteriormente.

Neste caso, em que as distâncias às barreiras passam a diferir em ordem de magnitude,

a derivada de pressões denota um regime de fluxo intermediário, denominado semi-radial,


99

correspondente ao segundo patamar97 da curva de derivada de pressões da figura em análise

no período compreendido entre 20 e 80h.

Figura 4.16 - Exemplo 2– NE= 4, NPGL=8. Persistindo-se com o mesmo exercício e dispondo o poço a apenas 50m da fronteira

mais próxima, finalmente se obtêm resultados espúrios. Com a representação do meio por

apenas 4 elementos o caráter esperado da resposta não é reproduzido. Veja-se a Figura

4.17.

Contudo, dispondo apenas 2 novos elementos de 10m nas imediações do poço os erros

voltam a faixa até agora considerada (veja-se a Figura 4.18).

97 Em valores numéricos, na proporção de 1 para 2, do primeiro patamar para o segundo.


100


Figura 4.18 - Exemplo 2– NE= 7 não uniformes, NPGL=8.


101

Exemplo 3: Reservatório Quadrado – 1 Fronteira Realimentada – Poço Descentralizado

Aplicando-se o aprendizado anterior e considerando-se, agora, uma das fronteiras como

realimentada, passa a ocorrer o comportamento exibido na Figura 4.19. Em termos práticos,

uma fronteira realimentada pode representar um possante aqüífero.

No caso em objeto são esperados três regimes de fluxo característicos: radial infinito,

semi-radial e regime permanente. Os dois primeiros já foram descritos anteriormente.

Quanto ao último, quando da influência da porção realimentada do reservatório, as quedas

de pressão no poço se tornam constantes e, como conseqüência, a derivada temporal passa a

tender a decrescer de forma abrupta.

Pela avaliação da Figura 4.19 podem-se constatar baixos erros relativos, sendo os

máximos sistematicamente ocorridos nos períodos correspondentes as transições entre

regimes característicos.

Resultado de semelhante qualidade é obtido ao se dividir o mesmo meio em 80

elementos uniformes. Veja-se a Figura 4.20.

Figura 4.19 - Exemplo 3 – NE=18 não uniformes, NPGL=8.


102

Figura 4.20 - Exemplo 3 – NE=80, NPGL=8.


103

Exemplo 4: Reservatório Quadrado – Poço Próximo a um Ângulo Reto

Aproximando-se o poço de um segundo limite também a 250m e lançando-se mão

do exposto anteriormente, com apenas 20 elementos de comprimento não uniforme, se

obtêm discrepâncias compatíveis com as até aqui exibidas.

Com essa configuração, o sistema atravessa três regimes distintos, conforme os

perpassados no exemplo 2. A única distinção está relacionada à proporção dos valores dos

dois patamares da curva de derivadas de pressão: desta vez, 1 para 4.

Figura 4.21 - Exemplo 4 – NE=20 não uniformes, NPGL=8. Exemplo 5: Reservatório Aberto– Poço Próximo Limites em Ângulo de 30º

Até aqui foram discutidos reservatórios de forma circular e quadrada – portanto,

com ângulos internos obtusos e retos. O intuito deste Exemplo 5 é avaliar como se

comporta a implementação do MEC quando da presença de ângulos agudos. Os parâmetros

de simulação são os apresentados na Tabela 4-5.


104

Tabela 4-5 – Exemplo 5: parâmetros considerados.

PARÂMETRO VALOR

k [mD] 300,0000

fi 0,2000

mi [cP] 1,0000

ct [cm²/kgf] 0,00010

h [m] 50,0000

q [m³/d] 500

rw [m] 0,10

A Figura 4.22 mostra que com elementos de extensão próxima a 2000m já é

possível representar a forma esperada da derivada de pressões, porém com erros máximos

absolutos da ordem de 11%.

Há dois regimes de fluxo característico estabelecidos, ambos de aspecto radial,

representados por dois patamares na proporção de 1 para 6.

Figura 4.22 – Exemplo 5 - NE=25, NPGL=8.


105

Duplicando-se a quantidade de elementos é possível trazer o erro máximo relativo

para a faixa com a qual se tem trabalhado até aqui.

Note-se que, neste caso, uma condição favorável é criada. Qual seja: o poço passa

estar localizado nas cercanias de uma extremidade de elemento e não mais nas

proximidades de um centro de elemento. Consulte-se o disposto no item 4.1.2.

Curiosamente, as discrepâncias também mudam de sinal.

Figura 4.23 – Exemplo 5 - NE=49, NPGL=8.


106

Exemplo 6: Reservatório Retangular – Poço entre Falhas Paralelas

Uma outra geometria de uso corrente e, por conseguinte, que deve ser bem

representada é a de fluxo em canal. Sistemas deposicionais turbidíticos típicos da costa

brasileira podem apresentar esta configuração.

Quando se trata de geometria acanalada duplamente simétrica98, e fechada em suas

extremidades são esperados três regimes característicos de fluxo: o radial infinito, o

linear99, de inclinação ½, e o pseudo-permanente, caso haja tempo suficiente para se

perceber todos os limites do reservatório.

Seja um reservatório retangular de lados 250m e 10.000m, conforme a Figura 4.24,

sem escala, tendo como parâmetros os especificados na Tabela 4-6.

Tabela 4-6 – Exemplo 6: parâmetros considerados.

PARÂMETRO VALOR

k [mD] 800,0000

fi 0,2000

mi [cP] 1,0000

ct [cm²/kgf] 0,00010

h [m] 50,0000

q [m³/d] 800

rw [m] 0,10

Em simulação via MEC, caso sejam adotados elementos uniformes de 2.000m na

direção x, o comportamento esperado não é reproduzido. Note-se que a distância do poço as

extremidades do canal é de 125m.

Já se dispondo de elementos de 500m ou de 250m o aspecto esperado da derivada

de pressões é alcançado e os erros são progressivamente reduzidos à faixa aceitável. Esta

afirmação pode ser comprovada através da Figura 4.25 e da Figura 4.26.

98 Simétrica em relação aos dois eixos coordenados. 99 Após o transitório atingir as fronteiras paralelas associadas ao canal as linhas de corrente equipotenciais tendem a se linearizar.


107

Figura 4.24 - Exemplo 6 – NE=12, NPGL=8.

Figura 4.25 - Exemplo 6 – NE=42, com comprimento uniforme na direção x, NPGL=8.


108

Figura 4.26 - Exemplo 6 – NE=84,com comprimento uniforme na direção x, NPGL=8.

Exemplo 7: Reservatório Infinito – Poço Fraturado

Uma outra característica dos programas computacionais desenvolvidos, que deve ser

explorada, é a de cômputo de mais de uma fonte de domínio. Para ilustrar esta

possibilidade, é escolhido o modelo de poço fraturado, sendo ele próprio representado por

um conjunto de fontes lineares100.

O fraturamento hidráulico é freqüentemente empregado na indústria de petróleo para

aumento de produtividade dos poços.

A Figura 4.27, a Figura 4.28 e a Figura 4.29 representam um reservatório de extensão

infinita sujeito a produção através de taxa constante de poço com fratura de 200m de

comprimento. Para os demais parâmetros considerados consulte-se a Tabela 4-7.

Transitórios correspondentes a poços fraturados tendem a apresentar, antes do regime

radial infinito, um regime de fluxo linearizado – de inclinação ½ no gráfico de diagnóstico.

100 Uma maneira mais elegante e mais exata para este cômputo seria a adoção da solução fundamental de fonte planar com extensão finita.


109

Os modelos, compatíveis com fluxo uniforme na admissão da fratura, são gerados com

21, 61 e 121 fontes, respectivamente.

Note-se que com poucas fontes, em tempos muito curtos, chega a ser desenvolvido um

regime radial a cada uma delas. À medida que a quantidade de fontes é aumentada este

regime vai desaparecendo. Com 121 fontes a resposta, se comparada à analítica

equivalente, pode ser considerada excelente.

O cálculo de fontes internas demanda pouquíssimo esforço computacional.

Tabela 4-7 - Exemplo 7: parâmetros considerados.

PARÂMETRO VALOR

k [mD] 50,0000 fi 0,2000

mi [cP] 1,0000 ct [cm²/kgf] 0,00010

h [m] 50,0000 q [m³/d] 800 Lf [m] 200,0000

Figura 4.27 – Exemplo 7 – NF=21, NPGL=8.


110

Figura 4.28 - Exemplo 7 – NF=61, NPGL=8.

Figura 4.29 – Exemplo 7 – NF=121, NPGL=8.


111

Seguindo-se esta mesma linha de raciocínio, poderia ter sido implementada solução

para poço horizontal101. Bastaria modificar a solução fundamental de Green de fonte linear

para fonte pontual e considerar as barreiras limitantes de topo e base de reservatório através

do método das imagens.

Exemplo 8: Modelo Composto – Duas Regiões Circulares Concêntricas

O modelo radial composto é definido por duas ou mais regiões circulares concêntricas,

sendo cada qual definida por k,µ,φ,ct e pelo seu raio externo. A razão entre os dois

primeiros parâmetros é denominada mobilidade, já o produto entre os dois últimos,

capacidade de estoque.

Este modelo é freqüentemente empregado em interpretação de testes em poços em

situações em que há presença de dois fluidos no meio poroso.

Note-se que deste ponto em diante, as respostas confrontadas com as soluções

analíticas advém do código preparado para lidar com domínios compostos de múltiplas

regiões homogêneas, ECDC.

Considere-se inicialmente apenas 2 regiões concêntricas, sendo a externa de extensão

infinita. Considere-se, também, que não haja contraste entre a capacidade de estoque das

regiões. Sejam os parâmetros básicos os apresentados na Tabela 4-8.

Tabela 4-8 - Exemplo 8: parâmetros considerados.

PARÂMETRO VALOR

k1/k2 [mD] 200/ VAR φ1/ φ2 15% /VAR

mi [cP] 1,0000 ct [cm²/kgf] 0,00010

h [m] 25,0000 q [m³/d] 800 rw [m] 0,1000

rI 50,0000

Deste modo, a derivada de pressões deve apresentar dois regimes de fluxo

característicos, ambos de aspecto radial – portanto, patamares. O primeiro relacionado à

mobilidade da região interna e o segundo relacionado à mobilidade da região externa.

101 Mesmo em se tratando de código preparado para lidar apenas com problemas bidimensionais.


112

O valor numérico de cada patamar é inversamente proporcional à mobilidade da região.

Em outras palavras, quanto mais baixo o patamar, melhor a permeabilidade da região

perpassada pelo transitório de pressões.

A Figura 4.30 mostra um caso em que a mobilidade do setor externo é 4 vezes inferior

a do setor interno.

A Figura 4.31 denota mobilidade do setor externo 2 vezes superior.

Para ambas as simulações são utilizados 16 elementos de interface e são obtidas

discrepâncias relativas máximas inferiores a 3%.

Figura 4.30 – Exemplo 8 - NE de interface = 16, NPGL=8.

k=200 mD

k=50 mD fi = 15%


113

Figura 4.31 – Exemplo 8 - NE de interface = 16, NPGL=8.

Adicionando-se aos dois exemplos anteriores um contraste de capacidade de estoque na

proporção de 10 para 1, ocorre acentuada modificação no período de transição entre os dois

patamares da curva de derivadas de pressão. O MEC as representa de forma fidedigna.

Vejam-se as 2 figuras seguintes.

O acréscimo de limite hexagonal selado à região externa, conforme Figura 4.33, gera

resultados com qualidade compatível com os até aqui manejados.

k=200 mD k=400 mD fi = 15%


114

Figura 4.32 - Exemplo 8 - NE de interface = 16, NPGL=8.

Figura 4.33 - Exemplo 8 - NE de interface = 16, NE de fronteira externa=6, NPGL=8.

k=200 mD

k=50 mD fi = 1,5%

k=200 mD

k=400 mD fi = 1,5%


115

Exemplo 9: Modelo Composto – Quatro Regiões Circulares Concêntricas

Neste exemplo são exploradas 4 regiões concêntricas ao invés de 2.

São adotados basicamente os mesmos parâmetros destacados na Tabela 4-8. Exceção

feita às distâncias radiais às interfaces – 50, 500 e 5000m – e às permeabilidades de cada

setor circular – 400, 200, 1600 e 800mD.

Na Figura 4.34 as interfaces são representadas por 8 elementos cada qual. Já na Figura

4.35 as duas interfaces internas são divididas em 16 elementos.

Os resultados das duas simulações são similares.

Figura 4.34 - Exemplo 9 - NE de interface = 8, NPGL=8.

k=400 mD

k=1600 mD

k=800 mD

k=200 mD


116

Figura 4.35 - Exemplo 9 - NE de interface = 2x16 + 8, NPGL=8. Exemplo 10: Modelo Linear Composto

Com o esgotamento dos modelos compostos disponíveis pelo PANSYSTEM102 e

tendo-se em mente que o modelo radial composto é caso particular do tipo de

implementação objeto desta dissertação, houve necessidade de se recorrer, primeiramente,

ao programa comercial de mesma finalidade denominado SAPHIR103, com vistas à

exploração do modelo linear composto.

Este modelo é ilimitado e possui uma interface semelhante à disposta no quadro menor

da Figura 4.36. São testadas duas relações entre propriedades – conforme a Tabela 4-9 - e

os resultados obtidos são muito bons.

Tabela 4-9 –Exemplo 10: parâmetros considerados

PARÂMETRO VALOR

k1/k2 (k2) [mD] 30,00/ 120,00 (10,00) fi 0,1000

mi [cP] 0,2500 ct [cm²/kgf] 0,00010

h [m] 25,0000 q [m³/d] 400 L [m] 300,0000

102 O PANSYSTEM não dispõe de modelo linear composto, com propriedades distintas entre regiões. 103 Versão 3.12. Para maiores informações, consultar o sítio: http://www.kappaeng.com/software/saphir.

k=400 mD

k=1600 mD

k=800 mD

k=200 mD


117

Figura 4.36 - Exemplo 10 - NE de interface = 25, sendo o menor com 10m. Exemplo 11: Modelo Composto por 3 regiões

Por último, o código computacional ECDC é testado contra o PANMESH, simulador

numérico baseado em elementos finitos, atrelado ao PANSYSTEM.

É avaliado o modelo disposto na Figura 4.37, sob as condições designadas na Tabela

4-10.

Tabela 4-10 –Exemplo 11: parâmetros considerados

PARÂMETRO VALOR

k1 /k2 /k3 [mD] 1000,00/ 2000,00 / 10,00 fi 0,20/ 0,20 / 0,05

mi [cP] 1,00 ct [cm²/kgf] 0,00010

h [m] 25,0000 q [m³/d] 500 L [m] 300,0000

Há um único poço, próximo a uma região de baixa permeabilidade (região 3).

Pode-se constatar que são obtidas respostas similares por ambos os métodos

numéricos.

K1=30mD K2=10mD

K2=120mD


118

Figura 4.37 - Exemplo 11 – modelo composto de 3 regiões.

4.4. Limitações do Algoritmo de Stehfest

Até este ponto, foi demonstrada a aptidão das implementações para lidar com todas as

características prometidas no item 1.3, com exceção da de produção das fontes internas ao

domínio por esquemas de vazões variáveis. Por outra, até aqui se considerou poços

produzindo a taxas constantes.

A produção de petróleo é geralmente regulada por estranguladores de fluxo. Esta

espécie de dispositivo é interposta à linha de produção, em superfície, a jusante da cabeça

de poço, com objetivo de regular pressões e vazões de todo o sistema produtivo.

Os estranguladores disponíveis têm a forma cilíndrica e seus diâmetros internos são

padronizados em determinadas bitolas. A substituição de um estrangulador por outro de

diferente diâmetro interno acarreta, portanto, mudança brusca (ou “em degrau”) na vazão

de produção.

MEF MEF


119

O algoritmo de STEHFEST (1970) não foi idealizado para o trato com funções

descontínuas104.

A Figura 4.38 consiste num exemplo ilustrativo: imagine-se que se queira inverter,

pelo algoritmo de STEHFEST (1970), do campo de Laplace para o campo real uma função

que sabidamente possui a natureza da representada pelos degraus em azul escuro.

Neste caso, após a ocorrência da primeira descontinuidade, em t=1, a resposta gerada

pelo algoritmo não é capaz de representar apropriadamente o segundo patamar.

Note-se, também, que o acréscimo progressivo no número de coeficientes de Stehfest

por passo de tempo tende a corresponder a um aumento de exatidão na resposta. No

entanto, a partir de NSteh (NS na figura) = 26, entra decisivamente em jogo a precisão da

própria máquina e a solução deixa de ser bem comportada.

Figura 4.38 - Exemplo sintético do uso do algorítmo de Stehfest. Caso de vazões variáveis.

Efetivamente, caso se queira avaliar a resposta de poços sujeitos a vazões variáveis, em

virtude da linearidade do problema, se pode obter a resposta isolada para cada poço

submetido à vazão constante, e aplicar, posteriormente, os Princípios da Superposição

dispostos nos itens 2.6 e 2.7.

104 De fato, isto só foi percebido durante o desenvolvimento do trabalho.


120

Naturalmente, não se trata de um procedimento desejado ou computacionalmente

eficiente, mas em certas circunstâncias pode ser administrado.

CAPÍTULO 5: CONCLUSÕES E SUGESTÕES

121

5.Conclusões e Sugestões

5.1. Conclusões

Em primeiro lugar, o estudo dos fundamentos envolvidos no desenvolvimento do

Método dos Elementos de Contorno abarca um útil e elegante ferramental matemático,

dilatador de horizontes. O postulante a dissecá-lo se depara com diversos conceitos

matemáticos estimulantes.

Adicionalmente, se pode afirmar que se trata de uma técnica viável para a solução do

tipo de problema estudado. Com alguns poucos cuidados, o MEC é capaz de representar o

caráter dos regimes de fluxo perpassados. Devem ser evitadas fontes muito próximas a

centros de elementos.

Relativamente poucos elementos são necessários para representar os simples problemas

apresentados. Até que se encontrem barreiras ou interfaces, a característica analítica da

solução é preservada. Os maiores erros se concentram nas transições entre regimes de fluxo

característicos. Nos diversos casos estudados, nestas porções foram obtidas discrepâncias

relativas inferiores a 6%, com poucos elementos. Neste sentido, a adoção da solução

numérica via espaço de Laplace é interessante na medida em que não acumula erros no

decorrer do tempo. O algoritmo de Stehfest é robusto e eficaz. Com 6 passos produz

excelentes resultados. Em contrapartida, não é apropriado para manusear históricos de

vazão variáveis.

Em relação a métodos analíticos, o MEC tem como vantagem a flexibilidade. Não são

esperadas restrições quanto à conformação das fronteiras e interfaces dos reservatórios.

Além disso, admite condições de contorno mais amplas e várias fontes no domínio.

Em relação a outros métodos numéricos (MDF e MEF), pode ser vantajoso, desde que

não haja não-linearidades. Em algumas circunstâncias tende a ser mais exato.


122

5.2. Sugestões

Entre os métodos numéricos mais difundidos - MDF, MEF e MEC - o MEC é aquele

cujo desenvolvimento é mais recente. Deste modo, constitui um campo ainda muito fértil

de trabalho.

Potenciais incrementos em relação ao que até aqui foi desenvolvido são:

Quanto à Exatidão da Resposta

• Pode-se implementar elementos de parâmetros lineares em lugar de elementos de

parâmetros constantes. Isto certamente se refletirá em ganhos de exatidão nos

cálculos dos coeficientes das matrizes H e G e dos valores em pontos internos.

Quanto à Aceleração da Solução

• Como se viu, a aceleração da solução de problemas da espécie analisada neste

trabalho – relativamente poucos elementos de contorno - está intimamente ligada ao

processo de montagem do sistema linear. Para outras circunstâncias, pode ser

interessante, no caso de múltiplas regiões homogêneas, que a matriz global nem

mesmo seja constituída e que seja operada alguma solução iterativa105, tal qual fez

DORS (2002).

• Uma outra vertente, para problemas com mais elementos e múltiplas regiões, pode

ser a reestruturação da matriz global para solução por método direto (veja-se

SANTIAGO, TELLES e VALENTIM, 1999).

Quanto ao Uso do Programa

• Especialmente para regiões compostas pode ser investido tempo na construção de

um pré-processador para que seja obtida maior facilidade na montagem do

problema.

Quanto ao método de inversão do campo de Laplace para o campo real

105 Por exemplo, via gradientes bi-conjugados.


123

• Apesar da escolha do algoritmo de STEHFEST (1970) ter sido considerada

acertada, do ponto de vista de sua robustez e da não acumulação de erros, resultou

insuficiente para cômputo de funções descontínuas. Uma promissora alternativa é o

algoritmo de Iseger (veja-se AL-AJMI et al, 2008).

Quanto à Inclusão de Novas Características Físicas e Geométricas no Simulador

Neste âmbito, pode-se buscar:

• Incluir de efeitos locais de poço – efeito de película106 e de estocagem107 – para

uma representação mais realista do problema prático;

• Estudar a adequação dos códigos para consideração de barreiras de extensão

limitada dentro do domínio, que possam representar falhas geológicas, condutivas,

selantes ou com imposição de capacidade de transmissão regulada;

• Reformular os programas para inclusão de escoamento de fluido muito

compressível – tipicamente, gás;

• Com conceito similar ao que foi utilizado para representação de poços fraturados,

manusear poços horizontais, mesmo com códigos preparados para problemas

bidimensionais. Basta substituir a solução fundamental de fonte linear pela de fonte

pontual e aplicar o conceito de imagens para representação de topo e base do

reservatório;

• Estender as aplicações para reservatórios estratificados, com cada substrato tendo

características petrofísicas distintas. Este tipo de reservatório ocasiona troca de

fluido entre camadas. É necessário avaliar a melhor forma de representá-la. Por

exemplo, se na própria solução fundamental ou com empilhamento de camadas.

Nesta segunda hipótese, fatalmente seria necessário dividir o domínio em

elementos;

106 A perfuração de poços de petróleo normalmente acarreta modificações das propriedades químico-físicas em suas cercanias. Estas alterações, por estarem restritas a uma região muito próxima à parede do poço e pela própria maneira de contabilizá-la nos cálculos, são embutidas nas soluções por uma parcela aditiva de queda de pressão, denominada efeito de película (ou skin). 107 Em instantes posteriores a abertura de poços de petróleo, antes mesmo da manifestação do reservatório de petróleo, ocorre descompressão dos fluidos no interior da tubulação. Esta descompressão é denominada efeito de estocagem.


124

• Expandir os conceitos aqui explorados para problemas tridimensionais e, em

decorrência, se avaliar, por exemplo, poços horizontais multifraturados;

• Representar escoamento bifásico e não linear água-óleo. Em campos de petróleo, a

injeção de água é largamente utilizada com a finalidade de manter a pressão de

reservatórios e de empurrar o fluido de interesse dos poços injetores em direção

aos produtores. É de fundamental importância que o MEC, mesmo que com algum

artifício, possa manusear este tipo de escoamento.

125

Referências Bibliográficas ABRAMOVITZ, M., STEGUN, I. A., 1972, Handbook of Mathematical Functions With

Formulas, Graphs, and Mathematical Tables, 10ª impressão, United States Department of

Commerce, National Bureau of Standards, Applied Mathematics Series 55.

AL-AJMI, N., AHMADI, M., OZKAN, E., KAZEMI H., 2008, “Numerical Inversion of

Laplace Transforms in the Solution of Transient Flow Problems With Discontinuities”,

SPE Annual Technical Conference and Exhibition, Denver, USA, Set.

BOURDET, D., WHITTLE, T.M., DOUGLAS, A.A., PIRARD, Y.M., 1983, “A New Set

of Type Curves Simplifies Well Test Analysis”, World Oil, Mai.

BOURDET, D., 2002, Well Test Analysis: The Use of Advanced Interpretation Models,

Handbook of Petroleum Exploration and Producion, nº 3, Amsterdam, Elsevier.

BREBBIA, C.A., 1980, The Boundary Element Method for Engineers, Pentech Press.

BREBBIA, C. A., TELLES, J. C. F., WROBEL, L. C., 1984, Boundary Element

Techniques, Berlin, Heidelberg, New York, Tokyo, Springer-Verlag.

CARSLAW, H.S., JAEGER, J.C., 1959, Conduction of Heat in Solids, London, Oxford

United Press.

CHURCHIILL, R.V., 1972, Operational Mathematics, New York, McGraw-Hill, 3rd

Edition.

DORS, C., 2002, Paralelização de Algoritmos MEC Via Subestruturação Baseada em

Solvers Iterativos – Aplicação a Problemas 3D Escalares e Vetoriais, Dissertação de

Mestrado, Programa de Pós-graduação em Engenharia Civil, UFOP, Escola de Minas,

Departamento de Engenharia Civil.

126

EARLOUGHER JR., R.C., 1977, Advances in Well Test Analysis, SPE Monograph Series,

Second Print, New York e Dallas, Society of Petroleum Engineers of AIME.

GRINGARTEN, A.C., RAMEY, H.J., 1973, “The Use of Source and Green’s Functions in

Solving Unsteady-Flow Problems in Reservoirs”, SPEJ, Out.

KIKANI, J., HORNE, R.N., 1992, “Pressure-Transient Analysis of Arbitrarily Shaped

Reservoirs With the Boundary-Element Method”, SPE Formation Evaluation, Mar, pp. 53-

60.

KIKANI, J., HORNE, R.N., 1993, “Modeling Pressure-Transient Behavior of Sectionally

Homogeneous Reservoirs by the Boundary-Element Method”, SPE Formation Evaluation,

Jun.

KRYUCHKOV, S., SANGER, S., BARDEN, R., 2001, “Implementation of the Boundary

Element Method in a Practical Reservoir Engineering Software Application”, Canadian

International Petroleum Conference, artigo 2001-001.

LIGGETT, J.A., LIU,P.L-F, 1979, “Unsteady Flow in Confined Aquifers: A Comparison of

Two Boundary Integral Methods”, Water Resources Res., Ago, 15, nº 4, pp 861-66.

MATSUKAWA, J., HORNE, R.N., 1988, “The Application of Boundary Integral Method

to Immiscible Displacement Problems”, SPERE, Ago, pp 1061-68.

NEWMAN, ALBERT B., 1936, “Heating and Cooling Rectangular and Cylindrical

Solids”, Ind. Eng. Chem., 28, no. 5, Feb.

OZKAN, E., RAGHAVAN, R., 1991a, “New solutions for Well-Test-Analysis Problems:

Part 1 – Analytical Considerations”, SPE Formation Evaluation, Set.

127

OZKAN, E., RAGHAVAN, R., 1991b, “New solutions for Well-Test-Analysis Problems:

Part 2 – Computational Considerations and Applications”, SPE Formation Evaluation, Set.

PISKUNOV, N., S/A, Cálculo Diferencial e Integral, Edições Cardoso.

RAGHAVAN, R., 1993, Well Test Analysis, New Jersey, PTR Prentice Hall.

RIZZO, F.J., SHIPPY, D.J., 1970, “A Method of Solution for Certain Problems of

Transient Heat Conduction”, AIAAJ, Nov, 8, nº 11, pp 2004-2009.

SANTIAGO, J.A.F., TELLES, J.C.F., VALENTIM, V.A., 1999, “An optimized block

matrix manipulation for boundary elements with subregions”, Advances in Engineering

Software, 30, pp 701-713.

STEHFEST, H., 1970, “Algorithm 368, Numerical Inversion of Laplace Transforms”,

Communications of ACM 13, N. 1, Jan, pp 47-49.

THOMAS, JOSÉ EDUARDO (org.), 2001, Fundamentos de Engenharia de Petróleo, Rio

de Janeiro, Editora Interciência Ltda.

VAN EVERDINGEN, A. F. e HURST, W., 1949, “The Application of the Laplace

Transformation to Flow Problems in Reservoirs”, Petroleum Transactions, AIME, Dez, pp

305-324.

WROBEL, L.C., 2002, The Boundary Element Method, Volume 1, Applications in

Thermo-Fluids and Acoustics, Brunel University, UK, John Wiley & Sons Ltda.

128

Apêndice A

Algoritmo de STEHFEST (1970) A inversão numérica de ( )sf no campo de Laplace para F(t) no campo real pode ser

aproximada pela expressão:

( )

≈ ∑=

it

fVt

tFNSteh

ii

2ln2ln

1

, (A.1)

onde os coeficientes Vi são definidos por:

( ) ( )( ) ( ) ( ) ( )( )[ ]

( )

∑+=

++

−−−−=

2/,min

2/1int2

2/12/

!2!!!2/

!21

ni

ik

nni

iikkikkn

kkV . (A.2)

Observações:

- NSteh deve ser par

- a soma de Vi deve ser aproximadamente igual a zero.

Tabela A-1: Coeficientes de Stehfest.

N Vi

2 6 10 12 1 2 1 0,083333 0,016667

2 -2 -49 -32,083333 -16,016667

3 366 1.279,000000 1.247,000000

4 -858 -15.623,666667 -27.554,333333

5 810 84.244,166667 263.280,833333

6 -270 -236.957,500000 -1.324.138,700000

7 375.911,666667 3.891.705,533333

8 -340.071,666667 -7.053.286,333333

9 164.062,500000 8.005.336,500000

10 -32.812,500000 -5.552.830,500000

11 2.155.507,200000

12 -359.251,200000

129

Apêndice B

Transformação de coordenadas de sistema ortotrópico para

sistema isotrópico equivalente

Seja o fluxo num meio poroso ortotrópico governado pela equação da difusão

t

u

y

u

k

k

x

u

k

k yx

∂∂=

∂∂+

∂∂

η1

2

2

2

2

, (B.1)

onde k é um valor de referência, obrigatoriamente o mesmo embutido em η .

O primeiro termo do primeiro membro pode ser expandido para

∂∂

∂∂

∂∂=

∂∂

x

x

x

u

xk

k

x

u

k

k xx ''2

2

, (B.2)

com x’ sendo o equivalente ao eixo original x para o sistema coordenado transformado.

Processando-se as derivadas parciais pela regra da cadeia para a multiplicação passa-se

a ter:

x

x

x

x

x

u

k

k

x

x

x

u

k

k

x

u

k

k xxx

∂∂

∂∂

∂∂+

∂∂

∂∂=

∂∂ ''

''

' 2

22

2

2

2

2

. (B.3)

Para que o termo em objeto mantenha o mesmo aspecto na equação governante, em

relação à equação original, é necessário que no sistema transformado:

1'

2

=

∂∂

x

x

k

kx e (B.4)

0''

2

2

=∂∂

∂∂

x

x

x

x

k

kx . (B.5)

130

De (B.4), portanto:

xk

k

x

x =∂∂ '

ou (B.6)

xk

kx

x

=' . (B.7)

Sendo x’ definido por (B.7), a condição (B.5) é automaticamente atendida.

Analogamente, para o eixo ortogonal y pode-se escrever:

yk

ky

y

=' . (B.8)

Deste modo, com a transformação de coordenadas definida por (B.7) e (B.8), um

sistema ortotrópico pode ser transformado num isotrópico equivalente.

Ademais, k pode assumir qualquer valor de referência. No entanto, por razões relativas

à interpretação de testes em poços de petróleo é conveniente que seja definido como:

yxkkk = . (B.9)

131

Apêndice C

Tabela de transformadas de Laplace utilizadas neste trabalho

Expressões extraídas de ABRAMOVITZ e STEGUN (1972).

Designação

(Dissertação/ ABRMVTZ et AL,1972) Função

Espaço de Laplace Função

Espaço Real

C.1 29.2.4

(diferenciação) ( ) ( )+− 0fsfs

dt

df

C.2 29.2.8

(convolução) ( ) ( )sfsf 21 ( ) ( ) τττ dftfff

t

2

0

121 ∫ −=∗

C.3 29.3.1 s

1 1

C.4 29.3.82 ske−

−

t

k

t

k

4exp

2

2

3π

C.5 29.3.83 skes

−1, 0≥k

t

kerfc

2

C.6 29.3.84 skes

−1, 0≥k

−

t

k

t 4exp

1 2

π

C.7 29.3.120 ( )skK0 , 0>k

−

t

k

t 4exp

21 2

132

Apêndice D

Solução fundamental: fonte instantânea planar

Seja uma fonte planar, de extensão infinita108, posicionada em ξ de forma ortogonal ao

eixo coordenado x, que submeta um meio igualmente infinito em todas as direções a uma

retirada instantânea de fluido, definida por um fluxo específico q~ .

Figura D.0.1 – Fonte Planar.

Se tal retirada é dominada por efeitos viscosos, então é regida pela equação da difusão

hidráulica, que para fluxo unidimensional pode ser expressa por:

02

2

=∂∂−

∂∂

t

u

x

uxη (D.1)

Considerando equilíbrio do meio prévio a retirada, se tem que:

( ) 00, ==txu (D.2)

Pela condição de espaço infinito, abordada na definição do problema:

( ) 0,lim =−∞→− txux ξξ (D.3)

108 Em y e z. Infinitesimal em x.

x

y

ξ

133

Considerando-se a lei de Darcy e valendo-se da simetria do problema proposto, a

condição de contorno interna pode ser expressa por:

wx

dx

ukq

w

Γ∂∂= ∫

Γ =ξµ2 (D.4)

Ademais, se a retirada se dá de forma igualmente distribuída pela fonte, então a

expressão (D.4), pode ser escrita como:

k

q

x

u

x 2

~µξ

=∂∂

=

(D.5)

A solução do problema descrito por (D.1) a (D.3) e (D.5) pode ser obtida a partir

de algumas maneiras distintas. Neste apêndice, como na maior parte do trabalho, é

explorada aquela que lança mão de transformadas de Laplace.

Procedendo-se assim e reescrevendo-se, portanto, o problema no espaço

transformado se tem:

02

2

=−∂∂

usx

uη (D.6)

Note-se que na equação da difusão, (D.6), já está embutida a condição inicial.

Nestes termos, a condição de contorno interna, no espaço transformado, passa a

ser:

sk

q

x

u

x

12µ

ξ

=∂∂

=

. (D.7)

As soluções da equação diferencial (D.6) tem a forma:

ηξηξ sxsx eCeCu −−− += 21 . (D.8)

Por conseguinte, a derivada espacial de (D.8) é:

ηξηξ ηη sxsx eCseCsx

u −−− −=∂∂

21 . (D.9)

134

Fazendo com que x se aproxime da posição da fonte planar, ξ , e pela substituição

de (D.7) em (D.9), intui-se que necessariamente C1=0 e C2 vale:

ssk

qC ηµ 1

2

~2 = . (D.10)

Substituindo-se (D.10) em (D.9), já com as considerações pregressas chega-se a:

ηξ

ηφsx

t

essc

qu −−=

2

1~ (D.11)

A expressão (D.11) corresponde à solução para a queda de pressão no espaço de

Laplace para uma fonte planar produzida a uma vazão específica constante.

Com o uso da equação (2.111), a solução para o impulso instantâneo é

prontamente obtida como:

xsx

xt

esc

qu ηξ

ηφ−−=

2

1~. (D.12)

Finalmente, de (C.6) se tem que a expressão equivalente a (D.12) no campo real,

que é:

( )t

x

xt

xetc

qu η

ξ

πηφ4

2

2

1~ −−= . (D.13)

Portanto:

( )( )

t

x

x

xet

txg ηξ

πηξ 4

2

2

1,,

−−= . (D.14)

135

Apêndice E Tabela de conversão de unidades

GRANDEZA SISTEMA

BRASILEIRO DE CAMPO

SISTEMA INTERNACIONAL

MULTIPLICAR POR

COMPRIMENTO metro [m] metro [m] 1,000 000 000

PERMEABILIDADE mili-Darcy [mD] metro quadrado

[m2] 9,869 232 667 E-16

PRESSÃO

quilograma-força por centímetro

quadrado [kgf/cm2]

Pascal [Pa] 9,806 652 000 E+4

TEMPO hora [h] segundo [s] 3,600 000 000 E+3

VAZÃO metro cúbico por

dia [m3/d] metro cúbico por segundo [m3/s]

1,157 407 407 E-5

VISCOSIDADE centi-Poise [cP] Pascal.segundo

[Pa.s] 1,000 000 000 E-3

136

Apêndice F Vetores Apontadores de ECDC Os arranjos globais em domínios compostos de múltiplas regiões são gerados NSteh

vezes por passo de tempo. No entanto, suas topologias, esparsas, são fixas.

No código computacional ECDC, a topologia dos arranjos globais é definida via 4

vetores apontadores (veja-se o Capítulo 3). Tais vetores são gerados da seguinte maneira:

e(j), reg1(j), reg2(j), j=1..et

‘entrada de dados.

‘e(j): número do elemento, reg1: 1ª região a que pertence o elemento, reg2: 2ª região a

que pertence o elemento, et: total de elementos.

‘se reg1=reg2 então elemento de contorno propriamente dito.

‘caso contrário, de interface.

rt = Max(reg1(j), reg2(j)), j=1..et

‘rt: número de regiões que possui o domínio composto.

k=0

Para j=1..rt

Para i=1..et

Se reg1(i)=j ou reg2(i)=j então

k=k+1

vapl(k)=e(i)

vapr(k)=j

fim

Próximo i

Próximo j

‘vapl(k): vetor apontador do elemento em que deve se encontrar o ponto-fonte.

‘vapr(k): vetor apontador da região sob avaliação.

137

z=0

para i=1..rt

para j=1..rt

para k=1 a et

se reg1(k)=i e reg2(k)=j então

z=z+1

vaph(z) = e(k)

se e(k)>et – eint então

z=z+1

vaph(z) = e(k)

fim

fim

Próximo k

Próximo j

Próximo i

‘vaph(k): vetor apontador de coluna da matriz H modificada.

‘eint: número total de elementos de interface.

k=0

para j=1..rt

para i=1..et

se reg1(i)=reg2(i) and reg2(i)=j então

k=k+1

vapg(k)=e(i)

fim

Próximo i

Próximo j

‘vapg(k): vetor apontador de coluna da matriz G modificada.

De posse de tais vetores, a montagem dos arranjos globais é imediata.

Documents

APLICAÇÃO DO MÉTODO DOS ELEMENTOS DE CONTORNO …objdig.ufrj.br/60/teses/coppe_m/ConradoKeidel.pdf · iii Keidel, Conrado Aplicação do Método dos Elementos de Contorno na Modelagem