MODELAGEM SÍSMICA EM MEIOS COMPLEXOS
Eldues Oliveira Martins
TESE SUBMETIDA AO CORPO DOCENTE DA COORDENAÇÃO DOS
PROGRAMAS DE PÓS-GRADUAÇÃO DE ENGENHARIA DA UNIVERSIDADE
FEDERAL
DO
RIO
DE
JANEIRO
COMO
PARTE
DOS
REQUISITOS
NECESSÁRIOS PARA A OBTENÇÃO DO GRAU DE MESTRE EM CIÊNCIAS EM
ENGENHARIA CIVIL.
Aprovada por:
Prof. Luiz Landau, D.Sc.
Dr. Djalma Manoel Soares Filho, D.Sc.
Prof. José Luis Drummond Alves, D.Sc.
Prof. Liliana Alcazar Diogo, D.Sc.
Prof. Wilson Mouzer Figueiró, D.Sc.
RIO DE JANEIRO, RJ - BRASIL
MAIO - 2003
MARTINS, ELDUES OLIVEIRA
Modelagem
Sísmica
em
Meios
Complexos [Rio de Janeiro] 2003.
XVI, 132p. 29,7 cm (COPPE/UFRJ,
M.Sc., Engenharia Civil, 2003)
Tese – Universidade Federal do Rio
de Janeiro, COPPE
1. Geofísica;
2. Modelagem Sísmica;
3. Anisotropia;
I. COPPE/UFRJ II. Título (série).
ii
“À minha querida e paciente esposa Luciana, que me deu meus muito amados
filhos Pedro, Mariana e João, pelo seu amor, suporte, atenção e compreensão
que foram fundamentais para a conclusão desta tese; e enfim,
a todos que amo e admiro.”
iii
AGRADECIMENTOS
Agradeço ao coordenador e orientador Prof. Luiz Landau e a todo corpo técnicoadministrativo do Laboratório de Métodos Computacionais em Engenharia – LAMCE pela
oportunidade da realização deste projeto.
Agradeço ao meu orientador Dr. Djalma M. Soares Filho pelo seu incansável esforço
no desenvolvimento da pesquisa científica, pela sua brilhante dedicação na orientação desta
tese e valorosa contribuição na minha formação profissional.
Agradeço ao geolólogo e amigo Ricardo Bedregal pelas proveitosas discussões sobre a
geologia do petróleo, pelos incentivos e apoio incondicional durante a realização desta tese.
Agradeço ao Prof. José Luis Drummond Alves e ao colega Dênis Araujo Figueiras de
Souza pelas contribuições na otimização visualização e desenvolvimento dos algoritimos em
Fortran.
Agradeço aos sempre amigos Josias J. da Silva, Alexandre Lopes, Ricardo Alencar,
Luis Fernando, Magda Almada, Mônica Caruso e a todos os outros que com certeza, sem eles,
a caminhada seria muito mais difícil.
Agradeço ao geofísico da PETROBRAS/CENPES/PDEXP/GEOF João Cláudio
Conceição da Silva pela colaboração para a realização deste trabalho.
iv
Resumo da Tese apresentada a COPPE/UFRJ como parte dos requisitos necessários para a
obtenção do grau de Mestre em Ciências (M. Sc.)
MODELAGEM SÍSMICA EM MEIOS COMPLEXOS
Eldues Oliveira Martins
Maio/2003
Orientadores: Djalma Manoel Soares Filho
Luiz Landau
Programa: Engenharia Civil.
Foram desenvolvidas e implementadas duas metodologias para simulações numéricas
de levantamentos sísmicos em meios complexos através da solução do sistema de equações
que regem o comportamento ondulatório em meios transversalmente isotrópicos
bidimensionais, usando o método das diferenças finitas. Para meios terrestres, foi
generalizado o algoritmo proposto por Zahradnik e Priolo [1994], no qual os componentes do
campo de deslocamento das partículas são propagados e os parâmetros elásticos são
introduzidos por integrações ao longo das linhas definidas pela malha, suposta regular.
Utilizam-se aproximações de segunda ordem para as derivadas parciais na geração do
operador de propagação. Para abordagens marinhas, modificou-se o algoritmo apresentado em
Levander [1986], no qual os campos de velocidade e tensão, assim como os parâmetros
elásticos, são definidos em malhas intercaladas regulares. Para este método, utilizam-se
aproximações de segunda ordem para as derivadas temporais e quarta ordem para as
espaciais, na geração dos operadores que realizam a propagação dos campos de velocidades e
tensão. Os métodos não estão restritos a qualquer tipo de geometria de aquisição (fontes e
receptores podem estar posicionados em qualquer lugar no modelo), assim como não há
qualquer limitação no que diz respeito ao modelo geológico, desde que as condições que
garantem a estabilidade e ausência de dispersões numéricas sejam satisfeitas. No que diz
respeito às aplicações, inicialmente, foi estudado o comportamento ondulatório em função do
grau de anisotropia TIV. Especificamente, verificou-se em instantâneos, obtidos em
modelagens para meios homogêneos, a influência dos parâmetros de Thomsen [1986] na
propagação de ondas compressionais e cisalhantes. Em seguida, foram obtidos sismogramas e
instantâneos em simulações sísmicas em dois modelos com altos contrastes de velocidades,
típicos de bacias brasileiras. Os resultados foram analisados e comparados com os obtidos em
simulações elásticas isotrópicas.
v
Abstract of Thesis presented to COPPE/UFRJ as a partial fulfillment of the requirements for
the degree of Master of Science (M.Sc.)
SEISMIC MODELING ON COMPLEX MEDIA
Eldues Oliveira Martins
Maio/2003
Advisors: Djalma Manoel Soares Filho
Luiz Landau
Departament: Civil Engineering.
Two finite differences type methods are introduced for seismic modeling on
transversely isotropic media. The first one is designed for seismic acquisition simulations on
on-shore geologic models and consists of a modification of the Zahradník and Priolo (1994)
proposal to accommodate seismic modeling on anisotropic media. In this method, the
components of the displacement vector are discretized in a regular grid and the elastic tensor
components are introduced through vertical and horizontal integrations along the grid lines.
Second order approximations are used for spatial and temporal derivatives and Vacuum
Formalism, as it is posed in Zahradník, Moczo and Heron (1993), is used as boundary
condition along the upper interface. This technique is not restricted to flat observation
surfaces, but it cannot be used on off-shore models, since it gives rise to instabilities where
shear modulus is small. The second method is based on Levander (1986) proposal technique
and is indicated to seismic simulations on off-shore geologic models. As it was done with the
Zahradník and Priolo (1986) method, we also adapt this second one to cope with transversely
isotropic media. In this case, the components of the velocity field and the elastic parameters
are defined in a staggered grid and the conditions for preventing spurious numerical
dispersions do not depend on shear modulus value. This fact makes off-shore seismic
simulations feasible using this approach. In relation to partial derivatives approximations, the
temporal derivatives are approximated in second order whereas the spatial ones are fourthorder approximated. Both methods are not restricted to any particular acquisition pattern
(sources and receivers may be located anywhere) and may be applied to any geologic case (as
far as the conditions for preventing spurious dispersions are satisfied). To illustrate, we show
snapshots and seismograms obtained on several simulations using the proposed methods and
indicate some of the characteristics of the wave field propagation on realist transversely
isotropic media.
vi
Índice
Página de Assinaturas . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .
i
Ficha Catalográfica . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .
ii
Dedicatória . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .
iii
Agradecimentos . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .
iv
Resumo . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .
v
Abstract . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .
vi
Índice . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . vii
Índice de Figuras . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .
x
Índice de Tabelas . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . xvi
1 Introdução
1
1.1 Elasticidade do Meio . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .
5
1.2 Anisotropia do Meio . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .
5
1.3 Heterogeneidade do Meio . . . . . . . . . . . . . . . . . . . . . . . . . . . .
8
2 Modelagem Sísmica para Meios Elásticos TIV: Caso Terrestre
10
2.1 Desenvolvimento . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 12
2.1.1
Discretização das Equações da Elastodinâmica para meios TIV . . . 12
vii
2.1.2
Condições de Contorno para as Bordas Laterais e Inferior . . . . . . 34
2.2 Termo Fonte . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 35
2.3 Definição da malha . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 37
2.3.1
Condições para Não Ocorrência de Dispersão Numérica . . . . . . . 37
2.3.2
Condições de Estabilidade Numérica . . . . . . . . . . . . . . . . . 38
2.4 Metodologia . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 38
2.5 Aplicações . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 42
2.5.1
Primeira Aplicação: Propagação de Onda em Meios Homogêneos
Anisotrópicos nos Modelos Propostos em Thomsen (1986) . . . . . 42
2.5.2
Segunda Aplicação: Modelo Formado por Duas Camadas Separadas
por Interface Curva . . . . . . . . . . . . . . . . . . . . . . . . . . . 50
2.5.3
Terceira Aplicação: Modelo com Altos Contrastes de Impedância . 58
3 Modelagem Sísmica para Meios Elásticos TIV: Caso Marinho
66
3.1 Desenvolvimento . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 67
3.1.1
Discretização pelo Método de Diferenças Finitas . . . . . . . . . . . 69
3.2 Metodologia . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 76
3.3 Aplicações . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 79
3.3.1
Primeira Aplicação: Modelagem em Meios Acústicos com operadores Elásticos . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 79
3.3.2
Segunda Aplicação: Modelagem Sísmica em Modelo Típico de Margens Passivas . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 81
viii
Conclusão
92
Referências
96
Apêndices
101
A Tensor de Elasticidade (Cijkl ) Expresso no Sistema de Coordenadas Principal
102
A.1 Componentes da tensão e deformação . . . . . . . . . . . . . . . . . . . . . 103
B Discretização por Diferenças Finitas
113
B.1 Malhas Intercaladas . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 113
ix
Lista de Figuras
1.1 Direção de polarização e propagação definida para um caso anisotrópico de
meio com finas camadas paralelas. . . . . . . . . . . . . . . . . . . . . . . .
6
1.2 Modelos litológicos equivalentes para diferentes tipos de simetrias. . . . . .
7
2.1 Linhas da malha onde são calculadas as médias geométricas dos parâmetros. 16
2.2 Modelo da malha do método PS2 de diferenças finitas para aproximação
das derivadas em relação a uma única variável para pontos da malha na
superfície. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 17
2.3 Modelo da malha do método PS2 de diferenças finitas para aproximação
das derivadas em relação a uma única variável para pontos da malha no
interior do modelo. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 18
2.4 Modelo da malha do método PS2 de diferenças finitas para aproximação das
derivadas em relação às variáveis x e z para pontos da malha na superfície
SE. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 25
2.5 Modelo da malha do método PS2 de diferenças finitas para aproximação das
derivadas em relação às variáveis x e z para pontos da malha na superfície
SW. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 26
x
2.6 Modelo da malha do método PS2 de diferenças finitas para aproximação das
derivadas em relação às variáveis x e z para pontos da malha na superfície
NE e NW. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 27
2.7 Modelo da malha do método PS2 de diferenças finitas para aproximação
das derivadas em relação às variáveis x e z para pontos da malha no interior
WN, WS, EN e ES. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 28
2.8 Gráfico da função fonte . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 36
2.9 Diagrama representando a 1a metodologia para modelagem em meios elásticos TIV. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 41
2.10 Instantâneos do componente vertical no tempo de 3 ms. a) alumínio-lucita.
b) folhelho anisotrópico. . . . . . . . . . . . . . . . . . . . . . . . . . . . . 44
2.11 Instantâneos do componente vertical no tempo de 3 ms. a) apatita. b)
biotita. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 44
2.12 Instantâneos do componente vertical no tempo de 3 ms. a) calcita. b)
arenito calcarenoso. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 45
2.13 Instantâneos do componente vertical no tempo de 3 ms. a) folhelho anisotrópicoss. b) folhelho argiloso. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 45
2.14 Instantâneos do componente vertical no tempo de 3 ms. a) calcáreo síltico.
b) arenito-folhelho. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 46
2.15 Instantâneos do componente vertical no tempo de 3 ms. a) arenito. b)
cristal de quartzo. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 46
2.16 Instantâneos do componente vertical no tempo de 3 ms. a) folhelho com
óleo. b) folhelho. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 47
xi
2.17 Instantâneos do componente vertical no tempo de 3 ms. a) folhelho anisotrópicols. b) calcáreo-folhelho. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 47
2.18 Instantâneos do componente vertical no tempo de 3 ms. a) siltito laminado.
b) arenito imaturo. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 48
2.19 Instantâneos do componente vertical no tempo de 3 ms. a) gelo b) material
de gipso intemperizado. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 48
2.20 Instantâneos do componente vertical no tempo de 3 ms. a) gas-arenitoágua-arenito. b) cristal. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 49
2.21 Instantâneos do componente vertical no tempo de 3 ms. a) cinza vulcânica. 49
2.22 Modelo de velocidade abordado na segunda aplicação. . . . . . . . . . . . . 51
2.23 Instantâneos do componente vertical do campo de onda, considerando fonte
explosiva na posição central do modelo homogêneo, VP =3000 m/s e VS =1732
m/s. a) Isotrópico. b) anisotrópico ε = δ = 0, 1. c) anisotrópico ε = δ =
0, 2. d) anisotrópico ε = δ = 0, 3. . . . . . . . . . . . . . . . . . . . . . . . 52
2.24 Instantâneos do componente horizontal do campo de onda, considerando
fonte explosiva na posição central do modelo homogêneo, VP =3000 m/s e
VS =1732 m/s. a) Isotrópico. b) anisotrópico ε = δ = 0, 1. c) anisotrópico
ε = δ = 0, 2. d) anisotrópico ε = δ = 0, 3. . . . . . . . . . . . . . . . . . . . 53
2.25 Sismogramas relativos ao componente vertical do campo de ondas. . . . . . 54
2.26 Sismogramas relativos ao componente horizontal do campo de ondas. . . . 55
2.27 Sismogramas relativos ao componente vertical do campo de ondas, sem a
presença das ondas superficiais. . . . . . . . . . . . . . . . . . . . . . . . . 56
xii
2.28 Sismogramas relativos ao componente horizontal do campo de ondas, sem
a presença das ondas superficiais. . . . . . . . . . . . . . . . . . . . . . . . 57
2.29 Modelo de velocidades da bacia sedimentar terrestre estudada. . . . . . . . 59
2.30 Distribuição dos parâmetros de Thomsen no modelo para o caso terrestre. . 59
2.31 Instantâneos relativos ao componente vertical para os modelos isotrópico
(à esquerda) e anisotrópico nos tempos de propagação 0,5; 0,75; 1,00 e 1,25 s. 61
2.32 Instantâneos relativos ao componente horizontal para os modelos isotrópico
(à esquerda) e anisotrópico nos tempos de propagação 0,5; 0,75; 1,00 e 1,25 s. 62
2.33 Sismogramas relativos ao componente vertical para os modelos isotrópico (à
esquerda) e anisotrópico (ε = δ = 0, 2, na camada que apresenta Vp = 2500
m/s). . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 63
2.34 Sismogramas relativos ao componente horizontal para os modelos isotrópico
(à esquerda) e anisotrópico (² = δ = 0, 2, na camada que apresenta Vp =
2500 m/s). . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 64
2.35 Ondas superficiais, especificamente as ondas diretas P e S, e as ondas
Rayleigh. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 65
3.1 Instantâneos relativos ao componente vertical, ilustrando os problemas de
dispersão numérica. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 67
3.2 Diagrama de velocidade. Malha intercalada de diferenças finitas e diagrama
espacial para as atualizações das velocidades. . . . . . . . . . . . . . . . . . 70
3.3 Diagrama de tensões. Malha intercalada de diferenças finitas e diagrama
espacial para as atualizações de tensões. . . . . . . . . . . . . . . . . . . . 72
xiii
3.4 Diagrama representando a 2a metodologia para modelagem em meios elásticos TIV. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 78
3.5 Instantâneos relativos aos componentes horizontal e vertical do campo de
velocidade (derivada do campo de deslocamento), U e V, aos componentes
do tensor de esforços τ xx , τ xz e τ zz , respectivamente, considerando um meio
acústico (µ = 0). . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 80
3.6 Modelo de velocidades da bacia sedimentar marítima estudada com os valores de velocidade da onda P em cada camada. . . . . . . . . . . . . . . . 82
3.7 Modelo de densidade no modelo para o caso marinho. . . . . . . . . . . . . 82
3.8 Distribuição dos parâmetros de Thomsen no modelo para o caso marinho. . 83
3.9 Instantâneos relativos ao componente vertical do campo de velocidades
obtidos na simulação sísmica no modelo para o caso marinho. . . . . . . . 84
3.10 Instantâneos relativos ao componente horizontal do campo de velocidades
obtidos na simulação sísmica no modelo para o caso marinho. . . . . . . . 85
3.11 Instantâneos relativos ao componente τ xx do campo de tensões obtidos na
simulação sísmica no modelo para caso marinho. . . . . . . . . . . . . . . . 86
3.12 Instantâneos relativos ao componente τ xz do campo de tensões obtidos na
simulação sísmica no modelo para o caso marinho. . . . . . . . . . . . . . . 87
3.13 Instantâneos relativos ao componente τ zz do campo de tensões obtidos na
simulação sísmica no modelo para o caso marinho. . . . . . . . . . . . . . . 88
3.14 Sismogramas relativos ao componente vertical do campo de velocidades
obtidos na simulação sísmica no modelo para o caso marinho. . . . . . . . 89
xiv
3.15 Sismogramas relativos ao componente horizontal do campo de velocidades
obtidos na simulação sísmica no modelo para o caso marinho. . . . . . . . 89
3.16 Sismogramas relativos ao componente τ xx do campo de tensões obtidos na
simulação sísmica no modelo para o caso marinho. . . . . . . . . . . . . . . 90
3.17 Sismogramas relativos ao componente τ xz do campo de tensões obtidos na
simulação sísmica no modelo para o caso marinho. . . . . . . . . . . . . . . 90
3.18 Sismogramas relativos ao componente τ zz do campo de tensões obtidos na
simulação sísmica no modelo para o modelo marinho. . . . . . . . . . . . . 91
xv
Lista de Tabelas
2.1 Parâmetros da malha usada nas modelagens na primeira aplicação para o
caso terrestre. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 42
2.2 Parâmetros de Thomsen usados na primeira aplicação para o caso terrestre. 43
2.3 Parâmetros da malha usada nas modelagens na segunda aplicação para o
caso terrestre. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 50
2.4 Parâmetros da malha usada nas modelagens na terceira aplicação para o
caso terrestre. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 60
3.1 Parâmetros da malha usada nas modelagens na primeira aplicação para o
caso marinho. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 79
3.2 Parâmetros da malha usada nas modelagens na segunda aplicação para o
caso marinho. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 81
xvi
Capítulo 1
Introdução
A modelagem sísmica tem sido uma das ferramentas mais importantes à disposição dos profissionais envolvidos nos processos de exploração e explotação de hidrocarbonetos através de métodos sísmicos. De fato, de posse de um conhecimento satisfatório
da resposta sísmica de cada componente do sistema petrolífero, intérpretes (geólogos,
geofísicos e engenheiros) podem trabalhar com mais segurança; através de simulações
numéricas de aquisições sísmicas, os geofísicos responsáveis pela definição dos parâmetros
de levantamento podem otimizar os custos envolvidos com o mesmo; e os profissionais
engajados no desenvolvimento de novas tecnologias de processamento e interpretação podem inferir a efetividade de cada técnica com aplicações em dados simulados, explorando
o conhecimento exato do modelo geológico que os geraram. Com relação às ferramentas
atualmente utilizadas no processamento de dados de campo, a modelagem encontra-se no
cerne dos algoritmos. Por exemplo, na migração pré-empilhamento em profundidade, a
solução da equação da onda é necessária nas etapas de depropagação do campo de onda
registrado e na geração das matrizes de tempo de trânsito.
1
O contínuo e acelerado desenvolvimento dos métodos de aquisição, processamento e interpretação sísmica têm permitido um aumento da recuperação e produção
das reservas petrolíferas. Na última década a exploração, explotação e caracterização de
reservas em estruturas e estratigrafias complexas e em condições especiais como fraturamentos, por exemplo, tem sido intensificado, requisitando assim a necessidade de se
desenvolver novas tecnologias capazes de se reconstruir modelos adequados de subsuperfície e determinando propriedades específicas, como a orientação e densidade de fraturas,
porosidade, permeabilidade, presença e saturação de fluidos, etc [1].
Uma propriedade fundamental de subsuperfície é a anisotropia, que atualmente na indústria raramente é levada em consideração em função do dispêndio computacional, que aumenta os custos de processamento, na inversão dos parâmetros elásticos e
migração. Se entende que a variação das propriedades físicas de acordo com a orientação
na qual se realiza a medida para a modelagem geológica é de grande importância, como
por exemplo, no fenômeno de birrefrigência das ondas de cisalhamento, no qual a onda
sísmica se divide em duas frentes de onda que se propagam com velocidades diferentes (o
movimento das partículas na onda mais rápida é paralela à direção de fraturamento do
meio).
Os dados de onda convertida implicam um novo paradigma de processamento, pois as reflexões sobre eventos em subsuperfície ocorrem de forma totalmente
assimétrica com respeito às reflexões de um só modo de propagação. Estas reflexões
dependem em grande parte do quociente de velocidade das ondas P e S, assim como
da profundidade do objetivo, no qual aumenta a complexidade do processamento para
eventos em subsuperfície.
2
O começo da sismologia usando diferenças finitas para resolver problemas
com propagação de onda tem aproximadamente 50 anos. Geologias complexas, provenientes de dobramentos e falhamentos decorrentes de variações de velocidades laterais
fortes e relevos de topografias severas são os principais responsáveis por dificultar o imageamento sísmico em subsuperfície fidedignamente (Soares Filho et al) [13].
Problemas modelando geofísica de exploração quase sempre dividem-se em
estruturas complexas ou estratigrafias complexas caracterizados por interfaces irregulares
separando unidades litológicas maiores. Então para ser efetivo, um algoritmo de modelagem numérica deve automaticamente aproximar as condições de contorno através destas
interfaces. Cada um algoritmo heterogêneo de diferenças finitas soluciona um sistema de
equações de segunda ordem em uma malha singular para malha normal, aplicado usualmente em modelagem terrestre, que foi introduzido por Alford (1974) [14] e estendido por
Kelly et al (1976) [15] ou quarta ordem em malha intercalada, aplicado usualmente em
modelagem marítma, que foi introduzido por Levander (1990) [4].
A técnica de diferenças finitas particularmente é atrativa para geometrias em
subsuperfícies estruturalmente complexas porque as condições de contorno dos contatos
litológicos são estimadas implicitamente pelo operador de diferenças finitas. O método de
diferenças finitas aproxima a equação da onda pela diferencial de segunda ordem, estas são
resolvidas iterativamente em cima de uma malha de discretização espacial em esquema
de pequenas extrapolações.
A utilização do método de diferenças finitas para criar sismogramas sintéticos não é novo. Cada vez mais surgem trabalhos de sismologia usando esta técnica
e discutindo em detalhes aspectos particulares. O programa para modelagem terrestre
3
em malha normal, desenvolvido neste trabalho, utiliza aproximação de segunda ordem
para ambos os intervalos, espacial e temporal, enquanto que o desenvolvido para a modelagem marítima usa aproximação de quarta ordem no espaço e segunda ordem no tempo
em malha intercalada. Em ambos os programas foram utilizadas condições de absorção
de energia. O método de diferenças finitas nos permite obter fotografias (“snapshots”)
do campo de ondas em qualquer tempo da propagação, (“common shot point”) com sismogramas em formato de dados padrões e ele tem a capacidade de recorrer a múltiplos
sismogramas sintéticos em um ponto comum em profundidade (“common depth point”)
simulando a aquisição dos dados.
Foi utilizado para aplicação de filtros e visualização de snapshots e sismogramas o pacote de processamento sísmico “Seismic Unix”, desenvolvido pelo “Center for
Wave Phenomena - Colorado School of Mines”.
Este pacote é livremente distribuído, juntamente com todos os códigos
fontes.
Simogramas sintéticos de modelos complexos são úteis na exploração, servindo para ilustrar como a técnica pode ajudar o intérprete. Como exemplo podemos citar
a definição dos vários parâmetros na implementação da modelagem por diferenças finitas,
que envolvem dispersão da malha, reflexões artificiais vindo da borda do modelo, definição
dos incrementos das amostras espaciais e temporais do modelo dentre outros,Kelly et al
(1976) [15].
4
1.1
Elasticidade do Meio
As perturbações que se propagam em meio elástico denominam-se ondas
elásticas. Este meio é constituído de qualquer material que tenda a preservar seu comprimento, forma e volume contra as forças externas.
Os materiais que possuem forças restauradoras que tendem a retornar o
material à sua condição original após a remoção das forças externas são meios elásticos.
A força restauradora é característica do material e tem origem nas forças de ligação entre
seus átomos ou moléculas individuais.
Refere-se como meio não-dispersivo aquele em que a forma da onda não
se altera à medida que a onda se propaga e sua velocidade é constante desde que sejam
fixadas as características de elasticidade e a densidade do meio, pois a velocidade depende
em geral, destes parâmetros. A onda sonora no ar é um exemplo de onda que não sofre
dispersão.
1.2
Anisotropia do Meio
A anisotropia é um termo geral que se refere às mudanças nas propriedades
físicas resultantes das mudanças na direção ou polarização (Figura 1.1). Anisotropia na
velocidade sísmica é uma característica comum em muitos materiais da Terra, mas, em
geral, não é levado em consideração no processamento de dados sísmicos.
Virtualmente, todos os atuais algoritmos de processamento em uso assumem
que a Terra é isotrópica, porém, na verdade, a Terra é anisotrópica em quase todas as
suas características.
5
Direção
da linha
Direção
da linha
Figura 1.1: Direção de polarização e propagação definida para um caso anisotrópico de
meio com finas camadas paralelas. As polarizações são definidas em termos da linha
de direção indicada e são consideradas ambas as direções de propagação horizontal e
vertical. Na Figura 1.1 (a) é considerado um acamamento horizontal, enquanto na Figura
1.1 (b) o acamamento é vertical. As Figuras mostram o desmembramento das diferentes
velocidades com diferentes direções de propagação e polarização. Figura adaptada de
Tatham and McCormack (1991) [16].
6
a)
b)
Isotropia Transversal
Vertical
Isotropia Transversal
Horizontal
c)
d)
Ortorrombica
Monoclinica
Figura 1.2: Modelos litológicos equivalentes para diferentes tipos de simetrias.
Vale a pena ressaltar que as considerações de isotropia são feitas já muito
cedo no desenvolvimento das definições da equação que rege a propagação do campo
de onda. Se uma anisotropia média for considerada, o problema se torna muito mais
complexo. Para se manter o problema controlado, algumas considerações sobre a simetria
e geometria da anisotropia são comumente feitas.
Neste trabalho foi usada isotropia transversalmente vertical (TIV), como
mostra a Figura 1.2 (a), nestes experimentos os materiais anisotrópicos assumiram isotropia
transversa (sem variação azimutal na velocidade) onde as ondas S são descritas pela polarização como ondas SV e SH. As ondas SV e SH podem ter diferentes velocidades.
7
1.3
Heterogeneidade do Meio
Na formulação para meios heterogêneos as propriedades elásticas podem ser
especificadas para cada ponto da malha de diferenças finitas, e as condições de contorno
são satisfeitas implicitamente, Kelly et al (1976) [15].
O meio elástico pode ser considerado como um conjunto de regiões litológicas
homogêneas, cada uma caracterizada por valores constantes de densidade e parâmetros
elásticos.
A propagação do campo de onda em cada região pode ser apresentada por
uma representação de diferenças finitas apropriada para as equações elásticas correspondentes. Por este método, todas as condições de contorno através das interfaces de separação de diferentes regiões podem ser satisfeitas explicitamente. As aproximações são
incorporadas às condições de contorno implicitamente, construindo representação de diferenças finitas para equações elásticas gerais que são válidas para regiões heterogêneas.
Estas aproximações tornam possível associar diferentes densidades e valores
de parâmetros elásticos a todo ponto da malha, trazendo flexibilidade na variação de
geometrias de superfícies complexas.
A dissertação aqui apresentada tem o objetivo de introduzir dois métodos
para modelagens sísmicas em meios transversalmente isotrópicos, compatíveis com as especificidades encontradas nos modelos geológicos abordados nos processos de exploração e
explotação de hidrocarbonetos. Para alcançar este intento foi preciso resolver diretamente
a equação diferencial que rege o comportamento ondulatório por estratégias que não envolvessem soluções particulares, tais como, as utilizadas nos métodos da refletividade e
8
traçamento de Raios. Especificamente, apresentamos soluções baseadas no método de
diferenças finitas, no qual as derivadas parciais espaciais e temporais foram discretizadas
com aproximações de segunda e quarta ordens. Para o caso de abordagens terrestres,
generalizamos o algoritmo PS2 de Zahradník e Priolo (1994) [2] e para o caso marinho,
implementamos o método apresentado em Levander (1990) [4], com modificações para
que o mesmo pudesse tratar casos anisotrópicos. Em ambos os métodos, desde que sejam
satisfeitas as condições necessárias para garantir a não existência de dispersões numéricas
e estabilidade, não existe qualquer limitação com respeito à geometria das camadas que
compõem o modelo. Assim como não há qualquer restrição à variação dos parâmetros
elásticos no interior da mesma.
Esta dissertação apresenta três capítulos. O primeiro capítulo diz respeito
a esta introdução; o segundo apresenta o método que generaliza o de Zahradník e Priolo
(1994) [2], indicado para modelagens terrestres; e o terceiro apresenta a modificação do
método de Levander (1990) [4], visando modelagens marinhas. No final do segundo e
terceiro capítulos, são apresentados exemplos de aplicação dos métodos propostos. Após
o terceiro capítulo, apresentamos as conclusões, comentários e propostas para trabalhos
futuros nesta área. Concluindo: apresentamos dois apêndices nos quais exibimos o desenvolvimento matemático empregado para expressar o tensor de elasticidade em meios
transversalmente isotrópicos em termos do sistema de coordenadas que melhor explora
esta simetria (simetria principal), e explicitamos a discretização em malhas intercaladas.
9
Capítulo 2
Modelagem Sísmica para Meios
Elásticos TIV: Caso Terrestre
Neste capítulo apresentamos uma generalização do algoritmo PS2 introduzido por Zahradník e Priolo (1994) [2] visando modelagens sísmicas 2-D por diferenças
finitas em meios transversalmente isotrópicos, com eixo de simetria vertical. O algoritmo
é indicado para simulações de aquisições sísmicas em modelos terrestres e pode ser aplicado nos casos de superfícies de observação irregulares. Quanto à complexidade estrutural e
estratigráfica dos modelos que podem ser abordados por essa técnica, não existe qualquer
restrição com relação ao número de camadas, que podem ser separadas por interfaces
arbitrariamente irregulares, assim como não há qualquer limitação quanto à variação dos
parâmetros no interior das camadas, desde que sejam satisfeitas as condições que garantem a inexistência de dispersões numéricas e instabilidades. O método proposto é capaz
de abordar modelos geológicos realistas, sendo, portanto, uma valiosa ferramenta a ser
usada por geólogos e geofísicos envolvidos nos processos exploratórios e explotatórios de
10
hidrocarbonetos.
Na formulação do método apresentado neste capítulo, temos:
(1) as equações que regem o campo de deslocamento das partículas são
discretizadas, numa malha regular, com aproximações de segunda ordem nas derivadas
parciais;
(2) os parâmetros são introduzidos por integrações ao longo das linhas da
malha;
(3) o termo fonte é dado pela segunda derivada da função gaussiana, como
apresentado em Cunha (1997) [6];
(4) as condições de contorno não reflexivas são as propostas por Emerman
e Stephen (1990) [10] combinadas com as bordas de absorção sugeridas por Cerjan et al
(1985) [7] e
(5) aplica-se a condição de contorno conhecida por Formalismo do Vácuo,
apresentada em Zahradník, Moczo e Hron (1993) [3], na superfície de observação.
O programa Fortran desenvolvido tem como dados de entrada a densidade ρ,
os parâmetros anisotrópicos ε e δ de Thomsen (1986) [5], e as velocidades das ondas compressionais VP e cisalhantes VS na direção vertical, ao longo do modelo. Os componentes
do tensor elástico Cijkl são computados a partir destes parâmetros e são integrados ao
longo das linhas da malha. Os dados de saída são instantâneos dos componentes do campo de deslocamento e sismogramas, especificamente os componentes vertical e horizontal
do campo de deslocamentos em todos os tempos registrados na superfície de observação.
A generalização para modelagens 3-D pode ser realizada explorando-se as facilidades de
paralelização do código.
11
Com relação às aplicações:
• Com o intuito de estudar o efeito da anisotropia na propagação do campo de ondas,
foram realizadas modelagens numéricas para simulação da propagação de onda em 23
meios dos presentes no artigo de Thomsen (1986). Como resultado destas simulações
serão apresentados instantâneos da componente vertical da velocidade para todos
os modelos abordados.
• Com o objetivo de estimar os efeitos decorrentes da anisotropia nas reflexões sísmicas
foram realizadas simulações em um modelo composto por duas camadas homogêneas
sendo a primeira anisotrópica, separadas por uma interface curva (“senoidal”).
• Para mostrar a capacidade do método proposto em modelagens em meios complexos,
foram obtidos instantâneos e sismogramas para os componentes vertical e horizontal
do campo de ondas num modelo com alta complexidade geológica.
2.1
2.1.1
Desenvolvimento
Discretização das Equações da Elastodinâmica para meios
TIV
Considere o sistema de equações que rege o comportamento ondulatório em
meios elásticos heterogêneos transversalmente isotrópicos (com eixo de simetria vertical),
12
bidimensional, isto é:
µ
¶
µ
¶
∂u
∂v
∂
∂
∂ 2u
C11
+
C13
+
ρ 2 =
∂t
∂x
∂x
∂x
∂z
µ
¶
µ
¶
∂u
∂v
∂
∂
C44
+
C44
,
+
∂z
∂z
∂z
∂x
µ
¶
µ
¶
∂2v
∂v
∂u
∂
∂
ρ 2 =
C44
+
C11
+
∂t
∂x
∂x
∂z
∂x
µ
¶
µ
¶
∂v
∂u
∂
∂
C33
+
C44
.
+
∂z
∂z
∂x
∂z
(2.1)
(2.2)
Onde u e v denotam, respectivamente, os componentes horizontal e vertical
do campo de deslocamento e C11 , C13 , C44 e C33 são os componentes do tensor de elasticidade, definidos nas Equações A.19, A.20, A.18 e A.17 apresentadas no apêndice, e ρ é a
densidade. Temos que: n−1, n e n+1 representam os índices referentes, respectivamente,
ao passo de tempo anterior, atual e posterior de valores de u e v de um dado ponto da
malha, e h representa a variação espacial.
Introduzindo os índices citados anteriormente i, j e n, teremos:
u(x, y, t) = uni,j = ui,j ,
(2.3)
n
v(x, y, t) = vi,j
= vi,j ,
(2.4)
e
13
com
n = 0, 1, 2, 3, ... .
(2.5)
Denotando o lado direito da Equação 2.1 por $1 (x, z), ou seja escrevendo
ρ
∂ 2u
= $1 (x, z),
∂t2
(2.6)
e aproximando em 2a ordem a derivada parcial temporal temos:
n−1
n
un+1
1 n
i,j − 2ui,j + ui,j
∼
$ (i, j) ,
=
2
DT
ρi,j 1
(2.7)
ou ainda,
un+1
i,j =
DT 2 n
$ (i, j) + 2uni,j − un−1
i,j ,
ρi,j 1
(2.8)
onde,
i = 1, ..., Nx;
j = 1, ..., Nz;
n = 1, ..., Nt;
e
DT é o intervalo de tempo usado na discretização a ser definido na Equação
2.76.
14
Analogamente, podemos escrever:
n+1
=
vi,j
DT 2
n−1
n
$2 (i, j) + 2vi,j
,
− vi,j
ρi,j
(2.9)
onde $2 (i, j) representa as derivadas espaciais da Equação 2.2.
Na discretização dos termos $1 (i, j) e $2 (i, j), consideramos dois tipos de
derivadas parciais espaciais, denominadas simples e mista. As derivadas simples envolvem
uma única variável espacial, isto é, são do tipo,
µ
¶
∂
∂f
a
,
∂x
∂x
∂
∂z
µ
¶
∂f
a
.
∂z
As derivadas mistas são do tipo,
∂
∂z
µ
¶
∂f
a
,
∂x
µ
¶
∂
∂f
a
.
∂x
∂z
Onde a e f são funções de x e z, no caso específico um dos parâmetros elásticos e um dos
campos de deslocamento, respectivamente.
Assim sendo serão discretizados esses dois tipos de derivadas. Inicialmente
considere a do tipo simples,
µ
¶
∂
∂f
a
.
∂x
∂x
15
(2.10)
i,j+1
N
i-1,j
E
W
i+1,j
S
i,j-1
Figura 2.1: Linhas da malha onde são calculadas as médias geométricas dos parâmetros.
A discretização de
∂
∂z
µ
¶
∂f
a
,
∂z
(2.11)
é, naturalmente, similar.
A derivada simples pode ser reescrita como
∂
g,
∂x
onde
g=a
∂f
.
∂x
16
(2.12)
h
h
(i-1/2,j)
(i+1/2,j)
(i,j)
Figura 2.2: Modelo da malha do método PS2 de diferenças finitas para aproximação das
derivadas em relação a uma única variável para pontos da malha na superfície.
Assim podemos aproximá-la por,
∂ ∼ g(i + 12 , j) − g(i − 12 , j)
g=
,
∂x
h
(2.13)
onde h é a distância horizontal e vertical entre dois pontos consecutivos na malha regular,
Figura 2.2, sendo:
1
∂f
g(i ± , j) ∼
| 1 .
=a
2
∂x (i± 2 ,j)
(2.14)
Os valores de g, Equação 2.12, nos pontos (i ± 12 , j) serão estimados a partir das seguintes
aproximações:
R (i+1,j) ∂f
dx f (i + 1, j) − f (i, j)
(i,j)
1
∂x
∼
g(i + , j) = R
=
,
R (i+1,j) 1
(i+1,j) 1
2
dx
dx
(i,j)
(i,j)
a
a
17
(2.15)
(i,j-1)
N
E
W
(i-1,j)
(i,j)
(i+1,j)
S
(i,j+1)
Figura 2.3: Modelo da malha do método PS2 de diferenças finitas para aproximação das
derivadas em relação a uma única variável para pontos da malha no interior do modelo.
ou
1
aE
g(i + , j) ∼
(fi+1,j − fi,j ) ,
=
2
h
(2.16)
onde aE é definido por
h
h
=R
,
aE = R
1
(i+1,j) 1
dx
dx
(i,j)
E a
a
(2.17)
na qual o símbolo E na integral denota integração ao longo do segmento horizontal que
vai de (i, j) a (i + 1, j), Figura 2.3. Analogamente, o valor de g em (i − 12 , j) é estimado
18
por
1
aW
(fi,j − fi−1,j ) ,
g(i − , j) ∼
=
2
h
(2.18)
no qual
aW = R
(i,j)
h
(i−1,j)
e W indica que a integração de
1
dx
a
h
,
1
dx
W
a
=R
(2.19)
1
no segmento horizontal que vai de (i − 1, j) a (i, j),
a
Figura 2.3.
Assim sendo, a aproximação usada no cálculo das derivadas parciais simples
em relação a variável x será:
µ
¶
∂
∂f ∼ 1
a
= 2 {aE (fi+1,j − fi,j ) − aW (fi,j − fi−1,j )}.
∂x
∂x
h
(2.20)
De forma análoga, pode-se encontrar a fórmula que aproximará as derivadas simples em
relação a variável z, isto é,
∂
∂z
µ
¶
∂f ∼ 1
a
= 2 {aS (fi,j+1 − fi,j ) − aN (fi,j − fi,j−1 )},
∂z
h
(2.21)
onde aN e aS são dados por:
h
h
=R
,
aS = R
1
(i,j+1) 1
dx
dx
(i,j)
S a
a
19
(2.22)
e
aN = R
(i,j)
h
(i,j−1)
h
=R
,
1
1
dx
dx
N
a
a
(2.23)
e S e N denotam integrações ao longo dos segmentos de reta verticais de (i, j) a (i, j + 1)
e (i, j − 1) a (i, j), respectivamente, Figura 2.3.
Com relação às derivadas mistas,
∂
∂z
µ
¶
∂f
a
,
∂x
µ
¶
∂
∂f
a
,
∂x
∂z
usaremos duas aproximações. A primeira, que será referida por forma curta, será usada
no cálculo do campo de deslocamento no interior do modelo. A segunda, desenvolvida
por forma cheia, será usada na superfície (j=1), em conjunto com o chamado Formalismo
do Vácuo, como sugerido em Zahradník e Priolo (1994).
Inicialmente, será deduzido a forma curta da aproximação usada para,
∂
∂z
µ
¶
∂f
a
.
∂x
Assim, mais uma vez usando
∂f
g∼
=a
∂x
20
e
∂g ∼ g(i, j + 1/2) − g(i, j − 1/2)
,
=
∂z
h
(2.24)
1
∂f
1
g(i, j ± ) ∼
|
=a
2
∂x (i,j± 2 ),
(2.25)
onde
R (i,j+1) ∂f
dz
(i,j)
∂x
g(i, j + 1/2) ∼
= R
(i,j+1) 1
dz
(i,j)
a
(2.26)
e
R (i,j)
∂f
dz
∂x ,
g(i, j − 1/2) ' R
1
(i,j)
dz
(i,j−1) a
(i,j−1)
(2.27)
ou ainda, usando as definições de aS e aN ,
aS
g(i, j + 1/2) '
h
Z
aN
g(i, j − 1/2) '
h
Z
(i,j+1)
(i,j)
∂f
dz
∂x
(2.28)
∂f
dz,
∂x
(2.29)
e
21
(i,j+1)
(i,j)
ou
g(i, j + 1/2) '
aS ∂f
|(i,j+1/2) h
h ∂x
(2.30)
g(i, j − 1/2) '
aN ∂f
|(i,j−1/2) h.
h ∂x
(2.31)
g(i, j + 1/2) ' aS
fi+1,j+1/2 − fi−1,j+1/2
2h
(2.32)
g(i, j − 1/2) ' aN
fi+1,j−1/2 − fi−1,j−1/2
.
2h
(2.33)
fi±1,j+1/2 '
fi±1,j + fi±1,j+1
2
(2.34)
fi±1,j−1/2 '
fi±1,j + fi±1,j−1
,
2
(2.35)
e
Assim,
e
Então, usando as aproximações
e
22
chega-se
g(i, j + 1/2) '
aS
[fi+1,j+1 + fi+1,j − fi−1,j+1 − fi−1,j ] ,
4h
(2.36)
g(i, j − 1/2) '
aN
[fi+1,j + fi+1,j−1 − fi−1,j − fi−1,j−1 ] ,
4h
(2.37)
portanto
∂
∂z
µ
¶
∂f
1
a
' 2 {aS [fi+1,j+1 + fi+1,j − fi−1,j+1 − fi−1,j ]
∂x
4h
− aN [fi+1,j + fi+1,j−1 − fi−1,j − fi−1,j−1 ]} .
Analogamente, a forma curta de aproximações de
∂
∂x
(2.38)
¡ ∂f ¢
a ∂z é dada por:
µ
¶
∂u
1
∂
a
' 2 {aE [fi+1,j+1 + fi,j+1 − fi+1,j−1 − fi,j−1 ]
∂x
∂z
4h
− aW [fi,j+1 + fi−1,j+1 − fi,j−1 − fi−1,j−1 ]} .
A forma cheia de aproximação da derivada
∂
∂z
(2.39)
¡ ∂f ¢
a ∂x é obtida de maneira similar, contudo
a aproximação inicial é um pouco diferente, ou seja,
∂g
1
'
[g(i − 1/2, j + 1/2) + g(i + 1/2, j + 1/2)]
∂z
2h
−[g(i − 1/2, j − 1/2) + g(i + 1/2, j − 1/2)] ,
onde
23
(2.40)
g(i + 1/2, j + 1/2) = a
∂f
|(i+1/2,j+1/2)
∂x
(2.41)
que é aproximada por
R (i+1,j+1/2) ∂f
Z
dx
(i,j+1/2)
aSE (i+1,j+1/2) ∂f
∂x
dx
=
g(i + 1/2, j + 1/2) ' R
(i+1,j+1/2) 1
h (i,j+1/2)
∂x
dx
(i,j+1/2)
a
(2.42)
onde aSE é definido por
h
h
=R
aSE ' R
1
(i+1,j+1/2) 1
dx
dx
(i,j+1/2)
SE
a
a
e SE indica que a integração de
(2.43)
1
é realizada no segmento de reta que vai de (i, j + 1/2)
a
a (i + 1, j + 1/2), como mostrado na Figura 2.4.
Assim sendo,
aSE
[fi+1,j+1/2 − fi,j+1/2 ] ,
h
(2.44)
aSE
[fi+1,j+1 + fi+1,j − fi,j+1 − fi,j ] .
2h
(2.45)
g(i + 1/2, j + 1/2) '
ou ainda
g(i + 1/2, j + 1/2) '
24
(i,j)
SE
(i+1/2,j+1/2)
Figura 2.4: Modelo da malha do método PS2 de diferenças finitas para aproximação das
derivadas em relação às variáveis x e z para pontos da malha na superfície SE.
De forma análoga, encontramos:
g(i − 1/2, j + 1/2) '
aSW
[fi,j+1 + fi,j − fi−1,j+1 − fi−1,j ] ,
2h
(2.46)
onde
h
aSW ' R
(i,j+1/2)
(i−1,j+1/2)
1
dx
a
h
1
dx
SW a
=R
(2.47)
em que SW indica que a integração é efetuada no segmento de reta que vai de (i−1, j+1/2)
a (i, j + 1/2) como mostrado na Figura 2.5.
g(i + 1/2, j − 1/2) '
aEN
[fi+1,j + fi+1,j−1 − fi,j − fi,j−1 ] ,
2h
25
(2.48)
(i,j)
SW
(i-1/2,j+1/2)
Figura 2.5: Modelo da malha do método PS2 de diferenças finitas para aproximação das
derivadas em relação às variáveis x e z para pontos da malha na superfície SW.
e
g(i − 1/2, j − 1/2) '
aW N
[fi,j + fi,j−1 − fi−1,j − fi−1,j−1 ] ,
2h
(2.49)
onde
h
h
aNE ' R
=R
,
1
1
(i+1,j−1/2)
dx
dx
(i,j−1/2)
NE a
a
(2.50)
e
h
aNW ' R
(i,j−1/2)
(i−1,j−1/2)
1
dx
a
h
,
1
dx
NW a
=R
(2.51)
onde NE e NW indicam que as integrações são realizadas nos segmentos apresentados
na Figura 2.6.
26
EN
WN
(i-1/2,j-1/2)
(i,j)
Figura 2.6: Modelo da malha do método PS2 de diferenças finitas para aproximação das
derivadas em relação às variáveis x e z para pontos da malha na superfície NE e NW.
Com isso, a forma cheia de
∂
∂z
∂
∂z
¡ ∂f ¢
a ∂x escreve-se:
µ
¶
1
∂f
a
'
{aSE (fi+1,j+1 + fi+1,j − fi,j+1 − fi,j )
∂x
4h2
+aSW (fi,j+1 + fi,j − fi−1,j+1 − fi−1,j )
−aNE (fi+1,j + fi+1,j−1 − fi,j − fi,j−1 )
−aNW (fi,j + fi,j−1 − fi−1,j − fi−1,j−1 )}.
(2.52)
Finalmente, de forma análoga pode-se escrever:
µ
¶
∂
1
∂f
a
'
{aES (fi+1,j+1 + fi,j+1 − fi+1,j − fi,j )
∂x
∂z
4h2
+aEN (fi+1,j + fi,j − fi+1,j−1 − fi,j−1 )
−aW S (fi,j+1 + fi−1,j+1 − fi,j − fi−1,j )
−aW N (fi,j + fi−1,j − fi,j−1 − fi−1,j−1 )}
27
(2.53)
WN
EN
(i,j)
WS
ES
Figura 2.7: Modelo da malha do método PS2 de diferenças finitas para aproximação das
derivadas em relação às variáveis x e z para pontos da malha no interior WN, WS, EN e
ES.
ES, EN, W S e W N denotam os segmentos de reta usados nas integrações que definem
os parâmetros aES , aEN , aW S e aW N , respectivamente, Figura 2.7.
h
,
aES ' R
(i+1/2,j+1) 1
dz
(i+1/2,j)
a
h
aEN ' R
(i+1/2,j)
1
dz
a
h
,
aW S ' R
(i−1/2,j+1) 1
dz
(i−1/2,j)
a
h
aW N ' R
(i−1/2,j)
1
dz
a
(i+1/2,j−1)
(i−1/2,j−1)
(2.54)
(2.55)
Agora podemos escrever ás Equações 2.8 e 2.9 para o cálculo das componemtes do campo de onda de deslocamento em todos os pontos da malha, especificamente no interior (i = 2,...., Nx-1) e (j = 2,....,Nz-1) e na superfície (borda superior, ou
seja, i = 2, ....Nx-1 e j = 1).
Assim, tendo em conta que $n1 (i, j) e $n2 (i, j) representam a soma das
28
derivadas parciais espaciais nas Equações 2.8 e 2.9, respectivamente, ou seja:
$n1 (i, j)
µ
¶
µ
¶
∂u
∂v
∂
∂
C11
(i, j) +
C13
(i, j)
=
∂x
∂x
∂x
∂z
µ
¶
µ
¶
∂v
∂u
∂
∂
C44
(i, j) +
C44
(i, j)
+
∂z
∂x
∂z
∂z
(2.56)
e
$n2 (i, j)
µ
¶
µ
¶
∂v
∂u
∂
∂
C44
(i, j) +
C44
(i, j)
=
∂x
∂x
∂x
∂z
µ
¶
µ
¶
∂u
∂v
∂
∂
C11
(i, j) +
C33
(i, j) ,
+
∂z
∂x
∂z
∂z
(2.57)
onde usando as aproximações das derivadas encontradas acima (usando a forma curta
para as derivadas mistas), temos:
µ
¶n
∂
1
∂u
C11
(i, j) ' 2 {C11E (i, j)[ui+1,j − ui,j ]
∂x
∂x
h
−C11w (i, j)[ui,j − ui−1,j ]}
(2.58)
µ
¶n
∂
1
∂v
C13
(i, j) '
{C13E (i, j)[vi+1,j+1 − vi,j+1
∂x
∂z
4h2
−vi+1,j−1 − vi,j−1 ]
−C13w (i, j)[vi,j+1 + vi−1,j+1
−vi,j−1 − vi−1,j−1 ]
29
(2.59)
∂
∂z
µ
¶n
1
∂v
C44
(i, j) '
{C44S (i, j)[vi+1,j+1 + vi+1,j
∂x
4h2
−vi−1,j+1 − vi−1,j ]
−C44N (i, j)[vi+1,j + vi+1,j−1
−vi−1,j − vi−1,j−1 ]
∂
∂z
(2.60)
µ
¶n
1
∂u
(i, j) ' 2 {C44S (i, j)[ui,j+1 − ui,j ]
C44
∂z
h
−C44N (i, j)[ui,j − ui,j−1 ]}
(2.61)
µ
¶n
1
∂
∂v
(i, j) ' 2 {C44E (i, j)[vi+1,j − vi,j ]
C44
∂x
∂x
h
−C44w (i, j)[vi,j − vi−1,j ]}
(2.62)
µ
¶n
∂
1
∂u
C44
(i, j) '
{C44E (i, j)[ui+1,j+1 + ui,j+1
∂x
∂z
4h2
−ui+1,j−1 − ui,j−1 ]
−C44w (i, j)[ui,j+1 + ui−1,j+1
−ui,j−1 − ui−1,j−1 ]
30
(2.63)
∂
∂z
µ
¶n
1
∂u
C11
(i, j) '
{C11S (i, j)[ui+1,j+1 + ui+1,j
∂x
4h2
−ui−1,j+1 − ui−1,j ]
−C11N (i, j)[ui+1,j + ui+1,j−1
−ui−1,j − ui−1,j ]
∂
∂z
(2.64)
µ
¶n
1
∂v
(i, j) ' 2 {C33S (i, j)[vi,j+1 − vi,j ]
C33
∂z
h
−C33N (i, j)[vi,j − vi,j−1 ]} .
(2.65)
Desta forma, as Equações 2.1 e 2.2 podem ser escritas como:
un+1
=
i,j
DT 2
{C11E (i, j)(ui+1,j − ui,j ) − C11E (i − 1, j)(ui,j − ui−1,j )
ρ(i, j)h2
+0, 25[C13E (i, j)(vi+1,j+1 − vi,j+1 − vi+1,j−1 − vi,j−1 )
−C13E (i − 1, j)(vi,j+1 + vi−1,j+1 − vi,j−1 − vi−1,j−1 )]
+0, 25[C44S (i, j)(vi+1,j+1 + vi+1,j − vi−1,j+1 − vi−1,j )
−C44S (i, j − 1)(vi+1,j + vi+1,j−1 − vi−1,j − vi−1,j−1 )]
+C44S (i, j)(ui,j+1 − ui,j ) − C44S (i, j − 1)(ui,j− ui,j−1 )}
+2ui,j − un−1
i,j
(2.66)
31
e
n+1
=
vi,j
DT 2
{C44E (i, j)(vi+1,j − vi,j ) − C44E (i − 1, j)(vi,j − vi−1,j )
ρ(i, j)h2
+0, 25[C44E (i, j)(ui+1,j+1 + ui,j+1 − ui+1,j−1 − ui,j−1 )
−C44E (i − 1, j)(ui,j+1 + ui−1,j+1 − ui,j−1 − ui−1,j−1 )]
+0, 25[C11S (i, j)(ui+1,j+1 + ui+1,j − ui−1,j+1 − ui−1,j )
−C11S (i, j − 1)(ui+1,j + ui+1,j−1 − ui−1,j − ui−1,j−1 )]
+C33S (i, j)(vi,j+1 − vi,j ) − C33S (i, j − 1)(vi,j− vi,j−1 )}
n−1
+2vi,j − vi,j
,
(2.67)
para i = 2,...., Nx-1; j = 2,....Nz-1 e k = 1,...., Nt.
Usando o fato que CIJW (i, j) = CIJE (i − 1, j) e CIJN (i, j) = CIJE (i − 1, j),
reduzimos pela metade o número de parâmetros envolvidos.
Para o cálculo dos componentes uni,j , ν ni,j em j = 1, utilizamos os operadores cheios nas aproximações das derivadas mistas, como sugerido em Zahradník e Priolo
(1994), além das hipóteses relativas ao chamado Formalismo do Vácuo (Zahradník, Moczo
e Hron (1993)), ou seja, para j = 0 assumimos:
n
uni,j = vi,j
= CIJ (i, j) = ρ(i, j) = 0 .
32
(2.68)
Assim, temos:
=
un+1
i,1
DT 2
{C11E (i, 1)(ui+1,1 − ui,1 ) − C11E (i − 1, 1)(ui,1 − ui−1,1 )
ρ(i, 1)h2
+0, 25[C13ES (i, 1)(vi+1,2 − vi,2 − vi+1,1 − vi,1 )
−C13ES (i − 1, 1)(vi,2 + vi−1,2 − vi,1 − vi−1,1 )]
+C44S (i, j)(ui,2 − ui,1 )
+0, 25[C44SE (i − 1, j)(vi,2 + vi,1 − vi−1,2 − vi−1,1 )
−C44SE (i, 1)(vi+1,2 + vi+1,1 − vi,2 − vi,1 )]
+2ui,1 − un−1
i,1
(2.69)
e
n+1
vi,1
=
DT 2
{0, 25[C44ES (i, j)(ui+1,2 + ui,2 − ui+1,1 − ui,1 )
ρ(i, 1)h2
−C44ES (i − 1, 1)(ui,2 + ui−1,2 − ui,1 − ui−1,1 )]
−C44E (i, j)(vi+1,1 − vi,1 ) − C44E (i − 1, 1)(vi,1 − vi−1,1 )
+0, 25[C13SE (i − 1, 1)(ui,2 + ui,1 − ui−1,2 − ui−1,1 )
+C13SE (i, 1)(ui+1,2 + ui+1,1 − ui,2 − ui,1 )]
+C33S (i, j)(vi,2 − vi,j )}
+2vi,j − vi,j ,
(2.70)
para i = 2,....,Nx-1, e k=1,....,Nt.
33
2.1.2
Condições de Contorno para as Bordas Laterais e Inferior
Serão aqui apresentados, os operadores envolvidos na estratégia usada para
atenuação das reflexões nas bordas do modelo nas modelagens elásticas, ou seja, as
equações introduzidas por Cerjan (1985), Emerman e Stephen (1983) e as equações não
reflexivas (equações da onda one-way).
As condições de bordas de absorção propostas por Cerjan (1985) [7] com os
operadores propostos por Emerman e Stephen (1990) [10] são implementadas para solucionar os problemas que mais surgem na aplicação de métodos de soluções discretizadas
para cálculo de propagação de onda, a presença de reflexões e envoltórias vindas das bordas da malha numérica. Estes eventos indesejados prejudicam o sinal sísmico real que se
propaga na região modelada.
Uma das soluções para evitar o efeito das bordas é aumentar a malha numérica, assim as reflexões das bordas e as envoltórias têm o mesmo tempo envolvido na modelagem. Obviamente esta solução aumenta consideravelmente o custo computacional.
As condições de contorno para bordas não reflexivas, introduzidas para o
método de diferenças finitas (Clayton and Enquist (1977) [11]; Reynolds (1978) [12]) são
baseadas na troca da equação de onda numa faixa da região das bordas por uma equação
unidirecional que não permite que a energia de propagação que vem da reflexão da borda
volte para a malha numérica.
34
2.2
Termo Fonte
As fontes de energia sísmica mais utilizadas são a dinamite e o vibrador em
levantamentos terrestres e canhões de ar comprimido em levantamentos marítimos. Cada
uma destas fontes emite um pulso característico conhecido como assinatura da fonte que
se propaga em todas as direções com objetivo principal de encontrar reservatórios de óleo
e gás pelas propriedades reflexivas das rochas no interior da Terra. Estes pulsos elásticos
ou detonações são de duração ou comprimento muito pequeno, da ordem de 200 ms, e se
refletem e se refratam em cada uma das camadas geológicas em profundidade, retornando
a superfície, no caso de sísmica de superfície, ou para o poço, no caso de sísmica interpoços,
com informações valiosas para pesquisa de petróleo.
Para iniciar o processo de propagação de ondas sísmicas, foi adicionado às
equações que regem o campo de deslocamento das partículas, um termo fonte, especificamente em termo da densidade de força, como sugerido em Cunha (1997) [6]:
£
¤
2
f (t) = 1 − 2π (πfc t)2 e−π(πfc t) ,
(2.71)
realizando a discretização, substituindo t → (n − 1)∆t − tf , fica:
2
¤
£
f ((n − 1)∆t) = 1 − 2π (πfc ((n − 1)∆t − tf ))2 e−π(πfc ((n−1)∆t−tf )) ,
(2.72)
onde fc, frequência central está em função da frequência de corte fcorte :
fcorte
fc = √ .
3 π
35
(2.73)
Função Fonte (Fc=60 Hz)
0,6
0,4
0,2
Amplitude
0,0
-0,2
-0,4
-0,6
-0,8
-1,0
-1,2
0,00
0,02
0,04
0,06
0,08
0,10
0,12
Tempo (s)
Figura 2.8: Gráfico da função fonte: Amplitude x Tempo
Para frequência de corte de 60 Hz e 0,00015 s de intervalo de tempo temos
o gráfico da funcão fonte apresentado na Figura 2.8.
O comprimento do pulso (Nf ) e o período do pulso (tf ) no domínio do
tempo são:
√
4 π
Nf =
,
Dtfmax
√
2 π
tf =
.
fcorte
(2.74)
A adequada simulação da fonte requer um tratamento especial para os
casos homogêneo e heterogêneo. Os receptores utilizados para registrar as reflexões destes
pulsos são basicamente de dois tipos: eletromagnético (geofones) para registro em terra,
e de pressão (hidrofones) para levantamentos na água.
Os hidrofones utilizam cristais piezoelétricos, que geram uma corrente elétri36
ca proporcional à variação de pressão produzida pela onda. Estes receptores, a exemplo
dos geofones, devem reproduzir o mais fielmente possível as vibrações mecânicas na forma
de oscilações elétricas.
Estas oscilações elétricas são transmitidas até ao sismógrafo, onde são digitalizadas, multiplexadas e registradas (ou retransmitidas via satélite para uma central de
computadores), após severo depuramento e amplificação eletrônica [25].
2.3
2.3.1
Definição da malha
Condições para Não Ocorrência de Dispersão Numérica
As dimensões da malha são de importância vital para o método das dife-
renças finitas. As figuras mostrando os detalhes do modelo, como as dimensões totais em
ambas as direções, as velocidades em cada camada e as velocidades máximas e mínimas
para cada modelo, serão apresentadas na seção aplicações.
Considera-se que a função velocidade ci,j é discretizada dentro de um valor
médio para cada quadrado da malha. Esta hipótese é válida desde que os espaçamentos
da malha sejam pequenos comparados com o comprimento de onda da propagação.
Uma relação entre a menor velocidade utilizada no modelo (cmin ) e a freqüência máxima (fmax ), limita o máximo valor do espaçamento da malha de forma a
não se ter excessiva dispersão de energia [31], lembrando que neste modelo utiliza-se
37
h = ∆x = ∆z:
h=
cmin
,
ω.fmax
(2.75)
onde ω representa o número máximo de amostras por comprimento de onda correspondente à freqüência máxima. O valor ótimo encontrado de maneira empírica para este
número é 5.
2.3.2
Condições de Estabilidade Numérica
Um problema muito importante que deve ser considerado é a estabilidade
numérica. Para solucionar este problema foi desenvolvida uma relação para controle
dos valores dos intervalos do tempo de amostragem para evitar que o sistema se torne
numericamente instável, onde cmax é a maior velocidade adotada no modelo e µ é uma
constante definida da mesma forma que na dispersão da malha [31]. O melhor valor
encontrado para esta constante é 5.
Dt =
2.4
h
,
µ.cmax
(2.76)
Metodologia
A primeira metodologia é definida nos ítens abaixo:
1- Definição da malha:
• Maior valor de velocidade das ondas compressionais e menor valor de velocidade das
38
ondas cisalhantes, no modelo fornecido pelo geofísico intérprete.
• Frequência máxima da fonte sísmica (fmax ).
• Obtenção de h e Dt empregando as Equações 2.75 e 2.76.
2- Entrada de Dados:
Modelo sísmico:
Considerando i = 1,...,Nx e j = 1,...,Nz:
• Campo de velocidades das ondas compressionais: Vp(i,j);
• Campo de velocidades das ondas cisalhantes: Vs(i,j);
• Campo de densidade: ρ(i,j);
• Campos dos parâmetros de Thomsen: ε(i,j) e δ(i,j);
Obs: A malha deve ser pelo menos dez vezes mais densa do que a encontrada
com a condição de dispersão, para a realização das integrações ao longo das linhas da
malha, Equação 2.75.
Fonte sísmica:
• Frequência de corte;
• Tipo de fonte (explosiva ou direcional);
• Posição da fonte em relação à malha.
39
Geometria de Aquisição:
• Número e posição dos receptores em relação à malha.
Intervalo de amostragem temporal para saída:
• Tempo de registro;
• Definição da malha.
3- Integração dos parâmetros:
• Obtenção dos parâmetros elásticos C11 (i,j), C13 (i,j), C33 (i,j) e C44 (i,j);
• Obtenção dos parâmetros C11E (i,j), C11W (i,j), C11S (i,j), C11N (i,j), C13E (i,j), C13W (i,j),
C33S (i,j), C33N (i,j), C44S (i,j), C44W (i,j), C44N (i,j) e C44E (i,j) através de integrações
ao longo das linhas horizontais e verticais;
• Obtenção dos parâmetros C13ES (i,j), C44SE (i,j);
Estes parâmetros são definidos em relação à malha definida.
4- Execução do programa principal:
Como mostra a Figura 2.9 os dados de entradas, de integração dos parâmetros servem como entrada para o programa principal.
5- Saídas:
• Sismogramas relativos aos componentes vertical e horizontal do campo de deslocamento das partículas;
• Instantâneos e filmes apresentando a propagação dos campos de deslocamento.
40
—
Figura 2.9: Diagrama representando a 1a metodologia para modelagem em meios elásticos
TIV.
41
2.5
2.5.1
Aplicações
Primeira Aplicação: Propagação de Onda em Meios Homogêneos Anisotrópicos nos Modelos Propostos em Thomsen (1986)
O objetivo desta aplicação é exibir a assinatura sísmica das frentes de onda
qP e qSV em meios transversalmente isotrópicos com eixo de simetria vertical encontrados
na natureza. Especificamente, escolhemos 23 modelos dos apresentados em Thomsen
(1986), listados na Tablela 2.2. Os parâmetros usados nas modelagens estão descritos na
Tabela 2.1. As Figuras 2.10 a 2.21 exibem os resultados.
Tabela 2.1: Parâmetros da malha usada nas modelagens na primeira aplicação.
h
Dt(s)
fcorte(Hz)
Nx
Nz
Ntotal
ixf
jzf
Na
jobs
5
0,0005
60
601
601
2000
301
301
50
3
distância entre pontos consecutivos nas direções x e z (malha regular)
tempo de amostragem de diferenças finitas
frequência de corte
número de pontos da malha na direção x
número de pontos da malha na direção z
número total de passos de tempo
índice relativo à posição da fonte explosiva na direção x
índice relativo à posição da fonte explosiva na direção z
número de colunas ou linhas da borda não-reflexiva
índice relativo à posição da superfície de observação (receptores)
42
Tabela 2.2: Parâmetros de Thomsen usados na primeira aplicação (obtidos em Thomsen,
1986).
ε
δ
γ
Amostra
Vp (m/s) Vs (m/s)
Arenito
3368
1829
0,110 -0,127 0,255
Folhelho
4529
2703
0,034 0,250 0,046
Arenito Imaturo
4476
2814
0,097 0,051 0,051
Calcáreo Síltico
4972
2899
0,056 -0,041 0,067
Folhelho argiloso
3928
2055
0,334 0,818 0,575
Siltito Laminado
4449
2585
0,091 0,688 0,046
Arenito Calcarenoso
5460
3219
0,000 -0,345 -0,007
Folhelho com Óleo
4231
2539
0,200 0,000 0,145
Cinza Vulcânica
4846
1856
0,200 -0,003 0,105
Cristal
4420
2091
1,120 -1,230 2,280
Cristal de Quartzo
6096
4481
-0,096 0,169 -0,159
Calcita
5334
3353
0,369 0,127 0,169
Biotita
4054
1341
1,222 -1,437 6,120
Apatita
6340
4389
0,097 0,257 0,079
Gelo
3627
1676
-0,038 -0,100 0,031
Composto Alumínio-Lucita
2868
1350
0,970 -0,890 1,300
Arenito-Folhelho
3009
1654
0,013 -0,010 0,035
Folhelho Anisotrópico-SS
3009
1654
0,059 -0,042 0,163
Calcáreo-Folhelho
3306
1819
0,134 -0,094 0,156
Folhelho Anisotrópico-LS
3306
1819
0,169 -0,123 0,271
Folhelho Anisotrópico
2745
1508
0,103 -0,073 0,345
Gás-Arenito-Água-Arenito
1409
780
0,022 -0,002 0,004
Material de Gipso intemperizado
1911
795
1,161 -1,075 2,781
43
0
100
200
300
400
500
600
0
0
0
100
100
200
200
300
300
400
400
500
500
600
100
200
300
600
a)
400
500
600
b)
Figura 2.10: Instantâneos do componente vertical no tempo de 3 ms. a) alumínio-lucita
ε = 0, 970, δ = −0, 890. b) folhelho anisotrópico ε = 0, 103, δ = −0, 073.
0
100
200
300
400
500
600
0
0
0
100
100
200
200
300
300
400
400
500
500
600
600
a)
100
200
300
400
500
600
b)
Figura 2.11: Instantâneos do componente vertical no tempo de 3 ms. a) apatita ε = 0, 097,
δ = 0, 257. b) biotita ε = 1, 222, δ = −1, 437.
44
000
100
200
300
400
500
600
000
0
0
100
100
200
200
300
300
400
400
500
500
600
100
200
300
400
500
600
600
a)
b)
Figura 2.12: Instantâneos do componente vertical no tempo de 3 ms. a) calcita ε = 0, 369,
δ = 0, 127. b) arenito calcarenoso ε = 0, 000, δ = −0, 345.
000
100
200
300
400
500
600
000
0
0
100
100
200
200
300
300
400
400
500
500
600
600
a)
100
200
300
400
500
600
b)
Figura 2.13: Instantâneos do componente vertical no tempo de 3 ms. a) folhelho
anisotrópico-ss ε = 0, 059, δ = −0, 042. b) folhelho argiloso ε = 0, 334, δ = 0, 818.
45
000
100
200
300
400
500
600
000
0
0
100
100
200
200
300
300
400
400
500
500
600
600
100
200
a)
300
400
500
600
b)
Figura 2.14: Instantâneos do componente vertical no tempo de 3 ms. a) calcáreo síltico
ε = 0, 056, δ = −0, 041. b) arenito-folhelho ε = 0, 013, δ = −0, 010.
000
100
200
300
400
500
600
000
0
0
100
100
200
200
300
300
400
400
500
500
600
600
a)
100
200
300
400
500
600
b)
Figura 2.15: Instantâneos do componente vertical no tempo de 3 ms. a) arenito ε = 0, 110,
δ = −0, 127. b) cristal de quartzo ε = −0, 096, δ = 0, 169.
46
000
100
200
300
400
500
600
000
0
0
100
100
200
200
300
300
400
400
500
500
600
600
100
200
a)
300
400
500
600
b)
Figura 2.16: Instantâneos do componente vertical no tempo de 3 ms. a) folhelho com óleo
ε = 0, 200, δ = 0, 000. b) folhelho ε = 0, 034, δ = 0, 250.
000
100
200
300
400
500
600
000
0
0
100
100
200
200
300
300
400
400
500
500
600
600
a)
100
200
300
400
500
600
b)
Figura 2.17: Instantâneos do componente vertical no tempo de 3 ms. a) folhelho
anisotrópico-ls ε = 0, 169, δ = −0, 123. b) calcáreo-folhelho ε = 0, 134, δ = −0, 094
47
000
100
200
300
400
500
600
000
0
0
100
100
200
200
300
300
400
400
500
500
600
100
200
300
400
500
600
600
a)
b)
Figura 2.18: Instantâneos do componente vertical no tempo de 3 ms. a) siltito laminado
ε = 0, 091, δ = 0, 688. b) arenito imaturo ε = 0, 097, δ = 0, 051.
000
100
200
300
400
500
600
000
0
0
100
100
200
200
300
300
400
400
500
500
600
600
a)
100
200
300
400
500
600
b)
Figura 2.19: Instantâneos do componente vertical no tempo de 3 ms. a) gelo ε = −0, 038,
δ = −0, 100. b) material de gipso intemperizado ε = 1, 161, δ = −1, 075.
48
000
100
200
300
400
500
600
000
0
0
100
100
200
200
300
300
400
400
500
500
600
600
100
200
a)
300
400
500
600
b)
Figura 2.20: Instantâneos do componente vertical no tempo de 3 ms. a) gas-arenito-águaarenito ε = 0, 022, δ = −0, 002. b) cristal ε = 1, 120, δ = −1, 230.
000
100
200
300
400
500
600
0
100
200
300
400
500
600
Figura 2.21: Instantâneos do componente vertical no tempo de 3 ms. a) cinza vulcânica
ε = 0, 200, δ = −0, 003.
49
2.5.2
Segunda Aplicação: Modelo Formado por Duas Camadas
Separadas por Interface Curva
Nesta aplicação exibimos instantâneos e sismogramas encontrados em si-
mulações sísmicas realizadas nos modelos definidos na Figura 2.22, visando estimar os
efeitos nas reflexões da interface que separa os dois meios decorrentes da anisotropia. Os
parâmetros usados nas modelagens estão descritos na Tabela 2.3.
As Figuras 2.23 e 2.24 mostram instantâneos relativos a propagação do campo de onda em meio homogêneo com Vp =3000 m/s e VS =1732 m/s. Neles observamos um
aumento gradual da elipsidade no que diz respeito à propagação das ondas compressionais
(quase compressionais, rigorosamente), com o aumento dos valores de ε e δ. Observamos
também um resquício de onda S originado por não usarmos uma fonte exploxiva exata.
Ao contrário da frente de ondas compressionais, as ondas cisalhantes exibem um caráter
mais isotrópico, independente dos valores de ε e δ assumidos nesta aplicação.
Tabela 2.3: Parâmetros da malha usada nas modelagens na segunda aplicação.
h
Dt(s)
fcorte(Hz)
Nx
Nz
Ntotal
ixf
jzf
Na
jobs
5
0,0005
60
400
400
2000
200
3
50
2
distância entre pontos consecutivos nas direções x e z (malha regular)
tempo de amostragem de diferenças finitas
frequência de corte
número de pontos da malha na direção x
número de pontos da malha na direção z
número total de passos de tempo
índice relativo à posição da fonte explosiva na direção x
índice relativo à posição da fonte explosiva na direção z
número de colunas ou linhas da borda não-reflexiva
índice relativo à posição da superfície de observação (receptores)
As Figuras 2.27 e 2.28 são o resultado da subtração das ondas superficiais
50
Superfície (m)
0
0
100
200
300
400
500
600
700
800
900
1000
1100
1200
1300
1400
1500
1600
1700
1800
1900
200
400
600
Profundidade (m)
800
1000
1200
1400
1600
1800
Figura 2.22: Representa os modelos sísmicos abordados na segunda aplicação. A primeira
camada apresenta VP =3000 m/s e VS = 1732 m/s; a segunda possui VP =3800 m/s e
VS =2196 m/s. Nos quatro casos abordados a segunda camada é isotrópica. No primeiro
caso a camada superior é isotrópica e nos demais transversalmente isotrópica com ε = δ
iguais a 0,1, 0,2 e 0,3, respectivamente. Em todos os casos a densidade foi considerada
invariável ao longo do modelo (2,4 g/cm3 ).
51
Componente V
Superfície (m)
500
1000
1500
2000
0
0
500
500
Profundidade (m)
Profundidade (m)
0
1000
1500
Superfície (m)
0
500
0
500
1500
2000
1500
2000
1000
1500
2000
2000
Superfície (m)
0
500
1000
Superfície (m)
1500
2000
0
0
500
500
Profundidade (m)
Profundidade (m)
1000
1000
1000
1000
1500
1500
2000
2000
Figura 2.23: Instantâneos do componente vertical do campo de onda, considerando fonte
explosiva na posição central do modelo homogêneo, VP =3000 m/s e VS =1732 m/s. a)
Caso isotrópico. b) Caso anisotrópico, com ε = δ = 0, 1. c) Caso anisotrópico, com
ε = δ = 0, 2. d) Caso anisotrópico, com ε = δ = 0, 3. Observe o crescimento da elipsidade
com o aumento do valor dos parâmetros de Thomsen, no caso das ondas compressionais. A
frente de ondas cisalhantes (mais lenta e de menor amplitude) é decorrente da aproximação
usada na introdução da fonte explosiva. Para os valores de anisotropia escolhidos as ondas
cisalhantes apresentam um caráter mais isotrópico em relação às compressionais.
52
Componente U
Superfície (m)
500
1000
1500
2000
0
0
500
500
Profundidade (m)
Profundidade (m)
0
1000
1500
Superfície (m)
0
500
0
500
1500
2000
1500
2000
1000
1500
2000
2000
Superfície (m)
0
500
1000
Superfície (m)
1500
2000
0
0
500
500
Profundidade (m)
Profundidade (m)
1000
1000
1000
1000
1500
1500
2000
2000
Figura 2.24: Instantâneos do componente horizontal do campo de onda, considerando
fonte explosiva na posição central do modelo homogêneo, VP =3000 m/s e VS =1732 m/s.
a) Caso isotrópico. b) Caso anisotrópico, com ε = δ = 0, 1. c) Caso anisotrópico, com
ε = δ = 0, 2. d) Caso anisotrópico, com ε = δ = 0, 3. Observe o crescimento da elipsidade
com o aumento do valor dos parâmetros de Thomsen, no caso das ondas compressionais. A
frente de ondas cisalhantes (mais lenta e de menor amplitude) é decorrente da aproximação
usada na introdução da fonte explosiva. Para os valores de anisotropia escolhidos as ondas
cisalhantes apresentam um caráter mais isotrópico em relação às compressionais.
53
Componente V
Superfície (m)
500
1000
1500
2000
0
0
0.1
0.1
0.2
0.2
0.3
0.3
Tempo de registro (s)
Tempo de registro (s)
0
0.4
0.5
Superfície (m)
0
500
0.7
0.7
0.8
0.8
0
500
1500
2000
0
0
0.1
0.1
0.2
0.2
0.3
0.3
Tempo de registro (s)
Tempo de registro (s)
1500
2000
Superfície (m)
Superfície (m)
1000
2000
0.5
0.6
500
1500
0.4
0.6
0
1000
0.4
0.5
1000
0.4
0.5
0.6
0.6
0.7
0.7
0.8
0.8
Figura 2.25: Sismogramas relativos ao componente vertical do campo de ondas. Em todos
os casos observa-se a predominância das ondas superficiais (ondas P e S diretas e ondas
Rayleigh) em relação as ondas refletidas na interface que delimita as duas camadas. a)
Caso isotrópico. b) Caso anisotrópico, com ε = δ = 0, 1. c) Caso anisotrópico, com
ε = δ = 0, 2. d) Caso anisotrópico, com ε = δ = 0, 3.
54
Componente U
Superfície (m)
500
1000
1500
2000
0
0
0.1
0.1
0.2
0.2
0.3
0.3
Tempo de registro (s)
Tempo de registro (s)
0
0.4
0.5
Superfície (m)
0
500
0.7
0.7
0.8
0.8
0
500
1500
2000
0
0
0.1
0.1
0.2
0.2
0.3
0.3
Tempo de registro (s)
Tempo de registro (s)
1500
2000
Superfície (m)
Superfície (m)
1000
2000
0.5
0.6
500
1500
0.4
0.6
0
1000
0.4
0.5
1000
0.4
0.5
0.6
0.6
0.7
0.7
0.8
0.8
Figura 2.26: Sismogramas relativos ao componente horizontal do campo de ondas. Em
todos os casos observa-se a predominância das ondas superficiais (ondas P e S diretas e
ondas Rayleigh) em relação as ondas refletidas na interface que delimita as duas camadas.
a) Caso isotrópico. b) Caso anisotrópico, com ε = δ = 0, 1. c) Caso anisotrópico, com
ε = δ = 0, 2. d) Caso anisotrópico, com ε = δ = 0, 3.
55
Componente V
Superfície (m)
500
1000
1500
2000
0
0
0.1
0.1
0.2
0.2
0.3
0.3
Tempo de registro (s)
Tempo de registro (s)
0
0.4
0.5
Superfície (m)
0
500
0.7
0.7
0.8
0.8
0
500
1500
2000
0
0
0.1
0.1
0.2
0.2
0.3
0.3
Tempo de registro (s)
Tempo de registro (s)
1500
2000
Superfície (m)
Superfície (m)
1000
2000
0.5
0.6
500
1500
0.4
0.6
0
1000
0.4
0.5
1000
0.4
0.5
0.6
0.6
0.7
0.7
0.8
0.8
Figura 2.27: Sismogramas relativos ao componente vertical do campo de ondas, sem a
presença das ondas superficiais. a) Caso isotrópico. b) Caso anisotrópico, com ε = δ =
0, 1. c) Caso anisotrópico, com ε = δ = 0, 2. d) Caso anisotrópico, com ε = δ = 0, 3.
Na parte (a) indicamos os quatro modos de onda refletida presente nesta modelagem.
Observamos variações consideráveis nos tempos de trânsito das ondas PP nos offsets
longos. As ondas SP também apresentam variações, porém inferiores às apresentadas nas
ondas PP. As ondas PS e SS praticamente não apresentam variações.
56
Componente U
Superfície (m)
500
1000
1500
2000
0
0
0.1
0.1
0.2
0.2
0.3
0.3
Tempo de registro (s)
Tempo de registro (s)
0
0.4
0.5
Superfície (m)
0
500
0.7
0.7
0.8
0.8
0
500
1500
2000
0
0
0.1
0.1
0.2
0.2
0.3
0.3
Tempo de registro (s)
Tempo de registro (s)
1500
2000
Superfície (m)
Superfície (m)
1000
2000
0.5
0.6
500
1500
0.4
0.6
0
1000
0.4
0.5
1000
0.4
0.5
0.6
0.6
0.7
0.7
0.8
0.8
Figura 2.28: Sismogramas relativos ao componente horizontal do campo de ondas, sem a
presença das ondas superficiais. a) Caso isotrópico. b) Caso anisotrópico, com ε = δ =
0, 1. c) Caso anisotrópico, com ε = δ = 0, 2. d) Caso anisotrópico, com ε = δ = 0, 3.
Na parte (a) indicamos os quatro modos de onda refletida presente nesta modelagem.
Observamos variações consideráveis nos tempos de trânsito das ondas PP nos offsets
longos. As ondas SP também apresentam variações, porém inferiores às apresentadas nas
ondas PP. As ondas PS e SS praticamente não apresentam variações.
57
dos sismogramas completos exibidos nas Figuras 2.25 e 2.26. Na parte (a) das figuras
indicamos os modos de onda relativos às reflexões. Observamos, como esperado, variações
consideráveis no tempo de trânsito das ondas PP nos offsets maiores. As ondas SP
também são modificadas com o valor crescente de ε e δ. As ondas PS e SS praticamente
não são alteradas com a variação dos parâmetros de Thomsen usados nesta aplicação.
2.5.3
Terceira Aplicação: Modelo com Altos Contrastes de Impedância
Nesta aplicação mostramos os instantâneos e sismogramas obtidos na simu-
lação sísmica realizada num modelo que apresenta alto grau de complexidade estrutural,
descrito nas Figuras 2.29 e 2.30. O objetivo é mostrar a capacidade do método de abordar modelos realistas, além de exibir os efeitos decorrentes da anisotropia. Os parâmetros
usados nas modelagens estão dispostos na Tabela 2.4 e os resultados são apresentados nas
Figuras 2.31 até 2.35.
Nas figuras relativas aos instantâneos, Figuras 2.31 e 2.32, observamos um
espalhamento lateral de energia superior na propagação na camada anisotrópica, também
constatado nos maiores offsets nas Figuras 2.33 e 2.34 as partes (a) e (b) destas figuras exibem os simogramas completos relativos aos componentes vertical e horizontal do campo
de ondas, respectivamente; enquanto que as partes (c) e (d) mostram os mesmos sismogramas subtraídos das superficiais (ondas P e S diretas e ondas Rayleigh). Nestas últimas
observamos o maior espalhamento lateral de energia no caso anisotrópico, como indicado
nos instantâneos. Tal fato tem o potencial de degradar o imageamento sísmico quando
58
Superfície (m)
0
1000
2000
3000
4000
5000
6000
7000
8000
9000
0
10000
Vp= 1700 m/s
250
500
Vp= 2500 m/s
Profundidade (m)
750
1000
Vp= 3200 m/s
1250
Vp= 3700 m/s
1500
Vp= 4200 m/s
1750
Vp= 4700 m/s
2000
Vp= 5200 m/s
2250
Vp= 5000 m/s
2500
Vp= 4000 m/s
2750
Vp= 5100 m/s
Vp= 6345 m/s
Figura 2.29: Modelo de velocidades da bacia sedimentar terrestre estudada com os valores
de velocidade da onda P em cada camada. Os√ valores de velocidade da onda S podem
ser obtidos dividindo nos valores de VP por 3. O valor de densidade é considerado
constante ao longo do modelo, igual a 2,4 g/cm3 .
Superfície (m)
0
1000
2000
3000
4000
5000
0
6000
7000
8000
9000
10000
250
500
Profundidade (m)
750
Vp=2500 m/s
Epsilon = delta = 0,2
1000
1250
1500
1750
2000
2250
2500
2750
Figura 2.30: Distribuição dos parâmetros de Thomsen no modelo para o caso terrestre.
Com exceção da camada que corresponde à velocidade das ondas compressionais igual a
2500 m/s, na qual ε = δ = 0, 2, temos valores nulos para estes parâmetros.
59
o caráter anisotrópico da propagação do campo de onda não for levado em consideração
na estratégia de processamento dos dados, como mencionado em Rosa Filho (2002) [29].
A Figura 2.35 mostra em detalhe a onda tipo Rayleigh encontrada com a aplicação do
Formalismo do Vácuo na superfície que delimita o modelo na parte superior. Observe a
dispersão características na parte inferior da figura. Com relação ao tempo de processamento gasto na simulação de um tiro no modelo em questão, utilizando um computador
pessoal com processador Pentium IV com velocidade de 1.8 GHz e 1.0 Gb de memória
RAM, foram necessários 146 minutos de processamento.
Finalmente, embora tenhamos considerado geometrias de aquisição com
fontes e receptores na superfície da Terra o método não está limitado a qualquer tipo
de levantamento terrestre.
Tabela 2.4: Parâmetros da malha usada nas modelagens na terceira aplicação.
h
Dt(s)
fcorte(Hz)
Nx
Nz
Ntotal
ixf
jzf
Na
jobs
5
0,00015
60
2000
560
20000
500
5
50
5
distância entre pontos consecutivos nas direções x e z (malha regular)
tempo de amostragem de diferenças finitas
frequência de corte
número de pontos da malha na direção x
número de pontos da malha na direção z
número total de passos de tempo
índice relativo à posição da fonte explosiva na direção x
índice relativo à posição da fonte explosiva na direção z
número de colunas ou linhas da borda não-reflexiva
índice relativo à posição da superfície de observação (receptores)
60
Componente V
Caso Isotrópico
Caso Anisotrópico
Superfície (m)
100
200
300
400
500
600
Superfície (m)
700
800
900
1000
0
0
250
250
500
500
750
Profundidade (m)
Profundidade (m)
0
1000
1250
1500
1750
2000
100
200
300
400
0
100
200
300
400
0
100
200
300
400
0
100
200
300
400
800
900
1000
700
800
900
1000
700
800
900
1000
700
800
900
1000
1500
1750
2000
2250
2500
2750
100
200
300
400
500
600
Superfície (m)
700
800
900
1000
0
0
250
250
500
500
750
Profundidade (m)
Profundidade (m)
700
1250
2500
0
1000
1250
1500
1750
2000
500
600
750
1000
1250
1500
1750
2000
2250
2250
2500
2500
2750
2750
Superfície (m)
0
100
200
300
400
500
600
Superfície (m)
700
800
900
1000
0
0
250
250
500
500
750
Profundidade (m)
Profundidade (m)
600
750
Superfície (m)
1000
1250
1500
1750
2000
500
600
750
1000
1250
1500
1750
2000
2250
2250
2500
2500
2750
2750
Superfície (m)
0
100
200
300
400
500
600
Superfície (m)
700
800
900
1000
0
0
250
250
500
500
750
Profundidade (m)
Profundidade (m)
500
1000
2250
2750
0
1000
1250
1500
1750
2000
500
600
750
1000
1250
1500
1750
2000
2250
2250
2500
2500
2750
2750
Figura 2.31: Instantâneos relativos ao componente vertical para os modelos isotrópico (à
esquerda) e anisotrópico nos tempos de propagação 0,5; 0,75; 1,00 e 1,25 s. Como resultado
da anisotropia presente na segunda camada, observamos um espalhamento lateral de
energia superior na propagação no modelo anisotrópico.
61
Componente U
Caso Isotrópico
Caso Anisotrópico
Superfície (m)
100
200
300
400
500
600
Superfície (m)
700
800
900
1000
0
0
250
250
500
500
750
Profundidade (m)
Profundidade (m)
0
1000
1250
1500
1750
2000
100
200
300
400
0
100
200
300
400
0
100
200
300
400
0
100
200
300
400
800
900
1000
700
800
900
1000
700
800
900
1000
700
800
900
1000
1500
1750
2000
2250
2500
2750
100
200
300
400
500
600
Superfície (m)
700
800
900
1000
0
0
250
250
500
500
750
Profundidade (m)
Profundidade (m)
700
1250
2500
0
1000
1250
1500
1750
2000
500
600
750
1000
1250
1500
1750
2000
2250
2250
2500
2500
2750
2750
Superfície (m)
0
100
200
300
400
500
600
Superfície (m)
700
800
900
1000
0
0
250
250
500
500
750
Profundidade (m)
Profundidade (m)
600
750
Superfície (m)
1000
1250
1500
1750
2000
500
600
750
1000
1250
1500
1750
2000
2250
2250
2500
2500
2750
2750
Superfície (m)
0
100
200
300
400
500
600
Superfície (m)
700
800
900
1000
0
0
250
250
500
500
750
Profundidade (m)
Profundidade (m)
500
1000
2250
2750
0
1000
1250
1500
1750
2000
500
600
750
1000
1250
1500
1750
2000
2250
2250
2500
2500
2750
2750
Figura 2.32: Instantâneos relativos ao componente horizontal para os modelos isotrópico
(à esquerda) e anisotrópico nos tempos de propagação 0,5; 0,75; 1,00 e 1,25 s. Como resultado da anisotropia presente na segunda camada, observamos um espalhamento lateral
de energia superior na propagação no modelo anisotrópico.
62
Componente V
Caso Anisotrópico
Caso Isotrópico
Superfície (m)
1000
2000
3000
4000
5000
6000
Superfície (m)
7000
8000
9000 10000
0
0
0.3
0.3
0.6
0.6
0.9
0.9
Tempo de registro (s)
Tempo de registro (s)
0
1.2
1.5
1.8
1000
2000
3000
4000
0
1000
2000
3000
4000
6000
7000
8000
9000 10000
7000
8000
9000 10000
1.5
1.8
2.1
2.4
2.4
2.7
2.7
3.0
Superfície (m)
0
1000
2000
3000
4000
5000
6000
Superfície (m)
7000
8000
9000 10000
0
0
0.3
0.3
0.6
0.6
0.9
0.9
Tempo de registro (s)
Tempo de registro (s)
5000
1.2
2.1
3.0
0
1.2
1.5
1.8
5000
6000
1.2
1.5
1.8
2.1
2.1
2.4
2.4
2.7
2.7
3.0
3.0
Figura 2.33: Sismogramas relativos ao componente vertical para os modelos isotrópico (à
esquerda) e anisotrópico (ε = δ = 0, 2, na camada que apresenta Vp = 2500 m/s). Os
sismogramas C e D não apresentam as ondas superficiais. Observe na figura D que as
condições de contorno nas bordas laterais não foram tão efetivas nos casos anisotrópicos.
63
Componente U
Caso Anisotrópico
Caso Isotrópico
Superfície (m)
1000
2000
3000
4000
5000
6000
Superfície (m)
7000
8000
9000 10000
0
0
0.3
0.3
0.6
0.6
0.9
0.9
Tempo de registro (s)
Tempo de registro (s)
0
1.2
1.5
1.8
1000
2000
3000
4000
0
1000
2000
3000
4000
6000
7000
8000
9000 10000
7000
8000
9000 10000
1.5
1.8
2.1
2.4
2.4
2.7
2.7
3.0
Superfície (m)
0
1000
2000
3000
4000
5000
6000
Superfície (m)
7000
8000
9000 10000
0
0
0.3
0.3
0.6
0.6
0.9
0.9
Tempo de registro (s)
Tempo de registro (s)
5000
1.2
2.1
3.0
0
1.2
1.5
1.8
5000
6000
1.2
1.5
1.8
2.1
2.1
2.4
2.4
2.7
2.7
3.0
3.0
Figura 2.34: Sismogramas relativos ao componente horizontal para os modelos isotrópico
(à esquerda) e anisotrópico (ε = δ = 0, 2, na camada que apresenta Vp = 2500 m/s).
Neste caso, a Figura D não apresenta reflexões laterais com a mesma intensidade das
observadas com relação ao componente vertical, (Figura 2.33 D).
64
Figura 2.35: Ondas superficiais, especificamente as ondas diretas P e S, e as ondas
Rayleigh. Os sismogramas c e d das Figuras 2.33 e 2.34 foram obtidos como resultado da
diferença entre os sismogramas completos (Figuras 2.33 e 2.34 a e b) e os que contém as
ondas superficiais exclusivamente. Observe na ampliação a dispersão das ondas Rayleigh.
65
Capítulo 3
Modelagem Sísmica para Meios
Elásticos TIV: Caso Marinho
Neste capítulo implementa-se uma generalização do método proposto por
Levander (1990) visando modelagem sísmica em meios elásticos anisotrópicos marinhos.
Neste método utiliza-se aproximações de 2a ordem para as derivadas temporais e de 4a
ordem para as derivadas que envolvem x e z. Os campos de velocidades (derivadas temporais dos campos de deslocamento) e as tensões τ xx , τ zz e τ xz , assim como os parâmetros do
modelo, são introduzidos através de malhas intercaladas, portanto, de acordo com o artigo
supracitado, a condição para inexistência da dispersão numérica apresenta dependência
com a velocidade das ondas compressionais exclusivamente.
Desta forma evita-se os problemas encontrados em modelagens em meios
acústicos usando o método apresentado no capítulo anterior.
Nas modelagens realizadas neste capítulo foram usadas a fonte, condição de
contorno de Cerjan (nas bordas laterais e inferior) iguais as do capítulo anterior e condição
66
Figura 3.1: Instantâneos relativos ao componente vertical, ilustrando os problemas de
dispersão numérica encontrados com o método apresentado no capítulo anterior. Num
caso considerando µ = 0 (meio acústico).
de superfície rígida.
A Figura 3.1 ilustra os problemas que surgem quando a constante de cisalhamento é nula em modelagens usando o método proposto no capítulo 2.
3.1
Desenvolvimento
Este algoritmo é uma modificação do apresentado no capítulo anterior, usan-
do o modelo de diferenças finitas em malha intercalada de Levander (1990) que considera
aproximações de derivada parcial de quarta ordem no espaço e de derivada parcial de
segunda ordem no tempo, para propagação da onda P-SV em meios heterogêneos 2D
baseado na formulação de malha intercalada de Madariaga-Vireux.
Este tipo de derivação é utilizada na modelagem com as equações de velocidades e tensões, nas quais uma dada quantidade pode ser definida nos pontos com posição
67
inteira, e a outra quantidade em posição fracionária. A vantagem evidente é a própria
proximidade dos pontos com maior peso em relação ao ponto de cálculo da derivada, o
que leva a uma maior precisão [26].
Implementando um sistema de equações de onda para meios TVI em termos das derivadas temporais, dos deslocamentos vertical e horizontal (Vz e Vx, respectivamente) e tensões (τ xx , τ xz e τ zz ) e avaliando Vx e a densidade nas posições inteiras
da malha (i,j), Vz e a densidade nas posições fracionárias da malha (i+1/2,j+1/2), os
esforços compressionais τ xx e τ zz e as constantes elásticas C11 , C13 e C33 são avaliados
em (i+1/2 ,j), o campo de esforços cisalhantes τ xz e a constante C44 são avaliados em
(i,j+1/2), conforme mostra as Figuras 3.2 e 3.3. Os parâmetros que são avaliados em
posições fracionárias são obtidos do modelo original por interpolações lineares simples
entre pontos vizinhos.
O modelo numérico foi desenvolvido apartir das Equações 3.1 e 3.2, que
compõem o sistema de equações elásticas hiperbólicas de primeira ordem, e as Equações
3.4, 3.5 e 3.6, leis constitutivas expressas em velocidades da partícula e tensões. As
duas componentes de velocidade não podem ser definidas no mesmo nó para a malha
intercalada completa. A condição de estabilidade e a curva de dispersão da velocidade
de fase da onda P não dependem da razão de Poisson, o comportamento da curva de
dispersão da velocidade de fase da onda S tem dependência fraca em relação a razão de
Poisson.
O método de malha intercalada tem a qualidade de poder corrigir modelos
com variações nas propriedades dos materiais, incluindo pequenas e grandes razões de
Poisson, com dispersão numérica mínima e anisotropia numérica. A análise de dispersão
68
indicada para um pequeno comprimento de onda no modelo precisa ser amostrado em 5
pontos da malha por comprimento de onda [20].
Este método pode ser usado para simular melhor a propagação de onda
em meios elásticos e acústicos alternados, tem bom comportamento em meios líquidos,
onde a velocidade da onda S vai para zero, e não é necessário tratamento especial para a
interface líquido-sólido, sendo ele ideal para simulação sísmica marítima, ou seja é estável
para contatos rocha-fluido [4].
3.1.1
Discretização pelo Método de Diferenças Finitas
As Equações 2.1 e 2.2 que descrevem o campo de deslocamento das partícu-
las em meios TIV que podem ser reescritos em termos dos campos de velocidades U e V
e tensões τ xx , τ xz e τ zz , temos:
∂ 2 u ∂τ xx ∂τ xz
=
+
∂t2
∂x
∂z
(3.1)
∂ 2w
∂τ zx ∂τ zz
+
,
=
∂t2
∂x
∂z
(3.2)
ρ
e
ρ
malha de velocidades
u(i, j) e w(i + 1/2, j + 1/2),
69
(3.3)
Figura 3.2: Diagrama de velocidade. Malha intercalada de diferenças finitas e diagrama espacial para as atualizações das velocidades. O quadrado cinza tem área h2 . Os
vértices do quadrado são pontos da malha (i,j), (i+1,j), (i+1,j+1) e (i,j+1). No diagrama
a velocidade horizontal é definida em (i,j), a velocidade vertical em índices fracionários
(i+1/2,j+1/2), a tensão normal em (i+1/2,j) e a tensão cisalhante em (i,j+1/2). As componentes de velocidades são definidas em intervalos de tempo (n+1/2)∆t e (n-1/2)∆t, entretanto as componentes de tensão são definidas em intervalos de tempo n∆t e (n+1)∆t.
Os diagramas espaciais para velocidade horizontal e vertical são mostrados, respectivamente, pela linha contínua e pontilhada, com os nós das tensões usados nas atualizações.
(Adaptado de Levander, 1988).
70
e as leis constitutivas
τ xx = C11
∂u
∂w
+ C13
,
∂x
∂z
(3.4)
τ zx = C44
µ
¶
(3.5)
∂w
∂u
+ C13
,
∂z
∂x
(3.6)
∂u ∂w
+
∂z
∂x
e
τ zz = C33
ou seja, se considerarmos os componentes U e V do campo de velocidade, temos o sistema
de equações parciais de 1a ordem,
ρ
∂τ xx ∂τ xz
∂U
=
+
,
∂t
∂x
∂z
(3.7)
ρ
∂V
∂τ zx ∂τ zz
=
+
,
∂t
∂x
∂z
(3.8)
∂τ xx
∂U
∂V
= C11
+ C13
,
∂t
∂x
∂z
(3.9)
µ
(3.10)
∂τ zx
= C44
∂t
∂U ∂V
+
∂z
∂x
71
¶
Figura 3.3: Diagrama de tensões. Malha intercalada de diferenças finitas e diagrama
espacial para as atualizações de tensões. O diagrama espacial para atualizações das tensões cisalhantes e normal, são mostrados, respectivamente, por uma linha pontilhada e
contínua, com os nós das tensões usados nas atualizações. (Adaptado de Levander, 1988).
e
∂τ zz
∂V
∂U
= C33
+ C13
.
∂t
∂z
∂x
(3.11)
Malha de tensões:
τ xx (i + 1/2, j),
τ zz (i + 1/2, j)
e
τ zx (i, j + 1/2).
Agora, considerando aproximações de 2a ordem nas derivadas parciais temporais e 4a ordem nas derivadas parciais com relação às variáveis espaciais, como em
72
Levander (1990), e tendo em vista as Figuras 3.2 e 3.3, podemos escrever:
ρ(i, j)
9
1 1 n
U n+1/2 (i, j) − U n−1/2 (i, j)
=
[ τ xx (i − 3/2, j) − τ nxx (i − 1/2, j)
Dt
Dx 24
8
1
1 1 n
9
[ τ (i, j − 3/2)
+ τ nxx (i + 1/2, j) − τ nxx (i + 3/2, j)] +
8
24
Dx 24 xz
1
9
9
(3.12)
− τ nxz (i, j − 1/2) + τ nxz (i, j + 1/2) − τ nxz (i, j + 3/2)],
8
8
24
V n+1/2 (i + 1/2, j + 1/2) − V n−1/2 (i + 1/2, +1/2j)
=
Dt
9
9
1 1 n
[ τ xz (i − 1, j + 1/2) − τ nxz (i, j + 1/2) + τ nxz (i + 1, j + 1/2)
Dx 24
8
8
9
1
1 1 n
[ τ zz (i + 1/2, j − 1) − τ nzz (i + 1/2, j)
− τ nxz (i + 2, j + 1/2)] +
24
Dz 24
8
1
9
+ τ nzz (i + 1/2, j + 1) − τ nzz (i + 1/2, j + 2)],
8
24
ρ(i + 1/2, j + 1/2)
(3.13)
n
τ n+1
1 1 n+1/2
xx (i + 1/2, j) − τ xx (i + 1/2, j)
(i − 1, j)
= C11 (i + 1/2, j)
[ U
Dt
Dx 24
1
9
9
− U n+1/2 (i, j) + U n+1/2 (i + 1, j) − U n+1/2 (i + 2, j)]
8
8
24
1 1 n+1/2
9
[ V
(i + 1/2, j − 3/2) − V n+1/2 (i + 1/2, j − 1/2)
+C13 (i + 1/2, j)
Dz 24
8
1
9
(3.14)
+ V n+1/2 (i + 1/2, j + 1/2) − V n+1/2 (i + 1/2, j + 3/2)],
8
24
73
n
τ n+1
1 1 n+1/2
zz (i + 1/2, j) − τ zz (i + 1/2, j)
(i − 1, j)
= C13 (i + 1/2, j)
[ U
Dt
Dx 24
1
9
9
− U n+1/2 (i, j) + U n+1/2 (i + 1, j) − U n+1/2 (i + 2, j)]
8
8
24
9
1 1 n+1/2
[ V
(i + 1/2, j − 3/2) − V n+1/2 (i + 1/2, j − 1/2)
+C33 (i + 1/2, j)
Dz 24
8
1
9
(3.15)
+ V n+1/2 (i + 1/2, j + 1/2) − V n+1/2 (i + 1/2, j + 3/2)]
8
24
e
n
τ n+1
1 1 n+1/2
xz (i, j + 1/2) − τ xz (i, j + 1/2)
= C44 (i, j + 1/2)
[ V
(i − 3/2, j + 1/2)
Dt
Dx 24
9
9
− V n+1/2 (i − 1/2, j + 1/2) + V n+1/2 (i + 1/2, j + 1/2)
8
8
1
1 1 n+1/2
[ U
(i, j − 1)
− V n+1/2 (i + 3/2, j + 1/2)] + C44 (i, j + 1/2)
24
Dz 24
9
1
9
(3.16)
− U n+1/2 (i, j) + U n+1/2 (i, j + 1) − U n+1/2 (i, j + 2)] .
8
8
24
Assim sendo, as fórmulas recursivas usadas na computação dos campos de velocidades e
tensões, são dadas por:
9
1
Dt
{[ τ nxx (i − 3/2, j) − τ nxx (i − 1/2, j)
hρ(i, j) 24
8
1
1
9
+ τ nxx (i + 1/2, j) − τ nxx (i + 3/2, j)] + [ τ nxz (i, j − 3/2)
8
24
24
1
9
9
− τ nxz (i, j − 1/2) + τ nxz (i, j + 1/2) − τ nxz (i, j + 3/2)]} + U n−1/2 (i, j),
8
8
24
U n+1/2 (i, j) =
74
(3.17)
Dt
1
{[ τ nxz (i − 1, j + 1/2)
hρ(i + 1/2, j + 1/2) 24
1
9
9
− τ nxz (i, j + 1/2) + τ nxz (i + 1, j + 1/2) − τ nxz (i + 2, j + 1/2)]
8
8
24
9
9
1
+[ τ nzz (i + 1/2, j − 1) − τ nzz (i + 1/2, j) + τ nzz (i + 1/2, j + 1)
24
8
8
1
− τ nzz (i + 1/2, j + 2)]} + V n−1/2 (i + 1/2, j + 1/2),
24
V n+1/2 (i + 1/2, j + 1/2) =
9
Dt
1
{C11 (i + 1/2, j)[ U n+1/2 (i − 1, j) − U n+1/2 (i, j)
h
24
8
1
9
+ U n+1/2 (i + 1, j) − U n+1/2 (i + 2, j)]
8
24
9
1
+C13 (i + 1/2, j)[ V n+1/2 (i + 1/2, j − 3/2) − V n+1/2 (i + 1/2, j − 1/2
24
8
1
9
+ V n+1/2 (i + 1/2, j + 1/2) − V n+1/2 (i + 1/2, j + 3/2)]} + τ nxx (i + 1/2, j),
8
24
(3.18)
τ n+1
xx (i + 1/2, j) =
9
τ n+1
1
Dt
zz (i + 1/2, j)
=
{C13 (i + 1/2, j)[ U n+1/2 (i − 1, j) − U n+1/2 (i, j)
Dt
h
24
8
1
9
+ U n+1/2 (i + 1, j) − U n+1/2 (i + 2, j)]
8
24
1
9
+C33 (i + 1/2, j)[ V n+1/2 (i + 1/2, j − 3/2) − V n+1/2 (i + 1/2, j − 1/2)
24
8
1
9
+ V n+1/2 (i + 1/2, j + 1/2) − V n+1/2 (i + 1/2, j + 3/2)]} + τ nzz (i + 1/2, j),
8
24
75
(3.19)
(3.20)
e
Dt
1
{C44 (i, j + 1/2)[ V n+1/2 (i − 3/2, j + 1/2)
h
24
9
9
− V n+1/2 (i − 1/2, j + 1/2) + V n+1/2 (i + 1/2, j + 1/2)
8
8
1
1
− V n+1/2 (i + 3/2, j + 1/2)] + C44 (i, j + 1/2)[ U n+1/2 (i, j − 1)
24
24
1
9
9
− U n+1/2 (i, j) + U n+1/2 (i, j + 1) − U n+1/2 (i, j + 2)]
8
8
24
τ n+1
xz (i, j + 1/2) =
+τ nxz (i, j + 1/2) .
3.2
(3.21)
Metodologia
A segunda metodologia é definida nos ítens abaixo:
1- Entrada de Dados:
Modelo sísmico:
• Campo de velocidades das ondas compressionais: Vp(i,j);
• Campo de velocidades das ondas cisalhantes: Vs(i,j);
• Campo de densidade: ρ(i,j);
• Campos dos parâmetros de Thomsen: ε(i,j) e δ(i,j);
Obs: A malha deve ser pelo menos duas vezes mais densa do que a encontrada com a condição de dispersão.
Fonte sísmica:
• Frequência de corte;
76
• Tipo de fonte (explosiva ou direcional);
• Posição da fonte em relação á malha;
Geometria de Aquisição:
• Número e posição dos receptores em relação á malha;
Intervalo de amostragem temporal para saída.
• Tempo de registro.
2- Compilação:
• Obtenção dos parâmetros elásticos C11 (i,j), C13 (i,j), C33 (i,j) e C44 (i,j);
• Obtenção dos parâmetros aE (i,j), aS (i,j) através de integrações ao longo das linhas
horizontais e verticais; Estes parâmetros são definidos em relação à malha definida.
3- Execução do programa principal:
Como mostra a Figura 3.4 os dados de entradas, de compilação e de condições
de bordas servem como imput para o programa principal.
4-Saídas:
• Sismogramas relativos aos componentes vertical e horizontal do campo de velocidades das partículas;
• Sismogramas relativos aos componentes do tensor de tensões, τ xx , τ zz e τ xz ;
• Instantâneos e filmes apresentando a propagação dos campos de velocidades e tensões.
77
Modelagem
Meios Elásticos
TIV
Método 2
•
•
•
•
•
Campos de Tensões
Campos de Velocidades
Vp(i,j), Vs(i,j)
Densidade ρ (i,j)
Epsilon ε (i, j ) , delta δ (i, j )
Cálculo dos
Parâmetros Elásticos
C11, C13, C33, C44
C11e, C13e, C33e, C44e
C11s, C13s, C33s, C44s
•
• Arranjo
Fonte/Receptor
Fonte Sísmica
Fcorte
•
•
Intervalo de
Amostragem Temporal
Tempo de registro
Programa
Principal
Saídas:
• Sismogramas
• Instantâneos
Figura 3.4: Diagrama representando a 2a metodologia para modelagem em meios elásticos
TIV.
78
3.3
3.3.1
Aplicações
Primeira Aplicação: Modelagem em Meios Acústicos com
operadores Elásticos
O objetivo desta aplicação foi mostrar a estabilidade do segundo método
em relação ao valor da constante de cisalhamento µ. Especificamente para modelos acústicos homogêneos, foram realizadas modelagens acústicas usando operadores recursivos,
Equações 3.17 e 3.21. Os parâmetros usados nas modelagens estão dispostos na Tabela
3.1 e os resultados são apresentados na Figura 3.5, partes (a), (b), (c), (d) e (e), as quais
correspondem aos instantâneos dos componentes U, V, τ xx , τ xz e τ zz , respectivamente.
Ao contrário do resultado com o método introduzido no capítulo anterior, não observamos
dispersões numéricas nestes resultados.
Tabela 3.1: Parâmetros da malha usada nas modelagens na primeira aplicação para o
caso marinho.
h
Dt(s)
fcorte(Hz)
Nx
Nz
Ntotal
ixf
jzf
Na
jobs
5
0,0002
60
400
400
1200
200
200
40
3
distância entre pontos consecutivos nas direções x e z
tempo de amostragem de diferenças finitas
frequência de corte
número de pontos do grid na direção x
número de pontos do grid na direção z
número total de passos de tempo
índice relativo à posição da fonte explosiva no direção x
índice relativo à posição da fonte explosiva no direção z
número de colunas ou linhas da borda não-reflexiva
índice relativo à posição da superfície de observação (receptores)
79
Superfície (m)
200
400
600
1000
1200
Superfície (m)
1400
1600
1800
2000
200
200
400
400
600
600
Profundidade (m)
0
800
1000
1200
1600
1800
1800
600
800
1000
1200
600
800
0
200
400
600
800
1400
1600
1800
2000
2000
Superfície (m)
400
400
1400
1600
1800
2000
0
200
200
400
400
600
600
1000
1200
Profundidade (m)
0
800
1200
1400
1600
1800
2000
1400
1600
1800
2000
Superfície (m)
1000
1200
800
1000
1200
1400
1400
1600
1600
1800
1800
2000
1000
1200
1400
200
200
800
1600
0
0
1000
1400
2000
Profundidade (m)
800
Profundidade (m)
Profundidade (m)
0
0
2000
Superfície (m)
0
200
400
600
800
1000
1200
0
200
400
Profundidade (m)
600
800
1000
1200
1400
1600
1800
2000
Figura 3.5: Instantâneos relativos aos componentes horizontal e vertical do campo de
velocidade (derivada do campo de deslocamento), U e V, aos componentes do tensor de
esforços τ xx , τ xz e τ zz , respectivamente, considerando um meio acústico (µ = 0). Observe
que não há dispersão numérica e o componente de tensão cisalhante τ xz é nulo, como era
de se esperar tratando-se de modelagem em meio acústico.
80
3.3.2
Segunda Aplicação: Modelagem Sísmica em Modelo Típico de Margens Passivas
Nesta aplicação mostramos a capacidade, do método introduzido neste capí-
tulo, para abordar modelos geológicos complexos, típicos de regimes distencionais. O modelo encontra-se descritos na Figura 3.6, os modelos de densidade e de distribuição dos
parâmetros de Thomsen estâo apresentados nas Figuras 3.7 e 3.8. Os parâmetros usados
na modelagem estão definidos na Tabela 3.2 e resultados são apresentados nas Figuras 3.9
à 3.18 . Constata-se a não existência de dispersão numérica e instabilidade na propagação
da onda ao longo do modelo. Os altos contrastes de impedância não são obstáculos para
este tipo de modelagem. Com relação ao tempo de processamento gasto na simulação
de um tiro no modelo em questão, utilizando uma estação de trabalho com processador
Intel Xeon com velocidade de 2.8 GHz e 1.0 Gb de memória RAM, foram necessários 126
minutos de processamento.
Tabela 3.2: Parâmetros da malha usada nas modelagens na segunda aplicação para o caso
marinho.
h
Dt(s)
fcorte(Hz)
Nx
Nz
Ntotal
ixf
jzf
Na
jobs
10
0,0004
30
3000
791
20000
1500
5
40
3
distância entre pontos consecutivos nas direções x e z
tempo de amostragem de diferenças finitas
frequência de corte
número de pontos do grid na direção x
número de pontos do grid na direção z
número total de passos de tempo
índice relativo à posição da fonte explosiva no direção x
índice relativo à posição da fonte explosiva no direção z
número de colunas ou linhas da borda não-reflexiva
índice relativo à posição da superfície de observação (receptores)
81
Superfície (m)
0
1000
2000
3000
4000
5000
6000
7000
8000
9000
10000
0
250
500
Profundidade (m)
750
1000
1250
1500
1750
2000
2250
2500
2750
Figura 3.6: Modelo de velocidades da bacia sedimentar marítima estudada com os valores
de velocidade da onda P em cada camada.
Superfície (m)
0
1000
2000
3000
4000
5000
6000
7000
8000
9000
10000
0
250
Densidade = 1 g/cm
500
3
Profundidade (m)
750
1000
1250
1500
1750
2000
Densidade = 2,2 g/cm
3
2250
2500
2750
Figura 3.7: Modelo de densidade no modelo para o caso marinho, Figura 3.6.
82
Superfície (m)
0
1000
2000
3000
4000
5000
6000
7000
8000
9000
0
10000
250
Epsilon = delta = 0
500
Profundidade (m)
750
1000
1250
1500
1750
2000
2250
2500
Epsilon = delta = 0,2
2750
Figura 3.8: Distribuição dos parâmetros de Thomsen no modelo para o caso marinho,
Figura 3.6.
83
Superfície (m)
200
400
600
Superfície (m)
800 1000 1200 1400 1600 1800 2000 2200 2400
0
50
100
150
200
250
300
Profundidade (m)
Profundidade (m)
0
350
400
450
500
550
600
650
700
0
200
400
600
800 1000 1200 1400 1600 1800 2000 2200 2400
0
200
400
600
800 1000 1200 1400 1600 1800 2000 2200 2400
0
200
400
600
800 1000 1200 1400 1600 1800 2000 2200 2400
0
200
400
600
800 1000 1200 1400 1600 1800 2000 2200 2400
0
50
100
150
200
250
300
350
400
450
500
550
600
650
700
Superfície (m)
200
400
600
Superfície (m)
800 1000 1200 1400 1600 1800 2000 2200 2400
0
50
100
150
200
250
300
Profundidade (m)
Profundidade (m)
0
350
400
450
500
550
600
650
700
0
50
100
150
200
250
300
350
400
450
500
550
600
650
700
Superfície (m)
200
400
600
Superfície (m)
800 1000 1200 1400 1600 1800 2000 2200 2400
0
50
100
150
200
250
300
Profundidade (m)
Profundidade (m)
0
350
400
450
500
550
600
650
700
0
50
100
150
200
250
300
350
400
450
500
550
600
650
700
Superfície (m)
200
400
600
Superfície (m)
800 1000 1200 1400 1600 1800 2000 2200 2400
0
50
100
150
200
250
300
Profundidade (m)
Profundidade (m)
0
350
400
450
500
550
600
650
700
0
50
100
150
200
250
300
350
400
450
500
550
600
650
700
Figura 3.9: Instantâneos relativos ao componente vertical do campo de velocidades obtidos
na simulação sísmica no modelo definido na Figura 3.6.
84
Superfície (m)
200
400
600
Superfície (m)
800 1000 1200 1400 1600 1800 2000 2200 2400
0
50
100
150
200
250
300
Profundidade (m)
Profundidade (m)
0
350
400
450
500
550
600
650
700
0
200
400
600
800 1000 1200 1400 1600 1800 2000 2200 2400
0
200
400
600
800 1000 1200 1400 1600 1800 2000 2200 2400
0
200
400
600
800 1000 1200 1400 1600 1800 2000 2200 2400
0
200
400
600
800 1000 1200 1400 1600 1800 2000 2200 2400
0
50
100
150
200
250
300
350
400
450
500
550
600
650
700
Superfície (m)
200
400
600
Superfície (m)
800 1000 1200 1400 1600 1800 2000 2200 2400
0
50
100
150
200
250
300
Profundidade (m)
Profundidade (m)
0
350
400
450
500
550
600
650
700
0
50
100
150
200
250
300
350
400
450
500
550
600
650
700
Superfície (m)
200
400
600
Superfície (m)
800 1000 1200 1400 1600 1800 2000 2200 2400
0
50
100
150
200
250
300
350
400
450
500
550
Profundidade (m)
Profundidade (m)
0
600
650
700
0
50
100
150
200
250
300
350
400
450
500
550
600
650
700
Superfície (m)
200
400
600
Superfície (m)
800 1000 1200 1400 1600 1800 2000 2200 2400
0
50
100
150
200
250
300
Profundidade (m)
Profundidade (m)
0
350
400
450
500
550
600
650
700
0
50
100
150
200
250
300
350
400
450
500
550
600
650
700
Figura 3.10: Instantâneos relativos ao componente horizontal do campo de velocidades
obtidos na simulação sísmica no modelo definido na Figura 3.6.
85
Superfície (m)
200
400
600
Superfície (m)
800 1000 1200 1400 1600 1800 2000 2200 2400
0
50
100
150
200
250
300
Profundidade (m)
Profundidade (m)
0
350
400
450
500
550
600
650
700
0
200
400
600
800 1000 1200 1400 1600 1800 2000 2200 2400
0
200
400
600
800 1000 1200 1400 1600 1800 2000 2200 2400
0
200
400
600
800 1000 1200 1400 1600 1800 2000 2200 2400
0
200
400
600
800 1000 1200 1400 1600 1800 2000 2200 2400
0
50
100
150
200
250
300
350
400
450
500
550
600
650
700
Superfície (m)
200
400
600
Superfície (m)
800 1000 1200 1400 1600 1800 2000 2200 2400
0
50
100
150
200
250
300
Profundidade (m)
Profundidade (m)
0
350
400
450
500
550
600
650
700
0
50
100
150
200
250
300
350
400
450
500
550
600
650
700
Superfície (m)
200
400
600
Superfície (m)
800 1000 1200 1400 1600 1800 2000 2200 2400
0
50
100
150
200
250
300
350
400
450
500
550
Profundidade (m)
Profundidade (m)
0
600
650
700
0
50
100
150
200
250
300
350
400
450
500
550
600
650
700
Superfície (m)
200
400
600
Superfície (m)
800 1000 1200 1400 1600 1800 2000 2200 2400
0
50
100
150
200
250
300
Profundidade (m)
Profundidade (m)
0
350
400
450
500
550
600
650
700
0
50
100
150
200
250
300
350
400
450
500
550
600
650
700
Figura 3.11: Instantâneos relativos ao componente τ xx do campo de tensões obtidos na
simulação sísmica no modelo definido na Figura 3.6.
86
Superfície (m)
200
400
600
Superfície (m)
800 1000 1200 1400 1600 1800 2000 2200 2400
0
50
100
150
200
250
300
Profundidade (m)
Profundidade (m)
0
350
400
450
500
550
600
650
700
0
200
400
600
800 1000 1200 1400 1600 1800 2000 2200 2400
0
200
400
600
800 1000 1200 1400 1600 1800 2000 2200 2400
0
200
400
600
800 1000 1200 1400 1600 1800 2000 2200 2400
0
200
400
600
800 1000 1200 1400 1600 1800 2000 2200 2400
0
50
100
150
200
250
300
350
400
450
500
550
600
650
700
Superfície (m)
200
400
600
Superfície (m)
800 1000 1200 1400 1600 1800 2000 2200 2400
0
50
100
150
200
250
300
Profundidade (m)
Profundidade (m)
0
350
400
450
500
550
600
650
700
0
50
100
150
200
250
300
350
400
450
500
550
600
650
700
Superfície (m)
200
400
600
Superfície (m)
800 1000 1200 1400 1600 1800 2000 2200 2400
0
50
100
150
200
250
300
Profundidade (m)
Profundidade (m)
0
350
400
450
500
550
600
650
700
0
50
100
150
200
250
300
350
400
450
500
550
600
650
700
Superfície (m)
200
400
600
Superfície (m)
800 1000 1200 1400 1600 1800 2000 2200 2400
0
50
100
150
200
250
300
Profundidade (m)
Profundidade (m)
0
350
400
450
500
550
600
650
700
0
50
100
150
200
250
300
350
400
450
500
550
600
650
700
Figura 3.12: Instantâneos relativos ao componente τ xz do campo de tensões obtidos na
simulação sísmica no modelo definido na Figura 3.6. Como era de se esperar a tensão
cisalhante é nula na parte acústica, 0,5 e 1,0s.
87
Superfície (m)
200
400
600
Superfície (m)
800 1000 1200 1400 1600 1800 2000 2200 2400
0
50
100
150
200
250
300
Profundidade (m)
Profundidade (m)
0
350
400
450
500
550
600
650
700
0
200
400
600
800 1000 1200 1400 1600 1800 2000 2200 2400
0
200
400
600
800 1000 1200 1400 1600 1800 2000 2200 2400
0
200
400
600
800 1000 1200 1400 1600 1800 2000 2200 2400
0
200
400
600
800 1000 1200 1400 1600 1800 2000 2200 2400
0
50
100
150
200
250
300
350
400
450
500
550
600
650
700
Superfície (m)
200
400
600
Superfície (m)
800 1000 1200 1400 1600 1800 2000 2200 2400
0
50
100
150
200
250
300
Profundidade (m)
Profundidade (m)
0
350
400
450
500
550
600
650
700
0
50
100
150
200
250
300
350
400
450
500
550
600
650
700
Superfície (m)
200
400
600
Superfície (m)
800 1000 1200 1400 1600 1800 2000 2200 2400
0
50
100
150
200
250
300
350
400
450
500
550
Profundidade (m)
Profundidade (m)
0
600
650
700
0
50
100
150
200
250
300
350
400
450
500
550
600
650
700
Superfície (m)
200
400
600
Superfície (m)
800 1000 1200 1400 1600 1800 2000 2200 2400
0
50
100
150
200
250
300
Profundidade (m)
Profundidade (m)
0
350
400
450
500
550
600
650
700
0
50
100
150
200
250
300
350
400
450
500
550
600
650
700
Figura 3.13: Instantâneos relativos ao componente τ zz do campo de tensões obtidos na
simulação sísmica no modelo definido na Figura 3.6.
88
Superfície (m)
0
2000
4000
6000
8000
10000
12000
14000
16000
18000
20000
22000
24000
0
0.8
1.6
Tempoderegistro(s)
2.4
3.2
4.0
4.8
5.6
6.4
7.2
8.0
Figura 3.14: Sismogramas relativos ao componente vertical do campo de velocidades
obtidos na simulação sísmica no modelo definido na Figura 3.6.
Superfície (m)
0
2000
4000
6000
8000
10000
12000
14000
16000
18000
20000
22000
24000
0
0.8
1.6
Tempoderegistro(s)
2.4
3.2
4.0
4.8
5.6
6.4
7.2
8.0
Figura 3.15: Sismogramas relativos ao componente horizontal do campo de velocidades
obtidos na simulação sísmica no modelo definido na Figura 3.6.
89
Superfície (m)
0
2000
4000
6000
8000
10000
12000
14000
16000
18000
20000
22000
24000
0
0.8
1.6
Tempoderegistro(s)
2.4
3.2
4.0
4.8
5.6
6.4
7.2
8.0
Figura 3.16: Sismogramas relativos ao componente τ xx do campo de tensões obtidos na
simulação sísmica no modelo definido na Figura 3.6.
Superfície (m)
0
2000
4000
6000
8000
10000
12000
14000
16000
18000
20000
22000
24000
0
0.8
1.6
Tempoderegistro(s)
2.4
3.2
4.0
4.8
5.6
6.4
7.2
8.0
Figura 3.17: Sismogramas relativos ao componente τ xz do campo de tensões obtidos na
simulação sísmica no modelo definido na Figura 3.6. Como era de se esperar a tensão
cisalhante é nula na parte acústica.
90
Superfície (m)
0
2000
4000
6000
8000
10000
12000
14000
16000
18000
20000
22000
24000
0
0.8
1.6
Tempoderegistro(s)
2.4
3.2
4.0
4.8
5.6
6.4
7.2
8.0
Figura 3.18: Sismogramas relativos ao componente τ zz do campo de tensões obtidos na
simulação sísmica no modelo definido na Figura 3.6.
91
Conclusão, Comentários e Trabalhos
Futuros
Foram desenvolvidos dois métodos para modelagens elásticas por diferenças
finitas visando abordagens em métodos sísmicos realísticos que apresentam anisotropia
do tipo TIV.
O primeiro método é uma generalização do algoritmo de Zahradník e Priolo
(1994), no qual os parâmetros elásticos são introduzidos por integrações ao longo das
linhas da malha, neste trabalho suposta regular. Este método é indicado para simulações
sísmicas terrestres e não está restrito a superfícies de observação horizontais.
O fato de a malha ser definida através da condição de não dispersão que
envolve o menor valor de velocidade de ondas cisalhantes, faz com que este método não
seja adequado para modelagens marinhas, pois esta velocidade é igual a zero na parte
líquida, o que implica uma distância nula entre os pontos consecutivos da malha.
O segundo método amplia as possibilidades da proposta apresentada em
Levander (1990) para a realização de modelagens elásticas em meios transversalmente
isotrópicos. Este método tem a conveniência de poder ser aplicado em meios marinhos,
92
pois a condição de não dispersão do mesmo envolve a velocidade das ondas compressionais
exclusivamente.
Quanto ao tipo de levantamento, ambos os métodos não fazem qualquer
restrição à geometria de aquisição, ou seja, todos os tipos de levantamentos empregados
atualmente, por exemplo, superfície-superfície com uma ou duas componentes, tomografia
interpoços, VSP, aquisições marinhas com cabos de fundo ou verticais, podem ser simulados facilmente.
Foram realizadas cinco aplicações.
A primeira aplicação objetivou o estudo do comportamento ondulatório em
meios transversalmente isotrópicos homogêneos, definidos com valores das velocidades das
ondas compressionais e cisalhantes e parâmetros de Thomsen (1986).
O comportamento da frente de onda P na direção horizontal em função da
maior velocidade de propagação nessa direção e a frente de onda S na direção vertical,
dependem de como varia o parâmetro de Thomsen ε nas litologias estudadas.
O parâmetro δ se traduz por uma “aceleração ”, para δ > 0, ou “desaceleração”, para δ < 0, da frente de onda P nas direções diagonais, enquanto que, para as
ondas S, ocorre uma “desaceleração”, para δ > 0 e “aceleração”, para δ < 0 (Rosa Filho
2002 [29]).
A segunda aplicação exibiu o comportamento de diversos tipos de ondas
refletidas numa interface curva que separa um meio transversalmente isotrópico, com
vários valores de parâmetros de Thomsen, de um meio isotrópico.
Observou-se um aumento gradual da elipsidade no que diz respeito à propagação das ondas compressionais (quase compressionais, rigorosamente), com o aumento
93
dos valores de ε e δ.
Observou-se também um resquício de onda S originado por não se usar uma
fonte exploxiva exata.
Ao contrário da frente de ondas compressionais, as ondas cisalhantes exibem
um caráter mais isotrópico, independente dos valores de ε e δ assumidos nesta aplicação.
Como esperado, ocorreram variações consideráveis no tempo de trânsito das
ondas PP nos offsets maiores.
As ondas SP também são modificadas com o valor crescente de ε e δ e as
ondas PS e SS praticamente não são alteradas com a variação dos parâmetros de Thomsen
usados nesta aplicação.
A terceira aplicação exibiu a capacidade do primeiro método em simulações
em modelos geológicos complexos. Instantâneos e sismogramas relativos aos componentes
do campo de deslocamento mostram as diferenças entre a propagação de ondas em meios
isotrópicos e moderadamente anisotrópicos (parâmetros de Thomsen iguais a 0,2).
Observou-se um espalhamento lateral de energia superior na propagação na
camada anisotrópica, também constatado nos maiores offsets. Tal fato tem o potencial de
degradar o imageamento sísmico quando o caráter anisotrópico da propagação do campo
de onda não for levado em consideração na estratégia de processamento dos dados.
Foi detectada a onda tipo Rayleigh encontrada com a implementação do
Formalismo do Vácuo na superfície que delimita o modelo na parte superior.
A quarta aplicação mostrou a estabilidade do segundo método em meios
acústicos e a quinta exibiu a capacidade deste método em modelagens marinhas complexas.
94
Este método pode ser usado para simular melhor a propagação de onda
em meios elásticos e acústicos alternados, tem bom comportamento em meios líquidos,
onde a velocidade da onda S vai para zero, e não é necessário tratamento especial para a
interface líquido-sólido, sendo ele ideal para simulação sísmica marítima, ou seja é estável
para contatos rocha-fluido [4].
Em relação a trabalhos futuros:
Comparar os dois métodos em modelagens terrestres (aplicar as mesmas
condições de contorno);
Generalizar as metodologias para abordagens em meios transversalmente
isotrópicos com eixo de simetria variável;
Estender os algoritmos para simulações tridimensionais e implementá-los em
ambientes de “cluster” de computadores pessoais;
Modificar as condições de contorno na superfície para incluir o operador
cheio encontrado em Zahradník e Priolo (1994) e comparar com a condição utilizada
neste tabalho;
Buscar estratégias para separação dos campos de onda compressional e cisalhante visando migrações sísmicas elásticas anisotrópicas por diferenças finitas.
95
Referências
[1] SAVASTA, Olga., KLÍE, H., TORO, W., et al. Estudios de anisotropía sísmica en
geofísica de exploración y producción. Vision Tecnologica, v.8, n.1, 2000.
[2] ZAHRADNÍK, J., PRIOLO, E. Heterogeneous formulatinos of elastodynamic equations and finite-difference schemes. Geophysics, Tulsa, July,1994.
[3] ZAHRADNÍK, J., MOCZO, P., HRON, F. Testing four elastic finite-difference
schemes for behavior at discontinuities. Bull. Seism. Soc. Am., v.83, p. 107-129.
[4] LEVANDER,
A.
R.
Fourth-order
finite-difference
P-SV
seismograms.
In:_.Numerical modeling of seismic wave propagation. Tulsa:
Society
of Exploration Geophysicists, 1990. p.190 - 201. ( Geophysics Reprint Series, n.13).
[5] THOMSEN, L. Weak elastic anisotropy. Geophysics, Tulsa, v.51, n.10, p.1954—1966,
oct.,1986.
[6] CUNHA, P. E. M. Estratégias eficientes para migração reversa no tempo préempilhamento 3-D em profundidade pelo método das diferenças finitas.
1997. 134f. Dissertação (Mestrado em Ciências em Geofísica). __ . Programa de
96
Pesquisa de Pós-Graduação em Geofísica, Universidade Federal da Bahia, Bahia,
1997.
[7] CERJAN, C., KOSLOFF, D., KOSLOFF, R., et al. A Nonreflecting bondary condition for discrete acoustic and elastic wave equations. In:_.Numerical modeling
of seismic wave propagation. Tulsa: Society of Exploration Geophysicists, 1990.
p.444 - 447. ( Geophysics Reprint Series, n.13).
[8] BOORE, D.M. Finite - difference methods for seismic waves in computational
physics, New York: Academic Press, 1972. v.11, p.1-37.
[9] AKI, K., RICHARDS, P. Quantitative seismology: theory and methods. W. H.
Freeman & CO, 1980.
[10] ERMERMAN, S. H., STEPHEN, R. A. Comment on absorbing boundary conditions
for acoustic and elastic wave equations. In:_.Numerical modeling of seismic
wave propagation. Tulsa: Society of Exploration Geophysicists, 1990. p.463 - 467.
( Geophysics Reprint Series, n.13).
[11] CLAYTON, R., ENGQUIST, B. Absorbing boundary conditions for acoustic and
elastic wave equations. In:_.Numerical modeling of seismic wave propagation.
Tulsa: Society of Exploration Geophysicists, 1990. p.448 - 459. ( Geophysics Reprint
Series, n.13).
[12] REYNOLDS, A. C. Boundary conditions for the numerical solution of wave propagation problems. Geophysics, v 43, p.1099-1110, 1978.
97
[13] SOARES FILHO, D. M., ROSALBA. J. F., ALMEIDA, M. A., et al. Seismic modeling with complex geology and topography. In: INTERNATIONAL CONGRESS
OF THE BRAZILIAN GEOPHYSICAL SOCIETY, 7 ,2001. Bahia. Resumos Expandidos... Bahia:. SBGf33899.
[14] ALFORD, R. M., KELLY, K. R., BOORE, D. M. Accuracy of finite-difference modeling of the acoustic wave equation. In:_.Numerical modeling of seismic wave
propagation. Tulsa: Society of Exploration Geophysicists, 1990. p.310 - 318. ( Geophysics Reprint Series, n.13).
[15] KELLY, K. R., WARD, R. W., TREITEL, S., et al. Synthetic seismograms: a finitediffrence approach. In:_.Numerical modeling of seismic wave propagation.
Tulsa: Society of Exploration Geophysicists, 1990. p.36 - 61. ( Geophysics Reprint
Series, n.13).
[16] TATHAM, R. H., McCORMACK, M. D. Multicomponent seismology in
petroleum exploration. Tulsa: Society of Exploration Geophysicists, 1991. (Investigations in Geophysics n.6).
[17] CRAMPIN, S. Suggestions for a consistent terminology for seismic anisotropy. Geophysical Prospecting, v.37, p.753—770, 1989.
[18] GEOLTRAIN, S. Propagation of elastic waves in transversely isotropic media. CWP074, Colorado, sept., 1988.
[19] ALTERMAN, Z., KARAL Jr, F. C. Propagation of elastic waves in layered media
by finite difference methods. In:_.Numerical modeling of seismic wave prop98
agation. Tulsa: Society of Exploration Geophysicists, 1990. p.4 - 35. ( Geophysics
Reprint Series, n.13).
[20] VIRIEUX, J. SH- wave propagation in heterogeneous media: Velocity-stress finitedifference method. In:_.Numerical modeling of seismic wave propagation.
Tulsa: Society of Exploration Geophysicists, 1990. p.167 - 176. ( Geophysics Reprint
Series, n.13).
[21] VIRIEUX, J. P-SV wave propagation in heterogeneous media: velocity-stress finitedifference method. In:_.Numerical modeling of seismic wave propagation.
Tulsa: Society of Exploration Geophysicists, 1990. p.177 - 189. ( Geophysics Reprint
Series, n.13).
[22] BAYLISS, A., JORDAN, K. E., LEMESURIER, B. J., et al. A Fourth-order accurate
finite-difference sheme for the computation of elastic waves. In:_.Numerical modeling of seismic wave propagation. Tulsa: Society of Exploration Geophysicists,
1990. p.202 - 219. ( Geophysics Reprint Series, n.13).
[23] MUFTI, I. R. Seismic modeling in the implicit mode. In:_.Numerical modeling
of seismic wave propagation. Tulsa: Society of Exploration Geophysicists, 1990.
p.220 - 258. ( Geophysics Reprint Series, n.13).
[24] STEPHEN, R. A. A Comparison of finite difference and reflectivity seimograms for
marine models. In:_.Numerical modeling of seismic wave propagation. Tulsa:
Society of Exploration Geophysicists, 1990. p.391- 409. ( Geophysics Reprint Series,
n.13).
99
[25] TRIGGIA, A. A; et al, Fundamentos de engenharia de petróleo. Rio de Janeiro:
Interciência/PETROBRAS, 2001
[26] PASSOS, G. F. Modelagem sísmica em meios efetivos simulando reservatório
de hidrocarbonetos. 2001. 116f. Dissertação (Mestrado em Ciências em Geofísica).
__ . Programa de Pesquisa de Pós-Graduação em Geofísica, Universidade Federal
da Bahia, Bahia, 2001.
[27] CHATTERJEE, R. Mathematical theory of continuum mechanics. Índia:
Narosa Publishing House, 1999.
[28] FARIA, E. L. Modeling migration and focusing analysis in transversely
isotropic media. University of Texas at Austin, Texas 1993.
[29] ROSA FILHO, J. C. Modelagem sísmica de ondas elásticas e migração reversa no tempo em meios transversalmente isotrópico. 2002. Dissertação (Mestrado em Engenharia Civil). __ . Coordenação dos Programas de Pós-Graduação em
Engenharia, Universidade Federal do Rio de Janeiro, Rio de Janeiro, 2002.
[30] NYE, J.F. Physical properties of crystals: Oxford Press, 1957.
[31] MUFTI, I. R., PITA, J. A. HUNTLEY, R. W. Finite-difference depth migration of
exploration-scale 3-D seismic data. Geophysics, v.61, p.776-794, 1996.
100
Apêndices
101
Apêndice A
Tensor de Elasticidade (Cijkl )
Expresso no Sistema de Coordenadas
Principal
As deformações em posições, orientações e forma dos elementos de um corpo
se relacionam com as forças superficiais aplicadas sobre os mesmos através das constantes
dos materiais (constantes elásticas) e são relacionados a partir da Lei de Hooke, Equação
A.1.
σ ij = bij + Cijkl εkl ,
onde σ ij = 0, quando todos εij = 0, então bij = 0.
102
(A.1)
A.1
Componentes da tensão e deformação
A teoria elástica clássica assume uma generalização da Lei de Hooke ex-
pressando cada componente do tensor de tensão (σ ij ) como uma combinação linear do
tensor de deformação (εkl ) [30]. Onde cada índice de direção pode assumir valores 1,2,3
(representando direções x, y, z), existem nove relações semelhantes, cada uma envolvendo
componentes de tensão e nove componentes de deformação. Estas nove equações podem
ser escritas compactamente como a equação A.2, com 81 coeficientes mas nem todos são
independentes. Os coeficientes (Cijkl ) são chamados de constantes elásticas ou tensor módulos elásticos de corpo e caracterizam completamente as propriedades elásticas do corpo,
Cijkl é um tensor simétrico de quarta ordem, (3x3x3x3), Equação A.2.
σ ij =
3 X
3
X
Cijkl εkl , (i, j, k, l = 1, 2, 3).
k=1 l=1
Forma matricial da Lei de Hooke:
103
(A.2)



































































σ 11  
 
 

σ 21 
 
 
 

σ 31 
 
 
 

σ 41 
 
 
 

σ 51 
=
 
 

σ 61 
 
 
 
 
σ 71  
 
 
 
σ 81  
 
 
σ 91

σ 11  
 
 

σ 21 
 
 
 

σ 31 
 
 
 

σ 41 
 
 
 

σ 51 
=
 
 

σ 61 
 
 
 
 
σ 71  
 
 
 
σ 81  
 
 
σ 91

C11 C12 C13 C14 C15 C16 C17 C18 C19  



C21 C22 C23 C24 C25 C26 C27 C28 C29 




C31 C32 C33 C34 C35 C36 C37 C38 C39 




C41 C42 C43 C44 C45 C46 C47 C48 C49 




C51 C52 C53 C54 C55 C56 C57 C58 C59 




C61 C62 C63 C64 C65 C66 C67 C68 C69 




C71 C72 C73 C74 C75 C76 C77 C78 C79  



C81 C82 C83 C84 C85 C86 C87 C88 C89  


C91 C92 C93 C94 C95 C96 C97 C98 C99

ε11 


ε21 



ε31 



ε41 



ε51 
,


ε61 




ε71 



ε81 


ε91

C11 ε11 C12 ε21 C13 ε31 C14 ε41 C15 ε51 C16 ε61 C17 ε71 C18 ε81 C19 ε91 


C21 ε11 C22 ε21 C23 ε31 C24 ε41 C25 ε51 C26 ε61 C27 ε71 C28 ε81 C29 ε91 



C31 ε11 C32 ε21 C33 ε31 C34 ε41 C35 ε51 C36 ε61 C37 ε71 C38 ε81 C39 ε91 



C41 ε11 C42 ε21 C43 ε31 C44 ε41 C45 ε51 C46 ε61 C47 ε71 C48 ε81 C49 ε91 



C51 ε11 C52 ε21 C53 ε31 C54 ε41 C55 ε51 C56 ε61 C57 ε71 C58 ε81 C59 ε91 
,


C61 ε11 C62 ε21 C63 ε31 C64 ε41 C65 ε51 C66 ε61 C67 ε71 C68 ε81 C69 ε91 




C71 ε11 C72 ε21 C73 ε31 C74 ε41 C75 ε51 C76 ε61 C77 ε71 C78 ε81 C79 ε91 



C81 ε11 C82 ε21 C83 ε31 C84 ε41 C85 ε51 C86 ε61 C87 ε71 C88 ε81 C89 ε91 


C91 ε11 C92 ε21 C93 ε31 C94 ε41 C95 ε51 C96 ε61 C97 ε71 C98 ε81 C99 ε91
104
onde a simetria dos tensores de tensão
σ ij = σ ji ,
Cijkl εkl = Cjikl εkl ,
(A.3)
então,
Cijkl = Cjikl .
Agora Cijkl pode ser escrita como uma soma
0
00
Cijkl = Cijkl + Cijkl ,
onde
1
0
0
Cijkl = (Cijkl + Cjikl ) = Cijlk ,
2
(A.4)
1
00
00
Cijkl = (Cijkl − Cjikl ) = −Cijlk ,
2
(A.5)
assim podemos reescrever σ ij como,
0
00
σ ij = Cijkl εkl + Cijkl εkl ,
105
(A.6)
usando álgebra indicial, temos,
00
00
Cijkl εkl = Cijlk εlk ,
(A.7)
00
= −Cijkl εlk ,
onde
00
00
Cijkl = −Cijlk ,
(A.8)
00
= −Cijkl εkl ,
e
εkl = εlk ,
(A.9)
portanto
00
Cijkl εkl = 0,
isto facilita a redução da equação A.6 para
0
σ ij = Cijkl εkl ,
106
(A.10)
onde
0
0
Cijkl = Cijlk ,
(A.11)
combinando a equação A.3 com a equação A.11,
σ ij = Cijkl εkl ,
(A.12)
Cijkl é simétrico em relação aos índices ij e kl.
Temos,
Cijkl = Cjikl = Cijlk = Cjilk .
(A.13)
A simetria dos tensores de deformação
(εkl = εlk ) ,
reduz para seis os termos da direita de cada equação independente A.2, assim, reduz de
nove para seis equações independentes. Portanto reduzindo para 36 constantes distintas
[27].
A forma mais geral do tensor de constante elástica de quarta ordem é,
Cijkl = λδ ij δ kl + µ (δ ik δ jl + δ il δ jk ) + ν (δ ik δjl − δ il δjk ) ,
107
(A.14)
onde λ, µ, ν são constantes elásticas.
A elasticidade pode ser representada mais compactamente com as trocas de
índices ij por α e kl por β, temos:
11 → 1,
22 → 2,
33 → 3,
32 = 23 → 4,
31 = 13 → 5
e
12 = 21 → 6,
assim o tensor Cijkl (3x3x3x3) pode ser representado como uma matrix Cαβ (6x6).
108
Pela reciprocidade, isto, é:
Cαβ = Cβα ,
assim reduzindo para 21 constantes.

Cαβ









=










C11 C12 C13 C14 C15 C16 


0 C22 C23 C24 C25 C26 



0
0 C33 C34 C35 C36 

.

0
0
0 C44 C45 C46 




0
0
0
0 C55 C56 


0
0
0
0
0 C66
Considerando um material elasticamente simétrico com relação a um eixo.
Se for aplicado uma rotação do sistema xyz sobre o eixo z com um ângulo qualquer, temos
que:
C14 = C15 = C24 = C25 = C34 = C35 = C46 = C56 = 0,
e
C16 = C26 = C36 = C45 = 0,
109
assim a matrix para Cαβ

Cαβ









=









C11 C11 − 2C66 C13
0
0
0
C11
C13
0
0
0
0
C33
0
0
0
0
0
C44
0
0
0
0
0
C44
0
0
0
0
0

0 


0 



0 

,

0 




0 


C66
(A.15a)
sendo que o número de constantes elásticas se reduz para 5.
C11 : controla a propagação de ondas P no eixo horizontal.
C33 : controla a velocidade vertical de ondas P.
C44 : controla a propagação vertical e horizontal das ondas SV e a propagação vertical das ondas SH.
C13 : controla a propagação oblíqua de ondas P e SV.
C66 : controla a propagação horizontal da onda SH.
Os parâmetros anisotrópicos de Thomsen em função dos parâmetros elásticos são:
² = (C11 − C33 ) /2C33 ,
¤
£
2
,
τ = 2 (C13 + C44 )2 − (C33 − C44 ) (C11 + C33 − 2C44 ) /2C33
γ = (C66 − C44 ) /2C44 ,
110
(A.16)
ou com os prâmetros elásticos em função dos parâmetros anisotrópicos:
C13 = −C44 ±
C33 = ρVP2 ,
(A.17)
C44 = ρVS2 ,
(A.18)
C11 = (2² + 1)C33,
(A.19)
q
2
2
C44
+ (2τ + 1)C33
− (2τ + 2)C33 C44 ,
(A.20)
C66 = (2γ + 1)C44 .
Os parâmetros anisotrópicos estão relacionados com a anisotropia para ângulos polares pequenos, com a dependência angular da componente de propagação horizontal da onda P e com a anisotropia da onda SH, respectivamente.
Quando τ = ² há a anisotropia elipsoidal da onda P.
As velocidades verticais e as de fase das ondas P e S que são
expressas por,
VP = (C33 /ρ)1/2 ,
111
(A.21)
VS = (C44 /ρ)1/2 ,
(A.22)
VP2 (θ) = α20 [1+ ∈ sin2 θ + D (θ)],
(A.23)
2
(θ) = β 20 [1 + α20 ∈ sin2 θ/β 20 − α20 D (θ) β 20 ],
VSV
(A.24)
2
VSH
(θ) = β 20 [1 + 2γ sin2 θ],
(A.25)
onde ρ é a densidade, θ o ângulo de fase e,
¡
¢
D(θ) = 1 − β 20 /α20 [1 + 4τ sin2 θ cos2 θ/(1 − β 20 /α20 )2
¡
¢2
+ 4(1 − β 20 /α20 + ²)² sin4 θ/ 1 − β 20 /α20 ]1/2 − 1.
(A.26)
As velocidades de fase (Thonsem, 1986) VP , VSV e VSH representam as
velocidades de propagação da frente de onda, que diferem da propagação de energia (velocidade de grupo).
112
Apêndice B
Discretização por Diferenças Finitas
B.1
Malhas Intercaladas
Para se utilizar a formulação de malha intercalada, com o intuito de aumentar a precisão
dos cálclulos e a estabilidade do algoritmo, expandi-se a função u (x + n∆x) em série de
Taylor com valores fracionários.
¡
¢
¡
¢
Assim para u x + 12 ∆x e u x − 12 ∆x , temos respectivamente,
µ
¶
¶2
µ
1 (2)
1
∆x (1)
∆x
u x + ∆x = u(x) +
u (x) +
u (x)
2
2
2
2!
¶3
µ
1 (3)
∆x
u (x) + ...,
+
2
3!
µ
µ
¶
¶2
1
∆x (1)
1 (2)
∆x
u x − ∆x = u(x) −
u (x) +
u (x)
2
2
2
2!
¶3
µ
∆x
1 (3)
u (x) + ...,
−
2
3!
113
(B.1)
(B.2)
¢
¡
¢
¡
Fazendo para u x + 32 ∆x e u x − 32 ∆x , temos,
µ
¶
¶2
µ
1 (2)
3
3∆x (1)
3∆x
u x + ∆x = u(x) +
u (x) +
u (x)
2
2
2
2!
¶3
µ
1 (3)
3∆x
u (x) + ...,
+
2
3!
µ
µ
¶
¶2
3
3∆x (1)
1 (2)
3∆x
u x − ∆x = u (x) −
u (x) +
u (x)
2
2
2
2!
¶3
µ
3∆x
1 (3)
u (x) + ...,
−
2
3!
(B.3)
(B.4)
Subtraindo-se B.2 de B.1,
µ
¶
¶
µ
1
1
u x + ∆x − u x − ∆x = ∆xu(1) (x)
2
2
¶3
µ
¡
¢
2 (3)
∆x
u (x) + O0 ∆x5 ,
+
2
3!
(B.5)
Subtraindo-se B.4 de B.3,
µ
¶
µ
¶
3
3
u x + ∆x − u x − ∆x = 3∆xu(1) (x) +
2
2
µ
¶3
¢
3∆x
2 (3)
00 ¡
u (x) + O ∆x5 ,
2
3!
(B.6)
Fazendo-se [((B.5) x 27) - B.6] e desprezando os termos de ordem superior a três,
· µ
¶
µ
¶¸ · µ
¶
µ
¶¸
1
3
1
3
27 u x + ∆x − u x − ∆x − u x + ∆x − u x − ∆x
2
2
2
2
= 27∆xu(1) (x) − 3∆xu(1) (x) ,
114
(B.7)
Isolando-se a primeira derivada u(1) (x) =
∂u(x)
∂x
em relação a x em B.7 tem-se uma aprox-
imação,
∂u(x)
1
=
∂x
∆x
½ · µ
¶
µ
¶¸
· µ
¶
µ
¶¸¾
9
1
1
1
3
3
u x + ∆x − u x − ∆x −
u x + ∆x − u x − ∆x
,
8
2
2
24
2
2
(B.8)
115