José Roberto Silvestre

Análise numérica de estabilidade de poços de petróleo

com relevância a produção de areia

Dissertação de Mestrado

Dissertação apresentada ao Programa de Pós-Graduação em Engenharia Civil da PUC-Rio como requisito parcial para obtenção do título de Mestre em Engenharia Civil. Área de Concentração: Geotecnia

Orientadores: Eurípedes do Amaral Vargas Jr. João Luiz Elias Campos

Rio de Janeiro, março de 2004.

José Roberto Silvestre

Análise numérica de estabilidade de poços de petróleo com relevância a produção de areia

Dissertação apresentada como requisito parcial para obtenção do título de Mestre pelo Programa de Pós-Graduação em Engenharia Civil da PUC-Rio. Aprovada pela Comissão Examinadora abaixo assinada.

Eurípedes do Amaral Vargas Jr. Orientador

PUC/Rio

João Luiz Elias Campos TecGraf

Luiz Fernando Campos Ramos Martha PUC/Rio

Aldo Durand Fárfan UENF

José Eugênio Leal Coordenador Setorial do Centro Técnico Científico - PUC-Rio

Rio de Janeiro, 12 de março de 2004

Todos os direitos reservados. É proibida a reprodução total ou parcial do trabalho sem autorização da universidade, do autor e do orientador.

José Roberto Silvestre Graduado em Engenharia Civil, pela UFJF - Universidade Federal de Juiz de Fora, em 2000.

Ficha Catalográfica

Silvestre, José Roberto

Análise numérica de estabilidade de poços de petróleo com relevância a produção de areia / José Roberto Silvestre; orientadores: Eurípedes do Amaral Vargas Jr; João Luiz Elias Campos. - Rio de Janeiro: PUC, Departamento de Engenharia Civil, 2004.

v., 133f.: il. ; 29,7 cm

Dissertação (mestrado) - Pontifícia Universidade Católica do Rio de Janeiro, Departamento de Engenharia Civil

Inclui referências bibliográficas.

1. Engenharia Civil - Teses. 2. Produção de areia. 3. Modelo elastoplástico Lade-Kim. 4. Mecanismos de ruptura. I. Vargas Jr., Eurípedes do Amaral. II. Campos, João Luiz Elias. III. Pontifícia Universidade Católica do Rio de Janeiro. Departamento de Engenharia Civil. IV. Título.

CDD: 6.24

A meus pais, A minha família,

por terem confiado em mim.

Agradecimentos

A meus pais Luís e Ana, pela atenção, carinho e ensinamentos sem os quais eu

não teria chegado aqui.

Ao Professor Eurípedes do Amaral Vargas Jr., pela paciência, confiança

depositada, orientação, estímulo e conhecimento transmitido neste convívio.

Ao meu co-orientador João Luiz pela participação em momentos decisivos da

dissertação.

A Armando Prestes, Raquel Veloso, Aldo, Rodrigo Peluci e Sergio Montoya que,

de alguma forma, contribuíram nesta dissertação.

Aos professores do Departamento de Engenharia Civil da PUC - Rio.

A todos os amigos, pela amizade e apoio nas horas difíceis, especialmente aos das

salas 317 e 609.

A todos os funcionários da PUC.

A PUC e a Capes pelos recursos financeiros a pesquisa.

A minha família pelo incentivo e apoio que me deram.

Resumo Silvestre, José Roberto; Vargas Jr., Eurípedes do Amaral; Campos, João Luiz Elias. Análise numérica de estabilidade de poços de petróleo com relevância a produção de areia. Rio de Janeiro, 2004. 133p. Dissertação de Mestrado - Departamento de Engenharia Civil, Pontifícia Universidade Católica do Rio de Janeiro.

Este trabalho visa simular situações que levam uma rocha a ruptura

considerando o fluxo de fluidos, em uma tentativa de correlacionar os

mecanismos envolvidos na ruptura ao processo conhecido como produção de areia

em rochas reservatório, ou seja, a produção de sólidos durante a extração de

hidrocarbonetos. Há algum tempo este assunto tem sido tratado do ponto de vista

da mecânica das rochas dada a similaridade do seu comportamento com outros

fenômenos estudados nesta área.

Como se trata de um processo de extração de fluidos, a sua influência deve

ser considerada no comportamento mecânico da rocha. Assim, simulações

numéricas, representativas de situações favoráveis a este fenômeno, foram

analisadas, utilizando-se, para tal, o programa comercial de elementos finitos

ABAQUS, que permite a simulação do processo de acoplamento fluido mecânico.

Dada a flexibilidade demonstrada pelo pacote numérico, possibilitou-se a

implementação de um modelo elastoplástico e a verificação de sua eficiência na

representação do comportamento da rocha reservatório neste processo. A

implementação numérica baseou-se no conhecido modelo de Lade-Kim. Este

modelo é apropriado à descrição do comportamento de materiais geomecânicos

granulares. Trabalhos experimentais em paralelo procuram determinar os

parâmetros do modelo de Lade – Kim. Estes parâmetros são usados na análise

numérica.

Palavras-chave Produção de areia; Modelo elastoplástico Lade – Kim; Mecanismos de

ruptura.

Abstract

Silvestre, José Roberto; Vargas Jr., Eurípedes do Amaral (advisor); Campos, João Luiz Elias (advisor). Numerical analysis of wells with relevance for sand production. Rio de Janeiro, 2004. 133p. MSc. Dissertation - Departamento de Engenharia Civil, Pontifícia Universidade Católica do Rio de Janeiro.

This work aims to simulate situations that take a rock the rupture

considering the fluid flow, in an attempt to correlate the involved mechanisms in

the rupture to the know process as sand production in rocks reservoir, that is, solid

production during hydrocarbon extraction.

Numerical simulations using the commercial finite element program

ABAQUS, reproduce fluid-mechanical coupling, and allowing the implementation

of an elastoplastic rock material behavior.

More precisely, the author implemented a numerical routine based on the

well-known Lade-Kim elastoplastic model, which can suitably describe the

behavior of geomaterials subjected to external stress-strain fields.

Palavras-chave Sand production; Lade – Kim elastoplastic model; Failure mechanisms.

Sumário

1 Introdução 22

2 Produção de areia em poços de petróleo 25

2.1. Introdução 25

2.2. Mecanismos de produção de areia 31

3 Análise numérica de problemas a poroelasticidade 36

3.1. Introdução 36

3.2. Teoria da poroelasticidade de Biot 36

3.2.1. Equações de equilíbrio 37

3.2.2. Relação deformação – deslocamento 37

3.2.3. Relação tensão – deformação 37

3.2.4. Equações governantes 39

3.2.5. Análise de problemas de poroelasticidade pelo programa

ABAQUS 40

3.2.6. Exemplos de validação 43

3.2.6.1. Adensamento unidimensional 43

3.2.6.2. Poço vertical em um estado de tensões não hidrostático 49

3.2.6.2.1. Carregamento modo 1 50

3.2.6.2.2. Carregamento modo 2 51

3.2.6.2.3. Carregamento modo 3 51

3.2.6.3. Solução para a inversa da transformada de Laplace 56

4 Modelo constitutivo implementado 58

4.1. Critérios da Plasticidade 58

4.1.1. Critério de escoamento 59

4.1.2. Leis de evolução para endurecimento–amolecimento 59

4.1.3. Lei de Fluxo 61

4.2. Modelo Lade - Kim 62

4.2.1. Comportamento elástico 62

4.2.2. Critério de ruptura 63

4.2.3. Função potencial plástica 64

4.2.4. Função de escoamento 67

4.3. Modelo de Mohr–Coulomb 70

4.3.1. Critério de ruptura 70

4.3.2. Função de escoamento 71

4.3.3. Função potencial plástica 72

5 Algoritmo para integração da relação tensão–deformação 74

5.1. Algoritmo do tipo explícito 76

5.2. Relações usadas no processo de integração 83

5.3. Detalhes da implementação numérica 87

5.4. Comparação entre ensaios de laboratório e simulação numérica 88

6 Análise da estabilidade de poços hipotéticos 93

6.1. Poço no arenito Rio Bonito 95

6.1.1. Estimativa da produção de sólidos 105

6.2. Poço hipotético no Calcário do Campo de Congro 107

6.2.1. Permeabilidade constante 107

6.2.2. Permeabilidade variável 115

6.2.3. Comparações entre os resultados obtidos utilizando os critérios de

Mohr–Coulomb e Lade-Kim 120

7 Conclusões e sugestões para trabalhos futuros 126

8 Referências Bibliográficas 129

Lista de figuras

Figura 2.01 – Produção de areia a nível microscópico (Dusseault e

Santarelli, 1989) 26

Figura 2.02 - Esquema da completação de um poço (Fjaer e

outros,1992) 27

Figura 2.03 – Tensão vertical efetiva próxima ao poço versus pressão de

colapso em um TWC (Veeken e outros, 1991). 30

Figura 2.04 - Dano de arenito brando durante o canhoneio (Dusseault e

Santarelli, 1989) 32

Figura 2.05 - Comparação de ensaios triaxiais com a ruptura de

cavidades (Morita,1994) 33

Figura 2.06 – Diagrama simplificado da envoltória de ruptura de cavidade

proposto por Morita (1987) 34

Figura 3.01 – Esquema da coluna poroelástica 44

Figura 3.02 – Deslocamento no topo da coluna com o tempo -

62.0=α 46

Figura 3.03 – Deslocamento no topo da coluna com o tempo -

1=α 46

Figura 3.04 – Excesso de poropressão na base da coluna com o tempo -

62.0=α 47

Figura 3.05 - Excesso de poropressão na base da coluna com o tempo -

1=α 47

Figura 3.06 – Excesso de poropressão ao longo da coluna -

62.0=α 48

Figura 3.07 - Excesso de poropressão ao longo da coluna - 1=α 48

Figura 3.08 – Esquema de um poço em um meio poroelástico 50

Figura 3.09 – malha de elementos finitos utilizada na simulação do

poço. 53

Figura 3.10 – Poropressão ao longo da direção 0=θ - 62.0=α 53

Figura 3.11 – Poropressão ao longo da direção 0=θ - 1=α 53

Figura 3.12 – Tensão tangencial ao longo da direção 0=θ -

62.0=α 54

Figura 3.13 – Tensão tangencial ao longo da direção 0=θ - 1=α 54

Figura 3.14 – Tensão radial ao longo da direção 0=θ - 62.0=α 55

Figura 3.15 – Tensão radial ao longo da direção 0=θ - 1=α 55

Figura 4.01 – Representação da superfície de escoamento no plano

desviador (Desai,1984). 60

Figura 4.02 – Superfície de ruptura do modelo Lade–Kim no plano triaxial

e octaédrico (Lade e Kim, 1988). 64

Figura 4.03 – Superfície potencial plástica do modelo Lade – Kim (Lade e

Kim, 1988) 65

Figura 4.04 – Dados experimentais e contornos de trabalho plástico

constante no plano octaédrico para Fuji River Sand (Lade e Kim,

1988). 68

Figura 4.05 – Função de escoamento no endurecimento e amolecimento

(Lade e Kim, 1988). 70

Figura 4.06 – Fluxo potencial plástico pela função de Menétrey e Willam

no plano desviador (Willam e Menetrey, 1995). 73

Figura 5.01 – Interpretação gráfica da solução proposta por Euler (Boyce

e Diprima, 1998) 77

Figura 5.02 – Curva tensão versus deformação axial para a fine sílica

sand. 89

Figura 5.03 - Curva deformação volumétrica versus deformação axial para

a fine sílica sand. 89

Figura 5.04 - Curva tensão versus deformação axial para o arenito Rio

Bonito. 90

Figura 5.05 - Curva deformação volumétrica versus deformação axial para

o arenito Rio Bonito. 90

Figura 5.06 - Curva tensão versus deformação axial para o arenito Vila

Velha. 91

Figura 5.07 – Curva deformação volumétrica versus deformação axial

para o arenito Vila Velha. 91

Figura 6.01 – Representação esquemática poço vertical a ser

estudado 93

Figura 6.02 - Malhas utilizadas na simulação, a esquerda a malha 1 e a

direita a malha 2. 94

Figura 6.03 – Campos de tensão principal maior na vizinhança do poço

para os casos 1, 2 e 3 com drawdown de 2 MPa. A esquerda

representa – se a escavação e a direita a utilização do poço. 97

Figura 6.04 – Campos de tensão principal menor na vizinhança do poço

para os casos 1, 2 e 3 e drawdown de 2 MPa. A esquerda representa

–se a escavação e a direita a utilização do poço. 98

Figura 6.05 – Distribuição da tensão principal maior na vizinhança do poço

para os casos 3 e 4, ilustrando os gráficos sobre a escavação como

indicado e os demais com o poço em produção. 99

Figura 6.06 – Campos de tensão principal menor na vizinhança do poço

para os casos 3 e 4. ilustrando os gráficos sobre a escavação como

indicado e os demais com o poço em produção. 100

Figura 6.07 - Zona rompida na vizinhança do poço para os casos 1, 2 e 3

e drawdown de 2 MPa. À esquerda está representada a escavação

do poço e à direita o poço em produção. 101

Figura 6.08 - Zona rompida na vizinhança do poço para os casos 3 e 4,

ilustrando os gráficos sobre a escavação como indicado e os demais

o poço em produção. 102

Figura 6.09 – Campos de deformação volumétrica na vizinhança do poço

para o caso 1. O primeiro gráfico representa a escavação como

indicado e o seguintes o poço em produção. 103

Figura 6.10 Curva tensão principal maior-raio adimensional e tensão

principal menor-raio adimensional. 104

Figura 6.11 - Indicador dos pontos de início da produção de

sólidos. 106

Figura 6.12 – Campo de tensão principal menor na vizinhança do poço

para os casos 6 e 8 com drawdown de 5 Mpa. À esquerda

representa–se a escavação e à direita a produção do poço. 108

Figura 6.13 – Campo de tensão principal menor na vizinhança do poço

para os casos 5 e 7, ilustrando os gráficos sobre a escavação como

indicado e os demais o poço em produção. 109

Figura 6.14 - Zona de ruptura na vizinhança do poço variando o estado de

tensão e drawdown de 2 MPa. À esquerda representa –se a

perfuração e à direita a utilização do poço. 110

Figura 6.15 - Zona rompida na vizinhança do poço variando o drawdown e

um carregamento aplicado de MPaeMPa hH 25'40' == σσ . O

primeiro gráfico representa a perfuração e os seguintes a utilização

do poço. 111

Figura 6.16 – Campo de deformação volumétrica na direção θ = 0 e θ = 45

sob uma tensão horizontal maior de 40 e 45 MPa (carregamentos 6 e

8) e um drawdown de 5 MPa. 112

Figura 6.17 – Campo de deformação volumétrica na vizinhança do poço

para os casos 6 e 8. À esquerda está representada a perfuração e à

direita a utilização do poço. 113

Figura 6.18 – Campo de deformação volumétrica na vizinhança do poço

para os casos 6 e 7, ilustrando os gráficos sobre a escavação como

indicado e os demais o poço em produção. 114

Figura 6.19 - Zona rompida na vizinhança do poço para os casos 6 e 7. À

esquerda está representada a escavação e à direita o poço em

produção. 116

Figura 6.20 – Campo de deformação volumétrica na vizinhança do poço

para os casos 6 e 7, ilustrando os gráficos sobre a escavação como

indicado e o demais o poço em produção. 117

Figura 6.21 - Distribuição de poro – pressão para o caso 7 com drawdown

de 10 MPa. A primeira ilustração considera a permeabilidade

constante como indicado e a segunda a permeabilidade como função

da porosidade. 118

Figura 6.22 – Curva deformação volumétrica-raio adimensional e

permeabilidade-raio adimensional em duas direções na etapa de

produção para o calcário do Campo de Congro 119

Figura 6.23 – Campo de deformação volumétrica na vizinhança do poço

para os modelos de Mohr - Coulomb e Lade - Kim para o caso 7 com

um drawdown de 10 MPa. O gráfico a esquerda corresponde à

escavação e a direita o poço em produção. 121

Figura 6.24 - Zona de ruptura na vizinhança do poço para o caso 7 e um

drawdown de 10 MPa. O gráfico a esquerda corresponde à

escavação e a direita o poço em produção. 122

Figura 6.25 – Distribuição do índice de vazios na vizinhança do poço para

os modelos de Mohr-Coulomb e Lade-Kim para o caso 7 com um

drawdown de 10 MPa. O gráfico a esquerda corresponde à

escavação e a direita a utilização do poço. 123

Figura 6.26 - excesso de poropressão na vizinhança do poço para os

modelos de Mohr-Coulomb e Lade-Kim com um drawdown de 10

MPa, mantendo um estado de tensão de MPaeMPa hH 2540 == σσ .

O primeiro gráfico corresponde ao modelo de Lade–Kim e o último o

de Mohr – Coulomb. 124

Figura 6.27 – Simulação de ensaios uniaxial e triaxial de compressão para

o modelo de Mohr – coulom e Lade – Kim para o calcário do Campo

de Congro. 125

Lista de Tabelas

Tabela 2.01 – Parâmetros que influenciam a produção de areia 28

Tabela 3.01 – Parâmetros poroelásticos do arenito de Berea (Detournay e

Cheng, 1993) 45

Tabela 5.01 – Parâmetros do modelo Lade-Kim para a fine sílica

sand 88

Tabela 5.02 – Parâmetros do modelo Lade-Kim para o arenito Vila

Velha 88

Tabela 5.03 – Parâmetros do modelo Lade-Kim para o arenito Rio

Bonito 88

Tabela 6.01 – Carregamento no contorno e pressão aplicada no poço e

na formação para o arenito Rio Bonito. 96

Tabela 6.02 – Carregamento no contorno e pressão aplicada no poço e

na formação para o arenito Rio Bonito. 107

Tabela 6.03 – Permeabilidade em função da porosidade (Soares,

2001) 115

Lista de símbolos

a Parâmetros do critério de ruptura do modelo Lade – Kim

A Constante da função de escoamento no amolecimento do

modelo Lade – Kim

B Coeficiente de Skempton

*B Constante da função de escoamento no amolecimento do

modelo Lade – Kim

B’ Matriz de acoplamento

c Coeficiente de difusividade

*c Coesão para o modelo de Mohr – Coulomb

C Parâmetro da função de escoamento do modelo Lade – Kim

1C

C

,C

e

2

3

Constantes auxiliares da solução de Detournay e Cheng

0|c Coesão no início da plastificação

D Matriz constitutiva elástica

'D Parâmetro da função de escoamento do modelo Lade – Kim

dd Drawdown

epD Matriz constitutiva elastoplástica

1D e 2D Constantes auxiliares da solução de Detournay e Cheng

pdw Incremento de trabalho plástico

εd Incremento de deformação total

edε Incremento de deformação elástica

pdε Incremento de deformação plástica

λd Parâmetro plástico

e Parâmetro relacionado a excentricidade

E Módulo de Young

nE^

Erro local

1e Excentricidade que descreve o contorno da função potencial

plástica

f Força de massa

''f Função de escoamento em termos do trabalho plástico

F Função de escoamento

Fa Valor aproximado da função real no campo de Laplace

pf Função de escoamento em termos de tensão

g Aceleração da gravidade

G Módulo de cisalhamento

pg Função potencial plástica no modelo Lade – Kim

h Parâmetro da função de escoamento do modelo Lade – Kim

H Constante poroelástica

h’ Parâmetro de endurecimento

H’ Matriz de fluxo

1I , e 2I 3IPrimeiro, segundo e terceiro invariantes de tensão

respectivamente

J Taxa de volume do meio da sua configuração corrente para a

configuração de referência

k Permeabilidade

K Módulo volumétrico do esqueleto sólido

^k Permeabilidade do meio poroso

'K Matriz de rigidez

sk Função que introduz a dependência da permeabilidade em

relação a saturação

sK Módulo volumétrico dos grãos sólidos

Kur Parâmetros elásticos do modelo de Lade – Kim

0K Função de Bessel modificada de ordem zero

1K Função de Bessel modificada de ordem 1

L’ Matriz de acoplamento

m Parâmetro do critério de ruptura do modelo Lade – Kim

n Parâmetro elástico do modelo Lade – Kim

'n Porosidade do meio

_n Vetor normal à superfície S1

:n t Volume de fluído absorvido pelo meio poroso por unidade de

volume

p poropressão

P Parâmetro da função de escoamento do modelo Lade – Kim

'p Tensão equivalente

*p Carregamento aplicado no topo da coluna

0p Poropressão na formação

0P Tensão hidrostática

pa Pressão atmosférica

wp Pressão no poço

q Parâmetro da função de escoamento do modelo Lade – Kim

Q Função potencial plástica

'q Tensão equivalente de Mises

_q

Constante que determina a variação do incremento no

algoritmo de integração

Q’ Vetor de fluxo

pQ Constante poroelástica

r Raio de um ponto qualquer na solução de Detournay e

Cheng

R Constante poroelástica

'r Terceiro invariante de tensão desviadora

'R Erro relativo

( )ϕ,ΘmcR Medida da tensão desviadora

mwR Função desviadora elíptica

rw Raio do poço

s Saturação

S Nível de tensão

'S Constante poroelástica

dS Tensor de tensões desviadoras

0S Tensão desviadora

t Tempo

t’ Tempo adimensional

'T Variável da função real usada na inversa de Stehfest

Ts Força de superfície

u Deslocamento

V Volume de fluído

wv Velocidade média do fluído em relação à parte sólida

pw Trabalho plástico

x Coordenada espacial

x’ Profundidade adimensional

α Coeficiente de Biot

'α Constante que determina o tipo de esquema de integração

*α Parâmetro da função de escoamento do modelo Lade - Kim

β Coeficiente de velocidade

∆ Incremento relacionado a variável que o acompanha

ijδ Delta de Kronecker

wuδ Campo variacional relacionado a poro – pressão

vδ Campo de velocidade virtual

εδ Taxa de deformação virtual

eσ∆ Incremento de tensão tentativa elástica

T∆ Subincremento de tempo usado no algoritmo de integração

ε Deformação

pε Deformação plástica

pl_ε Deformação plástica equivalente

φ Carga piezométrica

η Constante poroelástica

1η Parâmetro da função de ruptura do modelo Lade – Kim

ϕ Ângulo de atrito

µ Parâmetro da função potencial plástica do modelo Lade –

Kim

ν Coeficiente de Poisson

uv Coeficiente de Poisson não – drenado

θ Ângulo do sistema de coordenada cilíndrica

Θ Ângulo polar desviador

'θ Variação do volume de fluído

wρ Massa específica do fluído

σ Tensão

nσ Tensão normal a um plano

n^σ Tensão dada pelo método de Runge – Kutta Dormmand Price

rrσ Tensão radial

θσ r Tensão cisalhante em coordenada cilíndrica

i,TWCσ Tensão de ruptura inicial de uma amostra

w,vσ Maior tensão vertical efetiva suportada pela parede de um

poço horizontal

θθσ Tensão circunferencial (ou tangencial)

τ Tensão cisalhante

ξ Parcela que representa a porção elástica durante um

incremento de tensão tentativa

ψ Ângulo de dilatância

1ψ Parâmetro da função potencial plástica no modelo Lade –

Kim

2ψ Parâmetro da função potencial plástica no modelo Lade –

Kim

1 Introdução

Há algum tempo é conhecido pela indústria do petróleo que problemas com

a estabilidade de poços podem conduzir a gastos excessivos na perfuração,

completação, redução da produtividade e até ao processo de produção de areia.

Para o caso da produção de areia vários métodos têm sido propostos a fim de

prever e selecionar dispositivos para o seu controle.

Dois processos normalmente são associados à produção de areia, o

mecânico e o químico. A ruptura do material por esforços de tração e compressão

recai no processo mecânico. A extração do material rochoso para a criação do

poço, gera um rearranjo nas tensões pré – existentes, concentrando – as na

vizinhança do poço e submetendo a rocha freqüentemente ao colapso. A rocha

antes intacta dá lugar a um material granular com pouca ou nenhuma

característica coesiva, que pode ser facilmente carreado.

Além do colapso devido à compressão ou tracionamento da rocha pela

concentração de tensão, o fluxo interfere também neste processo quer seja pela

variação de tensão que ele irá causar ou pela ação erosiva, dado que à medida que

a rocha fragmenta o fluído terá força necessária para carrear algumas partículas.

Revela – se aqui o caráter do acoplamento fluído – mecânico existente neste

problema.

Ao longo da exploração do poço, a rocha reservatório é exposta a alguns

fluídos, tais como a lama de perfuração ou a água oriunda do próprio reservatório,

caracterizando o fluxo bifásico dado o fato inevitável da depletação. Estes fluídos

em contato com a rocha podem reagir quimicamente com os constituintes

minerais dos grãos ou do material cimentante. A água age fortemente, por

exemplo, em rochas como calcário e o dolomito, destruindo o material cimentante

e até mesmo o esqueleto sólido.

A ação dos fluídos não se limita a reações químicas, a resistência adicional

oferecida pela capilaridade também é afetada, observando que intrusão de água no

rocha provoca a sua diminuição.

23

Vários métodos tem sido desenvolvidos para a previsão da produção de

areia, considerando aspectos mecânicos ou químicos. Estes métodos podem ser

baseados em resultados de campo, laboratório ou modelos teóricos. Entretanto os

dois primeiros apresentam a desvantagem de não considerar várias situações a que

um poço é submetido e por vezes serem adequados apenas a um determinado

local. Os modelos teóricos teriam, portanto uma condição mais favorável para a

análise de diversas situações.

As soluções numéricas são as mais favoráveis deste ponto de vista, pois

conseguem agregar vários eventos que podem intervir na produção de areia.

Todavia, a proximidade deste tipo solução com a realidade depende de modelos

que consigam representar com mais fidelidade um dado fenômeno.

O enrijecimento da rocha quando submetida ao confinamento, o

“amolecimento” do material após a ruptura, a magnitude das deformações sofridas

pelo material e seu comportamento dilatante que desempenham um importante

papel no fluxo de fluídos, dada a relação com a permeabilidade são características

que um modelo constitutivo deve conseguir representar com esta fidelidade.

Esta pesquisa tem como objetivo aplicar apenas para o problema mecânico

um modelo elastoplástico que seja capaz de reproduzir o comportamento

mecânico da rocha, não somente associando mecanismos que levem a rocha à

ruptura com a produção de sólidos, mas também as deformações que

desempenham um papel fundamental no fluxo de meio poroso e

conseqüentemente na produção.

O modelo elastoplástico será introduzido no programa ABAQUS através de

um algoritmo do tipo explícito que permitirá a integração da relação tensão –

deformação.

O trabalho descrito aqui está dividido em sete capítulos. O capítulo 2

caracteriza a produção de areia, o que ela é e seus principais mecanismos.

O capítulo 3 descreve a poroelasticidade e mostra algumas validações com o

programa ABAQUS para o problema de acoplamento hidro – mecânico.

O capítulo 4 faz uma breve descrição dos principais pontos na teoria da

plasticidade, assim como descreve os modelos de Lade – Kim, o qual será

implementado e sucintamente o de Mohr – Coulomb utilizado pelo ABAQUS.

O capítulo 5 trata a respeito de esquemas para integração das relações

tensão–deformação elastoplásticas. E também do algoritmo implementado no

24

programa ABAQUS. Apresenta – se comparações entre a simulação de ensaios

triaxiais e resultado de ensaios de laboratório.

O capítulo 6 apresenta os resultados obtidos para a simulação de um poço

submetido a diferentes carregamentos.

2 Produção de areia em poços de petróleo

2.1 Introdução

A produção de partículas durante a extração de fluídos de uma rocha

reservatório é designada por produção de areia. A produção de areia ou produção

de sólidos referenciada normalmente a arenitos não consolidados, engloba além

destes materiais, rochas como calcário e arenitos de resistência média a elevada,

também susceptíveis a este fenômeno. Os danos mais comuns devido à produção

de areia incluem abrasão e desgaste de equipamentos, perda da produtividade do

poço, colapso de revestimento e disposição dos resíduos impregnados por

hidrocarbonetos.

Em uma escala microscópica, a produção de areia é um processo

desencadeado pela força de arraste associada ao fluído em movimento atuando em

um conjunto de partículas próximas a uma superfície livre (parede do poço,

cavidade do canhoneado). O transporte de partículas ocorre, quando esta força,

originada na diferença de pressão estabelecida entre o poço e a formação,

eventualmente se torna maior do que as forças estabilizadoras provenientes da

coesão do material e da tensão de arco desenvolvida entre as partículas, como

ilustrado na figura 2.01.

No intuito de compreender como a perfuração e a completação do poço

influem na produção de areia, descreve-se sucintamente tais etapas.

A perfuração de um poço de petróleo realizada através de uma sonda do

tipo rotativa, consiste no avanço obtido pela rotação e peso aplicado por uma

broca existente na extremidade de uma coluna de perfuração. Os fragmentos da

rocha perfurada são continuamente removidos por um fluido de perfuração. Este

fluido tem a função de limpar o fundo do poço dos fragmentos produzidos,

estabilizar a parede do poço, resfriar e lubrificar a coluna de perfuração e a broca.

Atingida a profundidade do reservatório, retira-se à coluna de perfuração e uma

coluna de revestimento de aço é colocada.

26

Figura 2.01 – Produção de areia a nível microscópico (Dusseault e Santarelli, 1989)

Concluída a perfuração, seguem-se às atividades de revestimento,

cimentação e caso necessário o canhoneio. A esta etapa designa-se o nome de

completação, ou seja, o conjunto de atividades destinadas a equipar um poço para

a produção de óleo ou gás.

Na figura 2.02, observa-se à união do revestimento à formação feita

através de um cimento e as perfurações (canhoneio) existentes através das quais

ocorrerá a produção de óleo ou gás. Essas perfurações realizadas com o uso de

pistolas se estendem desde o revestimento até a formação, suas formas dependem

do tipo e quantidade de carga utilizada, apresentando um aspecto fino e alongado

inicialmente. À medida que o poço é explorado, elas crescem e tornam - se

cavidades com possibilidade de virem a fundir. O crescimento destas perfurações

implica na perda de material sólido da formação, o qual é carreado pelo fluxo de

fluido, desencadeando o processo de produção de areia. Uma forma conveniente

de verificar a alteração da geometria destas cavidades é monitorar o volume

acumulado de areia produzido em um determinado intervalo.

Veeken (1991) sugere a classificação da produção de areia para uma

melhor interpretação e comparação dos eventos que levam um poço a produzir

27

sólidos. Esta classificação baseada em observações de campo permite a distinção

de três modos:

- Produção de areia transiente: refere-se ao declínio da produção de areia

com o tempo sob condições de produção constante do poço. Este tipo de

fenômeno é freqüentemente observado durante a limpeza da parede do poço após

a perfuração das cavidades.

- Produção de areia contínua: na maioria dos campos de exploração,

observa-se a produção contínua de areia. A limitação da produção de areia deve

obedecer a restrições operacionais como capacidade de separação, disposição de

areia, localização do poço, etc. Limites típicos são da ordem de 6 a 600 g/m3 para

formações de óleo e de 16 g/103 m3 para formações de gás. Dependendo da

capacidade do fluido em transportar partículas e da concentração, eventualmente

ocorrerá a obstrução de parte do intervalo produtor.

- Produção de areia catastrófica: causa a suspensão das operações de um

poço, pode ocorrer devido a uma produção maciça de areia que preenche e obstrui

o poço.

revestimento

canhoneadoformação

Figura 2.02 - Esquema da completação de um poço (Fjaer e outros, 1992)

Diante da necessidade de se minimizar a níveis toleráveis a produção de

areia, um mecanismo de controle deverá ser utilizado. Estes mecanismos atuam

diretamente nos parâmetros de produção do poço ou diretamente na produção de

areia. O controle feito sobre os parâmetros de produção consiste em manipular o

drawdown (diferença entre a pressão no poço e a poropressão no reservatório) e a

taxa de produção. O segundo tipo de controle é efetuado por técnicas como:

- Gravel packing;

- Consolidação química;

- Pré–consolidação da formação;

- Seleção das perfurações;

28

- Frac–Pack.

Destas técnicas, provavelmente a mais utilizada é o gravel packing, que

consiste na utilização de um filtro sólido granular, para impedir a produção de

grãos na formação. A consolidação química promove um aumento na resistência

da formação pela injeção de resina, todavia leva a uma redução da permeabilidade

e conseqüentemente da produção de um poço.

A seleção dos locais de perfuração é um outro meio para minimizar a

produção de areia, a idéia consiste em perfurar regiões onde a formação possui

uma boa resistência mecânica, e com isto uma pequena possibilidade de produzir

areia.

O frac-pack é uma técnica adequada para reservatórios onde a taxa de

produção de óleo é elevada, cuja convergência do fluxo radial pode se tornar um

agravante para a produção de sólidos. O método consiste na criação de uma

pequena fratura condutora, cuja função é transformar o fluxo radial em fluxo

linear através do poço, reduzindo o gradiente de pressão.

Apesar de controlarem a produção de areia, estas técnicas apresentam

como inconveniente o alto custo de instalação e manutenção, assim como a

redução da produção do poço. A escolha de qual técnica e quando utilizá-la torna

necessária a introdução de um método para a previsão da produção de areia.

As técnicas de previsão de produção de areia existentes são baseadas em

observações de campo, ensaios de laboratório ou modelos teóricos (Morita, 1989).

A previsão que utiliza a observação de campo consiste no estabelecimento de uma

correlação entre dados de produção de areia de um poço e parâmetros

operacionais. A precisão deste método é proporcionalmente dependente da

quantidade de parâmetros utilizados. Dada à dificuldade de monitorar e armazenar

parâmetros, poucos são incluídos em uma análise. A tabela 2.01 apresenta alguns

parâmetros que influenciam a produção de areia.

Formação Completação Produção

Resistência da rocha Diâmetro e orientação do

poço Taxa de fluxo

Tensões vertical e horizontal

in – situ

Perfuração – tamanho, fase,

densidade Drawdown

Heterogeneidade Dimensão dos tubos Depletação

Tabela 2.01 – Parâmetros que influenciam a produção de areia

29

Ensaios de laboratório permitem reproduzir a produção de areia em

condições controladas e desta forma observar separadamente a influência de cada

parâmetro. Além disto, possibilitam a validação de um modelo teórico e também

podem ser utilizados como ferramentas de previsão. Em seu trabalho,

Veeken(1991) cita a utilização do Thick Walled Cylinder (TWC), relacionando a

maior tensão vertical efetiva suportada pela parede de uma perfuração horizontal

( wv , )σ com a tensão de ruptura inicial de uma amostra ( )i,TWCσ , através da

expressão:

i,TWCw,v σ=σ (2.01)

a tensão vertical efetiva é definida empiricamente como: ddvwv += σσ , (2.02)

onde é o drawdown. Como a tensão inicial de ruptura não é facilmente

identificada, estabeleceu-se uma relação com a tensão de colapso

dd

( )TWCσ baseada

em uma série de ensaios em arenitos friáveis e consolidados, .

TWCw,v 86.0 σ×=σ (2.03)

Na figura 2.03 valores fornecidos pela equação acima são comparados

com dados de produção de areia obtidos em campo. Nota-se os valores

conservativos quanto ao início da ruptura, mas que podem ser utilizados com certa

confiança.

Modelos teóricos agrupam-se em soluções fechadas ou numéricas. Apesar

da solução fechada ter um uso atraente, por muitas vezes elas simplificam os

parâmetros envolvidos na produção de areia, como por exemplo, a trajetória de

tensão, geometria do poço e canhoneado. A solução numérica torna-se uma

valiosa ferramenta para a simulação deste processo.

A confiabilidade da solução numérica na representação dos mecanismos

de ruptura envolvidos na produção de areia dependerá principalmente do modelo

constitutivo adotado. Mclean (1991) discute sobre modelos que consideram o

caráter tridimensional do critério de ruptura e a natureza não linear do material,

que são fatores importantes na caracterização do comportamento volumétrico e de

endurecimento do material geológico.

30

Figura 2.03 – Tensão vertical efetiva próxima ao poço versus pressão de colapso em um

TWC (Veeken e outros, 1991).

Morita (1987) cita como obstáculos à utilização de códigos de elementos

finitos comerciais a pouca flexibilidade para representar o fluxo de óleo através de

um meio poroso, onde suas propriedades são complexas e por muitas vezes o

fluxo não é regido pela lei de Darcy.

Segundo Dusseault e Santarelli (1989), um modelo teórico deverá prever a

produção de areia tanto a nível microscópico partículas, fenômenos de arcos de

areia, desintegração de grãos, efeitos da força viscosa do fluido e variações da

resistência proporcionada pela coesão, assim como a nível macroscópico a análise

da plastificação e do caráter viscoso do fluido como funções do campo de tensões

e gradientes de pressão. O modelo a ser utilizado deverá considerar a envoltória

de plastificação e critério de ruptura do material, devendo ser validado com

ensaios de laboratório ou dados de campo.

A implementação destes modelos na produção de areia é caracterizada por

um alto grau de incerteza devido ao envolvimento de vários parâmetros

relacionados a ela. Assim como, a redução de parâmetros a serem utilizados

conduz a uma simplificação do problema e a desconfiança dos resultados.

31

2.2 Mecanismos de produção de areia

A produção de areia durante a exploração de hidrocarbonetos pode ser

relacionada a dois mecanismos fundamentais de desestabilização. O primeiro

relacionado ao comportamento mecânico, que leva a ruptura e a plastificação

localizada da rocha. O segundo está associado ao fluxo, que pode conduzir a

instabilidades hidrodinâmicas como o arraste de partículas.

O comportamento mecânico é influenciado pela ruptura por compressão e

por tração. A rocha fragmentada pela compressão, devido à concentração de

tensão na vizinhança do poço, disponibiliza partículas que podem ser “arracandas”

pela força de percolação do fluido, caracterizando a produção de areia por ruptura

de tração. Este tipo de produção de areia é comum em arenitos pouco

consolidados, produzindo baixa quantidade de areia e em geral de forma

esporádica (Morita e Boyd, 1991).

O mecanismo hidrodinâmico inicia-se com o cisalhamento do cimento que

une os grãos da rocha. Os grãos desagregados pela desintegração do cimento

mineral ficam susceptíveis ao carreamento pelas forças de percolação do fluido. A

remoção de partículas aumenta a cavidade do poço, redistribuindo as tensões, o

que favorece o primeiro processo de desestabilização. Observa-se nesta situação a

interação dos dois mecanismos.

A produção de sólidos normalmente não é decorrente apenas de um

evento, mas de uma seqüência. Portanto a história da formação deve ser conhecida

desde a perfuração até o momento em que o poço é utilizado. Dusseault e

Santarelli (1989) propuseram um modelo conceitual para a produção de sólidos.

Baseado neste modelo ilustra-se alguns mecanismos envolvidos neste fenômeno.

Anteriormente citou-se a perfuração de um poço como uma atividade a ser

conhecida, sua influência na produção de sólidos está relacionada à estabilização

das paredes do poço pelo fluido de perfuração. A utilização de um fluido não

penetrante, colocado a uma pressão maior do que a poropressão na formação,

confere à rocha na parede do poço um aumento de resistência, o que normalmente

impede a plastificação. Entretanto, se houver plastificação na parede, a perfuração

eventualmente origina um material granular, facilmente removido pelo fluido de

32

perfuração. Em arenitos de média a alta resistência, o cimento mineral pode ser

destruído gerando um material de resistência reduzida e usualmente granular.

Na completação, muitos cimentos utilizados para unir revestimento e

formação retraem, possibilitando a tensão radial na parede do poço ser menor do

que na etapa de perfuração. O descarregamento gerado propicia a formação de um

material granular, iniciada à produção de hidrocarbonetos, este material é

carreado. A cavidade produzida pelo canhoneamento apresenta regiões de coesão

reduzida e grãos triturados, o que leva mais sólidos a serem carreados. A figura

2.04 ilustra esta situação.

Figura 2.04 - Dano de arenito brando durante o canhoneio (Dusseault e Santarelli, 1989)

Uma analogia entre a região perfurada pelo canhoneio e um modelo

elastoplástico permite visualizar o comportamento da rocha. Envolvendo a

cavidade, existe um material granular comportando-se de forma totalmente

plástica. Próximo à cavidade, o arenito se comporta de forma elastoplástica, sob

regimes de “endurecimento” ou “amolecimento”. Afastando-se da cavidade, a

rocha se comporta de forma elástica, praticamente sem danos.

A geometria da perfuração é outro item a ser considerado. Há vinte anos

quando o sistema de perfuração era menos eficiente, cavidades de pequena

dimensão e uma baixa densidade de perfurações no intervalo produtor levavam a

problemas relacionados à ruptura por tração (Morita, 1987).

33

Durante a vida produtiva de um poço, a migração de finos e o

tamponamento de poros podem conduzir a redução da permeabilidade da rocha,

aumentando desta forma a concentração de fluxo e conseqüentemente as forças

viscosas do fluido, levando rochas com baixa resistência a uma ruptura por tração.

Morita correlacionou a ruptura das cavidades de perfuração a um ensaio triaxial, a

figura 2.05 ilustra a situação.

Figura 2.05 - Comparação de ensaios triaxiais com a ruptura de cavidades (Morita,1994)

Sob um estado de tensão desviador pequeno, a cavidade apresentará

ruptura apenas se houver um grande fluxo, pois provavelmente ela não se

plastificou (ponto A). À medida que a tensão desviadora aumenta, a rocha

plastifica. Neste caso, duas situações podem vir a se desenvolver. A primeira

caracteriza o acréscimo de carga até o ponto de ruptura (C) seguido de um

descarregamento de um descarregamento (ponto B), desenvolvendo carregamento

fissuras na rocha de uma forma controlada, em geral paralela à direção da tensão

desviadora. Essas fissuras possibilitam o material, antes “integro”, ser carreado

por um pequeno fluxo.

A outra situação é o carregamento além do pico de ruptura (ponto D), onde

as fissuras já não se desenvolveriam de forma controlada. O desenvolvimento

descontrolado de fissuras possibilita o surgimento de fraturas, desagregando e

reduzindo a resistência do material. O cenário referente ao descarregamento pode

ser interpretado como os disparos do canhoneio em que a rocha é solicitada

rapidamente, atingido o estado de ruptura.

Uma envoltória de ruptura para a cavidade é proposta por Morita (1989)

considerando os dois tipos de ruptura, como ilustrado na figura 2.06. Os

34

parâmetros de entrada nesta envoltória são o drawdown e o gradiente de pressão

normalizado, com este diagrama é possível estabelecer uma região segura de

operação do poço. Nota-se na figura que para um estado de tensão desviador

pequeno (proporcionado por baixos valores de drawdon), a cavidade deixará de

ser estável, caso aplique-se um alto gradiente de pressão. Aumentando-se o

drawdon, a ruptura por cisalhamento é favorecida e pequenos valores de gradiente

de poropressão são necessários para que ocorra a instabilidade da cavidade.

drawdown

Gradiente de poro-pressão

Ruptura por tração

Ruptura por cisalhamento

Operação segura

Figura 2.06 – Diagrama simplificado da envoltória de ruptura de cavidade proposto por

Morita (1987)

A depletação do reservatório é outro fator que influencia a produção de

areia, pelo simples fato de que a extração de hidrocarbonetos acarretará na

diminuição da poropressão na formação. Esta mudança na poropressão leva a um

aumento na tensão efetiva, proporcionando condições favoráveis a ruptura por

cisalhamento. Um procedimento para se evitar esta situação e também com o

objetivo de aumentar a produção é a injeção de água no reservatório.

Na literatura, relata-se o início ou aumento na produção de areia logo após

o water–cut, que é um fenômeno decorrente do processo de injeção. Algumas

causas ligadas ao water–cut que influenciam a produção de sólidos são

enumeradas abaixo:

- perda da pressão capilar que auxilia a união dos grãos ocasionada pela

injeção de água no sistema.

- baixa poropressão existente na formação no início do water-cut na época

da injeção, o que teria tornado favorável o desenvolvimento de tensões de

cisalhamento no reservatório.

35

- a tentativa de se manter os níveis de produção antes da injeção induz a

um aumento na taxa de fluxo, o que leva a uma redução na pressão do poço e a

um aumento no gradiente de pressão na superfície da cavidade.

- migração de finos devido ao alto fluxo imposto leva a uma maior taxa de

tamponamento.

- a possibilidade de fusão entre as cavidades aumenta devido ao seu

crescimento.

- ação química da água de injeção em determinados cimentos.

Compreende-se que as diferentes etapas da vida de um poço como a

perfuração, completação e produção influem de diversas formas na produção de

sólidos, que é vista então, não como um processo estacionário, mas evolutivo ao

longo de uma seqüência de eventos.

Atividades como a perfuração, concentram tensões ao redor do poço,

conduzindo por muitas vezes a rocha a ruptura por compressão.Assim como, o

aumento na taxa de produção ou a depletação do reservatório conduzem a mesma

situação. Em decorrência do aumento de fluxo de fluído, o gradiente de

poropressão próximo à parede do poço, favorecendo a ruptura por tração. Esta

pesquisa visa, então, analisar os mecanismos de ruptura por compressão e tração

da rocha reservatório descritos acima e associá-los a produção de areia.

3 Análise numérica de problemas a poroelasticidade

3.1 Introdução

Qualquer tipo de escavação em uma rocha leva a um descarregamento das

tensões pré-existentes, quer seja a abertura de túneis, shafts para minas, poços e

outros. Na exploração de poços, a variação da pressão do fluído que preenche os

poros da rocha, interage com o campo de tensões, à medida que a poropressão

varia ao redor do poço contribuindo de modo significativo no comportamento

mecânico.O interesse deste trabalho é associar os mecanismos de ruptura

decorrentes da alteração do estado de tensões in situ ao redor de um poço de

petróleo com processos de produção de areia.

Esta dependência do processo difusivo e relação tensão-deformação levam à

consideração do acoplamento fluido–mecânico no estudo da produção de areia.

Neste capítulo descreve-se sucintamente o processo de acoplamento baseado na

teoria poroelasticidade de Biot (1941). Além deste tópico, encontra-se uma

discussão sobre a solução adotada pelo programa de elementos finitos ABAQUS e

a validação de resultados obtidos pelo programa através da solução analítica

proposta por (Detournay e Cheng, 1988).

3.2 Teoria da poroelasticidade de Biot

O meio poroso descrito em sua teoria é elástico linear, isotrópico e

considera os poros totalmente ocupados por um fluído.A síntese feita neste tópico

basea - se no trabalho de Biot (1941).

As equações governantes do problema poroelástico provêem das equações

de equilíbrio, compatibilidade deformação–deslocamento, relação tensão-

deformação da teoria da elasticidade e da lei de Darcy.

37

3.2.1 Equações de equilíbrio

Estas equações são as mesmas dadas pela teoria da elasticidade, entretanto

considera-se que a tensão normal (σ ) a um plano é constituída por duas parcelas,

uma representando a tensão no esqueleto sólido e outra a poropressão no fluído, f

é a força de massa.

0fij,ij =+σ (3.01)

Neste tipo de notação, índices repetidos significam soma e o sinal de vírgula

a derivada.

3.2.2 Relação deformação – deslocamento

Estabelecem uma função entre deformação e deslocamento.

( )i,jj,iij uu21

+=ε (3.02)

onde ε é a deformação e u representa o deslocamento.

3.2.3 Relação tensão – deformação

Biot inclui na sua formulação uma variável adicional para descrever a

quantidade de fluído que ocupa os vazios do meio poroso. Esta variável representa

o incremento de volume de fluido por unidade de volume de material e é

designada por variação do volume de fluído (θ ’), que está relacionada a poro-

pressão designada por p.

A relação entre deformação e tensão é expressa por:

H3p

D ijkl

1ijklij

δ+σ=ε − (3.03)

onde D e matriz constitutiva do esqueleto sólido, descrita pelas constantes

elásticas do módulo de cisalhamento (G), módulo de Young (E) e o coeficiente de

Poisson υ ; H é uma constante relacionada ao fluído e δ é o delta de Kronecker.

A parcela Hp

3 introduz a poropressão e é adicionada somente as tensões

normais por não produzir qualquer tensão cisalhante. Seu efeito é igual nas três

38

componentes de deformação normal ao plano de referência devido à hipótese de

isotropia. A expressão (3.03) relaciona deformação à tensão e poropressão,

entretanto necessita-se de uma relação entre a variação do volume de fluído com a

tensão e a poropressão. Uma relação geral é dada por:

cpai += σθ ' (3.04)

sendo c e constantes, a σ o vetor de tensões e i um índice variando de 1 a 6.

Considerando novamente a hipótese de isotropia, uma mudança de sinais

na tensão cisalhante não deverá influenciar na variação do volume de fluído,

portanto as constantes que estão multiplicadas pelas tensões cisalhantes terão

valores nulos, o que permite reduzir o intervalo i de 1 a 3.

Biot desenvolve a expressão (3.04) na seguinte forma:

( ) pRH zyx1

31'

1+++= σσσθ (3.05)

Biot na sua formulação considera a existência de uma energia potencial, que

permite estabelecer a igualdade entre e (expressão 3.03). A partir de então

somente será referenciado. Expressando a tensão em função da deformação

através da expressão (3.03), obtêm-se:

1H H

H

pG ijkk

ijij αδν

νεεσ −

−

+=21

2 (3.06)

onde:

( )( ) H

Gυυα213

12−+

= (3.07)

Somando as três componentes de tensão normal que podem ser obtidas pela

expressão (3.06) e introduzindo essa soma na expressão (3.07), tem–se para a

expressão (3.05) a seguinte forma:

pkk Q

p+=αεθ ' (3.08)

onde

HRQp

α−=

11 (3.09)

As constantes elásticas presentes na expressão (3.03) estão relacionadas ao

esqueleto sólido. As constantes relacionadas ao fluído, presentes nas expressões

(3.08) e (3.09) podem ser interpretadas através de um simples exemplo de uma

39

amostra de solo envolvida por uma fina membrana, tal que as tensões aplicadas

sejam desprezíveis. Introduzindo um pequeno tubo que acompanhará a membrana

e submetendo uma poropressão negativa p, uma certa quantidade de água é

drenada. Pela expressão (3.05) tem-se:

pR1' −=θ (3.10)

Somando as três componentes normais de deformação dadas pela expressão

(3.03), a deformação volumétrica é dada por:

Hp

kk −=ε (3.11)

Logo, as constantes H1 e

R1 representam respectivamente a

compressibilidade do solo para uma variação na poropressão e a mudança no

volume de água para uma dada mudança de poropressão.

3.2.4 Equações governantes

O acoplamento fluído-mecânico é um processo transiente, procede-se a

seguir a caracterização das expressões descritas anteriormente como função do

tempo.

As tensões dadas na expressão (3.06) devem satisfazer a equação de

equilíbrio (3.01), utilizando a relação deformação-deslocamento (3.02) tem-se:

021

2 =∂∂

−∂∂

−+∇

ii

kki x

px

GuG αε

υ (3.12)

A expressão 3.12 descreve o comportamento mecânico do meio poroso,

observa-se a semelhança com a equação de Navier estudada na elasticidade.

Entretanto, ainda é necessária uma relação para o comportamento difusivo, esta

relação será obtida do balanço de massa. Considerando-se que um fluído

incompressível atravesse um cubo de dimensões infinitesimais, a taxa de fluído

que atravessa uma área unitária em um tempo t deverá ser igual à variação de

volume de fluido no cubo no mesmo tempo t.

O volume de fluído que atravessa o cubo é dado pela lei de Darcy

40

ii x

pkV∂∂

−= (3.13)

tem-se:

ii

ii

xV

t ∂∂

−=∂∂ 'θ (3.14)

Substituindo a expressão (3.13) e (3.14) em (3.08) obtêm-se:

tp

Qtpk

p

kk

∂∂

+∂∂

=∇12 ε

α (3.15)

As equações (3.12) e (3.15) formam o conjunto de equações governantes da

poroelasticidade. Apesar da formulação proposta por Biot explicar fenômenos da

área geomecânica como a subsidência devido à drenagem de um fluído ou a

ruptura por tração induzida pela pressurização de um poço, ela apresenta o

inconveniente de seus parâmetros não permitirem uma fácil interpretação física.

Rice e Cleary (1976) colocaram a formulação de Biot em função de parâmetros

usuais da mecânica dos solos e das rochas. Uma descrição mais aprofundada desta

formulação e outros trabalhos como o de Risnes (1992) é feita por Ferreira (1996).

3.2.5 Análise de problemas de poroelasticidade pelo programa ABAQUS

A complexidade das equações governantes da poroelasticidade torna a

geração de soluções analíticas uma tarefa difícil. A técnica numérica é então o

meio mais apropriado para a obtenção de resultados. O programa ABAQUS foi

selecionado pela sua potencialidade em resolver problemas diversos de

engenharia, no presente caso, o acoplamento fluido-mecânico.

O acoplamento fluído–mecânico como visto, consiste na solução de um

sistema de equações diferenciais de equilíbrio e balanço de massa de um meio

poroso. De acordo com o manual do usuário (ABAQUS - Theory Manual), o meio

poroso é considerado no programa como um meio multi-fásico constituído por

matéria sólida e seus vazios preenchidos por um líquido e um gás. As equações de

equilíbrio e continuidade discretizadas são solucionadas através do método de

Newton.

41

A condição de equilíbrio é expressa pelo princípio do trabalho virtual para

um determinado volume em um tempo t qualquer como:

( )∫ ∫ ∫∫ ⋅++⋅+⋅=v S V vwtV vvs dVgnsndVfdSTdV δρδδδεσ ': (3.16)

sendo:

:vδ campo de velocidade virtual

:εδ taxa de deformação virtual

:sT força de superfície

:f força de massa

:σ tensão

:wρ massa específica do fluído

:g aceleração da gravidade

:n t volume de fluído absorvido pelo meio poroso por unidade de volume

:'n porosidade do meio

s: saturação

A parcela referente à absorção fluído pelo sólido não será considerada neste

trabalho. O balanço de massa é obtido da mesma forma do item 3.2.4, entretanto o

fluxo é regido pela lei de Forchheimer, que é dada por:

( )x

kvvvsn www ∂∂

−=+φβ

^1' (3.17)

onde:

:vw velocidade média do fluido em relação à parte sólida (velocidade de

percolação)

:β coeficiente de velocidade

:k^

permeabilidade do meio poroso

:φ carga piezométrica

Nota-se que se a velocidade do fluído for baixa a lei de Forchheimer se

reduz à lei de Darcy, a mesma condição ocorre se o coeficiente de velocidade for

nulo. Assim, o balanço de massa usando a lei de Forchheimer é:

42

( ) ( )

( )∫

∫ =∆+

−∂∂

∂∂

+∆

+

+−+

V

S ww

ww

www

www

s

tt

w

wt

w

ww

dSvnsnut

dV

gx

uk

xu

vvgk

t

nsnJJ

nsnu

1

0'

1

'1'

1

_

0

0

00

ρρ

δ

ρδ

βρ

ρρ

ρρ

δ

(3.18)

sendo:

:u wδ campo variacional relacionado a poro-pressão

:0wρ massa específica em uma configuração de referência

0dVdVJ = : taxa de volume do meio da sua configuração corrente para a

configuração de referência

:t∆ incremento de tempo

:k s função em termo de saturação que introduz uma dependência da

permeabilidade em função da saturação

( ) :,exk permeabilidade em função da coordenada espacial e do índice de

vazios.

:n_

vetor unitário normal à superfície S1

Na equação (3.18) a permeabilidade k ficou definida pelo produto da

função por , onde para um meio saturado é igual a 1. Segundo o manual

de teoria do ABAQUS as equações governantes do processo difusão de fluído e

deformação são:

^

sk k sk

equação de equilíbrio: FpLuK =+ ''

equação de fluxo: ''' QpHdtduB =+

onde K é a matriz de rigidez, 'H é matriz de fluxo, e 'L 'B são matrizes

permitem o acoplamento.

Existem dois tipos de aproximação para resolver este sistema de equações.

Uma aproximação seria solucionar primeiro um conjunto de equações, depois com

o resultado da primeira equação resolver a segunda. Com o resultado da segunda

equação retorna-se para a primeira e verifica-se a variação de resultados. Se a

43

variação é mínima, a solução por este processo iterativo termina. Caso contrário, o

processo iterativo continua até que a variação de resultados seja mínima. O

segundo tipo de solução é resolver as duas equações ao mesmo tempo, este modo

é o adotado pelo programa ABAQUS. A validação do programa é feita no

próximo tópico.

3.2.6 Exemplos de validação

A validação do programa ABAQUS será feita com dois exemplos. O

primeiro se refere a uma coluna poroelástica e o segundo a um poço vertical

submetido a um estado de tensões não hidrostático.

3.2.6.1 Adensamento unidimensional

A situação em questão é a de um horizonte de solo de espessura L

repousando sobre uma camada rígida e impermeável, cuja superfície sofre a

aplicação de um carregamento sob condições drenadas. As condições de contorno

são e p = 0 em x = 0, para x = L tem-se )(* tHpxx −=σ 0=xu e 0=∂∂

xp .

Desde que o carregamento seja constante, a solução para este problema recai

na equação de difusão.

02

2=

∂∂

−∂∂

xpc

tp (3.19)

onde c é o coeficiente de difusividade, dado por

( )( )( ) ( )u

u

vvvvvkG

c−−

−−=

12112

22α (3.20)

44

Figura 3.01 – Esquema da coluna poroelástica

k é a permeabilidade, G módulo de cisalhamento, coeficiente de Poisson não

drenado, coeficiente de Poisson e

uv

v α é o coeficiente de Biot, definido aqui como

função do módulo volumétrico do esqueleto sólido ( )K e dos grãos ( )sK .

sKK

−= 1α (3.21)

A solução analítica para o caso da coluna poroelástica é dada por Detournay

e Cheng (1993). A poropressão e o deslocamento são dados em função das

variáveis adimensionais x’ e t’. A expressão do excesso de poropressão é dada

pela expressão (3.22).

( )','(1* txFGSpp −=

η (3.22)

onde:

Lxx =' (3.23)

24'

Lctt = (3.24)

( )∑∞

=

−

−=

,..3,1

22 'exp2

'sen41)','(m

tmxmm

txF πππ

(3.25)

( )vv

−−

=12

21αη (3.26)

u

u

vv

GBS

+−

=113' η (3.27)

45

( )( ) KKK

Bf

f

φαφαα

+−−=

1 (3.28)

sendo módulo volumétrico do fluído, fK φ a porosidade, B o coeficiente de

Skempton e p* o carregamento. A expressão para o deslocamento é dada por

uuu uxx ∆+= (3.29)

sendo:

( )( ) ( )'112

21*x

vGvLp

uu

uux −

−−

= (3.30)

( )( )( ) ( )','

112*

2 txFvvGvvLp

uu

ux −−

−=∆ (3.31)

( )[ ]∑∞

=

−−

=

,...3,1

22222 'exp1

2cos8

mtmxm

mF ππ

π (3.32)

onde é o deslocamento inicial elástico na condição não- drenada e é o

incremento de deslocamento dependente do tempo.

uxu xu∆

Apresenta-se a seguir a comparação de resultados da simulação numérica

com a solução analítica para o arenito de Berea, considerando a situação dos

materiais constituintes serem incompressíveis ou não. Esta condição é

caracterizada pelo coeficiente de Biot ( )α , que para a situação incompressível

assume o valor unitário. O carregamento p é igual a 1 MPa e os parâmetros

referentes ao material são listados na tabela 3.01. G 6000 MPa

υ 0.20

uυ 0.33

B 0.62

c 1.6 m2/s η 0.30

α 0.79

sK 36000 MPa

fK 3300 MPa

'n 0.19

k 1.9 x 102 mD

Tabela 3.01 – Parâmetros poroelásticos do arenito de Berea (Detournay e Cheng, 1993).

46

Por uma imposição do programa de elementos finitos, a reposta da pressão

de fluido é definida em excesso de poropressão. A comparação da resposta

analítica e numérica será feita obedecendo à convenção do programa.

4.00E-04

4.50E-04

5.00E-04

5.50E-04

6.00E-04

6.50E-04

0 100 200 300 400 500 600 700 800Tempo (s)

Des

loca

men

to (m

) Detournay e ChengAbaqus

Figura 3.02 – Deslocamento no topo da coluna com o tempo - 62.0=α

0.00E+00

1.00E-04

2.00E-04

3.00E-04

4.00E-04

5.00E-04

6.00E-04

7.00E-04

0 100 200 300 400 500 600 700 800Tempo (s)

Des

loca

men

to (m

)

Detournay e ChengAbaqus

Figura 3.03 – Deslocamento no topo da coluna com o tempo - 1=α

47

0

0.05

0.1

0.15

0.2

0.25

0.3

0.35

0.4

0.45

0 100 200 300 400 500 600 700 800Tempo (s)

Exce

sso

de p

orop

ress

ão (M

Pa)

Detournay e ChengAbaqus

Figura 3.04 – Excesso de poropressão na base da coluna com o tempo - 62.0=α

0.00

0.20

0.40

0.60

0.80

1.00

1.20

0 100 200 300 400 500 600 700 800Tempo (s)

Exce

sso

de p

oro-

pres

são

(MPa

)

Detournay e ChengAbaqus

Figura 3.05 - Excesso de poropressão na base da coluna com o tempo - 1=α

48

0

0.05

0.1

0.15

0.2

0.25

0.3

0.35

0.4

0.45

0.5

0.00 2.00 4.00 6.00 8.00 10.00 12.00Posição na coluna (m)

Exce

sso

de p

oro-

pres

são

(MPa

)

Detournay e ChengAbaqus

t = 507 s

t = 307 s

t = 107 s

t = 10 s

Figura 3.06 – Excesso de poropressão ao longo da coluna - 62.0=α

0.00

0.20

0.40

0.60

0.80

1.00

1.20

0.00 2.00 4.00 6.00 8.00 10.00 12.00Posição na coluna (m)

Exce

sso

de p

oro-

pres

são

(MPa

)

Detournay e Cheng

Abaqus

t = 527 s

t = 327 s

t = 127 s

t = 27 s

Figura 3.07 - Excesso de poropressão ao longo da coluna - 1=α

As figuras 3.02 e 3.03 mostram o deslocamento do topo da coluna com o

tempo, a resposta obtida pela simulação numérica possui boa concordância com a

solução analítica. É interessante notar também o deslocamento no início do

processo de adensamento, segundo a mecânica dos solos para constituintes

49

incompressíveis ( 1= )α o processo de adensamento não produz qualquer

deslocamento (situação não – drenada), ao passo que para materiais compressíveis

( 62.0= )α isto não se verifica.

Em relação às figuras 3.04 e 3.05 a concordância entre os resultados se

repete. Como no caso do deslocamento há uma diferença de comportamento entre

materiais compressíveis ou não. De acordo com a teoria de Terzaghi (1923), para

a situação de materiais incompressíveis, no instante do carregamento o fluído

absorve toda a carga aplicada ao meio poroso, como consequência o excesso de

poropressão é igual ao carregamento. O resultado mostrado na figura 3.05 está de

acordo com a teoria. Entretanto, isto não ocorre quando os materiais são

compressíveis (figura 3.04), onde parte da carga é suportada pelo esqueleto sólido

no instante do carregamento, este fato se reflete no deslocamento do topo da

coluna poroelástica como citado anteriormente.

As figuras 3.06 e 3.07 mostram o excesso de poropressão ao longo da

coluna onde também houve boa concordância de resultados. Este exemplo mostra

que o programa ABAQUS possui uma boa capacidade em representar o

comportamento de materiais compressíveis ou não.

3.2.6.2 Poço vertical em um estado de tensões não hidrostático

Detournay e Cheng (1989) propuseram a solução analítica para a escavação

de um poço vertical em uma formação saturada sujeita a um estado de tensões não

hidrostático. A solução analítica para este problema é originada no campo de

Laplace, assumindo um estado de deformação plana no plano perpendicular ao

eixo do poço e escavação instantânea; a solução da transformada é feita

numericamente.

A figura (3.08) esquematiza o exemplo do poço, onde estão representadas as

tensões nas direções y ( yyσ ) e x ( xxσ ), tensão hidrostática na formação ( ),

tensão desviadora ( , pressão na formação ( ), raio de um ponto qualquer(r),

raio do poço(r

0P

)0S 0p

w) e o ângulo ( . )θ

50

ooxx SP −=−σ

ooyy SP +=−σ

r

θ

wr

op

Figura 3.08 – Esquema de um poço em um meio poroelástico

A perfuração é simulada removendo no instante t = 0 as tensões atuantes na

parede do poço e impondo no poço um valor nulo de poropressão. O

carregamento foi dividido em três modos, uma parcela considerando a tensão

hidrostática in situ, a desviadora in situ e a poropressão atuante na formação. A

soma dos três efeitos com as tensões in situ reproduzem o efeito da perfuração.

3.2.6.2.1 Carregamento modo 1

Neste modo, o poço é submetido a um estado de tensão hidrostático. As

condições de contorno para este modo são:

01 Prr =σ (3.33)

01 =θσ r (3.34)

01 =p (3.35)

A solução para este modo é dada por:

2

2

0

1

rr

Pwrr −=

σ (3.36)

2

2

0

1

rr

Pw=θθσ

(3.37)

51

3.2.6.2.2 Carregamento modo 2

Neste modo é considerada apenas à ação da poropressão na formação. As

condições de contorno para este modo são:

02 =rrσ (3.38)

02 =θσ r (3.39)

02 pp = (3.40)

A solução é dada no campo de Laplace

( )( )βξ

0

0

0

~

KK

pps

−= (3.41)

( )( )

( )( )

−−=

βββ

ββξ

ησ

0

12

2

0

1

0

~

2KK

rr

KK

rr

ps wwrr (3.42)

( )( )

( )( )

( )( )

+−=

βξ

βββ

ββξ

ησθθ

0

0

0

12

2

0

1

0

~

2KK

KK

rr

KK

rr

ps ww (3.43)

O símbolo ~ representa a variável no campo de Laplace, é a função de

Bessel modificada de segundo tipo de ordem zero, é a função de Bessel

modificada de segundo tipo de ordem 1, s é a variável da função núcleo de

transformação, c é o coeficiente de difusividade definido anteriormente pela

expressão (3.20),

0K

1K

csr=ξ e

csrw=β .

3.2.6.2.3 Carregamento modo 3

O poço, neste modo, é submetido a um estado de tensão desviador. A

variável θ introduzida neste modo refere-se ao ângulo medido a partir do eixo x

no sentido anti – horário. As condições de contorno para este modo são:

)2cos(03 θσ Srr −= (3.44)

)2sen(03 θσ θ Sr = (3.45)

52

03 =p (3.46)

A solução para o campo de tensão e poropressão é dada por

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

( ) ( )θξ 2cos11

31911

2

2

221

22

0

~

−+

+−−

+−=

rr

CvvBKC

vvvvvB

SPs w

u

u

uu

u (3.47)

( )( ) ( ) ( ) ( )θξ

ξξ

ξ2cos3

1161

131

4

4

32

2

222110

~

−

−−

+

−+

=rr

Crr

Cv

KKCvvB

SSs ww

uu

urr (3.48)

( )( ) ( ) ( ) ( )θξ

ξξ

ξθθ 2cos3611

131

4

4

322110

~

+

++

−+

−=rr

CKKCvvB

SSs w

u

u (3.49)

( )( ) ( ) ( ) ( ) ( θξ

ξξ

ξθ 2sen3

12131

1312

4

4

32

2

222110

~

−

−−

+

−+

=rr

Crr

Cv

KKCvvB

SSs ww

uu

ur ) (3.50)

onde:

( )( )( )( )12

1 1112

DDvBvvv

Cu

uu

−+−−

=β

(3.51)

( )12

22

14DD

DvC u

−−

= (3.52)

( ) ( ) ( )( )12

2123

8DD

KvvDDC u

−−++

=β

ββ (3.53)

( ) ( )β11 2 KvvD u −= (3.54)

( ) ( )ββ 22 1 KvD −= (3.55)

A técnica utilizada para a inversão das expressões no campo de Laplace será

exposta no próximo tópico.

Para verificar a capacidade do programa em simular este problema, simula-

se o exemplo de um poço vertical de raio 0.1 m submetido a um estado de tensão

de MPaxx 2=σ e MPayy 4=σ , com uma poropressão na formação de 1 . A

figura 3.09 mostra a malha utilizada na simulação, o contorno está situado a uma

distância de 50 vezes o raio do poço. O material utilizado é o arenito de Berea.

Por uma imposição do programa a tensão é colocada como efetiva, tal como no

caso da coluna poroelástica, obedece-se à convenção do programa. A convenção

de sinais adotada para esforços mecânicos neste trabalho será a mesma do

programa ABAQUS, onde o esforço de compressão é negativo.

MPa

53

Figura 3.09 – malha de elementos finitos utilizada na simulação do poço.

0.00

0.20

0.40

0.60

0.80

1.00

1.20

0 10 20 30 40 50Raio adimensional

Poro

- pr

essã

o (M

Pa)

60

analítica

abaqus

t = 0.9

t = 12.9 sregime permanente

Figura 3.10 – Poropressão ao longo da direção 0=θ - 62.0=α

0.00

0.20

0.40

0.60

0.80

1.00

1.20

0 10 20 30 40 50Raio adimensional

Poro

- pr

essã

o (M

Pa)

60

analíticaabaqus

t = 0.8 s

t = 12.8 sregime permanente

Figura 3.11 – Poropressão ao longo da direção 0=θ - 1=α

54

-12.00

-11.00

-10.00

-9.00

-8.00

-7.00

-6.00

-5.00

-4.00

-3.00

-2.000 10 20 30 40 50 6

Raio adimensional

σ θθ

(MPa

)

0

analíticaabaqus

t = 0.9 s

t = 12.9 s

regime

Figura 3.12 – Tensão tangencial ao longo da direção 0=θ - 62.0=α

-12.00

-11.00

-10.00

-9.00

-8.00

-7.00

-6.00

-5.00

-4.00

-3.00

-2.000 10 20 30 40 50

Raio adimensional

σ θθ (

MPa

)

60

analítica

abaqus

t = 0.8 s

t = 12.8 s

regime t

Figura 3.13 - – Tensão tangencial ao longo da direção 0=θ - 1=α

55

-2.50

-2.00

-1.50

-1.00

-0.50

0.000 10 20 30 40 50 6

Raio adimensional

σ rr (

MPa

)

0

analíticaabaqus

t = 0.9 s

t = 12.8 s

regime permanente

Figura 3.14 - – Tensão radial ao longo da direção 0=θ - 62.0=α

-2.50

-2.00

-1.50

-1.00

-0.50

0.000 10 20 30 40 50 60

Raio adimensional

σ rr (

MPa

) analítica

abaqus

t = 0.8 s

t = 12.8 s

regime permanente

Figura 3.15 - Tensão radial ao longo da direção 0=θ - 1=α

Os resultados obtidos pela simulação numérica para a poropressão (figuras

3.10 e 3.11) mostram uma boa concordância com a solução analítica, mesmo para

56

a região próxima ao poço onde o gradiente de poro-pressão é elevado. Porém, à

medida que se afasta do poço e para tempos maiores, o resultado dado pela

simulação numérica se afasta um pouco da analítica, isto pode ser em virtude da

condição de contorno colocada na simulação numérica não representar o infinito.

Uma solução para este problema seria aumentar mais a malha ou aplicar

elementos que representassem o meio infinito, mas esta opção não existe no

programa para o caso do acoplamento fluido mecânico. Simulações com malhas

cujo contorno estava a 10, 20, 30 e 50 vezes o raio do poço foram feitas.

Observou-se pouca variação de resultados para malhas com o contorno acima de

20 vezes o raio do poço.

A mesma observação da poropressão não é feita para as tensões (figuras

3.12, 3.13, 3.14 e 3.15), uma possível razão é que os resultados para a tensão

foram analisados nos pontos de integração do elementos, os quais não se

localizam na direção 0=θ como é feita para a solução analítica. Uma tentativa de

interpolar os resultados de tensão para os nós foi feita, entretanto o programa não

permitiu identificar a que nó um valor de tensão estava associado. Todavia, o

comportamento das curvas de tensão versus deformação são muito semelhantes às

obtidas pela solução analítica o que garante uma certa confiabilidade.

3.2.6.3 Solução para a inversa da transformada de Laplace

A inversa da transformada de Laplace tem sido amplamente estudada na

literatura e vários métodos de solução foram propostos para sua resolução, que

levam em conta a natureza mal condicionada da sua solução. A inversa da

transformada de Laplace é usada nas expressões analíticas dadas no item anterior,

o método selecionado para a inversão foi o de Stehfest (1970). Este método está

baseado na amostragem de dados de acordo com uma série de delta. A solução

aproximada no tempo é dada pela seguinte expressão

( )∑=

=

N

iia i

TPV

TF

1 ')2ln('

'2ln (3.56)

57

( ) ∑+

=

+

−−−−−=

)2/,min(

21

2/2/

)!2()!()!1(!)!2/()!2(1'

Ni

ik

NiN

ikkikkkNkkVi (3.57)

onde:

Fa: valor aproximado da função real

' : é a variável da função real T

N: número de termos da série

P: é a função real no campo de Laplace, onde a variável da função núcleo de

transformação representada convencionalmente por s é igual a iT

)2ln( .

A expressão (3.57) difere daquela exposta no trabalho de Stehfest pelo

termo 2/NK no numerador, que substitui o termo 12/ +NK no trabalho original. É

importante notar também o valor que K assume na expressão (3.57). Como a

expressão (3.57) utiliza o fatorial, o domínio desta função exige que K seja

sempre um número inteiro. Portanto, K será aproximado para um valor inteiro,

quando i for igual a um numero par. Para que todos os termos participem do

somatório, K deverá ser truncado, do contrário erros surgirão durante a inversão

numérica.

O valor de N, número de termos da série, deverá assumir valores pares pelo

mesmo motivo do domínio da função, variando numa faixa entre 8 e 20.

Entretanto o uso de altos valores para N podem conduzir a erros na inversão. Um

pequeno teste variando o valor de N é recomendado, a fim de verificar a variação

dos resultados obtidos e na escolha do melhor valor de N.

4 Modelo constitutivo implementado

O conhecimento das condições que levam uma rocha reservatório a produzir areia

constitui um fator importante no uso e escolha dos métodos de controle. A previsão de

como a rocha responderá sob condições de operação de um poço necessariamente

envolve a utilização de uma relação constitutiva. Materiais geológicos mesmo

submetidos a pequenos carregamentos costumam se deformar de tal maneira que não

recuperam mais sua forma inicial, o que envolve além da deformação elástica a

consideração de deformações plásticas no modelo constitutivo.

Dois comportamentos ligados à deformação se destacam na rocha quando

submetida a carregamento. A dilatância, ou seja, o aumento de volume quando o estado

de tensão se aproxima do nível de ruptura da rocha e o colapso de poros que é o

fechamento de poros devido à ação da parcela hidrostática do carregamento. Exemplos

típicos deste comportamento são o fechamento de poços e a redução da permeabilidade.

Portanto, modelos que consigam simular tais situações são de grande valia na análise da

produção de areia e mais genericamente de processos mecânicos ao redor do poço.

4.1 Critérios da Plasticidade

A implementação e o desenvolvimento de modelos elastoplásticos na

geomecânica computacional apóiam-se fortemente em um número finito de critérios de

escoamento, funções de potencial plástico e leis de evolução (Jeremic, 2002). Alguns

aspectos deste embasamento serão abordados neste tópico e posteriormente utilizados na

implementação do algoritmo de integração do tipo explícito.

59

4.1.1 Critério de escoamento

O critério de escoamento define o limite entre o comportamento elástico e plástico

de um material, especificando se um determinado estado de tensões provoca ou não a

plastificação do material. É expresso normalmente em termos de tensão ou de seus

invariantes, como exemplificado na expressão (4.01):

0),,,,,( 231312332211 == σσσσσσfF (4.01)

Admitindo-se que o material seja homogêneo e isotrópico, F poderá ser

simplificada e expressa em termos das tensões principais:

0),,( 321 == σσσfF (4.02)

Como F fornece um escalar e as tensões principais permitem uma representação

gráfica a partir de um sistema de três eixos, o critério de escoamento pode ser

interpretado como uma superfície. Portanto, estados de tensão localizados dentro desta

superfície caracterizam o regime elástico e os que repousam nela caracterizam o plástico.

4.1.2 Leis de evolução para endurecimento–amolecimento

A superfície de escoamento não retrata apenas se um material está ou não no

regime elástico, mas também caracteriza a sua evolução ao longo da história de

carregamento. O termo evolução é empregado como uma mudança nas propriedades do

material e variações na sua resistência, relacionando-se aos fenômenos de

“endurecimento” e “amolecimento”.

A expressão do critério de escoamento (4.02) passa ser expressa não somente em

função das tensões ou invariantes, mas também de um parâmetro de endurecimento, o

critério de escoamento assume a forma:

( ) 0'')( =−= hffF σ (4.03)

onde σ representa o estado de tensão atuante em um corpo e o parâmetro de

endurecimento.

'h

60

Uma interpretação geométrica desta mudança é visualizada no plano desviador,

em que a superfície de escoamento ao longo de um carregamento pode transladar, variar

a forma e/ou tamanho, ou uma combinação destes. A figura (4.01) ilustra estes casos.

De acordo com o comportamento da superfície de escoamento, caracterizam -se as

situações de endurecimento isotrópico, cinemático ou misto. A figura (4.01-a) ilustra uma

superfície de escoamento estática, sem variação na forma ou tamanho, caracterizando um

material perfeitamente plástico. Na figura (4.01-b) ocorre o endurecimento isotrópico,

onde a variação de tamanho da superfície representa um aumento nas tensões que iniciam

a plastificação. O deslocamento da superfície de escoamento visto na figura (4.01-c)

representa o endurecimento cinemático.

(c) Figura 4.01 – Representação da superfície de escoamento no plano desviador (Desai,1984).

Para o endurecimento isotrópico, o parâmetro de endurecimento pode ser

relacionado à deformação plástica, a expressão (4.03) assume:

( ) 0')( =−= pffF εσ (4.04)

61

Outra forma para representar o parâmetro de endurecimento seria colocá-lo em

termos do trabalho plástico, a expressão (4.03) adquire a seguinte forma:

( ) 0')( =−= pwffF σ (4.05)

4.1.3 Lei de Fluxo

Na teoria clássica da plasticidade, o incremento de deformação total é decomposto

em uma parcela elástica e outra plástica. A parcela elástica está relacionada à lei de

Hooke como:

ijeij dDd σε 1−= (4.06)

sendo:

D : matriz constitutiva elástica do material edε : vetor incremento de deformação elástica

ijdσ :- vetor incremento de tensão

O incremento de deformação plástica é determinado através de uma lei de fluxo,

como na expressão (4.07). Esta lei é expressa por:

ij

pij

Qddσ

λε∂∂

= (4.07)

onde:

λd : escalar positivo

Q : função potencial plástica

A função potencial plástica (Q) é definida assumindo a existência de uma

superfície perpendicular ao incremento de deformação plástica.

A expressão (4.07) é conhecida também como condição de normalidade. Quando

a função potencial plástica é igual à função de escoamento, diz-se que o fluxo é

associado. Para os materiais geotécnicos é comum que a função de escoamento não

reproduza a condição de normalidade. Neste caso, adota-se outra função, o fluxo é

classificado como não associado.

62

A constante λd conhecida como parâmetro plástico, pode ser obtida através da

condição de consistência, a qual estabelece que um carregamento atuando sobre um

corpo sob estado de tensões plástico levará a outro estado de tensões plástico. Esta

condição exige que o critério de escoamento seja sempre satisfeito. A obtenção deste

parâmetro será feita no item que descreve as relações usadas no algoritmo de integração.

4.2 Modelo Lade-Kim

Este modelo tem sido relatado na literatura (Jeremic e outros, 1998) como capaz

de reproduzir com boa acurácia o comportamento de materiais granulares, representando

características como superfície de plastificação dependente dos três invariantes de tensão,

dilatância ligada a uma lei de fluxo não associativa, dependência do módulo de Young

em relação ao estado de tensão e o comportamento de endurecimento - amolecimento

relacionado às tensões atuantes.

Baseado na revisão e avaliação de dados experimentais, Lade e Kim (1988)

desenvolveram um modelo constitutivo para materiais friccionais como o concreto, areia

e rocha. O comportamento elástico é descrito através da lei de Hooke e o comportamento

plástico pelo critério de escoamento, função potencial plástica e um critério de ruptura,

expresso em termos dos invariantes de tensão. O modelo utiliza doze parâmetros,

incluindo a coesão presente no material.

4.2.1 Comportamento elástico

Como mencionado, o regime elástico obedece à lei de Hooke. O módulo de

Young é colocado em função da tensão principal menor (Lade e Kim, 1988), como segue

abaixo:

63

n

papaKurE

⋅⋅= 3σ

(4.08)

onde e são constantes adimensionais determinadas através de ensaio triaxial sob

vários valores de tensão confinante

Kur n

( )3σ , pa é a pressão atmosférica expressa na mesma

unidade de E e 3σ .

4.2.2 Critério de ruptura

Baseado no fato de que a rocha é um material friccional com muitas

características similares ao solo e o concreto, Lade e Kim ampliaram para rochas o

critério de ruptura desenvolvido para solos com uma envoltória curva (Lade e Duncan,

1977). O critério formulado em termos do primeiro e terceiro invariantes de tensão é

dado por:

11

3

31 27 η=

−

m

paI

II (4.09)

onde:

zyxI σσσ ++=1 (4.10)

( )yxxyzxzzxyzyyzxxzzyyxzxyzxyzyxI ττσττσττσττττττσσσ ++−++=3 (4.11)

sendo e o primeiro e o terceiro invariante de tensão respectivamente, m e 1I 3I 1η são os

parâmetros do modelo. O parâmetro m está relacionado a curvatura da superfície de

ruptura no plano triaxial (figura 4.02-a) e o parâmetro 1η com a seção transversal da

superfície de ruptura no plano octaédrico (figura 4.02 -b).

A resistência à tração e a coesão presente em uma rocha são incluídas através da

translação dos eixos de tensões principais ao longo do eixo hidrostático. Essa translação

inclui nas tensões normais uma parcela que representa a coesão do material dada por:

axx pa ⋅+= σσ_

(4.12)

64

ayy pa ⋅+= σσ_

(4.13)

azz pa ⋅+= σσ_

(4.14)

sendo a o parâmetro que representa a coesão.

Esta adição da parcela que representa a coesão deve ser feita antes da utilização das

expressões do critério de ruptura.

1σ

eixo hidrostático

1σ

32σ

(a)

1σ

3σ2σ

(b)

Figura 4.02 – Superfície de ruptura do modelo Lade–Kim no plano triaxial e octaédrico (Lade e

Kim, 1988).

4.2.3 Função potencial plástica

Lade e Kim propuseram uma lei de fluxo não associada para descrever o

comportamento das deformações plásticas. A função Q descrita na expressão (4.07) é

colocada no modelo Lade - Kim como: µ

ψψ

+−=

paI

II

IIg p 1

22

21

3

31

1

(4.15)

onde:

( )xzzyyxxzzxzyyzyxxyI σσσσσσττττττ ++−++=2 (4.16)

1ψ : fator que pondera a forma circular (do termo ) e a forma triangular (do termo

). Lade e Kim relacionaram este parâmetro ao parâmetro m do critério de ruptura por:

2I

3I

65

27.11 00155.0 m=ψ (4.17)

2ψ : fator que controla a interseção da superfície potencial plástica com o eixo

hidrostático.

µ : coeficiente que determina a curvatura dos meridianos (traço com o plano

triaxial).

2I : segundo invariante de tensão

Reescrevendo a expressão (4.07)

ij

ppij

gdd

σλε∂

∂=

(4.18)

A figura (4.03) mostra a superfície potencial plástica no espaço de tensões e os

vetores de incremento de deformação plástica permitindo identificar que quanto mais o

estado de tensão se afasta de um carregamento hidrostático e se dirige a ruptura, mais a

variação volumétrica deixa de ser compressiva e passa a ser dilatante.

eixo hidrostático

a

1p

σ

a

2p

2σ

Figura 4.03 – Superfície potencial plástica do modelo Lade – Kim (Lade e Kim, 1988)

Substituindo a expressão (4.15) em (4.18), tem-se:

66

( ) ( )

( ) ( )

( ) ( )

( )

( )

( )

−−

−−

−−

−−+−

−−+−

−−+−

=

23

31

122

21

23

31

122

21

23

31

122

21

23

312

122

21

1

23

312

122

21

1

23

312

122

21

1

1

22

22

22

II

II

II

II

II

II

II

IIG

II

IIG

II

IIG

paId

dddddd

xyzzxyzxy

zxyyzxyzx

yzxzxxyyz

xyyxyx

zxzzxz

yzzyzy

pxy

pxz

pyz

pz

py

px

τσττψτ

τσττψτ

τσττψτ

τσσψσσ

τσσψσσ

τσσψσσ

λ

εεγεεε

µ

(4.19)

onde:

( ) ( )1

23

123

21

11123II

IIIG µψµµψ ++−+=

(4.20)

Um fato relacionado à ocorrência da deformação plástica é o cumprimento da

condição de irreversibilidade de Prager, esta condição impõe:

0≥= pijijp ddw εσ (4.21)

sendo dwp o incremento de trabalho plástico

A função potencial plástica deve ser convexa em relação à origem do eixo de

tensões. Lade e Kim propõem que está condição é obedecida mediante a escolha de

parâmetros apropriados. Utilizando a expressão (4.18) e considerando que a função

potencial plástica é uma função homogênea de ordem µ , a equação (4.21) leva a:

ppp dgdw λµ= (4.22)

como pdλ >0, a condição de irreversibilidade exige que:

0≥pgµ (4.23)

Lade mostra em seu trabalho que negativo não é uma escolha adequada,

portanto

pg

µ e devem ser positivos. Como pg 273

31 ≥

II

e 32

21 ≥

−II

são sempre

positivos a condição de irreversibilidade conduz a:

67

0>µ (4.24)

( )327 12 +−> ψψ (4.25)

Em relação a materiais coesivos, o mesmo procedimento de translação dos eixos

descrito no critério de ruptura deve ser adotado.

4.2.4 Função de escoamento

Vários critérios de plastificação baseados em experimentos combinados ou não

com suposições do comportamento plástico exibido pelo material (fluxo associado ou não

associado) têm sido propostos. A identificação de pontos de plastificação em curvas de

tensão-deformação é uma destas propostas, porém não é um processo fácil. Uma das

dificuldades encontradas é o fato da plastificação ser um processo contínuo para materiais

geológicos. Não existindo, portanto, um ponto de plastificação distinto na curva tensão-

deformação.

Lade e Kim propuseram como alternativa analisar os contornos dos parâmetros de

endurecimento em um plano octáedrico e compará-los com dados experimentais da

superfície de plastificação. Adotaram como parâmetro de endurecimento o trabalho

plástico e associaram o valor do trabalho plástico a uma superfície de plastificação. Este

procedimento evitou o uso de trajetórias complicadas de tensão e a determinação de

pontos de plastificação na curva tensão-deformação. Além de representar a plastificação

em termos das deformações cisalhantes e volumétricas.

Através de ensaios com carregamento hidrostático estabeleceu-se uma função

monotônica para o trabalho plástico: P

p paIpaCw

= 1 (4.26)

sendo c e p parâmetros adimensionais do modelo.

Analisando o comportamento de curvas de mesmo trabalho plástico no plano

octáedrico (figura 4.04), a seguinte expressão para a superfície de escoamento foi

elaborada:

68

=pf qh

epaI

II

II

− 1

2

21

3

31

1ψ (4.27)

sendo 1ψ e h parâmetros adimensionais do modelo, q depende do nível de tensão e

assume os seguintes valores:

0=q - durante carregamento hidrostático

10 << q - durante o endurecimento

1=q - durante a ruptura

p1d,1 εσ

p2d,2 εσ

p3d,3 εσ−

Figura 4.04 – Dados experimentais e contornos de trabalho plástico constante no plano

octaédrico para Fuji River Sand (Lade e Kim, 1988).

Na figura (4.04) várias superfícies de plastificação com o mesmo valor de são

plotadas utilizando a expressão (4.27) a qual não está relacionada diretamente com o

trabalho plástico, contudo a função de escoamento representa razoavelmente bem o

contorno do trabalho plástico descrito pelos pontos. Outro fato é que os dados

experimentais indicam um comportamento anisotrópico enquanto a expressão (4.27) um

comportamento de expansão isotrópico, se adequando bem a um modelo de

endurecimento isotrópico.

pf

Como a função de escoamento e o parâmetro de endurecimento se relacionam de

forma única, através da manipulação das expressões (4.26) e (4.27) obtêm-se:

ρρ

11

'1''

=

paw

Df p

(4.28)

69

A expressão (4.28) coloca a função de escoamento em termos do trabalho plástico

para o regime de endurecimento, sendo válida para qualquer estado de tensão.

Recorrendo a expressão (4.05) o critério de escoamento pode ser definido como:

0'' =−= ffF p (4.29)

Lade e Kim propuseram uma função de escoamento em termos do trabalho plástico

para o regime de amolecimento. A expressão idealizada é uma exponencial decrescente

representando um regime de amolecimento isotrópico dada por:

−

= paw

B p

eAf*

' (4.30)

sendo A e B* constantes positivas. Estas constantes são obtidas sem a introdução de

outros parâmetros no modelo. Considerando que a tangente da curva de amolecimento é

igual à tangente negativa da curva de endurecimento e que o valor da função de

escoamento em termos do trabalho plástico no endurecimento é igual à função de

escoamento no amolecimento no ponto de ruptura (figura 4.05), tem-se:

( ) ruptura

p

paw

B

ruptura efA

−

=*

' (4.31)

( )rupturap

hardeningruptura

p

p

fpaw

d

dfB

'1''

,

*

= (4.32)

O mesmo tratamento dado a materiais coesivos no critério de ruptura é dado na

função de escoamento. Entretanto é apontado no modelo que inconsistências na condição

de irreversibilidade podem surgir. Estas inconsistências são tratadas considerando a

existência de uma superfície de escoamento inicial que se estende desde a origem do eixo

de tensões real até uma região próxima a superfície de ruptura. Isto implica que somente

deformações elásticas ocorrem na maior parte de trajetórias de tensão de tração.

Pequenos desvios no comportamento linear elástico são observados quando o material se

aproxima da ruptura por tração caracterizando a plastificação.

70

endurecimento

amolecimento

( )ap

pwd

pdf

( )ap

pwd

pdf−

appw

fp

fp (rup)

)rup(ap

pw

Figura 4.05 – Função de escoamento no endurecimento e amolecimento (Lade e Kim, 1988).

4.3 Modelo de Mohr–Coulomb

Um dos critérios mais utilizados na área geotécnica é o critério de Mohr - Coulomb

e é um dos modelos adotado pelo programa de elementos finitos ABAQUS.

4.3.1 Critério de ruptura

O critério de Mohr - Coulomb assume que a ruptura em um material ocorre quando

a tensão cisalhante atuante em um plano ultrapassar a tensão definida pela expressão:

ϕστ tan*nc −= (4.33)

onde τ é a tensão cisalhante e nσ é a tensão normal atuante em um plano, c é a coesão e

ϕ é o ângulo de atrito interno. O sinal negativo na expressão indica a tensão de

compressão, convenção que é adotada pelo programa.

Em termos das tensões principais o critério é colocado na forma:

( ) ( ) ( ) ( )ϕϕσσσσ cossen21

21

3131 cFrup −+−−= (4.34)

Portanto, diferente do critério do Lade, o critério de Mohr - Coulomb assume que a

ruptura é independente da tensão intermediária.

71

4.3.2 Função de escoamento

Neste modelo a função de escoamento é idêntica a função de ruptura. A função de

ruptura dada anteriormente pode ser reescrita em termos de invariantes de tensão. Esses

três invariantes são:

- tensão equivalente:

−=

^

31' σtraçop (4.35)

- tensão equivalente de Mises:

( )dd SSq :23'= (4.36)

- terceiro invariante de tensão desviadora:

31

:29'

⋅= ddd SSSr (4.37)

sendo neste caso o tensor de tensões, é o tensor de tensões desviadoras dado por

,

^σ

^+σ

dS

IpSd '= I é a matriz identidade. Em relação as operações descritas acima, a

seguinte convenção é adotada pelo ABAQUS, um ponto ( )⋅ representa multiplicação

entre vetores ou matriz, dois pontos ( ): representa a multiplicação das componentes

conjugadas correspondentes de duas matrizes aos pares onde o produto somado.

A função de escoamento é dada por:

*)tan('' cpqRF mc −−= ϕ (4.38)

sendo a medida de tensão desviadora definida como: mcR

( ) ϕππϕ

ϕ tan3

cos31

3sen

cos31,

+Θ+

+Θ=ΘmcR (4.39)

e o ângulo polar desviador dada pela expressão: Θ

( )3

''3cos

=Θqr (4.40)

72

4.3.3 Função potencial plástica

A seguinte lei de fluxo é assumida no programa ABAQUS:

σεε

∂∂

='

'

_G

gdd

pl

pl

(4.41)

onde:

σ∂∂

σ=Q:

'c1'g

(4.42)

( ) ( ) ψ−+ψ= tanpqRtan|ecQ 2mw

20 (4.43)

( )

⋅

−+Θ−−+Θ−

−+Θ−=Θ

ψπ ,3

45cos)1(4)12(cos)1(2

)12(cos)1(4,1

21

2211

21

21

221

1

mc

mw

R

eeeee

eeeR (4.44)

sendo:

Q : é a função potencial plástica

e : parâmetro relacionado à excentricidade que define a taxa com que a função G

se aproxima de uma assíntota (a função potencial tende a uma linha reta à medida que a

excentricidade se aproxima de zero).

ψ : ângulo de dilatância medido no plano p - qRmw

mwR : função desviadora elíptica

1e : excentricidade que descreve o contorno da função potencial plástica no plano

octáedrico.

0|c : coesão no início da plastificação

pl_

ε : deformação plástica equivalente

A função desviadora elíptica usada por Menétrey e Willam (1995) possibilita a

função potencial plástica ser contínua e suave no plano octáedrico garantindo que exista

apenas uma direção de fluxo, como ilustrado na figura (4.06).

73

1e1/2 ≤<

Figura 4.06 – Fluxo potencial plástico pela função de Menétrey e Willam no plano desviador

(Willam e Menetrey, 1995).

5 Algoritmo para integração da relação tensão–deformação

A análise de materiais com comportamento elastoplástico por elementos

finitos é feita de forma incremental e iterativa. A cada estágio do processo de

solução, incrementos de força são aplicados e os respectivos incrementos de

deslocamento calculados pela solução de equações de equilíbrio. As deformações

e tensões são computadas nos pontos de integração de cada elemento usando as

relações deformação–deslocamento e por uma lei elastoplástica de tensão–

deformação.

Caso ocorra a plastificação no material e um modelo de endurecimento

isotrópico esteja sendo utilizado, a solução destas leis elastoplásticas é obtida

resolvendo-se um sistema de equações da forma:

( )••

= εσσ hDep , (5.01)

ph•

=• _

ε ou pwh••

= (5.02)

onde dependerá do parâmetro utilizado para representar o endurecimento. •h

Nestas expressões σ representa o vetor de tensão, ε o vetor de deformação,

o parâmetro de endurecimento e a matriz elastoplástica, o sinal de ponto

acima das variáveis representa a derivada em relação ao tempo no qual a carga é

aplicada. A obtenção da matriz elastoplástica é mostrada no tópico 5.1.2.

h epD

Existem dois tipos de algoritmos para o cálculo das tensões e do parâmetro

de endurecimento. Para a compreensão do processo de integração mostra-se de

forma sintetizada uma expressão geral destes algoritmos.

Considerando a deformação total constituída pelas parcelas elástica e

plástica, a tensão na superfície de escoamento em um ponto qualquer de um

material que tenha se plastificado pode ser escrita como:

( )pnnn D 111 +++ −= εεσ (5.03)

onde o índice n se refere ao incremento, para a deformação:

75

11 ++ += nnn dεεε (5.04)

pn

pn

pn d 11 ++ += εεε (5.05)

introduzindo as expressões (5.04) e (5.05) em (5.03),

( )pnn

pnnn ddD 111 +++ −+−= εεεεσ (5.06)

o que permite escrever: pnnnn DdDd 111 +++ −+= εεσσ (5.07)

fazendo igual a , tem-se a parcela denominada de tensão

tentativa elástica, como mostrado na expressão (5.08):

1++ nn Ddεσ *1+nσ

pnnn Dd 1

*11 +++ −= εσσ (5.08)

A segunda parcela do lado direito da expressão (5.08) representa o

comportamento plástico. Esta parcela é que permite conduzir o estado de tensão

localizado fora da superfície de escoamento para a superfície de escoamento.

Durante o fluxo plástico o gradiente da função potencial plástico varia ao

longo da trajetória de deformação incremental de modo não conhecido. Então,

alguma hipótese deverá ser adotada para possibilitar a integração da relação

constitutiva (5.08) (Nogueira, 1998). Uma forma de integração é a regra do ponto

médio, que adota a hipótese da variação linear das tensões e do parâmetro de

endurecimento ao longo do incremento de deformação, portanto:

( ) 1' ''1 ++ +−= nnn σασασ α (5.09)

( ) 1' ''1 ++ +−= nnn hhh ααα (5.10)

sendo 1+nσ e os valores desconhecidos da tensão e do parâmetro de

endurecimento no fim do incremento,

1+nh

'α é uma constante com valor entre 0 e 1.

Assim, quando a deformação plástica na expressão (5.08) é substituída pela

lei de fluxo, as tensões utilizadas serão as dadas pela expressão (5.09).

Quando 'α é igual a 1 tem–se o algoritmo de backward Euler, classificado

como do tipo implícito, pois trabalha somente com as tensões do fim do

incremento. A vantagem deste tipo de algoritmo é que tensão resultante

automaticamente satisfaz o critério de escoamento e não exige, caso ocorra uma

mudança do estado elástico para um plástico, o cálculo da interseção da trajetória

de tensão com a superfície de plastificação.

76

Apesar de ser um método poderoso, é um algoritmo de difícil

implementação para modelos complexos, pois exige o cálculo da derivada de

segunda ordem da função potencial e de escoamento. Além disto, o algoritmo

pode apresentar divergências para superfícies de escoamento com vértices ou

rápidas mudanças na curvatura (Sloan, 2001).

Para 0=α tem – se o algoritmo do tipo explícito, que trabalha somente

com as tensões e parâmetro de endurecimento do início do incremento. O

esquema de forward Euler é um dos mais utilizados neste tipo de algoritmo. A

precisão deste tipo de algoritmo depende do tamanho do incremento, sendo usual

dividir o incremento em subincrementos de igual tamanho para melhorar a

precisão da integração. A vantagem deste método é a fácil implementação, pois

trabalha apenas com derivadas de primeira ordem. A implementação de um destes

algoritmos será discutida no próximo tópico.

5.1 Algoritmo do tipo explícito

A desvantagem do método forward Euler citado anteriormente é que não

existe uma regra para se estabelecer o número de subincrementos e também não se

pode garantir que as tensões resultantes satisfaçam o critério de escoamento no

fim do incremento. Outra desvantagem é a necessidade de se conhecer a

interseção do vetor tensão com a superfície de plastificação, quando o estado de

tensão muda do elástico para o plástico.

Sloan (1987) propôs um algoritmo do tipo explícito, cuja característica

chave é o controle de erro no processo de integração pela seleção automática do

tamanho do subincremento à medida que a integração procede. Este trabalho está

fundamentado neste algoritmo.

O método de Runge–Kutta–Dormand-Price é utilizado neste trabalho para

resolver as expressões (5.01) e (5.02) em alternativa ao método de Euler. A

dedução das expressões do método de Runge–Kutta–Dormand–Price não será

mostrada neste trabalho, mas pode ser encontrada em Dormand and Price (1980).

Todavia, para entender o seu mecanismo, será mostrado de forma resumida o

método de Euler, pois ambos tem a mesma concepção.

77

As expressões (5.01) e (5.02) definem um sistema de equações diferenciais

que regem o comportamento mecânico de um corpo, onde a princípio é conhecido

o estado de tensão e deformação atuantes em um corpo, as propriedades do

material e a matriz constitutiva ao nível do elemento. Portanto este problema pode

ser estudado como um problema de valor inicial.

Euler propôs que a solução deste tipo problema poderia ser obtida

aproximando a solução exata através do cálculo de uma série de pontos dentro de

um determinado intervalo. A figura (5.01) ilustra a solução exata (linha cheia) e a

aproximada (pontos) para uma equação diferencial em função de uma variável. A

solução das expressões 5.01 e 5.02 será feita por esta aproximação.

1w

nw1nw +

1v nv 1nv + v

w

Figura 5.01 – Interpretação gráfica da solução proposta por Euler (Boyce e Diprima,

1998)

Segundo Euler, cada valor da variável dependente (w) na curva (figura 5.01)

pode ser obtido pela soma do seu valor anterior com um incremento relacionado a

ela. Este incremento é o produto do incremento da variável independente (v) com

a derivada da tangente (w’) no ponto de valor inicial. A seguinte expressão foi

proposta:

( )nnnn1n w,v'wvww ⋅∆+=+ (5.11)

Nota – se uma semelhança com a expressão (5.08) onde w seria σ , seria

o incremento

v∆

εd e representaria a matriz constitutiva .

Entretanto, a convergência deste método é muito sensível ao incremento da

variável independente. A convergência é obtida, então, por dois modos. O

( nnn w,v'w ) D

78

primeiro é o uso de incrementos cada vez menores e o segundo pela obtenção de

um método mais eficiente. A obtenção de um método mais eficiente é feita

normalmente modificando-se o valor dado por ( )nnn w,v'wv ⋅∆ .

1nw +

nw +

( nnn w,v'w

(( )vwvv nn ∆+∆+ ,

∆

5 201

220351 σ +∆+

5 201

220351 hh ∆+∆+

54 565

8881 σσ +∆+

54 565

8881 hh +∆+∆

Na literatura um método citado é o de Euler aprimorado (Boyce e Di Prima,

1998) ou também conhecido por Euler Modificado. Este método determina

.através de dois estágios, o primeiro estágio utiliza o método de Euler para

determinar , o segundo estágio usa o valor de calculado anteriormente

no lugar de na expressão (5.11). Assim o valor de não dependeria

somente do valor no início do incremento mas também de algum ponto situado

dentro dele, então, pode–se dizer que a derivada

1nw +

1nw +

nw 1

) é avaliada duas

vezes para cada incremento. A expressão desenvolvida tem a seguinte forma:

( ) wwvwvww nnnnn +∆+=+ ','21

1 (5.12) )

O erro de truncamento local por este método, ou seja, o erro devido somente

ao método é proporcional a . O erro de truncamento global que reflete o erro

do truncamento local acrescido do erro dos parâmetros de entrada é proporcional a

. O método é classificado de acordo com o erro de truncamento global,

portanto o método de Euler Modificado é denominado de segunda ordem.

3v∆

2v∆

O método de Runge–Kutta-Dormand-Price utilizado neste trabalho tem o

mesmo princípio de Euler Modificado, entretanto a parcela ( )nnn w,v'wv ⋅ é

avaliada seis vezes, o que conduz a um erro de truncamento global de 5ª ordem. O

erro local para este método, expresso agora em termos de tensão é obtido a partir

das seguintes expressões:

64311 108145

297190

54031 σσσσσσ ∆∆−∆+∆+= −nn (5.13)

64311 108145

297190

54031 hhhhh nn ∆−∆+∆+= − (5.14)

631

^

1

^

216125

079.21

21619 σσσσσ ∆∆−∆+∆+= −nn (5.15)

631

^

1

^

216125

079.21

21619 hhhhh nn ∆−∆+∆+= − (5.16)

onde:

79

niiepi hD εσσ ∆

=∆ , (5.17)

iiniii

Qhhσ

σεσλ∂∂

∆∆=∆ ,, (5.18)

εε ∆∆=∆ nn T (5.19)

para i variando de 1 a 6, têm-se então:

11 −= nσσ

11 −= nhh (5.20)

11251 σσσ ∆+= −n

11251 hhh n ∆+= −

(5.21)

2113409

403 σσσσ ∆+∆+= −n

21

^

13^

409

403 hhhh n ∆+∆+= −

(5.22)

3211456

109

103 σσσσσ ∆−∆+∆+= −n

3211456

109

103 hhhhh n ∆−∆+∆+= −

(5.23)

43211572955

729880

2725

729226 σσσσσσ ∆+∆−∆−∆+= −n

43211572955

729880

2725

729226 hhhhhh n ∆+∆−∆−∆+= −

(5.24)

543211655

1892791

297226

25

270181 σσσσσσσ ∆+∆+∆−∆−∆−= −n

543211655

1892791

297226

25

270181 hhhhhhh n ∆+∆+∆−∆−∆−= −

(5.25)

A expressão para o erro local é obtida subtraindo as expressões (5.15) e

(5.13):

65431

^

28011

4027

7255

6310

36011 σσσσσ ∆+∆−∆+∆−∆=nE (5.26)

80

A mesma expressão é válida para o parâmetro de endurecimento

O processo de integração inicia assumindo que o subincremento de

deformação é igual ao próprio incremento de deformação. A tensão é calculada

somando a tensão do incremento anterior com o incremento de tensão dada pela

expressão (5.17). Em seguida é determinado o erro local e o erro relativo ,

como mostrado na expressão (5.27):

( )'R

n

nER

^σ

^

'= (5.27)

Se o erro relativo for menor do que a tolerância especificada, a tensão é

atualizada e termina o processo de integração. Caso o erro seja maior do que a

tolerância é feita a subdivisão do incremento.

O incremento de deformação passa a ser um somatório de subincrementos

de deformação. Para facilitar a operação de integração é associado ao

subincremento de deformação um subincremento de tempo adimensional (T), tal

que o incremento de deformação seja representando pelo valor unitário de T. O

subincremento de deformação é dado por:

kk T∆∆=∆ εε (5.28)

onde kε∆ é o subincremento de deformação e kT∆ é o subincremento de tempo.

O valor de ∆ é calculado baseado em uma extrapolação: kT

1

_

−∆=∆ kk TqT (5.29)

_q é um valor real positivo, o valor deste número dependerá se o subincremento de

deformação foi integrado com sucesso ou não. Sucessivos subincrementos são

calculados até que a soma deles se iguale ao incremento de deformação ou T seja

igual a 1.

Cabe ainda, definir a interseção do vetor tensão com a superfície de

escoamento, quando o estado de tensão muda do estado elástico para o plástico.

Supondo um estado de tensão elástico ( aσ ), tal que ( ) 0, <hF aσ , então, para um

incremento de tensão σ∆ que resulte na função de escoamento:

( ) 0, >∆+ hF a σσ (5.30)

81

deve-se encontrar a parcela de σ∆ responsável pelo carregamento elástico.

Considerando uma constante ξ tal que:

( ) 0, =hF σ (5.31)

onde:

σξσσ ∆+= a (5.32)

O valor de ξ varia entre 0 e 1. Para o valor 0 não há deformação elástica e

para 1 a deformação é totalmente elástica.

Visto que as equações (5.31) e (5.32) definem uma equação não-linear da

forma ( ) 0=ξF , métodos iterativos podem ser utilizados na obtenção de ξ ,como

a técnica de Newton Raphson, regula-falsi, secante e o Pegasus. Este último foi o

selecionado para determinar a interseção.

A seguir mostra-se o esquema do algoritmo de integração da relação tensão-

deformação.

1 – entrar com a tensão e o parâmetro de endurecimento inicial, o

incremento de deformação e a tolerância para o erro relativo durante a atualização

das tensões;

2 – calcular a tensão tentativa elástica ( )*σ a partir da soma da tensão

inicial ( )0σ e o incremento de tensão elástico ( )eσ∆

εσ ∆=∆ ee D (5.33)

et σσσ ∆+= 0 (5.34)

3 – verificar o critério de escoamento, caso ( ) FtolF −<*σ ou

( )( ) FtolFabs t <σ então o estado de tensão é elástico e as tensões podem ser

atualizadas, terminando o algoritmo;

4 – Caso a condição não seja verificada, determina-se se há a passagem do

estado elástico para o plástico pela condição ( ) 00 <σF e 0)( >tF σ , se a

condição é verificada determina-se a constante ξ e a parcela elástica do

incremento de tensão. Do contrário 0=ξ ;

– Atualiza-se a tensão inicial com a parcela de incremento elástico de

tensão

82

eσξσσ ∆+= 00 (5.35)

5 - Determina-se à quantidade de incremento de tensão que será corrigido

através de:

( ) ee σξσ ∆−=∆ 1 (5.36)

6 – Adotar T = 0 e 1=∆T (tamanho do subincremento);

7 – Determinar o incremento de tensão e parâmetro de endurecimento pelas

expressões (5.01) e (5.02) para i variando de 1 a 6.

iiei DbT λσσ ∆−∆∆=∆ (5.37)

it

iii bh σλ∆=∆ (5.38)

onde:

∂∂

∂∂

−∂∂

∂∂

∆∆∂∂

=∆ 0,max

i

ti

ii

ei

i QhFQDF

TF

σσ

σσ

σσ

λ (5.39)

iσ é a tensão dada pelas expressões de (5.20) a (5.25);

8 – Estimar o erro local de acordo com a expressão (5.26). Atualizar a

tensão e o parâmetro de endurecimento de acordo com as expressões (5.15) e

(5.16). Calcular o erro relativo de acordo com a expressão (5.27) e compará-lo

com a tolerância especificada;

9 – Se o erro relativo é maior do que a tolerância; o incremento de tempo

T∆ deverá ser reduzido. Esta redução é obtida multiplicando-se o incremento

pelo fator : _q

{ }1.0,'/9.0max_

Rtolq = (5.40)

e retornar com este incremento para o passo 6.

Caso o erro relativo seja menor que a tolerância, o incremento de tempo é

aumentado pelo mesmo fator q, mas com o seguinte valor

{ }1.1,'/9.0min_

Rtolq = (5.41)

O período de tempo T, as tensões bem como o parâmetro de endurecimento

é atualizado por:

83

TTT ∆+= (5.42)

^σσ =∆+ Tt (5.43)

^hh Tt =∆+ (5.44)

10 – verificar se o período de tempo T é maior do que 1, se for a integração

para o incremento termina, senão ela retorna para o passo 6 e prossegue até T se

igualar a 1.

5.2 Relações usadas no processo de integração

Para implementar o modelo de Lade-Kim no algoritmo explícito é

necessário deduzir a derivada da função de escoamento e obter o parâmetro

plástico a fim de se conhecer a matriz elastoplástica do elemento. A derivada da

função potencial plástica foi dada na expressão (4.19). A seguir são colocadas as

expressões da derivada da função de escoamento e do parâmetro plástico.

- Derivada da função de escoamento em termos de tensão: a derivada da

função de escoamento, expressão (4.27), é obtida através da regra da cadeia:

ij

p

ij

p

ij

p

ij

p IIfI

IfI

Iff

σσσσ ∂∂

∂

∂+

∂∂

∂

∂+

∂∂

∂

∂=

∂

∂ 3

3

2

2

1

1 (5.45)

A derivada da função de escoamento em relação aos invariantes é dada pelas

expressões:

qh

pp e

paI

IIf

Iq

Ih

If

+

∂∂

++

=∂

∂ 1

2

1

111

3 (5.46)

qh

p epaI

II

If

=

∂

∂ 122

21

2 (5.47)

qh

pp e

paI

II

Iqf

If

−

∂∂

=∂

∂ 123

31

133

ψ (5.48)

_q é uma variável em função do critério de ruptura, portanto:

( )( )

+

−−=

∂∂

m

paI

II

ImS

SIq 1

3

21

1

12*

1

*

1

311

ηαη

α (5.49)

84

( )( )m

paI

II

SIq

−−

=∂∂ 1

23

31

2*13 11 αη

α (5.50)

- Derivada da função de escoamento em relação ao trabalho plástico: em

relação às expressões (4.26) e (4.28):

- Endurecimento

( )11

11' −=

∂

∂ρ

ρρp

p

p wDpaw

f (5.51)

- Amolecimento

paBw

p

pp

epaAB

wf −

−=∂

∂ ' (5.52)

- Derivada dos invariantes em relação às tensões: as derivadas serão

colocadas em chaves obedecendo a seguinte seqüência de tensões:

11σ , 22σ , 33σ , 12σ , 13σ e 23σ

=∂∂

000111

1

σI (5.53)

( )( )( )

+−+−+−

=∂∂

12

31

23

2211

1133

3322

2

222

σσσ

σσσσσσ

σI (5.54)

( )( )( )

−−−−−−

=∂∂

12332331

23111223

23113112

2122211

2311133

2233322

3

222

σσσσσσσσσσσσ

σσσσσσσσσ

σI

(5.55)

O parâmetro plástico é obtido a partir da consideração do critério de

escoamento, ou seja, para qualquer estado de tensão localizado sobre a superfície

85

de escoamento um carregamento ocasionará um incremento nulo na função de

escoamento, tem-se:

0=dF (5.56)

colocando a expressão (5.56) em termos do trabalho plástico e das tensões:

0=∂∂

+∂∂

= pp

dwwFdFdF σ

σ (5.57)

o incremento de trabalho plástico é dado por:

ptp ddw εσ= (5.58)

O incremento de tensão pode ser obtido da relação do incremento de

deformação total com as parcelas elástica e plástica. O incremento de deformação

elástica está relacionado à lei de Hooke:

σε dDd e 1−= (5.59)

como:

pe ddd εεε += (5.60)

isto permite colocar o incremento de tensão em função do incremento de

deformação total e plástico, multiplicando a expressão (5.59) por D

( )pddDd εεσ −= (5.61)

substituindo (5.61) e (5.58) em (5.57)

( ) 0=∂∂

+−∂∂ pt

p

p dwFddDF εσεε

σ (5.62)

usando a lei de fluxo e simplificando a expressão (5.62), obtém-se o parâmetro

plástico

σσσ

εσλ

∂∂

∂∂

−∂∂

∂∂

∂∂

=QD

wFQDF

DdF

d

p

(5.63)

A matriz elastoplástica é deduzida a partir da expressão (5.61). Usando

novamente a lei de fluxo, tem-se:

∂∂

−=σ

λεσ QddDd (5.64)

86

Introduzindo a expressão (5.63) em (5.64) e multiplicando D por ambas as

parcelas dentro do parêntese, tem-se:

ε

σσσ

σσεσ dQD

wFQDF

DFQDDdd

p ∂∂

∂∂

−∂∂

∂∂

∂∂

∂∂

−= (5.65)

fatorando a expressão (5.65):

ε

σσσ

σσσ dQD

wFQDF

DFQDDd

p

∂∂

∂∂

−∂∂

∂∂

∂∂

∂∂

−= (5.66)

Comparando-se a expressão (5.66) com a (5.01), a matriz elastoplástica

é dada por:

epD

σσσ

σσ

∂∂

∂∂

−∂∂

∂∂

∂∂

∂∂

−=QD

wFQDF

DFQDDDep

p

(5.67)

Esta matriz é valida para qualquer modelo constitutivo que obedeça a uma

lei de endurecimento isotrópico.

A matriz constitutiva elástica é dada por:

D =( )

( )( )υυν

2111

−+−E

−−

−−

−−

−−

−−

−−

)1(221

)1(221

)1(221

11

11

11

000000000000000000100010001

νν

νν

νν

νν

νν

νν

νν

νν

νν

(5.68)

e a matriz inversa por:

87

ED 11 =−

( )( )

( )

++

+−−

−−−−

υυ

υυυ

υυυυ

120000001200000012000000100010001

(5.69)

5.3 Detalhes da implementação numérica

A implementação do algoritmo do tipo explícito com modelo Lade–Kim no

programa de elementos finitos ABAQUS foi feita utilizando a opção que o

programa oferece de incluir leis constitutivas através de uma subrotina Fortran.

Esta subrotina é denominada por UMAT (subrotina do usuário para definir o

comportamento mecânico de um material).

O programa ABAQUS permite ao usuário apenas escrever o código. O

compilamento e a linkagem do código são executados pelo próprio ABAQUS. A

subrotina é acionada com o fornecimento das variáveis utilizadas no modelo

constitutivo pelo ABAQUS,. As variáveis são processadas pela subrotina em cada

ponto de integração do elemento. Uma vez atualizadas, as variáveis são

retransmitidas ao ABAQUS. Esta operação é feita a cada incremento.

As variáveis ou parâmetros utilizados foram à tensão, a deformação, a

matriz constitutiva do elemento, o vetor PROPS (cujos valores são os parâmetros

do material e definidos no arquivo inp) e o vetor statev. O vetor statev permite que

qualquer variável criada no código seja armazenada. Fornecidas estas variáveis, o

algoritmo explícito pode corrigir as tensões e calcular a matriz constitutiva do

elemento.

Apesar da inclusão de subrotinas oferecer certa flexibilidade na análise de

materiais, deve-se observar que o programa não disponibiliza a “debugagem” da

subrotina implementada, ou seja, o usuário por precaução deve testá-la fora do

ABAQUS antes de utilizar o programa.

88

5.4 Comparação entre ensaios de laboratório e simulação numérica

A implementação foi verificada com três materiais. Um deles é descrito no

trabalho de Lade–Kim (1988). Os outros são arenitos cujas características

geomecânicas e parâmetros associados ao modelo de Lade-Kim foram

determinados por Barroso (2001). Nestes exemplos a linha cheia representa a

solução dada pelo ABAQUS e os pontos os resultados do ensaio. Para o arenito de

Vila Velha é colocada além da análise com o ABAQUS, a resposta obtida também

simplesmente usando a implementação que é representada por uma linha.

O material selecionado do trabalho de Lade–Kim (1988) é a fine silica sand,

pois ele permite verificar tanto a curva tensão versus deformação como a

deformação volumétrica versus deformação axial. Os parâmetros do modelo

Lade-Kim para esses materiais são:

Kur n m 1η c P 2ψ µ h α a υ

1170 0,53 0,1 24,7 0,324e-3 1,25 -3,69 2,26 0,355 0,515 0 0,2

Tabela 5.01 – Parâmetros do modelo Lade–Kim para a fine sílica sand

Kur n m 1η c P 2ψ µ h α a υ

3823,6 0,560 0,6648 17465 1,54e-3 1,00 -2,795 2,370 2,696 0,001 12 0,18

Tabela 5.02 – Parâmetros do modelo Lade–Kim para o arenito Vila Velha

Kur n m 1η c P 2ψ µ h α a υ

4003,8 0,5848 0,5356 2725 0,219e-4 1,2272 -0,1 0,9343 4,6696 0,1 5,17 0,2

Tabela 5.03 – Parâmetros do modelo Lade–Kim para o arenito Rio Bonito

Os três exemplos ilustrados reproduzem um ensaio triaxial convencional.

Na fase de amolecimento, alguns estados de tensão apresentaram dificuldade para

serem simulados. Esta dificuldade pode ser associada ao aparecimento de

autovalores negativos na matriz de rigidez, o que indicaria movimento de corpo

rígido. O programa apresenta como alternativas para minimizar esta instabilidade

a utilização do método de comprimento de arco, a introdução de forças viscosas

para reduzir os deslocamentos e o controle por deslocamento. O controle de

deslocamento consiste reduzir o incremento de força aplicado. Destes mecanismos

89

o que obteve melhor êxito foi o controle de deslocamento. Os resultados da

simulação são mostrados nas figuras 5.02 a 5.07.

Adotou-se como convenção de sinais a do programa ABAQUS, ou seja,

esforços de compressão negativos e de tração positivos.

-0.00

-200.00

-400.00

-600.00

-800.00

-1000.00

-1200.00

-1400.00

-0.00 -5.00 -10.00 -15.00 -20.00 -25.00 -30.00 -35.00

ε11(%)

σ 11

- σ33

(Mpa

)

lade - kimabaqus

101,325 KPa

202,65 KPa

506,625 KPa

Figura 5.02 – Curva tensão versus deformação axial para a fine sílica sand.

-1.00

-0.50

0.00

0.50

1.00

1.50

2.00

-0.00 -5.00 -10.00 -15.00 -20.00 -25.00 -30.00 -35.00

ε11(%)

ε vol

(%)

Lade - Kim

abaqus

101,325 KPa

202,65 KPa

506,625 KPa

Figura 5.03 - Curva deformação volumétrica versus deformação axial para a fine sílica

sand.

90

-0.00

-10.00

-20.00

-30.00

-40.00

-50.00

-60.00

-70.00

-0.00 -0.20 -0.40 -0.60 -0.80 -1.00 -1.20

ε11(%)

σ 11 -

σ33 (M

Pa)

ensaioabaqus

uniaxial

2 MPa

5 MPa

10 MPa

Figura 5.04 - Curva tensão versus deformação axial para o arenito Rio Bonito.

-0.50

0.00

0.50

1.00

1.50

2.00

2.50

-0.00 -0.20 -0.40 -0.60 -0.80 -1.00 -1.20ε11(%)

ε vol

(%)

ensaioabaqus

uniaxial

2 MPa

5 MPa

10 Mpa

Figura 5.05 - Curva deformação volumétrica versus deformação axial para o arenito Rio

Bonito.

91

-0.00

-20.00

-40.00

-60.00

-80.00

-100.00

-120.00

-0.00 -0.20 -0.40 -0.60 -0.80 -1.00 -1.20 -1.40 -1.60 -1.80 -2.00

ε11(%)

σ 11 -

σ33

(MPa

)

ensaioalgoritmoabaqus

uniaxial

2 MPa

5 MPa

10 MPa

Figura 5.06 - Curva tensão versus deformação axial para o arenito Vila Velha.

-1.00

-0.50

0.00

0.50

1.00

1.50

2.00

2.50

-2.50-2.00-1.50-1.00-0.500.00

ε11(%)

ε vol

(%)

ensaioalgoritmoabaqus

uniaxial

2 MPa

5 MPa

10 MPa

Figura 5.07 – Curva deformação volumétrica versus deformação axial para o arenito Vila

Velha.

92

A implementação do modelo Lade–Kim dentro do programa obteve bons

resultados quando comparados ao fine sílica sand (figuras 5.02 e 5.03) estudado

por Lade–Kim, reproduzindo com boa acurácia o comportamento das curvas

tensão versus deformação e deformação volumétrica versus deformação axial.

Uma boa concordância foi obtida para o arenito de Vila Velha (figuras 5.06 e

5.07) estudado por Barroso (2001). Entretanto para o arenito Rio Bonito a mesma

concordância não se verificou (figuras 5.04 e 5.05), uma possível razão para isso

estaria na representatividade do comportamento do material pelos parâmetros

obtidos, já que para os outros dois materiais a concordância foi satisfatória. Outro

fato observado é o comportamento dilatante da deformação próximo à ruptura

para boa parte das tensões aplicadas as amostra.

6 Análise da estabilidade de poços hipotéticos

Após comparar os resultados da simulação numérica com os ensaios

triaxiais, procede-se à avaliação da condição de estabilidade de um poço vertical

hipotético com a finalidade de se conhecer a zona de ruptura na sua vizinhança e

os processos de produção de areia.

Os materiais utilizados nesta análise foram o arenito Rio Bonito e o Calcário

do Campo de Congro descritos no capítulo anterior. O caso estudado é um poço

vertical hipotético submetido a um estado plano de deformação. O fluido de

perfuração é considerado penetrante, ou seja, o fluido não contribui para a

resistência mecânica da parede do poço. Os valores de carregamento no contorno

e o drawdown aplicados são modificados através da análise, a fim de verificar

qual etapa na vida de um poço tem maior influência na sua estabilidade.

A rocha para ambos materiais é considerada um meio poroso, contínuo e

isotrópico. Adicionalmente considerou-se no calcário a dependência da

permeabilidade com o estado de tensão. A figura 6.01 mostra a representação

esquemática do carregamento aplicado. O ângulo θ é medido no sentido anti–

horário a partir da direção de hσ .

r

mrw 1.0=

Hσ

hσθ

Figura 6.01 – Representação esquemát

ica poço vertical a ser estudado

94

Os casos a serem analisados adotam a poro–pressão na formação de 20 MPa

e a pressão no poço terá os valores de 10, 15 e 18 MPa. Os valores de

permeabilidade foram obtidos dos trabalhos de Barroso (2001) e Soares (2001),

para o arenito Rio Bonito a permeabilidade tem o valor de 1534 mD e para o

calcário do Campo de Congro 1,5 mD.

Duas malhas foram utilizadas, a primeira malha com 1057 elementos e 1153

nós simula o comportamento do poço com o arenito Rio Bonito e a segunda malha

mais refinada para o calcário, possuindo 2379 nós e 2280 elementos (figura 6.02).

malha 1

malha 2

Figura 6.02 - Malhas utilizadas na simulação, a esquerda a malha 1 e a direita a malha 2.

95

A malha 1 utilizada no arenito Rio Bonito, apesar de ter um refinamento

menor que a malha 2, mostrou obter resultados praticamente semelhantes a da

malha 2. Para o caso do calcário, o interesse em verificar a dependência da

permeabilidade com o estado de tensão exige um melhor conhecimento das

deformações em todo o reservatório, tornando necessária o uso de uma malha

mais refinada.

A simulação do poço é realizada em duas etapas, a primeira simula o

processo de escavação e a segunda a produção do poço. A escavação é feita

aplicando-se um carregamento na fronteira externa (contorno), enquanto, na

parede do poço não é aplicado qualquer carregamento ou restrição de

deslocamento. Na etapa de produção, o valor da diferença de pressão entre o poço

e a poro-pressão na formação, drawdown, é aplicado na fronteira interna (poço).

Por uma imposição do programa de elementos finitos o carregamento é aplicado

na forma de tensão efetiva e a pressão do fluido na forma de excesso de pressão.

6.1. Poço no arenito Rio Bonito

O objetivo deste caso é determinar sob que condições o poço vertical será

levado à ruptura. Para isto, plota–se os gráficos de isofaixa de tensão, deformação

volumétrica e a região em que ocorre a ruptura do material.

A plotagem da região do poço que sofreu ruptura foi possível através da

utilização da variável S (nível de tensão) do modelo Lade-Kim. S é definida como

a razão entre a função de ruptura (definida no tópico 4.2.2) e o parâmetro 1η . A

variável S pode assumir valores entre 0 e 1, onde valores abaixo de 1 indicam que

o material não atingiu a ruptura.

O carregamento no contorno e a pressão aplicada no poço para o arenito de

Rio Bonito estão descritos na tabela 6.01, o sinal de apostrofo indica tensão

efetiva.

96

Carregamento no contorno (MPa) Pressão (MPa) Carregamento

Hσ 'Hσ hσ '

hσ wp op

1 40 20 30 10 18 20

2 45 25 30 10 18 20

3 48 28 30 10 18 20

4 48 28 30 10 10 20

Tabela 6.01 – Carregamento no contorno e pressão aplicada no poço e na formação

para o arenito Rio Bonito

Na figura 6.03 mostra-se o campo de tensões para os carregamentos 1, 2 e 3.

Nota–se que não há diferença significativa na tensão principal maior entre os

processos de escavação e produção do poço. O acréscimo do carregamento no

contorno induz a uma concentração de tensão na vizinhança do poço, mas não se

nota o avanço desta concentração para o interior da formação com o aumento do

carregamento.

A figura 6.04 descreve a tensão principal menor para os carregamentos 1, 2

e 3. A influência do drawdown (etapa de produção) é pequena como no caso da

figura 6.03. Diferente do que é descrito na solução elástica de Kirsch, nota–se

uma concentração de tensão radial na vizinhança do poço. Esta concentração,

provavelmente, deve-se ao módulo de Young ser dependente do estado de tensão.

Esta dependência aumentaria a rigidez do material próximo à região de ruptura, o

que é desconsiderado na teoria elástica.

Para a figura 6.05, os carregamentos 3 e 4 mostram um aumento da

concentração da tensão principal maior na vizinhança do poço na direção de

0=θ com a aplicação de drawdown maiores. Entretanto, não houve avanço da

concentração de tensão para o interior da formação. Nota–se um descarregamento

de tensão na parede do poço maior para o caso 4, provavelmente, uma

conseqüência do amolecimento do material.

97

Os carregamentos 3 e 4 ilustrados na figura 6.06 mostram a tensão principal

menor tornando-se mais compressiva à medida que o drawdown aumenta. Isto

contraria a idéia de que o fluxo levaria a tensão radial para uma situação de baixa

tensão de compressão ou de tração na vizinhança do poço.

1σ

Carregamento 1: MPaddMPaMPa hH 2,10',20' === σσ

1σ

Carregamento 2: MPaddMPaMPa hH 2,10',25' === σσ

1σ

Carregamento 3: MPaddMPaMPa hH 2,10',28' === σσ

Figura 6.03 – Campos de tensão principal maior na vizinhança do poço para os casos 1,

2 e 3 com drawdown de 2 MPa. A esquerda representa – se a escavação e a direita a

utilização do poço.

98

3σ

Carregamento 1: MPaddMPaMPa hH 2,10',20' === σσ

3σ

Carregamento 2: MPaddMPaMPa hH 2,10',25' === σσ

3σ

Carregamento 3: MPaddMPaMPa hH 2,10',28' === σσ

Figura 6.04 – Campos de tensão principal menor na vizinhança do poço para os casos 1,

2 e 3 e drawdown de 2 MPa. A esquerda representa –se a escavação e a direita a

utilização do poço.

99

1σ

Carregamento 3 e 4 (escavação) - MPaMPa hH 10',28' == σσ

1 σ

Carregamento 3 (produção): MPaddMPaMPa hH 2,10',28' === σσ

1σ

Carregamento 4 (produção): MPaddMPaMPa hH 10,10',28' === σσ

Figura 6.05 – Distribuição da tensão principal maior na vizinhança do poço para os casos

3 e 4, ilustrando os gráficos sobre a escavação como indicado e os demais com o poço

em produção.

100

3σ

Carregamento 3 e 4 (escavação): - MPaMPa hH 10',28' == σσ

3σ

Carregamento 3 (produção): MPaddMPa hH 2,10',28' === σσ

3σ

Carregamento 4 (produção): MPaddMPa hH 10,10',28' === σσ

Figura 6.06 – Campos de tensão principal menor na vizinhança do poço para os casos 3

e 4. ilustrando os gráficos sobre a escavação como indicado e os demais com o poço em

produção.

101

S

Carregamento 1: MPaddMPa hH 2,10',20' === σσ

S

Carregamento 2: MPaddMPa hH 2,10',25' === σσ

S

Carregamento 3: MPaddMPa hH 2,10',28' === σσ

Figura 6.07 - Zona rompida na vizinhança do poço para os casos 1, 2 e 3 e drawdown de

2 MPa. À esquerda está representada a escavação do poço e à direita o poço em

produção.

102

S

Carregamento 3 e 4 (escavação) - MPaMPa hH 10',28' == σσ

S

Carregamento 3 (produção): MPaddMPa hH 2,10',28' === σσ

S

Carregamento 4 (produção): MPaddMPaMPa hH 10,10',28' === σσ

Figura 6.08 - Zona rompida na vizinhança do poço para os casos 3 e 4, ilustrando os

gráficos sobre a escavação como indicado e os demais o poço em produção.

103

vε

Carregamento 1 (escavação): - MPaMPa hH 10',20' == σσ

vε

Carregamento 2 (produção): MPaddMPaMPa hH 2,10',20' === σσ

Figura 6.09 – Campos de deformação volumétrica na vizinhança do poço para o caso 1.

O primeiro gráfico representa a escavação como indicado e o seguintes o poço em

produção.

A zona de ruptura descrita na figura 6.07 não muda significativamente sua

forma e tamanho da etapa de escavação para a de produção do poço. Para os

carregamentos 1, 2 e 3, mesmo a razão entre h

Hσ

σ variando de 2 a 2,8, não

promove alterações perceptíveis na zona rompida. Esta dificuldade da propagação

da zona rompida além da vizinhança do poço pode ser associada à resistência

mecânica descrita no item 5.4.

Os ensaios descritos na figura 5.04 mostram que uma pequena mudança da

tensão de confinamento no arenito aumenta em muito a resistência mecânica da

rocha. Na situação descrita pela figura 6.07, à medida que se afasta da parede do

poço, a tensão confinante ( 3σ ) aumenta, o que por sua vez contribui para o

aumento da resistência, justificando assim a pequena propagação da zona de

ruptura.

104

O pequeno avanço da concentração de tensão descrito nas figuras 6.03 e

6.05 pode ser explicado pelo mesmo motivo, uma vez que a concentração de

tensão principal maior representa é um indicador da região de ruptura.

-90

-80

-70

-60

-50

-40

-30

-20

-10

00 1 2 3 4 5 6 7 8 9 10

Raio adimensional

σ 1 (M

pa)

EscavaçãoProdução - dd= 2 MPaProdução - dd = 10 MPa

-25.00

-20.00

-15.00

-10.00

-5.00

0.00

5.00

0 1 2 3 4 5 6 7 8 9

Raio adimensional

σ 3

10

EscavaçãoProdução - dd = 2MPaProdução - dd = 10 MPa

Figura 6.10 Curva tensão principal maior-raio adimensional e tensão principal menor-raio

adimensional.

O gráfico da tensão principal maior versus raio adimensional ilustrado na

figura 6.10, permite observar que o acréscimo do valor da tensão principal maior

da etapa de escavação para a de produção é praticamente idêntico ao drawdown

105

aplicado. Constata-se através dos picos de tensão indicados nesta curva, que não

houve avanço da concentração de tensão com o aumento do drawdown, como foi

descrito na figura 6.05. O amolecimento do material mencionado anteriormente

também é verificado pelo decréscimo da tensão principal maior desde o pico de

tensão até à parede do poço.

Pelo gráfico de tensão principal menor versus raio adimensional da figura

6.10, nota-se o aumento da tensão principal menor com o aumento de drawdown,

o que contribui para aumento da resistência mecânica do arenito. O que justifica a

pequena propagação da zona de ruptura devido ao aumento do confinamento.

Como observado anteriormente e constatado neste gráfico, o desenvolvimento de

tensões de tração é mínimo.

As mesmas observações da figura 6.07 se aplicam à figura 6.08. Para ambas

as figuras, observou–se que a zona de ruptura não é alterada em sua forma ou

tamanho com a aplicação de um drawdown maior, o que indica a escavação como

a etapa que mais influencia na estabilidade do poço. A extensão da zona de

ruptura pode ser observada também pela distancia da parede do poço até o pico de

tensão principal maior.

A figura 6.09 mostra a deformação volumétrica para o caso 1 próximo à

parede do poço assumindo um caráter dilatante. A extensão da região que sofreu

dilatação é pequena e parece assumir a mesma forma que a região rompida.

6.1.1 Estimativa da produção de sólidos

Até o momento foi apresentado apenas o modo como o arenito responde

diante das condições a que foi submetido. Entretanto ainda não se mostrou uma

forma de como identificar o processo de produção de sólidos. Na literatura

referente à produção de sólidos (Dusseault, 1989), um mecanismo identificado

como o iniciador do processo é a redução da coesão devido à ruptura por

compressão do material. Embora haja um consenso quanto a este mecanismo,

ainda não existe um estudo quantitativo que indique o quanto à rocha perde da sua

coesão.

106

Apoiando neste fato, a fim de verificar a produção de sólidos, arbitrou – se a

situação de que a coesão do arenito na região de ruptura é de 20% da coesão

inicial, o que pelos parâmetros do modelo é da ordem de 0,0517 MPa. Esta

verificação consiste na simples comparação do valor da tensão radial com o valor

arbitrado para a coesão, caso a tensão radial ultrapasse este valor ocorrerá à

produção de sólidos. Para representar tal situação foi escolhido o carregamento 3

( MPaMPa hH 10',28' == σσ e MPadd 2= ), pois como mencionado

anteriormente as tensões de tração tem sido verificadas apenas para baixos valores

de drawdown.

A figura 6.11 ilustra a situação da produção de sólidos, os pontos negros

representam o local onde a produção de sólidos é iniciada.

Figura 6.11 - Indicador dos pontos de início da produção de sólidos.

A reduzida região que produz sólido é decorrente da pequena zona rompida

(figuras 6.07 e 6.08). A descontinuidade da região que produz sólidos pode ser

associada ao fato de que a região submetida a maior tensão tangencial (próxima de

º0=θ ) teve a coesão reduzida de forma mais rápida, o que permitiria o aumento

das tensões de tração além do limite estabelecido. Para a região próxima de

º90=θ , o descarregamento a que é submetida facilitaria o surgimento de tensões

de tração de valores mais elevados, ultrapassando desta forma a coesão presente

no material.

107

6.2. Poço hipotético no Calcário do Campo de Congro

O comportamento de um poço vertical escavado em um calcário será

avaliado para as situações de permeabilidade constante com a porosidade,

permeabilidade em função da porosidade, visando conhecer o comportamento do

fluxo diante de tais situações. Uma comparação qualitativa entre Mohr–Coulomb

e Lade–Kim também será feita.

Os carregamentos aplicados no contorno e a pressão aplicada no poço e na

formação estão descritos na tabela 6.02

Carregamento no contorno (MPa) Pressão (MPa) Carregamento

Hσ 'Hσ hσ '

hσ wp op

5 60 40 45 25 18 20

6 60 40 45 25 15 20

7 60 40 45 25 10 20

8 65 45 45 25 15 20

9 65 45 45 25 10 20

Tabela 6.02 - Carregamento no contorno e pressão aplicada no poço para o Calcário do

Campo de Congro.

6.2.1. Permeabilidade constante

Neste item ilustra–se apenas a tensão principal menor, já que está

desempenha o papel da tensão de confinamento, influenciando o comportamento

da zona rompida tanto na escavação como na produção de um poço.

A figura 6.12 ilustra a tensão principal resultante dos carregamentos 6 e 8,

onde o aumento no carregamento no contorno não produz alterações significativas

no campo de tensão principal menor.

108

Os carregamentos 5 e 7 ilustrados na figura 6.13 mostram que a aplicação

do drawdown induz a tensão principal menor a se tornar mais compressiva, o que

aconteceu de forma idêntica ao descrito na figura 6.06.

As figuras 6.14 e 6.15 mostram uma pequena zona de ruptura. O aumento

do carregamento no contorno (figura 6.14) não modificou a zona de ruptura. A

colocação de valores mais elevados para o drawdown (figura 6.15) também não

alterou significativamente a zona rompida. Como mencionado para o arenito, a

ruptura é mais influenciada pela etapa de escavação.

A figura 6.27 mostra a curva tensão desviadora versus deformação axial

para o calcário. Assim como para o arenito de Rio Bonito, o calcário apresenta

uma boa resistência mecânica, o que justifica a dificuldade de propagação da zona

rompida para o interior da formação.

3σ

Carregamento 6: MPaddMPaMPa hH 5,25',40' === σσ

3σ

Carregamento 8: MPaddMPaMPa hH 5,10',45' === σσ

Figura 6.12 – Campo de tensão principal menor na vizinhança do poço para os casos 6 e

8 com drawdown de 5 Mpa. À esquerda representa–se a escavação e à direita a

produção do poço.

109

3σ

Carregamento 5 e 7 (escavação): MPaMPa hH 25',40' == σσ

3σ

Carregamento 5 (produção): MPaddMPaMPa hH 2,25',40' === σσ

3σ

Carregamento 7 (produção): MPaddMPaMPa hH 10,25',40' === σσ

Figura 6.13 – Campo de tensão principal menor na vizinhança do poço para os casos 5 e

7, ilustrando os gráficos sobre a escavação como indicado e os demais o poço em

produção.

110

S

Carregamento 6: MPaddMPaMPa hH 5,25',40' === σσ

S

Carregamento 8: MPaddMPaMPa hH 5,25',45' === σσ

Figura 6.14 - Zona de ruptura na vizinhança do poço variando o estado de tensão e

drawdown de 2 MPa. À esquerda representa –se a perfuração e à direita a utilização do

poço.

111

S

Carregamento 5 e 6 (escavação): MPaddMPaMPa hH 2,25',40' === σσ

S

Carregamento 5 (produção): MPaddMPaMPa hH 5,25',40' === σσ

S

Carregamento 6 (produção): MPaddMPaMPa hH 10,25',40' === σσ

Figura 6.15 - Zona rompida na vizinhança do poço variando o drawdown e um

carregamento aplicado de MPaeMPa hH 25'40' == σσ . O primeiro gráfico

representa a perfuração e os seguintes a utilização do poço.

A figura 6.17 mostra o campo de deformação volumétrica aplicando-se os

carregamentos 6 e 8. Além do comportamento dilatante durante a ruptura indicada

pelas cores vermelha, amarela e verde, nota–se uma região mais comprimida

próxima ao poço (indicada pela cor azul escura), diferente do que ocorreu com o

112

arenito Rio Bonito. O aumento do carregamento no contorno aumenta

significativamente esta região.

Soares (2000) descreve que o calcário próximo à tensão hidrostática de 40

MPa apresenta o fenômeno conhecido como colapso de poros. Este fenômeno

caracteriza-se pelo aumento sensível da deformação volumétrica sem grandes

acréscimos no estado de tensão hidrostático. O comportamento obtido pela

simulação para este material parece estar relacionado a este evento.

Aplicando-se diferentes drawdown, a deformação volumétrica não se

mostrou tão sensível a esta mudança como ilustrado na figura 6.18. Uma possível

razão para isso pode ser que a diferença de pressão aplicada não tenha sido

suficiente para promover uma alteração na tensão hidrostática na região de maior

deformação compressiva.

Na figura 6.16, plota-se o gráfico de deformação volumétrica versus raio

adimensional na direção de teta igual a zero. Observa-se que a deformação

volumétrica apresenta um valor máximo de compressão próximo ao poço em

seguida assume um comportamento dilatante para uma pequena faixa situada

entre o poço e o valor Maximo de compressão, confirmando o comportamento

descrito nas figuras 6.17 e 6.18.

-6.00

-4.00

-2.00

0.00

2.00

4.00

6.00

8.00

10.00

0 2 4 6 8 10 12 14

Raio adimensional

ε vol

16

Escavação - 40 MPaProdução - 40 MPa, dd = 5 MPaEscavação - 45 MpaProdução - 45 MPa, dd = 5 MPa

Figura 6.16 – Campo de deformação volumétrica na direção θ = 0 e θ = 45 sob uma

tensão horizontal maior de 40 e 45 MPa (carregamentos 6 e 8) e um drawdown de 5

MPa.

113

vε

MPaddMPaMPa hH 5,25',40' === σσ

vε

MPaddMPaMPa hH 5,25',45' === σσ

Figura 6.17 – Campo de deformação volumétrica na vizinhança do poço para os casos 6

e 8. À esquerda está representada a perfuração e à direita a utilização do poço.

114

vε

Carregamento 6 e 7 (escavação): MPaMPa hH 25',40' == σσ

vε

Carregamento 6 (produção): MPaddMPaMPa hH 5,25',40' === σσ

vε

Carregamento 7 (produção): MPaddMPaMPa hH 10,25',40' === σσ

Figura 6.18 – Campo de deformação volumétrica na vizinhança do poço para os casos 6

e 7, ilustrando os gráficos sobre a escavação como indicado e os demais o poço em

produção.

115

6.2.2. Permeabilidade variável

Nesta simulação, a permeabilidade é colocada em função do índice de

vazios, e assim indiretamente com a deformação volumétrica.

A permeabilidade em função do índice de vazios foi obtida utilizando–se os

dados de Soares (2001) relativos ao calcário do Campo de Congro. Com estes

valores ajustou–se a expressão dada por Raghavan e Chin (2002).

Porosidade

(%)

Permeabilidade

(mD)

23,2 0,9

25,1 1,9

31,1 3,9 Tabela 6.03 - Permeabilidade em função da porosidade (Soares, 2001)

A relação de Raghavan e Chin é expressa por: 1

''

j

ii nn

kk

= (6.1)

onde k é a permeabilidade, n’ é a porosidade, j1 é um fator obtido através de

ajuste de curvas de dados experimentais, o subíndice i indica uma referência. O

valor de j1 que mostrou produzir valores próximos dos ensaios realizados por

Soares foi de 3,3. Os resultados ilustrados são a deformação volumétrica, excesso

de poropressão, região de ruptura e gráficos comparativos.

Como ilustrado na figura 6.19, a mudança na permeabilidade não alterou a

forma e a extensão da zona de ruptura em relação à situação de permeabilidade

constante. A ruptura como no item 6.2.1 parece ser governada pela etapa de

escavação e a zona de ruptura tem a mesma dificuldade para se propagar devido à

resistência mecânica do calcário.

O comportamento dilatante ilustrado na figura 6.20 para os carregamentos 6

e 7 limita–se à parede do poço. O aumento da região com maior deformação

volumétrica compressiva também não é significativo com a variação de

drawdown, idêntico ao que foi descrito no item 6.2.1.

116

S

MPaddMPaMPa hH 5,25',40' === σσ

S

MPaddMPaMPa hH 10,25',40' === σσ

Figura 6.19 - Zona rompida na vizinhança do poço para os casos 6 e 7. À esquerda está

representada a escavação e à direita o poço em produção.

117

vε

Carregamento 6 e 7 (escavação): MPaMPa hH 25',40' == σσ

vε

Carregamento 6 (produção): MPaddMPaMPa hH 5,25',40' === σσ

vε

Carregamento 7 (produção): MPaddMPaMPa hH 10,25',40' === σσ

Figura 6.20 – Campo de deformação volumétrica na vizinhança do poço para os casos 6

e 7, ilustrando os gráficos sobre a escavação como indicado e o demais o poço em

produção.

118

Excesso de poro pressão

teconskMPaddMPaMPa hH tan,10,25',40' ==== σσ

Excesso de poropressão

iávelékMPaddMPaMPa hH var,10,25',40' === σσ

Figura 6.21 - Distribuição de poro – pressão para o caso 7 com drawdown de 10 MPa. A

primeira ilustração considera a permeabilidade constante como indicado e a segunda a

permeabilidade como função da porosidade.

A figura 6.21 mostra a rápida variação do excesso de poropressão na direção

de θ igual a zero, o que pode ser associado à diminuição da permeabilidade nesta

região, dado que ocorre uma maior deformação volumétrica compressiva e com

isto a diminuição do índice de vazios.

Para visualizar essa associação de permeabilidade com a deformação

volumétrica, plotam-se as curvas de deformação volumétrica e permeabilidade nas

direções de teta igual a 0º e 45º na figura 6.22. Observa-se que a permeabilidade

apresenta um comportamento análogo à deformação volumétrica, reduzindo seu

valor próximo à região com maior deformação volumétrica compressiva e

119

aumentando na região de deformação volumétrica dilatante, refletindo o aumento

dos poros da rocha durante a ruptura.

Outro ponto a ser observado no gráfico de deformação volumétrica é a

pequena extensão na direção radial que apresenta comportamento dilatante,

conseqüente da pequena zona de ruptura formada ao redor do poço.

-6.00

-4.00

-2.00

0.00

2.00

4.00

6.00

8.00

0 2 4 6 8 10 12 14

Raio adimensionalε v (%

)

16

θ = 0θ = 45

0.00E+00

5.00E-09

1.00E-08

1.50E-08

2.00E-08

2.50E-08

0 2 4 6 8 10 12 14Raio adimensional

Perm

eabi

lidad

e (m

/s)

16

θ = 0 θ = 45

Figura 6.22 – Curva deformação volumétrica-raio adimensional e permeabilidade-raio

adimensional em duas direções na etapa de produção para o calcário do Campo de

Congro

120

6.2.3. Comparações entre os resultados obtidos utilizando os critérios de Mohr–Coulomb e Lade-Kim

Os parâmetros de Mohr–Coulomb para o calcário do Campo de Congro são

obtidos do trabalho de Coelho (2001). Em relação a estes parâmetros, existe a

dúvida se a amostra utilizada na determinação dos parâmetros do critério de

Mohr–Coulomb provêm do mesmo local da amostra que permitiu a determinação

dos parâmetros do modelo Lade–Kim no Campo de Congro. Por isto está análise

tem um caráter mais qualitativo. Os parâmetros utilizados nesta analise foram o

ângulo de atrito de 34º e a coesão de 7,5 MPa.

A simulação utilizou o carregamento 7 do item 6.2.2 para a comparação. Os

resultados ilustrados são a deformação volumétrica, a zona de ruptura, a

distribuição do índice de vazios e o excesso de poro-pressão.

Nota–se pela figura 6.23 que a região dilatada obtida pelo modelo de Mohr–

Coulomb é mais pronunciada que a do modelo Lade–Kim tanto pelo seu valor

como pela sua extensão que ela atinge no poço. Mas a forma do campo de

deformação volumétrica é praticamente idêntica para ambos os modelos.

Embora seja uma análise qualitativa, a figura 6.24 mostra uma diferença

significativa da extensão da zona de ruptura entre os dois critérios. Uma razão

provável para isso é a inclusão da tensão intermediária no modelo de Lade–Kim.

Outro fato observado é que a zona rompida não se altera da etapa de escavação

para a de produção, o que é semelhante ao modelo Lade-Kim.

O índice de vazios representado na figura 6.25 mostra uma distribuição

uniforme para o critério de Mohr – Coulomb, o que não compatibiliza com a

deformação volumétrica (figura 6.23). Isto indicaria que a deformação

volumétrica descrita pelo modelo de Mohr – Coulomb não consegue representar o

colapso de poros.

A distribuição uniforme do excesso de poropressão descrita na figura 6.26

para o modelo de Mohr–Coulomb pode ser explicada exatamente pela não

representação do colapso de poros. Para o modelo Lade–Kim a distribuição do

índice de vazios é semelhante à deformação volumétrica. Assim a região onde

ocorre o menor índice de vazios tem a variação mais rápida do excesso de

poropressão.

121

vε vε

Modelo Lade–Kim: iávelékMPaddMPaMPa hH var,10,25',40' === σσ

vε vε

Modelo Mohr-Coulomb: iávelékMPaddMPaMPa hH var,10,25',40' === σσ

Figura 6.23 – Campo de deformação volumétrica na vizinhança do poço para os modelos

de Mohr - Coulomb e Lade - Kim para o caso 7 com um drawdown de 10 MPa. O gráfico

a esquerda corresponde à escavação e a direita o poço em produção.

122

S

Modelo Lade–Kim: iávelékMPaddMPaMPa hH var,10,25',40' === σσ

S

Modelo Mohr-Coulomb: iávelékMPaddMPaMPa hH var,10,25',40' === σσ

Figura 6.24 - Zona de ruptura na vizinhança do poço para o caso 7 e um drawdown de

10 MPa. O gráfico a esquerda corresponde à escavação e a direita o poço em produção.

123

Indice de vazios

Modelo Lade–Kim: iávelékMPaddMPaMPa hH var,10,25',40' === σσ

Indice de vazios

Modelo Mohr-Coulomb: iávelékMPaddMPaMPa hH var,10,25',40' === σσ

Figura 6.25 – Distribuição do índice de vazios na vizinhança do poço para os modelos de

Mohr-Coulomb e Lade-Kim para o caso 7 com um drawdown de 10 MPa. O gráfico a

esquerda corresponde à escavação e a direita a utilização do poço.

124

Excesso de poropressão

Modelo Lade – Kim: iávelékMPaddMPaMPa hH var,10,25',40' === σσ

Excesso de poropressão

Modelo Moh-Coulomb: iávelékMPaddMPaMPa hH var,10,25',40' === σσ

Figura 6.26 - excesso de poropressão na vizinhança do poço para os modelos de Mohr-

Coulomb e Lade-Kim com um drawdown de 10 MPa, mantendo um estado de tensão de

MPaeMPa hH 2540 == σσ . O primeiro gráfico corresponde ao modelo de Lade–

Kim e o último o de Mohr – Coulomb.

125

-80

-70

-60

-50

-40

-30

-20

-10

00.00 2.00 4.00 6.00 8.00 10.00 12.00 14.00

ε11(%)

σ 11

- σ33

(MPa

)

Lade uniaxialLade 10 MPaMohr - Coulomb uniaxialMohr - Coulomb 10 MPa

Figura 6.27 – Simulação de ensaios uniaxial e triaxial de compressão para o modelo de

Mohr – coulom e Lade – Kim para o calcário do Campo de Congro.

7 Conclusões e sugestões para trabalhos futuros

O presente trabalho buscou associar mecanismos de ruptura a produção de

areia de rochas reservatório. Para isto, implementou–se o modelo elastoplástico de

Lade–Kim no programa de elementos finitos ABAQUS, que é capaz de simular o

acoplamento hidro–mecânico.

A capacidade do programa em simular o acoplamento hidro–mecânico foi

verificada através da comparação com a solução analítica de Detournay e Cheng

(1993) para a coluna poroelástica e Detournay e Cheng (1988) para o caso de um

poço vertical. Como mostrado no item (3.2.6.1), o programa ABAQUS conseguiu

reproduzir satisfatoriamente o comportamento da coluna poroelástica tanto para o

meio com constituintes compressíveis como para incompressíveis.

Para o caso do poço, a reprodução não teve a mesma eficiência. Porém o

próprio autor da solução analítica afirma que existe uma falha, no caso para uma

das tensões. Mas como pode ser visto no item 3.2.6.2, para as demais variáveis

testadas a simulação produziu resultados próximos da solução analítica. Isto

permitiu verificar com certa segurança o efeito da dependência da permeabilidade

com o estado de tensões, além da influência do fluxo de fluido no comportamento

mecânico da rocha.

O algoritmo do tipo explícito utilizado para integrar a relação tensão versus

deformação foi eficiente, reproduzindo satisfatoriamente os ensaios de Barroso

(2001) fora do ambiente do ABAQUS.

Introduzido o modelo Lade–Kim no programa ABAQUS, a simulação do

poço apresentou problemas de convergência na etapa de amolecimento para o

arenito Rio Bonito. Durante a análise, instabilidades como autovalores negativos

na matriz de rigidez foram observados, o que não permitia o término da

simulação.

Duas alternativas foram utilizadas para contornar tal problema. A primeira,

reduzir o incremento de carga, o que onerou o esforço computacional. A segunda,

127

restringir o valor da variável do modelo na função de escoamento que controla o

decaimento da resistência, o que pode ter influenciado no decréscimo da coesão

do material.

Os resultados do capítulo 6 mostraram uma zona rompida na vizinhança do

poço insignificante comparada aos breakouts descritos na literatura (Goodman

(1989), Cook (1990)). Assim como a extensão da região dilatante é ínfima.

A utilização de rochas com boa resistência à compressão dificultou a

analogia do mecanismo de ruptura com a produção de sólidos, normalmente

associada a rochas brandas. Outro fato que prejudicou a analogia e a que não se

chegou a uma conclusão foi à tensão principal menor se tornar mais compressiva

com o aumento do drawdown, o que limitou o efeito de arraste descrito no

processo de produção de sólidos. Portanto, a analogia não obteve o êxito esperado.

A introdução da permeabilidade dependente do estado de tensão mostrou a

influência do acoplamento hidro–mecânico. Embora não se tenha a informação da

vazão, que permitiria verificar facilmente as mudanças no fluxo, ela permitiu

visualizar a influência da deformação volumétrica no excesso de poro–pressão.

A discussão entre qual modelo é mais apropriado para simular o processo de

produção de sólidos é qualitativa. A primeira causa é não dispor de parâmetros

para os dois modelos de um mesmo material. O segundo fato é que para o modelo

Mohr–Coulomb uma série de modificações pode ser feita e não estão presentes no

programa ABAQUS, por exemplo, o módulo de Young dependente do estado de

tensão.

Independente disto notou–se para o caso do calcário, que a zona de ruptura

dada pelo critério de Mohr–Coulomb é significativamente maior do a do critério

Lade–Kim, o que colocaria uma maior área susceptível a produzir sólidos. Isto

possivelmente está relacionado à influência da tensão intermediária na resistência

mecânica.

Outro ponto de interesse é a “aparente” capacidade do modelo Lade–Kim

em considerar o colapso de poros, o que influenciaria no fluxo de fluído. Evento

esse que não foi descrito pelo modelo de Mohr–Coulomb. Deste ponto de vista, o

modelo Lade–Kim está mais capacitado a simular situações que levem em conta o

acoplamento fluído–mecânico.

128

Para trabalhos futuros sugere–se uma forma de se colocar a resistência

residual no modelo Lade–Kim, talvez pela introdução de um novo parâmetro. Mas

que esta iniciativa esteja em conjunto com um trabalho experimental.

O esforço computacional grande devido à dificuldade em simular os ensaios

dentro do ambiente do ABAQUS deixa uma questão, há como diminuí–lo? A

aplicação de um algoritmo do tipo implícito é válida no sentido de se responder a

essa questão. Assim como verificar a eficiência de ambos os algoritmos.

Os parâmetros do modelo Lade – Kim obtidos através de ensaios triaxiais

devem ser verificados através de ensaios em células triaxiais a fim de verificar se

conseguem reproduzir o comportamento dos materiais em questão para qualquer

estado de tensão.

A análise pelo método de elementos finitos não conseguiu reproduzir o

efeito de arraste, uma possibilidade seria então considerá – la como uma análise

preliminar, e em seguida aplicar a resposta obtida em uma análise por meio

discreto.

8 Referências Bibliográficas

ABAQUS–Theory Manual, version 6.2.

Abass, H. H.; Habbtar, A. H.; Shebatalham, A. Sand control during drilling,

perforation, completion and production. SPE 81492, Society of Petroleum

Engineers, Bahrain, june, 2003.

Barroso, E. V. Avaliação de um modelo elastoplástico para estudos de

processo de produção de areia em rochas produtoras de petróleo. Rio de

Janeiro, 2002. 188p. Tese de Doutorado. DEC, Puc, Rio.

Biot, M. A. General solution of the equations of elasticity and consolidation

for a porous material. Applied Mechanics Division, paper 55–A7, 1955.

Biot, M. A. General Theory of three – dimensional consolidation. Journal of

applied Physics. Vol. 12(155), p. 155–164, 1941.

Boyce, W; Diprima R. C. Equacões diferenciais elementares e problemas de

valore de contorno. Livros Técnicos e Científicos. Rio de Janeiro, RJ, 1998.

532p.

Coelho, L. C. Modelos de ruptura de poços de petróleo. Rio de Janeiro, 2001.

178 p. Tese de Doutorado. COPPE/UFRJ.

Cook, N. G. W., Ewy, R. T. Deformation and Fracture around cylindrical

openings in rock – I. Observations and analysis of deformations.

International Journal of Rock Mechanics and Mining Sciences & Geomechanics

Abstracts, vol. 27(5), p. 387–407, 1990.

130

Desai, C. S.; Siriwardane, H. J. Constitutive laws for engineering materials

with emphasis on Geologic Materials. Prentice–Hall, Inc., Englewood Cliffs,

New Jersey, 1984. 468 p.

Detournay, E.; Cheng, H–D. Fundamentals of poroelasticity. In

Comprehensive Rock Engineering (Principles, practice and projects),

Pergamon Press, Great Britain, 1993, vol. (2), p. 113–171.

Detournay, E.; Cheng, H–D. Poroelastic Response of a borehole in a

hydrostatic stress field. International Journal of Rock Mechanics and Mining

Sciences & Geomechanics abstracts. Vol. 25 (3), p. 171–182, 1988.

Dormand, J.R.; Prince, P. J. A family of embedded Runge-Kutta formulae.

Journal of Computation and Applied Mathematics. Vol. 6, p. 19–26, 1980.

Dusseault, M. B.; Santarellli, F. J. A conceptual model for massive solids

production in poorly – consolidated sandstones. Rock at Great Depth, Maury

& Fourmaintrauux (eds), Balkema, Rotterdam, 1989.

Ferreira, F. H. Uma implementação numérica para a solução de problemas

de poroelasticidade. Rio de Janeiro, 1996. 155 p. Dissertação de Mestrado.

DEC, PUC, Rio.

Fjaer, E. et al. Petroleum related rock mechanics. Elsevier Science

Publishers, N. Y. 1992. p. 257-268.

Jeremic, B.; Runesson, K.; Sture, S. A model for elastic–plastic pressure

sensitive materials subjected to large deformations. International Journal of

Solids and Structures, vol. 37 (22), p. 3079–3100, 2000.

Jeremic, B.; Yang. Z. Template elastic–plastic computations in

geomechanics. International Journal for Numerical and Analytical Methods in

Geomechanics, vol. 26, p. 1407–1427, 2002.

131

Lade, P. V., Kim, M. K. Single hardening constitutive model for soil, rock

and concrete. International Journal of Solids Structures, vol. 32 (14), p. 1963–

1978, 1995.

Lade, P. V.; Kim, M. K. Modeling rock strength in three dimensions.

International Journal of Rock Mechanics and Mining Sciences & Geomechanics

Abstracts, vol. 21 (1), p. 21-33, 1984.

Lade, P. V.; Kim, M. K. Single Hardening constitutive model for frictional

materials, I–Plastic potential function. Computers and Geotechnics, vol. 5, p.

307–324, 1988.

Lade, P. V.; Kim, M. K. Single Hardening constitutive model for frictional

materials, II–Yield criterion and plastic work contours. Computers and

Geotechnics, vol. 6, p. 13-29, 1988.

Lade, P. V.; Kim, M. K. Single Hardening constitutive model for frictional

materials, III–Comparisons with experimental data. Computers and

Geotechnics, vol. 6, p. 31 - 47, 1988.

Mclean, M. R.; Addis, M. A. Wellbore stability analysis: A review of current

methods of analysis and their field application. SPE 19941, Society of

Petroleum Engineers, Houston, Texas, February, 1990.

Melo, L. T. B. Utilização de um modelo elastoplástico para análise de

deformações em solos. Rio de Janeiro, 1995. 99 p. Dissertação de Mestrado.

DEC, PUC, Rio.

Morita, N. Field and laboratory verification of sand production prediction

models. SPE 27341, Society of Petroleum Engineers, Lafayette, February 1994.

Morita, N.; Boyd, P. A. Typical sand production problems: case studies and

strategies for sand control. SPE 22739, Society of Petroleum Engineers,

132

Dallas, Texas, October 1991.

Morita, N; Whitfill, D. L. Realistic Sand Production Prediction Numerical

Approach. SPE 16989, Society of Petroleum Engineers, Dallas, Texas,

September, 1987

.

Nogueira, C. L. Análise não linear de escavações e aterros. Rio de Janeiro,

1998. 250p. Tese de Doutorado. DEC, PUC, Rio.

Raghavan, R.; Chin, L. Y. Productivity Changes in Reservoirs with stress–

dependent permeability. SPE 77535, Society of Petroleum Engineers, San

Antonio, Texas, October 2002.

Rice, J. R. & Cleary, M. P. Some basics diffusion solutions for fluid–

saturated elastic porous media with compressible constituents. Reviews of

Geophysics and Space Physics, vol. 14 (2), p. 227–241, 1976.

Sloan, S. W. Substepping schemes for the numerical integration of

elastoplastic stress–strain relations. International Journal for Numerical

Methods in Engineering, vol. 24, p. 893–911, 1987.

Sloan, S. W.; Abbo, A.J. Sheng, D. Refined explicit integration of

elastoplastic models with automatic error control. Engineering

Computations, vol. 18(1/2), p.121-151, 2001.

Soares, A. C. Avaliação do colapso de Poros no Campo de Congro. Rio de

Janeiro: PDP/TEP, 2001. 41p. Relatório Técnico.

Soares, A. C. Um estudo experimental para definição do colapso de poros

em rochas carbonáticas. Rio de Janeiro, 2000. 185 p. Dissertação de

Mestrado, COPPE/UFRJ, RJ.

Stehfest, H. Numerical Inversion of Laplace Transform. Communications of

the ACM. vol. 13, p. 47 – 49, 1970.

133

Thomas, J. E. Fundamentos da Engenharia do Petróleo. Rio de Janeiro, RJ,

2001, 271 p.

Tronvoll, J. et al. Hydrodinamic erosion: a potential mechanism of sand

production in weak sandstones. International Journal of Rock Mechanics &

Mining Sciences, vol. 34 (3/4), paper 292.

Tronvoll, J.; Skjaerstein A.; Papamichos, E. Sand Production: Mechanical

Failure or hydrodynamic erosion? International Journal of Rock Mechanics

& Mining Sciences, vol. 34(3/4), paper 291.

Veeken, C. A. M. et al. Sand Production Prediction Review: Developing and

Integrate Approach. SPE 22792, Society of Petroleum Engineers, Dallas,

Texas, October, 1991.

Willam, K. J.; Menétrey, Ph. Triaxial failure criterion for concrete and its

generalization. ACI Structural Journal, vol. 92 (3), p. 311–318, 1995.