TEORIA ALTERNATIVA PARA ESTABILIZAÇÃO DA PERDA DEINJETIVIDADE

FERNANDO DIOGO DE SIQUEIRA

UNIVERSIDADE ESTADUAL DO NORTE FLUMINENSE - UENF

MACAÉ - RJ

FEVEREIRO - 2010

TEORIA ALTERNATIVA PARA ESTABILIZAÇÃO DA PERDA DEINJETIVIDADE

FERNANDO DIOGO DE SIQUEIRA

Dissertação apresentada ao Centro de Ciênciase Tecnologia da Universidade Estadual do NorteFluminense, como parte das exigências paraobtenção do título de Mestre em Engenharia deReservatório e de Exploração de Petróleo.

Orientador: Prof. Pavel Bedrikovetsky

UNIVERSIDADE ESTADUAL DO NORTE FLUMINENSE - UENF

MACAÉ - RJ

FEVEREIRO - 2010

TEORIA ALTERNATIVA PARA ESTABILIZAÇÃO DA PERDA DEINJETIVIDADE

FERNANDO DIOGO DE SIQUEIRA

Dissertação apresentada ao Centro de Ciênciase Tecnologia da Universidade Estadual do NorteFluminense, como parte das exigências paraobtenção do título de Mestre em Engenharia deReservatório e de Exploração de Petróleo.

Aprovada em 19 de Fevereiro de 2010.

BANCA EXAMINADORA

Prof. Dr. Carlos Enrique Pico Ortiz, D. Sc. - LENEP/CCT/UENF - Co-Orientador

Dr. Alexandre Sérvulo Lima Vaz Júnior, D. Sc. - UENF/CCT/LENEP

Dr. Andre Leibsohn Martins, D. Sc. - Petrobras/CENPES

Dr. Ronaldo Oliveira de Paiva, D. Sc. - Petrobras/UNBC

Prof. PhD. Pavel G. Bedrikovetsky - ASP/University of Adelaide/Austrália - Orientador

À Ana Catarina

ii

AGRADECIMENTOS

À Deus, por ter guiado meus passos e me abençoado.

À minha família, por ter me estimulado a seguir em frente. À minha ex-namorada Ana Catarina

pela compreensão, carinho, companherismo...

Ao Prof. Pavel Bedrikovetsky, que mais do que me orientar, me ensinou o significado da pesquisa

científica.

Aos membros do LENEP, especialmente ao prof. Adolfo por ter me ajudado sempre que necessário,

à Alexandre Sérvulo por diversas sugestões ao trabalho, e aos colegas pelo incentivo e pela amizade.

Ao CENPES pela disponibilização dos testes e pela colaboração no decorrer do trabalho, espe-

cialmente à Antônio Luíz Serra de Souza e Cláudio Alves Furtado.

À Banca Examinadora, especialmente ao Dr. André Leibsohn Martins, pela disponibilidade do

dia da defesa (na sexta-feira logo após o carnaval).

À ANP pela concessão da bolsa do programa de recursos humanos MCT-ANP-PRH20 durante

todo o mestrado.

SUMÁRIO

LISTA DE TABELAS vii

LISTA DE FIGURAS viii

SIMBOLOGIA xi

RESUMO xv

ABSTRACT xvi

1. INTRODUÇÃO 1

1.1 Objetivos . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 2

1.2 Organização do Documento . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 2

2. REVISÃO BIBLIOGRÁFICA 4

2.1 O Modelo Clássico de Filtração Profunda com Liberação de Partículas . . . . . . . . 5

2.1.1 Balanço de Massa . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 7

2.1.2 Cinética de Captura . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 8

2.1.3 Mecanismos de Captura . . . . . . . . . . . . . . . . . . . . . . . . . . . . 9

2.1.4 Crescimento do Diferencial de Pressão . . . . . . . . . . . . . . . . . . . . . 10

2.1.5 Adimensionalização do Sistema (1), (7), (10) . . . . . . . . . . . . . . . . . 11

2.2 Outros Modelos de Filtração Profunda . . . . . . . . . . . . . . . . . . . . . . . . . 12

2.3 Forças em Micro Escala . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 13

2.3.1 Força Eletrostática . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 13

2.3.2 Força de Arraste . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 14

2.3.3 Força de Elevação (“lifting”) . . . . . . . . . . . . . . . . . . . . . . . . . . 16

2.3.4 Força Gravitacional . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 16

3. METODOLOGIA 17

3.1 Revisão das Equações de Filtração Profunda . . . . . . . . . . . . . . . . . . . . . . 17

3.1.1 O Conceito de Capacidade de Estocagem . . . . . . . . . . . . . . . . . . . 17

3.1.2 Equações de Filtração Profunda Modificadas . . . . . . . . . . . . . . . . . 20

3.2 Bases Microscópicas do Conceito de Capacidade de Estocagem . . . . . . . . . . . . 21

3.3 Análise Dimensional de σcr . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 22

4. DESENVOLVIMENTO 25

4.1 Geometria Linear . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 25

4.1.1 Um Mecanismo . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 25

4.1.2 Dois mecanismos . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 35

4.2 Geometria Radial . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 39

5. RESULTADOS E DISCUSSÃO 43

5.1 Resultados para o Caso com um Mecanismo de Captura em Testemunhos . . . . . . 43

5.2 Resultados do Teste Vazão Alternada . . . . . . . . . . . . . . . . . . . . . . . . . . 46

5.3 Curvas Teóricas para σcr . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 48

5.4 Previsão de Impedância para Testemunhos com Dois Mecanismos de Captura . . . . 63

5.5 Previsão de Impedância para Poço . . . . . . . . . . . . . . . . . . . . . . . . . . . 65

6. CONCLUSÕES 67

Trabalhos Publicados . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 67

7. RECOMENDAÇÕES 68

REFERÊNCIAS BIBLIOGRÁFICAS 69

APÊNDICE A. CASO LINEAR A DOIS MECANISMOS COM EROSÃO 74

A.1 Aplicação do Lema de Agregação . . . . . . . . . . . . . . . . . . . . . . . . . . . 74

v

A.2 Obtenção da Equação Análoga à (175) em Função de σD Somente . . . . . . . . . . 75

A.3 Solução do Sistema (175)-(180) Utilizando o Método das Características . . . . . . . 77

A.3.1 Região 0 . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 77

A.3.2 Região I . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 78

A.3.3 Região II . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 79

A.3.4 Região III . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 80

A.4 Impedância . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 81

APÊNDICE B. CASO RADIAL COM EROSÃO 82

B.1 Obtenção da Equação Análoga à (225) em Função de σ Somente . . . . . . . . . . . 82

B.2 Solução do Sistema (225)-(229) pelo Método das Características . . . . . . . . . . . 83

B.2.1 Região 0 . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 83

B.2.2 Região I . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 84

B.2.3 Região II . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 85

B.2.4 Região III . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 85

B.3 Impedância . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 86

APÊNDICE C. FERRAMENTAS MATEMÁTICAS 87

C.1 Spline Cúbica Interpolante . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 87

C.2 Runge-Kutta de 4a ordem . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 88

vi

LISTA DE TABELAS

Tabela 1: Dados do teste e parâmetros de filtração . . . . . . . . . . . . . . . . . . . . . . . 44

Tabela 2: Dados do teste de vazão alternada . . . . . . . . . . . . . . . . . . . . . . . . . . 47

Tabela 3: Valores das constantes usadas em escala microscópica . . . . . . . . . . . . . . . . 49

Tabela 4: Parâmetros da simulação de impedância para o caso linear à dois mecanismos . . . 64

Tabela 5: Parâmetros usados na simulação de impedância para poços . . . . . . . . . . . . . 65

LISTA DE FIGURAS

Figura 1: Elemento de Meio Poroso (adaptado de Herzig et al. (1970)). . . . . . . . . . . . . 7

Figura 2: Gráficos de potencial total (a) e força elétrica (b) . . . . . . . . . . . . . . . . . . . 15

Figura 3: Esquema da "última" partícula do reboco interno . . . . . . . . . . . . . . . . . . 17

Figura 4: Concentração máxima retida em função da velocidade . . . . . . . . . . . . . . . . 19

Figura 5: Históricos esquemáticos de (a) velocidade e (b) concentração no efluente em exper-

imento hipotético . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 19

Figura 6: Idealização do meio poroso como um feixe de capilares em perfil (a) e em secção

transversal (b) . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 21

Figura 7: Transporte de suspensão em meio poroso com um mecanismo de captura e liberação

de partículas retidas: (a) movimento das frentes de concentração no plano (xD, tD),

com os pontos P = (xD, tD) e P’ = (xcrD (ηD) ,ηD); (b) perfil de retenção em quatro

momentos e (c) perfil de concentração da suspensão nesses quatro momentos . . . 26

Figura 8: Estágios de filtração profunda com erosão; histórico de concentração na saída do

testemunho (a) e curva de impedância (b) . . . . . . . . . . . . . . . . . . . . . . 30

Figura 9: Teste típico de vazões alternadas: (a) esquema de velocidade, (b) histórico de queda

de pressão, (c) histórico de impedância e (d) histórico de concentração no efluente . 32

Figura 10: Diagrama Caraterístico (a), perfil de retenção (b) e de suspensão (c) de partículas . 34

Figura 11: Diagrama (xD, tD) para o sistema linear com dois mecanismos . . . . . . . . . . . 37

Figura 12: Diagrama (r, t) para a geometria radial . . . . . . . . . . . . . . . . . . . . . . . . 40

Figura 13: Curvas de impedância para o primeiro (a) e segundo (b) testes. Os losangos denotam

os pontos experimentais e a linha contínua o ajuste realizado . . . . . . . . . . . . 45

Figura 14: Curvas de concentração em suspensão reduzida na saida do testemunho para o

primeiro (a) e segundo (b) testes . . . . . . . . . . . . . . . . . . . . . . . . . . . 46

Figura 15: Dados de pressão do teste de vazão alternada . . . . . . . . . . . . . . . . . . . . . 47

Figura 16: Ajuste dos dados do teste de vazão alternada: (a) ajuste da curva prevista (contínua)

contra os pontos laboratoriais e (b) concentração prevista da suspensão na saída do

testemunho . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 48

Figura 17: Curvas de concentração máxima retida em função da velocidade (a) e do número de

erosão (b), com análise de sensibilidade em relação à constante de Hammaker para

o primeiro teste . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 50

Figura 18: Curvas de concentração máxima retida em função da velocidade (a) e do número

de erosão (b), com análise de sensibilidade em relação ao raio da partícula para o

primeiro teste . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 51

Figura 19: Curvas de concentração máxima retida em função da velocidade (a) e do número de

erosão (b), com análise de sensibilidade em relação à constante de elevação χ para

o primeiro teste . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 52

Figura 20: Curvas de concentração máxima retida em função da velocidade (a) e do número de

erosão (b), com análise de sensibilidade em relação à constante de elevação ω para

o primeiro teste . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 53

Figura 21: Curvas de concentração máxima retida em função da velocidade (a) e do número de

erosão (b), com análise de sensibilidade em relação à constante de Hammaker para

o segundo teste . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 54

Figura 22: Curvas de concentração máxima retida em função da velocidade (a) e do número

de erosão (b), com análise de sensibilidade em relação ao raio da partícula para o

segundo teste . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 55

Figura 23: Curvas de concentração máxima retida em função da velocidade (a) e do número de

erosão (b), com análise de sensibilidade em relação à constante de elevação χ para

o segundo teste . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 56

Figura 24: Curvas de concentração máxima retida em função da velocidade (a) e do número de

erosão (b), com análise de sensibilidade em relação à constante de elevação ω para

o segundo teste . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 57

Figura 25: Curvas de concentração máxima retida em função da velocidade (a) e do número de

erosão (b), com análise de sensibilidade em relação à constante de Hammaker para

o teste de vazões alternadas . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 58

Figura 26: Curvas de concentração máxima retida em função da velocidade (a) e do número de

erosão (b), com análise de sensibilidade em relação ao raio da partícula para o teste

de vazões alternadas . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 59

ix

Figura 27: Curvas de concentração máxima retida em função da velocidade (a) e do número de

erosão (b), com análise de sensibilidade em relação à constante de elevação χ para

o teste de vazões alternadas . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 60

Figura 28: Curvas de concentração máxima retida em função da velocidade (a) e do número de

erosão (b), com análise de sensibilidade em relação à constante de elevação ω para

o teste de vazões alternadas . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 61

Figura 29: Impedância para o caso de dois mecanismos com erosão em testemunhos conforme

os parâmetros da Tabela 4) . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 64

Figura 30: Impedância para um poço sob as condições da Tabela 5 . . . . . . . . . . . . . . . 65

x

SIMBOLOGIA

Alfabeto Latino

A123 Constante de Hammaker [J]

Cm Concentração Molar [mol/L]

c Concentração de partículas em suspenção

cD Concentração de partículas em suspenção reduzida

c+D Concentração de partículas em suspenção reduzida à direita da frente de descontinuidade

(erosão)

c−D Concentração de partículas em suspenção reduzida à esquerda da frente de descontinuidade

(erosão)

D Constante Dielétrica, Velocidade da frente de descontinuidade

F Força [N]

f Densidade de distribuição [1/m]

H Espessura do canal do poro [m]

h Distância de separação entre grão e partícula [m]

J Impedância

J∗ Pseudo-impedância

k0 Permeabilidade inicial [m2, mD]

kdet Taxa de liberação de partículas [s−1]

k(σ) Fator de redução da permeabilidade adimensional devido à retenção de partículas

kB Constante de Boltzmann [kg m−2 s−2 K−1]

L Comprimento (do testemunho) [m]

m Coeficiente angular

n Concentração numérica de poros na rocha, [m−2]

p Pressão [Pa]

q Vazão de injeção (em poços) [m−3/dia]

r Raio [m]

T Temperatura Absoluta [K]

t Tempo [s]

tD Tempo adimensional [VPI]

t∗D Tempo adimensional de troca de vazão [VPI]

U Velocidade Darcy [m/s]

U Velocidade Intersticial [m/s]

V Energia de Interação [J]

x Coordenada cartesiana [m]

xD Coordenada cartesiana adimensional

Z Razão entre a separação grão-poro e o raio da partícula

z Valência

Alfabeto Grego

β Coeficiente de Dano à Formação

ε Razão entre forças de arraste e elétrica, número de erosão

ε0 Permissividade do vácuo [C−2 J−1 m−1]

ξ Intercepto da Característica com o eixo x

ξD Intercepto da Característica com o eixo xD

κ Inverso do comprimento de Debye, [m−1]

κC Coeficiente de atrito seco de Coulomb

xii

λ Coeficiente de filtração profunda, [1/m]

µ Viscosidade dinâmica [Pa.s]

v Concentração numérica de íons [m−3]

ρ Massa específica do Fluido [kg/m3]

ρpart Massa específica da Partícula [kg/m3]

∆ρ Diferença entre as massas específicas da partícula e da fluido [kg/m3]

η Valor temporal do intercepto da característica [s]

ηD Valor temporal adimensional do intercepto da característica

Π j Parâmetro adimensional referente à quantidade j

σ Concentração de partículas retidas

~σ Vetor de concentrações retidas por cada um dos mecanismo

σcr Concentração de partículas retidas máxima

σD Concentração de partículas retidas reduzida

σcrD Concentração de partículas retidas reduzida máxima

σ+D Concentração de partículas retidas reduzida à direita da frente de descontinuidade (erosão)

∆σcrD Diferença entre a Concentração máxima retida reduzida atual e a aquela da vazão anterior

σLJ Diâmetro de colisão atômico [m]

Φ Função

φ Porosidade

χ Coeficiente de elevação

ψ Potencial da superfície [mV]

ω Coeficiente da força de arraste

xiii

Outros Símbolos

∇ Operador Nabla

Subscritos

BR Born (energia potencial)

c Reboco

cr Máxima (concentração retida)

cf Arraste

D Adimensional

DLR Dupla camada elétrica repulsiva (energia potencial)

e Elétrico/eletrostático; Externo; erosão

l Elevação

g Gravidade

i Índice para íons

LVA London-van der Waals (energia potencial)

n Normal (para força)

s Partículas suspensas

p Poro

0 Inicial

det Liberação

123 Interação entre material 1 com 3 no meio 2

LJ Referente a Lennard-Jones (Potêncial de Born)

xiv

RESUMO

A liberação de partículas da matriz durante o transporte de suspensões em meios porosos é um

fenômeno amplamente observado em testes laboratoriais e em reservatórios naturais. Um novo mod-

elo para liberação de partículas é proposto, baseado no equilíbrio mecânico da partícula posicionada

no reboco interno dentro do poro. Este equilíbrio, que envolve os torques das forças de arraste,

eletrostática e de elevação agindo sobre a partícula na matriz, determina a concentração máxima de

retenção do processo de captura, e é dependente da razão adimensional entre as forças de arraste

e normal atuantes na partícula na superfície do reboco interno. A concentração máxima retida em

função do número de erosão, que é a razão entre as forças de arraste e elétrica, completa o sistema de

equações básicas para a filtragem de colóides com liberação de partículas.

O sistema resultante para o caso unidimensional linear com um mecanismo de captura admite

solução analítica desde que sejam constantes o coeficiente de filtração e a porosidade. Estas fórmulas

explícitas permitem o cálculo dos parâmetros do modelo (concentração máxima retida, coeficientes

de filtração e dano à formação) a partir de dados de diferencial de pressão ao longo do testemunho du-

rante a injeção da suspensão. Os valores de concentração retida encontrados por ajuste de testes labo-

ratoriais usando essas fórmulas foram comparados com o balanço de torques em escala microscópica,

para validar a teoria proposta. Soluções semi-analíticas para geometria linear na presença de dois

mecanismos e a solução do caso radial (com apenas um mecanismo) também serão apresentadas

aqui.

Palavras chave: Suspensões; Colóides; Liberação de Partículas; DLVO; Injetividade; Filtração;

Escoamento em Meio Poroso; Escoamento de Suspensões; Transporte de Fluidos com Sólidos

ABSTRACT

Particle detachment from matrix during suspension transport in porous media is a widely observed

phenomena in laboratory and natural reservoirs. A new model for particle detachment is proposed,

based on mechanical equilibrium of particle located on pore internal cake. This equilibrium, which

involves drag, electrostatic and lifting force torque acting on matrix particle, determines the maximum

retention concentration of capture process, and is function of the dimensionless ratio between drag

and normal forces acting on internal cake surface particle. The maximum retained concentration as

function of dislodging number, which is the ratio among drag and electrical forces, closes governing

equation system for colloid transport with particle release.

The resulting system for linear coreflood with one capture mechanism admits analytical solution

under the assumption of constant porosity and filtration coefficient. Those explicit formulae allow

model parameters (maximum retention concentration, filtration and formation damage coefficients)

calculation from pressure drop data across the core during suspension injection. The values of re-

tention concentration found in laboratory tests using those formulae are compared to torque balance

in micro scale, in order to validate theory. Semi-analytical solutions for the linear geometry, two

mechanisms case and radial geometry (one mechanism only) are also presented.

Keywords: Suspension; Colloids; Particle detachment; DLVO; Injectivity; Filtration; Flow in

Porous Media; Suspension Flow; Fluid Transport with Solids

1

1. INTRODUÇÃO

Na maioria dos campos de petróleo se faz necessário complementar a energia natural do reser-

vatório, o que normalmente é feito através da injeção de fluidos. Neste sentido, o método de recu-

peração mais comum é a injeção de água para deslocamento imiscível do óleo. Embora a eficiência

máxima desse deslocamento na escala do reservatório seja um problema relativamente bem conhecido

(BEDRIKOVETSKY, 1992), o fenomêno da Perda de Injetividade ainda requer alguma elucidação.

Perda de Injetividade é o aumento progressivo da resistência hidráulica (redução da permeabili-

dade) que se dá no decorrer da injeção. Este aumento se deve à captura das diversas partículas líquidas

e/ou sólidas presentes na água. Ao serem capturadas, estas partículas obliteram os poros conectados,

reduzindo seu diamêtro efetivo e por conseguinte aumentando o diferencial de pressão requerido para

sustentar uma dada vazão.

Obviamente, o aumento da pressão de injeção com o tempo requer o dimensionamento tão pre-

ciso quanto possível do aparato relacionado. Em campos offshore, onde os recursos são críticos, o

diagnóstico correto sobre as principais causas da Perda de Injetividade é fundamental, não só para

a adequada especificação dos equipamentos, mas para o projeto de injeção como um todo. Dados

recentes (PLUEN, 2007) mostram que o custo associado à produção do óleo é praticamente o mesmo

que o custo de manejo da água. Não obstante, a legislação ambiental, em virtualmente todos os

lugares do mundo, está ficando cada vez mais restritiva ao descarte da água produzida. Esta ex-

igência ambiental, bem como a disponibilidade, fazem da água produzida uma forte candidata para

a reinjeção. Algumas companhias inclusive consideram a possibilidade de injeção dessa água em

formações depletadas/abandonadas e até mesmo em aquíferos.

Por essas razões se fazem necessários modelos capazes de prever o comportamento da pressão

de injeção. Com isso em mente foi construído o assim chamado modelo matemático clássico de

perda de injetividade (ALTOÉ et al., 2004; CIVAN, 2007; HERZIG et al., 1970). De fato, o modelo

clássico, ou suas variantes, são capazes de explicar uma boa parte dos dados de campo, não só das

bacias brasileiras como do resto do mundo. Recentemente, porém, surgiram dados que não podem

ser explicados, sequer qualitativamente, por este modelo. Pesam contra ele os argumentos teóricos

descritos a seguir.

Um dos argumentos teóricos contra o modelo clássico é que ele envolve mecanismos demais. Em

outras palavras, o modelo clássico agrega diferentes formas de dano a formação num único corpo

que não é tão coeso quanto deveria. Nesse sentido, postula-se que ao invés desses mecanismos, deve

2

haver alguma mudança estrutural do meio poroso que explique a Perda de Injetividade.

Além disso, foram observados testes laboratoriais (FURTADO et al., 2002, 2003) nos quais a

injeção de água culmina na estabilização da permeabilidade. Porém, ao contrário do que ocorre em

poços, não foi observado reboco externo nestes testes. Ainda que isto ocorresse, deveria haver erosão

neste reboco para que a queda de injetividade se estabilizasse, o que não pode ocorrer nestes testes

porque não há fluxo paralelo ao reboco.

Comportamento ainda mais peculiar ocorre caso, após a estabilização, seja alterada a vazão de

injeção. Caso ocorra o aumento da vazão, é observado um aumento na concentração de partículas

a jusante do testemunho, além de uma diminuição na pressão de injeção. Caso ocorra diminuição

na vazão, ocorre novamente filtragem das partículas pelo meio poroso, com consequente perda de

injetividade até posterior estabilização. Com efeito, isto é bastante diferente do esperado levando-se

em conta o modelo clássico.

1.1 Objetivos

O objetivo deste trabalho é construir uma teoria para a estabilização da injetividade quando o

reboco externo não se forma. Para tal, postula-se que a concentração máxima de partículas retidas

é função de um tipo de velocidade de fluxo adimensional, o número de erosão. Esta função, que

depende também do meio poroso e da suspensão a ser filtrada pelo meio poroso, pode ser entendida a

partir do balanço de torques das forças atuantes nas partículas da superfície do reboco interno.

Com a introdução da função descrita no parágrafo anterior, se faz possível atacar analiticamente

o problema unidimensional linear à um mecanismo de captura. É possível também resover, semi-

analiticamente, os casos unidimensional linear com dois mecanismos e o radial com um mecanismo,

desde que o coeficiente de filtração e a porosidade sejam constantes. A partir das soluções encontradas

para o caso linear à um mecanismo é possível analisar dados experimentais, obtendo-se assim os

parâmetros do modelo apenas com dados de diferencial de pressão. A concentração máxima retida

calculada dos experimentos pode então ser comparada com o balanço de torques das forças atuantes

nas partículas do reboco interno, assumindo-se uma geometria de meio poroso simplificada. Este

último procedimento é suficiente para validar, ao menos qualitativamente, a teoria proposta.

1.2 Organização do Documento

No Capítulo 2, “REVISÃO BIBLIOGRÀFICA”, apresenta-se uma breve exposição dos modelos

clássicos de filtração profunda e das forças em micro-escala.

No Capítulo 3, “METODOLOGIA”, apresenta-se a base matemática e a formulação dos conceitos

3

fundamentais à teoria proposta nesta dissertação.

No Capítulo 4, “DESENVOLVIMENTO”, apresenta-se a solução do sistema exposto em "METODOLO-

GIA" para os casos radial e linear à um e dois mecanismos

No Capítulo 5, “RESULTADOS”, apresentam-se os resultados mais pertinentes do trabalho,

como os ajustes de dados laboratoriais de pressão e as estimativas de concentração retida máxima.

No Capítulo 6, “CONCLUSÕES”, será feita uma digressão sobre os fatos mais importantes ap-

resentados pelo modelo.

No Capítulo 7, “RECOMENDAÇÕES”, serã feitas sugestões para trabalhos futuros.

Apresenta-se a seguir as “REFERÊNCIAS BIBLIOGRÁFICAS” e os “APÊNDICES”:

O APÊNDICE A, detalha a solução para o sistema Linear à Dois Mecanismos com Erosão.

O APÊNDICE B, mostra as minúncias da solução do problema Radial com Erosão.

O APÊNDICE C, exibe os processos matemáticos auxiliares usados nas soluções semi-analíticas

empregadas nesta dissertação.

4

2. REVISÃO BIBLIOGRÁFICA

O transporte de suspensões e emulsões em meios porosos ocorre em inúmeros processos das en-

genharias ambiental, química, petrolífera e civil, entre outros ramos industriais. O fluxo de partículas

sólidas e líquidas e sua captura e liberação ocorrem na propagação de vírus, bactérias e gotículas

de óleo em aquíferos, filtragem industrial de líquidos e gases, cromatografia, injeção de água do

mar ou produzida em campos de petróleo, invasão de filtrado em formações quando da sua per-

furação, migração de finos em reservatórios de petróleo e aquíferos, entre outros processos (ENTOV;

MIRZHADZHANZADE, 1990; ELIMELECH et al., 1995; KHILAR; FOGLER, 1998; CIVAN,

2007; MAYS; HUNT, 2005). A retenção e liberação de partículas são fundamentalmente importantes

para processos ambientais, onde a concentração de partículas não deve exceder um determinado lim-

ite, enquanto a alteração da permeabilidade é importante para a produção de petróleo, devido à seus

efeitos na produtividade e injetividade de poços. A modelagem matemática da filtração profunda,

considerando a captura e liberação de partículas, é parte essencial no projeto de todos os processos

mencionados anteriormente.

Para melhor compreender a natureza dos fluidos injetados nos meios porosos com respeito ao ma-

terial particulado nele disseminado, se faz necessário classifica-los quanto ao tamanho dessas partícu-

las, o que é usualmente feito dividindo-os em duas classes, segundo Elimelech et al. (1995): as

soluções coloidais e as suspensões. As soluções coloidais apresentam partículas menores que 1 mi-

crometro, enquanto nas suspensões as partículas são maiores que esse valor. As razões pelas quais o

valor de 1 micrometro foi adotado como limite também estão esmiuçadas em Elimelech et al. (1995).

Neste trabalho, porém, não será feita distinção entre estas duas, e usar-se-ão ambos os termos como

sinônimos.

Pode-se avaliar a migração, retenção e liberação de colóides, tanto em laboratório quanto em es-

cala de campo, a partir de uma equação de balanço de massa, com termo de sumidouro modelando

a retenção e um termo de fonte para considerar a liberação (LOGAN, 2001; FOPPEN; SCHIJVEN,

2006; FRIMMEL et al., 2007). O termo de captura é proporcional ao fluxo advectivo da partícula,

sendo a constante de proporcionalidade chamanda de coeficiente de filtração. O termo de liberação é

proporcional a concentração de partículas retidas, com seu respectivo coeficiente denominado de taxa

de liberação. A teoria avançada da dependência do coeficiente de filtração com a difusão Browniana,

sedimentação gravitacional e velocidade da suspensão já se encontra desenvolvida Tufenkji e Elim-

elech (2004a), enquanto o coeficiente de liberação é determinado, na maioria dos casos, a partir de

5

dados experimentais.

Outra limitação desse modelo com cinética de liberação de partículas é a estabilização assintótica

da concentração retida (e por conseguinte da permeabilidade) em respeito ao tempo, o que contradiz

a observações feitas em migração de finos devido ao incremento no diferencial de pressão ou sob

mudança na salinidade, que mostram estabilização praticamente instantânea (KHILAR; FOGLER,

1998; MIRANDA; UNDERDOWN, 1993). Testes em testemunhos mostram incremento imediato de

permeabilidade após aumento abrupto da vazão (FURTADO et al., 2005).

Há muito foi reconhecido que a liberação de partículas ocorre se o equilíbrio mecânico da partícula

retida no reboco interno não é satisfeito. As forças atuantes nessa partícula são eletrostáticas, de ar-

raste, de elevação e gravitacionais. Alguns autores consideram o balanço de forças entre a força de

arraste, gerada pelo fluxo do fluido, e a força de fricção com um coeficiente de Coulomb empírico

(CIVAN, 2007). Outra abordagem considera o momento dessas forças (JIAO; SHARMA, 1994; FRE-

ITAS; SHARMA, 2001). Ainda que ambas as abordagens sejam equivalentes, a equação advectiva-

difusiva com termo cinético de liberação não reflete o equilíbrio mecânico da partícula, já que esse

termo não é afetado pelo equilíbrio em questão. Nenhum modelo em macro-escala, até onde vai o

conhecimento do autor, é encontrado na literatura considerando a influência do equilíbrio mecânico

na partícula na liberação das partículas em meios porosos.

A seguir serão descritos os modelos de filtração profunda. Há de se fazer uma exposição detal-

hada do modelo de filtração profunda clássico, englobando todos os diversos mecanismos que ele

contempla, bem como a modelagem usual de liberação de partículas. Após, serão expostas as teorias

mais recentes sobre filtração profunda. Por último, serão descritas as forças atuantes na partícula,

relevantes à explicação da teoria em micro escala.

2.1 O Modelo Clássico de Filtração Profunda com Liberação dePartículas

O modelo de filtração profunda clássico consiste em duas equações, correspondentes ao balanço

de partículas e à cinética de captura (HERZIG et al., 1970; PAYTAKES et al., 1974). Esse modelo

contem o coeficiente de filtração profunda, que é função das concentração de partículas capturadas.

Várias formas para esse coeficiente já foram propostas na literatura, conforme o(s) mecanismo(s) de

captura em questão (HERZIG et al., 1970; KUHNEN et al., 2000), as quais, por vezes, possibilitaram

a solução analítica.

Usualmente ocorre a medida da concentração de efluente no decorrer dos testes laboratoriais para

fluxo de suspensões em testemunhos, o que permite calcular o coeficiente de filtração (SOO et al.,

1986; PANG; SHARMA, 1997; BOLSTER et al., 1998; FOPPEN et al., 2005; FOPPEN; SCHIJVEN,

6

2006). De fato, a solução inversa ao problema da filtração profunda em meios porosos é bem posta

(ALVAREZ et al., 2006), ou seja, existe, é única e estável. A solução direta, desde que o coeficiente

de filtração seja conhecido, permite prever não só o perfil de retenção de partículas como também

a curva de concentração do efluente (curva de breaktrhough). O tempo de breaktrough é aquele a

partir do qual é possível detectar a mudança na concentração do efluente ocasionada pela injeção da

suspensão.

Recentemente, porém, foram publicados uma série de trabalhos que mostram que os perfis de

retenção previstos a partir do coeficiente de filtração, obtido pela inversão da curva de efluente,

não correspondem ao perfis encontrados experimentalmente (BRADFORD et al., 2002, 2003, 2004;

BRADFORD; BETTAHAR, 2006; TUFENKJI; ELIMELECH, 2004b; AL-ABDUWANI et al., 2004;

AL-ABDUWANI, 2005). Em outras palavras, esses autores reportam que é impossível ajustar simul-

taneamente o perfil de partículas retidas e a concentração delas no efluente através de optimização

sobre o coeficiente de filtração. Essa contradição torna necessária a análise dos fundamentos do mod-

elo clássico para transporte de suspensões em meios porosos, incluindo a derivação em micro escala

das equações, bem como a formulação de possíveis generalizações.

Outra motivação para a generalização do modelo clássico é a diferença entre os tempos de

breakthrough previstos e observados experimentalmente. Em testes laboratoriais de filtração pro-

funda esses tempos exibem valores entre 0,4 a 100 volumes porosos, enquanto a teoria clássica prevê

que este tempo é igual à um volume poroso apenas (HARVEY; GARABEDIAN, 1991; ROQUE et

al., 1995; KRETZSCHMAR et al., 1997; CHAUVETEAU et al., 1998; CAMESANO et al., 1999;

VEERAPEN et al., 2001; MASSEI et al., 2002). A introdução da dependência da porosidade com a

retenção de partículas no modelo clássico consegue explicar tempos de breakthrough menores que a

unidade (CORAPCIOGLU; CHOI, 1996).

A taxa de retenção, não só para captura por exclusão pelo tamanho, como também na captura

por forças elétricas, depende dos tamanhos de poros e de partículas, que variam devido a natureza

das rochas e das suspensões. Como o modelo clássico lida apenas com médias, ele não considera

características microscópicas como as distribuições de tamanho, energia superficial e carga. Alguns

modelos considerando esses problemas serão abordados brevemente na secção 2.2.

Esta seção se dedicará a descrição do modelo clássico, com especial atenção à filtração profunda.

Seguir-se-á neste trabalho a abordagem presente em Herzig et al. (1970), deduzindo-se a equação

do balanço de massa e expondo a equação de cinética de captura, com uma breve explanação sobre

os mecanismos de captura em seguida. Logo após, seram apresentadas as equações de crescimento

do diferencial de pressão. Em seguida, tratar-se-á da adimensionalização do sistema, dos modelos

avançados de filtração profunda e por fim das forças atuantes sobre a partícula depositada em micro

escala.

7

2.1.1 Balanço de Massa

Para o modelo clássico de filtração profunda assumem-se as seguintes premissas: não há agre-

gação de partículas, a difusão é pequena, a água e as partículas são incompressíveis e seus volumes

são aditivos na suspensão (lei de Amagat). Além disso, considera-se o fluido newtoniano com vis-

cosidade independente da concentração da suspensão. Todas essas hipóteses são razoáveis caso a

concentração de partículas em suspensão seja pequena.

Seja um meio poroso com área de seção transversal variável A(x), saturado por apenas um fluido,

no qual se deseja injetar uma determinada suspensão, que será filtrada. A conservação da massa

postula que:

{Tx. de Acumulação de Partículas}= {Fluxo de Partículas Entrando}−{Fluxo de Partículas Saindo} .

Seja c a concentração da suspensão (em vol./vol. da suspensão), σ a concentração de partículas

retidas (em vol./vol. do meio poroso), U a velocidade (Darcy) da suspensão e φ a porosidade (ver

Figura 1). Então escreve-se a equação do balanço de volumes unidimensional como:

∆xA(x)[(φc+σ)t+∆t− (φc+σ)t

]∆t

= (AUc)x− (AUc)x+4x .

Rearranjando obtém-se:

A(x)[(φc+σ)t+∆t− (φc+σ)t

]∆t

=(AUc)x− (AUc)x+4x

∆x.

Fazendo-se o limite ∆x,4t→ 0:

A(x)∂

∂ t(φc+σ)+

∂

∂x[A(x)Uc] = 0. (1)

É possível estender a equação (1) para o caso tridimensional. Esta fórmula vem dada por (BAREN-

BLATT et al., 1990):

Figura 1: Elemento de Meio Poroso (adaptado de Herzig et al. (1970)).

8

∂

∂ t(φc+σ)+∇ · (Uc) = 0. (2)

Usualmente as condições iniciais e de contorno são, para a geometria linear unidimensional:

c(x,0) = σ (x,0) = 0, (3)

c(0, t) = c0, (4)

nas quais c0 é a concentração da suspensão injetada, assumida constante. Já para o caso radial unidi-

mensional:

c(r,0) = σ (r,0) = 0 (5)

c(rw, t) = c0 (6)

2.1.2 Cinética de Captura

A Eq. (1) tem duas funções dependentes, fazendo-se necessário então incluir mais uma equação.

Neste contexto, introduz-se a cinética de captura:

∂σ

∂ t= λ (σ)Uc− kdet (U)σ , (7)

em que λ (σ) é a função de filtração e kdet (U) é a taxa de liberação de partículas. Caso λ (σ)

seja constante, λ é denominado de coeficiente de filtração. Em Herzig et al. (1970) são discutidas

diversas formas para λ (σ). Neste sentido, fica demonstrado nesta referência que caso a deposição

seja pequena, λ (σ) pode ser considerado constante. A justificativa básica disto é que se quando

há pouca deposição, o número de "vagas" para partículas se depositarem permanece essencialmente

o mesmo. Como λ depende essencialmente da quantidade de vagas disponíveis para deposição de

partículas, ele deve ser constante neste caso. Outra forma comum para λ (σ) é a linear, que não será

adotada neste trabalho porque o tempo de estabilização da permeabilidade não é finito, em contraste

com os testes de laboratório a serem apresentados mais a frente.

Muitas vezes se faz necessário ou é mais natural conceber a filtragem das suspensões em meios

porosos como dependentes de 2 ou mais mecanismos de captura. Sejam n mecanismos de captura,

cada qual com uma taxa de filtração λi (~σ), que é responsável pela deposição da concentração σi,

onde ~σ = (σ1,σ2, . . . ,σn). Fazendo σ = ∑ni=1 σi em (1) e escrevendo:

∂σi

∂ t= λi (σ)Uc, (8)

9

constrói-se o modelo clássico de filtração profunda à mais de um mecanismo. No apêndice A será

detalhado o procedimento de solução para o sistema a dois mecanismos com liberação de partículas

segundo a teoria proposta, aglutinando os dois mecanismos de captura em um só, conforme descrito

em Guedes et al. (2006).

O modelo proposto neste trabalho não se utilizará do parâmetro kdet na equação (7), uma vez

que considera a liberação de partículas instantânea, ao passo que qualquer valor finito positivo deste

parâmetro faz com que a quantidade máxima retida seja alcançada apenas assintóticamente, ou seja,

requer tempo infinito para estabilização da concentração retida. Então, sempre que essa equação for

referida, entenda-se que kdet = 0.

2.1.3 Mecanismos de Captura

Segundo Herzig et al. (1970) existem essencialmente cinco tipos de mecanismos de captura

(chamados por essa referência de “processos de captura”):

1. Sedimentação: Caso as partículas apresentem massa específica distinta a do fluido que as man-

tem em suspensão, elas se tornam suscetíveis à ação da gravidade e não possuem mais a mesma

velocidade que o fluido, podendo ir de encontro ao meio filtrante por sedimentação, onde são

capturadas.

2. Inércia: Devido ao peso aparente do fluido, as partículas podem se desviar das linhas de fluxo

quando nelas ocorre uma mudança abrupta. Não podendo acompanhar essa linhas, a partícula

pode entrar em contato com os grãos do meio poroso.

3. Efeitos Hidrodinâmicos: Devido a não uniformidade do campo de cisalhamento e a não esfe-

ricidade das partículas, alguns efeitos hidrodinâmicos podem ocorrer, ocasionando migração

lateral das partículas suspensas, com possível encontro delas com os grãos.

4. Interceptação Direta: Ainda que a partícula e o fluido possuam a mesma massa específica,

as partículas, devido ao seu tamanho, podem não serem capazes de acompanhar as menores

tortuosidades das linhas de fluxo do fluido carreador e por isso elas colidem com as paredes das

áreas convergentes dos poros. Este mecanismo é algumas vezes denominado de exclusão pelo

tamanho.

5. Difusão por Movimento Browniano: As partículas se difundem e pode acessar áreas que não

são usualmente irrigadas pela suspensão, e ficam aí retidas. Para que o movimento Browniano

seja importante, Kim e Karrila (2005) afirmam que o tamanho da partícula deve ser menor que

0,1 micrometros.

10

2.1.4 Crescimento do Diferencial de Pressão

Para modelar o aumento do diferencial de pressão no meio poroso devido ao acumulo de partícu-

las, é utilizada a lei de Darcy modificada (BARENBLATT et al., 1990):

U =−k0k (σ)

µ∇P, (9)

em que k0 é a permeabilidade inicial, µ é a viscosidade do fluido, P é a pressão e k (σ) é uma relação

que prevê como se dá a redução percentual da permeabilidade com a deposição. Como todos os casos

tratados nesta dissertação são unidimensionais, reescreve-se (9) como:

U =−k0k (σ)

µ

∂P∂x

. (10)

A função k (σ) é escrita como:

k (σ) =1

1+βσ. (11)

na qual o coeficiente β é o coeficiente de dano à formação.

Neste trabalho, β será considerado constante, mas em geral depende das características das

partículas depositadas e em suspensão, bem como das características do meio poroso. Em caso de

atuação de dois mecanismos, utiliza-se:

k (σ1,σ2) =1

1+β1σ1 +β2σ2. (12)

A equação (10) pode ser utilizada tanto no caso linear quanto no caso radial. Considere um teste

laboratorial de injeção de suspensão em testemunho (um exemplo de fluxo em meio poroso linear).

Seja ∆P o diferencial de pressão necessário para sustentar uma determinada vazão U , possivelmente

variável no tempo. Para calcular ∆PU manipula-se (10), do que resulta:

∆P(t)U (t)

=µ

k0

L∫0

dxk (σ (x, t))

, (13)

na qual L é o comprimento do testemunho. A importância da equação (13) é o relacionamento entre

as quantidades controláveis do teste (velocidade de Darcy e pressão) com os parâmetros da filtração

(resumidos pela função k (σ (x, t))).

Considere agora o caso axi-simétrico radial. Neste caso, a velocidade de Darcy U do escoamento

depende da posição radial r:

11

U (r) =q

2πrh, (14)

onde q é a vazão de injeção e h é a espessura. Assim, a equação (10) deve ser alterada para:

q2πrh

=−k0k (σ)

µ

∂P∂ r

, (15)

cuja manipulação resulta em:

∆P =qµ

2πhk0

re∫rw

drrk (σ (r, t))

, (16)

em que re é o raio externo, usualmente tomado como a metade da distância entre o poço injetor e

produtor.

2.1.5 Adimensionalização do Sistema (1), (7), (10)

Para adimensionalizar o sistema clássico de filtração profunda em geometria linear, definem-se

as seguintes variáveis adimensionais:

cD =cc0, σD =

σ

φc0, λD = λL, xD =

xL, tD =

∫ t0 U (η)dη

φL. (17)

A variável tD nada mais é do que o número de volumes porosos injetados. Com a substituição de (17)

em (1) e (7) obtêm-se para a geometria linear:

∂

∂ tD(cD +σD)+

∂cD

∂xD= 0, (18)

∂σD

∂ tD= λDcD. (19)

As condições (3) e (4) se tornam respectivamente:

cD (xD,0) = σD (xD,0) = 0 (20)

cD (0, tD) = 1 (21)

Para mensurar a evolução da resistência à injeção, define-se impedância J (tD), que está rela-

cionada tanto com a pressão quanto com a velocidade. Para a geometria linear:

12

J (tD) =∆P(tD)U (0)∆P(0)U (tD)

, (22)

que pode ser manipulada usando (13), resultando em:

J (tD) = 1+βφc0

1∫0

σDdxD. (23)

Embora não seja adimensionalisado neste trabalho o sistema de filtração profunda para a geome-

tria radial, utilizar-se-á a impedância nele também. Caso a vazão de injeção seja mantida constante, a

definição dessa quantidade em coordenadas radiais é:

J (tD) =∆P(t)∆P(0)

, (24)

e combinando esta equação com (16), conclui-se que:

J (tD) = 1+β

lnre/rw

re∫rw

σ (r, tD)dr. (25)

Entretanto, é preciso observar que a definição do tempo adimensional difere daquela do sistema linear

e vem dada por:

tD =

∫ t0 q(η)dη

φπ (r2e − r2

w)h≈∫ t

0 q(η)dη

φπr2eh

, (26)

onde re é o raio externo.

2.2 Outros Modelos de Filtração Profunda

Muitos estudos têm sido propostos para considerar variações locais na eficiência da captura

(ELIMELECH et al., 1995; KHILAR; FOGLER, 1998). A função de distribuição de densidade de

carga superficial foi ajustada conforme os dados experimentais (BAYGENTS et al., 1998; SIMONI

et al., 1998; BOLSTER et al., 1998, 1999). Em particular, a distribuição bimodal ajusta com precisão

considerável diversos dados simultâneos de perfil de retenção e concentração de efluente nos quais o

modelo clássico falha (TUFENKJI et al., 2003; TUFENKJI; ELIMELECH, 2005).

Sharma e Yortsos (1987a, 1987b, 1987c) deduziram as equações básicas de balanço de particulado

no transporte de suspensões em meios porosos. O modelo leva em consideração a variação nas dis-

tribuições de poro e partículas nos diversos mecanismos de captura. É assumido que o todo o espaço

poroso é acessível para partículas e que a população de partículas se move com a velocidade média

do fluido que a carreia. A introdução dos fatores de redução ao fluxo e à acessibilidade na equação do

13

balanço de populações é apresentada em Santos e Bedrikovetsky (2006). Em Bedrikovetsky (2008) é

incluída a estocagem de partículas em poros não-conectados e a variação na porosidade.

2.3 Forças em Micro Escala

Nesta secção, tratar-se-á das forças em micro escala cujo balanço de momentos determina a con-

centração máxima retida conforme a teoria aqui proposta (ver detalhes no próximo capítulo). Basica-

mente, pode-se dividir as forças atuantes sobre a partícula conforme suas naturezas em: eletrostática,

de arraste, de elevação e gravitacional. Nas próximas subsecções elas serão explanadas em detalhe.

2.3.1 Força Eletrostática

Tendo em vista a teoria DLVO (DERJAGIN; LANDAU, 1941; GREGORY, 1981; ELIMELECH

et al., 1995; KHILAR; FOGLER, 1998; ISRAELACHVILI, 2006), referente as interações inter-

partícula, é possível descrever a força eletrostática como resultantes da soma de três potenciais de

interação distintos: o potencial de London - van der Waals, o potencial de dupla camada e o potencial

de Born.

O potencial de London - van der Waals gera uma força atrativa:

VLVA =A132

6

[2(1+Z)Z (2+Z)

+ ln(

Z2+Z

)], (27)

em que Z = hrs

é a razão entre a separação entre a rocha e a partícula e A132 é a constante de Hammaker.

O potencial de Dupla Camada, cuja força é repulsiva ou atrativa a depender se as cargas da

partícula e do meio poroso são iguais ou distintas respectivamente:

VDLR =ε0Drs

4

[2ψ01ψ02 ln

(1+ exp(−κh)1− exp(−κh)

)−(ψ

201 +ψ

202)

ln(1− exp(−2κh))], (28)

em que ε0 é a permissividade do espaço livre (vácuo), D é a constante dielétrica, ψ01 e ψ02 são os

potenciais da partícula e do meio poroso, e κ é o inverso do comprimento de Debye, calculado pela

seguinte expressão:

κ =

√e2 ∑viz2

iε0DkBT

, (29)

em que kB é a constante de Boltzmann, vi é a concentração média de íons (em íons por unidade de

volume), zi é a valência do i-ésimo íon e e é a carga do elétron. Para soluções aquosas sob temperatura

normal, (29) pode ser escrita como (ELIMELECH et al., 1995):

14

κ = 2,32.109√(

∑Cmiz2i), (30)

em que Cmi é a concentração molar do composto iônico (em mols/L) e κ está em m-1.

O potencial de Born, cuja força é repulsiva e que vem dado por:

VBR =A132

7560

(σLJ

rs

)6[

8+Z

(2+Z)7

](31)

no qual σLJ é o diâmetro atômico de colisão no potencial de Lennard-Jones (LANDAU; LIFSHITZ,

1980).

A força elétrica, que tem valor líquido negativo e por isso é atrativa, vem dada por:

Fe =−∂VT

∂h, (32)

onde VT =VLVA +VDLR +VBR. A Figura 2 mostra os gráficos da energia potencial total e força versus

a separação entre partícula e superfície de grão.

O mínimo de energia em h= h1 corresponde ao ponto onde a força é nula (derivada primeira igual

à zero). O ponto de inflexão h = h2 corresponde ao zero da derivada segunda, que é a força mínima.

Para construir essa Figura, foram usados os seguintes valores para os parâmetros: A132 = 2,00zJ,

ψ01 =−30mV, ψ02 =−50mV, ε0 = 8,854.10−12C2J−1m−1, D = 78,0 (constante dielétrica da água

segundo Khilar e Fogler (1998)), diâmetro atômico de colisão no potencial de Lennard-Jones é σLJ =

0,5nm (KHILAR; FOGLER, 1998); κ foi calculado segundo a fórmula (30) com Cm = 0,513M e

z = 1, que é a valência do cloreto e do sódio no Cloreto de Sódio NaCl.

2.3.2 Força de Arraste

Uma expressão geral para a força de arraste atuante sobre uma partícula esférica em fluxo do tipo

Hele-Shaw entre duas placas paralelas (LANDAU; LIFSHITZ, 1987) é dada por:

Fc f =ωπµr2

sU ′

H, (33)

na qual ω é um coeficiente empírico no intervalo [10;60], µ é a viscosidade, rs é o raio da partícula,

U ′ é a velocidade média numa dada secção transversal e H é a altura do canal (O’NEIL, 1968; ALT-

MANN; RIPPERGER, 1997; AL-ABDUWANI et al., 2005).

15

(a)

(b)

Figura 2: Gráficos de potencial total (a) e força elétrica (b)

16

2.3.3 Força de Elevação (“lifting”)

A forma geral para a força de elevação sobre uma partícula esférica em fluxo do tipo Hele-Shaw

entre duas placas paralelas é dada por:

Fl = χr3s

√ρµU ′3

H3 , (34)

na qual H é a altura do canal e χ é o coeficiente de elevação, dado por 89,5 segundo Kang et al. (2004),

enquanto Altmann e Ripperger (1997) estima um valor de 1190. Expressões similares à (34) foram

usadas por Akhatov et al. (2008). Todos os autores citados nesta secção derivaram suas equações

baseados em Saffmann (1965, 1968).

2.3.4 Força Gravitacional

Nesta seção considerar-se-á não a força gravitacional bruta apenas, mas também o empuxo gerado

pela imersão completa da partícula no fluido. Assim, concatenamos ambas as forças como:

Fg =4πr3

s3

g∆ρ, (35)

onde g é a aceleração da gravidade e ∆ρ = ρpart−ρ , ou seja, é a diferença entre as massas específicas

da partícula e do fluido.

Todavia, a gravidade tende a liberar partículas da porção superior do poro e depositar partículas

na porção inferior dele. Aliando-se a isso a pouca relevância dessa força comparada com aquela de

origem eletrostática para os tamanhos de partícula deste trabalho (ver Capítulo 5, Tabela 3), torna-se

possível desprezar a força gravitacional nas estimativas de concentração máxima retida em função da

velocidade.

17

3. METODOLOGIA

Neste capítulo formular-se-á o problema matemático subjacente à estabilização da injetividade

devido à erosão do reboco interno, construindo a Teoria Alternativa de Perda de Injetividade. Primeiro,

serão reescritas as equações de captura de partículas em meios porosos levando-se em conta o con-

ceito de capacidade de estocagem. Em seguida, serão descritas as bases microscópicas do conceito

de capacidade de estocagem, de forma que ele possa ser colocado validado, bem como tornar pos-

sível a previsão das curvas de concentração máxima retida. Finalmente, será feita uma breve análise

dimensional da concentração máxima retida.

3.1 Revisão das Equações de Filtração Profunda

Nesta secção serão reescritas as equações de filtração profunda em meios porosos, passando

primeiro pela enunciação do conceito de capacidade de estocagem e depois incorporando este princí-

pio nas equações deduzidas anteriormente.

3.1.1 O Conceito de Capacidade de Estocagem

Considere as partículas aderidas à parede do poro confome mostra a Figura 3.

Figura 3: Esquema da "última" partícula do reboco interno

Nessa figura pode-se notar uma "última" partícula depositada, sobre a qual atuam a forças elétrica

18

Fe, a força de arraste Fc f , a força de lifting Fl e a força da gravidade Fg, que pode ser desprezada

frente às outras três forças. Conforme a deposição se processa, o espaço disponível ao fluxo do poro

diminui e, caso a velocidade do fluido seja constante, as forças de arraste Fc f e de lifting Fl aumentam

até o momento em que vencem as forças elétricas e a partícula é simplesmente levada com o fluxo, o

que ocasiona o fim da captura e o estabelecimento da concentração máxima de retenção.

Seguindo a abordagem unidimensional apresentada em Jiao e Sharma (1994) e em Freitas e

Sharma (2001), assume-se que o equilíbrio da partícula na superfície do reboco interno é determi-

nado pela igualdade de torques para as forças normal e de arraste:

Fc f lc f > Fnln, (36)

onde Fn = Fe−Fl é a força normal resultante; lcr e ln são os braços das forças de arraste e normal

respectivamente. O critério de torque (36) tem a mesma forma que o critério de balanço de forças:

Fc f > κCFn, (37)

onde κC é o coeficiente de atrito seco de Coulomb.

Resta agora discutir qual seria o valor mais adequado para a força eletrostática no equilíbrio da

partícula. Considere o último momento no qual a partícula ainda esta alojada na superfície do reboco

interno logo antes de ser liberada pela força de arraste. A Figura 3 mostra como a partícula gira em

torno do ponto tangente com a partícula vizinha neste instante, onde a partícula acaba de deixar o

estado de equilíbrio. Durante este movimento, a distância entre a partícula e o reboco aumenta de um

valor mínimo até um valor onde a interação é desprezível, o que corresponde à liberação da partícula.

Se em alguma distância intermediária a desigualdade (36) não é satisfeita, a partícula para naquela

posição e não deixa o reboco (não é liberada), o que implica que o torque da força de arraste deve

vencer o maior torque possível da força normal, que é atingido quando a força eletrostática é máxima,

ou seja, deve-se utilizar Fe (h2) (ver Figura 2.b) como valor da força eletrostática. Os braços das

forças de arraste e normal (lc f e ln) tem a mesma ordem de magnitude, tendo sua razão igual à√

3

para o caso de partículas esféricas de mesmo tamanho (ver Figura 3).

Define-se matematicamente o conceito de capacidade de estocagem da seguinte forma: existe

uma função

σcr := σcr(Fe,Fl,Fc f , ln, lc f ,φ ,φc

), (38)

denominada função capacidade de estocagem tal que as partículas são filtradas à uma taxa constante

pelo meio poroso enquanto σ ≤ σcr (ε). Para os casos onde existem dois (ou mais) mecanismos de

captura, devera se considerar que apenas um deles está sujeito a erosão, o que significa aplicar a

19

restrição na concentração retida em apenas um dos mecanismos (neste trabalho a restrição é aplicada

no primeiro mecanismo). Caso a suspensão a ser injetada mantenha suas características constantes

(salinidade, tipo e tamanho das partículas entre outras), a concentração máxima retida será função da

velocidade da injeção apenas, como mostra a Figura 4.

Figura 4: Concentração máxima retida em função da velocidade

Seja o seguinte experimento para determinação da função máxima de retenção σcr (U). Primeiro,

injeta-se uma suspensão com concentraçao em suspenção c0 à velocidade U1 num meio poroso com

espessura desprezível frente ao comprimento médio de captura das partículas 1/λ , para que o perfil de

retenção nesse meio poroso possa ser considerado uniforme, como mostra a Figura 5.a.

(a)

(b)

Figura 5: Históricos esquemáticos de (a) velocidade e (b) concentração no efluente em experimentohipotético

20

A concentração no efluente c(L, t), que é mostrada na Figura 5.b, é menor que aquela injetada

devido à retenção de partículas, permanecendo constante até o momento t1, quando a concentração

de partículas depositadas atinge o valor máximo σcr (ε1). A partir desse instante, a concentração no

efluente aumenta abruptamente até a c0. Logo após, em t2, injeta-se água sem partículas com a mesma

velocidade U1, o que causa o deslocamento da suspensão pela água limpa, e a concentração reduz à

zero durante um curto período após a troca da suspensão pela água “limpa” (Figura 5.b).

Depois aumenta-se a velocidade de injeção para U2 no momento t3, o que resulta em liberação

instantânea de algumas partículas capturadas que deixam o testemunho e criam um banco de concen-

tração no efluente. O volume do banco ∆σ1 (área sob a curva de concentração, Figura 4) corresponde

ao número de partículas liberadas da rocha durante o incremento de velocidade de U1 para U2 (Figura

5).

Aumentos subseqüentes da velocidade de fluxo da água sem partículas resultam em diversos

bancos de partículas no efluente; ∆σk é o número de partículas liberado na mudança da velocidade

de Uk+1 para Uk+2 (Figura 5.a). A medida dos bancos de concentração no efluente permite obter a

função de retenção máxima (Figura 4).

3.1.2 Equações de Filtração Profunda Modificadas

De acordo com o exposto anteriormente, deve-se adequar o sistema de equações de filtração

profunda, modificando sua cinética de captura para:

∂σ

∂ t= λcU, σ < σcr (U) , (39)

σ = σcr (U) caso contrário,

e para o caso com 2 mecanismos, reescrevemos a cinética de captura de σ1 como:

∂σ1

∂ t= λ1cU, σ1 < σcr (U) , (40)

σ1 = σcr (U) , caso contrário,

e o balanço de massa (1) se mantém. Para verificar se o sistema proposto acima de fato é aplicável,

ajustar-se-ão dados de laboratório e de campo.

Para a geometria linear, escreve-se (39) e (40) em termos adimensionais respectivamente como:

21

∂σD

∂ tD= λDcD σD < σcrD (U) , (41)

σD = σcrD (U) caso contrário,

∂σ1D

∂ tD= λDcD σ1D < σcrD (U) , (42)

σ1D = σcrD (U) caso contrário.

3.2 Bases Microscópicas do Conceito de Capacidade de Estocagem

Para testar a hipótese da capacidade de estocagem, concebe-se que o meio poroso seja com-

posto por n capilares quadrados paralelos de secção transversal de lado H (Figura 6). Além disso,

considera-se que as partículas depositadas e em suspensão são todas esféricas e com o mesmo raio.

Assumindo que o fluxo em cada um desses capilares segue a equação de Poiseülle, a porosidade e a

permeabilidade inicial desse meio poroso são escritas como:

φ = nH2,

k0 =nH4

8π

. (43)

Rearranjando e colocando H e n em evidência temos:

H =√

8πk0φ

,

n = φ 2

8πk0

. (44)

Essa decomposição do meio poroso é importante pois torna possível encontrar os parâmetros necessários

à avaliação das forças.

a) b)

Figura 6: Idealização do meio poroso como um feixe de capilares em perfil (a) e em secção transversal(b)

Efetuamos o cálculo teórico da concentração máxima retida a partir do balanço de torques sobre a

última partícula depositada (seguindo o mesmo critério que (CIVAN, 2007)): as forças Fl e Fe atuam

sobre o braço ln enquanto a força Fc f atua sobre o braço lc f . Em linguagem matemática:

22

Fn =√

3Fc f . (45)

Substituindo as expressões para as forças (equações (32), (33) e (34)) obtém-se:

Fe−χr3s

√√√√ ρµU3

(H−2hc)3 −√

3ωπµr2sU

H−2hc= 0, (46)

em que hc é a altura do reboco interno, U é a velocidade intersticial e Fe é o módulo da força elétrica

(atrativa). A velocidade intersticial se relaciona com a velocidade Darcy por:

U =Uφ. (47)

Resta agora escrever hc em função de σcr, o que pode ser feito através de:

σcr =Vdepositado

Vrocha=

4φ (1−φc)

H2

(Hhc−h2

c), (48)

ou alternativamente:

hc =H(

1−√

1− σcrφ(1−φc)

)2

(49)

Substituindo (47),(49) em (46) obtém-se a expressão final implícita que determina a função ca-

pacidade de estocagem através de:

Fe−χr3s

√√√√√ ρµU3

φ 3(

H√

1− σcrφ(1−φc)

)3 −√

3ωπµr2sU

φH√

1− σcrφ(1−φc)

= 0. (50)

Para verificação serão usados os ajustes realizados nos dados de laboratório para se obter pontos

(U,σcr) da curva da função capacidade de estocagem. Esses pontos servirão como base para ajustar

a curva teórica descrita pela equação (50), na qual os parâmetros devem estar no intervalo descrito

na literatura. Para a obtenção dessas curvas, utilizaremos o método de Newton-Raphson na equação

(50), visando ajustar os dados laboratoriais disponíveis. O método de Newton-Raphson utilizado está

descrito em Ruggiero e Lopes (1996).

3.3 Análise Dimensional de σcr

A Figura 3 mostra que existem cinco parâmetros que influênciam a concentração máxima retida

σcr: as forças elétrica Fe, de arraste Fc f e de elevação Fl , juntamente com os braços ln e lc f , referentes

23

as forças normais e de arraste respectivamente. Além deles, devem ser incluídas as porosidades da

formação e do reboco interno, já que o limite superior para a concentração máxima retida é φ (1−φc).

Para proceder com a análise dimensional segundo Barenblatt (1996), defini-se então, como a classe

de sistema de unidades, a classe MLT e usamos a força elétrica Fe e o braço ln como as grandezas

linearmente independentes quanto a função dimensão. A partir da equação (38), é possível definir

então 3 números adimensionais:

ΠFl =Fl

Fe, (51)

ΠFc f =Fc f

Fe, (52)

Πlc f =lc f

ln, (53)

que juntamente com φ e φc são os parâmetros adimensionais que regem a concentração máxima

retida. Embora σcr seja adimensional, é fácil perceber que a quantidade máxima retida depende das

porosidades da rocha e do reboco, fato que torna a comparação direta de σcr entre dois meios porosos

muito difícil. Por isso, se faz necessário representar a concentração máxima retida por outro número

adimensional:

Σcr =σcr

φ (1−φc), (54)

que é diretamente comparável entre meios porosos distintos, variando entre zero e um. Com essas

definições, é possível escrever o teorema dos Π’s de Buckingham da seguinte forma:

Σcr = Φ

(ΠFc f ,ΠFl ,Πlc f ,φ ,φc

), (55)

em que Φ é uma relação (função) qualquer relacionando os parâmetros adimensionais.

Nesta dissertação, já foi considerada uma forma específica de relação para σcr (45). Para essa

relação específicamente, escreve-se:

ε = ΠFc f , (56)

Πlc f =√

3, (57)

ΠFl = 1− ε√

3, (58)

24

e o enunciado do teorema dos Π’s de Buckingham se torna:

Σcr = Φ

(ε,1− ε

√3,√

3)= Φε (ε) , (59)

em que Φε é uma função que só depende de ε . Nesse caso especial, por haver apenas um parâmetro

adimensional independente, chamar-se-á ε de número de erosão.

25

4. DESENVOLVIMENTO

Este capítulo objetiva resolver as equações de filtração profunda modificadas no Capítulo 3.

Primeiro resolver-se-á o sistema de geometria linear com um e dois mecanismos, seguindo-se com a

geometria radial de um mecanismo.

4.1 Geometria Linear

4.1.1 Um Mecanismo

Para um mecanismo, temos o sistema enunciado no capítulo 3 (equações (1), (39), (10), (11),

(3), (4) ). Esse sistema admite solução semi-analítica sob as premissas de porosidade, coeficientes de

filtração e dano à formação constantes. As condições sobre as quais essas premissas são verificadas

encontram-se na literatura, conforme descrito a seguir.

Seja o caso de captura elétrica com alta concentração de partículas retidas, formando poli-camadas

na parede dos poros. A acumulação de partículas retidas resulta no decremento da secção transver-

sal porosa e num aumento da velocidade intersticial. Por um lado, o aumento da velocidade causa

o aumento do fluxo carreador de partículas da solução para as áreas proximais do grão onde a in-

terceptação de partículas e captura pela dupla camada elétrica ocorre (ELIMELECH, 1994), o que

acarreta um aumento do coeficiente de filtração. Por outro lado, o aumento de velocidade causa o de-

créscimo do tempo de permanência e conseqüente redução da probabilidade e colisão entre partícula

e grão (TUFENKJI; ELIMELECH, 2004a), o que diminui o coeficiente de filtração. A existência de

dois fatores competitivos pode causar uma pequena variação no coeficiente de filtração em função

da velocidade (JOHNSON et al., 2007). No caso da deposição da partículas por camadas, o efeito

da espessura do reboco interno (concentração retida) no coeficiente de filtração também é pequeno.

A modelagem numérica de Dabrowski (1988) mostra que a diminuição da porosidade não acarreta

nenhuma modificação significativa na filtração profunda de suspensões.

Enfim, para solução do sistema (18), (41), (10), (3) e (4), recorre-se ao método das características.

A Figura 7.a mostra um esboço das principais características do sistema, juntamente com a frente de

erosão xcrD .

26

(a)

(b)

(c)

Figura 7: Transporte de suspensão em meio poroso com um mecanismo de captura e liberação departículas retidas: (a) movimento das frentes de concentração no plano (xD, tD), com os pontos P =(xD, tD) e P’ = (xcrD (ηD) ,ηD); (b) perfil de retenção em quatro momentos e (c) perfil de concentraçãoda suspensão nesses quatro momentos

Da análise deste gráfico, é possível depreender que a única parte da solução que não esta ime-

diatamente disponível é aquela pertencente a região II. Isso acontece porque as caracteristícas que

passam por um ponto genérico P = (xD, tD) não mais interceptam o eixo do tempo adimensional, mas

sim a frente de erosão xcrD no ponto P’ = (xcrD (ηD) ,ηD), no qual ηD é o valor de tD no encontro entre

essas duas linhas. Por isso, deve-se escrever:

27

cD = σD = 0 Região 0 (60)

cD = exp(−λDxD) Região I (61)

σD = λD (tD− xD)exp(−λDxD) Região I (62)

cD = exp [−λD (xD− xcrD (ηD))] Região II (63)

σD = σcrD exp [−λD (xD− xcrD (ηD))] Região II (64)

cD = 1 Região III (65)

σD = σcrD Região III (66)

As linhas limítrofes entre as regiões 0-I e I-II vem dadas respectivamente por:

xD = tD, (67)

tD = xD +σcrD

λD. (68)

O perfil destas soluções é apresentado nas Figuras 7.b e 7.c.

O próximo passo é encontrar xcrD (tD). Para isso, deve-se escrever a condição de continuidade no

fluxo das partículas nessa frente (LANDAU; LIFSHITZ, 1987):

c+D (1−D)−σ+D D = c−D (1−D)−σcrDD, (69)

em que D é a velocidade da descontinuidade (frente de erosão), c+D e c−D são as concentrações (reduzi-

das) em suspenção referenciadas à direita e à esquerda da descontinuidade respectivamente, e σ+D é a

concentração reduzida retida à direita da frente. A partir de (39) conclui-se que σD é contínua, logo

σ+D = σcrD , e:

(c+D− c−D

)(1−D) = 0, (70)

o que por sua vez implica que ou a velocidade da frente é 1 ou a concentração em suspensão é

contínua. Como a velocidade é sempre menor que 1, c+D = c−D . Com isso em mente, já se pode iniciar

o procedimento para a obtenção de xcrD (tD). Em primeiro lugar:

28

σD (xcrD (tD) , tD) = σcrD. (71)

Diferenciando (71) com respeito a tD:

∂σD

∂ tD

∣∣∣∣(xcrD(tD),tD)

+dxcrD

dtD

∂σD

∂xD

∣∣∣∣(xcrD(tD),tD)

= 0. (72)

As equações (39), (70) e (A.2) (ver apêndice A) permitem escrever:

∂σD

∂ tD

∣∣∣∣(xcrD(tD),tD)

= λD, (73)

∂σD

∂xD

∣∣∣∣(xcrD(tD),tD)

=−λD (1+σcrD) . (74)

Usando (73) e (74) em (72) conclui-se que:

dxcrD

dtD=

11+σcrD

. (75)

Para resolver (75), basta incluir como condição inicial o tempo a partir do qual a erosão começa (na

posição xD = 0):

xcrD

(σcrD

λD

)= 0, (76)

o que resulta em:

xcrD (tD) =

(tD−

σcrDλD

)1+σcrD

. (77)

Por fim, deve-se encontrar o valor do intercepto ηD na região II. Para tal, utiliza-se a equação da

característica P’P:

xD− xcrD (ηD) = tD−ηD, (78)

substituindo (77) em (78), obtêm-se:

ηD =(1+σ

−1crD

)(tD− xD)−λ

−1D . (79)

Com (79), é possível reescrever (63) e (64) como:

29

cD = exp{−λD

[− tD

σcrD

+(1+σ

−1crD

)xD +

1λD

]}Região II (80)

σD = σcrD exp{−λD

[− tD

σcrD

+(1+σ

−1crD

)xD +

1λD

]}Região II (81)

A solução apresentada nesta secção permite restabelecer a seguinte figura física e estrutura de

fluxo, durante a injeção de suspensão no meio poroso com liberação de partículas. A Figura 7.a

apresenta as frente de concentração em suspensão juntamente com as frentes frontais e distais no plano

(xD, tD). A frente de concentração da suspensão se propaga com a velocidade do fluido carreador; as

concentrações retida e em suspensão a frente dessa fronteira são nulas (zona 0). Atrás dessa frente

(zona I), a concentração em suspensão atingiu o estado estacionário e a concentração retida cresce

proporcionalmente ao tempo. No momento tD = σcrD/λD, quando a captura de partículas acaba na face

de injeção do testemunho, aparece a frente de erosão com velocidade (1+σcrD)−1. Não há retenção

de partículas atrás da frente de erosão, na zona III, e a concentração em suspensão neste momento é

aquela injetada com a água .

Ocorre filtração profunda com coeficiente de filtração constante na zona II, conforme as equações

(80) e (81).

Perfis de concentração retida e em suspensão são apresentados nas Figuras 7.b e 7.c. Quatro mo-

mentos relevantes foram selecionados para os perfis: dois momentos antes da liberação de partículas,

tD1 e tD2 , e dois posteriores ao início da erosão, tD3 e tD4 .

O perfil correspondente a tD1 nas Figuras 7.b,c foram construídos antes das partículas serem

observadas no efluente, tD1 < 1, com alguma retenção de partículas. Esse perfil aumenta no momento

tD2 , quando o breakthrough já ocorreu; a concentração de partículas em suspensão entra em estado

estacionário. Esses perfis não mudam qualitativamente até o início da erosão, tDc = σcrD/λD. Durante

o período tDc até tD3 , a zona de erosão se propaga pelo testemunho. A concentração retida se iguala

à máxima na zona de erosão, e a concentração em suspensão se torna unitária.

No momento tD4 a frente de erosão já atingiu o final do testemunho e o fluxo da suspensão se

estabiliza. O histórico de concentração em suspensão a jusante do testemunho é apresentado na

Figura 8.a. Antes do primeiro volume poroso, essa concentração é nula, mudando abruptamente

neste momento para exp(−λD), ainda constante até o fim da zona I, momento a partir do qual cresce

monotonicamente até a unidade.

30

(a)

(b)

Figura 8: Estágios de filtração profunda com erosão; histórico de concentração na saída do teste-munho (a) e curva de impedância (b)

Para o cálculo da impedância, utilizam-se as equações a seguir: caso tD < 1:

J (tD) = 1+φβc0

(1

λD+ exp(−λDtD)+ tD

), (82)

para 1≤ tD < σcrD/λD (ver simplificação em Bedrikovetsky et al. (2001)):

J (tD)≈ 1+φβc0 (1− exp(−λD)) tD. (83)

Pouco depois do início da propagação da frente de erosão no meio poroso (logo após seu início,

tD > 1+ σcrD/λD)

J (tD) = 1+φβc0

λD (σcrD +1)

[λDσcrDtD−σ

2crD

exp(

λD (tD−σcrD−1)σcrD

−1)]

, (84)

e outro caracterizado pelo fim da filtragem de partículas (passagem da frente de erosão por todo o

meio poroso, tD > 1+σcrD (1+ 1/λD)):

31

J (tD) = 1+φβc0σcrD . (85)

Uma curva esquemática de impedância com o tempo é mostrada na Figura 8.b, mostrando o cresci-

mento não-linear antes de um volume poroso injetado. Após isso o crescimento é linear até o mo-

mento em que a erosão se inicia, a partir daí o crescimento é não-linear novamente. Quando a con-

centração máxima retida é atingida em todo o testemunho, a impedância torna-se constante.

Vazões Alternadas

Um tipo de teste recorrente na literatura é aquele realizado à múltiplas vazões (como na Figura

9). Neste caso também se faz possível resolver analiticamente todas as equações. Na prática, para

resolver este sistema de equações translada-se o eixo do tempo de maneira que o tempo de mudança

da vazão seja a origem do eixo transladado, ou seja utiliza-se t∗D = tD− tD0 , onde tD0 é o momento

no qual se da a mudança de vazão. Para simplificar, será omitido daqui em diante o asterisco em

t∗D, sendo a coordenada de tempo adimensional já transladada denotada por tD simplesmente. Além

disso, seram considerados constantes e iguais para todas as vazões os coeficientes de dano à formação

β e de filtração λ . Assim, neste novo sistema de coordenadas, se a vazão for reduzida, resolvem-se

as equações (18), (41) e (10) com as seguintes condições iniciais e de contorno:

cD (xD,0) = 1, (86)

σD (xD,0) = σcr0D , (87)

cD (0, tD) = 1, (88)

que são as condições fisicamente correspondentes ao testemunho saturado pelo fluido injetado (com

partículas) e com a concentração retida igual a máxima concentração retida anterior. Para conceber

esta configuração, assume-se que a erosão já tenha atingido todo o testemunho durante a vazão an-

terior. Caso a vazão seja aumentada, ocorre a liberação de partículas ao invés de captura segundo o

conceito de capacidade de estocagem, escrevendo-se as seguintes condições iniciais e de contorno:

cD (xD,0) = 1+(σcrD−σcr0D) , (89)

σD (xD,0) = σcr0D, (90)

cD (0, tD) = 1, (91)

A solução com diminuição da vazão (condições (86)-(88)) vem dada por (ver Figura 10, Regiões

IV, V, VI e VII), com ∆σcrD = σcrD−σcr0D:

32

(a)

(b)

(c)

(d)

Figura 9: Teste típico de vazões alternadas: (a) esquema de velocidade, (b) histórico de queda depressão, (c) histórico de impedância e (d) histórico de concentração no efluente

33

σD = σcr0D +1− e−λDtD Região IV (92)

cD = e−λDtD Região IV (93)

σD = σcr0D +1+[λD (tD− xD)−1]e−λDxD Região V (94)

cD = e−λDxD Região V (95)

σD = σcr0D +1+(∆σcrD−1)exp(

λDtD−λD∆σcrDxD−1∆σcrD−1 +1

)Região VI (96)

cD = exp(

λDtD−λD∆σcrDxD−1∆σcrD−1 +1

)Região VI (97)

σD = σcrD Região VII (98)

cD = 1 Região VII (99)

e as linhas limítrofes entre as regiões IV-V, V-VI são dadas por:

xD = tD, (100)

tD = xD +∆σcrD

λD. (101)

Para encontrar a frente de erosão, utiliza-se o mesmo equacionamento utilizado para teste em

apenas uma vazão (equações (72)-(77)). A única diferença é a equação diferencial parcial em σD

(equação (A.2)), que passa a ser:

∂σD

∂ tD+

∂σD

∂xD= λD (1+σcr0D−σ) . (102)

Com esta modificação, a frente de erosão é dada por:

xcrD (tD) =tD

∆σcrD

− 1λD

. (103)

Resta agora enunciar a impedância para os diversos momentos da deposição: para tD < 1:

J (tD) = 1+φβc0

[σcr0D +1− 2

λD+ tD +

(2

λD−1+ tD

)e−λDtD

], (104)

para 1≤ tD <∆σcrD

λD:

34

(a)

(b)

(c)

Figura 10: Diagrama Caraterístico (a), perfil de retenção (b) e de suspensão (c) de partículas

35

J (tD) = 1+φβc0

[σcr0D +1+

(1+

2λD

)e−λD− 2

λD+(

1− e−λD)

tD

], (105)

para ∆σcrDλD

+1≤ tD < ∆σcrD

(1+ 1

λD

):

J (tD) = 1+φβc0

{1+σcr0D +

1−∆σcrD

λD+

(1− 1

∆σcrD

)tD

+

(∆σcrD−2+ 1/∆σcrD

λD

)exp(

λD−2∆σcrD +11−∆σcrD

)(106)

−(

∆σcrD−2+ 1/∆σcrD

λD

)exp(

1− λD (tD−∆σcrD−1)1−∆σcrD

)},

e após a erosão, tD ≥ ∆σcrD

(1+ 1

λD

):

J (tD) = 1+φβc0σcrD . (107)

Já para o caso de aumento na vazão de injeção, utilizam-se as condições (89)-(91) que, somadas

ao fato de não ocorrer mais deposição, acarretam na seguinte solução:

cD (xD, tD) = 1+(σcrD−σcr0D) , Região VIII, (108)

σD (xD, tD) = σcr0D, Região VIII, (109)

cD (xD, tD) = 1, Região IX, (110)

σD (xD, tD) = σcr0D, Região IX, (111)

Já a impedância para este caso resulta em:

J = 1+φβc0σcr0D , (112)

sendo constante.

4.1.2 Dois mecanismos

Para o sistema de dois mecanismos, considera-se que apenas o primeiro sofre erosão e resolve-se

o seguinte sistema de equações adimensionais (equações (18), (42), (19), (10), (12), (20) e (21)):

36

∂

∂ tD(cD +σ1D +σ2D)+

∂cD

∂xD= 0, (113)

∂σ1D∂ tD

= λ1DcD, σ1D < σcrD (ε) ,

σ1D = σcrD (ε) , caso contrário,(114)

∂σ2D

∂ tD= λ2DcD, (115)

U =− k0

µL(1+φc0β1σ1D +φc0β2σ2D)

∂ p∂xD

, (116)

cD (x,0) = σ1D (x,0) = σ2D (x,0) = 0, (117)

cD (0, t) = 1. (118)

Seguindo a mesma idéia aplicada ao caso com um mecanismo, constrói-se a Figura 11. Dessa

forma, torna-se possível dividir a solução em regiões e encontrar a frente de erosão, como feito no

caso precedente. Os detalhes da solução, que usa o lema da agregação conforme descrito em Guedes

et al. (2006), estão no apêndice A.

A solução é:

σ1D = σ2D = cD = 0, Região 0 (119)

cD = e−(λ1D+λ2D)xD, Região I (120)

σ1D = λ1D (tD− xD)e−(λ1D+λ2D)xD , Região I (121)

σ2D = λ2D (tD− xD)e−(λ1D+λ2D)xD , Região I (122)

ηD +(xD− xcrD (ηD))− tD = 0, Região II (123)

cD = e−λ2D xDe−λ1D(xD−xcrD(ηD)), Região II (124)

σ1D = σcrDe−(λ1D+λ2D)(xD−xcrD(ηD)), Região II (125)

σ2D =λ2D

λ1D

σcrDe−(λ1D+λ2D)(xD−xcrD(ηD)), Região II (126)

cD = e−λ2DxD, Região III (127)

σ2D =λ2D

λ1D

σcrD +λ2D (tD− tcrD (xD))e−λ2DxD, Região III (128)

σ1D = σcrD. Região III (129)

37

Figura 11: Diagrama (xD, tD) para o sistema linear com dois mecanismos

A equação (123) representa a coordenada no eixo tD da interseção da característica referente ao

ponto (xD, tD) com a curva da frente de erosão. A definição das duas linhas retas que delimitam as

regiões 0-I e I-II são respectivamente:

tD = xD, (130)

tD = xD +σcrD

λ1D

. (131)

Agora que a solução é conhecida tanto na região II quanto na III, é possível proceder como na

secção anterior e encontrar a frente de erosão. Para isso, escreve-se (ver ((175))-((177))):

∂σD

∂ tD+

∂σD

∂xD=−(λ1D +λ2D)σD, σD ≤

(1+

λ2D

λ1D

)σcrD. (132)

A quantidade(

1+λ2Dλ1D

)σcrD nada mais é do que a quantidade total de partículas depositadas quando

σ1D = σcrD . Essa quantidade é obtida a partir do lema de agregação (detalhes no apêndice A, seção

A.1).

Com tudo isso, a derivação da frente de erosão segue diretamente: o primeiro passo é escrever

dσD (x, tcr (x))dxD

=∂σD

∂ tD

dtcrD

dxD+

∂σD

∂xD. (133)

O primeiro membro da equação (133) é 0 pois:

38

dσD (xD, tcrD (xD))

dxD=

dσD

dσ1D

dσ1D (xD, tcrD (xD))

dxD= 0. (134)

A derivada do segundo membro com respeito ao tempo vem dada pela cinética de captura (114) e

(115), sendo:

(∂σD

∂ tD

)xD,tcrD(xD)

= (λ1D +λ2D)cD (xD, tcrD) = (λ1D +λ2D)e−λ2DxcrD (135)

e a derivada na coordenada xD é obtida através de (132):

(∂σD

∂xD

)xD,tcrD(xD)

=−(λ1D +λ2D)

[e−λ2D xcrD +

(1+

λ2D

λ1D

)σcrD

]. (136)

Usando (134), (135) e (136) em (133) e rearranjando convenientemente obtêm-se:

dtcrD

dxD=

(1+

λ2D

λ1D

)σcrDeλ2DxD +1, (137)

o que forma um problema de valor inicial caso inclua-se:

tcrD (0) =σcrD

λ1D

. (138)

Resolver (137), (138) recai sobre a frente de erosão

tcrD (xD) =

(1

λ1D

+1

λ2D

)σcrD

(eλ2D xD−1

)+ xD +

σcrD

λ1D

, (139)

com a usual definição de xcrD:

xD = xcrD (tcrD (xD)) . (140)

No tocante a impedância, a integração procede segundo:

J (tD) = 1+φc0

β1

1∫0

σ1DdxD +β2

1∫0

σ2DdxD

, (141)

o que não pode ser feito analiticamente para todos os tempos. Foi utilizada então a integração

numérica por integração direta da Spline Cúbica Interpolante (RALSTON; RABINOWITZ, 2001);

que será descrita no apêndice C. Todavia, é possível determinar algumas informações sobre a curva

de impedância analiticamente: a reta que rege o início da deposição (que ainda pode ser modelada

pela teoria clássica) até tD =σcrDλ1D

vem dada por

39

J (tD)≈ 1+ c0φβ1λ1D +β2λ2D

λ1D +λ2D

(1− e−(λ1D+λ2D)

)tD, (142)

e a reta após a erosão em todo o testemunho é

J (tD) = 1+ c0φ

{β1σcrD +β2

[λ2D +1

λ2D

e−λ2D − 1λ2D

− σcrD

λ2D

(λ2D−1+ e−λ2D

)+(

1− e−λ2D

)tD

]}.

(143)

4.2 Geometria Radial

Para a filtração profunda na vizinhança de um poço injetor, onde a velocidade decresce rapida-

mente dentro da zona de retenção, o coeficiente de filtração é assumido constante (PANG; SHARMA,

1997; WENNBERG; SHARMA, 1997). O sistema de equações que modela a injeção de suspensões

em geometria radial, pode ser deduzido a partir das equações (1), (5), (6), (7) e (10) e vem dado por:

φ∂c∂ t

+q

2πrh∂c∂ r

=−∂σ

∂ t, (144)

∂σ

∂ t = cU (r)λ , σ ≤ σcr (U (r)) ,

σ = σcr (U (r)) , caso contrário,(145)

U =− k0

µ (1+βσ)

∂ p∂ r

, (146)

c(rw, t) = c0, (147)

c(r,0) = σ (r,0) = 0. (148)

Há de se lembrar que:

U (r) =q

2πrh, (149)

o que explica porque U (e conseqüentemente σcr) em (145) depende de r. Embora seja possível adi-

mensionalizar o sistema (144)-(148), é mais cômodo mantê-lo dimensional para permitir comparação

direta com dados de campo, a partir de um simulador construído para tal.

Novamente então recorre-se ao diagrama no plano (r, t) (Figura 12) para resolver o sistema de

equações (144)-(148) (detalhes no apêndice B), com η sendo o intercepto em t da característica que

passa pelo ponto (xD, tD) com a frente de erosão. A solução desse sistema é:

40

Figura 12: Diagrama (r, t) para a geometria radial

c(r, t) = σ (r, t) = 0 Região 0 (150)

σ (r, t) =c0qλ

2πhr

[t− φπh

q

(r2− r2

w)]

e−λ (r−rw) Região I (151)

c(r, t) = c0e−λ (r−rw) Região I (152)

0 = t−η− φπhq

(r2− r2

cr (η))

Região II (153)

σ (r, t) =rcr (η)σcr (rcr (η))

re(−λ (r−rcr(η))) Região II (154)

c(r, t) = c0e−λ (r−rcr(η)) Região II (155)

σ (r, t) = σcr (r) Região III (156)

c(r, t) = c0, Região III (157)

onde rcr (t) representa a posição (distância radial) da frente de erosão. Uma diferença qualitativa

importante deste caso para o linear é que as equações que determinam as regiões 0, I e II são agora

parábolas, dadas pelas seguintes equações:

t =φπh

q

(r2− r2

w), (158)

t =φπh

q

(r2− r2

w)+

2πrwhqλc0

σcr (rw) . (159)

41

Resta agora encontrar a frente de erosão. Para isso procede-se como nos casos anteriores: primeiro

escreve-se a equação (144) em termos de σ apenas (detalhes no apêndice B):

2φπrhλq

∂σ

∂ t+

1λ

∂σ

∂ r=−

(1+

1rλ

)σ , (160)

e depois escreve-se a derivada sobre a frente de erosão como:

(dσ

dr

)(r,tcr(r))

=

(∂σ

∂ r

)(r,tcr(r))

+

(∂σ

∂ t

)(r,tcr(r))

dtcr

dr. (161)

O primeiro membro da equação precedente é simplesmente:

(dσ

dr

)(r,tcr(r))

=dσcr (r)

dr, (162)

com a derivada temporal do segundo membro dada por (usando (145), (149)):

(∂σ

∂ t

)(r,tcr(r))

=qλ

2πrhc0 (163)

e a derivada espacial do mesmo membro escrita, com auxílio de (160), como:

(∂σ

∂ r

)(r,tcr(r))

=−(

λ +1r

)σcr (r)−λφc0. (164)

Uma diferença do caso radial para os casos lineares é que o primeiro membro de (161) admite uma

derivada não-nula a depender da curva σcr (r). Além disso, a citada curva aparece na expressão (164).

Isso impossibilita a integração analítica a não ser que a curva seja determinada de antemão, integrável

e a razão σcr(r)r também seja integrável. Embora este trabalho determine curvas teóricas de σcr (U),

a derivação a seguir não assumirá qualquer tipo especial de curva, mas imporá a condição de que a

derivada σcr (r) seja sempre não-negativa, o que é bem razoável de acordo com os experimentos de

vazão alternada.

Agora, faz-se uso de (162), (163) e (164) em (161) e obtêm-se a relação diferencial que define a

frente de erosão:

dtcr

dr=

2πrhqλc0

[σ′cr (r)+

(λ +

1r

)σcr (r)+λφc0

], (165)

onde σ ′cr (r) é a derivada da curva σcr (r) que forma um problema de valor inicial com a seguinte

relação:

tcr (rw) =2πrwhqλc0

σcr (rw) . (166)

42

Da imposição sobre a derivada, conclui-se que (165) será sempre não negativa e que tcr é monótona,

admitindo uma inversa, que é justamente rcr (t), ou seja:

rcr (tcr (r)) = r. (167)

A curva definida por (165) e (166) pode ser encontrada através de métodos numéricos de solução

de equações diferenciais ordinárias para virtualmente qualquer curva σcr (U) que se conheça. No

apêndice, descrever-se-á o método utilizado neste trabalho, que foi um Runge-Kutta de 4a ordem

(RUGGIERO; LOPES, 1996).

A impedância do caso radial é calculada numericamente pela expressão:

J (tD) = 1+β

ln(re/rw)

re∫rw

σdr, (168)

embora seja possível aproxima-la até antes da erosão (t < 2πrwhqλc0

σcr (rw)) pela reta enunciada pela

teoria clássica (ver em Altoé et al. (2004)).

43

5. RESULTADOS E DISCUSSÃO

Nesta secção serão apresentados os resultados obtidos a partir da teoria desenvolvida. Primeiro,

apresentar-se-ão os ajustes de dados obtidos a partir do modelo com um mecanismo de captura em

testemunhos, inclusive para o caso de vazão alternada. Após isso, será possível o ajuste de curvas

teóricas de concentração máxima retida em função da velocidade. Seguir-se-ão alguns resultados para

o caso a dois mecanismos também em testemunhos e finalmente algumas previsões de impedância

em poços.

5.1 Resultados para o Caso com um Mecanismo de Captura emTestemunhos

Nesta secção, os três parâmetros empíricos do modelo para transporte de colóides com liberação

de partículas - o coeficiente de filtração, o de dano à formação e a concentração máxima retida - são

calculados a partir da impedância obtida de testes laboratoriais realizados por Furtado et al. (2002,

2003, 2005). A litologia das amostras não é informada nos testes. Os testemunhos em questão foram

encapsulados em silicone, sem pressão confinante. A água da formação foi simulada por solução

salina com cloretos; concentração do ânion Cl− de 90000 ppm, concentração de Na+ de 31000 ppm,

além de outros cátions metálicos (Ca2+, K+, Mg2+, Sr+ e Ba+). A salinidade total foi de 148000

ppm e o pH foi de 6,78.

Foi injetada nos testemunhos água de formação até a estabilização da permeabilidade, para evitar

os efeitos da migração de finos no resultado dos testes. Em seguida, foi injetado óleo até a estabiliza-

ção da permeabilidade efetiva ou seja, até a permanência de água residual apenas. Posteriormente, o

óleo foi deslocado por água do mar filtrada até estabilização da permeabilidade efetiva, para simular

o efeito do óleo residual. A partir desse ponto, injetou-se a suspensão em questão, à uma vazão con-

stante, monitorando-se a diferença de pressão no testemunho apenas. A velocidade e outros dados do

testemunho estão expostos na Tabela 1. É possível obter os parâmetros λ , β e σcr a partir apenas de

dados de pressão, ao contrário daquilo que ocorre com o modelo clássico (BEDRIKOVETSKY et al.,

2001). Para isso, basta mensurar três características da curva de impedância: seu coeficiente angular

no período anterior à erosão:

m = φβc0 (1− exp(−λD)) , (169)

44

o tempo a partir do qual a impedância deixa de ser linear, ou seja, o tempo

tDc =σcrD

λD, (170)

e o tempo a partir do qual a impedância se estabiliza,

tDst =σcrD

λD+1+σcrD. (171)

O sistema (169)-(171) permite solução explícita:

σcrD = tDst − tDc−1, λD =tDst − tDc−1

tDc

, β =m[

φc0

(1− exp

(− tDst−tDc−1

tDc

))] . (172)

Tabela 1: Dados do teste e parâmetros de filtração

Parâmetro 10 teste 20 teste

Porosidade (%) 18,6 7,7

Permeabilidade Inicial (mD) 118 78

Concentração de Injeção (ppm) 50 50

Comprimento do Testemunho (m) 0,048 0,052

Diâmetro do Testemunho (m) 0,038 0,038

Coeficiente de Dano à Formação 400 170

Coeficiente de Filtração (1/m) 207 57,7

Ponto (U = 0(m/s),σcr = φ (1−φc)) (0; 0,1674) (0; 0,0693)

Ponto (Uteste (m/s),σcrteste) (0,001; 0,002) (0,0001; 0,0058)

m 0,003720 0,0006219

tDc 21 500

tDst 232 2001

A possibilidade de caracterizar completamente o sistema de dano à formação com erosão a partir

apenas da curva de diferença de pressão em testes de laboratório é uma vantagem da teoria proposta,

visto que, como ja foi dito, isto não pode ser feito no modelo clássico.

As Figuras 13.a e 13.b apresentam as impedâncias medidas e previstas pelo modelo para os testes

1 e 2 respectivamente, com a curva contínua correspondendo à curva prevista. As Figuras 14.a e 14.b

apresentam a curva de breaktrough prevista para cada um dos testes. A partir da Tabela 1 é possível

observar que ambos os testes exibem coeficientes de filtração e dano à formação dentro da faixa usual

de valores (PANG; SHARMA, 1997; BEDRIKOVETSKY et al., 2001) e que as concentrações máx-

45

imas retida se mantém a valores consideravelmente inferiores à porosidade, demonstrando o sentido

físico desse parâmetro para esses testes.

Em contraponto ao modelo clássico, a estabilização da perda de injetividade se dá de maneira

suave. Neste sentido, é possível afirmar que a teoria presente neste trabalho é um complemento ao

modelo clássico, devendo ser aplicado sempre que a perda de injetividade se der de forma gradual e

que nos casos de perda abrupta, aplica-se o usual sistema clássico. Em outras palavras, este modelo

não suplanta nem contém o modelo clássico, pois são aplicados em contextos distintos. Como os

parâmetros do modelo podem ser obtidos com apenas os testes de pressão, a comfirmação da aplica-

bilidade do modelo pode ser feita através do monitoramento da concentração suspensa à jusante do

testemunho, ou seja, da curva de breakthrough.

(a)

(b)

Figura 13: Curvas de impedância para o primeiro (a) e segundo (b) testes. Os losangos denotam ospontos experimentais e a linha contínua o ajuste realizado

46

(a)

(b)

Figura 14: Curvas de concentração em suspensão reduzida na saida do testemunho para o primeiro(a) e segundo (b) testes

5.2 Resultados do Teste Vazão Alternada

Para utilizar os dados do teste de vazão alternada (Figura 15), foi necessário assumir que existem

dois mecanimos de captura atuantes. Além disso, assumiu-se que é possível descontar a influência

do segundo mecanismo subtraindo dos dados de diferencial de pressão uma reta, referente ao período

de maior vazão, pertinente ao segundo mecanismo (não sujeito à erosão). Finalmente, considerou-se

que a concentração máxima retida do primeiro mecanismo é zero na vazão mais elevada, ou seja,

não há captura pelo primeiro mecanismo nesta vazão. A subtração da reta descrita no parágrafo

anterior dos dados de entrada gerou os pontos experimentais da Figura 16.a, que foi utilizada para

ajustar a impedância, sendo esta quantidade chamada de pseudo-impedância J∗. Os parametros do

ajuste, bem como alguns outros dados, estão na Tabela 2. A Figura 16.b mostra a concentração de

efluente prevista para este teste se apenas um mecanismo fosse atuante. Para aproximar a presença do

segundo mecanismo, multiplica-se a curva de breakthrough por exp(−λ2D), que é a concentração da

47

suspenção filtrada à taxa constante pelo 2 mecanismo.

Tabela 2: Dados do teste de vazão alternada

Parâmetro Valor

Porosidade (%) 29,96

Permeabilidade Inicial (mD) 1040

Concentração de Injeção (ppm) 600

Comprimento do Testemunho (m) 0,0519

Diâmetro do Testemunho (m) 0,037

Coeficiente de Dano à Formação 1750

Coeficiente de Filtração (1/m) 48,2

Ponto (U = 0(m/s), σcr = φ (1−φc)) (0; 0,2696)

Ponto (Uteste (m/s),σcrteste) (0,00155; 0,009348)

m 0,2888

tDc 20,8

tDst 73,8

Este teste revela também uma velocidade para qual a concentração máxima retida é nula, o que

permite impor uma condição extra sobre o ajuste de qualquer curva teórica (como a equação (56))

para σcr com a velocidade: este ajuste deve mostrar que, para uma velocidade igual ou menor aquela

corrspondente a vazão mais alta, σcr = 0. Esta condição é honrada no ajuste apresentado neste tra-

balho, e esse ajuste é exibido na próxima secção.

Figura 15: Dados de pressão do teste de vazão alternada

48

(a)

(b)

Figura 16: Ajuste dos dados do teste de vazão alternada: (a) ajuste da curva prevista (contínua) contraos pontos laboratoriais e (b) concentração prevista da suspensão na saída do testemunho

5.3 Curvas Teóricas para σcr

A partir dos parâmetros recolhidos nos testes das duas seções anteriores, é possível construir, para

cada exemplo, curvas teóricas de concentração máxima retida em função da velocidade e do número

de erosão, desde que se assuma um valor de porosidade do reboco interno (fixado em 10% em todos

os cálculos). Ao fixar-se um valor de porosidade do reboco interno, obtêm-se para cada teste dois

pontos: um obtido através do par de vazão - concentração máxima retida experimentais e outro que

é a concentração máxima retida à vazão nula (estimada a partir da porosidade do testemunho e do

reboco interno).

Com esses dois pontos, foi possível ajustar uma curva conforme aquela prevista pela equação

(50), manipulando-se os parâmetros pertinentes as forças em micro-escala (ver Tabela 3).

49

Tabela 3: Valores das constantes usadas em escala microscópica

Parâmetro 10 Teste 20 Teste Vazão Alternada

Constante de Hammaker (J) 2,00.10−21 1,82.10−21 1,95.10−21

Inverso do comprimento de Debye (1/m) 2,35.109 2,35.109 2,35.109

Potencial do Grão (V) -0,03 -0,03 -0,03

Potencial da Partícula (V) -0,05 -0,05 -0,05

Raio da Partícula (m) 4,45.10−7 7,00.10−7 1,25.10−6

Constante de Elevação 1190 1190 600

Constante de Arraste 60 60 40

Viscosidade (cP) 1,00 1,00 1,00

Densidade do Fluido (kg/m3) 1000 1000 1000

Densidade da Partícula (kg/m3) 2750 2750 2750

Força Eletrostática (N) 9,343.10−11 4,504.10−11 2,087.10−10

Força de Arraste (N) 5,089.10−11 2,500.10−11 1,114.10−10

Força de Elevação (N) 5,256.10−12 1,758.10−12 1,563.10−11

Gravidade (N) 6,337.10−15 2,467.10−14 1,405.10−13

Os dados utilizados para o ajuste dessas curvas se encontram na Tabela 3, bem como o cálculo das

forças no ponto experimental de cada teste. Os valores usados na Tabela foram estimados, pois seus

valores não foram reportados. Todavia, como os valores estimados desses parâmetros estão dentro da

faixa usual, é possível afirmar que as curvas teóricas de concentração retida validam qualitativamente

a teoria proposta.

A partir destas considerações, foi possível grafar as curvas teóricas dadas pela equação (50) com

os valores da Tabela 3 para cada teste, tanto em função da velocidade quanto em função do número

de erosão ε , equação (56), e realizar uma análise de sensibilidade para as quatro variáveis que se

revelaram mais influentes no comportamento de σcr: a constante de Hammaker, o raio da partícula

injetada e as constantes de elevação e de arraste. Estes gráficos são dados pelas Figuras (17)-(28).

50

(a)

(b)

Figura 17: Curvas de concentração máxima retida em função da velocidade (a) e do número de erosão(b), com análise de sensibilidade em relação à constante de Hammaker para o primeiro teste

51

(a)

(b)

Figura 18: Curvas de concentração máxima retida em função da velocidade (a) e do número de erosão(b), com análise de sensibilidade em relação ao raio da partícula para o primeiro teste

52

(a)

(b)

Figura 19: Curvas de concentração máxima retida em função da velocidade (a) e do número de erosão(b), com análise de sensibilidade em relação à constante de elevação χ para o primeiro teste

53

(a)

(b)

Figura 20: Curvas de concentração máxima retida em função da velocidade (a) e do número de erosão(b), com análise de sensibilidade em relação à constante de elevação ω para o primeiro teste

54

(a)

(b)

Figura 21: Curvas de concentração máxima retida em função da velocidade (a) e do número de erosão(b), com análise de sensibilidade em relação à constante de Hammaker para o segundo teste

55

(a)

(b)

Figura 22: Curvas de concentração máxima retida em função da velocidade (a) e do número de erosão(b), com análise de sensibilidade em relação ao raio da partícula para o segundo teste

56

(a)

(b)

Figura 23: Curvas de concentração máxima retida em função da velocidade (a) e do número de erosão(b), com análise de sensibilidade em relação à constante de elevação χ para o segundo teste

57

(a)

(b)

Figura 24: Curvas de concentração máxima retida em função da velocidade (a) e do número de erosão(b), com análise de sensibilidade em relação à constante de elevação ω para o segundo teste

58

(a)

(b)

Figura 25: Curvas de concentração máxima retida em função da velocidade (a) e do número deerosão (b), com análise de sensibilidade em relação à constante de Hammaker para o teste de vazõesalternadas

59

(a)

(b)

Figura 26: Curvas de concentração máxima retida em função da velocidade (a) e do número de erosão(b), com análise de sensibilidade em relação ao raio da partícula para o teste de vazões alternadas

60

(a)

(b)

Figura 27: Curvas de concentração máxima retida em função da velocidade (a) e do número deerosão (b), com análise de sensibilidade em relação à constante de elevação χ para o teste de vazõesalternadas

61

(a)

(b)

Figura 28: Curvas de concentração máxima retida em função da velocidade (a) e do número deerosão (b), com análise de sensibilidade em relação à constante de elevação ω para o teste de vazõesalternadas

62

A análise geral dessas Figuras permite concluir que, dentre todos os parâmetros, a constante

de Hammaker é aquele que tem maior peso para a concentração máxima retida como função da

velocidade. Porém, se tratarmos σcr como função do número de erosão a constante de arraste é que

assume a maior importância, embora acompanhada de perto pela constante de Hammaker. O raio

da partícula mostra influência razoável em σcr tanto no número de erosão quanto na velocidade. A

constante de lifting, por sua vez, quase não afeta a dependência da concentração máxima retida com

a velocidade, mas mostra importante contribuição para o número de erosão.

A análise dimensional da secção 3.3 antecipa que não é possível comparar as figuras referentes

à diferentes testes pois eles exibem diferentes porosidades. Mas mesmo dentro de um mesmo teste,

observa-se que a análise dimensional proposta não está completa: se estivesse, todas as curvas calcu-

ladas deveriam coincidir nos gráficos de σcr contra ε . Uma possível explicação para essa incoerência

é a diferença na escala na qual a análise foi feita daquela dos parâmetros sobre os quais foi testada a

sensibilidade. Em outras palavras, a sensibilidade testada nas Figuras (17)-(28) se referia a paramet-

ros que estão em uma escala menor que aquela apresentada pelas forças, que são a base da análise

apresentada.

Outra feição curiosa observada nesses gráficos é que o número de erosão ε aumenta com a

diminuição da velocidade, em contrapartida à proporcionalidade direta da força de arraste com U

apresentada na equação (33). Essa aparente contradição é explicada pela diminuição da área aberta

ao fluxo devido a deposição das partículas, o que pode ser contemplado no terceiro membro das

equações (46) e (50):

Fc f =ωπµr2

sUφH√

1−Σcr. (173)

De maneira geral o aumento da constante de Hammaker ocasiona o aumento da concentração

máxima retida (Figuras (17), (21) e (25)). Este aumento é justificado, pois este parâmetro esta rela-

cionado ao potêncial atrativo de London - van der Waals. Embora o potêncial de Born, que é repul-

sivo, também esteja relacionado à constante de Hammaker, sua interação é de pequeno alcance, não

influênciando muito no módulo da força máxima, e sim na localização (distância entre interfaces, h)

dessa força. No que tange ao número da erosão, o aumento da constante de Hammaker leva à redução

dessa quantidade, o que é explicado em parte pela presença da força elétrica no denominador. Por

outro lado, a diminuição desse número indica o aumento, por mais tímido que seja, da influência da

força de elevação. Na Figura 25 ocorre uma “inversão” da tendência descrita aqui, e a linha referente

a 50% do valor da constante de Hammaker aparece um pouco acima da curva referente aos 90%. Esta

inversão parece ter sido ocasionada pelo baixo valor de A123.

O raio da partícula é inversamente proporcional à concentração máxima retida (Figuras (18), (22)

e (26)). Isto provavelmente se deve ao aumento de área da partícula para a ação da força de arraste,

63

enquanto a força elétrica aumenta apenas linearmente com o aumento do raio (ELIMELECH et al.,

1995). As curvas contra o número de erosão parecem corroborar esta interpretação.

A constante de elevação tem praticamente nenhuma influência sobre a concentração máxima

retida (Figuras (19), (23) e (27)). Os valores maiores dessa constante levam a números de erosão

menores, o que é justificado para manter o balanço de torques da equação (50) (ver equaçao (58)).

A constante de arraste tem influência moderada sobre a concentração máxima de partículas reti-

das (Figuras (20), (24) e (28)). É possível concluir disso que ela é a mais importante das forças

hidrodinâmicas. Os curvas referentes ao número de erosão são também muito influênciadas por esse

parâmetro, de maneira que quanto maior a constante de arraste, maiores são os números de erosão

associados, o que corresponde ao comportamento esperado.

A forma dos gráficos nas Figuras 17.a-28.a lembra uma parábola, especialmente nas baixas ve-

locidades. Este comportamento das menores velocidades também não é inexperado: como a força de

arraste tem ordem 1 na velocidade e a força de elevação é de ordem 1,5, é esperado que essa última

decaia mais rapidamente com o decremento de U . A partir das Figuras 17.b-28.b é possível concluir

que à veocidades pequenas, a força de elevação pode ser desprezada. Ao se desprezar a força de

elevação na equação (50), é possível relacionar explicitamente a concentração máxima retida:

σcr = φ (1−φc)

(1− 3ω2π2µ2r4

sφ 2H2F2

eU2), (174)

ou seja, a dependência é quadrática, como na parábola. Como a correção imposta pela força de

elevação não é excessivamente severa (no caso mais agudo, a diferença relativa entre o valor ideal de

ε ≈ 0,57 e o menor valor das curvas de número de erosão ε ≈ 0,47, é menor que 20%), a forma da

curva ainda lembra a parábola. Outra conclusão importante é que, no limite de baixa deposição que é

utilizado na prática, a força de elevação não pode ser desconsiderada dos cálculos.

5.4 Previsão de Impedância para Testemunhos com Dois Mecan-ismos de Captura

Como ilustração do caso linear à dois mecanismos, plota-se o gráfico de impedância (Figura 29),

cujos parâmetros do modelo usado estão dados na Tabela 4. Observa-se neste caso que a estabilização

da injetividade tardará muito mais a ocorrer, visto que o segundo mecanismo diminui a quantidade de

partículas disponível para a filtração conforme a distância a partir da face de injeção do testemunho.

64

Tabela 4: Parâmetros da simulação de impedância para o caso linear à dois mecanismos

Parâmetros Valor

λ1 (1/m) 10

λ2 (1/m) 5

β1 100

β2 50

Concentração máxima retida 0,020

Porosidade (%) 20

Velocidade de Injeção (mm/s) 2,50

Comprimento do Testemunho (m) 0,075

Concentração da suspensão injetada 100 ppm

Figura 29: Impedância para o caso de dois mecanismos com erosão em testemunhos conforme osparâmetros da Tabela 4)

Além disso, é possível observar duas linhas retas, uma no início da injeção e outra após um tempo

bem mais longo. A linha reta dos tempos mais curtos é interpretada como a atuação simultânea dos

dois mecanismos de captura sem erosão, o que de fato se reduz ao sistema clássico, que exibe essse

comportamento linear. O comportamento linear observado nos tempos mais longos corresponde ao

fim da deposição do primeiro mecanismo com a atuação ininterrupta do segundo, que não é afetado

pela liberação de partículas. A premissa de coeficiente de filtração constante determina que o sistema

se comporte, no que tange ao segundo mecanismo, de maneira similar aquela prevista pelo modelo

clássico após o final da erosão, o que explica o comportamento linear observado.

65

5.5 Previsão de Impedância para Poço

Como exemplo de Impedância para o caso radial (poço), foi realizada uma simulação com os

parâmetros expostos na Tabela 5. Para essa simulação, foi considerada uma concentração máxima

retida média ao invés de considerar a dependência σcr (U). A Figura (30) apresenta os resultados

dessa simulação.

Tabela 5: Parâmetros usados na simulação de impedância para poços

Parâmetro Valor

Porosidade (%) 29,96

Concentração da suspensão injetada 100 ppm

Diâmetro do Poço (in) 7,5

Raio externo (m) 250

Espessura do intervalo (m) 50

Lambda (1/m) 48,2

Beta 50

Concentração máxima retida 0,02996

Vazão de injeção (m3/dia) 1000

Figura 30: Impedância para um poço sob as condições da Tabela 5

Novamente observa-se que a impedância tem crescimento gradual em oposição ao modelo clás-

66

sico de Perda de Injetividade. Nessa simulação, porém, a estabilização não seria alcançada no decorrer

da vida útil do poço. Todavia, a depender de σcr (U) e λ , é possível que a estabilização possa ocorrer

dentro da vida produtiva do poço, permitindo melhor avaliar a necessidade de intervenção para esse

poço.

67

6. CONCLUSÕES

1. A estabilização do diferencial de pressão pode ser explicada pelo balanço de torques das forças

que agem sobre a partícula depositada.

2. A solução analítica da geometria linear mostra que há uma frente de erosão de partículas, atrás

da qual não ocorre mais deposição.

3. O modelo proposto prevê uma tendência de estabilização de permeabilidade gradual; não sub-

stituindo o modelo clássico mas devendo ser aplicado em outros cenários, os quais apresentam

estabilização gradual da Perda de Injetividade.

4. A partir de testes laboratoriais, foram ajustadas curvas teóricas de concentração máxima retida,

o que validou o modelo proposto qualitativamente.

5. A partir do ajuste de dados laboratoriais foi possível determinar que a força de arraste é a força

hidrodinâmica mais importante no processo de libreação de partículas.

6. Esses testes laboratoriais mostram também que, para o limite de pequena deposição, a força de

elevação não pode ser desprezada do balanço de torques.

Trabalhos Publicados

BEDRIKOVETSKY, P.; SIQUEIRA, F.D.; FURTADO, C.A.; SOUSA, A.L.S. Quantitative The-

ory for fines Migration and Formation Damage, In: SPE International Symposium and Exhibition on

Formation Damage Control, SPE128384, 2010

BEDRIKOVETSKY, P.; SIQUEIRA, F. D.; FURTADO, C. A.; SOUZA, A. L. S. Modified parti-

cle detachment model for colloidal transport in porous media Transport in Porous Media, 2009

BEDRIKOVETSKY, P.; SIQUEIRA, F.D.; FURTADO, C.A.; SOUSA, A.L.S. Teoria Alternativa

à Estabilização da Perda da Injetividade em Poços Injetores, In: 5º Congresso Brasileiro de P & D

em Petróleo e Gas, 2009

SIQUEIRA, F.D.; BEDRIKOVETSKY, P.; ORTIZ, C.E.P. Teoria Alternativa de Perda de Inje-

tividade - Verificação do Conceito de Capacidade de Estocagem, In: Reunião Anual de Avaliação

PRH-ANP 2009, 2009

68

7. RECOMENDAÇÕES

• Validar a teoria em micro escala quantitativamente realizando experimentos cuidadosos, de

modo a reproduzir as condições por ela postuladas

• Reformular o critério de balanço de torques em outras topologias de meio poroso e/ou para

partículas distribuídas segundo tamanho, usando teoria da percolação

• Avaliar detalhatamente a física das soluções do caso radial e linear a dois mecanismos

• Aplicar a teoria aqui proposto na área de Migraçao de Finos.

69

REFERÊNCIAS BIBLIOGRÁFICAS

AKHATOV, I.; HOEY, J. M.; SWENSON, O.; SCHULS, D. L. Aerosol focusing on micro-capillaries:Theory and experiment. Journal of Aerosol Science, v. 39, p. 671–709, 2008.

AL-ABDUWANI, F. Filtration and External Filter Cake Build-up in Sandstones. Tese (Doutorado)— Delft University of Technology, The Netherlands, 2005.

AL-ABDUWANI, F.; BEDRIKOVETSKY, P.; FARAJZADEH, R.; BROEK, W. van den; CURRIE,P. External filter cake erosion: Mathematical model and experimental study. In: SPE paper 94635presented at the SPE 6th European Formation Damage Conference, Scheveningen, The Netherlands.[S.l.: s.n.], 2005.

AL-ABDUWANI, F. A. H.; ZWART, B.-R. de; FARAJZADEH, R.; BROEK, W. M. G. T.van den; CURRIE, P. K. Utilising static filtration experiments to test existing filtration theories forconformance. In: Produced Water Workshop, Aberdeen, UK. [S.l.: s.n.], 2004.

ALTMANN, J.; RIPPERGER, S. Particle deposition and layer formation at the cross flow microfiltration. Journal of Menbrane Science, v. 124, p. 119–128, 1997.

ALTOÉ, J. E.; BEDRIKOVETSKY, P.; GOMES, A. C. A.; SIQUEIRA, A. G.; SOUZA, A. L. S.Effects of oil-water mobility on injectivity impairment due to suspended particles. In: SPE paper88501 presented at the SPE Asia Pacific Oil and Gas Conference and Exhibition held in Perth,Australia. [S.l.: s.n.], 2004.

ALVAREZ, A. C.; BEDRIKOVETSKY, P.; HIME, G.; MARCHESIN, D.; RODRÍGUEZ, J. R. Afast inverse solver for the filtration function for flow of water with particles in porous media. Journalof Inverse Problems, v. 22, p. 69–88, 2006.

BARENBLATT, G. I. Scaling, self-similarity, and intermediate asymptotics. [S.l.]: CambridgeUniversity Press, 1996.

BARENBLATT, G. I.; ENTOV, V. M.; RYZHIK, V. M. Theory of Fluid Flows Through naturalRocks. [S.l.]: Kluwer Academic Publisher, 1990.

BAYGENTS, J. C.; GLYNN, G. R.; ALBINGER, O.; BIESEMEIER, B. K.; OGDEN, K. L.;ARNOLD, R. G. Variation of surface charge density in monoclonal bacterial populations:Implications for transport through porous media. Environmental Sci Technol, v. 32, p. 1596–1603,1998.

BEDRIKOVETSKY, P. Upscaling of stochastic micro model for suspension transport in porousmedia. Transport in Porous Media, v. 75, p. 335–369, 2008.

BEDRIKOVETSKY, P. G. Mathematical Theory of Oil and Gas Recovery. [S.l.]: Kluwer AcademicPublisher, 1992.

BEDRIKOVETSKY, P. G.; MARCHESIN, D.; CHECAIRA, F.; SERRA, A. L.; RESENDE, E.Characterization of deep bed filtration system from laboratory pressure drop measurements. Journalof Petroleum Science and Engineering, v. 64, n. 3, p. 167–177, 2001.

70

BOLSTER, C.; HORNBERGER, G.; MILLS, A. A method for calculating bacterial depositioncoefficients using the fraction of bacteria recovered from laboratory columns. Environ. Sci. Technol.,v. 32, p. 1329–1332, 1998.

BOLSTER, C. H.; MILLS, A. L.; HORNBERGER, G. M.; HERMAN, J. S. Spatial distribution ofbacteria experiments in intact cores. Water Resources Res., v. 35, p. 1797–1807, 1999.

BRADFORD, S.; BETTAHAR, M. Concentration dependent transport of colloids in saturatedporous media. Journal of Contaminant Hydrology, v. 82, p. 99–117, 2006.

BRADFORD, S.; BETTAHAR, M.; SIMUNEK, J.; GENUCHTEN, M. van. Straining andattachment of colloids in physically heterogeneous porous media. Vadose Zone Journal, v. 3, p.384–394, 2004.

BRADFORD, S.; SIMUNEK, J.; BETTAHAR, M.; GENUCHTEN, M. van; YATES, S. R. Modelingcolloid attachment, straining, and exclusion in saturated porous media. Environmental Science andTechnology, v. 37, p. 2242–2250, 2003.

BRADFORD, S.; YATES, S.; BETTAHAR, M.; SIMUNEK, J. Physical factors affecting thetransport and fate of colloids in saturated porous media. Water Resources, v. 38, p. 1327, 2002.

CAMESANO, T.; UNICE, K.; LOGAN, B. Blocking and ripening of colloids in porous media andtheir implications for bacterial transport. Colloids and Surfaces A: Physicochemical and EngineeringAspects, v. 160, p. 291–307, 1999.

CHAUVETEAU, G.; NABZAR, L.; COSTE, J.-P. Physics and modeling of permeability damageinduced by particle deposition. In: spe paper 39463. [S.l.: s.n.], 1998.

CIVAN, F. Reservoir Formation Damage: Fundamentals, Modeling, Assesment, and Mitigation. 2nded.. ed. [S.l.]: Gulf Professional Publishing, 2007.

CORAPCIOGLU, M. Y.; CHOI, H. Modeling colloid transport in unsaturated porous media andvalidation with laboratory column data. Water Resour. Res., v. 32, p. 3437–3449, 1996.

DABROWSKI, W. Consequences of the mass balance simplification in modelling deep filtration.Water Resources, v. 22, n. 10, p. 1219–1227, 1988.

DERJAGIN, B. V.; LANDAU, L. D. Theory of the stability of strongly charged lyophobic sols andof the adhesion of strongly charged particles in solutions of electrolytes. Acta physicochim. URSS,v. 14, n. 6, p. 633–662, 1941.

ELIMELECH, M. Particle deposition onideal collectors from dilute flowing suspensions:Mathematical formulation, numerical solution, and simulations. Journal of Separation Technology,v. 4, p. 186–212, 1994.

ELIMELECH, M.; JIA, X.; GREGORY, J.; WILLIAMS, R. Particle Deposition & Aggregation:Measurement, Modelling and Simulation (Colloid & surface engineering). [S.l.]: Butterworth-Heinemann, 1995.

ENTOV, V. M.; MIRZHADZHANZADE, A. K. Hidrodinâmica em Perfuração. [S.l.]: Nedra (emRusso), 1990.

FOPPEN, J.; MPOROKOSO, A.; SCHIJVEN, J. Determining straining of escherichia coli frombreakthrough curves. Journal of Contaminant Hydrology, v. 76, p. 191–210, 2005.

71

FOPPEN, J. W. A.; SCHIJVEN, J. F. Evaluation of data from the literature on the transport andsurvival of escherichia coli and thermotolerant coliforms in aquifers under saturated conditions.Water Resources, v. 40, p. 401–426, 2006.

FREITAS, A. M.; SHARMA, M. M. Detachment of particles from surfaces: An afm study. Journalof Colloid and Interface Science, v. 233, p. 73–82, 2001.

FRIMMEL, F. H.; KAMMER, F. von der; FLEMMING, F. C. Colloidal Transport in Porous Media.[S.l.]: Springer-Verlag Berlin Heidelberg, 2007.

FURTADO, C. J. A.; CARVALHO, M. R.; RONDINI, S. B. Avaliação da Injeção de Água Produzidano Campo de Pampo. [S.l.], 2002.

FURTADO, C. J. A.; CARVALHO, M. R.; RONDINI, S. B. Avaliação e Previsão da Reinjeção emPampo. [S.l.], 2003.

FURTADO, C. J. A.; CARVALHO, M. R. de; RONDININI, S. B.; SOUZA, A. L. S. de. Advancedexperimental tests for evaluation of water injection and re-injection. In: Annual Seminar on Reservesand Reservoirs, Rio de Janeiro, Brazil. [S.l.: s.n.], 2005.

GREGORY, J. Approximate expressions for retarded van der waals interaction. Journal of Colloidand Interface Science, v. 83, p. 138–145, 1981.

GUEDES, R. G.; AL-ABDUWANI, F.; BEDRIKOVETSKY, P.; CURRIE, P. Injectivity declineunder multiple particle capture mechanisms. In: SPE International Symposium and Exhibition onFormation Damage Control held in Lafayette, Lousiana. [S.l.]: SPE98623, 2006.

HARVEY, R. W.; GARABEDIAN, S. P. Use of colloid filtration theory in modeling movement ofbacteria through a contaminated sandy aquifer. Environ. Sci. Technol., v. 25, p. 178–185, 1991.

HERZIG, J. P.; LECLERC, D. M.; GOFF, P. le. Flow of suspensions through porous media -application to deep bed filtration. Industrial and Engineering Chemistry, v. 65, n. 5, p. 8–35, 1970.

HOFFMAN, J. D. Numerical Methods for Engineers and Scientists. [S.l.]: Marcel Dekker, 2001.

ISRAELACHVILI, J. Intermolecular and Surface Forces. [S.l.]: Academic Press, 2006.

JIAO, D.; SHARMA, M. M. Mechanism of cake build-up in cross flow filtration of colloidalsuspensions. Journal of Colloid and Interface Science, v. 162, p. 454–462, 1994.

JOHNSON, W. P.; LI, X.; ASSEMI, S. Deposition and re-entrainment dynamics of microbes andnon-biological colloids during non-perturbed transport in porous media in the presence of an energybarrier to deposition. Advances in Water Resources, v. 30, p. 1432–1454, 2007.

KANG, S.-T.; SUBRAMANI, A.; HOEK, E.; DESHUSSES, M.; MATSUMOTO, M. Directobservation of biofouling in cross-flow microfiltration: mechanisms of deposition and release.Journal of Membrane Science, v. 244, p. 151–165, 2004.

KHILAR, K.; FOGLER, S. Migration of Fines in Porous Media. [S.l.]: Kluwer Academic Publisher,1998.

KIM, S.; KARRILA, S. J. Microhydrodynamics: Principles and Selected Applications. [S.l.]: DoverPublications, 2005.

72

KRETZSCHMAR, R.; BARMETTLER, K.; GROLIMUND, D.; YAN, Y.-d.; BORKOVEC, M.;STICHER, H. Experimental determination of colloid deposition rates and collision efficiencies innatural porous media. Water Resources Research, v. 33, p. 1129, 1997.

KUHNEN, F.; BARMETTLER, K.; BHATTACHARJEE, S.; ELIMELECH, M.; KRETZSCHMAR,R. Transport of iron oxide colloids in packed quartz sand media: Monolayer and multilayerdeposition. Journal of Colloid and Interface Science, v. 231, p. 32–41, 2000.

LANDAU, L. D.; LIFSHITZ, E. M. Statistical Physics, Part 1. [S.l.]: Pergamon Press, Oxford, 1980.

LANDAU, L. D.; LIFSHITZ, E. M. Fluid Mecanics. 2nd edition. ed. [S.l.]: Pergamon Press, 1987.(Course on Theoretical Physics, v. 6).

LOGAN, D. J. Transport Modeling in Hydrogeochemical Systems. [S.l.]: Springer, 2001.

MASSEI, N.; LACROIX, M.; WANG, H. Q.; DUPONT, J. Transport of particulate material anddissolved tracer in a highly permeable porous medium: comparison of the transfer parameters.Journal of Contaminant Hydrology, v. 57, p. 21–39, 2002.

MAYS, D.; HUNT, J. R. Hydrodynamic aspects of particle clogging in porous media. Journal ofEnvironmental Science and Technology, v. 39, p. 577–584, 2005.

MIRANDA, R. M.; UNDERDOWN, D. R. Laboratory measurement of critical rate: A novelapproach for quantifying fines migration problems. In: SPE paper 25432 presented at the ProductionOperations Symposium held in Oklahoma City. [S.l.: s.n.], 1993.

O’NEIL, M. E. A sphere in contact with a plane wall in a slow linear shear flow. ChemicalEngineering Science, v. 23, p. 1293–1298, 1968.

PANG, S.; SHARMA, M. M. A model for predicting injectivity decline in water-injection wells.SPEFE - 28489, v. 12, p. 194–201, 1997.

PAYTAKES, A. C.; RAJAGOPALAN, R.; TIEN, C. Application of porous medium models to thestudy of deep bed filtration. The Canadian Journal of Chemical Engineering, v. 727, p. 52–56, 1974.

PLUEN, P. Oilfield water management: Total’s current practices & perspectives. In: Seminar onOilfield Water Management, Rio de Janeiro, October 17th 2007. [S.l.: s.n.], 2007.

RALSTON, A.; RABINOWITZ, P. A First Course in Numerical Analysis. [S.l.]: Dover Publications,2001.

ROQUE, C.; CHAUVETEAU, G.; THIBAULT, G.; BOUTECA, M. Mechanisms of formationdamage by retention of particles suspended in injection water. In: SPE paper 30110, EuropeanFormation Damage Conference, The Hague. [S.l.: s.n.], 1995.

RUGGIERO, M. A. G.; LOPES, V. L. R. Cálculo Numérico: Aspectos Teóricos e Computacionais.[S.l.]: Pearson Education, 1996.

SAFFMANN, P. G. The lift on a small sphere in a slow shear flow. Journal of Fluid Mechanics,v. 22, p. 385–400, 1965.

SAFFMANN, P. G. Correction to "the lift on a small sphere in a slow shear flow". Journal of FluidMechanics, v. 31, p. 624, 1968.

SANTOS, A.; BEDRIKOVETSKY, P. A stochastic model for particulate suspension flow in porousmedia. Journal of Transport in Porous Media, v. 13, p. 30–52, 2006.

73

SHARMA, M. M.; YORTSOS, Y. C. Fines migration in porous media. AIChE J., v. 33, p.1654–1662, 1987.

SHARMA, M. M.; YORTSOS, Y. C. A network model for deep bed filtration processes. AIChE J.,v. 33, p. 1644–1653, 1987.

SHARMA, M. M.; YORTSOS, Y. C. Transport of particulate suspensions in porous media: Modelformulation. AIChE J., v. 33, p. 1636, 1987.

SIMONI, S. F.; HARMS, H.; BOSMA, T. N. P.; ZEHNDER, A. J. B. Population heterogeneityaffects transport of bacteria through sand columns at low flow rates. Environ. Sci. Technol., v. 32, p.2100–2105, 1998.

SOO, H.; WILLIAMS, M.; RADKE, C. J. A filtration model for flow of dilute, stable emulsions inporous media - 2. parameter evaluation and estimation. Chem Eng Sci, v. 41, p. 261–272, 1986.

TUFENKJI, N.; ELIMELECH, M. Correlation equation for predicting single-collector efficiencyin physicochemical filtration in sarurated porous media. Journal of Environmental Science andTechnology, v. 38, p. 529–536, 2004.

TUFENKJI, N.; ELIMELECH, M. Deviation from classical colloid filtration theory in the presenceof repulsive dlvo interactions. Langmuir, v. 20, p. 10818–10828, 2004.

TUFENKJI, N.; ELIMELECH, M. Breakdown of colloid filtration theory: Role of the secondaryenergy minimum and surface charge heterogeneities. Langmuir, v. 21, p. 841–852, 2005.

TUFENKJI, N.; REDMAN, J. A.; ELIMELECH, M. Interpreting deposition patterns of microbialparticles in laboratory-scale column experiments. Environ. Sci. Technol., v. 37, p. 616–623, 2003.

VEERAPEN, J.; NICOT, B.; CHAUVETAU, G. In-depth permeability damage by particle depositionat high flow rates. In: SPE paper 68962 presented at presented at the SPE European FormationDamage Conference to be held in The Hague, The Netherlands 21-22 May. [S.l.: s.n.], 2001.

WENNBERG, K. E.; SHARMA, M. M. Determination of the filtration coefficient and the transitiontime for water injection wells. In: Paper SPE 38181 Presented At The SPE European FormationDamage Conference. [S.l.: s.n.], 1997.

74

APÊNDICE A. CASO LINEAR A DOIS

MECANISMOS COM EROSÃO

O sistema a ser resolvido é:

∂

∂ tD(cD +σ1D +σ2D)+

∂cD

∂xD= 0, (175)

∂σ1D

∂ tD= λ1DcD, σ1D < σcrD (ε) , (176)

∂σ2D

∂ tD= λ2DcD, (177)

U =− k0

µL(1+φc0β1σ1D +φc0β2σ2D)

∂ p∂xD

, (178)

cD (x,0) = σ1D (x,0) = σ2D (x,0) = 0, (179)

cD (0, t) = 1. (180)

Os passos do detalhamento serão:

1.Utilizar o Lema de Agregação de Guedes et al. (2006) para reduzir o sistema à um mecanismo.

2.Obter a equação análoga à (175) mas apenas em função de σ .

3.Utilizar o método das características para resolver o sistema de equações (175)-(180).

4.Obter a equação da impedância (equação (141)).

A.1 Aplicação do Lema de Agregação

Para efetuar a agregação, escreve-se:

75

dσ1D

dσ2D

=

∂σ1D∂ tD

∂σ2D∂ tD

=H(σcrD−σ1D)λ1D

λ2D

, (181)

onde H(x) é a função passo unitário de Heavside.

Assim, se σcrD > σ1D , têm-se:

dσ1D

dσ2D

=λ1D

λ2D

,

σ1D =λ1D

λ2D

σ2D, (182)

σD = σ1D +σ2D,

σ1D = σD(1+

λ2Dλ1D

) σ2D = σD(1+

λ1Dλ2D

), (183)

e se σcrD = σ1D , obtêm-se:

σD = σcrD +σ2D . (184)

A equação (183) torna possível a correspondência entre a quantidade máxima retida do primeiro

mecanismo e a quantidade equivalente total, ou seja:

σ (σ1D = σcrD) =

(1+

λ2D

λ1D

)σcrD. (185)

Dessa forma construímos a função de taxa de captura como a soma das taxas de cada mecanismo em

função de σ , ou seja:

λD = λ1DH((

1+λ2D

λ1D

)σcrD−σ

)+λ2D. (186)

A.2 Obtenção da Equação Análoga à (175) em Função de σD So-mente

Primeiramente explicita-se c na cinética de captura:

76

cD =1

λD (σ)

∂σD

∂ tD. (187)

Agora introduz-se a variável Φ definindo-a como:

Φ(σD) =

σD∫0

dyλD (y)

. (188)

A partir da qual expressa-se cD como:

cD =∂Φ(σD)

∂ tD. (189)

Substituindo (189) em (175) obtêm-se:

∂

∂ tD

(∂Φ(σD)

∂ tD+σ

)+

∂

∂xD

(∂Φ(σD)

∂ tD

)= 0,

que após a inversão da ordem de derivação do segundo termo do primeiro membro fica:

∂

∂ tD

(∂Φ(σD)

∂ tD+σ

)+

∂

∂ tD

(∂Φ(σD)

∂xD

)= 0 (190)

Agora integra-se (190) com respeito à tD e obtêm-se:

[∂Φ(σD)

∂ tD+σD

]tD

0+

[∂Φ(σD)

∂xD

]tD

0= 0,

0 =1

λD (σD)

∂σD

∂ tD+σD− cD (xD,0)−σD (xD,0)+

+1

λD (σD)

∂σD

∂xD− 1

λD (σ (xD,0))dσD (xD,0)

dxD,

aplicando as condições iniciais e de contorno concluí-se que:

∂σD

∂ tD+

∂σD

∂xD=−λD (σD)σD. (191)

A equação (191) é a equação correspondente à (175) em termos de apenas σ . Como a ordem da

derivada em (190) foi invertida, (191) vale até onde a frente de erosão ainda não passou. Em outras

palavras, a validade dessa expressão é mantida sempre que a função de filtração for contínua, e desde

que se mantenham as condições iniciais e de contorno.

77

A.3 Solução do Sistema (175)-(180) Utilizando o Método das Car-acterísticas

Para resolver o sistema (175)-(180) é necessário empregar o método das características. A

primeira análise que se deve fazer no diagrama (xD, tD) (Figura 11) é dividi-lo em duas partes: uma

região diretamente influenciada pelas condições iniciais (região 0) e outra diretamente afetada pelas

condições de contorno (regiões I, II e III).

A.3.1 Região 0

Na região 0, escreve-se:

dσD

dtD=

∂σD

∂ tD+

dxD

dtD

∂σD

∂xD, (192)

e a comparação com (191) determina que:

dσD

dtD=−(λ1D +λ2D)σD (193)

dxD

dtD= 1. (194)

A equação (194) implica, por integração direta, que:

xD−ξD = tD, (195)

onde ξ se refere a intersecção da característica do ponto (xD, tD) com o eixo xD. De (193) segue que:

σD (xD, tD) = σD (ξD,0)e−(λ1D+λ2D)tD ,

e incluíndo (179) conclui-se que:

σD (xD, tD) = 0⇒ σ1D (xD, tD) = σ2D (xD, tD) = 0, (196)

e de (187) conclui-se que:

cD = 0. (197)

Para finalizar, vamos definir a divisa entre a região 0 e I. Para isso, basta utilizar (195) com ξ = 0,

o que ocasiona em (130):

78

xD = tD. (198)

A.3.2 Região I

Na região I utiliza-se o método das características como segue:

dσD

dxD=

∂σD

∂xD+

dtDdxD

∂σD

∂ tD, (199)

e a comparação de (199) com (191) permite aferir que:

dσD

dxD=−(λ1D +λ2D)σD (200)

dtDdxD

= 1. (201)

De (201) conclui-se que:

tD−ηD = xD. (202)

E de (200):

σD (xD, tD) = σD (0,ηD)e−(λ1D+λ2D)xD , (203)

com σD (0,ηD) sendo dado pela integração de (176) e (177) somadas na cordenada xD = 0:

σD (0,ηD)−σD (0,0) =

η∫0

(λ1D +λ2D)dtD,

σD (0,ηD) = (λ1D +λ2D)ηD. (204)

A utilização de (204) em (203) resulta em:

σD (xD, tD) = (λ1D +λ2D)(tD− xD)e−(λ1D+λ2D)xD . (205)

Finalmente a aplicação de (183) em (205) fornece:

σ1D = λ1D (tD− xD)e−(λ1D+λ2D)xD, (206)

σ2D = λ2D (tD− xD)e−(λ1D+λ2D)xD. (207)

79

Para obter cD basta somar (176) e (177) para usar com (205):

cD = e−(λ1D+λ2D)xD. (208)

A.3.3 Região II

Nesta região ainda valem as relações (200) e (201). Porém, deve-se levar em conta que a frente de

erosão já se propaga no meio poroso, e que sua presença muda sobremaneira o processo de deposição

de partículas. Em primeiro lugar, avaliar-se-á a integração de (201):

tD−η = xD− xcrD (ηD) ,

tD− xD + xcrD (η)−ηD = 0, (209)

ou seja, a coordenada de intersecção entre a característica e a frente de erosão ηD é definida implici-

tamente por (209). Agora que se passe para o caso da integração de (200):

σD (xD, tD) = σD (xcrD (ηD) ,ηD)e−(λ1D+λ2D)(xD−xcrD(ηD)), (210)

por (183), (182) e pelo fato de que sobre xcrD , σD (σ1D = σcrD) =(

1+λ2Dλ1D

)σcrD , segue:

σD (xD, tD) =(

1+λ2D

λ1D

)σcrDe−(λ1D+λ2D)(xD−xcrD(ηD)), (211)

σ1D (xD, tD) = σcrDe−(λ1D+λ2D)(xD−xcrD(ηD)), (212)

σ2D (xD, tD) =λ2D

λ1D

σcrDe−(λ1D+λ2D)(xD−xcrD(ηD)). (213)

Para encontrar cD desta vez é melhor voltar para (175), substituir (176), (177) e integrar a ex-

pressão resultante pelo método das características. Este procedimento mostrará que a velocidade da

característica resultante é a mesma que em (201), mas o diferencial em cD vem dado por:

dcD

dxD=−(λ1D +λ2D)cD, (214)

o que implica em:

80

cD (xD, tD) = cD (xcrD (ηD) ,ηD)e−(λ1D+λ2D)(xD−xcrD(ηD)), (215)

e com (127) (ver subsecção A.3.4) obtêm-se finalmente:

cD (xD, tD) = e−[λ2DxD+λ1D(xD−xcrD(ηD))]. (216)

A.3.4 Região III

Nesta região a frente de erosão já passou, o que significa que só ha deposição pelo segundo

mecanismo de captura. Assim, escreve-se (175) como:

∂cD

∂ tD+

∂cD

∂xD=−∂σ2D

∂ tD=−λ2DcD. (217)

Aplicando o método das características novamente obtêm-se:

dcD

dxD=−λ2DcD (218)

dtDdxD

= 1⇒ tD−ηD = xD (219)

A integração de (218) resulta em (com uso de (219)):

cD = e−λ2D xD. (220)

Para obter σ2D é necessário integrar (217) de tcrD (xD) até tD, mantendo xD fixo:

σ2D (xD, tD)−σ2D (xD, tcrD (xD)) =

tD∫tcrD(xD)

λ2De−λ2DxDdη

usando (185):

σ2D (xD, tD) =λ2D

λ1D

σcrD +λ2D (tD− tcrD (xD))e−λ2DxD. (221)

A solução σ1D é trivial:

σ1D = σcrD . (222)

81

A.4 Impedância

Para obter a fórmula da impedância, recorre-se à uma de suas definições:

J (tD) =∆p(tD)∆p(0)

. (223)

A partir de (178) escreve-se:

∆p(tD) =UµL

k0

1∫0

(1+φc0β1σ1D +φc0β2σ2D)dxD

=UµL

k0

1+φc0β1

1∫0

σ1DdxD +φc0β2

1∫0

σ2DdxD

.em tD = 0, σ1D = σ2D = 0, logo:

∆p(0) =UµL

k0,

o que acarreta:

J (tD) = 1+φc0β1

1∫0

σ1DdxD +φc0β2

1∫0

σ2DdxD. (224)

82

APÊNDICE B. CASO RADIAL COM

EROSÃO

As equações básicas pertinentes a geometria radial com erosão são:

φ∂c∂ t

+q

2πrh∂c∂ r

=−∂σ

∂ t, (225)

∂σ

∂ t= cU (r)λ (σ ;σcr (r)) , (226)

U =− k0

µ (1+βσ)

∂ p∂ r

, (227)

c(0, t) = c0, (228)

c(r,0) = σ (r,0) = 0. (229)

Seguir-se-ão agora os mesmos passos realizados no apêndice A, exceto o primeiro, visto que este caso

trata de um mecanismo apenas.

B.1 Obtenção da Equação Análoga à (225) em Função de σ So-mente

Em primeiro lugar, substitui-se a expressão para velocidade da geometria radial (149) em (226)

para obter:

∂σ

∂ t=

qcλH(σcr (r)−σ)

2πrh, (230)

agora é explicitar c em (230) para obter (assumindo que σ (r, t)< σcr (r)):

c =2πrhλq

∂σ

∂ t. (231)

83

A substituição de (231) em (225) resulta em:

2φπrhλq

∂ 2σ

∂ t2 +1λ

∂ 2σ

∂ r∂ t=−

(1+

1rλ

)∂σ

∂ t. (232)

Invertendo a ordem da derivada no segundo termo do primeiro membro de (232), e integrando com

respeito do tempo no intervalo [0; t] obtêm-se:

2φπrhλq

∂σ

∂ t+

1λ

∂σ

∂ r=−

(1+

1rλ

)σ . (233)

B.2 Solução do Sistema (225)-(229) pelo Método das Caracterís-ticas

Procede-se nesta secção da mesma maneira que na Secção A.3 (ver Figura 12).

B.2.1 Região 0

Para a região 0, faz-se uso da seguinte expressão:

dcdt

=∂c∂ t

+∂c∂ r

drdt

, (234)

o que acarreta em (quando aplicada em (225)):

drdt

=q

2φπrh, (235)

dcdt

=qcλ

2πrh. (236)

A integração de (235) resulta em:

r2−ξ2 =

qφπh

t, (237)

ou seja, as características no espaço (r, t) são parábolas. Integrando (236) obtêm-se:

c(r, t) = c(ξ ,0)exp

λqφπh

ξ −√

φπhq t +ξ 2

2

, (238)

o que acaba, com uso de (229) e (237), culminando em:

84

c(r, t) = σ (r, t) = 0. (239)

A linha que delineia a região 0, juntamente com o eixo r é dada por (237) com ξ = 0, resultando

em (158).

B.2.2 Região I

Nesta região escreve-se:

dσ

dr=

∂σ

∂ r+

∂σ

∂ tdtdr

, (240)

donde concluí-se que:

dtdr

=2φπrh

q, (241)

dσ

dr=−

(λ +

1r

)σ . (242)

A integração de (241) resulta em:

t−η =πφh

q

(r2− r2

w), (243)

e de (242):

σ (r, t) =rwσ (rw,η)

rexp [−λ (r− rw)] . (244)

Para determinar a quantidade σ (rw,η), procede-se integrando a cinética de captura (230) em

r = rw:

σ (rw,η) =qλc0

2πrwhη =

qλc0

2πrwh

[t− φπh

q

(r2− r2

w)]

, (245)

a substituição de (245) em (244) resulta em:

σ (r, t) =c0qλ

2πhr

[t− φπh

q

(r2− r2

w)]

e−λ (r−rw), (246)

e usando a expressão (231) concluímos que:

c = c0e−λ (r−rw). (247)

85

Para encontrar a linha limítrofe entre as Regiões I e II, basta substituir η por tcr (0) em (243). A

determinação de tcr (0) é feita a partir de (245), com σ (rw, t) = σcr (rw):

tcr (0) =2πrwhqλc0

σcr (rw) , (248)

e a substituição de (248) em (243) resulta em (159).

B.2.3 Região II

Para esta região, ainda são válidas as relações (241) e (242). Porém, deve-se alterar (243) e (244),

devido à mudança no limite inferior de integração em r, que agora jaz sobre a frente de erosão. O

resultado final disso é:

t−η =φπh

q

(r2− r2

cr (η)), (249)

σ (r, t) =rcr (η)σ (rcr (η) ,η)

rexp [−λ (r− rcr (η))] . (250)

como σ (rcr (η) ,η) = σcr (rcr (η)):

σ (r, t) =rcr (η)σcr (rcr (η))

rexp [−λ (r− rcr (η))] , (251)

com a concentração dada por:

c = c0e−λ (r−rcr(η)). (252)

A equação é obtida aplicando (214) à (225).

B.2.4 Região III

Nesta Região, a frente de erosão já passou e não há mais deposição alguma. Logo:

σ = σcr (r) , (253)

c = c0. (254)

86

B.3 Impedância

Para obter a fórmula básica da impedância, o primeiro passo é manipular (227) para obter:

∆p(t) =− µq2πhk0

lnre

rw+β

re∫rw

σ (r, t)dr

, (255)

e pela aplicação de (255) em (223) concluí-se que:

J (tD) = 1+β

ln(re/rw)

re∫rw

σ (r, tD)dr. (256)

O tempo adimensional tD (volumes porosos injetados) é definido por:

tD =qt

φπh(r2e − r2

w). (257)

87

APÊNDICE C. FERRAMENTAS

MATEMÁTICAS

Neste apêndice serão descritos os artifícios numéricos usados na implementação dos simuladores.

A Secção C.1 tratará das Splines Cúbicas. Em C.2 será tratado o método de Runge-Kutta de 4a ordem,

usado na solução de (165)-(166).

C.1 Spline Cúbica Interpolante

Em alguns momentos os simuladores de perda de injetividade utilizados neste trabalho foi necessário

interpolar um certo conjunto de dados. O solução utilizada nestes casos foi a Spline Cúbica In-

terpolante. Para a implementação desse módulo, foram usadas as seguintes equações (RALSTON;

RABINOWITZ, 2001):

hi = ai−ai−1, i = 1,2, · · · ,n−1, (258)

S (x) = Mi−1(ai− x)3

6hi+Mi

(x−ai−1)3

6hi

+

(yi−1−

Mi−1h2i

6

)(ai− x)+

(yi−

Mih2i

6

)(x−ai−1) , (259)

onde x ∈ Ii = [ai−1,ai] e os n ai‘s representam as abcissas das n ordenadas yi da curva que se pretende

interpolar. Os coeficientes Mi são definidos pela solução do sistema abaixo:

2M0 +ω0M1 = d0, (260)

τiMi−1 +2Mi +ωiMi+1 = di, i = 1,2, · · · ,n−2, (261)

τn−1Mn−2 +2Mn−1 = dn, (262)

88

com as seguintes definições inclusas:

θi =yi− yi−1

hi, i = 1,2, · · · ,n−1, (263)

ωi =hi+1

hi +hi+1, τi = 1−λi, di =

6(θi+1−θi)

hi +hi+1, i = 1,2, · · · ,n−2. (264)

Agora, é possível construir pelo menos dois tipos de Spline: caso se conheça as derivadas da curva a

ser interpolada(y′0 ey′n−1

), faz-se:

ω0 = µn−1 = 1, (265)

dn−1 =6

hn−1

(y′n−1−

yn−1− yn−2

hn−1

), (266)

d0 =6h1

(y1− y0

h1− y′0

). (267)

Caso as derivadas não possam ser utilizadas, ainda é possível construir a Sline Cúbica Natural, cuja

premissa básica é S′′ (a0) = S′′ (an−1) = 0:

λ0 = d0 = µn−1 = dn−1 = 0. (268)

O sistema de equações lineares (260)-(262) apresenta-se com matriz dos coeficientes tridiago-

nal. Esta forma especial admite um algoritmo bastante eficiente para sua solução, conhecido como

Algoritmo de Thomas (ver Hoffman (2001)).

C.2 Runge-Kutta de 4a ordem

Para resolver o problema de valor inicial:

dydx

= f (x,y) , (269)

y(a) = y0, (270)

nos pontos (xi,yi) , i = [0,n−1], espaçados de hi = xi+1−xi nas abcissas e com x ∈ [a,b], utiliza-se o

seguinte esquema (RALSTON; RABINOWITZ, 2001):

89

yi+1 = yi +0,17476028k1 +0,55148066k2 +1,20553560k3 +0,17118478k4, (271)

k1 = hi f (xi,yi) ,

k2 = hi f (xi +0,4hi,yi +0,4k1) ,

k3 = hi f (xi +0,45573725hi,yi +0,2969776k1 +0,15875964k2) ,

k4 = hi f (xi +hi,yi +0,21810040k1−3,05096516k2 +3,832864k3) .

a vantagem deste esquema em vista dos outros de 4a ordem é que ele minimiza o coeficiente de erro

na quarta ordem (γ4 em Ralston e Rabinowitz (2001)).