Aspectos da modelagem acústica usando elementos finitos espectrais

Prévia do material em texto

UNIVERSIDADE FEDERAL DO PARÁ
INSTITUTO DE GEOCIÊNCIAS
PROGRAMA DE PÓS-GRADUAÇÃO EM GEOFÍSICA
DISSERTAÇÃO DE MESTRADO
Aspectos da modelagem acústica usando elementos
finitos espectrais
LUCAS DE CASTRO COSTA
Belém - Pará
2022
LUCAS DE CASTRO COSTA
Aspectos da modelagem acústica usando elementos
finitos espectrais
Dissertação apresentada ao Programa de Pós-Graduação
em Geofísica do Instituto de Geociências da Universidade
Federal do Pará para obtenção do título de Mestre em
Geofísica.
Área de Concentração: Métodos Sísmicos
Linha de pesquisa: Modelagem e inversão de dados
geofísicos
Orientador: Prof. Jessé Carvalho Costa
Belém - Pará
2022
Dados Internacionais de Catalogação na Publicação (CIP) de acordo com ISBD
Sistema de Bibliotecas da Universidade Federal do Pará
Gerada automaticamente pelo módulo Ficat, mediante os dados fornecidos pelo(a) autor(a)
C837a Costa, Lucas de Castro.
 Aspectos da modelagem acústica usando elementos finitos
espectrais / Lucas de Castro Costa. — 2022.
 37 f. : il. color.
 Orientador(a): Prof. Dr. Jessé Carvalho Costa
 Dissertação (Mestrado) - Universidade Federal do Pará,
Instituto de Geociências, Programa de Pós-Graduação em
Geofísica, Belém, 2022.
 1. Elementos finitos espectrais. 2. Polinômio de Legendre.
3. Teorema da reciprocidade. 4. Função de Green. I. Título.
CDD 550
Powered by TCPDF (www.tcpdf.org)
LUCAS DE CASTRO COSTA
Aspectos da modelagem acústica usando elementos finitos
espectrais
Dissertação apresentada ao Programa de Pós-Graduação
em Geofísica do Instituto de Geociências da Universidade
Federal do Pará para obtenção do título de Mestre em
Geofísica.
Data de aprovação: 30 de Setembro de 2022
Banca Examinadora:
 
AGRADECIMENTOS
Quero agradecer a minha família, especialmente minha mãe, a minha namorada, meus
avós e aos meus irmãos por terem me apoiado ao longo dessa minha trajetória pelo
mestrado. Agradeço ao meu orientador, Jessé Carvalho Costa, por ter me auxiliado e
ajudado bastante para a conclusão, mesmo com todos os problemas de tempo por causa
da pandemia. Gostaria de agradecer, também, ao Bruno e a Natiê por me auxiliarem e
serem bastante pacientes para tirarem dúvidas quanto a programação e a escrita.
ABSTRACT
The spectral elements’ method (SEM) allows the numerical computation of wave
propagation in very complex models. It can accurately represent boundary conditions,
ocean floor relief, and topography. This dissertation investigates some aspects of modeling
acoustic wave propagation with SEM in the context of seismic exploration. Through
numerical experiments, we evaluate the accuracy of SEM to simulate wave propagation for
seismic application and its computation costs with the objective to assess the adequacy of
SEM for seismic modeling, imaging, and inversion. The computational complexity of SEM
scales very unfavorably compared to the traditional finite-differences methods. However,
for the generation of accurate synthetic data sets in very complex models, and for the
representation of ocean floor relief in nodes or OBC acquisitions in low-frequency band
FWI, it can be a very competitive method.
Keywords: Spectral finite elements. Polynomial. Reciprocity Theorem. Degree.
Green’s Function.
RESUMO
O método dos elementos espectrais permite o cálculo numérico da propagação de ondas
em modelos muito complexos. Ele pode representar com precisão as condições de contorno,
relevo do fundo do oceano e topografia. Esta dissertação investiga alguns aspectos da
modelagem da propagação de ondas acústicas com método dos elementos espectrais no
contexto da exploração sísmica. Por meio de experimentos numéricos, avaliamos a precisão
do método dos elementos espectrais para simular a propagação de ondas para aplicação
sísmica e seus custos computacionais com o objetivo de avaliar a adequação do método dos
elementos espectrais para modelagem sísmica, imageamento e inversão. A complexidade
computacional do método dos elementos espectrais escala muito desfavoravelmente em
comparação com os métodos tradicionais de diferenças finitas. No entanto, para a geração de
conjuntos de dados sintéticos precisos em modelos muito complexos e para a representação
do relevo do fundo oceânico em nós ou aquisições OBC na FWI de banda de baixa
frequência, pode ser um método competitivo.
Palavras-chaves: Elementos finitos espectrais. Polinômio de Legendre. Teorema da
reciprocidade. Função de Green.
LISTA DE FIGURAS
2.1 Mapeamento de um elemento Ω𝑒 no domínio quadrangular de referência, de
dimensões [−1, 1] × [−1, 1], no método de elementos finitos espectral. Os
campos de onda dentro de cada elemento são determinados por interpolação
de Lagrange a partir do valor dos campos em cada nó. A coordenada
de cada nó é obtida através de pontos da quadratura de Gauss-Lobatto-
Legendre para um dado grau do polinômio interpolador. A figura representa
a configuração de nós para interpolação com polinômios de Lagrange de
grau 2. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 4
3.1 Seção sísmica para o campo de pressão calculada analiticamente, à esquerda;
seção para o campo de pressão obtida usando elementos finitos espectral de
ordem 6, a direita. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 13
3.2 Seções sísmicas com o resíduo entre a solução analítica e calculada usando
elementos finitos espectral de: ordem 6 (à esquerda’) e ordem 7 (à direita). 13
3.3 Modelo de velocidade, escala de cores em m/s. . . . . . . . . . . . . . . . 14
3.4 Modelo de densidade, escala de cores em kg/m3. . . . . . . . . . . . . . . 14
3.5 Instantâneo do campo de onda logo após a ativação da fonte no meio com
velocidade e densidade indicadas nas Figuras 3.3 e 3.4, respectivamente.
Observa-se o padrão de radiação da fonte em uma região homogênea. . . . 14
3.6 Instantâneo do campo de onda após a reflexão da frente de onda na superfície
livre. A velocidade e a densidade do modelo estão indicadas nas Figuras
3.3 e 3.4, respectivamente. . . . . . . . . . . . . . . . . . . . . . . . . . . . 15
3.7 Instantâneo do campo de onda após a frente de onda ter passado pela base
do modelo. As fronteiras absorventes foram eficazes para suprimir reflexões
espúrias. A velocidade e a densidade do modelo estão indicadas nas Figuras
3.3 e 3.4, respectivamente. . . . . . . . . . . . . . . . . . . . . . . . . . . . 15
3.8 Seções sísmicas calculadas usando elementos finitos espectral de: ordem 4
(à esquerda) e ordem 7 (à direita); diferença entre as duas seções(abaixo). . 16
3.9 Modelo de velocidade Marmousi, escala de cores em m/s. . . . . . . . . . 17
3.10 Modelo de densidade Marmousi, escala de cores em kg/m3. . . . . . . . . 17
3.11 Instantâneo do campo de onda após 1.4 s de propagação. Fonte de injeção
de volume nas coordenadas (4500 m, 15 m), com superfície livre. . . . . . 17
3.12 Instantâneo do campo de onda após 1.4 s de propagação. Fonte de injeção
de volume nas coordenadas (4500 m, 15 m), sem superfície livre. . . . . . 18
3.13 Instantâneo do campo de onda após 1.4 s de propagação. Fonte de injeção
de volume, localizada nas coordenadas (1500 m, 100 m), com condição de
superfície livre no topo. Observa-se a qualidade das condições de absorção
CPML, não há efeitos de borda significativos na lateral e na base do modelo. 18
3.14 Campo de pressão: sem superfície livre(topo), com superfície-livre (base). . 19
3.15 Gradiente vertical do campo de pressão: sem superfície livre(topo), com
superfície-livre (base). . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 20
3.16 Gradiente horizontal do campo de pressão: sem superfície livre(topo), com
superfície-livre (base). . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 21
3.17 Relação de reciprocidade entre fontes monopolares(topo);relações de re-
ciprocidade entre fonte monopolar e dipolo vertical(meio); relações de
reciprocidade entre fonte comopolar e dipolo horizontal(base). Simulação
com superfície livre. . . . . . . . . .. . . . . . . . . . . . . . . . . . . . . . 22
3.18 Relação de reciprocidade entre fontes do tipo dipolos verticais(topo);relações
de reciprocidade entre dipolo vertical e dipolo horizontal(meio); relações de
reciprocidade entre dipolos horizontais(base). Simulação com superfície livre. 23
SUMÁRIO
1 INTRODUÇÃO 1
2 METODOLOGIA 3
2.1 ELEMENTOS FINITOS ESPECTRAL . . . . . . . . . . . . . . . . . . . 3
2.1.1 Matrizes de compliância e flutuabilidade . . . . . . . . . . . 7
2.1.2 Injeção das fontes . . . . . . . . . . . . . . . . . . . . . . . . . 8
2.1.3 Construção do sistema linear global . . . . . . . . . . . . . . 9
2.1.4 Condições de fronteira . . . . . . . . . . . . . . . . . . . . . . 10
3 RESULTADOS E DISCUSSÃO 12
4 CONCLUSÕES 24
REFERÊNCIAS 25
APÊNDICE A - RELAÇÕES DE RECIPROCIDADE PARA O CAMPO
DE ONDA 26
1 INTRODUÇÃO
A teoria que descreve a propagação de uma onda mecânica é a base da sismologia e
sísmica de exploração. Essa teoria está baseada na resolução das equações diferenciais
parciais desde o modelo mais simples – onda acústica em um meio com densidade constante
– até o modelo com uma complexidade maior. Saber a solução da equação da onda para um
dado conjunto ou distribuição de fontes em um dado meio é crucial em várias aplicações
dentro da sísmica de exploração: estudos de iluminação para um levantamento sísmico,
amarração poço-sísmica, migração RTM, FWI, por exemplo.
Porém, nem sempre esta solução pode ser obtida facilmente. Na verdade, a não ser em
casos mais simples de meios homogêneos ou meios com complexidade muito baixa, a solução
analítica nunca está ao alcance de nossas mãos. Neste contexto, a procura por soluções
numéricas se torna relevante, o que é comumente conhecido por modelagem numérica. Esta
usa modelos matemáticos para descrever as condições físicas de cenários geológicos usando
equações. Com modelos numéricos, existem técnicas, como métodos de diferenças finitas e
elementos finitos, para aproximar as soluções dessas equações. Experimentos numéricos
podem, então, ser realizados nesses modelos produzindo os resultados que podem ser
interpretados no contexto do processo geológico (Ismail-Zadeh, 2010).
Neste trabalho será solucionada a equação acústica escalar da onda, para modelar o
campo de pressão. Para isso, utilizamos do método de elementos espectrais e diferenças
finitas de segunda ordem para fazer a evolução do campo de pressão.
O método de elementos espectrais é baseada em uma formulação fraca das equações
da onda e combina a flexibilidade do método de elementos finitos com a acurácia de um
método pseudo-espectral (Fichtner (2010)). Aplicado ao problema de propagação da onda,
isto resulta em uma maior precisão na determinação do valor de campo em cada elemento.
No entanto, este método exige um custo computacional elevado se comparado ao custo
exigido pelo método de diferenças finitas.
Neste trabalho, será feita a implementação da equação acústica com densidade variável
utilizando o método de elementos espectrais. Para a análise da precisão de modelagem,
será feita a variação do polinômio de Legendre, que são aplicados para evitar que o campo
de pressão nas bordas do intervalo da célula apresente oscilações devido ao fenômeno de
Runge, provocado quando o grau do polinômio de Legendre possui ordem elevada, como
dito por Fichtner (2010).
Para fazer a comparação da precisão de modelagem acústica, foi considerado o teorema
da reciprocidade, neste é dito que o mesmo sismograma deve ser registrado se as localizações
da fonte e do geofone forem trocadas(Fokkema, 2013).
Por fim, é necessário analisar a precisão e a estabilidade do modelador e se ele está
possuindo um bom funcionamento. Para isso, é importante o conhecimento da função
1
2
de Green para estas equações a fim de construir a função "semi-analítica". Será baseado
nos cálculos realizados em Snieder (1998), desenvolvendo a função de Green para nossas
equações a partir da função de Green 2D para o meio acústico de segunda ordem.
2 METODOLOGIA
2.1 ELEMENTOS FINITOS ESPECTRAL
A solução numérica da equação da onda através do método de elementos finitos
espectrais apresenta maior acurácia que métodos de diferenças finitas na implementação de
superfície livre, no tratamento estável de variações topográficas, no cálculo da amplitude
das reflexões em interfaces que separam domínios com forte contraste de propriedades
físicas (Komatitsch, 2005; Trinh, 2017). Nesta etapa da dissertação implementamos o
método de elementos finitos espectral em malhas cartesianas regulares. O custo desta
implementação é superior ao método de diferenças finitas.
Implementamos a solução numérica para a equação de onda com densidade variável:
𝜅(x)𝜕
2𝑝(x, 𝑡)
𝜕𝑡2
− ∇ ·
(︃
1
𝜌(x)∇𝑝(x, 𝑡)
)︃
− 𝑞(x, 𝑡) + ∇ ·
(︃
1
𝜌(x)f(x, 𝑡)
)︃
= 0 (2.1)
na qual, x representa o vetor posição e 𝑡 o tempo, 𝑝(x, 𝑡) representa o campo de pressão,
𝜅(x) indica o módulo de compressibilidade e 𝜌(x) a densidade do meio; 𝑞(x, 𝑡) indica a
taxa de injeção de volume e f(x, 𝑡) a densidade de fontes dipolares.
No método de elementos finitos espectral o domínio de propagação Ω é particionado em
subdomínios quadrangulares, Ω𝑒. Para representar modelos com geometria complexa, cada
elemento Ω𝑒 é mapeado em um domínio de referência definido pelo quadrado (𝜉1, 𝜉2) ∈
[−1, 1] × [−1, 1], como ilustra a Figura 2.1. Em cada elemento, o campo de onda é
amostrado em pontos definidos pela quadratura de Gauss-Lobatto-Legendre (Komatitsch,
2005). O campo de onda dentro de cada elemento é representado por interpolação de
Lagrange (Deville, 2002),
𝑝(x, 𝑡)|Ω𝑒 ≈
𝑁∑︁
𝑖=0
𝑁∑︁
𝑗=0
𝑝𝑖𝑗(𝑡) 𝑙𝑖(𝑥1(𝜉1)) 𝑙𝑗(𝑥2(𝜉2)) =
𝑁∑︁
𝑖=0
𝑁∑︁
𝑗=0
𝑝𝑖𝑗(𝑡) 𝑙𝑖(𝜉1)𝑙𝑗(𝜉2) , (2.2)
em que (𝑥1, 𝑥2) representam as coordenadas de um ponto no elemento da malha de
discretização do modelo e (𝜉1, 𝜉2) representam as coordenadas locais em que este ponto é
mapeado no domínio de referência; 𝑙𝑖(𝜉1) e 𝑙𝑗(𝜉2) indicam polinômios de Lagrange 1D em
cada direção coordenada e 𝑁 indica o grau do polinômio interpolador.
Os polinômios de Lagrange são especificados em função das coordenadas dos nós,
(𝜉𝑘1 , 𝜉𝑘2 ), no domínio de referência. Por exemplo,
𝑙𝑖(𝜉1) =
𝑁∏︀
𝑘=0, 𝑘 ̸=𝑖
(𝜉1 − 𝜉𝑘1 )
𝑁∏︀
𝑘=0, 𝑘 ̸=𝑖
(𝜉𝑖1 − 𝜉𝑘1 )
. (2.3)
3
4
Observa-se desta definição que
𝑙𝑖(𝜉𝑘1 ) =
⎧⎨⎩ 1, se 𝑘 ̸= 𝑖0, caso contrário , (2.4)
e analogamente para 𝑙𝑗(𝜉2).
Ω
e
1
4
5
6
7
2
3
8
9 1
2
3
4
5
6
7
8
9
Figura 2.1: Mapeamento de um elemento Ω𝑒 no domínio quadrangular de referência, de
dimensões [−1, 1] × [−1, 1], no método de elementos finitos espectral. Os campos de
onda dentro de cada elemento são determinados por interpolação de Lagrange a partir do
valor dos campos em cada nó. A coordenada de cada nó é obtida através de pontos da
quadratura de Gauss-Lobatto-Legendre para um dado grau do polinômio interpolador. A
figura representa a configuração de nós para interpolação com polinômios de Lagrange de
grau 2.
A aproximação numérica da solução representada pelo campo nos nós de cada elemento,
𝑝𝑖𝑗(𝑡), é obtida usando o método de Galerkin(Deville, 2002). Esta abordagem exige que
resíduo da aproximação numérica da solução seja ortogonal a todas as funções base,
𝑙𝑖(𝑥1)𝑙𝑗(𝑥2), ou seja,
∫︁
Ω
𝑑x𝑙𝑖(𝑥1)𝑙𝑗(𝑥2)
[︃
𝜅(x)𝜕
2𝑝(x, 𝑡)
𝜕𝑡2
− ∇ ·
(︃
1
𝜌(x)∇𝑝(x, 𝑡)
)︃
− 𝑞(x, 𝑡) + ∇ ·
(︃
1
𝜌(x)f(x, 𝑡)
)︃]︃
= 0(2.5)
para todas as funções base utilizadas na representação da solução em cada elemento Ω𝑒.
Observando que ∫︁
Ω
𝑑x 𝜑(𝑥1, 𝑥2) ≡
𝑁𝑒∑︁
𝑒
∫︁
Ω𝑒
𝑑x 𝜑(𝑥1, 𝑥2) (2.6)
e que as funções base são linearmente independentes em cada elemento Ω𝑒, segue que
∫︁
Ω𝑒
𝑑x𝑙𝑖(𝑥1)𝑙𝑗(𝑥2)
[︃
𝜅(x)𝜕
2𝑝(x, 𝑡)
𝜕𝑡2
− ∇ ·
(︃
1
𝜌(x)∇𝑝(x, 𝑡)
)︃
− 𝑞(x, 𝑡) + ∇ ·
(︃
1
𝜌(x)f(x, 𝑡)
)︃]︃
= 0 .(2.7)
5
Substituindo a aproximação polinomial para o campo de onda em cada elemento na
equação 2.7
𝑁∑︁
𝑟=0
𝑁∑︁
𝑠=0
{︃[︂∫︁
Ω𝑒
𝑑x 𝜅(𝑥1, 𝑥2) 𝑙𝑖(𝑥1)𝑙𝑗(𝑥2) 𝑙𝑟(𝑥1)𝑙𝑠(𝑥2)]︂
𝑑2𝑝𝑟𝑠(𝑡)
𝑑𝑡2
+[︃∫︁
Ω𝑒
𝑑x 𝑙𝑖(𝑥1)𝑙𝑗(𝑥2) ∇ ·
(︃
1
𝜌(x)∇ [𝑙𝑟(𝑥1)𝑙𝑠(𝑥2)]
)︃]︃
𝑝𝑟𝑠(𝑡)
}︃
−
[︂∫︁
Ω𝑒
𝑑x 𝑙𝑖(𝑥1)𝑙𝑗(𝑥2)𝑞(x, 𝑡)
]︂
+
[︃∫︁
Ω𝑒
𝑑x 𝑙𝑖(𝑥1)𝑙𝑗(𝑥2)∇ ·
(︃
1
𝜌(x)f(x, 𝑡)
)︃]︃
= 0 , (2.8)
para todas as bases de interpolação 𝑙𝑖(𝜉1)𝑙𝑗(𝜉2) no elemento Ω𝑒. Esta expressão pode ser
reescrita como um sistemas de equações diferenciais ordinárias para os valores os do campo
onda em cada nó da malha, 𝑝𝑖𝑗(𝑡),
𝑁∑︁
𝑟=0
𝑁∑︁
𝑠=0
[︃
𝐾𝑖𝑗𝑟𝑠
𝑑2𝑝𝑟𝑠(𝑡)
𝑑𝑡2
+𝐵𝑖𝑗𝑟𝑠𝑝𝑟𝑠(𝑡)
]︃
− 𝑞𝑖𝑗(𝑡) + 𝑑𝑖𝑗(𝑡) = 0 , (2.9)
na qual
𝐾𝑖𝑗𝑟𝑠 ≡
∫︁
Ω𝑒
𝑑x𝜅(𝑥1, 𝑥2)𝑙𝑖(𝑥1)𝑙𝑗(𝑥2)𝑙𝑟(𝑥1)𝑙𝑠(𝑥2) , (2.10)
𝐵𝑖𝑗𝑟𝑠 ≡
∫︁
Ω𝑒
𝑑x𝑙𝑖(𝑥1)𝑙𝑗(𝑥2)∇ ·
(︃
1
𝜌(x)∇ [𝑙𝑟(𝑥1)𝑙𝑠(𝑥2)]
)︃
, (2.11)
𝑞𝑖𝑗(𝑡) ≡
∫︁
Ω𝑒
𝑑x𝑙𝑖(𝑥1)𝑙𝑗(𝑥2)𝑞(x, 𝑡) , (2.12)
𝑑𝑖𝑗(𝑡) ≡
∫︁
Ω𝑒
𝑑x𝑙𝑖(𝑥1)𝑙𝑗(𝑥2)∇ ·
(︃
1
𝜌(x)f(x, 𝑡)
)︃
. (2.13)
A expressão para o arranjo 𝐵𝑖𝑗𝑟𝑠 pode ser reescrita aplicando-se integração por partes,
𝐵𝑖𝑗𝑟𝑠 =
∫︁
𝜕Ω𝑒
𝑑𝑆
1
𝜌(x) 𝑙𝑖(𝑥1)𝑙𝑗(𝑥2)𝜈 · ∇ [𝑙𝑟(𝑥1)𝑙𝑠(𝑥2)]
−
∫︁
Ω𝑒
𝑑x
1
𝜌(x)∇ [𝑙𝑖(𝑥1)𝑙𝑗(𝑥2)] · ∇ [𝑙𝑟(𝑥1)𝑙𝑠(𝑥2)] . (2.14)
As integrais em cada elemento devem ser somadas de acordo com a equação 2.6. Observando
que as direções normais, 𝜈, à fronteira de elementos adjacentes são opostas, a soma das
integrais na fronteira de cada elemento se anula se as propriedades físicas são contínuas
através da fronteira; caso contrário, estas integrais devem ser avaliadas e contribuem para
a solução.
Para prosseguir a implementação do algorítmo é necessário especificar: a) a represen-
tação das propriedades físicas 𝜌(x) e 𝜅(x); b) o mapeamento entre as coordenadas do
domínio Ω𝑒 e o quadrado de referência. Uma possibilidade é especificar as propriedades
físicas em cada nó de um elemento e usar o mesmo polinômio de Lagrange para interpolar
6
a densidade e o módulo de compressibilidade em cada elemento,
𝜅(x)|Ω𝑒 ≈
𝑁∑︁
𝑖=0
𝑁∑︁
𝑗=0
𝜅𝑖𝑗𝑙𝑖(𝜉1)𝑙𝑗(𝜉2) , (2.15)
𝜌(x)|Ω𝑒 ≈
𝑁∑︁
𝑖=0
𝑁∑︁
𝑗=0
𝜌𝑖𝑗𝑙𝑖(𝜉1)𝑙𝑗(𝜉2) . (2.16)
Em nossa implementação atual optamos por manter 𝜌 e 𝜅 constantes em cada elemento e
assim reduzir o custo de armazenamento associado à representação do modelo.
O mapeamento do domínio Ω𝑒 no domínio de referência também admite escolhas. Para
representar malhas mais complexas, o mapeamento também pode ser feito a partir das
coordenadas dos nós em cada elemento (Komatitsch, 2005),
x(𝜉1, 𝜉2)|Ω𝑒 =
𝑁∑︁
𝑖=0
𝑁∑︁
𝑗=0
x𝑖𝑗 𝑙𝑖(𝜉1)𝑙𝑗(𝜉2) . (2.17)
Esta representação é bastante geral para acomodar elementos com fronteiras curvilíneas.
Por outro lado, a avaliação das integrais na equação 2.7 requer o cálculo do Jacobiano
destas transformações em cada elemento. Alternativamente, para simplificar a avaliação
do Jacobiano na integração numérica da equação 2.7, podemos exigir que os elementos
tenham lados retilíneos. Em nossa implementação atual optamos por malhas cartesianas
em que cada elemento é retangular, desta forma o mapeamento de Ω𝑒 para o domínio de
referência é ainda mais simples
𝑥1(𝜉1) =
𝑥01 + 𝑥𝑁1
2 + 𝜉1
𝑥𝑁1 − 𝑥01
2 , (2.18)
𝑥2(𝜉2) =
𝑥02 + 𝑥𝑁2
2 + 𝜉2
𝑥𝑁2 − 𝑥02
2 , (2.19)
nas quais 𝑥𝑘1 e 𝑥𝑘2 indicam coordenadas dos nós em cada elemento. Definindo
Δ1 ≡
𝑥𝑁1 − 𝑥01
2 , (2.20)
Δ2 ≡
𝑥𝑁2 − 𝑥02
2 , (2.21)
o Jacobiano deste mapeamento é simplesmente
𝐽𝑒 = Δ1 Δ2. (2.22)
A avaliação numérica das integrais na equação 2.7 utiliza a quadratura de Gauss-
Lobatto-Legendre (Deville, 2002). Em 1D a aproximação numérica para a integração por
7
quadratura de uma função 𝑓(𝜉) no intervalo de referência [−1, 1] é dada por
1∫︁
−1
𝑑𝜉𝑓(𝜉) ≈
1∫︁
−1
𝑑𝜉
𝑁∑︁
𝑖=0
𝑓(𝜉𝑖)𝑙𝑖(𝜉) =
𝑁∑︁
𝑖=0
𝑤𝑖𝑓(𝜉𝑖) , (2.23)
na qual a escolha dos pontos de avaliação da função 𝜉𝑖 e os pesos da quadratura 𝑤𝑖 é o
que definem a quadratura de Gauss-Lobatto-Legendre de ordem 𝑁 . Para esta quadratura
os pontos 𝜉𝑖 são as raízes da equação
(1 − 𝜉2) 𝑑
𝑑𝜉
𝐿𝑁−1(𝜉) = 0 , (2.24)
na qual 𝐿𝑁(𝜉) representa o polinômio de Legendre de ordem 𝑁 . Os pesos 𝑤𝑖 são obtidos
da expressão (Deville, 2002):
𝑤𝑖 =
⎧⎨⎩
2
𝑁(𝑁−1) , se 𝜉𝑖 = ±1
2
𝑁(𝑁−1)[𝐿𝑁−1(𝜉𝑖)]2
, se − 1 < 𝜉𝑖 < 1 .
(2.25)
Os polinômios de Legendre são definidos no intervalo 𝜉 ∈ [−1, 1] e podem ser calculados
por recursão (Arfken, 2005):
𝐿0(𝜉) = 1 , (2.26)
𝐿1(𝜉) = 𝜉 , (2.27)
𝐿𝑁(𝜉) =
1
𝑁
[(2𝑁 − 1)𝜉𝐿𝑁−1(𝜉) − (𝑁 − 1)𝐿𝑁−2(𝜉)] , 𝑁 ≥ 2 . (2.28)
As derivadas dos polinômios de Lagrange avaliadas nos pontos da quadratura também
são necessárias para avaliar as integrais na equação 2.31. Para polinômios de grau par,
pode-se verificar (Deville, 2002) que estas derivadas são:
𝑙′𝑖(𝜉𝑗) ≡
𝑑
𝑑𝜉
𝑙𝑖(𝜉𝑗) =
⎧⎪⎪⎨⎪⎪⎩
𝑁∑︀
𝑘=0,𝑘 ̸=𝑖
1
𝜉𝑖−𝜉𝑘
, se 𝑗 = 𝑖 ,
1
𝜉𝑗−𝜉𝑖
𝐿𝑁 (𝜉𝑗)
𝐿𝑁 (𝜉𝑖) , se 𝑗 ̸= 𝑖 .
(2.29)
Para 𝑁 = 2, os pontos para avaliação da integral por quadratura são 𝜉0 = −1, 𝜉1 = 0
e 𝜉2 = 1, e os pesos correspondentes são 𝑤0 = 1/3, 𝑤1 = 4/3 e 𝑤2 = 1/3.
2.1.1 Matrizes de compliância e flutuabilidade
Nossa escolha para o mapeamento de Ω𝑒 no domínio de referência, equações 2.18 e
2.19 , e a atribuição de propriedades físicas homogêneas para cada elemento simplifica a
avaliação dos arranjos 𝐾𝑖𝑗𝑟𝑠 e 𝐵𝑖𝑗𝑟𝑠. O arranjo de compliância 𝐾𝑖𝑗𝑟𝑠 em cada elemento é
8
dado pela expressão:
𝐾𝑒𝑖𝑗𝑟𝑠 =
⎧⎪⎪⎨⎪⎪⎩
𝜅𝑒𝐽𝑒
[︃
1∫︀
−1
𝑑𝜉1𝑙
2
𝑖 (𝜉1)
]︃ [︃
1∫︀
−1
𝑑𝜉2𝑙
2
𝑗 (𝜉2)
]︃
, se 𝑖 = 𝑟, 𝑗 = 𝑠
0, caso contrário ,
(2.30)
em que 𝜅𝑒 representa o módulo de compressibilidade do elemento. Para arranjo 𝐵𝑖𝑗𝑟𝑠,
associado a flutuabilidade (inverso da densidade) do elemento, obtém-se
𝐵𝑒𝑖𝑗𝑟𝑠 =
1
𝜌𝑒
𝐽𝑒
⎧⎨⎩ 1Δ21
⎡⎣ 1∫︁
−1
𝑑𝜉1𝑙
′
𝑖(𝜉1)𝑙′𝑟(𝜉1)
⎤⎦⎡⎣ 1∫︁
−1
𝑑𝜉2𝑙
2
𝑗 (𝜉2)
⎤⎦
+ 1Δ22
⎡⎣ 1∫︁
−1
𝑑𝜉2𝑙
′
𝑗(𝜉2)𝑙′𝑠(𝜉2)
⎤⎦ ⎡⎣ 1∫︁
−1
𝑑𝜉1𝑙
2
𝑖 (𝜉1)
⎤⎦ ,
⎫⎬⎭ (2.31)
em que 𝜌𝑒 é a densidade do elemento e 𝑙′(𝜉) indica a derivada em relação ao argumento.
Todas as integrais nas equações 2.30 e 2.31 são avaliadas por quadratura de Gauss-Lobatto-
Legendre.
2.1.2 Injeção das fontes
Para uma fonte puntual de injeção de volume,
𝑞(x, 𝑡) = 𝑞(𝑡)𝛿(x − x𝑠), (2.32)
na qual 𝑞(𝑡) é a taxa de injeção de volume e x𝑠 o ponto de injeção da fonte; a aplicação
da fonte em cada nó (m,n) do elemento que contém a fonte, Ω𝑠, pode ser avaliada
analiticamente,
∫︁
Ω𝑠
𝑞(𝑡)𝛿(x − x𝑠)𝑙𝑚(𝑥1)𝑙𝑛(𝑥2)𝑑𝑥1𝑑𝑥2 = 𝑞(𝑡)𝑙𝑚(𝑥𝑠1)𝑙𝑛(𝑥𝑠2) . (2.33)
Analogamente, para uma fonte unidirecional puntual, polarizada na direção n ≡ (𝑛1, 𝑛2),
com intensidade 𝑓(𝑡),
f(x, 𝑡) = 𝑓(𝑡)n𝛿(x − x𝑠), (2.34)
9
a injeção da fonte nos nós da malha requer a avalição da integral
∫︁
Ω𝑠
∇ ·
(︃
1
𝜌(x)𝑓(𝑡)n𝛿(x − x𝑠)
)︃
𝑙𝑚(𝑥1)𝑙𝑛(𝑥2)𝑑𝑥1𝑑𝑥2
= 𝑓(𝑡)
∫︁
Ω𝑠
[︃
∇ ·
(︃
1
𝜌(x)n𝛿(x − x𝑠)𝑙𝑚(𝑥1)𝑙𝑛(𝑥2)
)︃
− 1
𝜌(x)𝛿(x − x𝑠)n · ∇(𝑙𝑚(𝑥1)𝑙𝑛(𝑥2))
]︃
𝑑𝑥1𝑑𝑥2
= − 𝑓(𝑡)
𝜌(x𝑠)
[𝑙′𝑚(𝑥𝑠1)𝑙𝑛(𝑥𝑠2)𝑛1 + 𝑙𝑚(𝑥𝑠1)𝑙′𝑛(𝑥𝑠2)𝑛2] ,
(2.35)
na qual estamos pressupondo que a fonte está contida dentro do elemento. Deve se observar
as coordenadas (𝑥𝑠1, 𝑥𝑠2) no elemento são mapeadas em (𝜉𝑠1, 𝜉𝑠2) no domínio de referência
para avaliação dos polinômios de Lagrange e suas derivadas.
2.1.3 Construção do sistema linear global
Após os cálculos dos arranjos 𝐾𝑖𝑗𝑟𝑠 e 𝐵𝑖𝑗𝑘𝑠 em cada elemento e dos termos fontes 𝑞𝑖𝑗 e
𝑑𝑖𝑗 é necessário fazer a soma sobre todos os elementos, como indicado pela equação 2.6.
Para poder efetuar esta soma sobre todo o domínio é necessário mapear os índices dos nós
locais em cada elemento Ω𝑒, os quais estão indexados por um par de índices (𝑖𝑗), em um
índice 𝐼 que referencia o nó de forma global na malha. Este mapeamento é implementado
através da função,
𝐼 ≡ local_global(𝑖, 𝑗, 𝑒) . (2.36)
Esta função é chamada dentro do laço que efetua a soma sobre todos os elementos do
domínio Ω para montar o sistema linear especificado pela indexação global dos nós da
malha.Efetuada a soma sobre todos os elementos da malha para cada nó com a indexação
global, obtém-se o sistema linear:
𝐾𝐼𝐼
𝑑2
𝑑𝑡2
𝑝𝐼(𝑡) −
𝑁𝐺∑︁
𝑅=1
𝐵𝐼𝑅 𝑝𝑅(𝑡) = 𝑞𝐼(𝑡) − 𝑑𝐼(𝑡) , 𝐼 ∈ {1, · · · , 𝑁𝐺} , (2.37)
em que, 𝑁𝐺 representa o número de nós em todo o domínio Ω. A matriz 𝐾𝐼𝐼 é diagonal,
em consequência da ortogonalidade do polinômios de Lagrange (Deville, 2002), conforme
indica a equação 2.30. Esta é uma das propriedades importantes desta metodologia em
relação a outras abordagens de elementos finitos (Komatitsch, 2010). Em consequência, o
sistema linear a ser resolvido pode ser reescrito explicitamente,
𝑑2
𝑑𝑡2
𝑝𝐼(𝑡) =
1
𝐾𝐼𝐼
⎡⎣ 𝑁𝐺∑︁
𝑅=1
𝐵𝐼𝑅 𝑝𝑅(𝑡) + 𝑞𝐼(𝑡) − 𝑑𝐼(𝑡)
⎤⎦ , 𝐼 ∈ {1, · · · , 𝑁𝐺}. (2.38)
10
Para integração no tempo, usamos o esquema de diferenças finitas de segunda ordem:
𝑝𝐼(𝑡+ Δ𝑡) = 2𝑝𝐼(𝑡) − 𝑝(𝑡− Δ𝑡) +
(Δ𝑡)2
𝐾𝐼𝐼
⎡⎣ 𝑁𝐺∑︁
𝑅=1
𝐵𝐼𝑅 𝑝𝑅(𝑡) + 𝑞𝐼(𝑡) − 𝑑𝐼(𝑡)
⎤⎦ . (2.39)
A análise teórica de condições de estabilidade para se determinar o passo de evolução
no tempo Δ𝑡 não é tão simples como para o caso de esquemas de diferenças finitas e não
há uma metodologia estabelecida na literatura (Komatitsch, 2005). Nesta implementação
adotamos a prescrição mais conservadora proposta por Komatitsch (2005):
Δ𝑡 < 0.3 Δℎ
𝑐𝑚𝑎𝑥
, (2.40)
em que Δℎ indica a menor razão entre o menor espaçamento entre pontos na malha e 𝑐𝑚𝑎𝑥
reprepresenta a velocidade máxima de propagação no modelo.
2.1.4 Condições de fronteira
O sistema linear na equação 2.39 é singular. Para garantir a unicidade da solução
é necessário especificar condições de contorno na fronteira de Ω. Na implementação
aplicamos a condição de superfície livre em toda a fronteira, ou seja,
𝑝(x, 𝑡)|𝜕Ω = 0. (2.41)
Esta condição é implementada atribuindo um valor arbitrariamente grande à entrada da
matriz 𝐾𝐼𝐼 associadas aos nós na fronteira do domínio (Deville, 2002).
Para evitar reflexões espúrias nas laterais e na base do modelo, adicionamos nestas
regiões com fronteiras absorventes do tipo CPML (Komatitsch, 2007). Há também a opção
de se aplicar fronteiras absorventes no topo do modelo, se a condição de superfície livre
não for desejada.
Condições de fronteira CPML são implementadas através do sistema de equações:
𝜅(x)𝜕
2𝑝(x, 𝑡)
𝜕𝑡2
− ∇ ·
(︃
1
𝜌(x)∇𝑝(x, 𝑡)
)︃
+
𝐷∑︁
𝜈=1
{︃
𝜕
𝜕𝑥𝜈
(︃
1
𝜌(x)𝜓𝜈(x, 𝑡)
)︃
+ Ψ𝜈(x, 𝑡)
}︃
= 𝑞(x, 𝑡) − ∇ ·
(︃
1
𝜌(x)f(x, 𝑡)
)︃
𝜕𝜓𝜈(x, 𝑡)
𝜕𝑡
+ 𝛾𝜈(x)𝜓𝜈(x, 𝑡) = 𝜎𝜈(x)
𝜕𝑝(x, 𝑡)
𝜕𝑥𝜈
(2.42)
𝜕Ψ𝜈(x, 𝑡)
𝜕𝑡
+ 𝛾𝜈(x)Ψ𝜈(x, 𝑡) = 𝜎𝜈(x)
𝜕
𝜕𝑥𝜈
(︃
1
𝜌(x)
(︃
𝜕𝑝(x, 𝑡)
𝜕𝑥𝜈
− 𝜓𝜈(x, 𝑡)
)︃)︃
,
nas quais os campos auxiliares 𝜓𝜈(x, 𝑡) e Ψ𝜈(x, 𝑡) existem apenas nas regiões com as fron-
teiras absorventes e 𝐷 indica a dimensão do domínio. Fora destas regiões, as propriedades
11
que controlam a absorção do campo de onda, 𝛾𝜈(x) e 𝜎𝜈(x), são nulas.
3 RESULTADOS E DISCUSSÃO
Neste capítulo avaliamos a acurácia das soluções numéricas da equação da onda acústica
usando a nossa implementação do métodos de elementos finitos espectrais. O algoritmo
de elementos finitos espectrais foi implementado em FORTRAN2003 (Metcalf, 2011). Para
acelerar a performance do código usamos apenas diretivas openMP (Chapman, 2008). A
seguir apresentamos alguns exemplos numéricos para solução da equação acústica com
densidade variável. Estes experimentos indicam as características da implementação em
relação a acurácia, inclusão de fontes de injeção de volume e fontes dipolares, estabilidade
da condição de fronteira livre e eficácia das fronteiras absorventes na supressão de efeitos
de borda.
O primeiro teste compara a solução analítica da equação da onda com a solução
numérica em um meio homogêneo. A solução analítica foi computada pela convolução de
um pulso Ricker com frequência pico de 8 Hz com a função de Green para a propagação
de onda em 2D,
𝐺(x, 𝑡; x𝑠) =
1
2𝜋𝑐2
𝐻(𝑡− ‖x − x𝑠‖/𝑐)√︁
𝑡2 − (‖x − x𝑠‖/𝑐)2
. (3.1)
As soluções numéricas foram computadas em um modelo com velocidade de propagação
de 2500 m/s. A malha de elementos finitos é uniforme, contém 251×767 elementos de
dimensão 12 m × 12 m; o campo de onda foi aproximado por polinômios interpoladores de
graus 6 e 7. A Figura 3.1 apresenta as seções sísmicas obtidas com a solução analítica e
com elementos finitos de ordem 6, as seções foram normalizadas pela máxima amplitude. A
Figura 3.2 apresenta o resíduo entre a solução analítica e as soluções de elementos finitos de
grau 6 e 7, utilizando a mesma escala de amplitudes que a 3.1. Estes resultados mostram
a boa qualidade das soluções calculadas por elementos finitos de alta ordem; também se
observa que a amplitude dos resíduos apresentam anomalias de maior amplitude para
alguns afastamentos, este efeito pode estar associado a avaliação numérica da função de
Green que apresenta uma singularidade para 𝑡 = ‖x − x𝑠‖/𝑐.
Para avaliar a implementação de fontes dipolares e o desempenho da condição de
fronteira livre utilizamos o modelo simples com velocidade e densidade mostradas nas
Figuras 3.3 e 3.4, respectivamente. Este modelo não apresenta contraste de impedância
entre a base da camada superior e as regiões subjacentes. A simulação foi feita com uma
fonte dipolar localizada diretamente acima da interface vertical que separa as regiões
subjacentes, nas coordenadas (4608 m, 750 m). A orientação do dipolo forma 45∘ com a
vertical. O modelo foi representado em uma malha com elementos quadrados uniformes
de 12 m de lado. Em cada elemento, o campo de onda é interpolado usando o produto
tensorial de polinômios de Lagrange de grau 4, como ilustra a Figura 2.1. Fronteiras
absorventes do tipo CPML foram aplicadas. Em todas as simulações estas condições
12
13
Figura 3.1: Seção sísmica para o campo de pressão calculada analiticamente, à esquerda;
seção para o campo de pressão obtida usando elementos finitos espectral de ordem 6, a
direita.
Figura 3.2: Seções sísmicas com o resíduo entre a solução analítica e calculada usando
elementos finitos espectral de: ordem 6 (à esquerda’) e ordem 7 (à direita).
foram implementadas estendendo-se as fronteiras laterais e da base do modelo com bordas
contendo 25 elementos. A fonte injeta um pulso Ricker com frequência pico de 12 Hz. A
seguir apresentamos três instantâneos do campo de onda desta simulação.
A Figura 3.5 mostra o campo de onda logo após o fim da injeção do pulso fonte.
Observa-se que o campo de onda apresenta o padrão de radiação correto para uma fonte
dipolar em um meio homogêneo. Este tipo de fonte seria mais difícil de simular com a
mesma acurácia em um esquema de diferenças finitas. O instantâneo do campo de onda na
Figura 3.6 destaca a acurácia da condição de superfície livre. Pode-se observar a inversão
de polaridade do campo de onda refletido na superfície em relação à onda incidente e a
variação da amplitude do campo ao longo deste evento resultante do padrão de radiação
da fonte dipolar. A Figura 3.7 mostra o campo de onda após a frente de onda atravessar a
base do modelo, podemos observar que a fronteira absorvente na base do modelo foi capaz
14
de suprimir reflexões espúrias nesta fronteira.
Figura 3.3: Modelo de velocidade, escala de cores em m/s.
Figura 3.4: Modelo de densidade, escala de cores em kg/m3.
Figura 3.5: Instantâneo do campo de onda logo após a ativação da fonte no meio com
velocidade e densidade indicadas nas Figuras 3.3 e 3.4, respectivamente. Observa-se o
padrão de radiação da fonte em uma região homogênea.
A seguir avaliamos a variação do campo de onda registrado para diferentes graus do
polinômio interpolador. A Figura 3.8 mostras as duas seções símicas obtidas usando
15
Figura 3.6: Instantâneo do campo de onda após a reflexão da frente de onda na superfície
livre. A velocidade e a densidade do modelo estão indicadas nas Figuras 3.3 e 3.4,
respectivamente.
Figura 3.7: Instantâneo do campo de onda após a frente de onda ter passado pela base
do modelo. As fronteirasabsorventes foram eficazes para suprimir reflexões espúrias. A
velocidade e a densidade do modelo estão indicadas nas Figuras 3.3 e 3.4, respectivamente.
polinômios interpoladores de grau 4 e 7, e a diferença entre as duas seções. Os valores de
amplitude estão na mesma escala com truncamento em 1% do máximo valor da amplitude
para solução de grau 7. Este resultado indica que há ganho de acurácia com o aumento
da ordem do polinômio, entretanto, o custo de armazenamento com campo de onda e o
número de operações de ponto flutuante para atualização do campo de onda escala com
(𝑁 + 1)2, em que 𝑁 indica o grau do polinômio; adicionalmente, o número de atualizações
do campo de onda para um mesmo tempo de simulação escala com 𝑁 , ou seja, o custo
total de operações de ponto flutuante escala com 𝑁(𝑁 + 1)2. Neste exemplo numérico, a
diferença entre as duas soluções não justifica o custo de se utilizar o polinômio de grau 7.
Modelo Marmousi
A seguir apresentamos simulações da propagação de ondas no modelo Marmousi. A
velocidade de propagação e a densidade para este modelo são apresentados, respectivamente,
16
Figura 3.8: Seções sísmicas calculadas usando elementos finitos espectral de: ordem 4 (à
esquerda) e ordem 7 (à direita); diferença entre as duas seções(abaixo).
nas Figuras 3.9 e 3.10. O modelo foi também representado por elementos quadrados
uniformes com 12 m de lado e o campo de onda interpolado com polinômios de Lagrange
de grau 4. Nos experimentos a seguir utilizamos um pulso fonte Ricker com frequência
pico igual a 15 Hz.
Inicialmente avaliamos a estabilidade do algoritmo na presença de superfície livre. O
tempo de simulação neste experimento foi de 4.5 s. Durante este intervalo de simulação
não houve ocorrência de instabilidade. Para destacar as características da simulação com
superfície livre, efetuamos uma segunda simulação sem a presença de superfície livre. As
Figuras 3.11 e 3.12 apresentam o instantâneo do campo de onda após 1.4 s de propagação
para as simulações com e sem a superfície livre, respectivamente. Comparando-se estas
Figuras, particularmente a frente de onda, observa-se a mudança da forma de onda causada
pela condição de superfície livre; também se verifica a maior variabilidade do campo de
onda causada pela interferência de eventos ascendentes e descendentes, particularmente
em profundidade abaixo de 1 km. Verificamos diferenças na amplitude dos eventos acima
de 1 km provavelmente causadas pelas diferenças no padrão de radiação da fonte dipolar
em relação à fonte de injeção de volume.
A Figura 3.13 destaca a qualidade das condições de absorção CPML. Nesta simulação a
17
fonte de injeção de volume está localizada próximo à fronteira lateral esquerda do modelo
nas coordenadas (1500 m, 100 m). A condição de superfície livre está sendo aplicada.
Novamente, não se observam eventos espúrios associados a reflexões na fronteira lateral ou
na base do modelo após 1.4 s de propagação.
Figura 3.9: Modelo de velocidade Marmousi, escala de cores em m/s.
Figura 3.10: Modelo de densidade Marmousi, escala de cores em kg/m3.
Figura 3.11: Instantâneo do campo de onda após 1.4 s de propagação. Fonte de injeção de
volume nas coordenadas (4500 m, 15 m), com superfície livre.
18
Figura 3.12: Instantâneo do campo de onda após 1.4 s de propagação. Fonte de injeção de
volume nas coordenadas (4500 m, 15 m), sem superfície livre.
Figura 3.13: Instantâneo do campo de onda após 1.4 s de propagação. Fonte de injeção de
volume, localizada nas coordenadas (1500 m, 100 m), com condição de superfície livre no
topo. Observa-se a qualidade das condições de absorção CPML, não há efeitos de borda
significativos na lateral e na base do modelo.
Finalmente, para o modelo Marmousi, avaliamos a acurácia com que as relações de
reciprocidade para o campo de pressão e seu gradiente, apresentadas no Apêndice 4, são
obedecidas pelas soluções numéricas. Para este conjunto de experimentos utilizamos um
pulso fonte Ricker com frequência pico de 8 Hz. As Figuras 3.14, 3.15 e 3.16 apresentam
seções sísmicas do campo de pressão e as componentes do gradiente do campo comparando
experimentos com e sem superfície livre. A fonte está localizada a 18 m de profundidade e
posição 4605 m, e 181 receptores estão a profundudade de 15 m a intervalos de 50 m. A
soluções foram computadas com polinômio interpolador de grau 4.
Para avaliar a acurácia das soluções numéricas em um meio heterogêneo, avaliamos as
relacões de reciprocidade para um traço sísmico com fonte localizada na posição 3000 m e
profundide 15 m e receptor nas coordenadas 6000 m e 10 m, com condição de superfície-
livre. A Figuras 3.17 e 3.18 apresentam os pares de traços para cada uma das relações
de reciprocidade envolvento fontes monopolares e fontes dipolares. A superposição dos
traços recíprocos mostra a acurácia com que a propagação e as condições de fronteira são
19
Figura 3.14: Campo de pressão: sem superfície livre(topo), com superfície-livre (base).
honradas pelo método de elementos finitos espectral, com polinômios interpoladores de
ordem 4. A acurácia com que estas relações de reciprocidade são obedecidas pode ser
ainda maior se aumentarmos o grau do polinômio interpolador. Contudo, considerando
como o custo computacional aumenta rapidamente com o grau do polinômio, avalio que
estas soluções tem acurácia suficiente para as aplicações em sismica de exploração.
O uso de malhas reglares em todos os experimentos numéricos permite uma melhor
comparação dos custos do SEM em relação ao método de diferenças finitas amplamente
utilizado pela indústria. Apesar do custo de armazenamento e custo computacional do
SEM escalar muito desfavoravelmente, há aplicações que justificam seu uso na sísmica de
exploração. Por exemplo, na aquisição sísmica com receptores no fundo do relevo oceânico;
neste caso o uso de malhas quadrangulares não uniformes, que representem melhor o
20
Figura 3.15: Gradiente vertical do campo de pressão: sem superfície livre(topo), com
superfície-livre (base).
relevo do fundo oceânico, elimina eventos espúrios de difração, causado pelos degraus na
representação do fundo do mar presente nas malhas de diferenças finitas. Estes eventos
espúrios podem comprometer seriamente a acurácia da inversão da forma de onda quando
os receptores estão no fundo do mar e as malhas de representação são mais esparsas devido
à menor banda de frequência do pulso sísmico (Trinh, 2017). Adiconalmente, na presença
de variações topográficas o SEM é bem mais adequado que a abordagem tradicional por
diferenças finitas. Finalmente, na presença de refletores com forte contraste de propriedades
físicas como, por exemplo, em flancos de domos de sal e falhas, as amplitudes do campo
de onda são computadas com muito mais acurácia em malhas irregulares que melhor
representem estas interfaces.
21
Figura 3.16: Gradiente horizontal do campo de pressão: sem superfície livre(topo), com
superfície-livre (base).
22
Figura 3.17: Relação de reciprocidade entre fontes monopolares(topo);relações de recipro-
cidade entre fonte monopolar e dipolo vertical(meio); relações de reciprocidade entre fonte
comopolar e dipolo horizontal(base). Simulação com superfície livre.
23
Figura 3.18: Relação de reciprocidade entre fontes do tipo dipolos verticais(topo);relações
de reciprocidade entre dipolo vertical e dipolo horizontal(meio); relações de reciprocidade
entre dipolos horizontais(base). Simulação com superfície livre.
4 CONCLUSÕES
Esta dissertação descreve a metodologia do método de elementos espectrias e os detalhes
de sua implementação em fortran90. A validação da implementação foi feita, inicialmente,
comparando os resultados da solução SEM com a solução analítica em meios homogêneos.
O resíduo entre a solução numérica computada com polinômio interpolador de grau 6
e a solução computada convento o pulso fontecom a função de Green 2D apresentam
valores aceitáveis para as aplicações em sísmica de exploração e podem ser atribuidos a
inacurácia na representaçãoda singularidade na função de Green. A seguir avaliamos,
em um meios heterogêneo a acurácia da representação de fontes dipolares e condições de
fronteira do tipo superfície livre, comparando o resíduo entre soluções computadas com
polinômio interpolador do campo de onda de graus 4 e 7. Observamos que o custo de
armazenamento escala com o quadrado do polinômio interpolador e o custo de operações
de ponto flutudante escala com o cudo da ordem do polinômio. Os exemplos numéricos
indicam que polinômios de ordem 4 apresentam um bom compromisso entre acurácia e
custo computacioal. Finalmente, para avaliar a acurácia do algoritmo em um modelo
heterogêneo mais complexo avaliamos a acurácia das todas as relações de reciprocidade
para o campo acústico e seu gradiente no modelo Marmousi com superfície livre. Os
resultados reforçam mais uma vez o bom compromisso entre acurácia e custo computacional
com polinômio interpolador de ordem 4. Devido ao maior custo computacional do SEM
em relação ao método de diferenças finitas, sua aplicação se justifica em contextos de
modelagem, imageamento e inversão na presença de variações topográficas, variação de
relevo do fundo oceânico em aquições tipo nodes e OBC, e para modelagem acurada da
amplitude do campo de onda na presença de interfaces com geometria irregular separando
domínios com forte contraste de propriedades físicas.
24
REFERÊNCIAS
Arfken, G.B., H.J. Weber, 2005, Mathematical methods for physicists international student
edition. Elsevier. 7
Chapman, B., G. Jost, R. Van Der Pas, 2008, Using OpenMP: portable shared memory
parallel programming. MIT press. 12
Deville, M.O., P.F. Fischer, E.H. Mund and others, 2002, High-order methods for incom-
pressible fluid flow, vol. 9. Cambridge university press. 3, 4, 6, 7, 9, 10
Fichtner, A., 2010, Full seismic waveform modelling and inversion. Springer Science &
Business Media. 1
Fokkema, J.T., P.M. van den Berg, 2013, Seismic applications of acoustic reciprocity.
Elsevier. 1
Ismail-Zadeh, A., P. Tackley , 2010, Computational methods for geodynamics. Cambridge
University Press. 1
Komatitsch, D., G. Erlebacher, D. Göddeke, D. Michéa, 2010, High-order finite-element
seismic wave propagation modeling with mpi on a large gpu cluster. Journal of compu-
tational physics, 229(20), 7692–7714. 9
Komatitsch, D., R. Martin, 2007, An unsplit convolutional perfectly matched layer improved
at grazing incidence for the seismic wave equation. Geophysics, 72(5), SM155–SM167.
10
Komatitsch, D., S. Tsuboi, J. Tromp, 2005, The spectral-element method in seismology.
Seismic earth: Array analysis of broadband seismograms, 205–227. 3, 6, 10
Metcalf, M., J. Reid, M. Cohen, 2011, Modern Fortran Explained. Oxford University Press.
12
Snieder, R.K., 1998, A guided tour of mathematical physics. Samizdat Press. 2
Trinh, P.T., R. Brossier, L. Métivier, L. Tavard, J. Virieux, P. Wellington, 2017, Efficient 3d
elastic fwi using a spectral-element method on cartesian-based mesh. In: SEG Technical
Program Expanded Abstracts 2017. Society of Exploration Geophysicists, 1533–1538. 3,
20
25
APÊNDICE A
RELAÇÕES DE RECIPROCIDADE PARA O CAMPO DE ONDA
As relações de reciprocidade expressão simetrias espacias e temporais do campo de
ondas. Nesta dissertação utilizamos as simetrias espaciais para avaliar a acurácia das
soluções numéricas em meios heterogêneos, para os quais não há soluções analíticas para a
equação da onda acústica.
As duas soluções fundamentais para a equação de onda em meios acústicos são:
i. a solução para uma fonte de injeção de volume puntual impulsiva (monopolo):
𝜅(x)𝜕
2𝐺(x, 𝑡; x′)
𝜕𝑡2
− ∇ ·
(︃
1
𝜌(x)∇𝐺(x, 𝑡; x
′)
)︃
− 𝛿(𝑡)𝛿(x − x′) = 0 (4.1)
ii. a solução para uma fonte unidirecional putual impulsiva (dipolo):
𝜅(x)𝜕
2Γ𝑛(x, 𝑡; x′)
𝜕𝑡2
−∇·
(︃
1
𝜌(x)∇Γ𝑛(x, 𝑡; x
′)
)︃
+∇·
(︃
1
𝜌(x)𝛿(𝑡)ê𝑛𝛿(x − x
′)
)︃
= 0 (4.2)
Considerando meios estacionários podemos representar as duas soluções acima no
domínio de Fourier,
𝐺(x, 𝑡; x′) = 12𝜋
∫︁ +∞
−∞
�̃�(x, 𝜔; x′)𝑒−𝑖𝜔𝑡𝑑𝜔 , (4.3)
Γ𝑛(x, 𝑡; x′) =
1
2𝜋
∫︁ +∞
−∞
Γ̃𝑛(x, 𝜔; x′)𝑒−𝑖𝜔𝑡𝑑𝜔 . (4.4)
Considerando ainda que
𝛿(𝑡) = 12𝜋
∫︁ +∞
−∞
𝑒−𝑖𝜔𝑡𝑑𝜔 (4.5)
e substituindo estas representações na equação de onda, seguem as correspondentes
equações de onda no domínio da frequência:
𝜅(x) − 𝜔2 �̃�(x, 𝜔; x′) − ∇ ·
(︃
1
𝜌(x)∇�̃�(x, 𝜔; x
′)
)︃
− 𝛿(x − x′) = 0 (4.6)
e
𝜅(x) − 𝜔2 Γ̃𝑛(x, 𝜔; x′) − ∇ ·
(︃
1
𝜌(x)∇Γ̃𝑛(x, 𝜔; x
′)
)︃
+ ∇ ·
(︃
1
𝜌(x) ê𝑛𝛿(x − x
′)
)︃
= 0 . (4.7)
As relações de reciprocidade espaciais são obtidas considerando as sequintes configura-
ções de aquisição: �̃�(x, 𝜔; x′), �̃�(x, 𝜔; x′′), Γ̃𝑛(x, 𝜔; x′) e Γ̃𝑚(x, 𝜔; x′′). Para um domínio
Ω limitado por uma fronteira 𝜕Ω considerar as seguintes interações:
26
27
i. ∫︁
Ω
[︁
�̃�(x, 𝜔; x′)(−𝜔2)�̃�(x, 𝜔; x′′) − �̃�(x, 𝜔; x′′)(−𝜔2)�̃�(x, 𝜔; x′)
]︁
𝑑Ω
,
ii. ∫︁
Ω
[︁
Γ̃𝑛(x, 𝜔; x′)(−𝜔2)�̃�(x, 𝜔; x′′) − �̃�(x, 𝜔; x′′)(−𝜔2)Γ̃𝑛(x, 𝜔; x′)
]︁
𝑑Ω
,
iii. ∫︁
Ω
[︁
Γ̃𝑛(x, 𝜔; x′)(−𝜔2)Γ̃𝑚(x, 𝜔; x′′) − Γ̃𝑚(x, 𝜔; x′′)(−𝜔2)Γ̃𝑛(x, 𝜔; x′)
]︁
𝑑Ω .
Observando que as integrais acima são identicamente nulas, utilizando as equações de
onda 4.6 e 4.7 segue que:
i.
∫︁
Ω
{︃
�̃�(x, 𝜔; x′)
[︃
∇ ·
(︃
1
𝜌(x)∇�̃�(x, 𝜔; x
′′)
)︃
+ 𝛿(x − x′′)
]︃
− �̃�(x, 𝜔; x′′)
[︃
∇ ·
(︃
1
𝜌(x)∇�̃�(x, 𝜔; x
′′)
)︃
+ 𝛿(x − x′)
]︃
𝑑Ω
}︃
= 0 ,
ii.
∫︁
Ω
{︃
Γ̃𝑛(x, 𝜔; x′)
[︃
∇ ·
(︃
1
𝜌(x)∇�̃�(x, 𝜔; x
′′)
)︃
+ 𝛿(x − x′′)
]︃
− �̃�(x, 𝜔; x′′)
[︃
∇ ·
(︃
1
𝜌(x)∇Γ̃𝑛(x, 𝜔; x
′′) − ∇ ·
(︃
1
𝜌(x) ê𝑛𝛿(x − x
′)
)︃)︃]︃
𝑑Ω
}︃
= 0 ,
iii.
∫︁
Ω
{︃
Γ̃𝑛(x, 𝜔; x′)
[︃
∇ ·
(︃
1
𝜌(x)∇Γ̃𝑚(x, 𝜔; x
′′)
)︃
− ∇ ·
(︃
1
𝜌(x) ê𝑚𝛿(x − x
′′)
)︃]︃
− Γ̃𝑚(x, 𝜔; x′′)
[︃
∇ ·
(︃
1
𝜌(x)∇Γ̃𝑛(x, 𝜔; x
′′) − ∇ ·
(︃
1
𝜌(x) ê𝑛𝛿(x − x
′)
)︃)︃]︃
𝑑Ω
}︃
= 0 ,
Aplicando o teorema de Green e a propriedade de filtragem da distribuição delta
i.
�̃�(x′′, 𝜔; x′) − �̃�(x′, 𝜔; x′′)+∫︁
𝜕Ω
𝑑(𝜕Ω)n̂ · 1
𝜌(x)
{︁
�̃�(x, 𝜔; x′)∇�̃�(x, 𝜔; x′′)
− �̃�(x, 𝜔; x′′)∇�̃�(x, 𝜔; x′)
}︁
= 0 ,
28
ii.
Γ̃𝑛(x′′, 𝜔; x′) −
1
𝜌(x′) ê𝑛 · ∇�̃�(x
′, 𝜔; x′′)+∫︁
𝜕Ω
𝑑(𝜕Ω)n̂ · 1
𝜌(x)
{︁
Γ̃𝑛(x, 𝜔; x′)∇�̃�(x, 𝜔; x′′)
− �̃�(x, 𝜔; x′′)∇Γ̃𝑛(x, 𝜔; x′)
}︁
= 0 ,
iii.
1
𝜌(x′′) ê𝑚 · ∇Γ̃𝑛(x
′′, 𝜔; x′) − 1
𝜌(x′) ê𝑛 · ∇Γ̃𝑚(x
′, 𝜔; x′′)+∫︁
𝜕Ω
𝑑(𝜕Ω)n̂ · 1
𝜌(x)
{︁
Γ̃𝑛(x, 𝜔; x′)∇Γ̃𝑚(x, 𝜔; x′′)
− Γ̃𝑚(x, 𝜔; x′′)∇Γ̃𝑛(x, 𝜔; x′)
}︁
= 0 .
Para condições de fronteira homogêneas: superfície livre (pressão nula); superfície
rígida (gradiente de pressão nulo); e condições de radiação de Sommerfeld, aproximadas
por condições de fronteira absorventes, todas as integrais de fronteira em i., ii. e iii. se
anulam. Consequentemente, com campo de onda acústico deve satisfazer as seguintes
relações de reciprocidade:
i.
�̃�(x′′, 𝜔; x′) = �̃�(x′, 𝜔; x′′); (4.8)
ii.
Γ̃𝑛(x′′, 𝜔; x′) =
1
𝜌(x′) ê𝑛 · ∇�̃�(x
′, 𝜔; x′′) ; (4.9)
iii.
1
𝜌(x′′) ê𝑚 · ∇Γ̃𝑛(x
′′, 𝜔; x′) = 1
𝜌(x′) ê𝑛 · ∇Γ̃𝑚(x
′, 𝜔; x′′) . (4.10)
	 Introdução
	 Metodologia
	Elementos finitos espectral
	Matrizes de compliância e flutuabilidade
	Injeção das fontes
	Construção do sistema linear global
	Condições de fronteira
	 Resultados e discussão
	 Conclusões
	REFERÊNCIAS
	APÊNDICE A - RELAÇÕES DE RECIPROCIDADE PARA O CAMPO DE ONDA